Determination of Paris parameters from S−N curve in the simulation of fatigue in Nomex honeycomb sandwich structures

Determination of S-N and Paris parameters form the two primary approaches as far as fatigue testing is concerned. With many economic factors and complexities involved in experimentation, these approaches are not related to each other and require implementation of separate test standards for both of them. In this work, an attempt is made on the interconvertibility of these experimental parameters and a numerical technique based on the virtual crack closure technique (VCCT) is proposed for the evaluation of Paris parameters from S−N curve, corresponding to four-point bending end notched flexural (4ENF) sandwich specimen made up of Nomex® aramid honeycomb core with glass fibre reinforced composite (GFRP) facesheets. VCCT was effectively utilised for the estimation of strain energy release rates (SERR), in both 2D plane stress/strain and 3D conditions, which were compared with each other, and the resultant was incorporated with the S−N data for the determination of parameters of the Paris curve. These determined parameters were employed for the simulation of fatigue delamination growth by using ‘Direct Cyclic’ step and compared with the results obtained from the literature. The comparison revealed a good agreement between the two, thereby validating the proposed technique.


1. Introduction
Fatigue life estimation is one of the principal tools in the fail-safe design of mechanical structures, in particular aerospace structures which are quite susceptible to failure due to cyclic loads. This fatigue life estimation is based on predominantly two concepts namely, cumulative fatigue damage (CFD), based on the stress-life approach i.e. S−N curves, and the fatigue crack growth (FCG) based on the advancement of fatigue crack with respect to the number of cycles, often referred to as a-N or Paris curves [1][2][3].
The stress-life approach is quite popular and is widely used to insight into the fatigue behaviour of structural components. This approach utilises a mathematical relation between the sustained range of stresses and the number of reported load cycles until failure or damage occurs in the structure [4]. Generally, the curve produced through such kind of mathematical relation is known as S−N curve, where S is the stress range of constant amplitude and is represented as a function of the number of load cycles to damage or failure, N, as seen in figure 1. This approach is very helpful for design engineers in order to explore the fatigue failure or damage behaviour of any structure, however, it does not provide sufficient information related to the onset or the propagation of fatigue cracks [5,6].
The fatigue crack growth (FCG) approach, on the other hand, utilises the correlation between fatigue and fracture, by developing material and fracture mode-specific dependent mathematical relations between fatigue crack growth rate (FCGR) and the reported number of load cycles, till the failure occurs. The curve so produced, referred to as the Paris Curve [7], depicts the rate of advancement of the crack / da dN in comparison to the change in strain energy release rate (SERR) as seen in figure 2.
Paris and Erdogan pioneered a technique to compare Irwin's stress intensity factor range with the rate of growth of crack per cycle [7] with the development of mathematical relations between the growth rate of cracks per cycle ( / da dN ) and the relative strain energy release rate (SERR), known as the Paris law as presented in the equation (1).
Where, C and m are Paris parameters and can be determined by the linear regression curve fit of experimental data and it can vary with stress ratio (R) and mode-mixity. With the final fracturing of a component, depicting the final failure under fatigue loading, therefore the S−N curves and the Paris curves reflect the same basic behavioural indication of the material or structure under investigation [8]. Although, both curves have their own methodologies involved, but there must exist some kind of relationship between their associated parameters [9]. Therefore, based on a representative curve of fatigue crack growth rate and the basic formats of an S−N curve, some mathematical relationships can be established between these two curves. This will facilitate the interconvertibility of test data in such a way that only one out of two types of curves are needed to be tested and the other one can be estimated from the existing or available test results, thereby reducing the financial burden of extra testing [10] and [11].
Among the pioneering works in this direction, a novel cohesive damage model for cyclic fatigue, based on idealisations of S−N data, was proposed by Carlos G. Davila in [12], both in mode-I and mixed-mode conditions using double cantilever beam (DCB) and mixed-mode bending (MMB) specimen respectively. Based on the assumption that the density of microcracks is represented by both stable tearing damage and damage due to cyclic loading, which thereby can be clubbed in a single damage variable can describe the state of damage and therefore, a tearing-based quasi-static cohesive law can be taken as the envelope of the fatigue damage.
Further, the difficulty in the implementation of Paris law, using cohesive elements, was identified, which is primarily due to the presence of a process zone, as Paris law is not a local measure of fatigue propagation rate. A  UMAT subroutine-based damage model was thereby developed using Abaqus ® , by adding fatigue damage accumulation based on Turon's [13] mixed-mode cohesive laws. The model demonstrated the capability of replicating the S−N curve which is exclusively based on the loading history and damage accumulation at the point, by predicting the crack propagation rate rather than what is obtained by Paris law, the algorithm of which is schematically summarised in figure 3. For the mixed-mode MMB specimen, the validation was made using the test data reported by Ratliffe [14], by fitting a quadratic polynomial for different crack lengths due to the unavailability of a closed-form expression for compliance calculations in MMB.
The present work is based on the four-point bending end notched flexural (4ENF) specimen of a honeycomb sandwich composite. For determining the mode-II fracture toughness of unidirectional laminated composites, 4ENF test was first described by Martin and Davidson in [15], and the test configuration was further analysed by beam theory and by finite element analysis. The study revealed the major advantages of the 4ENF test over other mode-II tests predominantly the promotion of stable crack propagation, reduction of frictional effects, and has the ability to perform experimental compliance calibration as the delamination advances. Further, a detailed comparative analysis between end-notched flexural (ENF) and 4ENF tests for the estimation of mode-II fracture toughness of bonded aluminium joints was conducted by De Oliveira in [16], where different adhesives and data  reduction techniques were utilised to compare the load-displacement curves and resistance curves (R-curves). the study concluded that the 4ENF test was more accurate and stable for the estimation of mode-II fracture toughness.

