Computation of energy release rate for interfacial crack in orthotropic bimaterials

<!-- metadata commented in wiki content


==Computation of energy release rate for interfacial crack in orthotropic bimaterials==

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
'''Salah Bouziane<sup>a,b</sup>, Hamoudi Bouzerd<sup>a,b</sup>'''</div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<sup>a</sup>Civil Engineering department, University of August 20, 1955, Skikda, Algeria. E-mails: [mailto:bouziane_21@yahoo.fr bouziane_21@yahoo.fr], [mailto:bib_ham@yahoo.fr bib_ham@yahoo.fr]</div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<sup>b</sup>Laboratory of Civil Engineering and Hydraulics, University of May 8, 1945, Guelma, Algeria.</div>
-->
==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

==1. Introduction==

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 Nemat-Nasser [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 Banks-Sills 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 T-spline basis functions [24]. The same method XIGA was used by Habib et al. [25] 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 two-dimensional 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 two-dimensional 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.

==2. Numerical modelling of interfacial crack==

The bimaterial has been discredized using a special mixed finite element RMQ-7 ('''R'''eissner '''M'''odified '''Q'''uadrilateral) as shown in [[#img-1|Figure 1]](a). The present mixed finite element used in this study is two-dimensional element with seven nodes: five displacement nodes and two stress nodes as shown in [[#img-1|Figure 1]](b). The node 5 coincides with the crack tip. This element was developed by Bouzerd [28], in the physical (<math> x,y </math>) plane, and was reformulated and validated by Bouziane et al. [29] in a natural (<math>\xi , \eta  </math>) plane.

<div id='img-1'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
|-
|  style="text-align: center;padding:10px;"| [[Image:Draft_Bouziane_473899494-image1.png|378px]] 
|  style="text-align: center;padding:10px;"| [[Image:Draft_Bouziane_473899494-image2.png|378px]]
|-
|  style="text-align: center;font-size: 75%;"|(a) Discretization of bimaterial
|  style="text-align: center;font-size: 75%;"|(b) RMQ-7 element
|- style="text-align: center; font-size: 75%;"
| colspan="2" style="padding:10px;"| '''Figure 1'''. Discretization of bimaterial and RMQ-7 element
|}


Displacement for the present mixed finite element can be given by

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>u=\sum _{i=1}^{5}{N}_{i}{u}_{i}</math>
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(1)
|}

where  <math display="inline">{N}_{i}</math> are the shape functions and  <math display="inline">{u}_{i}</math> is the nodal displacement corresponding to node <math> i </math>. For the present element, the shape functions are given as follows

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{N}_{1}=-\frac{1}{4}\left( 1- \xi \right) \left( 1-\eta \right) \xi</math>, <math>\quad {N}_{2}=\frac{1}{4}\left( 1+\xi \right) \left( 1-\eta \right) \xi</math>, <math>\quad {N}_{3}=\frac{1}{4}\left( 1+\xi \right) \left( 1+\eta \right)</math>, 
|-
|<math>\quad {N}_{4}=frac{1}{4}\left( 1-\xi \right) \left( 1+\eta \right)</math>,  <math>\quad {N}_{5}=</math><math>\frac{1}{2}\left( 1-{\xi }^{2}\right) \left( 1-\eta \right)</math>
|}
|  style="text-align: right;"|(2)
|}

The element stress component is approximated by

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math display="inline">\left\{ \sigma \right\} =\left[ M\right] \left\{ \tau \right\}</math> 
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(3)
|}

where <math>[M]</math> is the matrix of interpolation functions for stresses and <math>\{\tau\}</math> is the vector of nodal stresses.

