Stress Analysis of Three-Dimensional Media Containing Localized Zone by FEM-SGBEM Coupling

This paper presents an efficient numerical technique for stress analysis of three-dimensional infinite media containing cracks and localized complex regions. To enhance the computational efficiency of the boundary element methods generally found inefficient to treat nonlinearities and non-homogeneous data present within a domain and the finite element method FEM potentially demanding substantial computational cost in the modeling of an unbounded medium containing cracks, a coupling procedure exploiting positive features of both the FEM and a symmetric Galerkin boundary element method SGBEM is proposed. The former is utilized to model a finite, small part of the domain containing a complex region whereas the latter is employed to treat the remaining unbounded part possibly containing cracks. Use of boundary integral equations to form the key governing equation for the unbounded region offers essential benefits including the reduction of the spatial dimension and the corresponding discretization effort without the domain truncation. In addition, all involved boundary integral equations contain only weakly singular kernels thus allowing continuous interpolation functions to be utilized in the approximation and also easing the numerical integration. Nonlinearities and other complex behaviors within the localized regions are efficiently modeled by utilizing vast features of the FEM. A selected set of results is then reported to demonstrate the accuracy and capability of the technique.


Introduction
A physical modeling of three-dimensional solid media by an idealized mathematical domain that occupies the full space is standard and widely used when inputs and responses of interest are only localized in a zone with its length scale much smaller than that of the body.Influence of the remote boundary of a domain on such responses is generally insignificant for this particular case and can, therefore, be discarded in the modeling without loss of accuracy of the predicted solutions.Such situations arise in various engineering applications such as the simulation of crack growth in hydraulic fracturing process in which the fracture is generally treated as an isolated crack in an infinite medium, the evaluation and assessment of service life of large-scale structures in which the influence of embedded initial defects can be characterized by small pre-existing flaws, and the determination of effective properties of materials possessing a microstructure such as cracks, voids, inclusions, and localized inelastic zones.Unlike the stress analysis of linear elasticity problems, complexity of a mathematical model can substantially increase when an infinite body contains additionally a line of singularity and/or a localized nonlinear region.The former situation arises naturally when a surface of displacement discontinuities e.g., cracks and dislocations is present whereas the latter may result from applications of high-intensity loads, complex constitutive laws, containment of small defects and inhomogeneities, and localized non-mechanical effects such as temperature change .Besides various practical applications, such present complexity renders the modeling itself theoretically and computationally challenging.
Various analytical techniques e.g., integral transform methods, methods based on stress and displacement representations, techniques related to potential theories, etc. have been proposed and used extensively in the stress analysis of solid media e.g., 1-5 .However, their applications are very limited to either two-dimensional boundary value problems involving simple data or three-dimensional problems with extremely idealized settings.Such limitation becomes more apparent when complexity of involved physical phenomena increases e.g., complexity introduced by the presence of material nonlinearities, inhomogeneities, and embedded singularities .For those particular situations, a more sophisticated mathematical model is generally required in order to accurately predict responses of interest, and, as a major consequence, an analytical or closed-form solution of the corresponding boundary value problem cannot readily be obtained and numerical techniques offer better alternatives in the solution procedure.
The finite element method FEM and the boundary element method BEM are two robust numerical techniques extensively used in the modeling of various field problems.Both techniques possess a wide range of applications, and there are various situations that favor the FEM over the BEM and vice versa.The FEM has proven to be an efficient and powerful method for modeling a broad class of problems in structural and solid mechanics e.g., 6-8 .In principle, the basis of the FEM is sufficiently general allowing both nonlinearities and inhomogeneities present within the domain to be treated.In addition, a final system of discrete algebraic equations resulting from this method possesses, in general, desirable properties e.g., symmetry, sparseness, positive definiteness of the coefficient matrix, etc. .Nevertheless, the conventional FEM still exhibits some major shortcomings and requires nontrivial treatments when applied to certain classes of problems.For instance, a standard discretization procedure cannot directly be applied to boundary value problems involving an infinite domain.A domain truncation supplied by a set of remote boundary conditions is commonly employed to establish an approximate domain of finite dimensions prior to the discretization.It should be noted that defining such suitable truncated surface and corresponding boundary conditions remains the key issue and it can significantly influence the quality of approximate solutions.Another limited capability of the method to attain adequately accurate results with reasonably cheap computational cost is apparent when it is applied to solve fracture problems.In the analysis, it generally requires substantially fine meshes in a region surrounding the crack front in order to accurately capture the complex singular field and extract essential local fracture information such as the stress intensity factors e.g., 9-11 .
The boundary element method BEM has been found computationally efficient and attractive for modeling certain classes of linear boundary value problems since, for a homogeneous domain that is free of distributed sources, the key governing equation involves only integrals on the boundary of the domain e.g., 12-23 .As a direct consequence, the discretization effort and cost are significantly reduced, when compared to the FEM, due to the reduction of spatial dimensions of the governing equation by one.Another apparent advantage of the method is associated with its simplicity in the modeling of an infinite media.In such particular situation, the remote boundary of the domain can basically be discarded without loss via an appropriate treatment of remote conditions e.g., 14, 17, 19, 23 .Among various strategies utilized to form the BEM, the weakly singular symmetric Galerkin boundary element method SGBEM has become a wellestablished and well-known technique and, during the past two decades, has proven robust for three-dimensional analysis of linear elasticity problems e.g., 15, 16 , linearly elastic infinite media containing isolated cracks e.g., 14, 19, 23 , and cracks in finite bodies e.g., 16,18,20,21,23 .Superior features of this particular technique over other types of the BEM are due mainly to that all kernels appearing in the governing integral equations are only weakly singular of O 1/r and that a final system of linear algebraic equations resulting from the discretization possesses a symmetric coefficient matrix.The weakly singular nature not only renders all involved integrals to be interpreted in an ordinary sense and evaluated numerically using standard quadrature but also allows standard C 0 interpolation functions to be employed in the approximation procedure.It has been also demonstrated that the weakly singular SGBEM along with the proper enrichment of an approximate field near the crack front yields highly accurate fracture data e.g., mixed-mode stress intensity factors even that relatively coarse meshes are employed in the discretization e.g., 18, 21, 23 .While the weakly singular SGBEM has gained significant success in the analysis of linear elasticity and fracture problems, the method still contains certain unfavorable features leading to its limited capability to solve various important classes of boundary value problems.For instance, the method either becomes computationally inefficient or experiences mathematical difficulty when applied to solve problems involving nonlinearity and nonhomogeneous media.As the geometry of the domain becomes increasingly complex and its size and surface to volume ratio are relatively large requiring a large number of elements to reasonably represent the entire boundary of the domain , the method tends to consume considerable computational resources in comparison with the standard FEM.Although the SGBEM yields a symmetric system of linear equations, the coefficient matrix is fully dense and each of its entries must be computed by means of a double surface integration.
In the past two decades, various investigators have seriously attempted to develop efficient and accurate numerical procedures for analysis of elasticity and fracture problems by exploiting positive features of both the BEM and the FEM.The fundamental idea is to decompose the entire domain into two regions and then apply the BEM to model a linearly elastic region with small surface-to-volume ratio and possibly containing the displacement discontinuities e.g., cracks and dislocations and the FEM to model the remaining majority of the domain possibly exhibiting complex behavior e.g., material nonlinearity and nonhomogeneous data .The primary objective is to compromise between the requirement of computational resources and accuracy of predicted results.Within the context of linear elasticity, there have been several investigations directed towards the coupling of the conventional BEM and the standard FEM e.g., 24-26 and the coupling of the strongly singular SGBEM and the standard FEM e.g., 27-30 .It should be emphasized that the former type of coupling procedure generally destroys the desirable symmetric feature of the entire system of linear algebraic equations whereas the latter type requires special numerical treatment of strongly and hyper singular integrals e.g., 31, 32 .Extensive review of various types of coupling between boundary integral equation methods and finite element techniques can be found in 33 .Among those existing techniques, a particular symmetric coupling strategy between the weakly singular SGBEM and the standard FEM has been found computationally efficient and has recently become an attractive alternative for performing comprehensive stress and fracture analysis.This is due primarily to i the symmetry feature of the SGBEM that leads to the symmetric coupling formulation and ii the weakly singular nature of all involved boundary integrals requiring simpler theoretical and numerical treatment in comparison with strongly singular and hypersingular integrals.Xiao 16 first presented such coupling formulation for cracks in isotropic, linearly elastic finite bodies; more precisely, a pair of weakly singular, weak-form displacement and traction integral equations was utilized along with the principle of virtual work and the proper enforcement of continuity conditions on the interface to establish the symmetric coupling formulation.Later, Frangi and Novati 34 successfully implemented Xiao's formulation to analyze cracked bodies subjected to pure traction boundary conditions.Besides its accuracy and robustness, the technique was still restricted to the conforming discretization of the interface between the two regions.Springhetti et al. 35 relaxed such limitation by allowing the weak enforcement of continuity across the interface and also generalized the technique to treat both potential and elastostatic problems.Nevertheless, their main focus is on uncracked bodies made of linearly isotropic materials.Recently, Rungamornrat and Mear 36 extended the work of Xiao 16 to enable the treatment of both material anisotropy and nonmatching interface.While this particular coupling scheme has been well-established for decades, on the basis of an extensive literature survey, applications of this technique to model a problem of an infinite space containing isolated cracks and localized complex zones have not been recognized.
In this paper, a numerical procedure based on the symmetric coupling between the weakly singular SGBEM and the standard FEM is implemented to perform three-dimensional stress analysis of an infinite medium containing displacement discontinuities and localized complex zones.Vast features of the FEM are exploited to allow the treatment of very general localized zones, for instance, those exhibiting material nonlinearity, material nonuniformity, and other types of complexity.The weakly singular SGBEM is utilized to readily and efficiently model the remaining unbounded region.A pair of weakly singular boundary integral equations proposed by Rungamornrat and Mear 22 is employed as a basis for the development of SGBEM, and this, therefore, allows the treatment of the unbounded region that is made of an anisotropic linearly elastic material and contains cracks.It is worth emphasizing that while the present study is closely related to the work of Frangi and Novati 34 , Springhetti et al. 35 , and Rungamornrat and Mear 36 , the proposed technique offers additional crucial capabilities to treat an infinite domain, material nonlinearity in localized zones, and general material anisotropy.Following sections of this paper present basic equations and the coupling formulation, essential components for numerical implementations, numerical results and discussions, and conclusions and useful remarks.

