Numerical and Experimental Assessment of Post Impact Fatigue Life of Glass-fiber-reinforced Aluminum Laminates

In this research, dynamic progressive failure of Glass-Fiber-Reinforced aluminum laminates under low-energy impact was modelled. Intralaminar damage models, strain-based damage evolution laws, Puck failure criteria were used in ABAQUS-VUMAT software for modelling. Bilinear cohesive model was used for interface delamination, and the JohnsonCook models were employed for aluminum layers. Damage evolution behaviours of this hybrid composite were calculated. After that, energy dissipation mechanisms were examined to identify the progressive failure and delamination of composite layers and plastic deformation of aluminum layers. In order to determine stress intensity at crack tip, the analytical model for constant-amplitude fatigue crack propagation according to Paris law was applied. Also, bridging stress along crack length in aluminum layer was investigated by correlation between the delamination growth rate and energy release rate in hybrid composite layers. The obtained findings indicated that the highest amount of peak low velocity impact force belonged to Glare 4 3/2. The presented numerical method based on bridging stress phenomena can successfully be used for predicting the post impact fatigue life of Glare.


Introduction
A fiber metal laminate (FML) such as GLARE is a hybrid composite consisting of thin aluminum layer and fiberreinforced epoxy composite layer. The combination of the various metals having different strengths and the composite layer with different orientation can lead to the improvement of aircraft structure design. The most widely used metal for FML has been aluminum, and the fibers have been generally Kevlar or glass. FMLs with glass fiber (GLARE) and Kevlar fibers (ARALL) are used in aircraft structures [1][2][3][4][5] since they have several outstanding good points like weight saving, good impact resistance, and damage tolerance. Furthermore, the maintenance of metallic aircraft structure such as easy machining, forming, and mechanical fastening ability and damage inspection can apply on FMLs too [6]. GLARE laminate is a good candidate material for lots of applications in future aircraft structure, i.e. the fuselage, pressure bulkhead, tails and wings etc. These outstanding features of GLARE laminate have been already proved by the commercial aircraft Airbus A380 [7]. Impact damage is one of the concerning topics in aerospace structures. The fact of being unable to visually identify the interior damage existing in the composite layers, which sometimes go further the affected area, has been still a major safety problem. Accordingly, predicting the accurate amount of internal impact damage to Glare is critical. Due to out of-plane loads, like impacts, FMLs may tolerate damage due to various mechanisms including: (i) plastic deformation of the metal layers; (ii) matrix cracking and fiber failure; (iii) delamination between composite plies; and (iv) debonding of the metal and composite layer. In general, considering the absorbed energy, GLARE will act better in comparison with aluminum for impact loading. Also, GLARE at the time of being exposed to impact also develops a noticeable dent, which looks like monolithic alloys. The appearance of this dent demonstrates that an impact has occurred. Considerable internal damages including little surface appearance are Glass-fiber-reinforced Aluminum Laminates formed in traditional composites when affected by the same impact energy levels [8][9][10][11][12][13]. Due to the scatter in the results, impact tests usually waste time, and large sample sizes are needed for them. Hence, it has been tried to model the response of FMLs. Modelling damage from in and out of plane loads delamination is just related to the finite element analysis. There have been some analytical methods created to determine the impact loading [10,12]. In a study, Wu et al. [14] experimentally investigated the damage tolerance of fiber metal laminates which have been empowered bidirectionally. Caprino et al. [15][16][17] represented a mechanistic model to prognosticate the reaction of square fiber glass-aluminum laminates under low velocity impact, for which, several experimental tests should be conducted. Liu et al. [18] developed a theoretical and numerical techniques to simulate the dynamical impact responses of carbon fiber composite laminates. In this research, their model was developed for analyzing the low energy impact of GLARE, and the Johnson-Cook model was used for aluminum layers. It should be noticed that the energy dissipation mechanisms of GLARE laminates have been connected to the damage evolution and plastic deformation. This research well explained the damage evolution behaviors and energy dissipation mechanisms of GLARE laminate due to different failure modes including fiber breakage, matrix cracking, shear failure and delamination. A numerical technique was implemented by ABAQUS-VUMAT (user dynamic material subroutine). The impact forcetime/displacement curves and the energy dissipation, numerical results with different impact energy by FEM have been compared with the experimental results.
GLARE has been known for greatly tolerating fatigue and damage. Some results about fatigue behaviour of GLARE are as follows [19][20][21][22][23]: (1) Fatigue crack growth rates in GLARE panels are significantly lower than those in the monolithic aluminum panel under equal loading scenarios; (2) at the major part of fatigue life, growth rates of fatigue cracks in GLARE laminates have been reliably fixed; (3) at larger crack lengths in GLARE laminates crack growth rate was in creased.
Fatigue crack growth of fiber metal laminates (FML) has been studied by many researchers. Marissen [24] and Guo and Wu [25,26] developed an analytical prediction model for crack growth in FML. Both models used the superposition principle for stress intensity factor (SIF) which means that SIF at crack tip was established on crack opening in the aluminum layers and crack closing because of bridging stress. Marissen assumed that the bridging stress is constant along the crack length [24,27]. This assumption proved to be incorrect by Guo and Wu [25,26]. In their point of view, in addition to shear deformation proposed by Marissen, there existsa persistent delamination shape increase (ellipse or triangular), which has been identified by the crack growth in the aluminum layers and the delamination growth at the crack center. The inaccuracy of this assumption was shown by all data reported on fatigue crack growth in FMLs [24][25][26][27][28][29][30]. Alderliesten [30,31] developed an analytical method for crack growth on FML considering linear elastic fracture mechanics, and declared that SIF at crack tip is a function of the far field opening stress and the crack closing bridging stress in the aluminum layers. In this research a new prediction model was developed based on physical reality and accuracy for fatigue crack growth and life prediction. To improve the prediction models of this phenomenon, several facilitator assumptions were needed to be redefined to describe the crack growth behaviour in Simple mathematical relations. The two main facilitator assumptions for the method presented here were: (1) Linear Elastic Fracture Mechanics may be implemented, and (2) the empirical relation between the crack growth rate da/dN and the effective SIF range ∆K eff for thin monolithic aluminum sheets is also true for the aluminum layers in Glare. As a result, the effective SIF at the crack tip has been used to delineate the crack growth in the aluminum layers in Glare. This effective SIF needs to be incorporated when investigating the impact of crack opening restraint, fiber bridging, and delamination.
In this research, it has beentried to provide a model for prognosticating the very low velocity post-impact fatigue life of GLARE using Abaqus VUMAT and mathematical relations by using bridging stress phenomena and comparing with both low velocity impact and tension-tension fatigue test results.