For the RMQ-7 element ([[#img-1|Figure 1]](b)), the shape functions <math>M_{i2}</math>, used to evaluate <math>\sigma_{12}</math> and <math>\sigma_{22}</math> [29] for nodes 6 and 7 are obtained by

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math display="inline">{M}_{i2}^{6}=\frac{1}{6}\left( 1-2\xi \right) \left( 1-\right. </math><math>\left. 2\eta \right)</math>,   <math>\quad {M}_{i2}^{7}=</math><math>\frac{1}{6}\left( 1+2\xi \right) \left( 1-2\eta \right) \quad ,\quad i=1,2</math>
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(4)
|}

The element stiffness  matrix ''[K<sub>e</sub>]'' is given by the following expression 

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math display="inline">\left[ {K}_{e}\right] =\left[ \begin{matrix}\left[ {K}_{\sigma \sigma }\right] &\left[ {K}_{\sigma u}\right] \\{\left[ {K}_{\sigma u}\right] }^{T}&\left[ 0\right] \end{matrix}\right]</math> '' ''
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(5)
|}

where the sub-matrices,  <math display="inline">\left[ {K}_{\sigma \sigma }\right]</math>  and  <math display="inline">\left[ {K}_{\sigma u}\right]</math> , are given by the following relations

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\left[ {K}_{\sigma \sigma }\right] =-t\int_{{A}_{e}}^{}{\left[ M\right] }^{T}\left[ S\right] \left[ M\right] d{A}^{e}</math> <math>\, \left[ {K}_{\sigma u}\right] =t\int_{{A}_{e}}^{}{\left[ M\right] }^{T}\left[ B\right] d{A}^{e}</math>
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(6)
|}

where <math>[S]</math> is the compliance matrix, <math>[M]</math> is the matrix of interpolation functions for stresses, <math>[B]</math> is the strain-displacement matrix of shape function derivatives, <math>t</math> is the thickness, <math>A^e</math> is the element area and <math>T</math> indicate the matrix transpose.

==3. Computation of energy release rate==

The virtual crack extension method, associated with the mixed finite element RMQ-7, is used to calculate the energy release rate <math>G</math> [28]. In this technique, the first calculation of the deformation energy   <math display="inline">{\Pi }_{1}</math> is carried out in the initial configuration of the crack. The crack is then moved an infinitesimal distance  <math display="inline">\delta a</math> in the direction of its axis. The deformation energy   <math display="inline">{\Pi }_{2}</math> is evaluated again in the second configuration, the energy released during this crack length variation is

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\delta \Pi ={\Pi }_{2}-\, {\Pi }_{1}</math>
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(7)
|}

The energy release rate <math>G</math> will be evaluated thereafter starting from the relation

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>G=\frac{\delta \Pi }{\delta a}</math>
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(8)
|}

Calculation by the virtual crack extension method requires two finite element analysis. The use of the RMQ-7 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 RMQ-7 element is associated to crack tip, and consequently the length of crack <math>a</math> can be increased by a quantity  <math display="inline">\delta a</math> 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   <math display="inline">u(a)</math> and   <math display="inline">u(a+\delta a)</math> obtained in the structure with a crack length <math>a</math> and in the same structure with a crack length  <math display="inline">a+\delta a</math> are as close as the disturbance  <math display="inline">\delta a</math> is small compared to dimensions of the crack element. We can thus write with a good approximation

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>u\left( a\right) =u(a+\delta a)</math>
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(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  <math display="inline">\delta a</math>, the energy release rate is calculated as follows:

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>G=-\frac{\Pi \left( a+\delta a\right) -\Pi (a)}{\delta a}</math>
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(10)
|}

where <math display="inline">\Pi \left( a+\right. </math><math>\left. \delta a\right)</math>  and   <math display="inline">\Pi (a)</math> represent respectively the deformation energy of the cracked structure in the configuration   <math display="inline">a+\delta a</math> and <math>a</math>.  

In its discretized form, the deformation energy is written

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\Pi =\frac{1}{2}\sum _{i=1}^{ne}{\left\{ u\right\} }_{i}^{T}{\left[ K\right] }_{i}{\left\{ u\right\} }_{i}</math>
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(11)
|}

where <math>ne</math> is the total number of elements in discretized structure, <math display="inline">{\left\{ u\right\} }_{i}</math> the vertical vector containing the nodal values of element <math>i</math>, <math display="inline">{\left[ K\right] }_{i}</math> the elementary matrix of element <math>i</math>, and the exponent <math>T</math> indicates the transposed vector.

