Progressive Damage Modelling of Transverse Impact Behavior of Adhesively Bonded Composite Joints

Corresponding Author: Konstantinos Tserpes Department of Mechanical Engineering and Aeronautics, Laboratory of Technology and Strength of Materials, University of Patras, Greece Email: kitserpes@upatras.gr Abstract: A three-dimensional dynamic numerical model was developed, using the commercially available finite element software LS-DYNA to simulate the Low Velocity Impact (LVI) case on bonded joints consisting of woven CFRP adherents. The adhesive failure simulated using the Cohesive Zone Method (CZM), while the failures in the composite material simulated using a progressive damage material model based on Hashin’s failure criteria. The numerical results were compared with the corresponding numerical and experimental data found in literature, showing good agreement regarding the damage area and the global structural response. However, the maximum contact force was underestimated by 18.2% due to premature fiber failure in the lower adherent as a result of tensile bending stresses. In the adhesive layer, the stress distribution was found to alter from tensile near the edge of the lower substrate to compressive close to the edge of the upper substrate. The cohesive failure initiated from the edge under mixed mode loading while it propagated along the overlap length under mode II dominated loading. Following, a study on the effect of the post damage parameters of the CFRP material model to the damage accumulation of the joint was conducted. These parameters were calibrated according to experimental data from literature. The modified model was used for the estimation of residual tensile strength of the bonded joint containing impact induced damage. The results show tensile strength decrease of 16.3% for the joint with 25.4 mm overlap length and impact energy of 10 J. Finally, the effect of overlap length and impact energy was studied, in the means of damage accumulation in the CFRP and adhesive material and residual tensile strength of the joint. It was found that increasing the overlap length lead to reduction of the disbonding area due to impact loading resulting in lesser decrease of the residual tensile strength of the joint. For low impact energy of 5 J the associated damage in the adherents and the adhesive is minimal and doesn’t affect the tensile strength of the joint. Increasing the impact energy resulted in significant increase in the damage accumulation of the joint leading to equivalent decrease in the tensile strength.


Introduction
In recent years, the implementation of adhesives in the joining method of composite materials in structural applications in aeronautics has increased due to the many advantages over the traditional mechanical joints. Mechanical joints using bolts/fasteners have to deal with the potential problems such as unevenly distributed load in threads (Zhou et al., 2015a), stress concentration in threads (Zhou et al., 2015b), shear damage of bolts (Zhou et al., 2016), heat shorts for pure metal bolts or brittle fracture for pure ceramic bolts (Zhang et al., 2014), residual stress for novel metal-ceramic composite bolts (Zhou et al., 2015c;2015d). The use of adhesive bonding comes in conjunction with the establishment of the use of composite materials, as it enables to use them more efficiently. Some features which make adhesive bonding attractive include improved appearance, improved aerodynamic surfaces, good sealing, high strength, low weight, low stress concentration, low cost, fatigue and corrosion resistance and cost efficiency (Breuer, 2016;Banea and da Silva, 2009).
The joints are generally designed to carry in-plane loads, although they are prone to transverse impact loading from tool drops, flying debris or fragments. Vaidya et al. (2006) studied experimentally and numerically the mechanical behaviour of single lap bonded joints subjected to transverse Low Velocity Impact (LVI). They noticed that the transverse loading causes greater peel stress concentration in the adhesive layer compared to in-plane loading, due to the considerable bending deflection. They also noted that stress distribution around the crack tip changes from mixed mode to primary mode II as the crack propagates through the adhesive.
Impacts are categorized into low and high velocity. High velocity impact response is dominated by stress wave propagation through the material, in which the structure does not have time to respond, leading to very localized damage. As a result, the most critical failure modes include the formation of shear plug, fiber breakage, penetration and delamination. On the other hand, in LVI, the dynamic structural response of the target is of utmost importance as the contact duration is long enough for the entire structure to respond to the impact. In consequence more energy is absorbed elastically and the damage ranges in a larger area. Typical failure modes in the composite materials during LVI are matrix mode which occurs parallel to the fibers, delamination mode, fiber mode and penetration. Normally, during a LVI, matrix and delamination failures occur in lower impact energies compared to fiber failure and penetration. The impact energy required for each damage mode initiation is affected a lot by the geometry, material properties and boundary conditions.
The events of low velocity impacts are very common in airplanes during maintenance works and it is a matter of concern due to the low visual inspectability of the induced damages and due to their impact on the structural integrity of the joints. In literature, many experimental and numerical studies have been conducted to investigate the mechanical behaviour of composite materials under impact load. However, there have been very limited studies regarding the transverse impact on bonded lap joints in composite structures and often the association between the damages in composite material and adhesive material is not accounted. This study focuses on the damage of composite bonded joints under LVI loading and its associated loss of tensile strength capability. Finite Element (FE) models are implemented for this purpose in LS-DYNA software, to analyze the LVI and the uniaxial tensile loading of the joint. The results obtained by the FE model were compared and validated by experimental LVI data conducted on CFRP joints from de Oliveira et al. (2012), leading to the understanding of the experimental results that can be correlated to numerical simulation (Arrigoni, 2020).
The present paper is divided in four sections. The introduction is followed by the Finite Element modeling section where the developed FE model is presented as well as the adopted theories concerning the composite damage and debonding initiation and propagation. In section 3 the numerical results are presented. In the last section of the paper, the conclusions concerning the work are summarized.

