Abstract
The interfacial crack in bimaterials is a very interesting problem for composite materials and which has received particular attention from several researchers. In this study, we will propose a numerical modeling of the interfacial crack between two orthotropic materials using a special mixed finite element. For the calculation of the energy release rate, a technique, based on the association of the present mixed finite element with the virtual crack extension method, was used. The numerical model proposed, in this work, was used to study a problem of interfacial crack in bimaterials. Two cases were treated: isotropic and orthotropic bimaterials. The results obtained, using the present element, were compared with the values of the analytical solution and other numerical models found in the literature.
Keywords: Interfacial crack, Mixed finite element, Virtual crack extension method, Energy release rate, Orthotropic bimaterials
The interfacial fracture is a complex phenomenon which is still badly understood, so it would already be enough to justify its study. Indeed, the interface located between two dissimilar materials is, on the mechanical level, a weak point: when these materials are requested by stresses, of thermal origin for example, the fracture of the interface is a mode usually observed. Moreover, one knows little about the mechanical conditions which lead to this fracture. A comprehension of the interfacial fracture thus represents a significant stake in the field of composite materials.
The problem of the interfacial crack in isotropic bimaterials has been treated by many researchers. We can cite, for example, the work of Williams [1], Erdogan [2,3], England [4], Rice and Sih [5], Hutchinson et al. [6], Rice [7] and Suo and Hutchinson [8].
The cracks along the interface between two anisotropic plates were initially treated by Gotoh [9]. The case of plane deformation of interfacial crack between two anisotropic materials was studied by Clements [10], Willis [11], Qu and Bassani [12], Suo [13] and Ni and NematNasser [14]. Bassani and Qu [15] have explicitly resolved the special case of Griffith's problem and the solution of the general problem has been found by Suo [13] and Qu and Li [16]. The crack path in the anisotropic medium was studied theoretically and numerically by Gao et al. [17], a weak plane model was adopted to characterize the anisotropic fracture toughness and the maximum energy release rate criterion was chosen to predict the crack path. The problem of interfacial cracks in anisotropic bimaterials was also treated by Wang et al. [18], Juan and Dingreville [19].
Based on anisotropic elasticity, Tanaka et al. [20] evaluate the energy release rate by the modified crack closure integral of the finite element method, and convert to the stress intensity factor for the cases of cracks on elastic symmetrical planes. Two approaches have been described by BanksSills and Ikeda [21] for considering an interface crack between two anisotropic materials. Both approaches have been used for orthotropic and monoclinic materials. The problem of cracked orthotropic bimaterial was also studied by Bouchemella et al. [22]. Fracture analysis of orthotropic cracked media was investigated by applying the recently developed Extended IsoGeometric Analysis (XIGA) [23] using the Tspline basis functions [24]. The same method XIGA was used by Habib et al. [25](2017) for the analysis of static fracture behaviour for a crack in orthotropic materials.
Khatir and Wahab [26] used an inverse algorithm based on Proper Orthogonal Decomposition (POD) and Radial Basis Functions (RBF) for single and multiple cracks identification in plate structures. The inverse analyses combine experimental fracture mechanics tests with numerical models based on the XIGA method. The eXtended IsoGeometric Analysis combined with Particle Swarm Optimization (PSO) have been used for crack identification in twodimensional linear elastic problems (plate) based on inverse problem [27].
In this paper, a numerical modeling has been proposed to study the interfacial crack between two orthotropic materials. This model uses a twodimensional mixed finite element developed in a natural plane. It is an element with 7 nodes: 5 displacement nodes and 2 stress nodes. The proposed model was used to calculate the energy release rate in a cracked orthotropic bimaterial using a technique that combines the present element with the virtual crack extension method. In this work, two cases of interfacial cracks were treated: an isotropic bimaterial and an orthotropic bimaterial. The results obtained, using the present mixed finite element, were compared with the values of the analytical solution and other numerical models found in the literature.
The bimaterial has been discredized using a special mixed finite element RMQ7 (Reissner Modified Quadrilateral) as shown in Figure 1(a). The present mixed finite element used in this study is twodimensional element with seven nodes: five displacement nodes and two stress nodes as shown in Figure 1(b). The node 5 coincides with the crack tip. This element was developed by Bouzerd [28], in the physical (x, y) plane, and was reformulated and validated by Bouziane et al. [29] in a natural (ξ, η) plane.
Displacement for the present mixed finite element can be given by

(1) 
where are the shape functions and is the nodal displacement corresponding to node i. For the present element, the shape functions are given as follows
,
, , 
(2) 
The element stress component is approximated by