Formulation
Consider a three-dimensional infinite medium, denoted by Ω, containing a crack and a localized complex zone as shown schematically in Figure 1 a .The crack is represented by two geometrical coincident surfaces S C and S − C with their unit outward normal being denoted by n and n − , respectively, and the localized complex zone is denoted by Ω L .In the present study, the medium is assumed to be free of a body force and loading on its remote boundary, and both surfaces of the crack are subjected to prescribed self-equilibrated traction defined by C .Now, let us introduce an imaginary surface S I to decompose the body Ω into two subdomains, an unbounded "BEM-region" denoted by Ω B and a finite "FEM-region" denoted by Ω F , as indicated in Figure 1 b .The surface S I is selected such that the localized complex zone and the crack are embedded entirely in the FEM-region and in the BEM-region, respectively i.e., S C ∪ S − C ⊂ Ω B and Ω L ⊂ Ω F , and, in addition, the BEM-region must be linearly elastic.To clearly demonstrate the role of the interface between the two subregions in the formulation presented below, we define {S BI , t BI , u BI } and {S FI , t FI , u FI } as the interface, the unknown traction, and the unknown displacement on the interface of the BEM-region Ω B and the FEM-region Ω F , respectively.It is important to emphasize that the interfaces S BI and S FI are in fact identical to the imaginary surface S I .While the formulation is presented, for brevity, only for a domain containing a single crack and a single localized complex zone, it can readily be extended to treat multiple cracks and multiple complex zones; in such particular case, multiple FEM-regions are admissible.