FE Model
The geometry and the loading scenario that are selected for simulation are based on the available experimental data found on literature. In this work, a single-lap joint configuration was modelled with each adherent containing 8 plies of plain weave CFRP material with quasi-isotropic stacking sequence [(0/90)/(45/-45)]2S. The adhesive material was an epoxy film with 0.21 mm thickness and 25.4 mm overlap. Impact experiments were conducted in a drop tower apparatus with cut out dimensions of 150100 mm in accordance with ASTM -D7136 standard. The striker has a mass of 1.5746 kg with hemispherical head with 10 mm diameter made of high strength aluminium alloy. Impact energy is 10 J that corresponds to 3.56 m/s impact velocity. More details about the problem simulated in the current work can be found in de Oliveira et al. (2012). The testing apparatus can be seen in Fig. 1. The modelled geometry with the basic dimensions are presented in Fig. 2.
The numerical simulation conducted in the LS-DYNA FE software, which uses an explicit time integration solution algorithm based on central difference method. The aforementioned solution algorithm favors the solution of short duration dynamic problems like an impact event with high non linearities occurring during material failure.
The composite adherents were modelled using 8noded solid element with 3 DOF per node with reduced integration Element Formulation (ELFORM 1). The film adhesive was modelled using 8-noded cohesive elements (ELFORM 19). The impactor was modelled using 8-noded solid elements with rigid material model (MAT 20) to reduce computational cost without any noticeable differentiation when compared to an elastic material. In the numerical simulation the available hourglass control algorithms provided by LS-DYNA were tested by their effectiveness preventing the hourglass deformation modes in the underintegrated solid hexahedral elements. Viscous formulations of hourglass control algorithms were generally very ineffective and considered inappropriate for low velocity impacts. The stiffness form of Flanagan and Belytschko (Hallquist, 2006) algorithm (IHQ 4) had better results in preventing the modes. However, it resulted in unrealistic stiffness increase of the structure and high hourglass energy dissipation. The assumed strain co-rotational stiffness form of Belytschko -Bindeman algorithm (IHQ 6) was very effective counteracting the hourglass modes, with low hourglass energy and low stiffness increase. Thus, this algorithm was used in the numerical simulations.
A mass weighted damping was implemented (DAMPING_PART_MASS) on the nodes of the laminate simulating the material damping, to reduce the oscillations during the impact event. Using this type of damping, an external force is applied on each node with magnitude proportional to the nodal velocity and mass and with direction opposite to the nodal velocity. By trial and error, a damping coefficient with value Ds = 800 is selected.
In the case of transverse impact of the joint, the 4 lower edges of the joint were constrained from lateral (parallel to the edge) and vertical translation, allowing sliding (perpendicular to the edge) and rotation. For the in-plane uniaxial tension case, on one end of the adherent, all translational degrees of freedom of the upper and lower nodes of the adherent were constrained simulating the clamping grips, while the other end was constrained from lateral and vertical translation. Load application was performed in the longitudinal direction on the same nodes via displacement control until the joint failure. Contact was modelled between the impactor and the adherent using the ERODING_SURFACE TO_SURFACE keyword with SOFT = 2 option which refers to a pinball segment based contact algorithm for the calculation of the contact stiffness.
A non-uniform mesh is used in the numerical calculation, as shown in Fig. 3. The joint contains two mesh refinement regions in the overlap area. The first is away from the point of impact, where the refinement is done only across the overlap length to consider the high stress gradient that occurs in that direction. The second is near the impact area, where the refinement is done in both directions to capture the damage progression of the composite material. On each ply, one element is used in through thickness direction for computational efficiency. The study of the mesh size variation showed significant influence on the dynamic structural response and the failure load of the joint. The strong mesh dependency of the progressive damage modelling method has been reported in other works in literature and referred as significant disadvantage of this method. This is explained by the softening post damage mechanical behavior of the material, which leads to strain localization during material failure resulting in energy dissipation, which depends by the elements' size (Lapczyk and Hurtado, 2007). A mesh convergence study was conducted to find an appropriate element Composite adherents Adhesive Impactor size recording the variation of the composite material peak failure load. It was noticed that the required mesh size to achieve convergence depends on the material properties. The baseline model with brittle post damage material properties showed sufficient convergence with 64.351 elements (Fig. 4). On the other side the modified model with more ductile post damage material properties needed further element size decrease with 112.831 elements and thus this mesh size was used for the following simulations.

