Analysis of Crack Problems in Multilayered Elastic Medium by a Consecutive Stiffness Method

Abstract

We propose a boundary-element-based method for crack problems in multilayered elastic medium that consists of a set of individually homogeneous strata. The method divides the medium along the slit-like crack surface so that the effects of the elements placed along one crack surface become distinguishable from those placed along the other surface. As a result, the direct method that cannot be directly applied for crack problems turns out to be applicable. After that, we derive a recursive formula that obtains a “stiffness matrix” for each layer by exploiting the chain-like structure of the system, enabling a sequential computation to solve the displacements on the crack surface in each layer “consecutively” in a descending order from the very top layer to the very bottom one. In our method, the final system of equations only contains the unknown displacements on the crack surface, ensuring the efficiency of the method. The numerical examples demonstrate better accuracy and broader applicability of our method compared to the displacement discontinuity method and more-acceptable efficiency of our method compared to the conventional direct method.

1. Introduction

Crack problems in multilayered elastic media have been found to be serve in many industrial applications, such as numerical stability analysis for safer construction [1,2,3,4,5], the characterization of fractured porous media for optimal transport efficiency [6,7,8,9,10], and the simulation of hydraulic fracturing for the economic development of unconventional reservoirs [11,12,13,14,15]. It is well accepted that the core component to model hydraulic fracturing is the ability to compute the elastic response of a pressurized crack intersecting a number of layers, each of which is assumed to be homogeneous and isotropic individually [15,16,17,18].

The boundary element method (BEM) is a widely-chosen method among the various numerical methods developed to calculate the opening displacement of cracks in crack problems [19]. Compared to methods that require the discretization of the whole domain of interest, for example, the finite element method (FEM), the BEM only discretizes the boundaries of a domain, leading to fewer elements in computation and enabling a more efficient computation consequently [20,21]. In addition, the BEM exploits analytical fundamental solutions, ensuring a potentially more-accurate computation [22]. As a result, determining appropriate analytical fundamental solutions becomes critical to the analysis of crack problems for ccurate construction of boundary element equations and efficient subsequent implementation.

The governing equations of multilayered media are a system of partial differential equations (PDEs) from the theory of elasticity [23]. Thus, it is natural to conceive a method that determines the analytical fundamental solution by reducing the PDEs to ordinary differential equations (ODEs) for easier computation [24]. Therefore, Fourier transform (FT) has been implemented in many works to reduce the PDEs of the system, in which all of the interfaces have to be parallel [16,23,24,25,26]. After the solutions of displacement and stress to the ODEs are determined in the FT field, a subsequent inverse FT on the solutions yields the analytical fundamental solutions in the spatial field. Improper integrals that require numerical evaluation, however, are usually contained in the analytical fundamental solutions due to the inverse FT, deteriorating the accuracy and efficiency of the following boundary-element implementation [27,28,29]. In addition, when the interfaces are not parallel, FT-based methods become inapplicable.

FT-based methods require the introduction of elements along the interfaces, while the displacement discontinuity method (DDM), a branch of the BEM that utilizes the physical analog of a slit-like crack and a series distribution of constant displacement discontinuities, does not need to take care of the interfaces [30], making itself perhaps the most efficient method to solve crack problems [23]. Unfortunately, the DDM only applies to homogeneous media or bi-materials (bounded by two half planes) since the analytical fundamental solutions to the DDM only exist in the two cases [30,31]. Shou et al. [28,29] attempted to extend the DDM to multilayered media. Their extension, however, turned out to be a first-order approximation instead of a full solution [32,33]. Apart from that, the application of DDM requires the loading on the two crack surfaces the same in magnitude but opposite in direction.

Unlike the DDM with its strict limitations, the conventional BEM, also known as the direct method (DM), is versatile for multilayered elastic media. The DM constructs the boundary-element equations for each layer using Kelvin’s solution and the reciprocal theorem as analytical fundamental solutions [34]. Then, it assembles all of the equations into a large system of algebraic equations by using the interface continuity. After that, all of the unknowns are solved simultaneously. It is obvious that the system of the final equations is going to grow to be tremendously huge with an increase in the number of layers, undermining the efficiency significantly [19,35]. As a result, Maier and Novati [36] proposed a transfer-matrix method that decomposes the large system of equations into smaller ones by exploiting the chain-like property of the multilayered system. The method, however, resulted in ill-conditioned matrices when dealing with thick layers. Therefore, the authors [37] developed a successive stiffness matrix method to remedy the issue. Unfortunately, neither of the two methods directly apply to elastic media containing a slit-like crack in which one crack surface coincides with the other effectively because the effects of elements placed along one crack surface become indistinguishable from those placed along the other surface [34]. Although some authors [22,34] have suggested that division of the media along the slit-like crack surface facilitates the application of the DM, no explicit procedure has been found in the literature, particularly when the DM is implemented with the chain-like structure of the multilayered media for better efficiency.

As a result, we propose a boundary-element-based method for crack problems in the multilayered elastic media. The method divides the media along the slit-like crack surface so that the effects of the elements placed along one crack surface become distinguishable from those placed along the other surface. Thus, the direct method that cannot be directly applied for crack problems turns out to be applicable. After that, we derive a recursive formula that obtains a “stiffness matrix” for each layer by utilizing the chain-like property of the system, enabling a sequential computation to solve the displacements on the crack surface in each layer “consecutively” in a descending order from the very top layer to the very bottom. Furthermore, the final system of equations in our method only contains the unknown displacements on the crack surfaces, ensuring a satisfying efficiency of the method. Therefore, our proposed method becomes decent in accuracy, acceptable in efficiency, and broad in applicability.

