Mixed Finite Element Computation of Energy Release Rate in Anisotropic Materials Based on Virtual Crack Closure-Integral Method

A BSTRACT . The material with anisotropic properties are becoming widely essential due to the ease to manipulate their mechanical properties to obtain a particular quality, insure safety or a specific behavior. Those kinds of materials are considered anisotropic because their characteristics and behavior are dependent on every direction of the material’s orientation. In this work, the virtual crack closure-integral technique is implemented to a mixed finite element, in addition to the stiffness derivative procedure, to evaluate the ERR of crack extension in anisotropic materials. A simulation of a cracked edge rectangular plate with anisotropic characteristics is taken for example. The results obtained are in good agreement with the analytical results, making the proposed technique a good model for fracture investigation and allow it to study more complicated cases in future works. are compared with the analytical ones giving for the stiffness derivative procedure as1.09% for 90°, 2.80% for 60°, 2.57% for 45°, 2.88% for 30°, 0.08% for 15° and 3.69% for 0°. And the error gap for the crack closure are 1.25% for 90°, 1.10% for 60°, 1.54% for 45°, 1.89% for 30°, 1.57% for 15° and 0.38% for 0°. the of of the are the closest to the analytical which proves their efficiency and


INTRODUCTION
n the study of fracture mechanics problems, the anisotropic media represent the most general material characteristics and the most hazardous ones, making them challenging to the procedures of studying the crack propagation problems. In nature, most materials are anisotropic and in the recent works on industrial fields tends to focus on those type of materials [1][2][3][4][5]. I Several works on the crack propagation in anisotropic material made their reliability, where it is investigated with the linear elastic fracture mechanics [6] and the J-Integral [7] but the phase-field modeling developed for the anisotropic surface energy [8] is getting popular for brittle anisotropic material [9,10] and in crystalline material [11,12]. The finite element method is also known for treating the fracture mechanics problems, the Stroh formalism is used with an enriched boundary element method to offer better results for any degrees of anisotropy [6]. The extended finite element method (XFEM) also was enriched to treat the anisotropy [13]. The recently developed edge-based smoothed finite element method ES-FEM [14] and the singular edge-based smoothed finite element method sES-FEM [15,16] present a significantly improved accuracy of the finite element method FEM [17,18] without changing much of the FEM settings. The rotating material axes investigate the anisotropy of material and help to approves the new finite element methods [19] such as the scaled boundary finite element method [20], fractal two-level finite element method [21], the element-free Galerkin method with the fractal finite element method [22] and the ES-FEM [16]. The deviation of the principal direction of a crack is defined as kinking of the crack which can be seen as the first step of the propagation, evaluated either with the ERR [23][24][25][26] or more often with the stress intensity factor [27]. It's taken into account for the fatigue crack growth [28] the effect of the T-stress [27] and different other cases. Among the diverse techniques for fracture investigation, the standard Irwin's crack closure integral is widely used for the computation of the ERR. Ronald made a historical approach and a review of the virtual crack closure technique and its parent technique the crack closure integral [29]. The study of the kinked crack phenomenon with the virtual crack closure technique gave also great results [30]. The same method is implemented on developed finite elements such as the cell-based smoothed finite element [31], the floating node method [32], or also for three-dimensional elements such as the tetrahedral finite element [33]. The main purpose of this current work is to solve fracture mechanics problems of anisotropic media with a special mixed finite element RMQ7 (Reissner's Modified Quadrilateral with 7-nodes). This element has been associated with the virtual crack closure-integral technique to evaluate the ERR for anisotropic materials. The present mixed finite element has the particularity to contain seven nodes, two stress nodes with two degrees of freedom (σ12,σ22) and five displacement nodes with also two degrees of freedom (u 1 ,u 2 ). The RMQ7 is created with the purpose to investigate mainly the fracture mechanics-related problems with enough efficiency and easier steps of the calculation, giving by that a significant gain of time. This element was, at first, created by Bouzerd [34] with a setup in a physical (x, y) plane. Bouziane et al. [35] redevelop the element beginning from the parent element in a natural (ξ, η) plane. That added a rationalization of the calculations and a huge benefit of modeling distinctive types of cracks with different orientations. The method gave great precision on the paperwork of Bouziane et al. [36], where a numerical example of a center cracked plan with uniaxial tension was investigated in two cases of homogenous materials and bi-materials. Also the study of Mohamed Ben Ali et al. [37] on sandwich structures with two cases of symmetrical double cantilever beam and asymmetrical double cantilever beam. For the work of Benmalek et al. [38], the exactitude of the results was highly remarkable.

REISSNER'S MODIFIED QUADRILATERAL ELEMENT
he RMQ7 is the final element proposed by Bouzerd [34] with 7 Nodes and 14 degrees of freedom ( Fig. 1.d). It is reduced from the RMQ11, an element with 11 nodes and 22 degrees of freedom ( Fig. 1.c), with a static condensation procedure, that merges the internal degrees of freedom by reducing the size of equations per elimination of a certain number of variables. Bouziane et al. [35] introduced the design of the finite element in a natural (ξ, η) plane. The RMQ11 itself is gotten from the parent element RMQ5 by re-localization of certain variables and by displacement of the static nodal unknown of the corners toward the side itself. In which, the RMQ5 is obtained by adding a displacement node to the Reissner mixed element to get an element with 5 nodes and 22 degrees of freedom ( Fig. 1.b). Finally, the Reissner element is a four-nodded element with five degrees of freedom at each ( Fig. 1.a); its formulation is based on Reissner's mixed variational principle Bouziane et al. [35]. Three of its sides are compatible with linear traditional elements and present a displacement node at each corner. The last side, furthermore on its two displacement nodes about corner (node 1 and node 2), offers three extra nodes: a median node (node 5) and two middle nodes in the medium on each half-side (nodes 6 and 7), presenting the components of the stress vector along with the interface. Bouziane et al. [35] have presented the formulation and the validation of the element. The element displacement component is estimated as: where:      The shape functions are: The stress field in any point is written where M is the matrix of interpolation functions for stresses and {τ} vector of nodal stresses. In the configuration of Fig.  1, the shape functions used to estimate are given by: The shape functions used to calculate and is given as follows: The nodal estimation of the displacement and stress fields is expressed by: where [B] is the strain-displacement transformation matrix.
The element matrix [K e ] is given by and [S] is the compliance matrix; [K σσ ] is a block on the reduced element stiffness matrix that depends on the stress interpolation functions only and [Kσu] is a block on the reduced element stiffness matrix that depends on the stress interpolation functions and the strain displacement transformation matrix. Otherwise, the relation to define the relationship between the forces and the displacement on the nodes is:

COMPUTATION OF ENERGY RELEASE RATE
he virtual crack closure-integral method, in addition to the stiffness derivative procedure, has been associated with the present mixed finite element to calculate the ERR for anisotropic materials. This section is reserved for the presentation of the methodology adopted for the calculation of the ERR.

Virtual crack closure-integral method (CCI)
The computation of the ERR can be evaluated with the crack closure integral technique [30,31,39] which, for a homogenous material, considers the strain energy released for a small crack extension is equal to the energy essential to close the crack. The ERR for this technique is expressed by Irwin [40] for Mode I, as: T where σy is the normal stress applied in the closed crack, Δl is the small crack extension. Δv is the displacement of two points on the crack lips along the y-axis (Fig. 2).
For Mode II: where Δu is the displacement of two points of the crack lips along the x-axis, and τ is the shear stress along the closed crack. The application of the Eqns. (12) and (13) can be executed by computing the stresses for closing the crack and calculate the displacements corresponding, which could be expressed by: (1) Considering the implementation of the two Eqns. (14) and (15) in the finite element RMQ7 proposed, it is mandatory to define the nodes that will help in getting the value of the ERR. First, the forces that are submitted on the nodes (5, 6, 5' and 7') which represent the first step of the crack closer integral technique and are determined by the Eqn. (11), then the second step is to determine the displacement on the nodes (4, 3', 7 and 6') ( Fig. 3)   Defining the forces and the displacement will allow the evaluation of the ERR for the modes I and II as: The addition gives the total ERR:

Stiffness derivative procedure (SDP)
The potential energy π of a cracked body containing a crack is given: The force F is due to external loads acting on the boundaries other than the crack faces, knowing that the body forces are absent and F does not vary with length proper to the cracked body; then, in conclusion, the ERR is obtained only by the derivative of the stiffness, as follow: Or as Bouziane et al. [36] states it: where: a : initial crack length.
nf: total number of elements concerned by the disturbance δa. {ν} f : vertical vector containing the nodal values of the element concerned by the disturbance.
The derivative δK is obtained from the two-state of the crack body (before and after crack extension), it is the difference between the stiffness of the elements around the crack tip in the configuration M 0 and the stiffness of the same elements of the configuration M 1 (Fig. 4).
The computation of the strain-ERR occurs with an extension crack-length δa very small and the nodal displacements are obtained from the configuration M 1 .

ELASTIC CONSTANTS FOR A NEW COORDINATE SYSTEM
or a given body, the elastic constants that are known in a certain system will be different in another system. Considering a body in a local Cartesian coordinate system (x,y) with its elastic constants known or defined; a different Cartesian coordinate system (x',y') would have a different value. In general, the principal elastic constant is given for an orthotropic media in the principal coordinate system and can be recalculated for another coordinate system. Such as, the two-dimensional strain-stress relation for an orthotropic media in the principal coordinate system is  where: Note that subscripts 1, 2, and 3 denote the components corresponding to the tension in the two main axes and the in-plane shear, respectively. The stresses for the coordinate system (x,y) and the stresses for the coordinate system (x',y') are related by the transformation matrix [T] such as: where the transformation matrix [T] is given: It is now possible to have the relation between the two compliance matrices of the two coordinate systems, such as: Moreover, the stress-strain relation of the same body in a coordinate system oriented at an angle φ is considered as anisotropic and defined as: here Ei',νi' and G'12 are the Young's Modulus, Poisson's coefficient, and shear modulus respectfully for the orientated axes plane, and η i ' is the coefficient of mutual influence of the first type [41][42] (i = 1; 2). These modules and coefficients for the new axes are well determined by reversing every element of the new stiffness matrix as follow: and so on for every elastic constant: 21 12 22 S S       (32) Note that both of the methods of Voyiadjis [41] and Lekhnitskii [43] were used to obtain the elastics constants for the oriented coordinate system, and the values were slightly different from the actual technique presented, which made a significant impact on the final results of the ERR.

NUMERICAL VALIDATION
n this part, a demonstration is required to approve the results proposed by the method mentioned above. A plate of 1mm width and 2 mm length with a central crack edge crack, is been submitted to a traction σ= 1.0 Pa pulling the plate upward while the plate is clamped at the bottom left corner and the displacement along the y axis is blocked for the rest of the bottom edge. The material of the plate is a graphite-epoxy composite (Fig. 6) with a Young's Modulus in the two main principal directions of the fiber E 1 = 144.8 GPa and E 2 = 11.7 GPa, Poisson's ratio of ν 12 =0.21and a shear modulus of G 12 = 9.66 GPa. The material is considered orthotropic when its orientation is at 0°. To get the anisotropic properties out from the same material, orientation with an angle φ varying from -90° to 90° changes its elastic constants giving E'1, E'2, ν'12, and G'12. (Tab. 1). Chen et al. [16] have presented different finite element method applications on this problem, in Tab. 2 FEM represents the Finite Element Method and ES FEM is the Edge Smoothed Finite Element Method. The analytical results are provided from the work of Bowie and Freese [44]and the results for the finite element methods proposed on the paper were compared to those last.

I
The medialization is made with an elaborated program in MATHLAB R2014b, the number of elements adopted is 450, and the number of nodes generated is 1201. The values of the ERR are calculated using the Eqn. (21) for the stiffness derivative procedure with the integration relative domain "δa/a" stable at 10 -6 to 10 -13 . Moreover, for the crack closure technique the Eqn. (17)    The results are compared with the analytical ones giving error gaps for the stiffness derivative procedure as1.09% for 90°, 2.80% for 60°, 2.57% for 45°, 2.88% for 30°, 0.08% for 15° and 3.69% for 0°. And the error gap for the crack closure integral are 1.25% for 90°, 1.10% for 60°, 1.54% for 45°, 1.89% for 30°, 1.57% for 15° and 0.38% for 0°. As the chart in (Fig. 7) shows, the results of both of the methods are the closest to the analytical results, which proves their efficiency and accuracy.

CONCLUSION
n this paper, a special mixed finite element RMQ7 is adopted to resolve discontinuity problems in anisotropic media. To simulate the anisotropy, an orthotropic media is oriented with an angle, which results in anisotropic elastic properties that are different from an oriented ax to another. In the computation of the ERR after a crack extension with the current element, two techniques were compared side by side. The crack closure integral gave closer values to the analytics results compared to the stiffness derivative method values, the implementation of those two techniques on the RMQ7 represent an accurate evaluation and computation efficiency of the ERR and outperform the standard FEM and the ES FEM. The RMQ7 present one mesh for the calculation of the ERR, which represent a considerable profit in computing times and set in data compared to the traditional methods, which use two meshes. The stiffness derivative procedure computation is dependent on the domain size, as it is considered between 10 -6 and 10 -13 of the crack length "δa/a".