Acessibilidade / Reportar erro

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

Abstract

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.

Keywords
3D XFEM; Crack propagation; Puck’s theory; composite joints; tapered structures

1 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. 2016MOU, H. L., T. C. ZOU, Z. Y. FENG, and J. XIE. 2016. “Crashworthiness Analysis and Evaluation of Fuselage Section with Sub-floor Composite Sinusoidal Specimens.” Latin American Journal of Solids and Structures 13 (6):1186-1202.). 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 and Soden (2002)HINTON, M. J., A. S. KADDOUR, and P. D. SODEN. 2002. “A Comparison of the Predictive Capabilities of Current Failure Theories for Composite Laminates, Judged Against Experimental Evidence.” Composites Science and Technology 62 (12):1725-1797. doi: 10.1016/S0266-3538(02)00125-2.
https://doi.org/10.1016/S0266-3538(02)00...
and Kaddour, Hinton and Soden (2004)KADDOUR, A. S., M. J. HINTON, and P. D. SODEN. 2004. “A Comparison of the Predictive Capabilities of Current Failure Theories for Composite Laminates: Additional Contributions.” Composites Science and Technology 64 (3):449-476. doi: 10.1016/S0266-3538(03)00226-4.
https://doi.org/10.1016/S0266-3538(03)00...
. 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)DONG, H., Z. LI, J. WANG, and B. L. KARIHALOO. 2016. “A New Fatigue Failure Theory for Multidirectional Fiber-Reinforced Composite Laminates with Arbitrary Stacking Sequence.” International Journal of Fatigue 87:294-300. doi: 10.1016/j.ijfatigue.2016.02.012.
https://doi.org/10.1016/j.ijfatigue.2016...
, who applied a modified version to analyze fatigue under uniaxial loading. Brunbauer, Gaier and Pinter (2015)BRUNBAUER, J., C. GAIER, and G. PINTER. 2015. “Computational fatigue life prediction of continuously fibre reinforced multiaxial composites.” Composites Part B: Engineering 80:269-277. doi: 10.1016/j.compositesb.2015.06.002.
https://doi.org/10.1016/j.compositesb.20...
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 efficient reproduction of mechanical behavior of this type of composite structures is relevant not only for academy, but also for industry.

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

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. 2015FREITAS, A., D. A. F. TORRES, P. T. R. MENDONÇA, and C. S. BARCELLOS. 2015. “Comparative analysis of Ck-and C0-GFEM applied to two-dimensional problems of confined plasticity.” Latin American Journal of Solids and Structures 12 (5):861-882.). 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)SUKUMAR, N., D. J. SROLOVITZ, T. J. BAKER, and J. H. PREVOST. 2003. “Brittle Fracture in Polycrystalline Microstructures with the Extended Finite Element Method.” International Journal for Numerical Methods in Engineering 56 (14):2015-2037. doi: 10.1002/nme.653
https://doi.org/10.1002/nme.653...
managed to integrate routines to DynELA for dynamic elastic XFEM. Later, following the same research line, Giner et al. (2009)GINER, E., N. SUKUMAR, J. E. TARANCON, and F. J. FUENMAYOR. 2009. “An Abaqus Implementation of the Extended Finite Element Method.” Engineering Fracture Mechanics 76 (3):347-368. doi: 10.1016/j.engfracmech.2008.10.015.
https://doi.org/10.1016/j.engfracmech.20...
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)SHI, J., J. LUA, H. WAISMAN, P. LIU, T. BELYTSCHKO, N. SUKUMAR, and Y. LIANG. 2008. “X-FEM Toolkit for Automated Crack Onset and Growth Prediction.” In 49th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference. American Institute of Aeronautics and Astronautics. 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)SHI, J., D. CHOPP, J. LUA, N. SUKUMAR, and T. BELYTSCHKO. 2010. “Abaqus Implementation of Extended Finite Element Method Using a Level Set Representation for Three-Dimensional Fatigue Crack Growth and Life Predictions.” Engineering Fracture Mechanics 77 (14):2840-2863. doi: 10.1016/j.engfracmech.2010.06.009.
https://doi.org/10.1016/j.engfracmech.20...
, 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)DEVECI, H. A., L. AYDIN, and H. S. ARTEM. 2016. “Buckling Optimization of Composite Laminates Using a Hybrid Algorithm Under Puck Failure Criterion Constraint.” Journal of Reinforced Plastics and Composites 35 (16):1233-1247. doi: 10.1177/0731684416646860.
https://doi.org/10.1177/0731684416646860...
coupled the theory as constraint for buckling analysis in a hybrid approach. Lee et al. (2015)LEE, C-S, J-H KIM, S-K KIM, D-M RYU, and J-M LEE. 2015. “Initial and Progressive Failure Analyses for Composite Laminates Using Puck failure criterion and Damage-Coupled Finite Element Method.” Composite Structures 121:406-419. doi: 10.1016/j.compstruct.2014.11.011.
https://doi.org/10.1016/j.compstruct.201...
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, 2014HOCHARD, C., S. MIOT, and Y. THOLLON. 2014. “Fatigue of laminated composite structures with stress concentrations.” Composites Part B: Engineering 65:11-16. doi: https://doi.org/10.1016/j.compositesb.2013.10.020.
https://doi.org/10.1016/j.compositesb.20...
) 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 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.

