Fracture and mode mixity analysis of shear deformable composite beams

To analyse delaminated composite beams with high accuracy under mixed-mode I/II fracture conditions first-, second-, third- and Reddy’s third-order shear deformable theories are discussed in this paper. The developed models are based on the concept of two equivalent single layers and the system of exact kinematic conditions. To deduce the equilibrium equations of the linearly elastic system, the principle of virtual work is utilised. As an example, a built-in configuration with different delamination position and external loads are investigated. The mechanical fields at the delamination tip are provided and compared to finite element results. To carry out the fracture mechanical investigation, the J-integral with zero-area path is introduced. Moreover, by taking the advantage of the J-integral, a partitioning method is proposed to determine the ratio of mode-I and mode-II in-plane fracture modes. Finally, in terms of the mode mixity, the results of the presented evaluation techniques are compared to numerical solutions and previously published models in the literature.

Beam elements with different delamination position over the beam thickness carried out on different type of pre-cracked specimens including mode-I [23], mode-II [8], mixed-mode I/II [58], mode-III [14,25], mixed-mode I/III [1], mixed-mode II/III [43,50] and mixed-mode I/II/III [57] fracture conditions. Finally, these experiments must be evaluated based on very simple mechanical structures where the delamination zone is known a priori, such as delaminated beam and plate structures.
To increase the accuracy of the evaluation techniques and to capture more precisely the complex fracture mechanical behaviour of these mechanical systems, the researchers have been inspired to introduce novel, higher-order beam and plate theories [33,53,54]. The first-order shear deformable beam theory (FSDT) assumes independent rotation around the bending axes [46]. As a next step, the second-order shear deformable theory (SSDT) describes the in-plane displacement as a quadratic function in terms of the through-thickness coordinate. The third-order shear deformable theory (TSDT) captures the longitudinal displacement in the form of cubic function [12]. And finally but not least, in order to satisfy dynamic boundary conditions at the top and bottom surfaces of the beam, Reddy's third-order shear deformable theory is also developed [41,47].
In the following, to analyse a delaminated composite beam under mixed-mode I/II fracture conditions, the basic idea of higher-order beam theories is utilised with the semi-layerwise approach [51]. The concept of the proposed description can be seen in Fig. 1. where beam elements with an artificial delamination are depicted. To perform fracture mechanical investigation, the J -integral with zero-area path is applied [42,50]. Moreover, a partitioning method is proposed to determine the mode-I and mode-II in-plane fracture modes without any semi-analytical considerations [29]. Finally, the obtained results of the mode partitioning are compared to previously published evaluation methods in the literature [4,21,34,56].
Applications already exist in the literature which uses the J -integral to evaluate delaminated composite beams [3,30,44] under pure mode-I or pure mode-II fracture conditions. The effect of reinforcement direction and delamination length on the energy release rate have already been investigated in numerous studies, but J -integral-based mode mixity evaluation, which undeniably requires complicated mechanical models, has not been extensively studied.
It is important to note that the methods described herein are only suitable for describing brittle composite materials. If the brittle carbon or glass fibres are protected by more damage tolerant resin system, the application of the linear elastic fracture mechanics is no longer appropriate. The nonlinear fracture process zone and the growth of multiple cracks at small scale should also be considered [55].