Numerical Method for GLARE Low
Energy Impact

Johnson-cook Model for Aluminum Layer
For predicting the mechanical behaviour of aluminum layer, the Johnson-Cook relation without temperature consideration was used [32]: Where A, B and C are material parameters, ̅ is the equivalent plastic strain, and n is the material constant. A damage parameter Ψ was used for simulation of ductile damage of aluminum layer.
Where d 1 , d 2 , d 3 and d 4 are material parameters and p/q is the hydrostatic pressure per deviatoric stress.

Composite Damage Model
The intralaminar damage model was proposed by Kachanov [33] and then Lemaitre and Chabocheextended the model [34]. It has been capable to predict intralaminar failure of a composite material, and the results from the relations were presented by Donadon et al. [35] and Appruzzese and Falzon [36]. Damage variables represented the local damage in a Representative Volume Element (RVE) of the composite material. In one dimension, by identifyinga damage variable d, the effective stress , or the stress, can be computed over the section of the RVE. For an undamaged pristine material d= 0, whiled=1 denotes complete failure. Lemaitre and Chaboche expressed that the principle of strain equivalence can be applied to obtain the constitutive law of damaged material [34].
The relationship existing between the stress tensor and true stress tensor was: Where D is damage matrix: The stress tensor and strain tensor were related as: The relationship between the stress increment tensor Δ and strain tensor Δ wasas follows: Where is inelastic strain, and Δ are the damage and incremental damage matricesas expressed by Donadon et al [34].

