A computational framework for predicting onset and crack propagation in composite structures via eXtended Finite Element Method XFEM

The eXtended Finite Element Method XFEM has been reliably used for analyzing crack growth in 3D structural elements over last years. In fact, many researchers have worked in this field, but it is scarce to find scientific contributions about 3D XFEM models applied to the failure of non-standard composite parts, such as tapered structures and thick laminated composites. Thus, a new computational framework is developed, which is based on a new enhanced golden section search algorithm and 3D Puck’s action plane principle in order to define the crack initiation direction. This in-formation is integrated into a XFEM and used to enrich elements, which have failed during analysis. Compared to the traditional algorithm, the new methodology has convergence one order higher than the traditional one; and it is 20 times more efficient computationally. Therefore, if more precision is needed, then higher gains are achieved combined to lower computational cost by using the proposed framework. Moreover, thick laminated composites with layers mainly oriented to 90o were simulated under tension and compression via the computational framework, displaying results as reported in the literature. Also, compact tension tests with 0°, 90° and 45° specimens were evaluated, and numerical results were qualitatively coherent with experimental data.


INTRODUCTION
There was undoubtedly a huge advance on the analysis of failure on composite structures in later years, nevertheless the concept of damage tolerance in such structures is still a field wide open for researches and poses many issues already unsolved.Even with modern approaches involving Fracture Mechanics and Representative Volume Elements, the intrinsic structural behavior is too complex and, therefore, micro-scale level approaches bring challenges to be applied widely in an industry environment such as aeronautics damage tolerance MOU et al. 2016 .Thus, one alternative consists on using meso-scale phenomenological damage models based, for example, on the Puck's theory.In fact, this theory is widely used and found in the literature due to easy implementation and good predictions.To corroborate this affirmation, it is enough to see that Puck's theory was one of the best evaluated in the first World Wide Failure Exercise WWFE I and showed good results in WWFE II, as well, "Exhibiting good predictive capability, none or one fundamental weakness and many relatively minor weaknesses" as stated by Hinton, Kaddour andSoden 2002 andKaddour, Hinton andSoden 2004 .WWFE II was the second edition of this benchmark, which tried to evaluate and compare the response of composite mechanics models.
Many researchers have evaluated and proposed modified and enhanced versions of the Puck's model for specific applications as Dong et al. 2016 , who applied a modified version to analyze fatigue under uniaxial loading.Brunbauer, Gaier and Pinter 2015 who, in a similar way, have studied multiaxial loadings.In opposition to World Wide Failure Exercises, which focused on plane stress problems, Puck's theory was used with successful for predicting translaminar failure, but in rare situations.In fact, translaminar failure mechanism involves a 3D stress state and can be found in parts like tapered structures and joints as shown by Figure 1, as well as thick laminated composites.These structures are widely used in rotorcraft industries as primary and vital parts.Therefore, correct and Angelo, Marcus Vinícius a* Ribeiro, Marcelo Leite b Tita, Volnei b efficient reproduction of mechanical behavior of this type of composite structures is relevant not only for academy, but also for industry.Regarding XFEM approach, it has been widely used on numeric simulation, with a huge diversity of enrichment functions, due to its intrinsic capability to model complex geometries independently of the mesh.It gets even more powerful when a previous knowledge of solution is available since it allows for specific enrichment functions to be used FREITAS et al. 2015 .Taking profit from the fact that the approach can map crack evolution on complex structures, many authors have tried to couple XFEM routines to non-or commercial programs.Sukumar et al. 2003 managed to integrate routines to DynELA for dynamic elastic XFEM.Later, following the same research line, Giner et al. 2009 included a model in a User Element subroutine UEL via AbaqusTM solver.Nevertheless, only 2D structures could be evaluated and the algorithm was based on plane stress assumption to be used for isotropic materials.Also, by using AbaqusTM, Shi et al. 2008 implemented a library, easily accessed from the software that allowed user to predict onset and propagation of cracks with-out further intervention from the user, but for only 2D application.From an aeronautical industry point of view, this integration to commercial programs is important, because it usually implies robustness and computational efficiency.Considering only proprietary routines, an extension to 3D space was presented by Shi et al. 2010 , where the authors used level-set and fast-marching methods to track the crack developed in isotropic materials.In addition, there were researchers, who used Puck's criterion along with or as basis for other approaches.Deveci, Aydin and Artem 2016 coupled the theory as constraint for buckling analysis in a hybrid approach.Lee et al. 2015 tried to make a damage-coupled FEM.However, an application of the action plane theory combining XFEM enrichment initiation as well as crack direction calculator has not been found in the literature by the present authors.
Based on the scenario pointed above, the present work aims to propose a new computational framework, combining XFEM and 3D Puck's failure theory for evaluating onset and propagation of failure evolution of cracks on aeronautical tapered structures made of composite materials.How-ever, preliminary studies performed by industry classified results and academy HOCHARD, MIOT and THOLLON, 2014 indicated that this type of structure shows translaminar failure mechanism, which tends to occur in the region close to not exactly on the section transition, on the side of the smaller thickness composite portion, where there are not resin pockets.Thus, the present work aims to predict onset/initialization and propagation of failure in thick laminated composites, which represents the behavior of translaminar failure mechanism in tapered composite structures.In other words, a new computational framework is developed, that is based on an enhanced golden section search algorithm and 3D Puck's action plane principle in order to define the crack initiation direction.This information is integrated into A computational framework for predicting onset and crack propagation in composite structures via eXtended Finite Element Method XFEM Latin American Journal of Solids and Structures, 2018, 15 11 Thematic Section , e70 3/14 XFEM and used to enrich elements that have failed during analysis.The combination of Puck's theory, enhanced golden section search algorithm and XFEM applied to fiber reinforced composite materials is new and makes the numerical solution extremely efficient computational-wise, which is very important to solve academy and industry problems.In terms of the computational aspects, the proposed framework is composed by two main modules: 1 onset/initialization and 2 propagation of failure, which will be discussed in present manuscript.