Governing Equations for Ω B
The total boundary of the BEM-region Ω B , denoted by S B , consists of the reduced crack surface S BC ≡ S C on which the traction is fully prescribed and the interface S BI where neither the traction nor the displacement is known a priori.Note again that the subscript "B" is added only to emphasize that those surfaces are associated with the BEM-region.To form a set of governing integral equations for this region, a pair of weakly singular, weak-form displacement and traction boundary integral equations developed by Rungamornrat and Mear 22 is employed.These two integral equations were derived from standard boundary integral relations for the displacement and stress along with a systematic regularization technique.The final form of completely regularized integral equations is well suited for establishing the symmetric formulation for the weakly singular SGBEM.Such pair of integral equations is given here, for convenience in further reference, by where t k and v k are sufficiently smooth test functions; D m • n i ε ism ∂ • /∂ξ s is a surface differential operator with ε ism denoting a standard alternating symbol; v j ξ u j ξ for ξ ∈ S BI and Δu j ξ for ξ ∈ S BC with Δu j ξ u j ξ − u − j ξ denoting the jump in Localized complex zone where δ pj is a standard Kronecker-delta symbol; E ijkl are elastic moduli; A tkoe mjdn and K ik jl ξ − y are defined by in which r ξ − y, r r , z, z jk z i E ijkl z l and the closed contour integral is to be evaluated over a unit circle z 1 on a plane defined by z • r 0. It is evident that the kernel H p ij ξ − y is given in an explicit form independent of material properties and the kernels Towards obtaining a system of governing integral equations for the BEM-region Ω B , the weak-form boundary integral equation for the traction 2.2 is applied directly to the crack surface S BC with the test function being chosen such that v 0 on S BI and to the interface S BI with the test function being chosen such that v 0 on S BC , and the weakform boundary integral equation for the displacement 2.1 is applied only to the interface S BI .A final set of three integral equations is given concisely by where { u BI , t BI } are test functions defined on the interface S BI and all involved bilinear integral operators are defined, with subscripts P, Q ∈ {I, C} being introduced to clearly indicate the surface of integration, by 12 13 It should be noted that the linear operator A PQ X, Y is in a symmetric form satisfying the condition A PQ X, Y A QP Y, X , and, as a consequence, it renders the left hand side of the system 2.9 being in a symmetric form.Although such symmetric formulation can readily be obtained, the right-hand side of 2.9 still contains the unknown traction on the interface t BI .The treatment of a term F BI u BI , t BI will be addressed once the formulation for the FEMregion Ω F is established.

Governing Equations for Ω F
Let us consider, next, the FEM-region Ω F .For generality, the entire boundary of this particular region can be decomposed into two surfaces: the interface S FI on which both the traction t FI and the displacement u FI are unknown a priori and the surface S FT on which the traction t FT is fully prescribed.The existence of the surface S FT is apparent for the case that the FEM-region contains embedded holes or voids.It is also important to emphasize that, in the development of a key governing equation for Ω F , the traction t FI is treated, in a fashion different from that for the BEM-region, as unknown data instead of the primary unknown variable.In addition, to be capable of modeling a complex localized zone embedded within the FEM-region, a constitutive model governing the material behavior utilized in the present study is assumed to be sufficiently general allowing the treatment of material nonlinearity, anisotropy, and inhomogeneity.The treatment of such complex material models has been extensively investigated and well established within the context of nonlinear finite element methods e.g., 6, 37, 38 , and those standard procedures also apply to the current implementation and will not be presented for brevity.Here, we only outline the key governing equation for the FEM-region and certain unknowns and necessary data connected to those for the BEM-region.
Following standard formulation of the finite element technique, the weak-form equation governing the FEM-region can readily be obtained via the principle of virtual work 6-8 and the final equation can be expressed in a concise form by where σ is a stress tensor; u is a suitably well-behaved test function defined over the domain Ω F ; u FI and u FT are the restriction of u on the interface S FI and boundary S FT , respectively; the integral operators are defined, with subscripts P ∈ {I, T}, by 16 in which ε ij y denotes the virtual strain tensor defined by ε ij y ∂ u i /∂y j ∂ u j /∂y i /2.Note again that a function form of the stress tensor in terms of the primary unknown depends primarily on a constitutive model employed.For a special case of the FEM-region being made of a homogeneous, linearly elastic material, the stress tensor can be expressed directly and explicitly in terms of elastic constants E ijkl and the strain tensor ε i.e., σ ij E ijkl ε kl , and, within the context of an infinitesimal deformation theory i.e., ε ij y ∂u i /∂y j ∂u j /∂y i /2 , the integral operator K FF can be expressed directly in terms of the displacement u as It should be remarked that the factor of one half in the definition 2.17 has been introduced for convenience to cast this term in a form analogous to that for F BP given by 2.13 , and this, as a result, leads to the factor of two appearing on the right-hand side of 2.15 .It is also worth noting that the first term on the right-hand side of 2.15 still contains the unknown traction on the interface t FI .