2 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 2010DEUSCHLE, M. H. 2010. “3d Failure Analysis of UD Fibre Reinforced Composites: Puck's Theory Within FEA.” Doctor in Engineering, Institute for Statics and Dynamics of Aerospace Structures - University of Stuttgart, Universitat Stuttgart (3-930683-99-7).). 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)MIOT, S., C. HOCHARD, and N. LAHELLEC. 2010. “A Non-Local Criterion for Modelling Unbalanced Woven Ply Laminates with Stress Concentrations.” Composite Structures 92 (7):1574-1580. doi: 10.1016/j.compstruct.2009.11.019.
https://doi.org/10.1016/j.compstruct.200...
. 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)ANGELO, M. V., A. C. GALUCIO, and V. TITA. 2012a. “Development and Validation of a PFA Model for Composite Aeronautical Structures with Stress Concentrators.” X World Conference on Computational Mechanics, São Paulo - Brasil. and Angelo, Galucio, and Tita (2012b)ANGELO, M. V., A. C. GALUCIO, and V. TITA. 2012b. “Parametric Analysis of Puck-Matzenmiller Theory based Damage Model for Composite Structures.” International Journal of Vehicle Structures & Systems 4 (4):152-157. doi: 10.4273/ijvss.4.4.07.
https://doi.org/10.4273/ijvss.4.4.07...
.

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 AbaqusTM, 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)ANGELO, M. V., A. C. GALUCIO, and V. TITA. 2012a. “Development and Validation of a PFA Model for Composite Aeronautical Structures with Stress Concentrators.” X World Conference on Computational Mechanics, São Paulo - Brasil.. 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 AbaqusTM 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.

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

2.1 Crack Initiation

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 2010DEUSCHLE, M. H. 2010. “3d Failure Analysis of UD Fibre Reinforced Composites: Puck's Theory Within FEA.” Doctor in Engineering, Institute for Statics and Dynamics of Aerospace Structures - University of Stuttgart, Universitat Stuttgart (3-930683-99-7).). Works have been done on optimizing this search such as Szlosarek et al. (2016)SZLOSAREK, R., T. KARALL, C. HAHNE, A. BERGER, and N. ENZINGER. 2016. “Fracture Angle Search with Puck's 3D Interfiber Fracture Criterion Using the Damped Newton's Method.” Composites: Mechanics, Computations, Applications: An International Journal 7 (4):319-339. doi: 10.1615/CompMechComputApplIntJ.v7.i4.40.
https://doi.org/10.1615/CompMechComputAp...
, but a definitive method is still an open problem. An example of this plane is shown in Figure 3.

Figure 3
Puck’s theory: Action plane (WIEGAND, PETRINIC, and ELLIOTT 2008WIEGAND, J., N PETRINIC, and B. ELLIOTT. 2008. “An Algorithm for Determination of the Fracture Angle for the Three-Dimensional Puck Matrix Failure Criterion for UD Composites.” Composites Science and Technology 68 (12):2511-2517. doi: 10.1016/j.compscitech.2008.05.004.
https://doi.org/10.1016/j.compscitech.20...
)