2. Fatigue crack growth in composites
Fatigue crack growth in composite and sandwich materials gained much importance, particularly in the last couple of decades. Among the pioneering works, failure due to fatigue, at loads much below the quasi-static strength and safe-fail tolerant approach was demonstrated in composite structures [14] and [15]. The analysis of fatigue crack growth is based on the evaluation of fracture mechanics parameters such as J-integral, stress intensity factor (SIF) or SERR, and crack growth per cycle ( / da dN ) under cyclic fatigue loading [17][18][19]. Such parameters were plotted on the log-log scale and a sigmoidal or S-shaped curve was achieved to study the FCGR for both isotropic as well as anisotropic materials [3].
Initiation and propagation of fatigue crack in metal to composite bonded joints, through experimental and numerical approaches, was done in [20] and fatigue life was predicted corresponding to the interfacial crack propagation [21][22][23]. FCG was analysed through the development of a flexbeam-shaped composite helicopter rotor hub and a review on fibre metal laminates 'Glare' were presented in [24] and [25] respectively. Similar work on crack propagation and crack front shape corresponding to repaired and unrepaired panels was reported in [26]. Fatigue behaviour of monolithic and sandwich composites was investigated in [27] and [28] respectively. Various approaches like fracture mechanics-based, strain or stress based and cohesive-based were comprehensively reviewed in [29] with mode-I fracture-based tests on the specimen like double cantilever beam (DCB) and tapered double cantilever beam (TDCB) in [19,30,31].
Interfacial fatigue behaviour of both damaged and undamaged honeycomb sandwich composite structures was investigated under a four-point bending mechanism [32]. The study revealed that the undamaged structure failed due to the compressive collapse of facesheets, whereas the damaged structure failed due to the collapse of honeycomb cell walls at the tip of the unbonded region. Analytical and experimental approaches were coupled in Table 1. Features and limitations of VCCT.

