A Boundary-Element Analysis of Crack Problems in Multilayered Elastic Media: A Review

: Crack problems in multilayered elastic media have attracted extensive attention for years due to their wide applications in both a theoretical analysis and practical industry. The boundary element method (BEM) is widely chosen among various numerical methods to solve the crack problems. Compared to other numerical methods, such as the phase ﬁeld method (PFM) or the ﬁnite elem ent method (FEM), the BEM ensures satisfying accuracy, broad applicability, and satisfactory eﬃciency. Therefore, this paper reviews the state -of-the-art progress in a boundary-element analysis of the crack problems in multilayered elastic media by concentrating on implementations of the two branches of the BEM: the displacement discontinuity method (DDM) and the direct method (DM). The review shows limitation of the DDM in applicability at ﬁrst and subsequently reveals the inap-plicability of the conventional DM for the crack problems. After that, the review outlines a pre-treatment that makes the DM applicable for the crack problems and presents a DM-based method that solves the crack problems more eﬃciently than the conventional DM but still more slow ly than the DDM. Then, the review highlights a method that combines the DDM and the DM so that it shares both the eﬃciency of the DDM and broad applicability of the DM after the pre -treatment, making it a promising candidate for an analysis of the crack problems. In addition, the paper presents numerical examples to demonstrate an even faster approximation with the combined method for a thin layer, which is one of the challenges for hydraulic-fracturing simulation. Finally, the review concludes with a comprehensive summary and an outlook for future study.