The action plane inclination in 2D problems can be defined analytically (PUCK, KOPP and KNOPS 2002PUCK, A., and H. SCHÜRMANN. 2002. “Failure Analysis of FRP Laminates by Means of Physically Based Phenomenological Models.” Composites Science and Technology 62 (12):1633-1662. doi: 10.1016/S0266-3538(01)00208-1.
https://doi.org/10.1016/S0266-3538(01)00...
), 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. 2011DASSAULT SYSTÈMES SIMULIA CORP. 2011. Abaqus Version 6.11 Documentation. Dassault Systems.).

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.1 Action 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.

IFF A = [ ( 1 Y t p t * S 12 ) σ 22 ] 2 + ( σ 12 S 12 ) 2 + p t * S 12 σ 22 (1)

η E / G = 1 η r E / G * 1 + c * E / G ( IFF A 1 ) ξ E / G * + η r E / G * (2)

E 22 = η E E 22 G 12 = η G G 12 (3)

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.

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

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.

2.1.2 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)WIEGAND, J., N PETRINIC, and B. ELLIOTT. 2008. “An Algorithm for Determination of the Fracture Angle for the Three-Dimensional Puck Matrix Failure Criterion for UD Composites.” Composites Science and Technology 68 (12):2511-2517. doi: 10.1016/j.compscitech.2008.05.004.
https://doi.org/10.1016/j.compscitech.20...
, 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):

φ = 1 + 5 2 (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 102. This value was taken since after 1/100 of degree, the IFF parameter is not affected.

Figure 5
Flowchart for implementing Golden Section search.

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)SCHIRMAIER, F. J., J. WEILAND, L. KÄRGER, and F. HENNING. 2014. “A New Efficient and Reliable Algorithm to Determine the Fracture Angle for Puck’s 3D Matrix Failure Criterion for UD Composites.” Composites Science and Technology 100:19-25. doi: 10.1016/j.compscitech.2014.05.033.
https://doi.org/10.1016/j.compscitech.20...
, 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.

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

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.

Table 1
Performance comparison.

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 2010TSAI, C. H., J. KOLIBAL, and M. LI. 2010. “The golden section search algorithm for finding a good shape parameter for meshless collocation methods.” Engineering Analysis with Boundary Elements 34 (8):738-746. doi: 10.1016/j.enganabound.2010.03.003.
https://doi.org/10.1016/j.enganabound.20...
), medical appliance (LU and HARGREAVES 2008LU, W., and B. A. HARGREAVES. 2008. “Multiresolution field map estimation using golden section search for water-fat separation.” Magnetic Resonance in Medicine 60 (1):236-244. doi: 10.1002/mrm.21544.
https://doi.org/10.1002/mrm.21544...
) and even aircraft control (PATRON, BOTEZ, and LABOUR 2012PATRON, R. F., R. BOTEZ, and D. LABOUR. 2012. “Vertical profile optimization for the Flight Management System CMA-9000 using the golden section search method.” IECON 2012 - 38th Annual Conference, Montreal - Canada.).

3 RESULTS AND DISCUSSION

As commented previously, there were preliminary studies performed by industry (classified results) and academy (HOCHARD, MIOT and THOLLON, 2014HOCHARD, C., S. MIOT, and Y. THOLLON. 2014. “Fatigue of laminated composite structures with stress concentrations.” Composites Part B: Engineering 65:11-16. doi: https://doi.org/10.1016/j.compositesb.2013.10.020.
https://doi.org/10.1016/j.compositesb.20...
), 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.

Table 2
Material properties for case studies.