Governing Equations for Ω
A set of governing equations of the entire domain Ω can directly be obtained by combining a set of weakly singular, weak-form boundary integral equations 2.9 and the virtual work equation 2.15 .In particular, the last equations of 2.9 and 2.15 are properly combined, and this finally leads to where E is given by From the continuity of the traction and displacement across the interface of the BEM-region and FEM-region i.e., t BI y t FI y 0 and u BI y u FI y for all y ∈ S I S BI S FI , the test functions u BI and u are chosen such that u BI y u FI y for all y ∈ S I S BI S FI and, as a direct consequence, E identically vanishes.It is therefore evident that the right-hand side of 2.19 involves only prescribed boundary data, and, in addition, if the integral operator K FF possesses a symmetric form, 2.19 constitutes a symmetric formulation for the boundary value problem currently treated.

Numerical Implementation
This section briefly summarizes numerical procedures adopted to construct approximate solutions of a set of governing equations 2.19 and to postprocess certain quantities of interest.The discretization of the BEM-region and the FEM-region is first discussed.Then, components essential for numerical evaluation of weakly singular and nearly singular double-surface integrals, evaluations of kernels, and determination of general mixed-mode stress intensity factors are addressed.Finally, the key strategy for establishing the coupling between the in-house weakly singular SGBEM code and the reliable commercial finite element package is discussed.

Discretization
Standard Galerkin strategy is adopted to construct an approximate solution of the governing equation 2.19 .For the BEM-region Ω B , only the crack surface S BC and the interface S BI need to be discretized.In such discretization, standard isoparametric C 0 elements e.g., 8-node quadrilateral and 6-node triangular elements are employed throughout except along the crack front where special 9-node crack-tip elements are employed to accurately capture the asymptotic field near the singularity zone.Shape functions of these special elements are properly enriched by square root functions, and, in addition, extra degrees of freedom are introduced along the element boundary adjacent to the crack front to directly represent the gradient of the relative crack-face displacement 18, 21, 23 .These positive features also enable the calculation of the mixed-mode stress intensity factors i.e., mode-I, mode-II, and mode-III stress intensity factors in an accurate and efficient manner with Mathematical Problems in Engineering use of reasonably coarse meshes.For the FEM-region Ω F , standard three-dimensional, isoparametric C 0 elements e.g., ten-node tetrahedral elements, fifteen-node prism elements, and twenty-node brick elements are utilized throughout in the domain discretization.
It is important to note that the BEM-region and the FEM-region are discretized such that meshes on the interface of the two regions conform i.e., the two discretized interfaces are geometrically identical .A simple means to generate those conforming interfaces is to mesh the FEM-region first and then use its surface mesh as the interface mesh of the BEMregion.With this strategy, all nodal points on both discretized interfaces are essentially coincident.The key advantage of using conforming meshes is that the strong continuity of the displacement, the traction, and the test functions across the interface can be enforced exactly, and, as a result, the condition E 0 is also satisfied in the discretization level.It should be emphasized also that nodes on the interface of the BEM-region contain six degrees of freedom i.e., three displacement degrees of freedom and three traction degrees of freedom while nodes on the FEM-region contain only three displacement degrees of freedom.

Numerical Integration
For the FEM-region, all integrals arising from the discretization of the weak-form equation contain only regular integrands, and, as a result, they can be efficiently and accurately integrated by standard Gaussian quadrature.On the contrary, numerical evaluation of integrals arising from the discretization of the BEM-region is nontrivial since it involves the treatment of three types of double surface integrals i.e., regular integrals, weakly singular integrals, and nearly singular integrals .The regular double surface integral arises when it involves a pair of remote outer and inner elements i.e., the distance between any source and field points is relatively large when compared to the size of the two elements .This renders its integrand nonsingular and well-behaved and, as a result, allows the integral to be accurately and efficiently integrated by standard Gaussian quadrature.
The weakly singular double surface integral arises when the outer surface of integration is the same as the inner surface.For this particular case, the source and field points can be identical and this renders the integrand singular of order 1/r.While the integral of this type exists in an ordinary sense, it was pointed out by Xiao 16 that the numerical integration by Gaussian quadrature becomes computationally inefficient and such inaccurate evaluation can significantly pollute the quality of approximate solutions.To circumvent this situation, a series of transformations such as a well-known triangular polar transformation and a logarithmic transformation is applied first both to remove the singularity and to regularize the rapid variation of the integrand.The final integral contains a nonsingular integrand well suited to be integrated by Gaussian quadrature.Details of this numerical quadrature can be found in 16, 39, 40 .
The most challenging task is to compute nearly singular integrals involving relatively close or adjacent inner and outer elements.Although the integrand is not singular, it exhibits rapid variation in the zone where both source and field points are nearly identical.Such complex behavior of the integrand was found very difficult and inefficient to be treated by standard Gaussian quadrature 16 .To improve the accuracy of such quadrature, the triangular polar transformation is applied first and then a series of logarithmic transformations is adopted for both radial and angular directions to further regularize the rapid variation integrand.The resulting integral was found well-suited for being integrated by standard Gaussian quadrature 16, 40-42 .