Progressive Damage Modeling of Composite Material
Numerical modelling of the composite adherents is done using LS-DYNA material model MAT 162. This material model incorporates strain-based Hashin-type failure criteria (Hashin, 1980) and damage mechanics approach of Matzenmiller et al. (1995) for the characterization of the strain softening behaviour after the damage initiation. It also accounts for strain rates effects in strength and elastic modulus which is useful for simulations of high velocity impacts.
ΜΑΤ 162 requires a total of 34 material properties and parameters to describe the full response of an orthotropic unidirectional or woven composite material under different damage modes. There are nine elastic constants and ten failure strength properties which can be determined from standard ASTM test methods, except for fiber crush Strength (SFC) and fiber shear Strength (SFS). The rest of material properties are determined by non-standard experimental techniques such low velocity impact experiments and quasi-static punch shear testing in conjunction with supporting numerical simulations. More details regarding the methologies for determination of the material parameters can be found in (Gama, 2015;Xiao et al., 2007;Jordan et al., 2014). A main drawback of this material model is the excessive material testing required to determine all the input parameters. In this work, the post damage parameters of the material model are calibrated by fitting to the experimental structural response and their influence is examined. The values of the required input parameters are presented in Table 1.
More theoretical details and constitutive equations regarding the material model are given below.
A set of quadratic failure functions (Table 2) is used to define the initiation of different damage mechanisms related to fiber fracture, fiber crush, fiber shear, in-plane matrix crack and delamination. Table 3 and 4, the description of the symbols of Equation (1) to Equation (7) is presented.
Fiber compression -direction a 2 22 Fiber compression -direction b 2 22 10 10 10  Fiber mode failure in a direction f8, f10 Fiber mode failure in b direction f11 Fiber crush failure mode f12 Perpendicular matrix failure mode f13 Delamination mode r7, r9 Damage threshold of fiber mode failure in a direction f8, f10 Damage threshold of fiber mode failure in b direction f11 Damage threshold of fiber crush failure mode f12 Damage threshold of perpendicular matrix failure mode f13 Damage threshold of delamination mode Quasi-static shear strengths SSRC Shear strength when c <0

Sab
Layer shear strength due to matrix shear failure After damage initiation, the progressive damage model assumes linear elastic response which is governed by the reduced stiffness matrix with the updated damage variables ϖi. Damage variable ϖi with i = 1,…,6 is used to relate the onset and growth of damage to stiffness reduction of the material. The compliance matrix [S] is related to the damage variables as: A damage variable ωi for an individual failure mode j is given by: The damage variables are coupled to the individual damage modes (j = 1,…,6) by the vector-value function qij of the form:

Cohesive Zone Modeling of Adhesive Material
Numerical modelling of the epoxy film adhesive is done using LS-DYNA material model MAT 138. This material model simulates the adhesive material with the Cohesive Zone Method (CZM). It includes a bilinear traction-separation law with a quadratic mixed mode delamination criterion. Theoretical details and constitutive laws regarding the CZM, can be found in (LSTC, 2014;Floros et al., 2015).

Results and Discussion
In the following subsections, the results of the FE models of LVI are presented. For the model in section 3.1, the elastic and strength properties used for the MAT 162 material model of CFRP found on literature from the corresponding author (de Oliveira et al., 2012;Mendes and Donadon, 2014), while for the rest parameters, baseline parameters are used, which was found in Deka et al. (2008). In section 3.2, the critical MAT 162 parameters are calibrated according to the experimental data of the LVI testings. Finally, the influence of the overlap length of the joint and of the impact energy are studied by means of damage accumulation in the joint materials and residual tensile strength of the joint. due to bending tensile stress. This localized failure instantly propagates through the upper plies inducing further fiber fracture and delamination of the adhesive layer, leading to a sudden decrease of the reaction force. The numerical model also predicts 5.22 J energy absorption, which is in good agreement (5.7% deviation) with the mean value of experimental data.