(3) 
where [M] is the matrix of interpolation functions for stresses and {τ} is the vector of nodal stresses.
For the RMQ7 element (see Figure 1(b)), the shape functions M_{i2}, used to evaluate σ_{12} and σ_{22}[29] for nodes 6 and 7 are obtained by

(4) 
The element stiffness matrix [K_{e}] is given by the following expression

(5) 
where the submatrices, and , are given by the following relations

(6) 
where [S] is the compliance matrix, [M] is the matrix of interpolation functions for stresses, [B] is the straindisplacement matrix of shape function derivatives, t is the thickness, A^{e} is the element area and T indicate the matrix transpose.
(a) Discretization of bimaterial  (b) RMQ7 element 
Figure 1. Discretization of bimaterial and RMQ7 element 
The virtual crack extension method, associated with the mixed finite element RMQ7, is used to calculate the energy release rate G [28]. In this technique, the first calculation of the deformation energy is carried out in the initial configuration of the crack. The crack is then moved an infinitesimal distance in the direction of its axis. The deformation energy is evaluated again in the second configuration, the energy released during this crack length variation is

(7) 
The energy release rate G will be evaluated thereafter starting from the relation

(8) 
Calculation by the virtual crack extension method requires two finite element analysis. The use of the RMQ7 element makes it possible to introduce one mesh for the calculation of the energy release rate, which represents a considerable profit in computing times and setting data compared to the traditional techniques which use two meshes [28].
Indeed the intermediate displacement node of the RMQ7 element is associated to crack tip, and consequently the length of crack "a" can be increased by a quantity while acting inside strict of the crack element by translation of the tip crack node without disturbing the remainder of the mesh.
With the assumption on materials and displacements (linear elastic behaviour and small displacements), the solutions and obtained in the structure with a crack length "a" and in the same structure with a crack length are as close as the disturbance is small compared to dimensions of the crack element. We can thus write with a good approximation

(9) 
Several calculations on simple examples enabled us to confirm the relation Equation (9), which is theoretically coherent and physically acceptable, considering the assumptions used.
If we consider that the external loading does not vary during the increase , the energy release rate is calculated as follows

(10) 
where and represent respectively the deformation energy of the cracked structure in the configuration and "a".
In its discretized form, the deformation energy is written

(11) 
with:
ne : total number of elements in discretized structure,
: vertical vector containing the nodal values of element i,
: elementary matrix of element i, and the exponent "T" indicates the transposed vector.
By substitution of Equation (11) in Equation (10), the expression of the energy release rate G becomes

(12) 
Taking account of Equation (9), the expression Equation (12) can be written in the following form

(13) 
and as only the crack element is disturbed, then G results more simply in the relation

(14) 
where the index "f" indicates that the matrix and vector used are those of the crack element.
The expression Equation (14) shows that only the crack element is concerned, and consequently it is enough to place in the mesh another RMQ7 element equivalent to that placed on the crack, in other words an element which has the same geometry and made up of same material. The energy release rate is calculated according to the relation Equation (14) with only one discretization starting from the difference of the elementary matrices of the element containing the crack and representing the state and its equivalent element representing the state "a". The expression Equation (14) can be written differently as follows

(15) 
In practice, we carry out the discretization of the cracked structure in the configuration , and we locate the element containing the crack like its equivalent element representing the configuration "a", in order to save their elementary matrices during the assembly operation and before the application of the boundary conditions.
After the resolution phase, the nodal values of the crack element are extracted, and a special module is used to evaluate the energy release rate according to the following formula

(16) 
4.1 Presentation of the example
The example treated, in this study, is the interfacial crack centered of a bimaterial plate. This example was studied by Chow et al. [30] with plane strain condition. This rectangular bimaterial is made of material #1 and #2 and subjected to a tension σ_{22}^{0} = 1 MPa. As shown in Figure 2, the dimensions of the bimaterial are the half crack length a=1mm, the width w=20a and the height h=20a. Two cases are treated in this example. In the first case it is assumed that the materials #1 and #2 are isotropic and in the second case the materials are considered to be orthotropic (carbon composites: AS4/35016) with layup angle of 0 and 90 degree. The material properties of the used materials are defined in Table 1.
Figure 2. Bimaterial plate

A stress σ^{0}_{11 }is applied to the side of the material #2. In the case of plane strain, this stress is expressed by