METHODS
To correctly reproduce the behavior of a composite structure, it is necessary to model with good accuracy polymer matrix behavior as well as out-of-plane stresses.In order to proper model it, the present work invokes Mohr's assumption, where the collapse of the matrix is dependent on the stress state acting on the plane in which failure occurs.On the other hand, since the fibers are the main responsible for providing stiffness in a composite structure, their failure strongly modifies the stress/strain field and quickly reduces its strength, which implies stress concentration and convergence problems on FEA Finite Element Analysis due to the fact that fiber failure mechanisms are assumed to be very brittle and abrupt.
In order to overcome singularities and computational issues that appear near the failure crack , solution via eXtended Finite Element Method XFEM was applied.Matzenmiller's theory, which is based on a Weibull distribution of fiber failures was used to simulate mechanisms in the longitudinal direction aligned to fibers , and Puck's theory was used to simulate failure caused by normal stresses in transverse directions normal to fibers and shear stresses DEUSCHLE 2010 .In addition, a non-local approach was adopted to avoid premature stops in calculus due to convergence problems when fiber fails as proposed by Miot, Hochard, and Lahellec 2010 .The material model was previously implemented as a User Material subroutine UMAT into AbaqusTM, and full details can be seen in Angelo, Galucio, and Tita 2012a and Angelo, Galucio, and Tita 2012b .
Initiation and propagation routines are organized as shown by the analysis flowchart in Figure 2, where different colors represent the sub-modules, which were implemented.First of all, state variables are initialized and load is applied.A linear elastic calculation is performed according to generalized Hooke's law applied to the composite structure.The user subroutine feeds the constitutive matrix to native XFEM solver from Abaqus TM , which performs the solution of the equilibrium equation.The UMAT routine for golden section algorithm is applied to define the action plane AP of failure in order to calculate the failure indexes according to Puck's theory and to update state variables.To enhance convergence performance, these calculations are carried out over space averaged mechanical properties obtained by Angelo, Galucio and Tita 2012a .A new linear elastic analysis is performed with the new set of properties updated properties .If failure indexes exceed failure threshold index, subroutine UDMGINI is called to determine fracture angle inclination and enrichment application.The UDMGINI is a subroutine implemented in FORTRAN, which runs at each end of step during Abaqus TM analysis.It must be coded so that it provides as output a boolean parameter for applying enrichment and the angle according to the element that will be enriched.This procedure is repeated until load is fully applied.The main difficult for implementing Puck's theory is the proper determination of the action plane, the basis of this theory, and it is defined as the plane in which the Puck's indexes are maximum DEUSCHLE 2010 .Works have been done on optimizing this search such as Szlosarek et al. 2016 , but a definitive method is still an open problem.An example of this plane is shown in Figure 3.The action plane inclination in 2D problems can be defined analytically PUCK, KOPP and KNOPS 2002 , but for 3D problems, due to transcendental equations, no analytical solution is possible making it necessary to use numerical methods in order to determine the angle.As previously stated, in the present work, User Defined Damage Initiation UDMGINI subroutine is implemented into AbaqusTM.According to this subroutine, the user must provide the angle and threshold for crack onset DASSAULT SYSTÈMES SIMULIA CORP. 2011 .
The default XFEM from AbaqusTM is able to proper monitor cracks intrinsically, but it cannot be used directly since it supports only Linear Elastic Fracture Mechanics LEFM for isotropic materials.To overcome this limitation, user defined subroutines were employed.For these subroutines, even if the elements are homogeneous, the effect of fibers are taken into account by the material model and, the fracture plane considers the presence of fibers even if they are not explicitly modeled.Considering some regions of the tapered structures or joints, this feature would also allow to model resin pockets and regions without fibers, but an isotropic criterion should be used in these regions 2.1.1Action plane inclination: Puck's approach According to Puck's theory, in order to proper describe the mechanisms of failure related to the polymer matrix inter-fiber failure , a new coordinate system is defined at each finite element mesh integration point regarding the Puck's "inclination" of the "action plane" in which the stresses are acting.Three modes of failures are defined concerning the sign and ratio between components of the stress tensor.Each mode has its own governing law and set of parameters that require tests and assumptions to be determined.
By averaging the stress tensor in the new orientation system, a failure index, which is called "stress exposure", is calculated according to Eq. 1 and used to quantify the amount of damage in the matrix at each integration point.This index will be argument for the determination of the damage variable for inter-fiber mechanism as shown by Eq. 2 .This will be the variable by which elastic moduli will be reduced following Eq. 3 .It is important to notice that it was used the same notation of Puck's works.The action plane angle may be determined through "brute force", i.e. a sweep made through all possible angles and the index of Inter-Fiber Failure IFF is calculated for each one.The highest value is used to determine the critical angle and the inclination angle, considering the respective failure index.This method is extremely computationally inefficient; even though it is assured that the right angle will be reached.Figure 4 illustrates how the failure index behaves when calculated for each inclination.To increase precision, it is strategic to choose smaller steps for sweeping the angles, always taking it implies on higher number of steps and amount of calculations required.Figure 4 displays failure index plotted against failure angle.A sweep from -90° to 90° is performed at each step.As analysis evolves, failure indexes assume higher values.In this example analysis, there is no crack and results are shown up to 90% of complete analysis.
The calculation of action angle is performed through the subroutine UDMGINI from Abaqus.This is used specifically to set a user defined damage initiation criterion for enriched elements.As shown in Figure 2, the routine is called at the end of each increment and during its execution.Thus, the user-defined code will check threshold failure index according to the failure function Equation 1 and predict the crack initiation direction.

