This paper describes the application of the complex variable meshless manifold method (CVMMM) to stress intensity factor analyses of structures containing interface cracks between dissimilar materials. A discontinuous function and the near-tip asymptotic displacement functions are added to the CVMMM approximation using the framework of complex variable moving least-squares (CVMLS) approximation. This enables the domain to be modeled by CVMMM without explicitly meshing the crack surfaces. The enriched crack-tip functions are chosen as those that span the asymptotic displacement fields for an interfacial crack. The complex stress intensity factors for bimaterial interfacial cracks were numerically evaluated using the method. Good agreement between the numerical results and the reference solutions for benchmark interfacial crack problems is realized.

1. Introduction

Composite materials are found in a variety of important structures such as adhesive joints, composite laminates, and electronic and optic components. The abrupt change in properties at the interface is one source of failure of composite materials; the failure behavior can be rationalized by stress analysis and fracture mechanics [1].

Stress intensity factors are important parameters in the investigations of fracture toughness of composite materials. The stress intensity factors have been well established for a linear elastic homogeneous material [2]. But for a bimaterial composite, the situation is more complicated since the field equations are complex and the calculation of opening and sliding mode parameters is far from obvious. A bimaterial is perfectly bonded with no specific material parameters describing the interface. Interface toughness strongly depends on the relative amount of crack surface shearing to opening; a complex stress intensity factor has to be employed for the interface crack.

Numerical methods are the primary methods of choice in the analysis of cracks in bimaterial [3, 4]. In these methods, the finite element method (FEM) is a popular numerical method for calculating the stress intensity factor, but the process of FEM cannot be completed easily even if the automatic mesh generation technique is used, because the finite element mesh should be regenerated whenever the location or length of the assumed crack is changed. Hence, meshless methods have been developed and used with increasing frequency over the last decade [5]. As compared with the standard FEM approach, meshless method can avoid the distortion of mesh when extreme large deformation is encountered [6] and provides an efficient means for addressing high gradient problems, such as that occurring in strain concentration [7].

Meshless methods can mainly be divided into two kinds; one is based on moving least-square (MLS) approximations and the other is based on the partition of unity. The first methods include element-free Galerkin method [8], meshless local Petrov-Galerkin method [9, 10], reproducing kernel particle method [11], wavelet particle method [12], radial basis functions [13], moving particle finite element method [14], and complex variable meshless method [1517]; the second include hp-clouds [18] and the partition of unity finite element method [19, 20], and the meshless manifold method [21]. Other meshfree methods include the boundary element-free method [22, 23], the particle in cell (PIC) method [24], and the generalized finite difference method [25]. A comprehensive review of mesh-free methods can be found in [26].

In practical applications, the structural analysis system based on meshless methods to evaluate the fracture mechanics parameters of practical models is developed. Once crack information such as shape, location, and number is added to the model, the stress analysis is performed to evaluate results such as stress, strain, and displacement distribution. Lin used partition of unity method to model the discontinuous material [27]. Cheng and Li presented complex variable moving least-square approximation in conjunction with meshless for stress intensity factor analysis of an interfacial crack [2]. Li et al. applied enriched meshless manifold method, which is a combination of the partition of unity method and the finite cover approximation theory, to two-dimensional crack modeling [28]. Sukumar et al. investigated partition of unity enrichment for bimaterial interface crack [29].

The complex variable numerical manifold method (CVMMM) uses complex variable moving least-squares (CVMLS) approximation to construct meshless approximations [30, 31]. The CVMLS, which can avoid forming the ill-conditional or singular equations and can improve the computational efficiency rapidly, first was presented by Cheng et al. [32]. Now various complex variable meshless methods have been undertaken that marked the beginning of intensive research efforts, such as the complex variable meshless method [33, 34], the complex variable boundary element-free method [35], the complex variable reproducing kernel particle method [36, 37], and the complex variable element-free Galerkin method [38, 39].

In this paper, we develop the CVMMM for bimaterial interface cracks. This work extends the capabilities of the CVMMM to the analysis of cracks that lie at the interface of two elastically homogeneous isotropic materials. A detailed description of the numerical implementation of the CVMMM for 2D crack modeling in isotropic media is given in [30]. In this study, interfacial crack modeling within the CVMMM framework is proposed, and a numerical study of the accuracy and robustness of the method is conducted.

2. Interfacial Cracks at Bimaterial

There is a fundamental difference between the analysis and interpretation of the stress intensity factors for interface cracks in bimaterials as compared to cracks in homogeneous materials. Bimaterials exhibit a coupling of tensile and shear effects, and the singular stresses are oscillatory in the vicinity of the interface crack tip. The stress field is characterised by a complex stress intensity factor, , together with the bimaterial constant, , relating the elastic properties of the two materials.

Consider bimaterial joined along the x1-axis as shown in Figure 1. Assuming traction free crack surfaces, the complex stress field in the vicinity of the crack tip is where the and are the values 1 and 2. The angular functions and equal unity along the interface ahead of the crack tip . For the angular functions can be found in [40]. The expression includes the bimaterial constant, , defined as where is one of the Dundurs’ parameters, which is defined as where is the corresponding shear modulus; a constant can be given as