Evaluation of Kernels
To further reduce the computational cost required to form the coefficient matrix contributed from the BEM-region, all involved kernels n i ξ H As a result, this task can readily be achieved via a standard procedure.For the last three kernels, the computational cost is significantly different for isotropic and anisotropic materials.For isotropic materials, such kernels only involve elementary functions and can therefore be evaluated in a straightforward fashion.On the contrary, the kernels U p i ξ − y , G p mj ξ − y , and C tk mj ξ − y for general anisotropy are expressed in terms of a line integral over a unit circle see 2.4 -2.6 , and 2.8 .Direct evaluation of such line integral for every pair of points ξ, y arising from the numerical integration is obviously computationally expensive.To avoid this massive computation, a well-known interpolation technique e.g., 21, 23, 36 is adopted to approximate values of those kernels.Specifically, the interpolant of each kernel is formed based on a two-dimensional grid using standard quadratic shape functions.Values of kernels at all grid points are obtained by performing direct numerical integration of the line integral 2.8 via Gaussian quadrature and then using the relations 2.4 -2.6 .The accuracy of such approximation can readily be controlled by refining the interpolation grid.

Determination of Stress Intensity Factors
Stress intensity factors play an important role in linear elastic fracture mechanics in the prediction of crack growth initiation and propagation direction and also in the fatigue-life assessment.This fracture data provides a complete measure of the dominant behavior of the stress field in a local region surrounding the crack front.To obtain highly accurate stress intensity factors, we supply the developed coupling technique with two crucial components, one associated with the use of special crack-tip elements to accurately capture the neartip field and the other corresponding to the use of an explicit formula to extract such fracture data.The latter feature is a direct consequence of the extra degrees of freedom being introduced along the crack front to represent the gradient of relative crack-face displacement.Once a discretized system of algebraic equations is solved, nodal quantities along the crack front are extracted and then postprocessed to obtain the stress intensity factors.An explicit expression for the mixed-mode stress intensity factors in terms of nodal data along the crack front, local geometry of the crack front, and material properties can be found in the work of Li et al. 18 for cracks in isotropic media and Rungamornrat and Mear 23, 36 for cracks in general anisotropic media.In the current investigation, both formulae are implemented.

Coupling of SGBEM and Commercial FE Package
To further enhance the modeling capability, the weakly singular SGBEM can be coupled with a reliable commercial finite element code that supports user-defined subroutines.The key objective is to exploit available vast features of such FE package e.g., mesh generation, userdefined elements, powerful linear and nonlinear solvers, and various material models, etc. to treat a complex, localized FEM-region and utilize the SGBEM in-house code to supply information associated with the majority of the domain that is unbounded and possibly contains isolated discontinuities.
In the coupling procedure, the governing equation for the BEM-region is first discretized into a system of linear algebraic equations.The corresponding coefficient matrix and the vector involving the prescribed data are constructed using the in-house code, and they can be viewed as a stiffness matrix and a load vector of a "super element" containing all degrees of freedom of the BEM-region.This piece of information is then imported into the commercial FE package via a user-defined subroutine channel and then assembled with element stiffness matrices contributed from the discretized FEM-region.Since meshes of both interfaces one associated with the BEM-region and the other corresponding to the FEM-region are conforming, the assembly procedure can readily be achieved by using a proper numbering strategy.Specifically, nodes on the interface of the BEM-region are named identical to nodes on the interface of the FEM-region associated with the same displacement degrees of freedom .It is important to emphasize that all interface nodes of the BEM-region possess six degrees of freedom i.e., three displacement degrees of freedom and three traction degrees of freedom but there are only three displacement degrees of freedom per interface node of the FEM-region.To overcome such situation, each interface node of the BEM-region is fictitiously treated as double nodes where the first node is chosen to represent the displacement degrees of freedom and is numbered in the same way as its coincident interface node of the FEM-region whereas the second node is chosen with different name to represent the traction degrees of freedom.With this particular scheme, the assembling procedure follows naturally that for a standard finite element technique.
Once the coupling analysis is complete, nodal quantities associated with the BEMregion are extracted from the output file generated by the FE package and then postprocessed for quantities of interest.For instance, the displacement and stress within the BEM-region can readily be computed from the standard displacement and stress boundary integral relations 17, 22 , and the stress intensity factors can be calculated using an explicit expression proposed in 18, 23 .

Numerical Results and Discussion
As a means to verify both the formulation and the numerical implementations, we first carry out numerical experiments on boundary value problems in which the analytical solution exists.In the analysis, a series of meshes is adopted in order to investigate both the convergence and accuracy of the numerical solutions.Once the method is tested, it is then applied to solve more complex boundary value problems in order to demonstrate its capability and robustness.For brevity of the presentation, a selected set of results are reported and discussed as follows.

Isolated Spherical Void under Uniform Pressure
Consider an isolated spherical void of radius a embedded in a three-dimensional infinite medium as shown schematically in Figure 2 a .The void is subjected to uniform pressure σ 0 .In the analysis, two constitutive models are investigated: one associated with an isotropic, linearly elastic material with Young's modulus E and Poisson ratio ν 0.3 and the other chosen to represent an isotropic hardening material obeying J 2 -flow theory of plasticity 43 .