Pros Cons
Calculates crack growth based on principles of linear elastic fracture mechanics by calculating nodal displacements and forces.
Not applicable to ductile fractures as it lacks the mechanism to consider deformations due to plasticity. One of the state-of-the-art techniques for the evaluation of SERR.
Requires pre-crack and crack path definitions. Requires very few inputs like fracture toughness, besides regular material parameters for the evaluation of SERR.
Mesh sensitivity is relatively high and requires finer mesh at the interfaces. Capable of being directly implemented in fatigue loading by means of numerical integration of Paris law. Supports modelling of delamination/debonding at interfaces between bonded surfaces to allow evaluation of structural strength and failure modes. Capable of simulating crack onset as well as propagation without re-meshing.

Pros Cons
Calculates crack growth based on principles of damage mechanics by means of constitutive traction separation laws.
Requires crack path definitions.
Applicable on both ductile as well as brittle fractures. Generally, SERR is not evaluated directly and requires the implementation of subroutines for the estimation of SERR from damage parameter. Do not require pre-crack definition.
Implementation in fatigue loading is limited and requires dedicated subroutines. Supports modelling of delamination/debonding at interfaces between bonded surfaces to allow evaluation of structural strength and failure modes Determination of constitutive traction separation laws is sometimes a tedious process depending upon the nature of the material. Capable of SERR evaluations as well as simulation of crack propagation.
[33], with the aim of investigating FCG in a DCB specimen under a temperature-specific environment, demonstrating the higher fracture toughness and lower FCGR under a colder environment.

Virtual crack closure technique (VCCT)
Based on Irwin's crack closure concept, which assumes that strain energy released due to the opening of a crack is equal to the strain energy required for the closure of the crack [34], VCCT is a LEFM based technique utilising nodal displacements of the opened surface prior to the crack tip and the reaction force at the crack tip, for the calculation of the SERR, which is thereby employed to estimate both the crack onset and crack propagation for static as well as fatigue problems. Typical crack advancement and nodal displacement schematics have been presented in figure 4. Initially proposed by Rybicki and Kanninen [35], VCCT emerged as an effective technique for the determination of SERR by evaluating the energy per unit area required to close a small increment of crack length.   The oscillatory behaviour of SERR was decoded through the stress zone in the vicinity of the crack-tip for the bidirectional specimen with the help of VCCT [36]. Since fibre composites and sandwich materials exhibit brittle fracture based on LEFM theory, therefore over years VCCT has been extensively employed in simulating failures in these materials and has become quite popular because of its ease of use and ability to classify SERR into the various fracture modes [37]. This was implemented for various materials, orientations, structural components and fracture modes in [38] - [39]. A comprehensive report about the history, approach, and applications of VCCT has been presented by Krueger [40], which was further implemented by many researchers for modelling delamination propagation in 2D as well as 3D simulations in [41]- [42].
Modelling strategies of VCCT, CZM and extended finite element method (XFEM) approaches in Abaqus ® were provided in [43], and a coupled approach utilizing XFEM and VCCT was proposed in [44], for the simulation of mode-I fatigue delamination in composites. The crack propagation behaviour of a post-buckled symmetrical and unsymmetrical laminated composite under compressive loading was studied in [45] through implicit finite element analysis (FEA) using VCCT. The FEM-VCCT and XFEM with phantom nodes (XFEMPN VCCT) algorithms were presented in [46] for the re-meshing-free FCG simulation and life estimation in Abaqus ® . Further, a Python script-based 'adaptive VCCT' (AVCCT) was introduced and implemented in Abaqus ® in the form of a new XFEMPN-AVCCT algorithm and concluded that the adaptive VCCT can significantly enhance the accuracy of FCG simulation and life estimation. As seen VCCT is a widely implemented technique, the features and limitations of which are presented in table 1.