3.1 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)PUCK, A., and H. SCHÜRMANN. 1998. “Failure Analysis of FRP Laminates by Means of Physically Based Phenomenological Models.” Composites Science and Technology 58 (7):1045-1067. doi: 10.1016/S0266-3538(96)00140-6.
https://doi.org/10.1016/S0266-3538(96)00...
, it is expected that the angle for failure occur around 54° (DEUSCHLE 2010DEUSCHLE, M. H. 2010. “3d Failure Analysis of UD Fibre Reinforced Composites: Puck's Theory Within FEA.” Doctor in Engineering, Institute for Statics and Dynamics of Aerospace Structures - University of Stuttgart, Universitat Stuttgart (3-930683-99-7).). 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 2002PUCK, A., J. KOPP, and M. KNOPS. 2002. “Guidelines for the Determination of the Parameters in Puck’s Action Plane Strength Criterion.” Composites Science and Technology 62 (3):371-378. doi: 10.1016/S0266-3538(01)00202-0.
https://doi.org/10.1016/S0266-3538(01)00...
). 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.

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

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

Figure 8
Failure index for each angle: Compression analyses.

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

Figure 9
Failure index for each angle: Tension analyses.

Table 3 shows comparison the current method to Puck’s method for the Case Study 1.

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

3.2 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 2015AMERICAN SOCIETY FOR TESTING AND MATERIALS. 2015. E1922-04(2015). In Standard Test Method for Translaminar Fracture Toughness of Laminated and Pultruded Polymer Matrix Composite Materials. Philadelphia: American Society for Testing Materials.) and it has been commonly used for investigating translaminar failure and crack propagation in laminated composites (LAFFAN et al. 2010LAFFAN, M. J., S. T. PINHO, P. ROBINSON, and L. IANNUCCI. 2010. “Measurement of the in situ Ply Fracture Toughness Associated with mode I fibre Tensile Failure in FRP. Part II: Size and lay-up effects.” Composites Science and Technology 70 (4):614-621. doi: 10.1016/j.compscitech.2009.12.011.
https://doi.org/10.1016/j.compscitech.20...
, LAFFAN et al. 2012LAFFAN, M. J., S. T. PINHO, P. ROBINSON, and A. J. MCMILLAN. 2012. “Translaminar Fracture Toughness Testing of Composites: A review.” Polymer Testing 31 (3):481-489. doi: 10.1016/j.polymertesting.2012.01.002.
https://doi.org/10.1016/j.polymertesting...
, TEIXEIRA, PINHO, and ROBINSON 2011TEIXEIRA, R. F., S. T. PINHO, and P. ROBINSON. 2011. “Translaminar Ply Fracture Toughness of Advanced Composites.” International Conference on Composite Materials, Jeju, 2011., FERNANDES 2014FERNANDES, M. F. da C. 2014. “Translaminar Fracture of Thin-Ply Composite Laminates.” Masters in Science, Faculdade de Engenharia da Universidade do Porto, Universidade do Porto.). 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.

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

Final displacement for each analysis is also presented in Table 4.

Table 4
Details of computational analyses for Case Study 2

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)JAMALI, J. 2014. “Mechanistic Failure Criterion for Unidirectional and Random Fibre Polymer Composites.” The University of Western Ontario. and Slepetz and Carlson (1976)SLEPETZ, J. M., and L. CARLSON. 1976. Fracture of composite compact tension specimens. Massachusets: Army Materials and Mechanics Research Center.. 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.

Figure 11
Compact tension experimental results: (A) 0° (JAMALI 2014JAMALI, J. 2014. “Mechanistic Failure Criterion for Unidirectional and Random Fibre Polymer Composites.” The University of Western Ontario.), (B) 90° (JAMALI 2014) and (C) 45° specimen (SLEPETZ and CARLSON 1976SLEPETZ, J. M., and L. CARLSON. 1976. Fracture of composite compact tension specimens. Massachusets: Army Materials and Mechanics Research Center.).

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 12
Compact tension model results: (A) 0°, (B) 90° and (C) 45° specimen.

4 CONCLUSIONS

The present work showed a new computational framework, combining XFEM and 3D Puck’s failure theory, for evaluating onset and propagation of failure (evolution of cracks) in composite structures. Besides, the proposed approach for determining crack inclination angle uses an enhanced golden section method, which has proven to be robust and precise, filling a gap on the available researches in the field. Compared to the algorithm proposed by Puck, the new methodology has convergence one order higher than the first method developed by Puck; and it is 20 times more efficient computationally for a 0.1° precision. Therefore, if more precision is needed, then higher gains are achieved combined to lower computational cost by using the proposed method.