Mechanical model: the method of two equivalent single layers (ESLs)
A delaminated beam with an arbitrary stacking sequence including orthotropic or transversely isotropic plies is depicted in Fig. 2. As it is illustrated by this figure, different mathematical models must be applied to capture the displacement field of the cracked and uncracked portion separately [50]. Firstly, to describe the undelaminated portion, the membrane displacement field must remain continuous along the thickness coordinate. While in the case of the delaminated portion, the presence of the delamination splits the functions into separated displacement fields resulting in two individual sub-laminates [48]. To handle this phenomenon properly, each and every sub-laminate is modelled as an equivalent single layer (ESL). This method is the so-called semilayerwise technique, which has already been introduced and used in numerous studies to describe delaminated composite plates under mixed-mode II/III fracture conditions [45]. Furthermore, thedisplacement field can be a b Fig. 2 Cross section and the assumed deformation of the undelaminated portion in the X −Z plane (a) and distributions of the transverse shear strain by different theories (b) using 2ESLs expressed by using the above-mentioned higher-order polynomials to increase the accuracy in terms of shear stress distribution: where i denotes the index of actual ESL, z (i) is the local thickness coordinate of the ith ESL and always coincides with the local midplane, u 0 is the global, u 0i is the local membrane displacement, θ denotes the rotation of the cross section about the Y axes, φ means the second-order, and λ expresses the third-order polynomial terms. Moreover, w (i) represents the separate transverse deflection functions. The displacement field of FSDT and SSDT can be obtained by substituting φ = 0 and λ = 0 into the Eq. (1), respectively. In the case of the undelaminated part, the kinematic continuity between ESLs can be imposed by the system of exact kinematic conditions [50]. Based on the above assumed displacement fields, the nonzero terms of the strain field become where ε x(i) is the normal strain of the ith ESL in the X -direction, and γ xz(i) is the shear strain of the ith layer in the X −Z plane. The in-plane strains can be decomposed into constant, linear and higher-order terms with respect to the through-thickness coordinate: and: By applying the constitutive equation of orthotropic materials the stress resultants can be introduced for each and: where N x(i) defines the normal force, M x(i) denotes the bending moment, Q x(i) is the transverse shear force, and L x(i) , P x(i) , R x(i) , S x(i) are the higher-order stress resultants [41,50]. By using the through-thickness coordinate dependence and the presented constitutive law, the relationship between the strain field and the stress resultants can be formulated as: and ⎛ where A pp(i) denotes the extensional, B pp(i) represents the coupling, D pp(i) is the bending, E pp(i) , F pp(i) , G pp(i) and H pp(i) are the higher-order stiffnesses. These terms can be calculated by [41,49]: where N l(i) is the number of orthotropic layers within the actual ith ESL, b defines the width of the beam, z i m and z i m+1 are the local bottom and top coordinates of the mth layer in the ith ESL [49]. Furthermore, C (i) is the stiffness matrix of the mth layer within the ith ESL expressed as:

Description of the undelaminated portion
In order to impose the continuity between two ESLs, the following condition can be formulated [50]: (u (1) , w (1) ) (2) , w (2) ) which imposes the displacement continuity and the identity of the transverse deflections. The second condition defines the position of the global reference plane: with respect to ESL1. As a first step, to develop semi-layerwise model only for FSDT solution, these equations are already enough. However, it is reasonable to increase the accuracy of the model further by using SSDT and TSDT. Thus, the shear strain continuity at the interface plane, referring to Fig. 2, can be also imposed: Furthermore, in accordance with the basic concept of Reddy-TSDT theory, even the traction-free boundary condition can be imposed on the bottom and top surfaces of the beam: Finally, by taking into account the previously discussed conditions, in each and every case the in-pane displacement functions can be expressed in the following invariant form: where the K i j matrices depend on the actual ESL thicknesses and the order of the applied theory, i refers to the ESL number, the summation index j defines the component in ψ j primary parameter vector, and w(x) expresses the transverse deflection of the beam [49].

Reddy's third-order beam theory
Using the above-mentioned exact kinematic conditions by Eqs. (11)-(14), 5 unknown displacement parameters can be eliminated from the original 10 parameters. The secondary parameters are: u 0i , φ i for i = 1..2 and θ 2 . The unknown terms and the vector of primary parameters become: The nonzero elements of the K i j and K (3) i j matrices in Eq. (15) can be found in Appendix A.1.1. It is important to highlight, by the application of the Reddy's third-order beam theory, the first derivate of the deflection becomes also a primary parameter. Later on we will see that it results much more complicated differential equation with respect to the transverse deflection than the other higher-order theories.