Cohesive zone method (CZM)
Initially proposed by D.S. Dugdale and G.I. Barenblatt in the early sixties [47], cohesive zone modelling (CZM) proved to be an effective and widely used method to examine damage under quasi-static and fatigue loading by means of a constitutive law which is often referred to as traction separation law (TSL). For the investigation of interfacial debonding, the energy potential function was utilised in [48] through the concept of CZM. The method has the advantage of simulating failure not only based on LEFM but also failure pertaining to elasticplastic fracture mechanics (EPFM), along with the independence of the existence of a pre-existing crack in the simulation. The pros, cons and predictive capability of the CZM were reviewed when applied to different materials such as concrete and steel [49]. Some exemplary works include investigation of the numerical simulation of crack growth based on decohesion elements placed between the layers of solid elements, such that the crack probably extends in laminated composite material through CZM [50]. Two-dimensional and three-dimensional fatigue modelling using CZM for the simulation of FCG in composite materials was performed for mode-I, mode-II, and mixedmode fracture on DCB, end-loaded split (ELS), mixed mode end notched flexural (MMENF) specimens  respectively [51][52][53][54][55]. Turon et al. and Kawashita-Hallett in [56], simulated mode-I fatigue damage model with a trilinear cohesive law, constructed by the superposition of two bilinear CZMs. The models were implemented in Abaqus ® using a user-defined element (UEL) subroutine and FE analyses were conducted under high-cycle fatigue (HCF) loading for the DCB specimens. The study revealed that the trilinear CZM's have more accuracy

3. Problem statement
Based on the literature excerpts, particularly related to sandwich materials [58][59][60], this study is focused on the prediction of Paris law parameters from the S−N curve of a honeycomb sandwich composite structure for which a finite element analysis (FEA) based numerical approach has been proposed thereby. Four-point bending end notched flexural (4ENF) specimen has been taken into consideration employing the VCCT in Abaqus ® for the evaluation of SERR which ultimately leads to the determination of Paris law parameters (C and m). The 4ENF specimen comprises of face-sheet made up of glass fibre reinforcement plastic (GFRP) of thickness 0.72 mm of each and Nomex aramid honeycomb core of thickness 11.26 mm, as seen in figure 5 along with boundary conditions. Geometrical dimensions are represented in table 3, and the mechanical properties of face-sheet material and honeycomb core material, that have been used in this study are respectively presented in tables 4 and 5.
A quasi-static load corresponding to 516 N was taken as per the load-displacement curve obtained from the maximum envelope of the hysteresis curve presented in [57], from which the maximum fatigue load (P max ) was also allocated for the fatigue analysis presented in the later part of the work.

4. Methodology
The main methodology involved utilises the S−N data, from which the number of cycles is extracted. At the same time, a numerical FEA model based on static VCCT in Abaqus ® is implemented for both two and three dimensions, utilising the material and geometric properties and peak fatigue load from experimental data. The two-dimensional model is made using VCCT for the evaluation of SERR with the known increments of crack length a, for both plane stress and plane strain conditions. This is compared with the three-dimensional model, revealing the relationship between the actual three-dimensional model with the idealistic two-dimensional model. These values are clubbed with the number of cycles obtained from the S−N curve, with the initial crack corresponding to the initial no. of cycles and the final one with the final no. of cycles.
Relative values of crack length, no. of cycles and SERR are taken into consideration and the rate of the advancement of crack with respect to no. of cycles is evaluated. Finally, Δa/ΔN is plotted against ΔG on a loglog scale and thereby the values of Paris parameters C and m are evaluated. The whole methodology is schematically depicted in detail in figure 6.

5. Modelling and simulation
In compliance with the proposed methodology, load-displacement and S−N data have been extracted from reported data in [57], as depicted in figures 7 and 8 respectively.