Tensile Failure Modes
To predict tensile fiber and matrix failure modes, the same damage initiation and evolution law was used. was defined as strain based failure index to discover tensile damage initiation and it should be noticed that = 1 applies to fiber failure and = 2 for matrix failure [37]: Where is the failure initiation strain in tension and is the ultimate tensile strength in longitudinal direction. Figure 1 shows intralaminar damage model behaviour for tensile failure modes. To calculate the damage parameter for respective failure mode, the following equations were used: is the maximum strain, is intralaminar fracture energy and * characteristic length. Introducing the characteristic length guaranteed that a persistent energy release rate per unitarea of crack was generated independently from the mesh that the FE model used [22][23][24].

Fiber Compressive Failure
A similar approach was used for Fiber compressive failure [38]. The damage initiation and evolution law for fiber compressive failure has been shown in Figure 2.  Where is the failure index, is the failure initiation strain, isthe maximum strain at failure and is the damage parameter.

Matrix Compressive Failure
The first damage model under compressive loading for failure in the transverse direction was proposed by Hash in [39] and developed by Puck and Schurmann [40].  Figure 3 displays Mohr-Coulomb behaviour as suggested by Puck and Schurmann for matrix compression, and Figure  4 shows Rotation of stresses on fracture plane rotated relative to the material coordinate system.

Interlayer Damage
The interlayer damage for delamination simulation by using the cohesive zone model is implemented. This method was constructed by supposing that there exists another material between the various composite layers with its own constitutive law. Regarding the cohesive crack model, the various structural elements or layers are connected by this interface. This theory states that around the crack tip, there is a cohesive damage area which connects the tractions for displacement jumps at the interface where the crack may appear. The damage initiation and the interfacial strength of the material have been correlated together, and the damage evolution has been related to the strain energy release rates. The remaining traction power was omitted, when the total area under the traction-separation curve became as big as the critical fracture toughness of the interface material. The most outstanding good point of this theory was that both the crack initiation and propagation were included in the same model, and this decreased th elaboration of the FE model. (Figure 5)

Material Properties
For simulating finite element model with intralaminar and inter-laminar damage together, a wide range of lamina properties was required to describe the material behaviour, and (Tables 1 and 2) show all the material properties used in this research.

Finite Element Modeling
Numerical modelling was conducted by ABAQUS-VUMAT for the evaluation of the impact related damage of GLARE specimens. A C3D8R solid hexagonal element was used. The Johnson-Cook model was implemented for aluminum layers. The dynamic penalty method was used to simulate the contact between the GLARE and impactor. Due to the excessive distortion of elements, an enhanced stiffness relaxation was employed to control hour glassing phenomena. Intra-laminar damage models with Puck failure criteria were implemented. Bilinear cohesive model was used to predict the delamination, and thefinite-thicknesscohesive element was established for the delamination interface between two layers.

Fatigue Crack Growth in FMLs
Mechanism of fatigue crack growth in FMLs has been a complex phenomenon. To characterize this phenomenonand to develop an analytical method, several simplifying assumptions were required. Though, the suggested method according to the real crack growth phenomenon, should bestill accurate. The SIF at crack tip was expressed as: bending tip farfield bridging The SIF expression based on the linear elastic theory for monolithic metals, caused by the far field stress applied in the Al layers was:

Secondary Bending Effect
If the laminate layout is unsymmetrical, a secondary bending will happen and leads to a shift of the neutral line. The thickness of a FML laminate is thin. Therefore, a beam theory could be used for the stress variation in a laminate caused by the secondary bending. Because of the eccentricity of the neutral lines, shown in Figure 6, the neutral axis shifted and led to a bending stress. Tocalculate a bending stress, at first, the displacement of the neutral line was needed to be obtained for the further consideration of the stress level variation in a laminate. By using the neutralline model [41], the bending moment, -6 , wasexpressed as: .z Where . is the applied force and 7 is the thickness of laminate. By using the beam theory for a thin laminate, bending moment was expressed as: for plane bending. A metallic surface of laminate should be divided into cracked ( @ AB , = 1 ) and uncracked (@ C , = 2) parts, respectively. Then, by putting eq. 19 into eq. 20, the differential equation yielded to: The general solution of the above differential equation was: Thus, the bending stress in thickness direction of FML could be calculated as:  The fibers in FMLs are insensitive to fatigue. They transfer a significant part of the load over the crack, and restrain the crack opening, as shown (Figure 7). Due to this restraining, the crack opening in GLARE was smaller as compared to the monolithic metal. This mechanism resulted in smaller cracktip SIF for equal crack length, and applied load compared to the monolithic metal.  The crackopening of a fatigue crack in an infinite FMLs at the location x along the crack can be expressed as [25,30,43]: x , δ f (x) and δ pp (x) represented the crack opening displacement caused by the far field stress, bridgingstress in aluminum layer, the elongation of the fibers over the delaminated length and the deformation of the glass/epoxy prepreg layer; respectively. The crack opening displacement far from crack caused by the far field stress was: 8 D , D and / are theelastic modulus of aluminum, far field stress of aluminum layerand crack length; respectively.
The bridged crack opening displacementE FB ;G? caused by bridgingstress in aluminum layer was calculated using series of crackopening displacements caused by the loads acting at apoint with positionG . G is the point load coordinates in the far field direction at the crack tip that the fiber stress affected it. , Where H;G? is the delamination shape, / I is start crack size after impact and FB,D is bridged stress at aluminum layer.
Index 1 and 2 at the above equations refer to the interface parallel and perpendicular to the loading direction. For this model, the expression of the bridging stress distribution was calculated numerically by assigning the crack length in Jelement as shown in Figure 8; and FB ;G? was obtained from: The crack closing bridging SIF caused by the bridging stresses can be obtained as: Heretofore, the SIF at the crack tip for a crack with an arbitrary delamination shape was determined; the crack growth rate can be obtained with the empirical Paris relation [30] for aluminum layer: Where v impact and b allistic are impactor velocity and ballistic limit of GLARE types which investigated by Seyed Yaghobi et al. [44] By using eq. 41, crack growth life of GLARE can be calculated as: 1 cg a n a cg eff

Specimens
Five types of GLARE -3 2/1, 4 2/1, 4 3/2, 5 2/1 and 5 3/2 r were fabricated with the thickness of aluminum layer of 0.4mm, and the thickness of composite layer was 2.2mm. It should be noted that, similar to a traditional composite, different GLARE lay-ups were possible, and to classify these laminates, a coding system was essential for design and production. The code for an arbitrary laminate was: GLARE A B/C-t where A defines the grade of the laminate as shown in (Table 3), B indicates the number of aluminum alloy layers, C indicates the number of glass/epoxy prepreg layers, and t indicates the thickness of aluminum layers.

Low Energy Impact Testing
The impact test was conducted by a drop-weight impact tower with a velocity range of 1.80 to 2.40m/s (Figure 9). After the first impact, a pneumatic breaking system was activated, and prevented the impactor to hit the specimen again. The specimen was clamped between the top and bottom steel plates of the fixture with 70×150 mm 2 with a 50 mm diameter circular opening at the center of each plate. The impactor was a10 mm diameter steel hemispherical with a mass of 6.29kg. Figure 9. Impact test machine.

Post-impact Fatigue Testing
The tension-tension fatigue testing was conducted by a SANTAM-50 machine at 2 levels of 117 and 350 MPa for 11 J impact energy and 80 MPa and 120 MPa for 18 J impact energy with a frequency of 10Hz as shown in Figure 10.