By substitution of Equation (11) in Equation (10), the expression of the energy release rate <math>G</math> becomes

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>G=-\frac{1}{2\delta a}\left[ \sum _{i=1}^{ne}{\left\{ u(a+\delta a)\right\} }_{i}^{T}{\left[ K(a+\delta a)\right] }_{i}{\left\{ u(a+\delta a)\right\} }_{i}-\right. </math><math>\left. \sum _{i=1}^{ne}{\left\{ u(a)\right\} }_{i}^{T}{\left[ K(a)\right] }_{i}{\left\{ u(a)\right\} }_{i}\right]</math> 
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(12)
|}

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

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>G=-\frac{1}{2\delta a}\sum _{i=1}^{ne}{\left\{ u(a+\delta a)\right\} }_{i}^{T}\left[ {\left[ K(a+\delta a)\right] }_{i}-{\left[ K(a)\right] }_{i}\right] {\left\{ u(a+\delta a)\right\} }_{i}</math>
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(13)
|}

and as only the crack element is disturbed, then <math> G </math> results more simply in the relation

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>G=-\frac{1}{2\delta a}{\left\{ u(a+\delta a)\right\} }_{f}^{T}\left[ {\left[ K(a+\delta a)\right] }_{f}-\right. </math><math>\left. {\left[ K(a)\right] }_{fi}\right] {\left\{ u(a+\delta a)\right\} }_{f}</math>
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(14)
|}

where the index <math> f </math> 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 RMQ-7 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 <math display="inline">a+\delta a</math>  and its equivalent element representing the state <math>a</math>.  The expression Equation (14) can be written differently as follows

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{G=-\frac{1}{2}\left\{ u\right\} }_{f}^{T}\frac{{\left[ \delta K\right] }_{f}}{\delta a}{\left\{ u\right\} }_{f}</math>
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(15)
|}

In practice, we carry out the discretization of the cracked structure in the configuration  <math display="inline">a+\Delta a</math>, and we locate the element containing the crack like its equivalent element representing the configuration <math>a</math>, 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

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{G=-\frac{1}{2}\left\{ u\right\} }_{f}^{T}\frac{{\left[ \Delta K\right] }_{f}}{\Delta a}{\left\{ u\right\} }_{f}</math>
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(16)
|}

==4. Numerical examples ==