Moreover, proposed framework performance for translaminar failure was evaluated through thick laminated composites with layers mainly oriented to 90o simulated under tension and compression, displaying results as reported in the literature. In addition, compact tension tests with 0°, 90° and 45° specimens were evaluated, and numerical results were qualitatively coherent with experimental data. Therefore the computational framework can be used to predict onset and crack propagation in composite structures with translaminar failure

5 ACKNOWLEDGEMENTS

The authors acknowledge the Consulat Général de France à São Paulo for supporting this Research. Volnei Tita acknowledges the financial support of the National Council for Scientific and Technological Development (CNPq process number: 401170/2014-4 and 310094/2015-1) and Air Force Office of Scientific Research - AFOSR (Grant FA9550-16-1-0222). Marcelo L. Ribeiro acknowledges the financial support of the São Paulo Research Foundation (FAPESP 2015/13844-8).

6 REFERENCES

  • AMERICAN SOCIETY FOR TESTING AND MATERIALS. 2015. E1922-04(2015). In Standard Test Method for Translaminar Fracture Toughness of Laminated and Pultruded Polymer Matrix Composite Materials. Philadelphia: American Society for Testing Materials.
  • ANGELO, M. V., A. C. GALUCIO, and V. TITA. 2012a. “Development and Validation of a PFA Model for Composite Aeronautical Structures with Stress Concentrators.” X World Conference on Computational Mechanics, São Paulo - Brasil.
  • ANGELO, M. V., A. C. GALUCIO, and V. TITA. 2012b. “Parametric Analysis of Puck-Matzenmiller Theory based Damage Model for Composite Structures.” International Journal of Vehicle Structures & Systems 4 (4):152-157. doi: 10.4273/ijvss.4.4.07.
    » https://doi.org/10.4273/ijvss.4.4.07
  • BRUNBAUER, J., C. GAIER, and G. PINTER. 2015. “Computational fatigue life prediction of continuously fibre reinforced multiaxial composites.” Composites Part B: Engineering 80:269-277. doi: 10.1016/j.compositesb.2015.06.002.
    » https://doi.org/10.1016/j.compositesb.2015.06.002.
  • DASSAULT SYSTÈMES SIMULIA CORP. 2011. Abaqus Version 6.11 Documentation. Dassault Systems.
  • DEUSCHLE, M. H. 2010. “3d Failure Analysis of UD Fibre Reinforced Composites: Puck's Theory Within FEA.” Doctor in Engineering, Institute for Statics and Dynamics of Aerospace Structures - University of Stuttgart, Universitat Stuttgart (3-930683-99-7).
  • DEVECI, H. A., L. AYDIN, and H. S. ARTEM. 2016. “Buckling Optimization of Composite Laminates Using a Hybrid Algorithm Under Puck Failure Criterion Constraint.” Journal of Reinforced Plastics and Composites 35 (16):1233-1247. doi: 10.1177/0731684416646860.
    » https://doi.org/10.1177/0731684416646860
  • DONG, H., Z. LI, J. WANG, and B. L. KARIHALOO. 2016. “A New Fatigue Failure Theory for Multidirectional Fiber-Reinforced Composite Laminates with Arbitrary Stacking Sequence.” International Journal of Fatigue 87:294-300. doi: 10.1016/j.ijfatigue.2016.02.012.
    » https://doi.org/10.1016/j.ijfatigue.2016.02.012
  • FERNANDES, M. F. da C. 2014. “Translaminar Fracture of Thin-Ply Composite Laminates.” Masters in Science, Faculdade de Engenharia da Universidade do Porto, Universidade do Porto.
  • FREITAS, A., D. A. F. TORRES, P. T. R. MENDONÇA, and C. S. BARCELLOS. 2015. “Comparative analysis of Ck-and C0-GFEM applied to two-dimensional problems of confined plasticity.” Latin American Journal of Solids and Structures 12 (5):861-882.
  • GINER, E., N. SUKUMAR, J. E. TARANCON, and F. J. FUENMAYOR. 2009. “An Abaqus Implementation of the Extended Finite Element Method.” Engineering Fracture Mechanics 76 (3):347-368. doi: 10.1016/j.engfracmech.2008.10.015.
    » https://doi.org/10.1016/j.engfracmech.2008.10.015
  • HINTON, M. J., A. S. KADDOUR, and P. D. SODEN. 2002. “A Comparison of the Predictive Capabilities of Current Failure Theories for Composite Laminates, Judged Against Experimental Evidence.” Composites Science and Technology 62 (12):1725-1797. doi: 10.1016/S0266-3538(02)00125-2.
    » https://doi.org/10.1016/S0266-3538(02)00125-2
  • HOCHARD, C., S. MIOT, and Y. THOLLON. 2014. “Fatigue of laminated composite structures with stress concentrations.” Composites Part B: Engineering 65:11-16. doi: https://doi.org/10.1016/j.compositesb.2013.10.020
    » https://doi.org/10.1016/j.compositesb.2013.10.020
  • JAMALI, J. 2014. “Mechanistic Failure Criterion for Unidirectional and Random Fibre Polymer Composites.” The University of Western Ontario.
  • KADDOUR, A. S., M. J. HINTON, and P. D. SODEN. 2004. “A Comparison of the Predictive Capabilities of Current Failure Theories for Composite Laminates: Additional Contributions.” Composites Science and Technology 64 (3):449-476. doi: 10.1016/S0266-3538(03)00226-4.
    » https://doi.org/10.1016/S0266-3538(03)00226-4
  • LAFFAN, M. J., S. T. PINHO, P. ROBINSON, and L. IANNUCCI. 2010. “Measurement of the in situ Ply Fracture Toughness Associated with mode I fibre Tensile Failure in FRP. Part II: Size and lay-up effects.” Composites Science and Technology 70 (4):614-621. doi: 10.1016/j.compscitech.2009.12.011.
    » https://doi.org/10.1016/j.compscitech.2009.12.011
  • LAFFAN, M. J., S. T. PINHO, P. ROBINSON, and A. J. MCMILLAN. 2012. “Translaminar Fracture Toughness Testing of Composites: A review.” Polymer Testing 31 (3):481-489. doi: 10.1016/j.polymertesting.2012.01.002.
    » https://doi.org/10.1016/j.polymertesting.2012.01.002
  • LEE, C-S, J-H KIM, S-K KIM, D-M RYU, and J-M LEE. 2015. “Initial and Progressive Failure Analyses for Composite Laminates Using Puck failure criterion and Damage-Coupled Finite Element Method.” Composite Structures 121:406-419. doi: 10.1016/j.compstruct.2014.11.011.
    » https://doi.org/10.1016/j.compstruct.2014.11.011
  • LU, W., and B. A. HARGREAVES. 2008. “Multiresolution field map estimation using golden section search for water-fat separation.” Magnetic Resonance in Medicine 60 (1):236-244. doi: 10.1002/mrm.21544.
    » https://doi.org/10.1002/mrm.21544
  • MIOT, S., C. HOCHARD, and N. LAHELLEC. 2010. “A Non-Local Criterion for Modelling Unbalanced Woven Ply Laminates with Stress Concentrations.” Composite Structures 92 (7):1574-1580. doi: 10.1016/j.compstruct.2009.11.019.
    » https://doi.org/10.1016/j.compstruct.2009.11.019
  • MOU, H. L., T. C. ZOU, Z. Y. FENG, and J. XIE. 2016. “Crashworthiness Analysis and Evaluation of Fuselage Section with Sub-floor Composite Sinusoidal Specimens.” Latin American Journal of Solids and Structures 13 (6):1186-1202.
  • PATRON, R. F., R. BOTEZ, and D. LABOUR. 2012. “Vertical profile optimization for the Flight Management System CMA-9000 using the golden section search method.” IECON 2012 - 38th Annual Conference, Montreal - Canada.
  • PUCK, A., J. KOPP, and M. KNOPS. 2002. “Guidelines for the Determination of the Parameters in Puck’s Action Plane Strength Criterion.” Composites Science and Technology 62 (3):371-378. doi: 10.1016/S0266-3538(01)00202-0.
    » https://doi.org/10.1016/S0266-3538(01)00202-0
  • PUCK, A., and H. SCHÜRMANN. 1998. “Failure Analysis of FRP Laminates by Means of Physically Based Phenomenological Models.” Composites Science and Technology 58 (7):1045-1067. doi: 10.1016/S0266-3538(96)00140-6.
    » https://doi.org/10.1016/S0266-3538(96)00140-6
  • PUCK, A., and H. SCHÜRMANN. 2002. “Failure Analysis of FRP Laminates by Means of Physically Based Phenomenological Models.” Composites Science and Technology 62 (12):1633-1662. doi: 10.1016/S0266-3538(01)00208-1.
    » https://doi.org/10.1016/S0266-3538(01)00208-1
  • SCHIRMAIER, F. J., J. WEILAND, L. KÄRGER, and F. HENNING. 2014. “A New Efficient and Reliable Algorithm to Determine the Fracture Angle for Puck’s 3D Matrix Failure Criterion for UD Composites.” Composites Science and Technology 100:19-25. doi: 10.1016/j.compscitech.2014.05.033.
    » https://doi.org/10.1016/j.compscitech.2014.05.033
  • SHI, J., D. CHOPP, J. LUA, N. SUKUMAR, and T. BELYTSCHKO. 2010. “Abaqus Implementation of Extended Finite Element Method Using a Level Set Representation for Three-Dimensional Fatigue Crack Growth and Life Predictions.” Engineering Fracture Mechanics 77 (14):2840-2863. doi: 10.1016/j.engfracmech.2010.06.009.
    » https://doi.org/10.1016/j.engfracmech.2010.06.009
  • SHI, J., J. LUA, H. WAISMAN, P. LIU, T. BELYTSCHKO, N. SUKUMAR, and Y. LIANG. 2008. “X-FEM Toolkit for Automated Crack Onset and Growth Prediction.” In 49th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference. American Institute of Aeronautics and Astronautics.
  • SLEPETZ, J. M., and L. CARLSON. 1976. Fracture of composite compact tension specimens. Massachusets: Army Materials and Mechanics Research Center.
  • SUKUMAR, N., D. J. SROLOVITZ, T. J. BAKER, and J. H. PREVOST. 2003. “Brittle Fracture in Polycrystalline Microstructures with the Extended Finite Element Method.” International Journal for Numerical Methods in Engineering 56 (14):2015-2037. doi: 10.1002/nme.653
    » https://doi.org/10.1002/nme.653
  • SZLOSAREK, R., T. KARALL, C. HAHNE, A. BERGER, and N. ENZINGER. 2016. “Fracture Angle Search with Puck's 3D Interfiber Fracture Criterion Using the Damped Newton's Method.” Composites: Mechanics, Computations, Applications: An International Journal 7 (4):319-339. doi: 10.1615/CompMechComputApplIntJ.v7.i4.40.
    » https://doi.org/10.1615/CompMechComputApplIntJ.v7.i4.40
  • TEIXEIRA, R. F., S. T. PINHO, and P. ROBINSON. 2011. “Translaminar Ply Fracture Toughness of Advanced Composites.” International Conference on Composite Materials, Jeju, 2011.
  • TSAI, C. H., J. KOLIBAL, and M. LI. 2010. “The golden section search algorithm for finding a good shape parameter for meshless collocation methods.” Engineering Analysis with Boundary Elements 34 (8):738-746. doi: 10.1016/j.enganabound.2010.03.003.
    » https://doi.org/10.1016/j.enganabound.2010.03.003
  • WIEGAND, J., N PETRINIC, and B. ELLIOTT. 2008. “An Algorithm for Determination of the Fracture Angle for the Three-Dimensional Puck Matrix Failure Criterion for UD Composites.” Composites Science and Technology 68 (12):2511-2517. doi: 10.1016/j.compscitech.2008.05.004.
    » https://doi.org/10.1016/j.compscitech.2008.05.004
  • Available Online: May 25, 2018

Publication Dates

  • Publication in this collection
    2018

History

  • Received
    23 July 2017
  • Reviewed
    11 Apr 2018
  • Accepted
    05 May 2018
Individual owner www.lajss.org - São Paulo - SP - Brazil
E-mail: lajsssecretary@gmsie.usp.br