Microsoft Word - numero_57_art_26_3096 S. Derouiche et alii, Frattura ed Integrità Strutturale, 57 (2021) 359-372; DOI: 10.3221/IGF-ESIS.57.26 359 Mixed Finite Element Computation of Energy Release Rate in Anisotropic Materials Based on Virtual Crack Closure-Integral Method Sami Derouiche Civil Engineering department, University of August 20, 1955, Skikda, Algeria. sami.derouiche25@gmail.com, http://orcid.org/0000-0002-3171-1265 Salah Bouziane, Hamoudi Bouzerd Civil Engineering department, University of August 20, 1955, Skikda, Algeria. Laboratory of Civil Engineering and Hydraulics, University of May 8, 1945, Guelma, Algeria. bouziane_21@yahoo.fr, http://orcid.org/0000-0002-0831-1167 bib_ham@yahoo.fr, http://orcid.org/0000-0002-8216-2980 ABSTRACT. 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. KEYWORDS. Fracture Mechanics; Finite Element Method; Anisotropy; Energy Release Rate; Virtual Crack Closure-Integral method. Citation: Derouiche, S., Bouziane, S., Bouzerd, H., Mixed finite element computation of Energy Release Rate in anisotropic materials based on virtual crack closure-integral method, Frattura ed Integrità Strutturale, 57 (2021) 359-372. Received: 12.05.2021 Accepted: 16.06.2021 Published: 01.07.2021 Copyright: © 2021 This is an open access article under the terms of the CC-BY 4.0, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. 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–5]. I https://youtu.be/xhp9MEwTcgQ S. Derouiche et alii, Frattura ed Integrità Strutturale, 57 (2021) 359-372; DOI: 10.3221/IGF-ESIS.57.26 360 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–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 (u1,u2). 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: T S. Derouiche et alii, Frattura ed Integrità Strutturale, 57 (2021) 359-372; DOI: 10.3221/IGF-ESIS.57.26 361     u N q (1) where:    1 1 2 2 3 3 4 4 5 51 2 1 2 1 2 1 2 1 2, , , , , , , , ,tq u u u u u u u u u u (2) {q}t is the vector of nodal displacements and [N] is the matrix of interpolation functions for displacements. Figure 1: Different stages for the mixed finite element formulation. The shape functions are: S. Derouiche et alii, Frattura ed Integrità Strutturale, 57 (2021) 359-372; DOI: 10.3221/IGF-ESIS.57.26 362                                           1 2 3 4 2 5 1 1 1 2 1 1 1 4 1 1 1 4 1 1 1 4 1 1 1 2 N N N N N (3) The stress field in any point is written       M (4) 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:                                 8 11 9 11 10 11 11 11 1 1 2 1 2 4 1 1 2 1 2 4 1 1 2 1 2 4 1 1 2 1 2 4 M M M M (5) The shape functions used to calculate 𝜎 and 𝜎 is given as follows:                                  6 2 7 2 8 2 9 2 1 1 2 1 2 6 1 1 2 1 2 6 1 1 2 1 3 1 1 2 1 4 1, 2 i i i i M M M M i (6) The nodal estimation of the displacement and stress fields is expressed by:                                       0 0 M qB (7) where [B] is the strain-displacement transformation matrix. S. Derouiche et alii, Frattura ed Integrità Strutturale, 57 (2021) 359-372; DOI: 10.3221/IGF-ESIS.57.26 363 The element matrix [Ke] is given by                     0 u e t u K K K K (8) where           e T e A K h M S M dA (9) and         e T e u A K h M B dA (10) [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:     eF K U (11) 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. Figure 2: The two configurations to consider for the crack closure-integral technique: a) the crack opening at the tip for a length of Δl. b) the stress needed to close the crack for the same length. The ERR for this technique is expressed by Irwin [40] for Mode I, as:     00 1 lim , 0 , 2 l I y l G l r v r dr l           (12) T S. Derouiche et alii, Frattura ed Integrità Strutturale, 57 (2021) 359-372; DOI: 10.3221/IGF-ESIS.57.26 364 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:     00 1 lim , 0 , 2 l II y l G l r u r dr l           (13) 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) ( 2 ) 0 1 2 l I yG v dx l       (14) ( 1) ( 2 ) 0 1 2 l II xyG u dx l       (15) 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) as: 43' 4 3' 76' 7 6' 43' 4 3' 76' 7 6' u u u u u u v v v v v v             (16) Figure 3: The two configurations of the crack closure integral technique occurred in RMQ7: a) Forces that act to extend the crack lips. b) Displacements resulted from the crack extension. S. Derouiche et alii, Frattura ed Integrità Strutturale, 57 (2021) 359-372; DOI: 10.3221/IGF-ESIS.57.26 365 Defining the forces and the displacement will allow the evaluation of the ERR for the modes I and II as:     5 43' 6 76' 5 43' 6 76' 1 2 1 2 I y y II x x G F v F v l G F u F u l           (17) The addition gives the total ERR: I IIG G G  (18) Stiffness derivative procedure (SDP) The potential energy π of a cracked body containing a crack is given: 1 2 T Tu Ku u F   (19) 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: 1 2 T KG u u l l        (20) Or as Bouziane et al. [36] states it:                1 1 2 nf T f f f f K G v v a (21) where: a : initial crack length. δa: crack length extension. 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 M0 and the stiffness of the same elements of the configuration M1 (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 M1. 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 F S. Derouiche et alii, Frattura ed Integrità Strutturale, 57 (2021) 359-372; DOI: 10.3221/IGF-ESIS.57.26 366 (22) Figure 4: The two configurations for the energy-release-rate calculation with the stiffness derivative procedure. Figure 5: An arbitrary body with its material properties in the (x,y) coordinate system and an (x',y') coordinate system inclined with φ from the principal axes with a change of the values of the material properties. where:                                  12 1 1 11 12 13 21 12 22 23 2 2 13 32 33 12 1 0 1 0 1 0 0 E ES S S S S S S E E S S S G   (23)       S S. Derouiche et alii, Frattura ed Integrità Strutturale, 57 (2021) 359-372; DOI: 10.3221/IGF-ESIS.57.26 367 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:  { } [ ]{ '}T (24) where the transformation matrix [T] is given:            2 2 cos ² sin ² 2 cos sin [ ] sin ² cos ² 2 cos sin cos sin cos sin cos sin T               (25) It is now possible to have the relation between the two compliance matrices of the two coordinate systems, such as:        1'S T S T (26) Moreover, the stress-strain relation of the same body in a coordinate system oriented at an angle φ is considered as anisotropic and defined as:      ' ' 'S (27) where                                            21 1 121 2 11 12 13 12 2 12 22 23 121 2 13 32 33 1 2 12 12 12 1 1 ' 1 GE E S S S S S S S GE E S S S G G G       (28) 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: 1 1 4 2 2 41 111 12 2 1 cos ( 2 )sin cos sin E E E ES G E             (29) 2 2 4 2 2 42 222 21 1 1 cos ( 2 )sin cos sin E E E ES G E             (30) 12 12 4 4 2 212 1233 12 1 2 1 2 2 sin cos 2( (1 2 ) 1)sin cos G G G GS E E               (31) S. Derouiche et alii, Frattura ed Integrità Strutturale, 57 (2021) 359-372; DOI: 10.3221/IGF-ESIS.57.26 368 and so on for every elastic constant: 2112 22 S S       (32) 12 121 12 11 2 S E S E           (33) 131 33 S S      (34) 232 33 S S      (35) 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 E1 = 144.8 GPa and E2 = 11.7 GPa, Poisson’s ratio of ν12 =0.21and a shear modulus of G12 = 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. Figure 6: An edge central cracked plate submitted to traction. I S. Derouiche et alii, Frattura ed Integrità Strutturale, 57 (2021) 359-372; DOI: 10.3221/IGF-ESIS.57.26 369 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) is used. ϕ E1(GPa) E2(GPa) ν1 G(GPa) η1 η2 0° 144.8 11.7 0.21 9.66 0 0 30° 35.6 14.84 0.1066 10.27 -0.38 -0.33 60° 14.84 35.6 0.0445 10.27 -0.33 -0.37 90° 11.7 144.8 0.016968 9.66 0 0 Table 1: The elastic constants for orientated axes by 30°, 60°, and 90°. Angle FEM [16] ES FEM [16] RMQ7 Analytical solution SDP CCI -90° 4.1599 4.2818 4.3033 4.296 4.3505 -60° 6.3409 6.7759 6.8659 6.966 7.0432 -30° 10.682 11.432 11.5916 11.604 11.827 0° 13.228 13.783 14.7764 14.197 14.251 15° 12.451 13.11 13.5575 13.355 13.568 30° 10.682 11.432 11.4859 11.604 11.827 45° 8.4684 9.1359 9.1907 9.288 9.4331 60° 6.3409 6.7759 6.8454 6.966 7.0432 90° 4.1599 4.2818 4.3033 4.296 4.3505 Table 2: ERR (10-2 J/m²) for the anisotropic plate under tension with different axes rotations. Figure 7: Chart of the ERR for the five finite element methods compared to the analytical ones. S. Derouiche et alii, Frattura ed Integrità Strutturale, 57 (2021) 359-372; DOI: 10.3221/IGF-ESIS.57.26 370 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”. ACKNOWLEDGMENT his work was supported by the Laboratory of Civil Engineering and Hydraulic (LGCH) of the University of May 8, 1945 –Guelma, Algeria. REFERENCES [1] Akinwande, D., Brennan, C.J., Bunch, J.S., Egberts, P., Felts, J.R., Gao, H., Huang, R., Kim, J.S., Li, T., Li, Y., Liechti, K.M., Lu, N., Park, H.S., Reed, E.J., Wang, P., Yakobson, B.I., Zhang, T., Zhang, Y.W., Zhou, Y., Zhu, Y. (2017). A review on mechanics and mechanical properties of 2D materials—Graphene and beyond, Extrem. Mech. Lett., 13, pp. 42–77, DOI: 10.1016/j.eml.2017.01.008. [2] Liu, N., Hong, J., Pidaparti, R., Wang, X. (2016). Fracture patterns and the energy release rate of phosphorene, Nanoscale, 8(10), pp. 5728–36, DOI: 10.1039/c5nr08682e. [3] Larijani, N., Brouzoulis, J., Schilke, M., Ekh, M. (2014). The effect of anisotropy on crack propagation in pearlitic rail steel, Wear, pp. 1–12, DOI: 10.1016/j.wear.2013.11.034. [4] Leitner, T., Hohenwarter, A., Pippan, R. (2019). Anisotropy in fracture and fatigue resistance of pearlitic steels and its effect on the crack path, Int. J. Fatigue, 124, pp. 528–536, DOI: 10.1016/j.ijfatigue.2019.02.048. [5] Toribio, J., González, B., Matos, J.C. (2015). Crack tip fields and mixed mode fracture behaviour of progressively drawn pearlitic steel, Frat. Ed Integrita Strutt., 9(33), pp. 221–228, DOI: 10.3221/IGF-ESIS.33.28. [6] Hattori, G., Alatawi, I.A., Trevelyan, J. (2017). An extended boundary element method formulation for the direct calculation of the stress intensity factors in fully anisotropic materials, Int. J. Numer. Methods Eng., 109(7), pp. 965– 981, DOI: 10.1002/nme.5311. [7] Tafreshi, A. (2017). An analytical expression for the J 2 -integral of an interfacial crack in orthotropic bimaterials, DOI: 10.1111/ffe.12587. [8] Eggleston, J.J., McFadden, G.B., Voorhees, P.W. (2001). A phase-field model for highly anisotropic interfacial energy, Phys. D Nonlinear Phenom., 150(1–2), pp. 91–103, DOI: 10.1016/S0167-2789(00)00222-0. [9] Li, B., Peco, C., Millán, D., Arias, I., Arroyo, M. (2014). Phase-field modeling and simulation of fracture in brittle materials with strongly anisotropic surface energy, DOI: 10.1002/nme. [10] Zhang, S., Jiang, W., Tonks, M.R. (2020). A new phase field fracture model for brittle materials that accounts for elastic anisotropy, Comput. Methods Appl. Mech. Eng., 358, pp. 112643, DOI: 10.1016/j.cma.2019.112643. [11] Nguyen, T.T., Réthoré, J., Baietto, M.C. (2017). Phase field modelling of anisotropic crack propagation, Eur. J. Mech. I T S. Derouiche et alii, Frattura ed Integrità Strutturale, 57 (2021) 359-372; DOI: 10.3221/IGF-ESIS.57.26 371 A/Solids, 65, pp. 279–288, DOI: 10.1016/j.euromechsol.2017.05.002. [12] Shanthraj, P., Svendsen, B., Sharma, L., Roters, F., Raabe, D. (2017). Elasto-viscoplastic phase field modelling of anisotropic cleavage fracture, J. Mech. Phys. Solids, 99, pp. 19–34, DOI: 10.1016/j.jmps.2016.10.012. [13] Hattori, G., Rojas-Díaz, R., Sáez, A., Sukumar, N., García-Sánchez, F. (2012). New anisotropic crack-tip enrichment functions for the extended finite element method, Comput. Mech., 50(5), pp. 591–601, DOI: 10.1007/s00466-012-0691-0. [14] Liu, G.R., Nguyen-Thoi, T., Lam, K.Y. (2009). An edge-based smoothed finite element method (ES-FEM) for static, free and forced vibration analyses of solids, J. Sound Vib., 320(4–5), pp. 1100–1130, DOI: 10.1016/j.jsv.2008.08.027. [15] Chen, L., Liu, G.R., Nourbakhsh-Nia, N., Zeng, K. (2010). A singular edge-based smoothed finite element method (ES- FEM) for bimaterial interface cracks, Comput. Mech., 45(2–3), pp. 109–125, DOI: 10.1007/s00466-009-0422-3. [16] Chen, L., Liu, G.R., Jiang, Y., Zeng, K., Zhang, J. (2011). A singular edge-based smoothed finite element method (ES- FEM) for crack analyses in anisotropic media, Eng. Fract. Mech., 78(1), pp. 85–109, DOI: 10.1016/j.engfracmech.2010.09.018. [17] Peng, F., Huang, W., Zhang, Z.Q., Fu Guo, T., Ma, Y.E. (2020). Phase field simulation for fracture behavior of hyperelastic material at large deformation based on edge-based smoothed finite element method, Eng. Fract. Mech., 238, DOI: 10.1016/j.engfracmech.2020.107233. [18] Surendran, M., Lee, C., Nguyen-Xuan, H., Liu, G.R., Natarajan, S. (2021). Cell-based smoothed finite element method for modelling interfacial cracks with non-matching grids, Eng. Fract. Mech., 242, pp. 107476, DOI: 10.1016/j.engfracmech.2020.107476. [19] Ayatollahi, M.R., Nejati, M., Ghouli, S. (2020). The finite element over-deterministic method to calculate the coefficients of crack tip asymptotic fields in anisotropic planes, Eng. Fract. Mech., 231, pp. 106982, DOI: 10.1016/j.engfracmech.2020.106982. [20] Song, C., Wolf, J.P. (2002). Semi-analytical representation of stress singularities as occurring in cracks in anisotropic multi-materials with the scaled boundary finite-element method, Comput. Struct., 80(2), pp. 183–197, DOI: 10.1016/S0045-7949(01)00167-5. [21] Su, R.K.L., Sun, H.Y. (2003). Numerical solutions of two-dimensional anisotropic crack problems, Int. J. Solids Struct., 40(18), pp. 4615–35, DOI: 10.1016/S0020-7683(03)00310-X. [22] Rajesh, K.N., Rao, B.N. (2010). Two-dimensional analysis of anisotropic crack problems using coupled meshless and fractal finite element method, Int. J. Fract., 164(2), pp. 285–318, DOI: 10.1007/s10704-010-9496-3. [23] Yang, Y., Chu, S.J., Huang, W.S., Chen, H. (2020). Crack growth and energy release rate for an angled crack under mixed mode loading, Appl. Sci., 10(12), DOI: 10.3390/app10124227. [24] Leguillon, D., Murer, S. (2008). A criterion for crack kinking out of an interface, Key Eng. Mater., 385–387, pp. 9–12, DOI: 10.4028/www.scientific.net/kem.385-387.9. [25] Gosz, M., Okyar, A.F., Nair, S. (2003). A crack driving force criterion for the prediction of interface crack kinking in thin-film composites, Interface Sci., 11(3), pp. 329–338, DOI: 10.1023/A:1025100507702. [26] Argatov, I.I., Nazarov, S.A. (2002). Energy release caused by the kinking of a crack in a plane anisotropic solid, J. Appl. Math. Mech., 66(3), pp. 491–503, DOI: 10.1016/S0021-8928(02)00059-X. [27] Cornetti, P., Sapora, A., Carpinteri, A. (2014). T-stress effects on crack kinking in finite fracture mechanics, Eng. Fract. Mech., 132(1), pp. 169–176, DOI: 10.1016/j.engfracmech.2014.10.011. [28] Spagnoli, A., Carpinteri, A., Vantadori, S. (2013). On a kinked crack model to describe the influence of material microstructure on fatigue crack growth, Frat. Ed Integrita Strutt., 7(25), pp. 94–101, DOI: 10.3221/IGF-ESIS.25.14. [29] Krueger, R. (2004). Virtual crack closure technique: History, approach, and applications, Appl. Mech. Rev., 57(2), pp. 109, DOI: 10.1115/1.1595677. [30] Xie, D., Waas, A.M., Shahwan, K.W., Schroeder, J.A., Boeman, R.G. (2004). Computation of energy release rates for kinking cracks based on virtual crack closure technique, C. - Comput. Model. Eng. Sci., 6(6), pp. 515–24. [31] Zeng, W., Liu, G.R., Jiang, C., Dong, X.W., Chen, H.D., Bao, Y., Jiang, Y. (2016). An effective fracture analysis method based on the virtual crack closure-integral technique implemented in CS-FEM, Appl. Math. Model., 40(5–6), pp. 3783– 800, DOI: 10.1016/j.apm.2015.11.001. [32] De Carvalho, N. V., Krueger, R. (2016).Modeling Fatigue Damage Onset and Progression in Composites Using an Element-Based Virtual Crack Closure Technique Combined with the Floating Node Method. American Society for Composites Thirty-First Technical Conference, 1102. [33] Okada, H., Kamibeppu, T. (2005). A virtual crack closure-integral method (VCCM) for three-dimensional crack problems using linear tetrahedral finite elements, C. - Comput. Model. Eng. Sci., 10(3), pp. 229–38, DOI: 10.3970/cmes.2005.010.229. S. Derouiche et alii, Frattura ed Integrità Strutturale, 57 (2021) 359-372; DOI: 10.3221/IGF-ESIS.57.26 372 [34] Bouzerd, H. (1992).Elémént fini mixte pour interface cohérente ou fissurée. Université de Claude Bernard, Lyon I, France. [35] Bouziane, S., Bouzerd, H., Guenfoud, M. (2009). Mixed finite element for modelling interfaces, Eur. J. Comput. Mech., 18(2), pp. 155–75, DOI: 10.3166/ejcm.18.155-175. [36] Salah, B., Hamoudi, B., Noureddine, B., Mohamed, G. (2014). Energy release rate for kinking crack using mixed finite element, Struct. Eng. Mech., 50(5), pp. 665–677, DOI: 10.12989/sem.2014.50.5.665. [37] Mohamed Ben Ali, A., Bouziane, S., Bouzerd, H. (2021). Computation of mode I strain energy release rate of symmetrical and asymmetrical sandwich structures using mixed finite element, Frat. Ed Integrità Strutt. Ed Integrità Strutt., 15(56), pp. 229–239, DOI: 10.3221/igf-esis.56.19. [38] Benmalek, H., Bouziane, S., Bouzerd, H. (2021). Mixed Finite Element for the Analysis of FGM Beams, Int. Rev. Mech. Eng., 15(1), pp. 36, DOI: 10.15866/ireme.v15i1.19680. [39] Rybicki, E.F., Kanninen, M.F. (1977). A finite element calculation of stress intensity factors by a modified crack closure integral, Eng. Fract. Mech., 9(4), pp. 931–938, DOI: 10.1016/0013-7944(77)90013-3. [40] Irwin, G.R. (1957). Analysis of stresses and strains near the end of a crack traversing a plate. J. Applied Mechanics, 24,(1957) 361–364. [41] Voyiadjis, G.Z., Kattan, P.I. (2005). Mechanics of composite materials with MATLAB, 1. [42] Hwu, C. (2010). Anisotropic Elastic Plates, vol. 01, Boston, MA, Springer US. [43] Lekhnitskii, S.G. (1968). Anisotropic plates, Gordon and Breach,. [44] Bowie, O.L., Freese, C.E. (1972). Central crack in plane orthotropic rectangular sheet, Int. J. Fract. Mech., 8(1), pp. 49– 57, DOI: 10.1007/BF00185197.