(17) 
where E is the Young's modulus and ν is the Poisson's ratio of the material.
Isotropic  Orthotropic (0 degree)  Orthotropic (90 degree) 
G_{#1} = 1GPa
ν_{#1} = ν_{#2} = 0.3 
E_{3} = 142 GPa
E_{1}/E_{3} =E_{2}/E_{3} = 6.91x10^{2} G_{12}/E_{3} = 2.68x10^{2} G_{13}/E_{3}=G_{23}/E_{3}=4.23x10^{2} ν_{31} = ν_{32} = ν_{12} = 0.3 
E_{1} = 142 GPa
E_{2}/E_{1} =E_{3}/E_{1} = 6.91x10^{2} G_{23}/E_{1} = 2.68x10^{2} G_{13}/E_{1}=G_{12}/E_{1}=4.23x10^{2} ν_{12} = ν_{13} = ν_{23} = 0.3 
In the example above, the authors (Chow et al. 1995) calculate and compare the stress intensity factors K_{1} and K_{2}, the energy release rate is calculated according to K_{1} and K_{2} by the expression given by Qu and Bassani [31]. The results are resumed in Table 2 for the two materials (isotropic and orthotropic).
Material  Exact solution  Hybrid element  Mutual integral  Extrapolation technique  
205 nodes  679 nodes  237 nodes  679 nodes  237 nodes  
Isotropic  G_{#2}/G_{#1}=1  10,988E04  11,290E04  11,302 E04  11,253 E04  13,132 E04  12,554E04 
G_{#2}/G_{#1}=5  06,453E04  06,606E04  06,614 E04  06,592 E04  07,649 E04  07,326E04  
G_{#2}/G_{#1}=50  05,353E04  05,460E04  05,461 E04  05,444 E04  06,287 E04  06,026E04  
Orthotropic

[0/0]  03,170E04  03,257E04  03,262 E04  03,247 E04  03,793 E04  03,540E04 
[90/90]  02,200E04  02,221E04  02,216 E04  02,221 E04  02,549 E04  02,480E04  
[0/90]  02,640E04  02,685E04  02,679 E04  02,675 E04  03,094 E04  03,021E04 
4.2 Results and discussions
The mixed finite element RMQ7 is now used to calculate the energy release rate of the cracked bimaterial plate. For this purpose three meshes (207, 237 and 677 nodes) are used in order to be able to compare the results of RMQ7 element with the other elements results by using approximately the same number of nodes. The results obtained are resumed in the Table 3.
Material  RMQ7 mixed finite element  
207 nodes  237 nodes  677 nodes  
Isotropic  G_{#2}/G_{#1}=1  11,272E04  11,205E04  11,126E04 
G_{#2}/G_{#1}=5  06,393E04  06,486E04  06,438E04  
G_{#2}/G_{#1}=50  05,274E04  05,278E04  05,297E04  
Orthotropic  [0/0]  03,225E04  03,237E04  03,167 E04 
[90/90]  02,260E04  02,293E04  02,168 E04  
[0/90]  02,691E04  02,764E04  02,617 E04 
According to the number of nodes, the numerical results of the energy release rate for different methods are listed in Table 4, 5 and 6 for both the isotropic bimaterial and anisotropic bimaterial. The difference with exact solution for the different methods are calculated and consigned in Tables 4, 5 and 6. This difference is expressed by the Error (%) calculated as follows

(18) 
Compared to the exact solution, the numerical results show the accuracy and efficiency of the RMQ7 element. The difference between the values of exact solution and those of the mixed finite element vary between 0,10% and 4,70%.
Material  Exact solution  RMQ7 element  Hybrid element  
207 nodes  Error %  205 nodes  Error %  
Isotropic  G_{#2}/G_{#1}=1  10,988E04  11,272E04  2,58  11,290E04  2,75 
G_{#2}/G_{#1}=5  06,453E04  06,393E04  0,93  06,606E04  2,37  
G_{#2}/G_{#1}=50  05,353E04  05,274E04  1,48  05,460E04  2,00  
Orthotropic

[0/0]  03,170E04  03,225E04  1,74  03,257E04  2,74 
[90/90]  02,200E04  02,260E04  2,73  02,221E04  0,95  
[0/90]  02,640E04  02,691E04  1,93  02,685E04  1,70 
For isotropic bimaterials, the RMQ7 element, for the same number of nodes, shows a clear superiority compared to the eight noded isoparametric displacement finite element (extrapolation technique), and more accurate results than those of the mutual integral method. For example, with the RMQ7 element, the Error passed from 0,93% to 2,58% with 207 nodes whereas the Error varied from 2,00% to 2,75 using the hybrid element with 205 nodes. For orthotropic bimaterials, the element RMQ7 shows its performance compared to the classical displacement element. It still gives results clearly closer to the exact solution. Compared to the mutual integral method the RMQ7 element gives very satisfactory results. Using RMQ7 element with 677, the difference varied between 0,10% and 1,45% whereas it is between 0,73% and 2,90% using mutual method with 679 nodes.
Material  Exact solution  RMQ7 element  Mutual integral  Extrapolation technique  
237 nodes  Error %  237 nodes  Error %  237 nodes  Error %  
Isotropic  G_{#2}/G_{#1}=1  10,988E04  11,205E04  1,98  11,253E04  2,41  12,554E04  14,25  
G_{#2}/G_{#1}=5  06,453E04  06,486E04  0,51  06,592E04  2,15  07,326E04  13,53  
G_{#2}/G_{#1}=50  05,353E04  05,278E04  1,40  05,444E04  1,70  06,026E04  12,57  
Orthotropic

[0/0]  03,170E04  03,237E04  2,11  03,247E04  2,43  03,540E04  11,67  
[90/90]  02,200E04  02,293E04  4,23  02,221E04  0,95  02,480E04  12,73  
[0/90]  02,640E04  02,764E04  4,70  02,675E04  1,33  03,021E04  14,43 
The results obtained, using the present mixed finite element, show the efficiency and accuracy of the proposed numerical model, which can give an acceptable solution with a few degrees of freedom from a unique mesh. It should be noted that during numerical calculation, the choice of the variation of the crack length is very significant. Indeed, it is necessary that this variation is sufficiently small so that the approximation Equation (9) has a justification, and not too small to avoid problems involved in the precision machine.
The results show also, that the current techniques of the finite elements analysis make it possible to find an effective numerical solution and a high precision to the problems of fracture mechanic.
Material  Exact solution  RMQ7 element  Mutual integral  Extrapolation technique  
677 nodes  Error %  679 nodes  Error %  679 nodes  Error %  
Isotropic  G_{#2}/G_{#1}=1  10,988E04  11,126E04  1,26  11,302E04  2,86  13,132E04  19,51 
G_{#2}/G_{#1}=5  06,453E04  06,438E04  0,23  06,614E04  2,49  07,649E04  18,53  
G_{#2}/G_{#1}=50  05,353E04  05,297E04  1,05  05,461E04  2,02  06,287E04  17,45  
Orthotropic

[0/0]  03,170E04  03,167 E04  0,10  03,262E04  2,90  03,793E04  19,65 
[90/90]  02,200E04  02,168 E04  1,45  02,216E04  0,73  02,549E04  15,86  
[0/90]  02,640E04  02,617 E04  0,87  02,679E04  1,48  03,094E04  17,20 
In this paper, a numerical modeling has been proposed to study the interfacial crack between two orthotropic materials. This model uses a special mixed finite element developed in a natural plane. It is a twodimensional element with seven nodes: five displacement nodes and two stress nodes.
The proposed numerical model was used to calculate the energy release rate in a cracked orthotropic bimaterial using a technique that combines the present element with the virtual crack extension method. Two cases were treated: isotropic and orthotropic bimaterials.
The accuracy and the efficiency of the proposed model has been evaluated by comparing the numerical solution with an available analytical solution or numerical ones obtained from others methods. Comparisons with existing analytical or numerical solutions for interfacial cracks in orthotropic bimaterials proved that the proposed model provide a good accuracy and efficiency.
[1] Williams M. L. The stresses around a fault or crack in dissimilar media. Bulletin of seismology society of America, 49: 199204, 1959.
[2] Erdogan F. Stress distribution in nonhomogeneous elastic plane with cracks. J. appl. Mech., 30: 232236, 1963.
[3] Erdogan F. Stress distribution in bonded dissimilar materials with cracks. J. appl. Mech., 32: 403409, 1965.
[4] England A. H. A crack between dissimilar media. J. appl. Mech., 32: 400402, 1965.
[5] Rice J. R., Sih G. C. Plane problems of cracks in dissimilar media. J. appl. Mech., 32: 418423, 1965.
[6] Hutchinson J. W., Mear M. , Rice J. R. Crack Paralleling an Interface Between Dissimilar Materials. ASME Journal of Applied Mechanics, 54: 828832, 1987.
[7] Rice J. R. Elastic Fracture Mechanics Concepts for Interfacial Cracks. ASME Journal of Applied Mechanics, 55: 98103,1988.
[8] Suo Z., Hutchinson J.W. Interface crack between two elastic layers. Int J Fract., 43: 1–18, 1990.
[9] Gotoh M. Some problems of bonded anisotropic plates with cracks along the bond. Int. J. Fract. Mech., 3: 253265, 1967.
[10] Clements D.L. A crack between dissimilar anisotropic media. Int. J. Engng. Sci., 9: 257–265, 1971.
[11] Willis J. R. Fracture mechanics of interfacial cracks. J. Mech. Phys. Solids, 19: 353368, 1971.
[12] Qu J., Bassani J. L. Cracks on bimaterial and bicrystal interfaces. J. Mech. Phys. Solids, 37(4): 417433, 1989.
[13] Suo Z. Singularities, interfaces and cracks in dissimilar anisotropic media. Proc. R. Soc. Lond., A 427: 331358, 1990.
[14] Ni L., NematNasser S. Interface crack in anisotropic dissimilar materials: An analytic solution. J. Mech. Phys. Solids, 39(1): 113144, 1991.
[15] Bassani J. L., Qu J. Finite crack on bimaterial and bicrystal interfaces. J. Mech. Phys. Solids, 37(4): 435453, 1989.
[16] Qu J., Li Q. Interfacial dislocation and its applications to interface cracks in anisotropic bimaterials. J. Elasticity, 26: 169195, 1991.
[17] Gao Y., Liu Z., Zeng Q., Wang T., Zhuang Z., Hwang KC. Theoretical and numerical prediction of crack path in the material with anisotropic fracture toughness. Engineering Fracture Mechanics, 180: 330347, 2017.
[18] Wang X., Zhou, K., Wu M.S. Interface cracks with surface elasticity in anisotropic bimaterials. Int. J. of Solids and Structures, 59: 110120, 2015.
[19] Juan P.A., Dingreville R. Mechanics of finite cracks in dissimilar anisotropic elastic media considering interfacial elasticity. J. of the Mechanics and Physics of Solids, 99: 118, 2017.
[20] Tanaka K., Oharada K., Yamada D., Shimizu K. Fatigue crack propagation in shortcarbonfiber reinforced plastics evaluated based on anisotropic fracture mechanics. Int. Journal of Fatigue, 92: 415425, 2016.
[21] BanksSills L. and Ikeda T. Stress intensity factors for interface cracks between orthotropic and monoclinic material. Int. J. Fract., 167(1): 4756, 2011.
[22] Bouchemella S., Bouzerd, H., Khaldi, N. Modélisation des interfaces fissurées des bimatériaux orthotropes. XIXème Congrès Français de Mécanique, Marseille, France, 2009.
[23] Ghorashi S. Sh., Valizadeh N., Mohammadi S. Extended isogeometric analysis (XIGA) for simulation of stationary and propagating cracks. Int. J. Numer. Methods Eng., 89: 1069 –1101, 2012.
[24] Ghorashi S. Sh., Valizadeh N., Mohammadi S. , Rabczuk T. Tspline based XIGA for fracture analysis of orthotropic media. Computers and Structures, 147: 138–146, 2015.
[25] Habib S. H., Belaidi I., Khatir S., Abdel Wahab M. Numerical Simulation of cracked orthotropic materials using extended isogeometric analysis. Journal of Physics: Conf. Series 842 012061, 2017.
[26] Khatir S., Wahab M. A. Fast simulations for solving fracture mechanics inverse problems using PODRBF XIGA and Jaya algorithm. Engineering Fracture Mechanics, 205: 285300, 2019.
[27] Khatir S., Wahab M. A., Benaissa B., Köppen M. Crack identification using eXtended IsoGeometric Analysis and particle swarm optimization. In Fracture, Fatigue and Wear, pp. 210222, Springer, Singapore, 2018.
[28] Bouzerd H. Elément fini mixte pour interface cohérente ou fissurée. Thèse de doctorat, Université de Claude Bernard (Lyon I), France, 1992.
[29] Bouziane S., Bouzerd H., Guenfoud M. Mixed finite element for modelling interfaces. European Journal of Computational Mechanics, 18(2): 155175, 2009.
[30] Chow W. T., Beom H. G., Alturi S. N. Calculation of stress intensity factors for an interfacial crack between dissimilar anisotropic media, using a hybrid element method and mutual integral: Computational Mechanics, 15(6): 546557, 1995.
[31] Qu J., Bassani J. L. Interfacial fracture mechanics for anisotropic bimaterial. Journal of Applied Mechanics, 60: 422431, 1993.
Published on 09/08/20
Licence: CC BYNCSA license
Are you one of the authors of this document?