The phase field method (PFM) has become a popular numerical method for the analysis of the crack problems in recent years, since its implementation is free of crack tracking algorithms, enabling relatively easy simulation of complicated crack patterns, such as nucleation, branching, and coalescence [25][26][27][28].The PFM, however, requires highly fine discretization due to the small length scale and becomes staggeringly computational-resource intensive as a result [29][30][31], severely limiting its application in actual industry.In addition, the PFM is usually incorporated with the finite element method (FEM) during implementation.Therefore, some computational issues of the FEM, such as sensitivity to geometry of the domain and size of discretization, also affect the PFM [31][32][33].
Among various numerical methods for the analysis of the crack problems to compute crack opening displacement, the boundary element method (BEM) is widely chosen, because it only discretizes boundaries of a domain of interest and consequently enables a more-efficient computation and a more-convenient implementation [34][35][36][37], while most of the other methods, such as the PFM or FEM, require discretization of the whole domain [31].Furthermore, the BEM exploits existed analytical fundamental solutions, ensuring a potentially more accurate computation [38][39][40][41].Thus, it is essential to determine analytical fundamental solutions appropriately for precise determination of boundary-element equations and efficient implementation to subsequently solve the crack problems.
The governing equation of the multilayered media is well known as a system of partial differential equations (PDEs) from the theory of elasticity [37,41].Therefore, it is natural to find the analytical fundamental solutions by reducing the PDEs to ordinary differential equations (ODEs) for simpler calculation [42].As a result, a variant of the BEM called the Fourier transform (FT) method that utilizes FT-type convolution has been proposed by many authors to reduce the PDEs of the system [19,40,[42][43][44][45].When the solutions of stress and displacement to the ODEs are solved in the FT field, a following inverse FT-type convolution leads to analytical fundamental solutions in the spatial field.Unfortunately, the convolution that has to be conducted twice for the method often results in improper integrals that require numerical evaluation, exacerbating accuracy and efficiency of the subsequent boundary-element implementation [46][47][48].In addition, the FT method requires all the interfaces of the media to be parallel, otherwise the convolution cannot reduce the PDEs to ODEs.
Unlike the FT method that requires elements along the interfaces, it is not necessary for the displacement discontinuity method (DDM), a branch of the BEM, to take care of interfaces [49][50][51][52], making itself perhaps the most efficient numerical method for an analysis of the crack problems [37].Therefore, after the statement of the crack problem in Section 2, the following Section 3 outlines the DDM at first, showing its high efficiency and strict limit for crack problems in general multilayered media, in which the number of layers is more than three.Then, Section 4 describes the direct method (DM) that is short for the direct boundary element method (DBEM), the other branch of the classical BEM, revealing its inapplicability for crack problems.After that, Section 5 presents a pre-treatment that turns the DM applicable for the crack problems and a DM-based method that solves the problems more efficiently than the conventional DM but still more slowly than the DDM.Having reviewed both the DDM and DM, the paper highlights a method that combines the DDM and DM in Section 6 so that it shares both the efficiency of the DDM and broad applicability of the DM after the pre-treatment, making it a promising candidate to solve the crack problems.Section 6 also provides numerical examples to demonstrate an even faster approximation with the combined method for a thin layer, which is one of the challenges for actual hydraulic-fracturing simulation.Finally, Section 7 summarizes the current advances of a boundary-element analysis of the crack problems and suggests future directions of the study.

Statement of Problem
Figure 1 shows the configuration of the crack problem in the paper.A pressurized crack under known pressure pc passes through a number of layers.Each layer, characterized by Young's modulus E and Poisson ratio v, is assumed homogeneous and isotropic individually [53].The objective is to determine the crack opening displacement D in the quasi-static problem.We shall compare the crack opening displacement obtained with a boundary-element method with those obtained with an analytical or semi-analytical solution, if any.The percentage error of the method to be tested is defined as where Dn is the crack opening obtained with the method to be tested, and Da represents the crack opening determined with the benchmark method.The labeling in the figure follows the convention of the works authored by Maier and Novati [54].The index (i) denotes the ith layer and ranges from 1 (the bottom layer) to N (the top layer).The superscript i represents the material properties and quantities of interest (traction and displacement) associated with the ith layer.The long dash line in the center implies the initial status of a slit-like crack characterized by two surfaces, one of which overlaps the other effectively.The two short dash lines at the two sides indicate the fictitious side boundaries, which are assumed in "undisturbed" regions and characterized by vanishing (zero) displacement, since the distance of them, h, is assumed to be infinite.
Matrices and column vectors are denoted by bold symbols.Thus, vectors of the nodal traction and displacement are indicated by p and u, respectively.The quantities pertaining to the lower (bottom), upper (top), side boundaries, and crack surfaces of each layer are denoted by the subscripts b, t, s, and c, respectively.
The boundary conditions of the multilayered system are usually of the type where 0 is the null matrix or vector, and p0 indicates the vector of the given nodal tractions on the very top surface of the system.The former of Equation ( 1) usually reflects an underlying rigid rock formation characterized by zero displacement beneath the bottom of the multilayered media [45,54].The boundary condition on the side boundaries is then As a result, the equations associated with the side boundaries can be split out of those pertaining to other boundaries [54].
The local coordinate defined in [41] is adopted in the paper, leading to the following interface continuity conditions (see the appendix for details).

Displacement Discontinuity Method (DDM)
The DDM, proposed by Crouch [49], exploits the physical analog of a slit-like crack and a distribution of a series of constant displacement discontinuities or a distribution of a series of edge-dislocation dipoles (Figure 2).As a result, the DDM only focuses on the elements on the crack surface by discretizing it as a unity.In other words, it is not necessary for the DDM to introduce elements along interfaces, if any, while other numerical methods have to take care of the interfaces, making the DDM perhaps the most efficient numerical method for an analysis of the crack problems [37].

Homogeneous Medium
The DDM was proposed under the assumption that (1) a crack lies in a homogeneous and isotropic medium and (2) equal loadings apply on the two crack surfaces (Figure 2a).The boundary-element equation of the DDM can be written in a concise form as where K is the coefficient matrix obtained from the DDM.Its explicit form is referred to in Ref. [49].It is worth noting that the K obtained from the theory of constant displacement discontinuities and that obtained from the theory of edge dislocation are mathematically the same.D is the nodal opening displacement of the elements on the crack surface.p is the known nodal pressure on the crack surface.Thus, the crack opening displacement is In addition, the analytical solution to the crack problem in a homogeneous medium is given as the following [55,56]: The percentage errors of the DDM with respect to the analytical solution with different numbers of the uniform crack element are illustrated in Figure 3, in which only the left-half part of the crack is displayed because of the symmetry.The abscissa is normalized by x/l, so that the coordinate −1 represents the crack tip. Figure 3b   Figure 3 demonstrates that a refined mesh leads to a more satisfying accuracy along the crack surface, except the tip, where the error of the DDM keeps around 26%.

The Higher-Order Tip-Element Approach
In order to improve the tip accuracy of the DDM, Shou and Crouch [57] developed the tip-element approach by replacing the constant displacement discontinuity (DD), also known as the zeroth-order element, at the tips with a higher-order quadratic DD (Figure 4).Constant DD is still assumed for elements other than at the tip, while the higher-order tip element assumes D~r 1/2 , where D is the opening displacement and r is the distance from a point in the tip element to its nearer tip point (x = −l for the left tip element and x = l for the right tip element).As a result, the approach captures the asymptotic behavior of the near-tip stress field due to the DD, leading to better tip accuracy.The tip-element approach yields a boundary-element equation in the same form as Equation ( 4), while the explicit form of the K is referred to in Ref. [57].The errors of the DDM with the tip element and the conventional DDM are displayed in Figure 5, which demonstrates significant improvement of the tip accuracy due to the higher-order tip element.The tip error decreases from around 26% to around 9%.The tip-element approach, however, is restricted to the homogeneous medium.Once the media becomes inhomogeneous, the approach fails.In addition, some attempts have been made to place higher-order elements on the whole crack surface for even better accuracy [58,59].Unfortunately, they are also limited to a homogeneous medium.

Two Bonded Half-Planes
Although the tip-element approach only applies to the homogeneous medium, the conventional DDM can be extended to the two bonded half-planes (Figure 6), since the Dundur's equation that characterizes the stress field of an edge dislocation in the two bonded half-planes can be used as an analytical fundamental solution [55,60].
The DDM for the two bonded half-planes results in a boundary-element equation in the same form as Equation ( 4), while the explicit form of the K is based on the Dundur's equation and can be referred to Refs.[55,60].A semi-analytical solution to the crack problem in the two bonded half-planes is available with an additional continuity condition [61,62]: In Figure 7, the abscissa is normalized using Thus, coordinate −3 and 1 represent the two crack tips.The result shows that the DDM matches well with the benchmark, except at the tip.Unfortunately, we cannot apply the tip-element approach here since the media becomes inhomogeneous.

Summary
The section has reviewed the DDM and its implementation for the crack problems in a homogeneous medium and two bonded half-planes.The method only requires elements on the crack surface and provides well-accepted accuracy, except at the tip.The tip-element approach improves the tip accuracy significantly, while it is not applicable to inhomogeneous media.In addition, the DDM is restricted in the cases of a homogeneous medium and two bonded half-planes since there are no existing analytical fundamental solutions to construct the boundary-element equations for the DDM in general multilayered media.Although some attempts [20,47,48] have been made to extend the DDM for general multilayered media, the extensions turn out to be first-order approximations, instead of full solutions [63,64].As a result, we turn to the DM, the other branch of the BEM, for a more applicable implementation.

Direct Method (DM)
Unlike the DDM that utilizes a physical analog of a slit-like crack and a distribution of a series of DDs or edge-dislocation dipoles, the traditional boundary-element method determines boundary-element equations directly by using the Kelvin's solution and the reciprocal theorem [65][66][67][68][69]. Thus, the method is also known as the direct method (DM) [41].The implementation of the DM for the multilayered media involves construction of a boundary-element equation for each layer and subsequent assembly of the equations with the continuity conditions along interfaces.Such a procedure usually generates an awfully large system of equations, of which all the unknowns are then solved simultaneously, undermining both accuracy and efficiency significantly [70].Thus, it becomes necessary to decompose the large final system of equations into smaller ones by using the chain-like properties, i.e., the continuity conditions on the interfaces of the multilayered media [53,71].Unfortunately, the DM cannot distinguish the effects of elements placed along one crack surface from those placed along the other surface.Thus, the DM cannot be applied to the crack problems directly [41].The implementation of the DM for the multilayered media that does not contain cracks, however, inspires subsequent studies on the DM for the crack problems in multilayered media.As a result, we shall still review the DM here in a relatively concise manner.

Transfer Stiffness Method (TSM)
Maier and Novati [53] developed a method to decompose the final large system of a boundary-element equation into smaller ones.The method sets up the boundary-element equation for each layer of the multilayered media with the DM in the form where H is the stiffness matrix obtained from the DM, of which the explicit form is referred to in [41].It expresses the tractions in terms of the displacements.u and p are the nodal displacement and traction of elements, respectively.The superscript i denotes the current layer (the ith layer).We can rewrite Equation ( 9) to ( 10) by partitioning the boundary elements according to the categorization of the boundaries: When the elements on the top and bottom surfaces are placed in pair (Figure 8a), Equation (10) can be rewritten to Equation (11) so that the quantities on the bottom interface can be expressed in terms of those on the top interfaces after some basic arrangements: where Equation (11) can be written into a more concise form as From Equation ( 12), we can solve the unknown displacements and tractions on the very top and bottom boundaries of the media using the chain-like property of the multilayered media, i.e., Equation (3), as a continuity condition on the interface: Equation ( 13) is a recursive formula.When i = N, we can solve pb 1 and ut N , since ub 1 and pt N are prescribed in Equation (1) as the boundary condition.Having obtained pb 1 and ut N , we can evaluate i from 2 to N−1 to obtain all the unknown displacements and tractions on the bottom and top interfaces of the layer that ranges from 2 to N−1.In other words, the TSM solves the unknowns by transferring the displacements and tractions on the bottom interface to those on the top interface of a layer with a sequential product of matrices from T 1 to T i .
Although the TSM manages to decompose the large final system of equations into smaller ones to improve efficiency, it suffers several computational issues that limit its application severely [45,54]: (1) It requires elements on the top and bottom interfaces in pair; (2) It categorizes the elements according to the type of unknown (displacement or traction).In practical applications, however, the displacement in mm and traction in kPa or MPa differs in magnitude of several orders.Therefore, an extra scaling that makes displacement and traction in the same order is necessary in case of ill-conditioning, undermining the efficiency of the method; (3) It results in a highly ill-conditioned matrix H when dealing with a thick layer.
In order to address the issues above, Maier and Novati [54] proposed the successive stiffness method (SSM) that can be implemented in a more robust manner.

Successive Stiffness Method (SSM)
The SSM follows the same path of the TSM until Equation (10), from which the SSM derives a recursive formula [54]: where Ĥ is the stiffness matrix of the ith layer.When i = 1, When i = 2~N, Based on Equation ( 14), the SSM derived another two recursive formulae for i = 2~N: Having obtained Equations ( 14), ( 17) and ( 18), we can firstly obtain ut N from Equation ( 14), since pt N is given as a boundary condition.Substitution of the solved ut N into Equations ( 17) and ( 18) with i = N solves both ub N and pb N .The continuity condition on the interface gives ut N−1 = −ub N and pt N−1 = pb N .Then, we can repeat the procedure to obtain ub N−1 and pb N−1 by replacing i = N by i = N − 1 in Equations ( 17) and (18).In such a recursive manner, we can compute all the unknowns successively from the very top layer to the very bottom one in a descending order, making the method quite convenient for numerical implementation.In addition, the method does not suffer the computational issues of the TSM.

Summary
The section has reviewed the DM and its implementation for multilayered media.The TSM decomposes the final large system of equations by utilizing the chain-like property (continuity condition on interfaces) of the media and categorizing the boundary element according to the type of unknown (displacement or traction).On the other hand, the SSM categorizes the boundary element by the categorization of the boundaries (bottom or top) and consequently overcomes the computational issues of the TSM.It is worth noting that neither the two methods are applicable for the crack problems directly, since the DM cannot distinguish the effects of elements placed along one crack surface from those placed along the other crack surface.The development of the SSM, however, inspires us to conceive a different approach to implement the DM so that it can be used for crack problems in multilayered media.

Pre-Treatment
The inapplicability of the DM to solve crack problems is attributed to the incapability of the DM to distinguish effects of the elements on the two overlapping crack surfaces [41].If we divide the media artificially along the crack surface, however, the media will be divided into two subdomains with a pair of imaginary interfaces denoted by I and I' (Figure 9).Then, we construct the boundary-element equation for each subdomain and assemble them with the continuity conditions on I and I'.The assembled equation turns out to be applicable to computing opening displacement of a slit-like crack [38,41].We can now solve the crack problem in a homogeneous medium by dividing the medium in the way shown in Figure 9 and obtain the assembled boundary-element equation to compute the crack opening displacement.The details of the assembly and the consequent discretization can be referred to in Ref. [72].The results obtained with the DM are compared to those obtained with the conventional DDM and the DDM with the tip element (Figure 10).
Figure 10 indicates that the DM shows the most satisfying accuracy except at the tip, since the integrals involved in the implementation of the DM can be evaluated analytically under the plane-strain condition [41].Although the higher-order tip-element approach yields the best tip accuracy, the DM that uses the zeroth-order element can still obtain a comparable accuracy at the tip.On the other hand, the requirement of the pre-treatment makes the DM much slower than the DDM (Table 1).All the computations in the paper were conducted with MATLAB 2021a on a desktop with an Intel(R) Core (TM) i9-13900KF 3.00GHz CPU and 32GB RAM.Table 1 clearly indicates that the DM is much slower than the DDM, not only because of the pre-treatment but also due to the necessity to take care of the elements on the resultant interfaces I and I', while the DDM only concentrates on the crack surface elements.
The pre-treatment shown in Figure 9 enables the application of the DM for crack problems in multilayered media.In Figure 11, the blue and red lines represent the two crack surfaces that effectively overlap each other, respectively.We cut the media along the crack surface (Figure 11a) artificially and acquire the resultant subdomains (Figure 11b).After that, we set up a boundary-element equation for each subdomain and assembled the equations with continuity conditions on all the interfaces.It is obvious that an increase in the number of layers leads to a more cumbersome implementation, as we need to set up the boundary-element equation for each subdomain and assemble them all together.The procedure finally results in an extraordinarily huge system of equations, further deteriorating the efficiency of the DM.In order to improve the efficiency, Long et al. [72] proposed a consecutive stiffness method to break the huge system of equations into smaller ones and enable a recursive formula for faster computation by utilizing the chain-like property, i.e., the interfacial continuity condition of the media.

Formulations and Algorithm
The CSM divides the media and obtains the resultant subdomain in the way shown in Figure 11.After setting up a boundary-element equation for each subdomain, the CSM assembles two subdomains of a given layer through the resultant interface I' and I'', if any.If the given layer is an inner layer, e.g., the ith layer in Figure 11b, the CSM simply patches the two systems of equations into one, since there is no resultant interface for the layer.Such a procedure eventually leads to a new system of equations for each layer, instead of a system of equations for the whole multilayered media.After that, the CSM derives recursive formulae to calculate unknowns for a faster computation by following the procedure of the SSM.In addition, the CSM derives a particular recursive formula that only focuses on the crack opening displacement uc so that all the other unknowns are eliminated during intermediate procedures, ensuring a further faster computation.Considering the limit of the page, we only present the four key equations of the method below.The detailed derivation of the equations is referred to in Ref. [72].
Equation ( 21) is the particular recursive formula in which i ranges from 2 to N. In Equations ( 19)-( 22), the matrices, such as Ĥ, H̃, H ̅ , A, C, etc., and the vectors at the righthand side, such as pc, p̂t, ûb, p c, etc., are all known.The details of those quantities are avail-able in Ref. [72].Having obtained the four equations, we can conduct a sequential computation to solve the displacements on the crack surfaces consecutively in each layer in a descending order from the very top layer to the very bottom one.The flowchart of the computation is shown in Figure 12.The calculated displacement of a crack surface in all the layers uc contains two parts: the displacement of the upper crack surfaces, denoted by uc + , and that of the lower crack surfaces, denoted by uc − .According to the local coordinate used in the work, the crack opening displacement Dc is the following (see the appendix for details):

Numerical Examples
The discretization of the CSM is referred to in Ref. [72].We compare the crack opening displacements obtained with the DDM and the CSM for the crack problem in two bonded half-planes (Figure 13).Figure 14 illustrates that both the DDM and the DM match the semi-analytical solution well along the crack surface, except at the tips.The CSM obtains tip errors that are lower than 20%, which is acceptable for practical engineering applications.Furthermore, the CSM leads to generally satisfactory results when the ratio of crack segments falls in the practical range, i.e., the value of b1/b2 falls in one order of magnitude.The durations of the DDM, CSM, and DM with the artificial subdivision are listed in Table 2.  2 indicates that the CSM is faster than the DM.On the other hand, the CSM that needs to take care of the elements on the interfaces is still slower than the DDM that only focuses on the elements on crack surfaces.The difference gets even larger with an increase in the number of elements.

Summary
The section has reviewed the CSM that implements the DM in a way that an artificial division of the multilayered media is made, enabling sequential computation to solve the displacements on the crack surfaces in each layer consecutively in a descending order from the very top layer to the very bottom one.Figures 10 and 14 illustrate acceptable accuracy of the DM even with the zeroth-order element because analytical evaluation is available for the integrals involved in the DM under the plane-strain condition [41], indicating that the DM-based CSM can be used as a convincing benchmark to test other methods for crack problems in multilayered media.Furthermore, Tables 1 and 2 have demonstrated the inefficiency of the CSM.
It is also worth noting that the inefficiency excludes the CSM for practical industrial applications, such as hydraulic-fracturing treatment that usually sets up the system of equations to conduct millions of time steps' computation because of an excessively small time step for numerical stability [73][74][75][76].For each time step, the newly introduced elements on the imaginary interfaces due to the artificial division result in a larger system of equations to manipulate [77][78][79], severely deteriorating the efficiency.In addition, when multiple cracks are involved, the CSM gets even slower and cumbersome for numerical implementation, since each crack requires an artificial division and consequent additional computation.Therefore, it is natural to conceive a method that shares both the efficiency of the DDM and broad applicability of the DM after the division.

Combined Boundary-Element Method
Ref. [19] incorporated the generalized Kelvin's solution, the fundamental analytical solution of the DM, to the DDM with the Fourier transform, illustrating combination of the two methods from the aspect of basic theory.On the other hand, Ref. [13] attempts to extend the DDM with a "dividing and patch" approach, showing a simpler combination.Inspired by the two works, Long and Xu [80] presented a method that sets up the stiffness matrix of each layer with the DDM and characterizes the effects of the interfaces with the DM.Thus, the method that combines the DDM and the DM shares both the efficiency of the DDM and broad applicability of the DM after the artificial division.

Formulations and Algorithm
Figure 15 schematically illustrates a given (the ith) layer after the deformation due to crack opening.The interfaces are sketched parallel for visual simplicity.Then, two artificial half-planes in the same material property marked in green and blue are patched with the layer on its original bottom and top interfaces (vertical dash lines), respectively.Having patched the two half-planes to the given layer, we apply tractions of the same magnitudes but in the opposite directions of the interfacial tractions.Thus, the interfaces of the two half-planes deform artificial nodal displacement vb i and vt i , respectively (Figure 15b).The nodal displacement discontinuities on the two interfaces are In addition, the resultant medium becomes homogeneous, enabling implementation of the tip-element approach for the resultant artificial crack.Thus, we set up the matrix K i for the elements at the crack tips marked by red circles in Figure 15b with the tip-element approach and construct K i for the non-tip elements with the conventional DDM: where the matrix K i can be partitioned by the categorization of the boundaries: where the matrices Kct i , etc., are determined with the configuration of the ith layer.Thus, they can be treated as known quantities.
On the other hand, we can also construct the boundary-element equation for the patched half-planes based on the DM: where the matrices Hbb i and Htt i are determined with the configuration of the patched halfplanes only.Therefore, they can also be treated as known quantities.Equations ( 25) and ( 27) are the starting point of the combined method.The method derives recursive formulae to obtain the final equation system that contains only the crack opening Dc for an even faster computation.Similar to the review of the CSM (Section 5.2), we only present the key equations of the combined method here, considering the limit of the page.The detailed explanation and derivation are referred to in Ref. [80].
When i = 1, the stiffness matrices of the combined method are where the matrices K̃, etc., are resultant matrices after fundamental calculations of the basic matrices K, etc., and H, etc., in Equations ( 25) and (27).When i = 2~N, the stiffness matrices of the combined method are where the matrices J, M, N, Q, etc., are also resultant matrices after fundamental calculations of the basic matrices K, etc., and H, etc., in Equations ( 25) and (27).The details of the matrices are referred to in Ref. [80].Having calculated K ̂cc N , K ̂ct N , K ̂tc N , and K ̂ttN from Equation (29), we can obtain the crack opening displacement D ̂cN : where D ̂cN contains the crack opening displacement of the elements from the bottom layer (i = 1) to the top layer (i = N).From Equations ( 28) and ( 29), we can conduct a sequential computation to compute the stiffness matrix K ̂ of each layer in an ascending order from the bottom layer to the top one.Having obtained the K ̂N, we can solve the crack opening displacement D ̂cN from Equation (30).The flowchart of the computation is shown in Figure 16.
Figure 16.Flowchart of the computation to solve crack problems in multilayered media with the combined method.

Two Bonded Half-Planes
The discretization of the combined method is referred to in Ref. [80].We firstly compare the percentage errors of the DDM, the CSM, and the combined method for the crack problem in the two bonded half-planes in Figure 17, while the parameters remain the same as in Section 3.3.The figure illustrates a satisfying match of the three methods to the semianalytical solution, except at the tips.The combined method incorporated with the tipelement approach gives the best tip accuracy.Furthermore, it does not affect the results of the combined method significantly when the crack length ratio remains in the practical range (up to 25) in different materials.In addition, the durations of the three methods are listed in Table 3, which indicates that the combined method improves the efficiency compared to the CSM, although the DDM that does not require elements on interfaces, if any, is still the most efficient.The second example is calculation of crack opening in general multilayered media, in which the number of layers is more than three.We chose five-layered media, in which a crack intersects the parallel interfaces perpendicularly (Figure 18).We set the crack segment as l1 = 12l2.The Poisson ratios of the softer and stiffer layers are 0.25 and 0.3, respectively.The pressure p0 remains a known constant.The percentage difference of the crack opening displacement obtained using the combined method and that obtained using the CSM with a different Young's modulus contrast and different numbers of the crack element on the thin layers are displayed from Figures 19-21, in which only half of the results are shown due to the symmetry.
All of the three figures clearly illustrate a "peak" difference on the thin-and-softer layers neighbored by thick-and-stiffer ones.Fortunately, we can alleviate the "peak" difference by refining the discretization.In addition, higher E contrast results in a higher "peak" difference and requires more refined elements to mitigate it consequently.When the discretization is refined sufficiently, the results obtained with the two methods are quite close to each other.

Approximation of Thin-Layer Scaling
In practical hydraulic-fracturing treatments, a fluid-driven crack sometimes passes through a few thin layers sandwiched by thicker ones.Such thin layers usually require awfully refined discretization for accurate simulation, slowing down the computation significantly.It is well accepted that efficiency and accuracy are both critical factors for practical industrial applications.A satisfying numerical simulation in industry should keep a balance between the two factors, instead of sacrificing one to satisfy the other.Thus, the thin layers might be scaled as a single relatively thick layer with average effective material properties as in Equation (31) for better efficiency with the combined method: where j is the index of the n thin layers, and lj is the thickness of the jth layer.Thus, the original crack problem in Figure 18 is approximated with that in Figure 22.The combined method placed only one crack element on the crack surface in the scaled layer characterized by E̅ and v ̅ for even faster computation.The calculated crack opening displacement is then compared to that obtained using the combined method with coarse mesh (one element on the thin layer) in the original five-layered media and that obtained using the CSM with refined mesh (five elements on the thin layer) in the original five-layered media.As indicated previously, the results obtained using the CSM with refined mesh can be used as convincing benchmarks.The comparisons with different contrasts in Young's modulus are shown from Figures 23 to 25.In the three figures, the crack opening displacement w0 is normalized using where w is the crack opening displacement and G1 is the shear modulus of the first layer.The three figures validate the approximation of the thin-layer scaling for more efficient computation when the contrast in Young's modulus remains up to 10-fold, since the approximated crack opening obtained using the combined method differs from the original crack opening obtained using the CSM with refined mesh by less than 20% (Figure 23).When the contrast in Young's modulus gets higher, the original crack opening obtained using the CSM with refined mesh on the thin-and-soft layers neighbored by a thickand-hard one reaches about twice of the approximated opening calculated using the combined method (Figures 24b and 25b).On the other hand, it is interesting to find the approximation to be valid for a contrast in Young's modulus up to 50-fold when a thin-andstiff layer is neighbored by thick-and-soft ones (Figures 24a and 25a).

Summary
The section has reviewed the combined method that constructs the stiffness matrix of each layer with the DDM and characterizes the effects of the interfaces with the DM.As a result, the method shares both the efficiency of the DDM and broad applicability of the DM after the division, as demonstrated in the numerical examples.In addition, the final system of equations only contains unknown crack opening displacements, leading to even faster computation.The combined method has also adopted the tip-element approach, ensuring a well-acceptable accuracy including at the tips, while the tip accuracy cannot be improved with refined discretization only.The numerical examples also validated an approximation of thin-layer scaling with the combined method under practical conditions for a more-efficient industrial simulation, for example, hydraulic-fracturing simulation.The crack intersects interface(s), if any, perpendicularly in the numerical examples, since some analytical or semi-analytical solutions are available in such conditions.In fact, the combined method can deal with more-complicated crack patterns, such as crack inclining to interfaces.Examples of using the method to simulate propagation of practical hydraulic fractures can be found in Refs.[12,51,78].

Summary and Outlook
The BEM is widely chosen for an accurate and efficient analysis of crack problems in multilayered media.The two branches of the BEM, the DDM, and the DM have demonstrated their strengths in efficiency and applicability, respectively.Although several attempts have been made to extend the DDM for general multilayered media and various approaches have been developed to accelerate computation conducted with the DM, the numerical examples have still revealed the weaknesses of the DDM and the DM in applicability and efficiency, respectively.Therefore, the review highlighted the combined method that constructs the stiffness matrix of each layer with the DDM and characterizes the effects of the interfaces with the DM so that it shares both the satisfying efficiency of the DDM and the broad applicability of the DM after the artificial division.In addition, the final system of equations of the combined method only contains unknown crack opening displacements, resulting in even faster computation.The tip-element approach was incorporated into the combined method, ensuring satisfactory accuracy including at the tips, while the tip accuracy can hardly be improved with refined discretization only.The numerical examples also validated an approximation of thin-layer scaling with the combined method to further improve the efficiency for practical industrial applications.
The BEM, however, is still limited currently, given challenges to solve crack problems accurately and efficiently in multilayered elastic media associated with more complicated conditions.As a result, pressing needs still exist for future study to explore the area more comprehensively.
1.The DDM is an ideal method for crack problems due to its high efficiency.The limited applicability, however, restricts its application severely.The attempts to extend it for general multilayered media by exploiting the method of images with a superposition scheme turn out to be first-order approximates.As a result, new approaches are encouraged for a full solution to extend the DDM. 2. Initiation and propagation of a fatigue crack make notable contribution to failure in industry.Computation of such phenomena in multilayered material continues pos-ing challenges for modeling and simulation.In addition, it gets even more complicated when material becomes brittle or quasi-brittle.As a result, it may be interesting to address such issues with the BEM. 3. It is getting popular to implement the PFM with the FEM for the analysis of complicated crack patterns in the latest decade.The former derives partial differential governing equations, which are then solved with the latter numerically.On the other hand, few studies have been conducted to incorporate the PFM into the BEM.Thus, it will be quite encouraging to combine the PFM with the BEM for a more efficient and more accurate analysis of more complex crack problems in multilayered media.
−|u + |, and u − = −|u − |.The opening displacement D is defined as the relative displacement of the lower element with respect to the upper one [41].Thus, Equations (A1) and (A2) correspond to Equations ( 4) and (23), respectively.The appendix only presents the essence of the local coordinate, the details of which can be found in Ref. [41].

Figure 2 .
Figure 2. Physical analog of (a) a slit-like crack and (b) a distribution of a series of constant displacement discontinuities or (c) a distribution of a series of edge-dislocation dipoles.
involves 100 crack elements, while only 11 elements are shown for visual simplicity.

Figure 3 .
Figure 3. Percentage errors of the DDM with (a) 5 elements and (b) 100 elements.

Figure 4 .
Figure 4.The DDM (a) without the tip-element approach and (b) with the tip-element approach.

Figure 5 .
Figure 5. Percentage errors of the conventional DDM and the DDM with the tip element: (a) 5 crack elements and (b) 100 crack elements.

Figure 6 .Figure 7 .
Figure 6.Physical analog of (a) a slit-like crack in two bonded half-planes and (b) a distribution of a series of edge-dislocation dipoles.The percentage error of the DDM with respect to the semi-analytical solution is shown in Figure7.In the computation, the material properties are set as E1 = 3.07 (GPa), v1 = 0.35, E2 = 68.95(GPa), and v2 = 0.3.The thickness of layer d is assumed to be 100 times of the larger one of b1 and b2, approximating the condition of the infinite d.In total, 10 uniform elements are placed on the shorter segment of the crack for a reasonably refined discretization.

Figure 8 .
Figure 8. Elements on the top and bottom interfaces are placed (a) in pair and (b) not in pair.

Figure 9 .
Figure 9.The artificial division: (a) division of the original media along the crack surface (dash line) and (b) the resultant subdomains and the imaginary interfaces I and I'.

Figure 10 .
Figure 10.Percentage errors of the conventional DDM, the DM, and the DDM with tip element: (a) 5 crack elements and (b) 100 crack elements.

Figure 11 .
Figure 11.Pre-treatment of the DM for crack problems in multilayered media: (a) original configuration and (b) resultant configuration after the division.

Figure 12 .
Figure 12.Flowchart of the computation to solve crack problems in multilayered media with the CSM.

Figure 13 .Figure 14 .
Figure 13.Schematic view of the crack problem in two bonded half-planes.The parameters remain the same as in Section 3.3.The percentage errors of the DDM and the CSM are displayed in Figure14.

Figure 15 .
Figure 15.Schematic view of the artificial patch: (a) the ith layer after the deformation and (b) the ith layer patched with two half-planes in the same material property.

Figure 18 .
Figure 18.Schematic view of the crack problem in five-layered media.

Figure 19 .Figure 20 .
Figure 19.The difference of the combined method and the CSM with 10-fold E contrast: (a) 1 crack element on the thin layer and (b) 5 crack elements on the thin layer.

Figure 21 .
Figure 21.The difference of the combined method and the CSM with 50-fold E contrast: (a) 1 crack element on the thin layer and (b) 5 crack elements on the thin layer.

Figure 22 .
Figure 22.Schematic view of the crack problem in the approximated three-layered media.

Table 1 .
Durations of the three methods for the case of homogeneous medium.

Table 2 .
Durations of the DDM, CSM, and DM for the case of two bonded half-planes.

Table 3 .
Durations of the three methods for the case of two bonded half-planes.