Third-order beam theory
Using again the above discussed exact kinematic conditions [(referring to Eqs. (11)-(13)], apart from the traction-free condition, 3 unknown displacement parameters can be eliminated from the original 10 parameters. The secondary parameters are: u 0i for i = 1..2 and λ 1 . The unknown terms and the vector of primary parameters are: The nonzero elements of the K i j and K (3) i j matrices in Eq. (15) are defined in Appendix A.2.1.

Second-order beam theory
In the case of SSDT, 3 unknown displacement parameters can be eliminated from the original 8 parameters by using Eqs. (11)- (13). The secondary parameters are: u 0i for i = 1..2 and φ 1 . The unknown terms and the vector of primary parameters become Obviously, K i j = 0 in this case. The elements of the K i j and K (2) i j matrices in Eq. (15) are specified by Appendix A.3.1.

First-order beam theory
As it was previously discussed, in the case of FSDT theory, only the continuity conditions by Eqs. (11)-(12) continuity can be imposed, meaning that 2 unknown displacement parameters are eliminated from the original 6 parameters. These secondary parameters are: u 0(i) for i = 1, 2. The vector of primary parameters becomes Using FSDT theory K i j = 0, the elements of the K i j and K (1) i j matrices in Eq. (15) can be found in Appendix A.4.1.

Equilibrium equations
The equilibrium equations of the undelaminated beam can be deduced from the virtual work principle [41], by setting the sum of coefficients for the virtual membrane displacement δu 0 : the primary parameters δψ j : and the deflection δw: δw : to zero. It is worth giving attention to Eq. (22) where the application of Reddy's third-order theory results a more complex equilibrium equation [41] compared to the other theories.

Description of the delaminated portion
As it was previously discussed, the presence of a delamination divides the beam into a top and a bottom sub-laminates. Thus, diverse mathematical functions are applied to the top and bottom parts: where notations are exactly the same as they were in Eq. (15). Although in this context, u 0b and u 0t express the global membrane displacement of the bottom and top sub-laminates, local membrane displacement parameters are not defined, and finally but not least w b describes the transverse deflection of the bottom beam and w t represents the deflection of the top beam [49].

Reddy's third-order beam theory
According to Reddy's third-order beam theory, the top and bottom beams are traction-free at their top and bottom boundaries, as well. But based on our computational experience, if these conditions are strictly imposed, we always get inconsistent and over-constrained shear strain distribution around the delamination tip. In order to avoid this phenomenon in this paper, the traction-free conditions are only imposed at the top boundary of the top sub-laminate and at the bottom boundary of the bottom sub-laminate. Referring to Fig. 2 and Eq. (13), 2 unknown displacement parameters can be eliminated from the original 10 parameters. The secondary parameters are: φ i for i = 1..2. The unknown terms and the vector of primary parameters become: where circles denote the auto-continuity parameters. The role of these parameters will be important when the displacement continuity between the delaminated and undelaminated portions is discussed. The nonzero elements of the matrices K i j and K (3) i j in Eq. (24) can be found in Appendix A.1.2. Similarly to the case of undelaminated portion, the derivative of the deflections become also primary parameters which results more complex differential equations in terms of the transverse deflections.

Third-order beam theory
Since no other condition can be imposed, the primary parameters become: where circle denotes the above-mentioned auto-continuity parameter [50]. The corresponding K i j matrices in Eq. (24) can be found in Appendix A.2.2.

Second-order beam theory
No other condition can be imposed against the displacement field in the case of SSDT, as well. We only have primary parameters: and K i j = 0. The nonzero elements of the K i j and K (2) i j matrices in Eq. (24) are given in Appendix A.3.2.

First-order beam theory
In the case of FSDT, the primary parameters are: The elements of the K i j and K (1) i j matrices in Eq. (24) can be found in Appendix A.4.2. By using this theory, there is no need to define auto-continuity parameter and K (3) i j = K (2) i j = 0.

Equilibrium equations
Using the principle of virtual work, it is possible to determine the equilibrium equations of the delaminated part also, by setting the sum of coefficients for the virtual membrane displacements δu 0b and δu 0t : the primary parameters δψ j : and the transverse deflection terms δw b and δw t : to zero. As we can see, the application of Reddy's third-order theory results much more complex equilibrium equations in terms of the transverse deflections than the other higher-order theories [41].

Example: built-in configuration of a delaminated beam
As an example a built-in beam with asymmetric delamination is considered. The geometry of the structure is depicted by Fig. 3, where l denotes the total length, a represents the length of the delamination, b is the beam width and 2h is the total thickness of the beam with [±30/0 2 / ± 30/0] S lay-up. The corresponding thicknesses of the sub-laminates are denoted by t 1 and t 2 . The material properties of the transversely isotropic and cross-ply layers are given by Table 1 [49,52]. The structure of the ELSs in each and every case is determined according to Fig. 1.

Continuity conditions of the displacement field
The continuity conditions between the regions can be imposed by using the discussed primary parameters. The conditions for the components of the membrane displacement and the transverse deflection can be ensured by: where q undel is the number of the primary parameters in ψ p vector. Independently of the order of the applied theory, these equations always represent four conditions. The continuity of the first-, second-, and third-order terms in the displacement functions can be specified through the element of the primary parameter vectors by connecting the first-order terms with first-order, the second-order terms with second-order, the third-order with the third-order terms and the deflection derivatives with the corresponding deflection derivatives. These represent q undel number of conditions between the undelaminated and delaminated portions. For the Reddy's TSDT theory, by referring to Eqs. (16) and (25), these conditions can be formulated as: and As the number of parameters in the ψ p vector for the undelaminated and delaminated parts is generally not equal to each other, i.e.: q del = q undel , the continuity of the remaining terms cannot be expressed directly. The remaining θ 2 and ∂w t ∂ x terms in the primary parameter vector of the delaminated part, which are indicated by circles in Eq. (25), can be defined as auto-continuity parameters [50]. According to the theorem of auto-continuity, the continuity of these terms can be ensured by using only the primary parameters of the undelaminated part. The equations can be written as: and representing q undel + 2 = 4 + 2 = 6 for Reddy's TSDT. Regarding TSDT, the situation is exactly the same. The inequality, by referring to Eqs. (17) and (26), can be handled in a similar way: and Fig. 4 Continuity of the stress resultants at the delamination tip meaning q undel + 1 = 5 + 1 = 6 number of conditions. In connection with SSDT, referring to Eqs. (18) and (27), the continuity can be written as: and representing q undel +1 = 3+1 = 4 conditions. Fortunately, in connection with FSDT, by referring to Eqs. (19) and (28), the number of terms in the primary parameter vector for undelaminated and delaminated portions are equal to each other q del = q undel . No further auto-continuity condition is required to be imposed. It can be expressed easily by: meaning q del = q undel = 2 conditions [50].

Continuity conditions of the stress resultants
As a result of arbitrary external forces, and P x(i) unknown stress resultants take place at the left end of the delaminated portion in each ESL. These stress resultants and their derivatives are transferred between the portions resulting equivalent normal force, equivalent shear force and equivalent bending moments on the right hand side of the undelaminated portion. The basic concept of the idea is depicted by Fig. 4, where the same arrows are used to sign R x(i) , S x(i) , and L x(i) , P x(i) . The form of the equations can be read out from the presented equilibrium equations [49]. Thus, the continuity of the normal forces, which means one condition for each theory, can be formulated as: whereN (undel) x represents the equivalent normal force. The next conditions between the portions are the continuity of the equivalent bending moments. These continuity conditions can be formulated as: where K (undel) i j collects the K i j matrix elements of the undelaminated part: denotes the equivalent bending moments imposing four conditions for the Reddy's TSDT, five for the TSDT, three for the SSDT and two for the FSDT. It is worth giving attention to Reddy's TSDT where we obtained an additional bending continuity condition in terms of the deflection derivative (Reddy: j = 4), referring to Eqs. (16) and (22). Finally, the shear force continuity of the problem can be imposed for Reddy's TSDT as:Q and for the other theories: resulting one more condition for each case.

Boundary conditions
In order to impose built-in boundary condition at the left end of the undelaminated part the displacement field must be zero. This condition can be easily prescribed by using the primary parameters: The external forces acting on the right hand side of the sub-laminates can be considered as equivalent stress resultants in accordance with Eqs. (29)-(30) and (31)- (33). In the lack of external forces in the axial direction, two boundary conditions can be imposed for each theory: In the lack of external bending moments, the caused equivalent bending moments can be formulated as: where K (del) i j arranges the K i j matrix elements of the delaminated part: And finally, according to Fig. 3, the top and bottom parts of the delaminated regions are subjected to F t and F b forces, causing the following two boundary conditions for the Reddy's TSDT: and causing two boundary conditions for the other theories:

Displacements and stress distributions
Two different delamination scenarios are considered in this section in order to investigate the performance of the higher-order shear deformable theories. In each end every case, according to the presented Fig. 1, the analysis is carried out based on the novel semi-layerwise beam models. Moreover, to verify the obtained analytical results finite element (FE) analyses with the so-called virtual crack closure technique (FE-VCCT) are carried out [31], as well.
Regarding "Case I" delamination scenario, when the top and bottom sub-laminates are subjected to F t = 10 N and F b = −10 N external forces, the u in-plane displacement fields and the w deflections are depicted by Fig. 5. As we can see, by using any kind of higher-order theories, the results show good agreement with the finite element solution in connection with u in-plane displacement field. In terms of the w deflections the agreement between the numerical and the analytical results is also quite good. The application of the higherorder beam theories results only a bit smaller function values than the FEA solution. As it is also illustrated by the highlighted part of the figure, the contribution of the higher-order theories to the deflection becomes smaller and smaller by increasing the order of theory. Thus, the deflection improvement of the SSDT, TSDT and Reddy's TSDT compared with the FSDT is negligible. From this point of view, application of the higherorder theories is not necessary. In order to present the essence of the higher-order beam theories, Fig. 6 depicts further results regarding the σ x normal and τ xz shear stresses. The stresses are provided at the delamination tip, separately at the undelaminated X = +0 [mm] and the delaminated X = −0 [mm] portions. As we can see, certain stress discontinuities take place at the crack tip which result different σ x normal and τ xz shear stress distributions on the sides. Based on Fig. 6a, apart from the singular nature of the numerical solution, σ x normal stresses are in good agreement with each other. Although, if we take a look at the delaminated part, the Reddy's TSDT solution shows a little perturbation. It can only be explained by the rigorously imposed traction-free conditions. The application of the higher-order beam theories becomes quite important only if the τ xz shear stresses are investigated. Based on Fig. 6b, by increasing the order of the solutions, the τ xz shear stress can be described with a more and more accurate way. If the exact description of the τ xz shear stress at the crack tip is important, the application of higher-order theories becomes inevitable. It is worth giving attention to τ xz shear stress distribution which is obtained by the Reddy's TSDT. As it was discussed previously, it is already able to satisfy the traction-free boundary conditions all along the top and bottom surfaces of the beam. shear stress distribution represents key role of the fracture mechanical investigation, the application of the higher-order theories is unavoidable. Although, if we investigate the delaminated portion of the structure, the perturbation caused by the Reddy's TSDT solution is no longer negligible. Unfortunately, the traction-free condition can significantly disturb the σ x normal stress distribution.

J-integral
In order to perform in-plane fracture mechanical investigation the so-called J -integral can be applied [42]. Based on the basic definition, considering any arbitrary counterclockwise C contour around the crack tip, the J -integral can be formulated as: where U is the strain energy density, σ i j are the components of the stress tensors, u i means the displacement vector components, ds is the length increment along the contour and n i represents the components of the outward unit vector in the given (X 1 − X 2 ) Cartesian coordinate system. The basic concept is depicted by Fig. 9a. Moreover, considering only linear elastic fracture mechanics, the fundamental property of this integral is the following [2]: where G T denotes the total energy release rate and J represents the value of the contour integral. Regarding delaminated beams, and based on the discussed semi-layerwise approach, the J -integral can be calculated as a zero-area path integral. The idea is depicted by Fig. 9b, where the n unit vector always remains parallel to the X -axis. Finally, taking into consideration the actual coordinate system (X 1 = −X and X 2 = z (i) ) and the actual number of ESLs, the total value of the energy release rate becomes Many models were developed in the literature to perform mode partitioning in beam-type fracture specimens, but generally these evaluation strategies, apart from the Suo-Hutchinson method [21], apply semi-analytical considerations [34] or arbitrary assumptions to make the separation feasible [4,56]. Fortunately, the application of the J -integral makes these type of considerations absolutely avoidable and provides a methodical separation technique. For this purpose, we only have to separate the strain and stress fields into symmetrical and asymmetrical function components with respect to the delamination plane: where the (sym) and (ant) subscripts indicate the symmetrical and asymmetrical function components, respectively [29]. To make it understandable, Fig. 10 represents an example based on a previously calculated stress distribution function (TSDT solution of "Case I"). Finally, knowing the symmetrical and asymmetrical function components of the field quantities [29], the mode-I fracture mode becomes and the mode-II fracture mode can be calculated as: First of all, to verify the proposed evaluation technique, the mode mixity of a symmetrically delaminated beam with unidirectional composite layers has to be investigated. The geometry of the beam is exactly the same as it is depicted by Fig. 3, but in this test problem the ESLs consist of transversely isotropic layers. The  [4,21,34,56], the developed higher-order evaluation techniques predict the mode mixity quite well . For "Case I" delamination scenario, the obtained results are depicted by Fig. 12. As can be seen, under different external loadings the FSDT, SSDT and TSDT solutions give results quite close to each other and the Euler-Bernoulli-based evaluation techniques. Furthermore, the agreement with the numerical FEA-VCCT solution is also quite good. Although, the application of Reddy's TSDT significantly decreases the ratio of the mode-I energy release rate. The difference can only be explained by the rigorously imposed traction-free boundary condition. The caused σ x perturbation, referring to Fig. 6, is no longer insignificant and it certainly disturbs the ratio of the symmetrical and asymmetrical function components, as well. Figure 13 gives result of the mode mixity evaluations for "Case II". In general, we can state that the FSDT, SSDT and TSDT are located between the Williams' curvature based and the Reddy's TSDT solutions. The FSDT and the Luo-Tong theory predict almost the same mode mixity values as the FEM-VCCT. Unfortunately, as in the delamination scenario discussed above, the application of Reddy's TSDT significantly decreases the ratio of the mode-I energy release rate, which can only be attributable to the σ x perturbation in Fig. 8.

Conclusion
To analyse brittle interlaminar fracture in composite structures under mixed-mode I/II fracture conditions first-, second-, third-and Reddy's third-order shear deformable theories were discussed in this paper with two ESLs. The displacement continuity between the ESLs was described and imposed by using the system of exact kinematic conditions. Using Reddy's TSDT theory, the traction-free boundary condition was also imposed along the top and bottom surfaces of the beam. As an example, a built-in configuration with different delamination scenario subjected to different external loads was presented. To solve the higher-order beam problems, the corresponding continuity and boundary conditions were discussed. Based on the obtained results, the mechanical fields at the delamination tip were provided and compared to each other. As a next step, to calculate the G T total energy release rate and to perform in-plane mode mixity analysis, the J -integral with zero-area path was introduced. By using symmetric and asymmetric decomposition of the mechanical fields, a well-ordered evaluation technique was proposed.
The results of different higher-order theories, apart form the Reddy's TSDT, were close to each other and the Luo-Tong solution. Furthermore, these evaluation methods gave results between the Williams' curvature based and the Bruno-Greco solutions, which defined the most extreme mode ratios. Unfortunately, regarding Reddy's TSDT it was not true. The rigorously imposed traction-free condition significantly decreased the ratio of mode-I energy release rate and it significantly disturbed the σ x normal stress distribution in the thinner sub-laminate.
Nevertheless, this paper is the first attempt to perform mode separation by combining the higher-order theories with the J -integral. The next step could be the semi-layerwise approach-based finite element development and the consideration of the transverse stretch modelling [52].  In this Appendix the nonzero K i j matrix elements for the second-order beam theory are presented.