Results of the FE Model with Baseline CFRP Material Properties
The damaged area in composite and adhesive material predicted by this study has similar shape but it is larger, compared to the numerical results of de Oliveira et al. (2012). This correlates with the premature fiber failure, which results in increase in size of damage to dissipate the remaining impact energy. In Fig. 7 the damage areas of the numerical models are shown in comparison.

Influence of Properties of CFRP Material on Damage Accumulation of the Lap Joint
For the LVI case, the critical material parameters are determined according to the observed failures and stress field in the FE model and according to the literature references about the typical damages in similar problems. In the next subsections, a study is conducted regarding the influence of these critical material parameters in the structural response and the damage accumulation of the bonded joint.
The critical parameters studied are about the damage softening property of the composite material for fiber failure (m1 and m2) and matrix failure (m4). Material properties related to strain rate and laminate penetration are estimated to have low influence since the impact velocity is very low and thus not investigated.

Strain Softening Parameter for Damage in Fiber (m1 and m2)
The strain softening parameter for damage in fiber has the biggest influence in this FE model, since this was the main damage, which lead to macroscopic failure of the bonded joint. The tested values are lower related to the baseline model ranging from m = 0.3 to m = 0.9. Lowering the values of parameter m, the total strain energy for failure is increasing making the material failure of CFRP more ductile. The woven CFRP material is balanced in the two principal material directions, so m1 = m2 = m is assumed. Figure 8 the force -time curves for the four parameter's values are plotted. The curves are identical in the initial loading region until failure. As expected, lowering the m value, the failure load is increasing. Similarly, delamination area in the adhesive is lowered, especially for the case of m = 0.3, which further indicated the association between the events of fiber rapture in the low laminate and the adhesive delamination. The absorbed energy for the case of m = 0.3 is 2.43 J, which is much lower from the other cases where it is over 5 J. That happens due to the decrease in damage accumulation of composite material, as shown in Fig. 9 where the area subjected to fiber damage and total fiber failure is represented.

Strain Softening Parameter for Damage in Matrix (m4)
Figure 10 the force -time curves for the parameter's values which simulate hardening, perfect plastic, ductile softening and brittle softening mechanical behavior of the matrix are plotted. In this graph, it is observed that the models with m = 0.1 (ductile softening), m = 0.001 (plastic) and m = -0.2 (hardening) have identical response and damage accumulation. On the contrary, on the models with more brittle matrix with m = 0.5 and m = 1 there is a big increase in delamination area, since the composite material has reduced resistance in delamination damage growth. In the cases with relatively brittle matrix material, the CFRP also suffered from fiber mode shear damage due to the high transverse shear stresses on the region adjacent to the impact point as shown in Fig. 11. The result of the increased damage accumulation is reflected on the force -time curve with decrease in stiffness on the region with force over 2.5 kN and with decrease in the peak force. However, this is not reflected in the absorbed energy of the joint or in the delamination area of the adhesive.

Damage Limit Parameter for Elastic Modulus Reduction (ωmax)
This material model parameter defines the upper limit of the damage variable ϖi and thus the minimum value of the reduced elastic modulus according to Equation. 12. The parameter ωmax affects the mechanical behavior of the material model subjected to high strain.
In the force -time graph (Fig. 12), the influence of ωmax appears after a main failure occurs in the structure, which results in sudden loss of contact force. At this stage, the damage propagates through the material and the parameter ωmax acts as a barrier in this occurrence, decreasing the drop in the contact force. The decreased damage growth with smaller values of the parameters is also reflected in the absorbed energy of the joint and the delamination area of the adhesive, which are also lower.

Results of the FE Model with Modified CFRP Material Properties
The resulted parameters for the CFRP material that are determined from the calibration procedure are m1 = m2 = 0.3, m4 = 0.5 and ωmax = 0.993. The modified model had very good correlation with the experimental results as shown in the graph of Fig. 13. The peak force and the stiffness of the modified model is near identical with the experimental data while the post damage response also shows good agreement.
Following the damage progression in the joint during the impact is briefly presented as individual events Fig. 13 When major failure propagation occurs, it is also reflected on the force -time curve with steep decrease of the contact force as in E, F and G. Figure 14 illustrates the progressive adhesive failure on different instants of impact. On moment D, the crack initiates under mixed mode loading with mixity ratio β = 1.2. Due to the significantly lower Gc in mode I loading, the out of plane traction governs the adhesive failure. On moment E, adhesive debonding initiates on the region adjacent to impact point under mode II (β = 2.4) governed loading, as the projectile tends close the adhesive crack due to compressive forces. On moment F, the debonding propagates from the overlap edge towards the center of the overlap under mode II (β = 3.1) loading. The plies with 0/90 orientation undergo more fiber damage as they are oriented in the bending direction. Concerning the delamination and the fiber mode shear failure, they are developed equally on the internal plies of the top laminate (2nd to 7th plies). In these plies, the delamination failure forms a circular disk with 9 mm diameter and the fiber shear failure forms a circular ring with 7 mm external diameter.