The remaining parts of the paper are presented as follows: Section 2 provides the statement of the crack problem; Section 3 derives the formulation of the consecutive stiffness method; Section 4 demonstrates the accuracy, efficiency, and applicability of the method under the plane-strain condition; and Section 5 summarizes the paper with conclusions.

2. Statement of Problem

The geometry and labeling of the crack problem in our work follows the conventions of the previous work (Figure 1) [37]. A multilayered elastic media consists of a number of layers, each of which is assumed to be homogeneous as well as isotropic and characterized by Young’s modulus E and the Poisson ratio v. A pressurized crack subjected to known pressure p_c passes through the layers. We need to determine the crack opening displacement D_c.

In Figure 1, the index (i) denotes the ith layer, ranging from 1 (the bottom layer) to N (the top layer). The superscript i indicates the material properties and quantities of interest (traction and displacement) pertaining to the ith layer. The long dashed line represents the initial status of the crack as a slit-like crack with two surfaces, one of which coincides with the other effectively. The two short dashed lines along the two sides are the fictitious side boundaries, the distance of which h is assumed to be infinite. As a result, they are assumed to be “undisturbed” and characterized by zero displacements.

The quantities associated with the lower (bottom), upper (top), side boundaries, and crack surfaces of each layer will be labeled by the subscripts b, t, s, and c, respectively. Matrices and column vectors are denoted by bold symbols. Vectors of the nodal traction and displacement in the problem are indicated by p and u, respectively.

The bottom surface of the system often represents an underlying stiff rock formation, characterized by zero displacement, beneath the bottom of the multilayered medium [16], while the very top surface can be subjected to any prescribed boundary conditions, usually known traction [16,37]. Therefore, the boundary conditions of the system are assumed to be traction-free on the top surface for the sake of simplicity:

u_b¹ = 0, p_t^N = 0

(1)

where 0 is the null matrix or vector.

The boundary condition on the side boundaries is then:

u_sⁱ = 0 (1 ≤ i ≤ N)

(2)

Therefore, the equations pertaining to the side boundaries can be split from those associated with other boundaries [37].

We have adopted the local coordinates defined in [34], resulting in the following interface continuity conditions (see Appendix A for details):

u_tⁱ = −u_bⁱ⁺¹, p_tⁱ = p_bⁱ⁺¹

(3)

3. Formulation of the Method

We present Figure 2a by rotating Figure 1 by 90 degrees clockwise and by removing the symbols of quantities for visual simplicity. The initial slit-like upper and lower crack surfaces, one of which coincides with the other effectively, are marked in blue and red, respectively, for visual distinguishment. We then divide the media by cutting it along the initial slit-like crack surfaces in order to distinguish the effects of the elements placed along one crack surface from those placed along the other in subsequent numerical implementation, leading to the resultant medium shown in Figure 2b.

We indicate the upper and lower parts of each layer after the cutting by the superscripts ’ and ‘’, respectively (Figure 3). The dividing also results in the horizontal interfaces, denoted by I, in the bottom and top layers, as shown in Figure 3a,c.

(1) For the bottom layer (i = 1):

We define α to represent the boundaries other than the interface I, i.e., the bottom, top, and crack boundaries. After that, we can formulate the boundary element equations based on the direct method for the two parts of the bottom layer:

K^1′u^1′ = p^1′, and K^1″u^1″ = p^1″

(4)

where K^1′ and K^1″ are the coefficient matrices of the upper and lower parts, respectively. They are determined by the configuration of the layer only. Their explicit forms can be found in the work of Crouch and Starfield [34].

We can rewrite Equations (4) to (5) by discretizing with boundary elements and partitioning them according to the categorization of the boundaries:

$$\left[\begin{array}{ll}{K}_{{\alpha }^{\prime }{\alpha }^{\prime }}^1& K_{{\alpha }^{\prime }{I}^{\prime }}^1\\ K_{I^{\prime }{\alpha }^{\prime }}^1& K_{I^{\prime }{I}^{\prime }}^1\end{array}\right]\left\{\begin{array}{l}{u}_{{\alpha }^{\prime }}^1\\ u_{I^{\prime }}^1\end{array}\right\}=\left\{\begin{array}{l}{p}_{{\alpha }^{\prime }}^1\\ p_{I^{\prime }}^1\end{array}\right\},\mathrm{and}\left[\begin{array}{ll}{K}_{{\alpha }^{″}{\alpha }^{″}}^1& K_{{\alpha }^{″}{I}^{″}}^1\\ K_{I^{″}{\alpha }^{″}}^1& K_{I^{″}{I}^{″}}^1\end{array}\right]\left\{\begin{array}{l}{u}_{{\alpha }^{″}}^1\\ u_{I^{″}}^1\end{array}\right\}=\left\{\begin{array}{l}{p}_{{\alpha }^{″}}^1\\ p_{I^{″}}^1\end{array}\right\}$$

(5)

Note that the continuity conditions defined by the local coordinates apply to interfaces I’ and I″:

u_I’ = −u_I″, p_I’ = p_I″

(6)

The substitution of Equation (6) into Equation (5) with subsequent algebraic arrangement leads to:

$$\left[\begin{array}{ll}{H}_{{\alpha }^{\prime }{\alpha }^{\prime }}^1& H_{{\alpha }^{\prime }{\alpha }^{″}}^1\\ H_{{\alpha }^{″}{\alpha }^{\prime }}^1& H_{{\alpha }^{″}{\alpha }^{″}}^1\end{array}\right]\left\{\begin{array}{l}{u}_{{\alpha }^{\prime }}^1\\ u_{{\alpha }^{″}}^1\end{array}\right\}=\left\{\begin{array}{l}{p}_{{\alpha }^{\prime }}^1\\ p_{{\alpha }^{″}}^1\end{array}\right\}$$