===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 <math>\sigma_{22}^0 = 1</math> MPa. As shown in [[#img-2|Figure 2]], the dimensions of the bimaterial are the half crack length a=1mm, the width w=20a and the height <math>h=20</math>a. 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/3501-6) with lay-up angle of 0 and 90 degree. The material properties of the used materials are defined in [[#tab-1|Table 1]].

<div id='img-2'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
|-
|style="padding:10px;"|  [[Image:Draft_Bouziane_473899494-image3.png|198px]] 
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding:10px;"| '''Figure 2'''. Bimaterial plate
|}


A stress <math>\sigma_{11}^0</math> is applied to the side of the material #2. In the case of plane strain, this stress is expressed by

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\sigma }_{11}^{0}=\left[ \frac{{\nu }_{12\#2}+{\nu }_{13\#2}{\nu }_{32\#2}}{1+{\nu }_{13\#2}{\nu }_{31\#2}}-\right. </math><math>\left. \frac{{\nu }_{12\#1}+{\nu }_{13\#1}{\nu }_{32\#1}}{1+{\nu }_{13\#2}{\nu }_{31\#2}}\left( \frac{{E}_{1\#2}}{{E}_{1\#1}}\right) \right] {\sigma }_{22}^{0}</math>
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(17)
|}

where <math>E</math> is the Young's modulus and <math>\nu</math> is the Poisson's ratio of the material. 

<div class="center" style="font-size: 75%;">'''Table 1'''. Material property</div>

<div id='tab-1'></div>
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;font-size:85%;width:auto;" 
|-style="text-align:center"
!'''Isotropic''' !! '''Orthotropic (0 degree)''' !! '''Orthotropic (90 degree)'''
|-style="text-align:center"
|  <math>G_{\#1}=1</math> GPa
|  <math>E_3 = 142</math> GPa
|  <math>E_1 = 142</math> GPa
|-style="text-align:center"
|<math>\nu_{\#1}=\nu_{\#2}=0.3</math>
|<math>E_1/E_3 =E_2/E_3 = 6.91\times 10^{-2}</math>
|<math>E_2/E_1 =E_3/E_1 = 6.91\times 10^{-2}</math>
|-style="text-align:center"
| 
|<math>G_{12}/E_3 = 2.68\times 10^{-2}</math>
|<math>G_{23}/E_1 = 2.68\times 10^{-2}</math>
|-style="text-align:center"
|
|<math>G_{13}/E_3=G_{23}/E_3=4.23\times 10^{-2}</math>
|<math>G_{13}/E_1=G_{12}/E_1=4.23\times 10^{-2}</math>
|-style="text-align:center"
|
|<math>\nu_{31} = \nu_{32} = \nu_{12} = 0.3</math>
|<math>\nu_{12} = \nu_{13} = \nu_{23} = 0.3</math>
|}


In the example above, the authors (Chow et al. 1995) calculate and compare the stress intensity factors <math>K_1</math> and <math>K_2</math>, the energy release rate is calculated according to <math>K_1</math> and <math>K_2</math> by the expression given by Qu and Bassani [31]. The results are resumed in [[#tab-2|Table 2]] for the two materials (isotropic and orthotropic).

<div class="center" style="font-size: 75%;">'''Table 2.''' Energy release rate in the numerical example</div>

<div id='tab-2'></div>
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;font-size:85%;width:auto;" 
|-style="text-align:center"
!  rowspan='2' colspan='2' | Material !! rowspan='2' |Exact solution !! Hybrid element !!  colspan='2'  |Mutual integral !!  colspan='2'  |Extrapolation technique
|-style="text-align:center"
!  205 nodes !! 679 nodes !! 237 nodes !! 679 nodes !! 237 nodes
|-style="text-align:center"
|  rowspan='3' style="text-align:left;" |Isotropic
|  <math>G_{\#2}/G_{\#1}=1</math>
|  10,988E-04
|  11,290E-04
|  11,302 E-04
|  11,253 E-04
|  13,132 E-04
|  12,554E-04
|-style="text-align:center"
| <math>G_{\#2}/G_{\#1}=5</math> 
|  06,453E-04
|  06,606E-04
|  06,614 E-04
|  06,592 E-04
|  07,649 E-04
|  07,326E-04
|-style="text-align:center"
|  <math>G_{\#2}/G_{\#1}=50</math> 
|  05,353E-04
|  05,460E-04
|  05,461 E-04
|  05,444 E-04
|  06,287 E-04
|  06,026E-04
|-style="text-align:center"
|  rowspan='3' |Orthotropic
|  [0/0]
|  03,170E-04
|  03,257E-04
|  03,262 E-04
|  03,247 E-04
|  03,793 E-04
|  03,540E-04
|-style="text-align:center"
| [90/90]
|  02,200E-04
|  02,221E-04
|  02,216 E-04
|  02,221 E-04
|  02,549 E-04
|  02,480E-04
|-style="text-align:center"
|  [0/90]
|  02,640E-04
|  02,685E-04
|  02,679 E-04
|  02,675 E-04
|  03,094 E-04
|  03,021E-04
|}

===4.2 Results and discussions===

The mixed finite element RMQ-7 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 RMQ-7 element with the other elements results by using approximately the same number of nodes. The results obtained are resumed in the [[#tab-3|Table 3]].

<div class="center" style="font-size: 75%;">'''Table 3'''. Energy release rate obtained using RMQ-7 element</div>

<div id='tab-3'></div>
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;font-size:85%;width:auto;" 
|-style="text-align:center"
!  rowspan='2' colspan='2' | Material !! colspan='3'  |'''RMQ-7 mixed finite element'''
|-style="text-align:center"
!  207 nodes !! 237 nodes !! 677 nodes
|-style="text-align:center"
|  rowspan='3' style="text-align:left;"|Isotropic
|  <math>G_{\#2}/G_{\#1}=1</math>
|  11,272E-04
|  11,205E-04
|  11,126E-04
|-style="text-align:center"
| <math>G_{\#2}/G_{\#1}=5</math>
|  06,393E-04
|  06,486E-04
|  06,438E-04
|-style="text-align:center"
|  <math>G_{\#2}/G_{\#1}=50</math>
|  05,274E-04
|  05,278E-04
|  05,297E-04
|-style="text-align:center"
|  rowspan='3' style="text-align:left;"|Orthotropic
|  [0/0]
|  03,225E-04
|  03,237E-04
|  03,167 E-04
|-style="text-align:center"
| [90/90]
|  02,260E-04
|  02,293E-04
|  02,168 E-04
|-style="text-align:center"
|  [0/90]
|  02,691E-04
|  02,764E-04
|  02,617 E-04
|}


According to the number of nodes, the numerical results of the energy release rate for different methods are listed in [[#tab-4|Tables 4]], [[#tab-5|5]] and [[#tab-6|6]] for both the isotropic bimaterial and anisotropic bimaterial. The difference with exact solution for the different methods are calculated and consigned in [[#tab-4|Tables 4]], [[#tab-5|5]] and [[#tab-6|6]]. This difference is expressed by the Error (%) calculated as follows

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>Error\, (\%)=\frac{G-{G}_{exact}}{{G}_{exact}}\times 100\, \%</math>
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(18)
|}


Compared to the exact solution, the numerical results show the accuracy and efficiency of the RMQ-7 element. The difference between the values of exact solution and those of the mixed finite element vary between -0,10% and 4,70%.

<div class="center" style="font-size: 75%;">'''Table 4'''.  Energy release rate for crack along bimaterial interface, Mesh 1: 207 nodes</div>

<div id='tab-4'></div>
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;font-size:85%;width:auto;" 
|-style="text-align:center"
!  rowspan='2' colspan='2' | Material !! rowspan='2' |Exact solution !! colspan='2'  |'''RMQ-7 element'''!!  colspan='2'  |Hybrid element
|-style="text-align:center"
! 207 nodes !!  Error % !!   205 nodes !!   Error %
|-style="text-align:center"
|  rowspan='3' style="text-align:left;"|Isotropic
|  <math>G_{\#2}/G_{\#1}=1</math>
|  10,988E-04
|  11,272E-04
|  2,58
|  11,290E-04
|  2,75
|-style="text-align:center"
|  <math>G_{\#2}/G_{\#1}=5</math>
|  06,453E-04
|  06,393E-04
|  -0,93
|  06,606E-04
|  2,37
|-style="text-align:center"
|  <math>G_{\#2}/G_{\#1}=50</math>
|  05,353E-04
|  05,274E-04
|  -1,48
|  05,460E-04
|  2,00
|-style="text-align:center"
|  rowspan='3' style="text-align:left;"|Orthotropic
|  [0/0]
|  03,170E-04
|  03,225E-04
|  1,74
|  03,257E-04
|  2,74
|-style="text-align:center"
|  [90/90]
|  02,200E-04
|  02,260E-04
|  2,73
|  02,221E-04
|  0,95
|-style="text-align:center"
|  [0/90]
|  02,640E-04
|  02,691E-04
|  1,93
|  02,685E-04
|  1,70
|}


For isotropic bimaterials, the RMQ-7 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 RMQ-7 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 RMQ-7 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 RMQ-7 element gives very satisfactory results. Using RMQ-7 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.

<div class="center" style="font-size: 75%;">'''Table 5'''. Energy release rate for crack along bimaterial interface, Mesh 2: 237 nodes</div>

<div id='tab-5'></div>
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;font-size:85%;width:auto;" 
|-style="text-align:center"
!  rowspan='2' colspan='2' | Material !! Exact solution !! rowspan='2' |RMQ-7 element !!  colspan='2' | Mutual integral !!  colspan='2' | Extrapolation technique
|-style="text-align:center"
!  237 nodes !! Error % !! 237 nodes !! Error % !! 237 nodes !! Error %
|-style="text-align:center"
|  rowspan='3' style="text-align:left;"|Isotropic
|  <math>G_{\#2}/G_{\#1}=1</math>
|  10,988E-04
|  11,205E-04
|  1,98
|  11,253E-04
|  2,41
|  12,554E-04
|  14,25
|-style="text-align:center"
| <math>G_{\#2}/G_{\#1}=5</math>
|  style="vertical-align: top;"|06,453E-04
|  style="vertical-align: top;"|06,486E-04
|  style="vertical-align: top;"|0,51
|  colspan='2'  style="vertical-align: top;"|06,592E-04
|  style="vertical-align: top;"|2,15
|  style="vertical-align: top;"|07,326E-04
|  style="vertical-align: top;"|13,53
|-style="text-align:center"
|  <math>G_{\#2}/G_{\#1}=50</math>
|  style="border-bottom: 1pt solid black;vertical-align: top;"|05,353E-04
|  style="border-bottom: 1pt solid black;vertical-align: top;"|05,278E-04
|  style="border-bottom: 1pt solid black;vertical-align: top;"|-1,40
|  colspan='2'  style="border-bottom: 1pt solid black;vertical-align: top;"|05,444E-04
|  style="border-bottom: 1pt solid black;vertical-align: top;"|1,70
|  style="border-bottom: 1pt solid black;vertical-align: top;"|06,026E-04
|  style="border-bottom: 1pt solid black;vertical-align: top;"|12,57
|-style="text-align:center"
|  rowspan='3' |Orthotropic
|  style="border-top: 1pt solid black;vertical-align: top;"|[0/0]
|  style="border-top: 1pt solid black;vertical-align: top;"|03,170E-04
|  style="border-top: 1pt solid black;vertical-align: top;"|03,237E-04
|  style="border-top: 1pt solid black;vertical-align: top;"|2,11
|  colspan='2'  style="border-top: 1pt solid black;vertical-align: top;"|03,247E-04
|  style="border-top: 1pt solid black;vertical-align: top;"|2,43
|  style="border-top: 1pt solid black;vertical-align: top;"|03,540E-04
|  style="border-top: 1pt solid black;vertical-align: top;"|11,67
|-
|  style="vertical-align: top;"|[90/90]
|  style="vertical-align: top;"|02,200E-04
|  style="vertical-align: top;"|02,293E-04
|  style="vertical-align: top;"|4,23
|  colspan='2'  style="vertical-align: top;"|02,221E-04
|  style="vertical-align: top;"|0,95
|  style="vertical-align: top;"|02,480E-04
|  style="vertical-align: top;"|12,73
|-
|  style="border-bottom: 1pt solid black;vertical-align: top;"|[0/90]
|  style="border-bottom: 1pt solid black;vertical-align: top;"|02,640E-04
|  style="border-bottom: 1pt solid black;vertical-align: top;"|02,764E-04
|  style="border-bottom: 1pt solid black;vertical-align: top;"|4,70
|  colspan='2'  style="border-bottom: 1pt solid black;vertical-align: top;"|02,675E-04
|  style="border-bottom: 1pt solid black;vertical-align: top;"|1,33
|  style="border-bottom: 1pt solid black;vertical-align: top;"|03,021E-04
|  style="border-bottom: 1pt solid black;vertical-align: top;"|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  <math display="inline">\Delta a</math> 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.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
'''Table 6.''' Energy release rate for crack along bimaterial interface, Mesh 3: 677 nodesy</div>

{| style="width: 100%;border-collapse: collapse;" 
|-
|  rowspan='2' colspan='2'  style="border-top: 1pt solid black;vertical-align: top;"|'''Material'''
|  rowspan='2' style="border-top: 1pt solid black;vertical-align: top;"|'''Exact solution'''
|  colspan='2'  style="border-top: 1pt solid black;vertical-align: top;"|'''RMQ-7 element'''
|  colspan='2'  style="border-top: 1pt solid black;vertical-align: top;"|'''Mutual integral'''
|  colspan='2'  style="border-top: 1pt solid black;vertical-align: top;"|'''Extrapolation technique'''
|-
|  style="border-bottom: 1pt solid black;vertical-align: top;"|677 nodes
|  style="border-bottom: 1pt solid black;vertical-align: top;"|Error %
|  style="border-bottom: 1pt solid black;vertical-align: top;"|679 nodes
|  style="border-bottom: 1pt solid black;vertical-align: top;"|Error %
|  style="border-bottom: 1pt solid black;vertical-align: top;"|679 nodes
|  style="border-bottom: 1pt solid black;vertical-align: top;"|Error %
|-
|  rowspan='3' style="border-top: 1pt solid black;vertical-align: top;"|Isotropic
|  style="border-top: 1pt solid black;vertical-align: top;"|G<sub>#2</sub>/G<sub>#1</sub>=1
|  style="border-top: 1pt solid black;vertical-align: top;"|10,988E-04
|  style="border-top: 1pt solid black;vertical-align: top;"|11,126E-04
|  style="border-top: 1pt solid black;vertical-align: top;"|1,26
|  style="border-top: 1pt solid black;vertical-align: top;"|11,302E-04
|  style="border-top: 1pt solid black;vertical-align: top;"|2,86
|  style="border-top: 1pt solid black;vertical-align: top;"|13,132E-04
|  style="border-top: 1pt solid black;vertical-align: top;"|19,51
|-
|  style="vertical-align: top;"|G<sub>#2</sub>/G<sub>#1</sub>=5
|  style="vertical-align: top;"|06,453E-04
|  style="vertical-align: top;"|06,438E-04
|  style="vertical-align: top;"|-0,23
|  style="vertical-align: top;"|06,614E-04
|  style="vertical-align: top;"|2,49
|  style="vertical-align: top;"|07,649E-04
|  style="vertical-align: top;"|18,53
|-
|  style="border-bottom: 1pt solid black;vertical-align: top;"|G<sub>#2</sub>/G<sub>#1</sub>=50
|  style="border-bottom: 1pt solid black;vertical-align: top;"|05,353E-04
|  style="border-bottom: 1pt solid black;vertical-align: top;"|05,297E-04
|  style="border-bottom: 1pt solid black;vertical-align: top;"|-1,05
|  style="border-bottom: 1pt solid black;vertical-align: top;"|05,461E-04
|  style="border-bottom: 1pt solid black;vertical-align: top;"|2,02
|  style="border-bottom: 1pt solid black;vertical-align: top;"|06,287E-04
|  style="border-bottom: 1pt solid black;vertical-align: top;"|17,45
|-
|  rowspan='3' style="border-top: 1pt solid black;border-bottom: 1pt solid black;vertical-align: top;"|Orthotropic


|  style="border-top: 1pt solid black;vertical-align: top;"|[0/0]
|  style="border-top: 1pt solid black;vertical-align: top;"|03,170E-04
|  style="border-top: 1pt solid black;vertical-align: top;"|03,167 E-04
|  style="border-top: 1pt solid black;vertical-align: top;"|-0,10
|  style="border-top: 1pt solid black;vertical-align: top;"|03,262E-04
|  style="border-top: 1pt solid black;vertical-align: top;"|2,90
|  style="border-top: 1pt solid black;vertical-align: top;"|03,793E-04
|  style="border-top: 1pt solid black;vertical-align: top;"|19,65
|-
|  style="vertical-align: top;"|[90/90]
|  style="vertical-align: top;"|02,200E-04
|  style="vertical-align: top;"|02,168 E-04
|  style="vertical-align: top;"|-1,45
|  style="vertical-align: top;"|02,216E-04
|  style="vertical-align: top;"|0,73
|  style="vertical-align: top;"|02,549E-04
|  style="vertical-align: top;"|15,86
|-
|  style="border-bottom: 1pt solid black;vertical-align: top;"|[0/90]
|  style="border-bottom: 1pt solid black;vertical-align: top;"|02,640E-04
|  style="border-bottom: 1pt solid black;vertical-align: top;"|02,617 E-04
|  style="border-bottom: 1pt solid black;vertical-align: top;"|-0,87
|  style="border-bottom: 1pt solid black;vertical-align: top;"|02,679E-04
|  style="border-bottom: 1pt solid black;vertical-align: top;"|1,48
|  style="border-bottom: 1pt solid black;vertical-align: top;"|03,094E-04
|  style="border-bottom: 1pt solid black;vertical-align: top;"|17,20
|}

==5. Conclusion==

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 two-dimensional 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.

==References==

[1] Williams M. L. The stresses around a fault or crack in dissimilar media. Bulletin of seismology society of America, 49: 199-204, 1959.

[2] Erdogan F. Stress distribution in nonhomogeneous elastic plane with cracks. J. appl. Mech., 30: 232-236, 1963.

[3] Erdogan F. Stress distribution in bonded dissimilar materials with cracks.  J. appl. Mech., 32: 403-409, 1965.

[4] England A. H. A crack between dissimilar media. J. appl. Mech., 32: 400-402, 1965.

[5] Rice J. R., Sih G. C. Plane problems of cracks in dissimilar media.  J. appl. Mech., 32: 418-423, 1965.

[6] Hutchinson J. W., Mear M. , Rice J. R. Crack Paralleling an Interface Between Dissimilar Materials. ASME Journal of Applied Mechanics, 54: 828-832, 1987.

[7] Rice J. R. Elastic Fracture Mechanics Concepts for Interfacial Cracks. ASME Journal of Applied Mechanics, 55: 98-103,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: 253-265, 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: 353-368, 1971.

[12] Qu J., Bassani J. L. Cracks on bimaterial and bicrystal interfaces.  J. Mech. Phys. Solids, 37(4): 417-433, 1989.

[13] Suo  Z. Singularities, interfaces and cracks in dissimilar anisotropic media. Proc. R.  Soc. Lond., A 427: 331-358, 1990.

[14] Ni L.,  Nemat-Nasser S.  Interface crack in anisotropic dissimilar materials: An analytic solution. J. Mech. Phys. Solids, 39(1): 113-144, 1991.

[15] Bassani J. L., Qu J. Finite crack on bimaterial and bicrystal interfaces. J. Mech. Phys. Solids, 37(4): 435-453, 1989.

[16] Qu J., Li Q. Interfacial dislocation and its applications to interface cracks in anisotropic bimaterials. J. Elasticity, 26: 169-195, 1991.

[17] Gao Y., Liu Z., Zeng Q., Wang T., Zhuang  Z., Hwang K-C. Theoretical and numerical prediction of crack path in the material with anisotropic fracture toughness. Engineering Fracture Mechanics, 180: 330-347, 2017.

<span id='bau010'><span id='bau015'>[18] [https://www.sciencedirect.com/science/article/pii/S0020768315000347 Wang] X., [https://www.sciencedirect.com/science/article/pii/S0020768315000347 Zhou], K., [https://www.sciencedirect.com/science/article/pii/S0020768315000347 Wu] M.S. Interface cracks with surface elasticity in anisotropic bimaterials. [https://www.sciencedirect.com/science/journal/00207683 Int. J. of Solids and Structures], 59: 110-120, 2015.

<span id='bau0005'><span id='bau0010'>[19] Juan [https://www.sciencedirect.com/science/article/pii/S0022509616305610 P.-A.],  [https://www.sciencedirect.com/science/article/pii/S0022509616305610 Dingreville] R.  Mechanics of finite cracks in dissimilar anisotropic elastic media considering interfacial elasticity.  [https://www.sciencedirect.com/science/journal/00225096 J. of the Mechanics and Physics of Solids], [https://www.sciencedirect.com/science/journal/00225096/99/supp/C 99]: 1-18, 2017.

[20] Tanaka K., Oharada K., Yamada D., Shimizu K. Fatigue crack propagation in short-carbon-fiber reinforced plastics evaluated based on anisotropic fracture mechanics. Int. Journal of Fatigue, 92: 415-425, 2016.

[21] Banks-Sills L. and Ikeda T. Stress intensity factors for interface cracks between orthotropic and monoclinic material. Int. J. Fract., 167(1): 47-56, 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. T-spline 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 POD-RBF XIGA and Jaya algorithm. Engineering Fracture Mechanics, 205: 285-300, 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. 210-222, 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): 155-175, 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): 546-557, 1995.

[31] Qu J., Bassani  J. L. Interfacial fracture mechanics for anisotropic bimaterial.  Journal of Applied Mechanics, 60: 422-431, 1993.