Mesh-1
Mesh-2 Mesh-3 Figure 3: Three meshes adopted in the analysis for FEM-region; meshes for BEM-region are identical to the interface mesh of FEM-region.
For the latter material, the uniaxial stress-strain relation is assumed in a bilinear form with E 1 and E 2 denoting the modulus in the elastic regime and the modulus of the hardening zone, respectively, and σ y and ε y denoting the initial yielding stress and its corresponding strain, respectively.
To test the coupling technique, we first decompose the body into two regions by a fictitious spherical surface of radius 5a and centered at the origin as shown by a dashed line in Figure 2 b .It is important to remark that such a surface must be chosen relatively large compared to the void to ensure that the inelastic zone that may exist for the second constitutive model is fully contained in the FEM-region.In the experiments, three different meshes are adopted as shown in Figure 3.Although meshes for the BEM-region are not shown, they can simply obtain from the interface meshes of the FEM-region.As clearly illustrated in the figure, mesh-1, mesh-2, and mesh-3 consist of 12, 32, 144 boundary elements and 24, 128, 1152 finite elements, respectively.

Results for Isotropic Linearly Elastic Material
For linear elasticity, this particular boundary value problem admits the closed form solution e.g., 44 .Since the problem is spherically symmetric, only the radial displacement u r and the normal stress components {σ rr , σ θθ , σ φφ } are nonzero and they are given explicitly by note that these quantities are referred to a standard spherical coordinate system {r, θ, φ} with its origin located at the center of the void

4.2
This analytical solution is employed as a means to validate the proposed formulation and the numerical implementation.Numerical solutions for the radial displacement obtained from the three meshes are reported and compared with the exact solution in Figure 4.As evident from this set of results, the radial displacements obtained from the mesh-2 and the mesh-3 are highly accurate with only slight difference from the exact solution while that obtained from the mesh-1 is reasonably accurate except in the region very near the surface of the void.The discrepancy of solutions observed in the mesh-1 is due to that the level of refinement is too coarse to accurately capture the geometry and responses in the local region near the surface of the void.We further investigate the quality of numerical solutions for stresses.Since all nonzero stress components are related by 4.2 , only results for the radial stress component are reported.Figure 5 shows the normalized radial stress obtained from the three meshes and the exact solution versus the normalized radial coordinate.It is observed that the mesh-3 yields results that are almost indistinguishable from the exact solution, whereas the mesh-1 and mesh-2 give accurate results for relatively large r, and the level of accuracy decreases as the distance r approaches a.It is noted by passing that the degeneracy of the accuracy in computing stress is common in a standard, displacement-based, finite element technique.
To demonstrate the important role of the SGBEM in the treatment of an unbounded part of the domain instead of truncating the body as practically employed in the finite element modeling, we perform another FE analysis of the FEM-region alone without coupling with the BEM-region but imposing zero displacement condition at its interface.The radial displacement and the radial stress obtained for this particular case using the mesh-3 are reported along with the exact solution and those obtained from the coupling technique in Figures 6 and 7, respectively.As evident from these results, numerical solutions obtained from the FEM with a domain truncation strategy deviate from the exact solution as approaching the truncation surface while the proposed technique yields almost identical results to the exact solution.The concept of domain truncation to obtain a finite body is simple, but it still remains to choose a proper truncation surface and boundary conditions to be imposed on that surface to mimic the original boundary value problem.This coupling technique provides an alternative to treat the whole domain without any truncation and difficulty to treat the remote boundary.

Isotropic Hardening Material
For this particular case, we focus attention to the material with no hardening modulus i.e., E 2 0 since the corresponding boundary value problem admits the closed form solution.For a sufficiently high applied pressure σ 0 , a layer close to the boundary of the void becomes inelastic and the size of such inelastic zone measured by the radius r 0 becomes larger as   σ 0 increases.By incorporating J 2 -flow theory of plasticity and spherical symmetry, the radial displacement and the radial stress can be obtained exactly as given below: where the Poisson ratio ν is taken to be 0.3 and r 0 ae σ 0 /2σ y − 1/3 is the radius of an inelastic zone.
In the analysis, the pressure σ 0 1.625σ y is chosen to ensure that the medium contains an inelastic zone; in fact, this selected applied pressure corresponds to r 0 1.615a.Numerical results obtained from the current technique are reported along with the exact solution in Figure 8 for the normalized radial displacement and in Figure 9 for the normalized radial stress.It can be concluded from computed solutions that they finally converge to the exact solution as the mesh is refined.In particular, results obtained from the mesh-3 are nearly indistinguishable from the benchmark solution.It should be pointed out that results obtained from the same level of mesh refinement for this particular case are less accurate than those x 1 obtained for the linear elastic case.This is due to the complexity posed by the presence of an inelastic zone near the surface of the void, and, in order to capture this behavior accurately, it requires sufficiently fine meshes.