Influence of Overlap Length on Damage Accumulation of the Lap Joint
As shown in the force -time curves in Fig. 15, increasing the overlap length leads to an increase of the stiffness of the joint. This results in a decrease of the total deflection and in an increase of peak force. Moreover, for the case of overlap length of 25.4 mm, delamination area is 312.7 mm 2 , while for length of 35 mm the delamination area decreases rapidly to the value of 69.3 mm 2 , which is located adjacent to the impact location and is caused due to local fiber failure. For the case of 15 mm overlap length a small increase of the delamination area (343 mm 2 ) is noticed. Figure 16 the force -displacement curves of the uniaxial tensile loading of the damaged and intact joint configurations are plotted. The joint with 35 mm length has an identical failure load of around 72 kN for the cases with and without damage, which is caused by fiber failure of the composite adherent near the edge of the overlap. The 25.4 mm joint has a residual tensile strength reduction of 16.3% with failure mode caused by delamination of the adhesive material. Finally, the 15 mm joint has reduction of 26.4% of the tensile strength. It is observed that the damaged joint with 15 mm overlap length doesn't exhibit sudden failure as the other cases as seen in Fig. 16. On the contrary, the delamination failure occurs in two stages since the propagation of the debonding is blocked temporary until the applied displacement increases, causing total joint failure.
In the uniaxial tension, the adhesive film is loaded by mode II dominated loading with mixity ratio of β ≈ 5 before the formation of a plastic zone on the joint overlap. Before the delamination of the joint, the adhesive film is loaded with mixed mode loading with mixity ratio of β ≈ 2.

Influence of Impact Energy on Damage Accumulation of the Lap Joint
Increasing the impact energy results in an increase of the peak force during impact, as shown in Fig. 17. The total deflection of the joint is increased but the impact duration is not affected, as it is influenced solely by the joint stiffness and projectile mass. On the 10 J and 15 J impact cases, a dramatic increase in damage accumulation of the joint is observed compared to the 5 J case (Fig. 18). This mechanical behavior indicates the presence of a threshold energy, which after reaching it, the damage propagation in the composite and adhesive material is greatly increased. Figure 19 the force -displacement curves of the uniaxial tensile loading of the damaged and intact joint configurations are plotted. The failure load for the joint with 5 J impact energy is identical with the intact joint. This is explained by the minimal damage formation of the joint during the trasnverse impact loading. On the cases with 10 J and 15 J impact energy, there is a gradual decrease on the residual tensile strength of the joint by 16.3% and 23.4% respectively. From the aformentioned numerical results, in case of high velocity impact, the phenomena would be more localized close to the impact area leading mainly to fiber breakage, penetration and delamination

Conclusion
Comparing the finite element model of this work with the literature's experimental data, they are in good agreement concerning the stiffness of the joint and the energy absorption (5.7% deviation) during the impact. The peak force during impact predicted by the model with the baseline material parameters was underestimated by 18.2% due to fiber failure on the bottom plies of the lower adherent. After the calibration of the CFRP material parameters, very good correlation was achieved regarding the mechanical response of the structure.
The calibration procedure resulted in a composite material with ductile post damage softening mechanical behavior. For a more reliable adjustment of the material parameters, more data are needed like experimental results with different impact energies and data from non destructive testings.
On the mesh sensitivity study, a strong dependence between the mesh size required for convergence and the post damage softening parameters of the composite material was observed. Specifically, for material models with more ductile softening mechanical behavior, the mesh size needs refinement to achieve convergence.
In general, the weak link of the joint studied on this work is the CFRP material due to fiber fracture, which leaded to partial debonding of the adhesive film. The damage of the adhesive film in the transverse impact case initiates from the overlap edge under mixed mode loading and propagates through the overlap length under mode II governed loading. The model predicted that the damage induced by impact loading with 10 J energy resulted in decrease of the tensile strength of the joint by 16.3%.
On the models with varying overlap length, the joint with increased overlap length had greater stiffness, while the impact induced damage is decreased. On the uniaxial tensile loading model, the joint with 35 mm overlap length failed due to adherent material tensile fiber failure, while the joints with 25.4 mm and 15 mm overlap length failed due to debonding of the adhesive film. The model with 5 J impact energy had significant decrease of damage accumulation in the joint compared to the models with 10 J and 15 J, resulting in no reduction of the tensile strength.