(7)

where $\{\begin{array}{l}{H}_{{\alpha }^{\prime }{\alpha }^{\prime }}^1=K_{{\alpha }^{\prime }{\alpha }^{\prime }}^1-K_{{\alpha }^{\prime }{I}^{\prime }}^1{\left(K_{I^{\prime }{I}^{\prime }}^1+K_{I^{″}{I}^{″}}^1\right)}^{-1}{K}_{I^{\prime }{\alpha }^{\prime }}^1\\ H_{{\alpha }^{\prime }{\alpha }^{″}}^1=K_{{\alpha }^{\prime }{I}^{\prime }}^1{\left(K_{I^{\prime }{I}^{\prime }}^1+K_{I^{″}{I}^{″}}^1\right)}^{-1}{K}_{I^{″}{\alpha }^{″}}^1\\ H_{{\alpha }^{″}{\alpha }^{\prime }}^1=K_{{\alpha }^{″}{I}^{″}}^1{\left(K_{I^{\prime }{I}^{\prime }}^1+K_{I^{″}{I}^{″}}^1\right)}^{-1}{K}_{I^{\prime }{\alpha }^{\prime }}^1\\ H_{{\alpha }^{″}{\alpha }^{″}}^1=K_{{\alpha }^{″}{\alpha }^{″}}^1-K_{{\alpha }^{″}{I}^{\prime }}^1{\left(K_{I^{\prime }{I}^{\prime }}^1+K_{I^{″}{I}^{″}}^1\right)}^{-1}{K}_{I^{″}{\alpha }^{″}}^1\end{array}$,

$\mathrm{and}$ we define a new coefficient matrix $H^1\equiv \left[\begin{array}{ll}{H}_{{\alpha }^{\prime }{\alpha }^{\prime }}^1& H_{{\alpha }^{\prime }{\alpha }^{″}}^1\\ H_{{\alpha }^{″}{\alpha }^{\prime }}^1& H_{{\alpha }^{″}{\alpha }^{″}}^1\end{array}\right]$.

As indicated in the statement of the problem u_b¹ = 0, the rows and columns pertaining to the bottom boundary can be extracted out of H¹, leading to a submatrix $\widehat{H}$¹ that contains only the quantities on the top and crack boundaries:

$$\left[\begin{array}{ll}{H}_{tt}{}^1& H_{tc}{}^1\\ H_{ct}{}^1& H_{cc}{}^1\end{array}\right]\left\{\begin{array}{l}{u}_t{}^1\\ u_c{}^1\end{array}\right\}=\left\{\begin{array}{l}{p}_t{}^1\\ p_c{}^1\end{array}\right\}$$

(8)

In Equation (8), the submatrices with the subscript without prime relate to the boundary elements on the corresponding boundaries of the two parts. For example, H_tc¹ is the submatrix representing the effects of the elements on the top boundaries of the two parts (both t′ and t″) on the elements on the crack surfaces of the two parts (both c′ and c″).

Equation (8) is in our preferable form, containing only the displacement of the crack surface, which is our target, and the quantities of the top boundary that links to the bottom boundary of upper layers via the continuity condition of Equation (3), enabling a recursive procedure to solve the displacement of the crack surfaces in each layer.

From Equation (8), we have:

u_c¹ = (H_cc¹)⁻¹(p_c¹ − H_ct¹u_t¹)

(9)

The substitution of Equation (9) and the continuity condition Equation (3) into (8) leads to:

$$\begin{array}{l}{p}_t{}^1=\left[{\left(H_{tt}{}^1\right)}^{-1}-H_{tc}{}^1{\left(H_{cc}{}^1\right)}^{-1}{H}_{ct}{}^1\right]u_t{}^1+H_{tc}{}^1{\left(H_{cc}{}^1\right)}^{-1}{p}_c{}^1\\ \hspace{1em} ={\widehat{K}}_{tt}{}^1{u}_t{}^1+{\widehat{p}}_t{}^1\\ \hspace{1em} =-{\widehat{K}}_{tt}{}^1{u}_b{}^2+{{\widehat{p}}_t}^1={\widehat{p}}_b{}^2\end{array}$$

(10)

The stiffness matrix of the bottom layer $\widehat{K}$_tt¹ = (H_tt¹)⁻¹ – H_tc¹(H_cc¹)⁻¹H_ct¹, and $\widehat{p}$_t¹ = H_tc¹ (H_cc¹)⁻¹ p_c¹. Equation (10) links the quantities of the bottom layer to those of the second layer.

(2) For the second layer (I = 2):

As shown in Figure 3b, the cutting results in two separate parts, each of which requires a system of linear algebraic equations:

K^2′u^2′ = p^2′, and K^2″u^2″ = p^2″

(11)

or in a more-explicit form:

$$\left[\begin{array}{lll}{K}_{b^{\prime }{b}^{\prime }}^2& K_{b^{\prime }{c}^{\prime }}^2& K_{b^{\prime }{t}^{\prime }}^2\\ K_{c^{\prime }{b}^{\prime }}^2& K_{c^{\prime }{c}^{\prime }}^2& K_{c^{\prime }{t}^{\prime }}^2\\ K_{t^{\prime }{b}^{\prime }}^2& K_{t^{\prime }{c}^{\prime }}^2& K_{t^{\prime }{t}^{\prime }}^2\end{array}\right]\left\{\begin{array}{l}{u}_{b^{\prime }}^2\\ u_{c^{\prime }}^2\\ u_{t^{\prime }}^2\end{array}\right\}=\left\{\begin{array}{l}{p}_{b^{\prime }}^2\\ p_{c^{\prime }}^2\\ p_{t^{\prime }}^2\end{array}\right\},\mathrm{and}\left[\begin{array}{lll}{K}_{b^{″}{b}^{″}}^2& K_{b^{″}{c}^{″}}^2& K_{b^{″}{t}^{″}}^2\\ K_{c^{″}{b}^{″}}^2& K_{c^{″}{c}^{″}}^2& K_{c^{″}{t}^{″}}^2\\ K_{t^{″}{b}^{″}}^2& K_{t^{″}{c}^{″}}^2& K_{t^{″}{t}^{″}}^2\end{array}\right]\left\{\begin{array}{l}{u}_{b^{″}}^2\\ u_{c^{″}}^2\\ u_{t^{″}}^2\end{array}\right\}=\left\{\begin{array}{l}{p}_{b^{″}}^2\\ p_{c^{″}}^2\\ p_{t^{″}}^2\end{array}\right\}$$

(12)

The explicit forms of the coefficient matrices K^2′ and K^2″ are also found in the work of Crouch and Starfield [34].

Since the two parts are separate, we can simply “patch” the two equations into one:

$$\left[\begin{array}{ll}{K}^{2^{\prime }}& 0\\ 0& K^{2^{″}}\end{array}\right]\left\{\begin{array}{l}{u}^{2^{\prime }}\\ u^{2^{″}}\end{array}\right\}=\left\{\begin{array}{l}{p}^{2^{\prime }}\\ p^{2^{″}}\end{array}\right\}$$

(13)

We define the new coefficient matrix for the second layer H²$\equiv \left[\begin{array}{ll}{K}^{2^{\prime }}& 0\\ 0& K^{2^{″}}\end{array}\right]$, which is then partitioned according to the categorization of the boundaries along the whole second layer, e.g., both t^2′ and t^2″ are categorized as t²:

$$\left[\begin{array}{lll}{H}_{bb}{}^2& H_{bc}{}^2& H_{bt}{}^2\\ H_{cb}{}^2& H_{cc}{}^2& H_{ct}{}^2\\ H_{tb}{}^2& H_{tc}{}^2& H_{tt}{}^2\end{array}\right]\left\{\begin{array}{l}{u}_b{}^2\\ u_c{}^2\\ u_t{}^2\end{array}\right\}=\left\{\begin{array}{l}{p}_b{}^2\\ p_c{}^2\\ p_t{}^2\end{array}\right\}$$

(14)

From the third row, we obtain:

u_c² = (H_cc²)⁻¹(p_c² − H_ct²u_t² − H_cb²u_b²)

(15)

The substitution of Equations (15) and (10) into the second row of Equation (14) leads to:

$$\begin{array}{l}{p}_b{}^2=H_{bt}{}^2{u}_t{}^2+H_{bb}{}^2{u}_b{}^2+H_{bc}{}^2{u}_c{}^2\\ \hspace{1em} =-{\widehat{K}}_{tt}{}^1{u}_b{}^2+{\widehat{p}}_t{}^1\\ \hspace{1em} =p_t{}^1\end{array}$$

(16)

Then, we can obtain:

$$u_b{}^2=-{(A^2-{\widehat{K}}_{tt}{}^1)}^{-1}{B}^2{u}_t{}^2+{(A^2-{\widehat{K}}_{tt}{}^1)}^{-1}[H_{bc}{}^2{\left(H_{cc}{}^2\right)}^{-1}{p}_c{}^2-{\widehat{p}}_t{}^1]$$

(17)

where A² = H_bc²(H_cc²)⁻¹ H_cb² − H_bb², and B² = H_bc²(H_cc²)⁻¹ H_ct² –H_bt².

The substitution of Equation (17) into (15) results in:

$$\begin{array}{l}{u}_c{}^2=-{\left(H_{cc}{}^2\right)}^{-1}[H_{ct}{}^2-H_{cb}{}^2{(A^2-{\widehat{K}}_{tt}{}^1)}^{-1}{B}^2]u_t{}^2\\ \hspace{1em} +{\left(H_{cc}{}^2\right)}^{-1}\{[I-H_{cb}{}^2{(A^2-{\widehat{K}}_{tt}{}^1)}^{-1}{H}_{bc}{}^2{\left(H_{cc}{}^2\right)}^{-1}]p_c{}^2+H_{cb}{}^2{(A^2-{\widehat{K}}_{tt}{}^1)}^{-1}{\widehat{p}}_t{}^1\}\end{array}$$

(18)

where I is the identity matrix.

The substitution of Equations (17) and (18) into the first row of (14) leads to the stiffness matrix equation of the second layer:

$$\begin{array}{l}{p}_t{}^2={\widehat{K}}_{tt}{}^2{u}_t{}^2+C^2{p}_c{}^2+D^{21}{\widehat{p}}_t{}^1\\ \hspace{1em} ={\widehat{K}}_{tt}{}^2{u}_t{}^2+{\widehat{p}}_t{}^2\\ \hspace{1em} =-{\widehat{K}}_{tt}{}^2{u}_b{}^3+{\widehat{p}}_t{}^2=p_b{}^3\end{array}$$

(19)

where $\{\begin{array}{l}{C}^2=H_{tc}{}^2{\left(H_{cc}{}^2\right)}^{-1}\left[I-H_{cb}{}^2{\left(A^2-{\widehat{K}}_{tt}{}^1\right)}^{-1}{H}_{bc}{}^2{\left(H_{cc}{}^2\right)}^{-1}\right]+H_{tb}{}^2{\left(A^2-{\widehat{K}}_{tt}{}^1\right)}^{-1}{H}_{bc}{}^2{\left(H_{cc}{}^2\right)}^{-1}\\ D^{21}=H_{tc}{}^2{\left(H_{cc}{}^2\right)}^{-1}{H}_{cb}{}^2{\left(A^2-{\widehat{K}}_{tt}{}^1\right)}^{-1}-H_{tb}{}^2{\left(A^2-{\widehat{K}}_{tt}{}^1\right)}^{-1}\\ {\widehat{p}}_t{}^2=C^2{p}_c{}^2+D^{21}{\widehat{p}}_t{}^1\\ {\widehat{K}}_{tt}{}^2=H_{tt}{}^2-H_{tb}{}^2{\left(A^2-{\widehat{K}}_{tt}{}^1\right)}^{-1}{B}^2-H_{tc}{}^2{\left(H_{cc}{}^2\right)}^{-1}\left[H_{ct}{}^2-H_{cb}{}^2{\left(A^2-{\widehat{K}}_{tt}{}^1\right)}^{-1}{B}^2\right]\end{array}$

Having obtained Equation (19), we can treat the remaining inner layers (i runs from 2 to N–1) in the same manner as we did for the second layer. Thus, we can obtain the recursive stiffness equation for the ith layer simply by replacing superscripts 1 and 2 in Equation (19) with i–1 and i, respectively.

$$p_t{}^i={\widehat{K}}_{tt}{}^i{u}_t{}^i+{\widehat{p}}_t{}^i$$

(20)

(3) For the top layer (i = N):

The cutting of the top layer results in an interface that is shared by the two resultant upper and lower parts, the same as the bottom layer. Thus, we chose an analogous way that is used from Equations (5) to (10) to treat the top layer, leading to the new coefficient matrix of the top layer:

$$H^N\equiv \left[\begin{array}{ll}{H}_{{\alpha }^{\prime }{\alpha }^{\prime }}^N& H_{{\alpha }^{\prime }{\alpha }^{″}}^N\\ H_{{\alpha }^{″}{\alpha }^{\prime }}^N& H_{{\alpha }^{″}{\alpha }^{″}}^N\end{array}\right]$$

(21)

where $\{\begin{array}{l}{H}_{{\alpha }^{\prime }{\alpha }^{\prime }}^N=K_{{\alpha }^{\prime }{\alpha }^{\prime }}^N-K_{{\alpha }^{\prime }{I}^{\prime }}^N{\left(K_{I^{\prime }{I}^{\prime }}^N+K_{I^{″}{I}^{″}}^N\right)}^{-1}{K}_{I^{\prime }{\alpha }^{\prime }}^N\\ H_{{\alpha }^{\prime }{\alpha }^{″}}^N=K_{{\alpha }^{\prime }{I}^{\prime }}^N{\left(K_{I^{\prime }{I}^{\prime }}^N+K_{I^{″}{I}^{″}}^N\right)}^{-1}{K}_{I^{″}{\alpha }^{″}}^N\\ H_{{\alpha }^{″}{\alpha }^{\prime }}^N=K_{{\alpha }^{″}{I}^{″}}^N{\left(K_{I^{\prime }{I}^{\prime }}^N+K_{I^{″}{I}^{″}}^N\right)}^{-1}{K}_{I^{\prime }{\alpha }^{\prime }}^N\\ H_{{\alpha }^{″}{\alpha }^{″}}^N=K_{{\alpha }^{″}{\alpha }^{″}}^N-K_{{\alpha }^{″}{I}^{\prime }}^N{\left(K_{I^{\prime }{I}^{\prime }}^N+K_{I^{″}{I}^{″}}^N\right)}^{-1}{K}_{I^{″}{\alpha }^{″}}^N\end{array}$

Unfortunately, the displacements of the bottom boundary of the top layer remain unknown. Thus, we cannot extract the equations associated with them out of the coefficient matrix, as we did for the bottom layer. As a result, we treated the coefficient matrix in the same manner as we did for the second layer to obtain explicit expressions of the displacement of the boundaries of the top layer. The results share the same forms of Equations (17)–(19) by replacing superscripts 1 and 2 with N − 1 and N, respectively, for the top layer:

$$u_b{}^N=-{(A^N-{\widehat{K}}_{tt}{{}^N}^{-1})}^{-1}{B}^N{u}_t{}^N+{(A^N-{\widehat{K}}_{tt}{{}^N}^{-1})}^{-1}[H_{bc}{}^N{\left(H_{cc}{}^N\right)}^{-1}{p}_c{}^N-{\widehat{p}}_t{{}^N}^{-1}]$$

(22)

$$\begin{array}{l}{u}_c{}^N=-{\left(H_{cc}{}^N\right)}^{-1}[H_{ct}{}^N-H_{cb}{}^N{(A^N-{\widehat{K}}_{tt}{}^1)}^{-1}{B}^N]u_t{}^N\\ \hspace{1em} +{\left(H_{cc}{}^N\right)}^{-1}\{[I-H_{cb}{}^N{(A^N-{\widehat{K}}_{tt}{{}^N}^{-1})}^{-1}{H}_{bc}{}^N{\left(H_{cc}{}^N\right)}^{-1}]p_c{}^N+H_{cb}{}^N{(A^N-{\widehat{K}}_{tt}{{}^N}^{-1})}^{-1}{\widehat{p}}_t{{}^N}^{-1}\}\end{array}$$

(23)

$$p_t{}^N={\widehat{K}}_{tt}{}^N{u}_t{}^N+C^N{p}_c{}^N+D^{\left(N\right)(N-1)}{\widehat{p}}_t{{}^N}^{-1}$$

(24)

where $\{\begin{array}{l}{A}^N=H_{bc}{}^N{\left(H_{cc}{}^N\right)}^{-1}{H}_{cb}{}^N-H_{bb}{}^N\\ B^N=H_{bc}{}^N{\left(H_{cc}{}^N\right)}^{-1}{H}_{ct}{}^N-H_{bt}{}^N\\ C^N=H_{tc}{}^N{\left(H_{cc}{}^N\right)}^{-1}\left[I-H_{cb}{}^N{\left(A^N-{\widehat{K}}_{tt}{}^{N-1}\right)}^{-1}{H}_{bc}{}^N{\left(H_{cc}{}^N\right)}^{-1}\right]+H_{tb}{}^N{\left(A^N-{\widehat{K}}_{tt}{}^{N-1}\right)}^{-1}{H}_{bc}{}^N{\left(H_{cc}{}^N\right)}^{-1}\\ D^{(N)(N-1)}=H_{tc}{}^N{\left(H_{cc}{}^N\right)}^{-1}{H}_{cb}{}^N{\left(A^N-{\widehat{K}}_{tt}{}^{N-1}\right)}^{-1}-H_{tb}{}^N{\left(A^N-{\widehat{K}}_{tt}{}^{N-1}\right)}^{-1}\\ {\widehat{p}}_t^N=C^N{p}_c{}^N+D^{(N)(N-1)}{\widehat{p}}_t{}^{N-1}\\ {\widehat{K}}_{tt}{}^N=H_{tt}{}^N-H_{tb}{}^N{\left(A^N-{\widehat{K}}_{tt}{}^{N-1}\right)}^{-1}{B}^N-H_{tc}{}^N{\left(H_{cc}{}^N\right)}^{-1}\left[H_{ct}{}^N-H_{cb}{}^N{\left(A^N-{\widehat{K}}_{tt}{}^{N-1}\right)}^{-1}{B}^N\right]\end{array}$

It is obvious from the stiffness matrix equation, Equation (24), that we can obtain u_t^N easily when p_t^N is prescribed as a boundary condition of the whole multilayered media since all of the other quantities on the right-hand side of Equation (24) are known:

$$u_t{}^N={({\widehat{K}}_{tt}{}^N)}^{-1}{p}_t{}^N-C^N{p}_c{}^N-D^{(N)(N-1)}{\widehat{p}}_t{{}^N}^{-1}$$

(25)

The substitution of the calculated u_t^N into Equations (22) and (23) solves u_b^N and u_c^N, respectively, since all of the other quantities on the right-hand sides of the two equations are known. Then, we can obtain u_t ^N−1 via the interface continuity condition u_t^N−1 = −u_b^N. After that, we can obtain u_b ^N−1 and u_c ^N−1 in the same manner as we calculated u_b^N and u_c^N, enabling a sequential computation that solves the displacements on the top boundary, crack surface, and bottom boundary of each layer in descending order from the very top layer to the very bottom one. It is worth noting that if u_t^N is given as a boundary condition of the whole system, then we can start the sequential computation directly without Equation (25). Therefore, our method is applicable for any prescribed boundary condition at the very top surface and are not limited to the traction-given boundary condition.

Because we are particularly interested in the displacement of the crack surface u_cⁱ, it is more efficient that the displacement u_bⁱ and u_tⁱ both be expressed in terms of u_cⁱ so that only u_cⁱ is involved in the computation. Equations (17) and (18) give the expressions of u_bⁱ and u_cⁱ by replacing superscripts 1 and 2 with i–1 and i, respectively, for the inner layer:

$$u_b{}^i=-{(A^i-{\widehat{K}}_{tt}{{}^i}^{-1})}^{-1}{B}^i{u}_t{}^i+{(A^i-{\widehat{K}}_{tt}{{}^i}^{-1})}^{-1}[H_{bc}{}^i{\left(H_{cc}{}^i\right)}^{-1}{p}_c{}^i-{\widehat{p}}_t{{}^i}^{-1}]$$

(26)

$$\begin{array}{l}{u}_c{}^i=-{\left(H_{cc}{}^i\right)}^{-1}[H_{ct}{}^i-H_{cb}{}^i{(A^i-{\widehat{K}}_{tt}{{}^i}^{-1})}^{-1}{B}^i]u_t{}^i\\ \hspace{1em} +{\left(H_{cc}{}^i\right)}^{-1}\{[I-H_{cb}{}^i{(A^i-{\widehat{K}}_{tt}{}^{i-1})}^{-1}{H}_{bc}{}^i{\left(H_{cc}{}^i\right)}^{-1}]p_c{}^i+H_{cb}{}^i{(A^i-{\widehat{K}}_{tt}{}^{i-1})}^{-1}{\widehat{p}}_t{}^{i-1}\}\end{array}$$

(27)

From Equation (27), we have:

$$u_t{}^i={\widehat{H}}_{tt}{}^i{H}_{tc}{}^i{H}_{cc}{}^i{u}_c{}^i-{\widehat{H}}_{tt}{}^i{H}_{tc}{}^i{\tilde{p}}_c{}^i$$

(28)

The substitution of Equation (27) into (25) leads to:

$$u_b{}^i={\widehat{H}}_{bc}{}^i{u}_c{}^i+{\widehat{u}}_b{}^i$$

(29)

where $\{\begin{array}{l}{\widehat{H}}_{tt}{}^i=-{\left\{H_{tc}{}^i\left[H_{ct}{}^i-H_{cb}{}^i{\left(A^i-{\widehat{K}}_{tt}{}^{i-1}\right)}^{-1}{B}^i\right]\right\}}^{-1}\\ {\tilde{p}}_c{}^i=\left[I-H_{cb}{}^i{\left(A^i-{\widehat{K}}_{tt}{}^{i-1}\right)}^{-1}{H}_{bc}{}^i{\left(H_{cc}{}^i\right)}^{-1}\right]p_c{}^i+H_{cb}{}^i{\left(A^i-{\widehat{K}}_{tt}{}^{i-1}\right)}^{-1}{\widehat{p}}_t{}^{i-1}\\ {\widehat{H}}_{bc}^i=-{\left(A^i-{\widehat{K}}_{tt}{}^{i-1}\right)}^{-1}{B}^i{\widehat{H}}_{tt}{}^i{H}_{tc}{}^i{H}_{cc}{}^i\\ {\widehat{u}}_b{}^i={\left(A^i-{\widehat{K}}_{tt}{}^{i-1}\right)}^{-1}\left[B^i{\widehat{H}}_{tt}{}^i{H}_{tc}{}^i{\tilde{p}}_c{}^i+H_{bc}{}^i{\left(H_{cc}{}^i\right)}^{-1}{p}_c{}^i-{\widehat{p}}_t{}^{i-1}\right]\end{array}$

The combination of Equation (28) and the interface continuity condition −u_tⁱ⁻¹ = u_bⁱ gives:

$${\mathbf{u}}_c{}^{i-1}=-{\left({\mathbf{H}}_{ct}{}^i{\widehat{\mathbf{H}}}_{tt}{}^{i-1}{\mathbf{H}}_{tc}{}^{i-1}{\mathbf{H}}_{cc}{}^{i-1}\right)}^{-1}{\mathbf{H}}_{ct}{}^{i-1}\left[{\widehat{\mathbf{H}}}_{bc}{}^i{\mathbf{u}}_c{}^i+{\widehat{\mathbf{u}}}_b{}^i+{\widehat{\mathbf{H}}}_{tt}{}^{i-1}{\mathbf{H}}_{tc}{}^{i-1}{{\tilde{\mathbf{p}}}_c}^{i-1}\right\}$$

(30)

Equation (30) is valid for the inner layers; thus, i ranges from 3 to N for Equation (30).

For the bottom layer, the combination of Equation (29) with i = 2 and Equation (9) leads to:

$$u_c{}^1={\left(H_{cc}{}^1\right)}^{-1}[p_c{}^1+H_{ct}{}^1{u}_t{}^1({\widehat{H}}_{bc}{}^2{u}_c{}^2+{\widehat{u}}_b{}^i)]$$

(31)

Once u_t^N is determined from Equation (25) or given directly as a prescribed boundary condition, we can solve u_c^N with Equation (23), enabling 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 with Equations (30) and (31). The solved crack surface displacement in all of the layers u_c contains two parts: the displacement on the upper crack surfaces, denoted by u_c⁺, and that on the lower crack surfaces, denoted by u_c⁻. According to the local coordinate defined in [34], the crack opening displacement D_c is (see Appendix A for details):

D_c = u_c⁻ + u_c⁺

(32)

4. Results and Discussion

We shall demonstrate the results of three numerical examples to validate the accuracy and efficiency of the consecutive stiffness method (CSM) that is proposed in the previous section. A slit-like crack subjected to known normal traction perpendicularly intersects the interfaces, if any, which are assumed to be parallel, so that other methods can be used as benchmarks. The contrasts in Young’s modulus and the thickness of the layer maintain an order of magnitude, which falls within the range of practical fracturing treatment [15,38,39].

We defined the percentage error of a method as follows:

$$\mathrm{Err}\%=\frac{D_1-D_0}{D_0}\times 100\%$$

(33)

where D₁ is the crack opening obtained by the method to be tested, and D₀ represents the crack opening determined by the benchmark method.

4.1. Discretization

Figure 4 illustrates the discretization of the lower part of the top layer in our numerical examples schematically. We divided the whole initial slit-like crack surface into elements of a uniform size e, which was chosen carefully so that there would be sufficient space for at least two elements that are uniform in size e on the resultant horizontal interface (Figure 4b). After that, we placed adaptive elements on the bottom boundary (interfaces of the multilayered elastic media) of the top layer. We made sure that ten uniform elements in the same size e were placed adjacent to the crack surface. Then, we placed a series of elements, of which an element is 1.2 times larger than its previous one in size. When the size of an element became more than 5e, we made the elements and all of the following ones 5e in size exactly, preventing computational issues due to sharp contrasts in element size. As a result, we acquired an interface of one part in a size that is only slightly larger than 50l, in which l is half the length of the initial crack surface, as shown in Figure 5a.

We placed adaptive elements on the horizontal interface in the similar way. If the interface is too short for 10 uniform elements in size e, we removed the last two elements and introduced two new uniform elements in size e_m by dividing the remaining part in two evenly (Figure 4b). When the interface was longer than 10 uniform elements in size e, we used the same approach to treat the elements at the corner of the horizontal interface and the top surface (Figure 4c), preventing sharp contrasts in element size. After finishing the discretization of the horizontal interface, we chose a similar approach to discretize the top surface of the multilayered system. We first placed three elements that were uniform in size e_m with subsequent adaptive elements of which the maximum size does not exceed 5e. When the discretization became close to the corner of the top surface and the side boundary, we used the treatment shown in Figure 4b,c to determine the size of the last two elements on the top surface. Having completed the discretization of the lower part, we discretized the upper part using symmetry with respect to the crack surface. As a result, the distance of the two side boundaries h becomes about 100l, as shown in Figure 5a, approximating the condition of remote side boundaries (h→∞).

The approach to discretize the top layer also applies to the bottom layer. As for the inner layers, if any, discretization becomes simpler since there is no horizontal interface due to the division. We can discretize the two vertical interfaces in the same manner that we discretized the bottom boundary of the top layer, as illustrated in Figure 5b.

4.2. Homogeneous Medium

The first example is the calculation of the opening displacement of a crack under a known constant normal loading p₁ in an infinite homogeneous and isotropic plane (Figure 6).

The analytical solution to crack opening in the homogeneous medium is [31]:

$$D_0=-\frac{4(1-v^2)}{E}{p}_1\sqrt{l^2-x^2}$$

(34)

Then, the percentage errors of the DDM and the DM with respect to the analytical solution with different numbers of crack elements are illustrated in Figure 7, in which only the left-half part of the crack has been presented due to symmetry. The abscissa is normalized by x/l so that the coordinate −1 represents the crack tip. Figure 7b involves 100 crack elements, while only 11 elements are shown for visual simplicity.

Figure 7 clearly demonstrates that more-refined discretization leads to better accuracy along the whole crack surface, except at the crack tip, where the DDM yields an error that is almost twice that of the DM. The DM shows more satisfying accuracy than the DDM, particularly at the tip, since the integrals involved in the implementation of the DM can be evaluated analytically under the plane-strain condition and because the singularity at the tip of the DDM is stronger than that of the DM [28,29,34].

4.3. Two Bonded Half Planes

The second numerical example is the computation of the opening displacement of a crack that passes through two bonded half planes under known normal loading (Figure 8), in which the abscissa is normalized by

$$\xi =\left\{\begin{array}{c}\frac{2x-b_1}{b_1},x<0\hfill \\ \frac{2x-b_2}{b_2},x\ge 0\hfill \end{array}\right.$$

(35)

We assume that the loading satisfies an additional continuity condition (Equation (36)) so that a semi-analytical solution can be used as a benchmark [40], although neither the DDM nor the DM requires such an extra condition for implementation.

$$\frac{1-v_1^2}{E_1}{p}_1=\frac{1-v_2^2}{E_2}{p}_2$$

(36)

In addition, the two half planes are characterized by E₁, v₁ and E₂, v₂, respectively. We set E₁ = 3.07(GPa), v₁ = 0.35; E₂ = 68.95(Gpa), and v₂ = 0.3 and present the percentage errors of the DDM and the consecutive stiffness method (CSM) with respect to the semi-analytical solution in Figure 9 with different b₁/b₂ values. The thickness of layer d is set to be equal to 100 times that of the larger one of b₁ and b₂, approximating the condition of infinite d. We placed 10 uniform elements on the shorter segment of the crack surface.

Figure 9 illustrates that both methods match the benchmark along the crack surface well, except at the tip, where the CSM yields 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 a practical range. The first and second numerical examples demonstrated the more-satisfactory accuracy of our CSM compared to the DDM.

4.4. General Multilayered Media

The third numerical example computes the opening displacement of a crack intersecting a few layers under constant loading p₀ in a general multilayered media in which the number of layers is equal to or larger than three. The DDM fails in such a case since there is no analytical fundamental solution for it. Thus, we compared the efficiency of the CSM and that of the conventional DM, which led to a large system of equations and solved all of the unknowns simultaneously. We set the crack being distributed uniformly in the system so that the length of each crack segment was l₁, as shown in Figure 10. Ten elements were placed on each crack segment. The thicknesses of the bottom and top layers were both set to 100l₁, approximating the condition of their infinite thicknesses. We set E₁ = 5(GPa), v₁ = 0.3; E₂ = 50(Gpa), and v₂ = 0.25.

The durations required by the CSM and the conventional DM to solve the crack opening displacement are listed in

Table 1. Durations of the CSM and the conventional DM.

Method	Duration for 3 Layers (s)	Duration for 5 Layers (s)
CSM	58	101
DM	194	406

. The computation was conducted with MATLAB 2014a on a laptop with an Intel(R) Core(TM) i7-2720QM CPU and 12GB RAM. The table clearly indicates that our CSM is more efficient than the conventional DM. The contrast becomes even larger with an increase in the number of layers, which leads to more elements being involved in the computation.

The three numerical examples clearly demonstrate the better accuracy and broader applicability of our method compared to the DDM, and more-acceptable efficiency of our method compared to the conventional DM.

5. Conclusions

We have proposed a numerical method based on the DM for crack problems in a multilayered elastic media consisting of a set of individually homogeneous strata. The method distinguishes the effects of elements placed along one crack surface from those placed along another crack surface by dividing the whole media along the initial slit-like crack surface. After that, we derived a recursive formula that obtains a “stiffness matrix” for each layer by exploiting the chain-like feature of the system, enabling a sequential computation to solve the displacements on the crack surface in each layer “consecutively” in a descending order from the very top layer to the very bottom one. In addition, the final system of equations in our method only contains the unknown displacements on the crack surfaces, ensuring the efficiency of the method. The numerical examples demonstrate the better accuracy and broader applicability of our method compared to the DDM and more-acceptable efficiency compared to the conventional DM, making it a convincing benchmark to test new methods for crack problems in multilayered elastic media. We are planning to conduct analysis of crack problems using our method in complex geometry, which fails the FT-based methods, and under complicated loading, which fails the DDM.