Static model
Two-dimensional (2D) and three-dimensional (3D) quasi-static finite element modelling of the 4ENF specimen have been implemented in Abaqus ® . Loading and boundary conditions, corresponding to the specimen, have been provided by means of two roller constraints along the bottom facesheet and two roller loads on the top facesheet, as seen in figure 9(a). The modelling has been specifically focused on the LEFM based VCCT method, involving a 5 mm pre-crack provided as the unbonded region between the interfaces of the bottom surface of the facesheet top (FST) and the top surface of the honeycomb (HC). The remaining portion between FST and HC has been bonded using bond nodes as seen in figure 9(b). Tie constraints have been used between the bottom surface of the honeycomb and the top surface of the facesheet bottom (FSB) as shown in figure 9(c). Node-to-surface interaction has been established between the bonded region of FST and HC, and hard contact normal behaviour along with tangential behaviour with friction formulation coefficient having a value of 0.8 [61,62]. Besides this, the unbonded region is again provided with a 'Penalty' constraint enforcement method to prevent any penetration of the facesheet in the honeycomb core. Facesheets and the honeycomb core have been modelled with continuum shell elements (SC8R) and 3D-Stress elements (C3D8R) respectively with reduced integration and hourglass control, as per Abaqus ® documentation [63]. The final mesh has been shown in figure 9(d), and the details of the meshing strategy have been given in table 6.  Once the results of these were obtained, they were incorporated with the available S−N curve data for the evaluation of Paris parameters which was further utilised in the simulation of fatigue crack propagation.

Fatigue model
Direct cyclic low fatigue step has been used to model the fatigue crack propagation of 4ENFspecimen and table 7 prescribed the parameters that have been used in the direct cyclic step through step module.
Sinusoidal periodic amplitude has been defined through the load module and 2 Hz cyclic load frequency as well as load ratio (R) = 0.1 has been also provided as per the literature. The parameters of amplitude definition that have been used in modelling are given in table 8.
Once the FEA modelling is completed, the required data have been extracted from the post-processing module of the FEA results to predict the required fatigue crack growth rate curve or Paris curve. This was finally validated by comparison of these FEA results ( / da dN versus DG) with the results obtained from the literature [57].

6. Results and discussion
Corresponding to the modelling strategies, numerical results were obtained using the Abaqus ® standard involving VCCT in 'Static General' for the determination of SERR with respect to the increment in crack length a. Paris parameters were thereby evaluated with this data and were further implemented in the simulation of fatigue crack growth involving VCCT in 'Direct Cyclic', revealing the advancement of crack length with the no. of cycles.

Static results
The relevant FEA results visualisation i.e., displacement distribution of specimen, SERR distribution of specimen and bond status between the bottom surface of top facesheet and top surface of honeycomb were evaluated for both 2D and 3D models. The distribution of displacements on the 4ENF specimen are depicted in figures 10(a)-(b). Displacement distribution in the direction of crack propagation (x-direction) can be seen in figure 10(a) depicting predominant Mode-II fracture behaviour due to the longitudinal movement of the  facesheet under the application of the bending load. Whereas Mode-I behaviour has been presented in figure 10(b) depicting the displacement distribution in the direction of flexural load (z-direction), together contributing to mixed-mode conditions of traditional 4ENF system. The SERR for 2D and 3D models of 4ENF specimen through the above-described FEA modelling technique have been presented in figures 11(a)-(c). The distribution of SERR between the interface of facesheet and honeycomb after the application of load have been presented for both 2D and 3D models with respect to crack propagation. The bond status between the interface of the top facesheet (FST) and the top surface of the honeycomb (HC) can be visualised in figures 12(a)-(c), depicting how the interface of FST and HC are getting deboned with respect to crack propagation. Here the blue colour represents the debonding after the application load and yellowish-red colour represents the bonded region (after the onset of crack propagation) between the interface of FST and HC, and the grey colour between the interface of FST and HC represents unbonded (precrack) region.
The SERR values so obtained for the plane stress (2DPS), and plane strain (2DPE) can be seen in figure 13. The SERR values so obtained for 3D solid stress element are in consonance with the values obtained for the 2D model corresponding to plane strain conditions.