Action plane inclination: Proposed approach
As shown by Figure 4, there are multiple points of maxima and minima locally.Therefore, a robust approach was proposed in order to increase efficiency and avoid this problem.
Following the works of Wiegand, Petrinic and Elliott 2008 , it was implemented the approach based on the so called "golden section algorithm".Thus, the interval of interest, in this case -90°, 90° , was subdivided as per the golden ratio shown in Eq. 4 : The subdivision follows the flowchart shown by Figure 5.At each iteration, the subinterval is reduced up to a desired minimum size that should contain the sought angle.In this example, the stop criterion å was taken to be equal 2 10  .This value was taken since after 1/100 of degree, the IFF parameter is not affected.The Golden Section method is efficient, but it shows some limitation.According to the distribution of maxima and minima, the method may be unable to find the real global maximum, returning a local extreme value.
In order to improve efficiency, in the present work, the Golden Section method was expanded taking into account the hypothesis from Schirmaier et al. 2014 , which states that no extreme value appears closer than 25° from each other.Under this hypothesis, a new approach is proposed that consists in a portion-wise application of Golden Section method.According to this new methodology, the whole interest interval is divided in 9 x 20° and, for each sub-interval, the Golden Section method is applied, and this new approach is defined as enhanced Golden Section method.Nine potential maximum values are obtained and compared.Hence, the maximum among them is considered the global maximum as shown by Figure 6.In order to exemplify the performance enhancement by using the proposed approach, Table 1 shows results for each level of precision.From the Table 1, one can see clearly the performance improvement in terms of precision level for the determination of action plane angle inclination due to the usage of enhanced golden section method, when in comparison with traditional methodology proposed in Puck's original work.The proposed approach is, for instance, more than 20 times more efficient computationally for a 0.1° precision than Puck's original method.The gains of applying the enhanced golden section method are not as high when in comparison to original golden section algorithm in terms of computational cost.Nevertheless, the balance between robustness and performance makes of it the best choice that could be further extended to other works that made use of this algorithm such as meshless colocation TSAI, KOLIBAL, and LI 2010 , medical appliance LU and HARGREAVES 2008 and even aircraft control PATRON, BOTEZ, and LABOUR 2012 .