The complex stress intensity factor is defined as where is the phase angle. The notations and instead of and are adopted here to emphasise that a bimaterial system is studied.

The stresses, and , at the interface directly ahead of the tip, at , are given as

The associated crack surface displacements, and , at a distance behind the tip are given as where

3. CVMMM for a Bimaterial Interface Crack

3.1. The CVMMM

Key techniques used in the construction of the CVMMM are the finite cover theory and the covering function in manifold method.

The finite cover used in the CVMMM is referred to as the physical cover and the mathematical cover, respectively. The physical cover is used to define element domains, whereas the mathematical cover is used to construct approximation covering functions. The mathematical cover is often constructed from a simple pattern such as triangle and circle. The construction of a physical cover involves the computations of intersections of the mathematical cover and the boundary of the problem domain. However a mathematical cover is made up of the influence domain of the node, whereas a physical cover is used to define the boundary and discontinuities in the problem domain in the meshless manifold method. And the mathematical cover of each node is a local cover, which is used to construct a local approximation by the CVMLS approximation. The physical cover only constitutes the solving domain of the problem.

The covering function is expressed by the complex variable moving least-squares (CVMLS) approximation. Through the introduction of complex number the CVMLS can reduce the degrees of freedom of basis function; the difference between CVMLS and MLS is that the CVMLS approximation is for a vector and the MLS approximation is for a scale function. The trial function can be written in the following form [32]: where is a basis function and are the corresponding unknown coefficients that are associated with node i and a specific geometric entity (such as crack).

The unknown coefficients in the CVMLS approximation are fewer than in the MLS approximation. For the linear basis of the CVMLS, the basis function needs two unknown coefficients. For the quadratic basis, the basis function needs three unknown coefficients. For an arbitrary point in the domain, we need fewer nodes with domains of influence that cover the point, and thus we also require fewer nodes in the whole domain.

3.2. Enriched CVMMM for a Bimaterial Interface Crack

For a bimaterial interface crack, the Cartesian components of near-tip asymptotic displacement fields can be obtained. The crack-tip displacement fields in the upper-half plane (replace by for the lower-half plane) are The terms in (10) can be described as where and is the distance from the point to the tip of the crack and is the angle from the tangent to the crack path at the crack tip.


In the enriched trial function method, a single crack in two dimensions is considered, the trial function is obtained from the trial function (9), and the terms that describe the displacement field near the tip of a crack are considered. We have where and are the approximations of the exact solution and and are unknown parameters associated with the interfacial crack.

Letwe have which are the near-tip enrichment functions to model an interfacial crack in bimaterial media. Detailed discussion on the theoretical need for CVMMM in the modeling of cracks is provided in [31].

3.3. Weak Form and Discrete Equations

Consider bimaterials labeled as and with domains and , respectively (Figure 1). The problem domain and is the interface. A crack is assumed to occupy part of . From the principle of virtual work, the weak form for the discrete problem is where and are the approximating trial functions, is the small strain tensor, is the Cauchy stress tensor, and is the first variation operator. is the prescribed tractions.

The CVMMM discrete space contains functions that are discontinuous across the crack . Following (16), the approximating trial functions used to model bimaterial interface cracks are written in the form where

Using the arbitrariness of nodal variations, the final linear algebra equations system is obtained as follows: where is the vector of nodal unknowns and and are the stiffness matrix and the load vector, respectively. For the specific forms of and , see [30, 31].

4. Numerical Examples

Numerical examples of bimaterial interface crack problems using CVMMM are presented in the following section. Results for a semi-infinite plate containing the bimaterial interface crack in the centre under tensile load are presented in Section 4.1. The stress intensity factors are evaluated and compared with the exact solution. Section 4.2 shows the numerical results for a plate containing the bimaterial interface crack in the centre under tensile load. Further, the stress intensity factors are evaluated and compared with reference solution by BEM [41]. Section 4.3 shows the numerical results for the test piece with a bimaterial interface crack under residual stress, and the results are compared with other solutions by Yuuki using BEM [42].

4.1. Bimaterial Interface Central Crack in a Semi-Infinite Plate under Tensile Load

In this example, a semi-infinite plate with a bimaterial linear central crack was subjected to a tensile load, as shown in Figure 2. The stress intensity factors are evaluated and the results are compared to the analytical solution of the problem of an isolated crack of length 2 at the interface between two remote stresses; the width of the plate is set to 2, and the crack length and plate width are set to 6 mm and 60 mm in the present calculation. Elastic moduli of both materials are set to 68.5 GPa and 137 GPa, and Poisson’s ratio of both materials is fixed at 0.3. The height of plate is set to 78 mm. The two load cases studied were normal tension and mixed mode loading under plane strain conditions. Figure 2 shows the geometry and loading, with and denoting remote normal tension and shear stress, respectively. For the centre crack tip, the exact solutions of stress intensity factor are given as