Low Energy Impact
( Figures 11 to 15) show the impact force-time history curves for GLARE 5 2/1 and 4 3/2 in 11 J 18 J impact energy. Some numerical oscillation appeared near the peak force in the curves, which arised mainly because of numerical impact and dynamic reaction. In the present investigation, numerical filtering technique was used, and the output sampling time was 0.005 ms. At about 2 ms and 2.2 ms time for GLARE 5 2/1 and GLARE 4 3/2 specimens, the impactor reached the lowest end and then started to rebound until complete separation between the impactor and laminate at about 3.7 ms and 3.5ms time; respectively. There was a relatively good consistency found between numerical and experimental results regarding the impact force-time curves. For both numerical and experimental results between GLARE 4 3/2 and GLARE 5 2/1, the higher peak force and less contact time belonged to GLARE 4 3/2. This was due to the number of aluminum layers which have higher stiffness and impact resistance. Figure 16 shows the maximum central displacement and permanent central displacement corresponded to the lowest end of the impactor and the zeroimpact force at the end of impact; respectively. Since GLARE 4 3/2 absorbed more energy due to the number of its aluminum and composite layers compared to GLARE 5 2/1, its displacement was lower. Therefore, it can be concluded that in response to a very slow velocity impact, with increasing number of composite and aluminum layers in GLARE, more energy was absorbed, and the displacement decreased. At impact force-time/displacement graphs, there was a good consistency between the numerical and experimental results.     The energy dissipation of laminates has been shown in Figure 16. The energy dissipation has been the result of the plastic deformation of aluminum layers, the delamination between different layers, the damage evolution of composite layers and the frictional contact between different layers. When the impactor was going to contact the Aluminum layer, the energy dissipation began to be released; and because of absorbing the kinetic energy of the impactor, the laminate was going to be deformed. When the Impactor reached to its lowest impact point all its kinetic energy transformed into strain energy and dissipated energy. Afterwards, rebounding the impactor and elastic unloading was started. The total unrestrained energy of laminates is was estimated by the area under the impact forcedisplacement curve. The compression between FEA and experimental results for 11 J and 18 J impact energy showed that the errors for the total dissipated energy for GLARE 5 2/1 were 6.62% and 5.73%, and for GLARE 4 3/2 were 4.78% and 4.40%; respectively.    Energy dissipation at the end of 11 J and 18 J impact are showed in Figure 19. The dissipated energy due to interlaminar damage and delamination from the impact in GLARE 4 3/2 For 11 J impact energy 8.77% and for 18 J impact energy 9.96% is less than GLARE 5 2/1 which expressed that GLARE 4 3/2 have smaller damage. By increasing the impact energy, the amount of plastic dissipation energy increases at the aluminum layers.

Tension-Tension Fatigue
The stress level-cycles diagram (S-N) for GLARE 4 3/2 and GLARE 5 2/1 is shown in Figure 20 at 18000cycles. S-N curves are linear and fatigue life of GLARE 4 3/2 is longer than GLARE 5 2/1.

Figure 20. Stress level vs Cycle for GLARE 5 2/1 and GLARE 4 3/2 after 11 J impact.
Bridging stress distribution for GLARE 5 2/1 and GLARE 4 3/2 at 117 MPa and 350 MPa stress level has been shown in Figure 21 and Figure 22. The overall shape ofthe bridging stress diagram was constant for GLARE 5 2/1 and GLARE 4 3/2, and the bridging stress was increased near th ecrack tip. As has been known, this phenomenon had a direct effect on the post impact fatigue life of GLARE, and increased it much compared with the monolithic aluminum (almost 16 times).    In Figure 23 and Figure 24 Crack growth rates and tension-tension fatigue stress for GLARE 5 2/1 and GLARE 4 3/2 were shown after 11 J impact energy at 350 MPa. As clearly shown in these figures, crack growth rate was almost constant, while the crack growth rate in aluminum always increased, and this was another reason why GLAREhad an advantage over aluminum.  (Table 4) shows the comparison between the experimental and predicted post impact fatigue lives of GLARE. As shown in (Table 4), the predicted post impact fatigue lives of GLARE 5 2/1 and GLARE 4 3/2 were consistent with the experimental tests' results.

Conclusion
GLARE is a hybrid composite consisting of thin aluminum layer and fiber-reinforced epoxy composite layer. Generally, the impact tests waste a lot of time, and due to scattering in the findings, they need large sample sizes. Hence, it has been tried to model the response of FMLs. Modelling damage from in and out of plane loads delamination has been just related to the finite element analysis. The mechanisms of energy dissipation in GLARE laminates have been connected to the damage evolution and plastic deformation. This research well explained the damage evolution behaviors and energy dissipation mechanisms of GLARE laminate due to the different failure modes including fiber breakage, matrix cracking, shear failure and delamination. Numerical technique was implemented by ABAQUS-VUMAT (user dynamic material subroutine). The impact forcetime/displacement curves and the energy dissipation, numerical results with different impact energy by FEM were compared with the experimental results. In this research, a new prediction model was developed based on the physical reality and accuracy for fatigue crack growth and life prediction. To improve the prediction models of this phenomenon, several facilitator assumptions were redefined to describe the crack growth behaviour in simple mathematical relations. The following results have been obtained: 1. GLARE 4 3/2 showed the highest amount of peak low velocity impact force because it had three layers of aluminum.