Isolated Penny-Shaped Crack in Infinite Medium
Consider, next, a penny-shaped crack of radius a which is embedded in a linearly elastic, infinite medium as shown schematically in Figure 10 a .The body is made of either an  isotropic material with Poisson's ratio ν 0.3 or zinc and graphite-reinforced composite.The last two materials are transversely isotropic with the axis of material symmetry directing along the x 3 -axis, and their elastic constants are given in Table 1.The crack is subjected to two types of traction boundary conditions: the uniform normal traction σ 0 i.e., t 1 t 2 0, t 3 σ 0 as shown in Figure 10 b and the uniform shear traction τ 0 along the x 1 -axis i.e., t 1 τ 0 , t 2 t 3 0 as shown in Figure 10 c .The first loading condition gives rise to a pure opening-mode problem with the mode-I stress intensity factor along the crack front being constant and independent of material properties, while the second loading condition yields nonzero mode-II and mode-III stress intensity factors that vary along the crack front.The analytical solutions for both cases can be found in the work of Fabrikant 4 .As a means to verify the coupling formulation and implementation, we choose the FEM-region to be a cube of dimensions 2a × 2a × 2a centered at 0, 0, 2a as illustrated in Figure 11 a .In the analysis, we generate three meshes for both the crack surface and the FEM-region as shown in Figure 11 b .
For the first loading condition, numerical solutions for the mode-I stress intensity factor normalized by the exact solution are reported in Table 2 for all three materials.Clearly from these results, the current technique yields highly accurate stress intensity factors with error less than 1.6%, 0.6%, and 0.1% for mesh-1, mesh-2, and mesh-3, respectively.The weak dependence of numerical solutions on the level of mesh refinement is due mainly to the use of special crack-tip elements to model the near-tip field and directly capture the gradient of relative crack-face displacement along the crack front.Relatively coarse meshes can therefore be employed in the analysis to obtain sufficiently accurate stress intensity factors.
For the second loading condition, the normalized mode-II and mode-III stress intensity factors K 2 and K 3 are shown in Figure 12 for isotropic material, zinc and graphitereinforced composite.Based on this set of results, it can be concluded again that numerical solutions obtained from the three meshes are in excellent agreement with the exact solution; in particular, a coarse mesh also yields results of high accuracy.It should also be remarked that for this particular loading condition, the material anisotropy plays a significant role on values of the mixed-mode stress intensity factors.

Infinite Medium Containing Both Penny-Shaped Crack and Spherical Void
As a final example, we choose to test the proposed technique by solving a more complex boundary value problem in order to demonstrate its capability.Let us consider an infinite medium containing a spherical void of radius a and a penny-shaped crack of the same radius as shown schematically in Figure 13.The medium is subjected to uniform pressure σ 0 on the surface of the void, whereas the entire surface of the crack is traction-free.In the analysis, two constitutive models are investigated: one associated with an isotropic, linearly elastic material with Young's modulus E and Poisson ratio ν 0.3 and the other corresponding to an isotropic hardening material with the bilinear uniaxial stress-strain relation similar to that employed in Section 4.1.The primary quantity to be sought is the mode-I stress intensity factor along the crack front induced by the application of the pressure to the void.In addition, influence of an inelastic zone induced in the high-load-intensity region on such fracture data is also of interest.
In the modeling, we first decompose the medium into the FEM-region and the BEMregion using a fictitious spherical surface of radius 4a centered at the same location as that of the void as shown in Figure 14 a .Three meshes are adopted in numerical experiments as shown in Figure 14 b .In particular, the FEM-region, the interface, and the crack surface consist of {24, 12, 8}, {128, 32, 16}, and {1024, 128, 64} elements for mesh-1, mesh-2, and mesh-3, respectively.It should be noted also that the mesh-1 is obviously very coarse; in particular, only eight elements are utilized to discretize the entire crack surface and only four relatively large crack-tip elements are used along the crack front.First, the analysis is carried out for the elastic material with Poisson ratio ν 0.3, and the computed mode-I stress intensity factors are normalized and then reported as a function of angular position along the crack front for all three meshes in Figure 15.This set of results implies that the obtained numerical solutions exhibit good convergence; in particular, results obtained from the mesh-2 and mesh-3 are of comparable quality while results obtained from the mesh-1 still deviate from the converged solution.As confirmed by this convergence study, only the mesh-3 is used to generate other sets of useful results.
Next, we consider a medium made of an isotropic hardening material.In the analysis, we choose the modulus E 1 E and Poisson ratio ν 0.3 for the linear regime and choose either E 2 E/3 or E 2 0 for the hardening regime.With this set of material parameters, the behavior in the linear regime for a small level of applied pressure is identical to that obtained in the previous case.To investigate the influence of the inelastic zone induced near Mathematical Problems in Engineering the surface of the void on the stress intensity factor along the crack front, we carry out various experiments by varying the applied pressure σ 0 .The distribution of the stress intensity factor along the crack front is reported in Figure 16 for a hardening material with E 1 E and E 2 E/3 under five levels of the applied pressure σ 0 ∈ {0.25σ y , 1.00σ y , 1.25σ y , 1.50σ y , 1.75σ y }.The body is entirely elastic at σ 0 0.25σ y , slightly passes the initial yielding at σ 0 1.00σ y , and possesses a larger inelastic zone as the pressure increases further.It is obvious from Figure 16  that the presence of an elastic zone significantly alters the normalized values of the stress intensity factor from the linear elastic solution and such discrepancy becomes more apparent as the level of applied pressure increases.The localized inelastic zone acts as a stress riser, that is, it produces the stress field of higher intensity around the crack, and this therefore yields the higher normalized stress intensity factor when compared with the linear elastic case.Figure 17 shows an additional plot between the maximum normalized stress intensity factors versus the normalized applied pressure for both an isotropic linearly elastic material and two isotropic hardening materials associated with E 2 0 and E 2 E/3 .Results for both materials are identical for a low level of the applied pressure since the entire body is still elastic , and, for a higher level of the applied pressure, the maximum stress intensity factor for the case of the hardening material is significantly larger than that for the linear elastic material.In addition, such discrepancy tends to increase as the hardening modulus decreases.