Figure 3 shows the CVMMM nodes in the vicinity of the crack. The length of the smallest nodes is chosen to be equal to one-tenth of the crack length. The employed meshless nodes are 861.

Material combinations defined by ranging from −0.1 to 0 were treated. This corresponds to a range of the second constant of Dundurs, , of approximately 0–0.25 and should well cover a wide range of material combinations, wherein is considered a reasonable range.

Figures 4 and 5 show the stress intensity factors obtained for the pure tension and the mixed mode loading cases, respectively, including the analytical results from (21).

Define the relative error as where is the analytical result for mode and is the numerical solution. ranges from 0.25% to 0.5% for pure tension loading and from 0.5% to 3% for mixed loading. For , the deviations from the analytically computed values were between 1% and 3% for both load cases. The error increases with decreasing . The accuracy is influenced by the length of the smallest nodes. The same length is used throughout, but the stress gradient across the interface increases with decreasing . It should also be noticed that an error still remains at . This is an effect of the CVMMM discretisation.

Figures 4 and 5 show the comparisons between the analytically obtained mode mixity and results from the present method. From Figures 4 and 5 we can find that the bimaterial method can work well for the range of bimaterial constants studied.

The mode mixities from these calculations are included in Figures 4 and 5. The results differ significantly from those obtained by the present method taking the material mismatch, that is, coupled stresses, into account.

The CVMMM formulation for the energy release rate would be the same as for a homogeneous material, but the relation for the phase angle (or the mode mixity) could not be found so easily. In fact, the phase angle would have to be calculated by a separate numerical procedure. This would have increased the complexity of the present method considerably, thus violating the simplicity of the herein presented results.

4.2. Bimaterial Interface Central Crack in a Plate under Tensile Load

In this example, a rectangular plate with a bimaterial interfacial central crack is subjected to a tensile load, as shown in Figure 3. The width 2 and length 2 of the plate are set to 100 and 200 mm, respectively. Three cases of crack length 2 (40, 60, and 80 mm) are considered. Young’s modulus of the lower material (material 1) is set to 205.8 GPa. The ratio of Young’s modulus of material 2 () to that of material 1 () is defined as and is varied from 1.0 to 100.0. For both materials 1 and 2, Poisson’s ratio is fixed at 0.3. The tensile load is fixed at 9.8 MPa.

The employed finite element mesh (number of nodes: 861 and number of elements: 800) is shown in Figure 6. CVMMM analysis is performed using these nodes under the assumption of a plane stress state, and conventional FEM analysis was performed to verify the results. In FEM analysis, the double nodes are distributed in the region of the interface crack, whereas, in CVMMM, the nodal property is changed so as to enrich the region in the vicinity of the interface crack.

The stress intensity factors and of the problem are calculated using the CVMMM and FEM. The results are divided by to obtain nondimensional values and . The results are shown in Table 1 and are also compared with the reference solution [33]. The results are more comparable with the reference solution.

4.3. Test Sample with Bimaterial Interfacial Crack Subjected to Residual Stress

In this example, an edge crack in jointed dissimilar materials (as shown in Figure 7) is subjected to a uniform change in temperature, and the stress intensity factors of the interfacial crack are evaluated. Two types of glass-type material are assumed, and the properties of materials are shown in Table 2. No mechanical load is considered, and only thermal loads induced by a temperature change are taken into account.

Four crack lengths a (1, 2, 3, and 4 mm) are considered. The nodes shown in Figure 3 are used in all cases, and a plane stress state is assumed. Moreover, the stress intensity factors and are calculated using the CVMMM. The results are normalized as follows: where and are the coefficients of thermal expansion of materials 1 and 2, respectively.

Table 3 shows the numerical results which are compared with those from two-dimensional BEM shown in [42]. As can be seen, appropriate results are obtained by the proposed analysis using CVMMM.

5. Concluding Remarks

In this paper, CVMMM is applied to the two-dimensional stress intensity factor analysis of the bimaterial interface crack problem. This work has extended the capabilities of the extended CVMMM to the analysis of cracks that lie at the interface of two elastically homogeneous isotropic materials. The two-dimensional near-tip asymptotic displacement functions are added to the CVMMM approximation using the framework of complex variable moving least-squares approximation. The crack-tip enriched functions are chosen as those that span the asymptotic displacement fields for an interfacial crack. The stress intensity factors and of the interface crack evaluated by the CVMMM provide results comparable with those by exact solution and the reference BEM solution. The calculated results are more accurate than those obtained by FEM and BEM.

The proposed modelling method is expected to be implemented in a developed structural analysis system based on CVMMM, and solution of additional practical problems is expected. Moreover, more realistic simulations of the crack extension from the debonded region at the bimaterial interface to internal region are expected. In such cases, enrichment by trial function to express junctions may be used.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.


The work is supported by the National Natural Science Foundation of China through Grant 11271234 and the Science Program of Shandong Provincial Education Department through Grant J13LJ51.