RESULTS AND DISCUSSION
As commented previously, there were preliminary studies performed by industry classified results and academy HOCHARD, MIOT and THOLLON, 2014 , which indicated that tapered composite structures show translaminar failure mechanism.Furthermore, onset/initialization and propagation of failure in thick laminated composites represents the same behavior of translaminar failure mechanism.Based on these aspects, in the present work, it is proposed to investigate two different case studies: 1 thick laminates with layers mainly oriented to 90o under tension and compression and 2 0°, 90° and 45° specimens of compact tension test.For both case studies, material properties can be found in Table 2, where the composite material T700/M21 consists on epoxy res-in reinforced by carbon fibers, and sub-index 1, 2 and 3 corresponds to the orthotropic directions.

Case Study 1: 90° layers under tension and compression
When transversal compression load is applied on unidirectional composite, considering experimental results presented by Puck and Schürmann 1998 , it is expected that the angle for failure occur around 54° DEUSCHLE 2010 .Experimental data was obtained from Puck's original work where a fiber reinforced torsional spring was employed to isolate compressive stresses PUCK and SCHÜRMANN 2002 .On the other hand, when the same structure is loaded under tension transversal loading, it is expected that the angle for failure occurs at 0°, i.e. parallel to the fibers aligned to 90o according to laminate coordinate system.
To evaluate the performance of the proposed framework, a computational model was developed in order to simulate those tests.Figure 7 shows the model for tension loading, which was simulated by using a prescribed displacement and hexahedral elements C3D8 -AbaqusTM organized on one element per layer.These solid elements with 8 nodes were chosen since this type of element is able to capture intrinsically out-of-plane phenomena by using XFEM approach via AbaqusTM.In addition, hourglass modes are controlled by default AbaqusTM algorithm.A symmetry boundary condition is used, and results are verified in the control element highlighted by the Figure 7.A mesh dependency study was carried out and, after verifying the convergence of the results, the final model, considered mesh independent, is presented in Figure 7. Final mesh has 24128 nodes and 18842 C3D8 elements.In order to optimize resources, refinement is applied near the edges with a coarser mesh in other regions.
A computational framework for predicting onset and crack propagation in composite structures via eXtended Finite Element Method XFEM Latin American Journal of Solids and Structures, 2018, 15 11 Thematic Section , e70 8/14 The control element works as a virtual strain-gage and its numerical results are acquired for comparison with experiments.For compression simulations, the applied loading was kept below buck-ling critical load, then a linear analysis was performed.It is important to remark that the model aimed to prove the effectivity of the proposed approach Results obtained at the control element for the compression analyses are shown in Figure 8.The sequence of curves demonstrates the application of traditional Puck's Method work through a sweeping of every possible angle of failure during each step of load application.The maximum value for each curve represents the sought value for the angle of failure.The red dots correspond to the angle obtained using the newly developed approach.It can be noticed from Puck's work that there are two global maxima values, which were found to be around 53.4° or -53.4°.The phenomenon may be regarded like instability in buckling, i.e. the answer is independent of the choice of the A computational framework for predicting onset and crack propagation in composite structures via eXtended Finite Element Method XFEM Latin American Journal of Solids and Structures, 2018, 15 11 Thematic Section , e70 9/14 path symmetrical behaviors .It is important to highlight that loading was simulated by using an imposed displacement up to the first failure i.e.IFF 1 , which happened at 3.7mm for this case.
In a similar manner, Figure 9 shows the results for the tension analyses.Again, the results are exactly as expected.Failure occurs at 0° i.e. in a plane parallel to the fibers and perpendicular to the loading.Filled dots show the results obtained by the proposed approach in opposition to the curve, which is again obtained through Puck's approach.Again, the imposed displacement was adopted to be 3.7mm, but as it can be noticed, for tension analyses, the total failure index IFF reached much higher values.Table 3 shows comparison the current method to Puck's method for the Case Study 1.