Conclusions and Remarks
The coupling procedure between a standard finite element method FEM and a weakly singular, symmetric Galerkin boundary element method SGBEM has been successfully established for stress analysis of a three-dimensional infinite medium.The proposed technique has exploited the positive features of both the FEM and the SGBEM to enhance the modeling capability.The vast and very general features of the FEM have been employed as a basis to treat a localized region that may embed a zone exhibiting complex behavior, whereas the SGBEM has been used specifically to model the majority of the medium that is unbounded and possibly contains the surface of displacement discontinuities such as cracks and dislocations.The coupling formulation has been based primarily on the domain decomposition technique along with the proper enforcement of continuity of the displacement and traction on the interface of the two regions one modeled by the SGBEM and the other by the FEM to form the coupling equations.For the FEM subdomain, the key formulation follows directly the well-known principle of virtual work, whereas, for the SGBEM subdomain, the governing equation is formulated based on a pair of weakly singular, weak-form boundary integral equations for the displacement and traction.The advantage of using the weakly singular integral equations is associated with the permission to apply a space of continuous interpolation functions in the discretization of primary unknowns on the SGBEM subdomain.
In the numerical implementation, various aspects have been considered in order to enhance the accuracy and computational efficiency of the coupling technique.For instance, special crack-tip elements have been employed to better approximate the near-tip field.Shape functions of these special elements have properly been enriched by a square root function such that the resulting interpolation functions can capture the relative crack-face displacement with sufficiently high level of accuracy.As a direct consequence, it allows relatively large crack-tip elements to be employed along the crack front while still yielding very accurate stress intensity factors.Another important consideration is the use of an interpolation strategy to approximate values of kernels for generally anisotropic materials; this substantially reduces the computational cost associated with the direct evaluation of the line integral.Finally, special numerical quadratures have been adopted to efficiently evaluate both the weakly singular and nearly singular double surface integrals.To demonstrate and gain an insight into the coupling strategy, the formulation has been implemented first in terms of an in-house computer code for linear elasticity boundary value problems.Subsequently, the weakly singular SGBEM has successfully been coupled with a reliable commercial finite element package in order to exploit its rich features to model more complex localized region such as inelastic zones and inhomogeneities.As indicated by results from extensive numerical experiments, the current technique has been found promising and, in particular, numerical solutions exhibit good convergence and weak dependence on the level of mesh refinement.
As a final remark, while the developed technique is still restricted to an infinite domain and to conforming interfaces, it offers insight into the SGBEM-FEM coupling strategy in terms of the formulation, the implementation procedure, and its performance.This coupling strategy can directly be generalized to solve more practical boundary value problems involving a half space, for example, cracks and localized complex zone near the free surface.Another crucial extension is to enhance the feature of the current technique by using the weak enforcement of continuity across the interface.This will supply the flexibility of mesh generation.

Figure 1 :
Figure 1: a Schematic of three-dimensional infinite medium containing crack and localized complex zone and b schematic of BEM-region Ω B and FEM-region Ω F .
p mj ξ −y , C tk mj ξ −y are singular only at ξ y of O 1/r see details in Xiao 16 for discussion of the singularity behavior of the kernel H p ij ξ − y .It should be remarked that for isotropic materials, the kernels U p i ξ − y , G p mj ξ − y , C tk mj ξ − y possess an explicit form in terms of elementary functions see 18, 22 .
p ij ξ − y , n i y H p ij ξ − y , U p i ξ − y , G p mj ξ − y , C tk mj ξ − y must be evaluated in an efficient manner for any pair of source and field points {ξ, y}.For the first two kernels n i ξ H p ij ξ − y and n i y H p ij ξ − y , they only involve the calculation of a unit normal vector n and the elementary function H p ij .

Figure 2 :
Figure 2: a Schematic of three-dimensional infinite medium containing spherical void and b schematic of BEM-region and FEM-region.

Figure 6 :
Figure 6: Normalized radial displacement versus normalized radial coordinate for isotropic, linearly elastic material with ν 0.3.Results are obtained from mesh-3 for both the coupling technique and the FEM with domain truncation.

Figure 7 :
Figure 7: Normalized radial stress versus normalized radial coordinate for isotropic, linearly elastic material with ν 0.3.Results are obtained from mesh-3 for both the coupling technique and the FEM with domain truncation.

Figure 8 :
Figure 8: Normalized radial displacement versus normalized radial coordinate for isotropic hardening material with E 2 0.

Figure 9 :
Figure 9: Normalized radial stress versus normalized radial coordinate for isotropic hardening material with E 2 0.

Figure 10 :
Figure 10: a Schematic of infinite medium containing penny-shaped crack, b crack under uniform normal traction σ 0 , and c crack under uniform shear traction τ 0 .

Figure 11 :
Figure 11: a Schematic of selected FEM-region and the remaining BEM-region and b three meshes adopted in the analysis.

Figure 13 :
Figure 13: Schematic of infinite medium containing spherical void of radius a and penny-shaped crack of radius a and subjected to uniform pressure at surface of void.

Figure 14 :Figure 15 :
Figure 14: a Decomposition of domain into BEM-region and FEM-region by a fictitious spherical surface of radius 4a and b three meshes adopted in analysis.

5 ]Figure 16 :
Figure 16: Normalized mode-I stress intensity factor of penny-shaped crack embedded within infinite medium containing spherical void under uniform pressure.Results are reported for isotropic hardening material with E 1 E and E 2 E/3.

E 1 = 5 ]Figure 17 :
Figure 17: Maximum normalized mode-I stress intensity factor versus applied pressure at surface of void.Results are reported for isotropic linearly elastic material with ν 0.3 and two isotropic hardening materials.

Table 1 :
Nonzero elastic constants for zinc and graphite-reinforced composite where axis of material symmetry is taken to direct along the x 3 -coordinate direction .

Table 2 :
Normalized mode-I stress intensity factor for isolated penny-shaped crack subjected to uniform normal traction.