Evaluation of Paris parameters
In the static FEA modelling, crack length is measured at the points corresponding to the maximum SERR for each increment of the crack, and for each cyclic increment in case of fatigue analysis. This is done by considering the field output variable 'CNAREA', which gives an insight into the contact area associated with each node in contact. The numerical difference of total area with this 'CNAREA' gives the magnitude of crack length corresponding to each increment. In consonance with results obtained in figure 13, the crack propagation rates obtained by the 3D model were also compared with the 2D models for both plane stress and plane strain  conditions against relative SERR on a log-log scale for the estimation of Paris parameters, as seen in figure 14. It was again found that 3D crack estimation was comparing well with the plane strain conditions.
Once the equation (1)      dominance of shear stress over the normal loads which are primarily dominant over the facesheet and the honeycomb core respectively. Direct cyclic procedure, adopted of low cycle fatigue (LCF), demonstrated a progressive release of bonded nodes per increment and the SERR values were evaluated, which were thereby utilised in the procedure, for the numerical integration of the Paris law. This progressive release of the bonded nodes can be visualised by means of contact status, at the different numbers of cycles as seen in figure 16.
To further elaborate on this, fatigue crack propagation along the bonded top face of a unit cell, at the crack front, in the honeycomb-facesheet interface is presented in figure 17, giving an insight into the way how the crack propagates in a bonded honeycomb interface with a monolithic composite facesheet. The crack approaching unit cell branches along the cell geometry, as seen in figures 17(a)-(b). It further merges with the crack of the adjacent cell and propagates along the edge which can be visualised in figures 17(c)-(d), followed up by approach and subsequent branching to the next unit cell, as depicted in figures 17(e)-(f ). The crack length is thereby estimated by manually taking the mean of all the opened surfaces of all the unit cells at the crack front, which was required for the establishment of results presented later.
The SERR variations with respect to crack length a and number of cycles N are presented in figures 18, and 19, respectively. The SERR values depict a gradually increasing trend till failure is achieved somewhere around 36 mm of crack length, evident from the sharp rise in the SERR values after this point. A similar trend is also seen  in figure 19, in which the failure is seen at N = 9.0 × 10 2 cycles and the crack propagation a with respect to the number of cycles N is seen in figure 20.

Validation of Paris parameters
The FEA results from 'Direct Cyclic' utilising Paris parameters estimated by the proposed approach were compared with the experimental or reported data [57] on a / da dN versus ΔG plot, as depicted in figure 21, revealing a good agreement between each other.

7. Conclusion and future scope
In this work, the S−N curve from the reported data has been implemented into the development of fatigue crack growth rate or Paris curve through FEA based numerical approach by implementing the virtual crack closure technique (VCCT) for the four-point bending end notched flexural (4ENF) honeycomb sandwich composite specimen. A relationship has been successfully established between the S−N curve and a representative curve of fatigue crack growth rate, and it was successfully demonstrated that only one type of curve is needed as far as experimental testing is concerned and others can be estimated from the proposed approach. SERR results obtained from FEA, with proper implementation of (VCCT) provide a reliable source in the evaluation of fracture energy, which is fundamental in this approach.
The Paris parameters obtained from the approach were successfully utilised in Abaqus ® , implementing the virtual crack closure technique (VCCT) in a direct cyclic step and a good agreement between the experimental or reported data from literature and obtained FEA results was demonstrated, thereby justifying the implementation of numerical techniques in the evaluation of Paris parameters from S−N data. The analysis further gave an insight into the propagation of crack in a unit honeycomb cell, at the interface of facesheet and honeycomb core, with merging and branching of the cracks from adjacent cells.
As far as future recommendations are concerned, the shape or type of crack, the geometry of the structure or specimen, and the material of the structure etc. are worth investigating for better reliability and better applicability of honeycomb sandwich structures in various applications.