Case Study 2: Compact tension specimen
In order to verify the implementation of the initiation routine in the XFEM via material model, com-pact tension tests under imposed displacement were simulated as shown by the Figure 10.Dis-placement was increased up to final failure of the structure.Results are reported at the last step before complete failure, since there is not solution for the model when the structure is completely separated.
The model was developed using 8805 C3D8 elements to simulate 16 composite layers.Con-strained and imposed displacements are applied in reference points RPs .One of the RPs is con-strained while the other one has displacement steps, which are gradually applied as shown by the Figure 10 A .
Boundary conditions are transferred to the inner portion of each hole by rigid links.It is important to highlight that the compact tension test is a standardized procedure AMERICAN SOCIETY FOR TESTING AND MATERIALS 2015 and it has been commonly used for investigating translaminar failure and crack propagation in laminated composites LAFFAN et al. 2010, LAFFAN et al. 2012, TEIXEIRA, PINHO, and ROBINSON 2011, FERNANDES 2014 .A computational framework for predicting onset and crack propagation in composite structures via eXtended Finite Element Method XFEM Latin American Journal of Solids and Structures, 2018, 15 11 Thematic Section , e70 10/14 Both phenomena, translaminar failure and crack propagation are expected to be found in tapered composite structures, mainly according to the stacking sequence of the laminate.Thus, three different layups were evaluated as shown in Table 4, and respective results for experiments shown in Figure 11 are displayed in Figure 12.Final displacement for each analysis is also presented in Table 4. Since the strength in fiber direction is much higher than in transversal, it is expected the failure represented by trespassing cracks in the structure following the alignment of fibers.
The consistency of the numerical results can be evaluated qualitatively by using experimental data found in the literature as shown by Figure 11, which were obtained from Jamali 2014 andSlepetz andCarlson 1976 .In these experiments compact tension specimens, made of unidirectional composite materials were stressed monotonically up to failure.
Each experiment was performed for unidirectional composites to check how crack growths.Material properties were the same as Case Study 1, and augmented Lagrange method was used as hard contact to model traction separation.Failure progression was modeled through XFEM in which enriched golden section and Puck criterion were used to determine whether crack goes through the element or not.In fact, the shape functions of finite elements are enriched through Heaviside approach, and a discontinuity of displacement field is achieved.After the enrichment, the displacement field is discontinuous in the enriched direction provided by enhanced golden section combined to Puck's theory , but the elements are still present.For instance, all properties not affected by discontinuity are kept.This is particularly important for analyzing fiber reinforced composite due to its orthotropic characteristic.The behavior shown by experiments was just as expected with the crack following the fiber direction and propagating parallel to it Figure 11 .It was exactly the behavior obtained from the computational models as it can be seen in Figure 12 A , B and C .It is important to highlight that the results for the boolean parameter related to the enrichment output of UDMGINI are presented, where "1"for enriched elements is red, and "0" for not enriched ones is blue

Figure 1 :
Figure 1: A Tapered structure and B typical double lap joint evident similarity

Figure 2 :
Figure 2: Calculations for the onset of failure -in red: UMAT, in blue: UDMGINI, in green: Abaqus Default

Figure 3 :
Figure 3: Puck's theory: Action plane WIEGAND, PETRINIC, and ELLIOTT 2008 for predicting onset and crack propagation in composite structures via eXtended Finite Element Method XFEM Latin American Journal of Solids and Structures, 2018, 15 11 Thematic Section ,

Figure 4 :
Figure 4: Example of calculation of action plane inclination as proposed by Puck's algorithm.

Figure 6 :
Figure 6: Schematic representation of enhanced Golden Section method, highlighted one subinterval IFF vs Action Plane Inclination

Figure 7 :
Figure 7: Model under tension transversal loading: Boundary conditions and control element.

Figure 8 :
Figure 8: Failure index for each angle: Compression analyses.

Figure 9 :
Figure 9: Failure index for each angle: Tension analyses.

Figure 10 :
Figure 10: Compact tension model: A 3D iso-view B Detailed sketch with dimensions units: mm .

Table 2 :
Material properties for case studies.

Table 3 :
Case Study 1 -Comparison the current method to traditional Puck's method.

Table 4 :
Details of computational analyses for Case Study 2