PROGRESSIVE FAILURE ANALYSIS OF LOW ENERGY IMPACT IN CARBON FIBER FILAMENT WINDING CYLINDERS

Composite material is very attractive for structural applications due to its inherent mechanical properties and low weight. The improvement in manufacture process allow composite materials be used even as primary structures in modern aircraft design such as Boeing 787 without loss of airworthiness. During service life composite structures can be damaged by collisions, dropping tools during assembly or maintenance, etc. Several impact studies were conducted for plates, but few regards curved geometries or cylinders. This study presents a progressive damage analysis of low energy impact on carbon fiber filament winding cylinders. Three different layups were used for experimental tests. A new material model based in continuum damage mechanics were implemented as a FORTRAN subroutine linked to finite element software ABAQUS for explicit dynamic analysis (VUMAT). Good correlations for force vs. time and displacement vs. time between the numerical model experimental test results were obtained.


NOTATION:
The application of composite materials in aeronautical industry have been increasing in the last decades, even in large civil transportation aircrafts, mostly due to composite high stiffness and low weight [22].Its intrinsic anisotropy allows achieve an optimal material performance, regarding the structure geometry and function.Composite materials can provide design lighter structures without loss of airworthiness, which is a very attractive characteristic for aeronautical industry.Thus, the application of composite material in aeronautical structures makes possible save more weight decreasing the fuel consumption or increasing the payload.
However, despite few new designs as Airbus 380 and Boeing 787, the application of composite in structures is still limited by the difficulty in predicting their service life [21].From the designer standpoint, the development of consistent failure criteria, which predicts with a reasonable accuracy the damage process, is essential to help the design process.
Although the great number of failure criteria and progressive failure models, the failure process and subsequent damage evolution are still underdeveloped.Also, the difficulty in predicting the structural failure modes requires a well-planned test program [14].However, due to the high specific stiffness and strength, advanced composites have been used as engineering materials for several years.On the other hand, not only the static behavior of composites is quite difficult to predict with consistency, but also the dynamic behavior caused by impact loads [25].This problem is so complex that some standards were developed in order to help the comprehension of the phenomena involved during impact tests.For example, ASTM D7139 (2007) [3] guides the impact test for composite flat coupons.
Several studies were conducted for composite flat plates under impact loads by several researches ( [6], [9], [23]), but only few studies were made regarding impact in curved geometry ( [4], [10]).Maybe, this can explain why the guidelines for design composite vessels was established more than 40 years ago and have recommended to apply high safety factors in order to avoid failure, mostly for pressurized vessels [10].
Despite the high strength in fiber direction, out-of-plane stresses caused by impact loads due to bird strike or dropping tool in a composite structure could lead to severe damage.In a metallic structure, this kind of damage is easier to detect, on the other hand for composite structures, this is not true [4].
Regarding the challenges shown earlier, the aim of the present work consists on performing a Progressive Failure Analysis (PFA), using a new material model based on Continuous Damage Mechanic [18].The PFA was carried out for carbon fiber filament winding cylinders under low energy impact load.The material model was implemented as a FORTRAN subroutine (VUMAT -user material), which was linked to the commercial finite element program ABAQUS TM .Explicit algorithm for time integration was used to perform the numerical simulations.Finally, the computational results were compared to low energy impact experimental tests in order to investigate the potential and limitation of the material model implemented for simulating the failure mechanisms of the composite cylinder under impact.

Test Coupons
There are several testing standards (tensile, compression, shear, bending, fatigue, impact tests) for composite flat coupons, but there is no standard for impact in composite cylinders.In order to overcome this limitation, in this work, the dimensions of the test machine as well as the available filament winding mandrel diameter were taken into account.Two dimensions of coupon are shown in Figure 1 Also, the impact tests performed aims to provide enough data to evaluate the accuracy of the material model implemented in order to simulate the behavior of carbon fiber filament winding cylinders under impact loads.Thus, three different lay-ups were manufactured (Table 1): type A with 3,49 mm thick; type B with 3,25 mm thick and type C with 3,54 mm thick (average).Those lay ups allow to assess how different levels of anisotropy affect the structure response under impact loads.Besides, according to the coupons manufacturer, the plies thickness was function of the orientation.Hence, the plies with fiber orientation equal 30 o are 31% thicker than 90 o plies, and 60 o plies are 15% thicker than 90 o .One important remark is that the material properties are classified data, thus the properties are not present in this study.

Drop Tower Apparatus
The impact test equipment provides the data of force and displacement when the impactor interacts to the coupon during an impact event (Figure 2(a)).It is also measured the strains in two different points of the cylinders using a bidirectional strain gages (Figure 2(b)).The impact tests were conducted in laboratories at Katholieke University of Leuven (Belgium).
The principle of impact test is very simple, because this test consists on a certain mass dropped from a certain height, and this mass hits the test coupon.So, the force and displacement are measured when the impactor and the coupon are in contact.
The impact apparatus (Figure 2(a)) consist on two guiding bars to drive the falling weight during the test.Regarding the used equipment, these guiding bars limits the impact height to 1.8 m.The test coupons are placed at the bottom of the machine (Figure 2(c)).
Several types of impactor head can be fixed onto the impactor frame allowing test different types of materials.For harder material, a more pointed impactor head is used, for the case for softer material a more blunted impactor head is used.The tests were performed using a round impactor head with 16 mm of diameter.The round head avoid, in certain level, penetration of the coupon caused by the use of sharp head.(c) Test schema using V-block in the base (30 J tests) A piezoeletric crystal, placed between the impactor head and impactor frame, is used as load cell in order to acquire the impact force.The displacement data was acquire using a light detector placed at the bottom of the apparatus, which measures the intensity of a Light Emitting Diode (LED) mounted on the impactor frame.Once the LED has a constant intensity, the distance measured is proportional to 2 1 ~d , where d is the distance between the LED and the light detector.
The displacement measurement system is simple, reliable and stable, once the light detector integrates the light intensity.Despite that, the measurements should be done in the most sensitive area of the calibrated displacement range.
At the bottom of the drop tower is the local to place the coupons to be tested.In this work, the cylindrical coupons are positioned in a flat surface for the 10 J impact tests ( Figure 2(c)) and in a "V" base for the 30 J impact test.

EXPERIMENTAL TESTS RESULTS
The experimental test results for 30 J impact energy, there is no data from the strain gages.Regarding experimental tests for 10 J impact energy, it can be seen the results for strain gages.However using C-Scan method, it was not detected any damage in the structure for all cylinders type impacted under 10 J (Figure 4(b)).
In general, the damage initiation is detected in impact test force vs. time history with a sudden force drop due to stiffness reduction caused by unstable damage growth [19].Also, observing the force vs. time graphics, it is possible to identify the delamination threshold by the first force sudden drop.During the damage process, matrix cracking is the first type of damage, which happens in a structure due to impact loading.As commented by the literature, this type of damage does not affect the laminate stiffness during impact [19].
For type A cylinders under 30 J impact energy, the force vs. time and displacement vs. time results are shown in Figure 3.The graphics show that the results for each coupon are nearly identical.Figure 3 shows that the force increases rapidly close to 0.0052 s and then a sudden force drop occurs.After that, the force increases again and a new sudden drop occurs.This trend repeats for close to 0.0087 s of the impact event.A similar behavior was presented by Minak et al. [16].This part of peaks and valleys could indicate the initialization of delaminations in several layers.After this first time, the unstable delamination propagation can cause the further oscillations in the force vs. time history, as observed by Schoeppner and Abrate [19].Also, Figure 3 shows that the maximum force level does not occur in the first peak and the maximum force value does not occur in the same time of the maximum displacement.
On the other hand, as mention before, the cylinders impacted under 10 J did not show any damage (Figure 4(a)) as observed by the C-Scan images.However, the force vs. time history shows a similar trend as shown by specimens impacted under 30 J (Figure 4

Force
This trend indicates that the peaks and valleys are not only caused by delaminations and the further unstable propagation also.The boundary conditions and the cylindrical geometry can affect very strongly the cylinder impact behavior.Furthermore, it is well known that the force vs. time history has many oscillations, which could be introduced by two sources: the impactor can load the structure in a natural frequency (impactor ringing) and the coupon can show flexural vibrations [3].Despite delaminations, this oscillatory behavior also could be explained by the wave propagation across the cylindrical structure and the cylinder natural frequency.Further finite element simulations showed that the delay in the basis reaction force possess an interesting correlation with the impactor force and the support reaction force (Figure 5). Figure 5 shows a finite element simulation results for force vs. time history for a cylinder under 10 J impact energy.In fact, the oscillatory behavior is not due to any kind of damage.When the impactor just hit the cylinder, the force increases, but there is no reaction in the support yet.After 0.0006 s, the reaction force in the support increases, but the force in the impactor decreases.The next peak in the impactor force corresponds to decrease of the support reaction force.This trend repeats until 0.003 s of the impact event, after this time, there are not clear correlations.It could indicate that the existence of an interference between the support reaction and the impactor leading to harmonic resonance in the force vs. time history.
It is possible to observe that all type A cylinders present matrix cracking and delaminations between several layers, also a small indentation mark (dent) is also observed.Despite the several damage detected, fiber breakages for all coupons are not observed.
For type B cylinders, the force vs. time and displacement vs. time results are shown in Figure 6.The graphics show that the results for each coupon do not possess a considerable dispersion.The same type A observations still valid for type B coupons, but the force peak intensity is higher for this cylinder type, due to the differences in thickness and lay-up.Delaminations, matrix damage and indentation marks are detected for type B cylinders, also.
Figure 7 shows the force vs. time and displacement vs. time results for type C cylinders.In this case, the dispersion in the displacement is more pronounced than for previous results.For example, the coupon 4 displacement results diverge even before the impactor reaches the maximum displacement.Type C cylinders present a higher load peak than other cylinders type, and Type B cylinders present a higher displacement.All the comments for Type A cylinders are also valid for Type C cylinders.Delaminations, matrix damage and indentation marks are detected for type C cylinders, also. Disp.

Force
In order to proceed with the investigation of how the damage evolves in composite cylinders under impact loads, the present work uses a progressive failure model proposed by Ribeiro et al. [18].This material model was implemented by Fortran as an user material subroutine (VUMAT), which was linked to a finite element dynamic explicit code (ABAQUS TM /Explicit).It is important to mention that the model regards plane stress state and fiber failures do not affect matrix failures and vise e versa.However, ply orientation affects the damage evolution.

Fiber behavior model
Under tensile load in fiber direction, the unidirectional composite laminae behavior is linear elastic with brittle fracture ( [2] [12]), which is detected using maximum stress criterion (eq.( 1)): After the failure detection, the damage variable in fiber direction, d 1 , is considered equal 1.On the other hand, the behavior of fiber under compression loads will be linear elastic until a specified value, after that, a non-linear elastic behavior should be considered.The linear elastic to non-linear elastic limit, 0 C X , is identified similarly to tensile failure as shown by eq.( 2).The non-linear elastic behavior is simulated using a secant modulus as shown by eq.( 3).

Matrix behavior model
The damage process in matrix is mainly caused by the stress components 22  and 12  .A non-linear behavior could be observed in some experimental tests in composite materials, mostly when the fibers and load are aligned.This non-linear behavior is due to inelastic strains and damage in matrix [17].To model the damage process in matrix, two internal damage variables, d 2 and d 6 , were used.Those variables range from 0, for undamaged material, to 1 for total damaged material.Based on Continuous Damage Mechanics (CDM), the hypothesis of effective stress links the damage variables to the stresses [7] as show by eq. ( 4).
According to CDM, microcracks and microvoids open if the matrix under tensile stress, then the damage in the material increases.On the other hand, under compression, the microcracks and microvoids are closed [7].To model this behavior, the damage variable d 2 only evolves when 22 0   .However, the damage parameter d 6 can evolve despite the sign of 12  .
Other important characteristic in the material model consists on simulating how the damage process evolves in matrix with the mutual influence of 22 ).The material model regards that under compression, the matrix possess a non-linear elastic behavior.Again, the secant modulus (eq.( 5)) was used to simulate the compression of the matrix: The initiation of damage in composite structures can be identified as deterioration of materials properties [13].This deterioration can be evaluated by performing cyclic tensile or compression tests.The model regards that the damage process starts when the stress vs. strain curve are no longer linear.
Based on CDM and using some adjustments for the Poisson's ratios in order to take account the damage effect, Matzenmiller et al. [15] proposed the compliance tensor shown by eq (6), which was adopted in the present work. Where In order to avoid the material self-healing, the damage parameters d1, d2 and d6 are the maximum value along the simulation.All the experimental tests for material characterization, the parameters identification as well as the damage variables evolution equation are described in details by Ribeiro et al. [18].

FINITE ELEMENT MODEL
The impact on composite cylinders were simulated using the finite element code ABAQUS TM /explicit.The finite element model is shown by Figure 8.The Figure 9(a) shows the boundary condition ( 0 Z U  ) applied in the place that the coupon touches the V-block (Figure 2(c)).The Figure 9(a) also presents the boundary conditions (all rotations -R x , R y , R z -and displacements U x and U y are restricted) and initial conditions for XX J impact energy applied on the impactor.Figure 8 (b) shows that the impactor is modeled using rigid triangular elements (R3D3).The mass of 3.24 Kg is applied in the "mass point".Also, the radius, the "reference point" (where the impact forces are obtained) and the mesh impactor.A four node reduced integration shell element with six degree of freedom per node (S4R) were used to model the cylinder.The reduced integration element was chosen in order to reduce the simulation time and to avoid numerical problems.Nevertheless, when performing a Progressive Failure Analysis (PFA), some elements could become excessively distorted, aborting the analysis.In order to overcome this situation, the mesh is refined until the severe element distortion was eliminated.Due to element distortion, different meshes were used to simulate the different cylinders lay-up.The mesh for Type A cylinder has 11168 elements, because it is not observed an excessive element distortion.However, for Type B and C cylinders, due to the lay-up and the geometry, it was necessary to use 19026 elements.
The contact between the impactor and the cylinder were modeled as hard contact for normal behavior and penalty for tangential behavior.It is important to mention that hard contact and exponential contact formulation (soft contact) can provide different results for flat panels as shown by [20].However, preliminary simulations showed that this did not occur for the cylinders investigated in this work.
Other important issue for simulating the impact on composite filament winding cylinders is to predict the dissipation of impact energy.Part of the impact energy is dissipated by irreversible process as damage (matrix and fiber damage and delaminations), but other part is dissipated by damping phenomenon.
Damping in composite materials is dependent of several factors, for example: fiber volume fraction, composite lay-up, environmental factors, force magnitude, etc. [24].Also, the structure geometry has an important influence in the impact response.
ABAQUS TM provides the Rayleigth's model for direct integration dynamic analysis in order to simulate energy dissipation mechanisms through damping [5].In fact, in finite element analysis, damping is treated as a matrix, which can be treated in two different ways, as a material property or as a numerical object to oppose the excitation forces [11].
Considering the equilibrium equation in dynamic analysis shown by eq. ( 7).
According to [24], the following transformation can be applied in eq.( 7), , and the resulting equation is multiplied by   After that, the eq.( 7) can be written as:: Eq. ( 8) is used to obtain the critical damping value, which is calculated by: Base on classical modal analysis       of eq.( 7), it can be shown that . Regarding the modal analysis of eq. ( 8), the modal matrix of eigenvalues vector is Thus, the modal fraction of critical damping is: Using the Principle of the orthogonality for     , eq.( 10) results in: The Rayleigth's model introduces damping in the structure as a linear combination of mass and stiffness system matrices [11] (eq. ( 12)), where the parameters  and  can be ob- tained using eq.( 11) and experimental data.
Therefore the damping is proportional to the mass and the stiffness.The mass contribution is related to the low frequencies vibrations and the stiffness contribution is related to the high frequencies vibrations [5].
It can be proved that for a mode i that the fraction of critical damping proportional to   M is: And the fraction of critical damping proportional to [K] is: In this work, a reversal analysis was performed in order to obtain the Rayleigth's parameters.However, in the future works, these parameters will be obtained by dynamic experimental analyses.

RESULTS AND DISCUSSIONS
Considering the experimental results for 10 J impact energy, only the cylinders impacted under 30 J were simulated.Thus, several computational simulations were performed for each cylinder lay-up in order to find the "best model configuration" for 30 J impact energy.The mesh refinement, element type and Rayleigth's parameters values were evaluated.
During the simulation, the transverse shear stiffness is constant and not vary as the damage evolves, this is a software limitation because in fact, the transverse shear stiffness change as the damage evolves, also for composites the 5  6 correction factor used in ABAQUS TM is not adequate for laminates.First, it was carried out the simulation results for type A cylinder without damping effects.The comparison between numerical and experimental results is shown by Figure 10.It is possible to observe that the simulation was able to predict the behavior of the impact tests capturing peaks and valleys, mainly up to 0.0087 ms.On the other hand, in the final part of simulation (around 0.0127 ms), a sudden loss of contact between the impactor and the cylinder causes a sudden reduction of the impactor force.The energy absorbed by the damage mechanisms in computational simulations without damping effects was not enough to dissipate the kinetic energy of the impactor.This same behavior was also shown in Figure 4 for low energy level (10 J) and no damage.In fact, the finite element model without damping of the cylinder works as just like a dynamic system with a mass and a spring, so during the rebound, the cylinder increases the impactor speed, leading to a sudden contact loss between the impactor and the coupon.It does not happen in the real structure, where the composite cylinders possess dissipative energy mechanism as damping and damage.Schoeppner and Abrate [19] commented that delaminations and fiber breakage are the type of damage, which dissipates more enegy, also indentation could absorb a considerable amount of energy [1].Thus, it is necessary to include the damping effects in the fnite element model.
After some simulations, the lower frequency Rayleigth's parameter ( ), which performs suitable simulations of the the impact in cylinders, was rather high (Table 2).The  parameter values for all cylinders are high due to the cylindrical geometry and the polymer matrix damping characteristics [8].The stiffness proportional parameters (for high frequency)  for all cylinders were set to zero in order to obtain the suitable fit between experimental and numerical results, once the impact velocity is low.Also, in order to not increase the simulation time, the stiffness proportional parameters were set to zero, because the stable time integration for explicit simulations could be affected [5].
The results for the force vs. time for the model with damping parameter and experimental test results are shown by Figure 11 The addition of damping improve the quality of the finite element model results in the time period after 0.0087 ms .On the other hand, the valleys and peaks in the beginning of simulation were not well captured.Despite the divergence in the beginning of the impact, the addition of damping allows the finite element model perform better predictions for the most of impact event.Also, the damping effects reduce the difference between the peaks and valleys, improving the high frequency part of numerical simulations.Figure 11(b) shows the displacement vs. time results, the differences between the finite element model and the experiments could be explained, in part, by the indentation phenomenon in the real structure.
The kinetic energy dissipated in the experiments were around 15.65 J (Table 3), regarding the finite element model, the damage and damping dissipate 13.37 J.  Regarding the damage and damping, the model did not dissipates enough energy and the difference could be due to delaminations, which are simulated by the material model.It is important to mention that the fiber breakage and delaminations dissipate more energy than matrix damage [19].
For type B cylinders, the force vs. time results for experimental tests and model are shown by Figure 13.For this configuration, the finite element model did not detect the maximum peak force value.Again for type B cylinder, the damping effects in simulations affect the capacity to fit the peaks and valleys, when the impact starts.Despite the divergence when the impact begins, for most of the impact event, the model can predicts the mechanical behavior for this cylinder configuration under impact loading.The same observations about Rayleigth's damping parameters made for Type A cylinders still valid for Type B cylinders.Without damping effects, the end of force vs. time would possess a sudden force drop.Figure 13(b) shows the displacement vs. time results, in this case, the finite element model predicted well the displacement, despite the indentation in the real structure.
Regarding the energy dissipated, the model dissipates more energy than the experimental tests ( Table 4 ).Even dissipating more energy than the experimental tests, the experimental impact tests last for more time (Figure 13).Even the finite element model dissipating more energy, the experimental tests last for more time.Figure 15 shows the force vs. time results for Type C cylinders.Again, the damping parameters improve the quality of the finite element model.However, to the model cannot capture the peaks and valleys in the initial part of the impact test, but for the rest of the impact event, the model provides good results.Figure 15(b) shows the displacement vs. time results, the differences between the finite element model and the experiments could be explained, in part, by the indentation in the real structure.
Table 5 shows the speed and the variation of the kinetic energy before and after the impact.As for Type B, the finite element model dissipates more kinetic energy than the experimental tests.
Although the model dissipates more energy than the experiments, the impact event duration for experiments is longer than for simulation.Figure 16(a) shows the simulation kinetic energy variation along the impact event.

CONCLUSIONS
Impact is a very complex event, with different factors can influence the structural behavior such as, material type, energy level, damage mechanism structure geometry, etc.This work showed that the cylindrical geometry adds complexity and it has a huge influence in the structural behavior.The damping effects, which can be negligible in flat panels, have an important influence in more complex geometries.
In order to confirm the real structure damping, some experimental tests must be performed in the future.
Also, a delamination criterion must be implemented.Despite the lack of a delamination criterion, the intra-ply model performs well in order to predict the damage and the type of damage in the cylinder structures.
Finally, the simulations showed that the damping effects possess a huge influence in the cylinder impact behavior (regarding the boundaries condition, energy level and cylinder thickness), once even without delamination and any other type damage (10 J impact), there are dissipative sources acting in the structure.

11 component in fiber direction 22 component in transverse direction 12 stress component 11 component in fiber direction 22 component in transverse direction 12 
Stress StressShear Strain Stress Shear (a).It is important to notice that the cylinders (Figure1(b)) were manufactured by CTM (Navy Technological Center-Brazil), using filament winding process for carbon fiber and epoxy resin.(a) (b) Figure 1: (a) Coupon dimensions; (b) Cylinder manufactured by filament winding process

Figure 3 :
Figure 3: Force vs. time and displacement vs. time for type A cylinders (30 J impact energy).
Figure3shows that the force increases rapidly close to 0.0052 s and then a sudden force drop occurs.After that, the force increases again and a new sudden drop occurs.This trend repeats for close to 0.0087 s of the impact event.A similar behavior was presented by Minak et al.[16].This part of peaks and valleys could indicate the initialization of delaminations in several layers.After this first time, the unstable delamination propagation can cause the further oscillations in the force vs. time history, as observed by Schoeppner and Abrate[19].Also, Figure3shows that the maximum force level does not occur in the first peak and the maximum force value does not occur in the same time of the maximum displacement.On the other hand, as mention before, the cylinders impacted under 10 J did not show any damage (Figure4(a)) as observed by the C-Scan images.However, the force vs. time history shows a similar trend as shown by specimens impacted under 30 J (Figure4(b)).

Figure 5 :
Figure 5: Force vs. time history for impactor and base for type B.

Figure 6 :
Figure 6: Force vs. time and displacement vs. time for type B cylinders (30 J impact energy).

Figure 7 :
Figure 7: Force vs. time and displacement vs. time for type C cylinders (30 J impact energy).

 and 12 
. To account this mutual influence, the damage variables d 2 and d 6 are function of the ply orientation angle (

Figure 8 :
Finite element model: (a) Cylinder geometry, boundary conditions and initial conditions for 30 J impact energy; (b) Impactor geometry, reference point, mass point and mesh.

Figure 9 :
Finite element mesh: (a) Type A; (b) Type B and C.

Figure 10 :
Figure 10: Type A cylinders: Simulation vs. Experimental results for (without damping effects).

Figure 11 :
Simulation vs. Experimental results for Type A cylinders: (a) force vs time (b) displacement vs. time.

Figure 12 (Figure 12 :
Figure 12(a) presents the kinetic energy along the impact and Figure 12(b) shows the damaged area due to impact.
Figure 13: (a) Force vs. time experimental and simulation results for Type B cylinders.(b) Displacement vs. time simulation and experimental results.

Figure 14 (Figure 14 :
a) shows the model kinetic energy loss and Figure 14(b) shows the damaged area.Simulation results for type B cylinder: (a) kinetic energy; (b) damaged area.Type C are thicker than Type B and Type A, these difference in thickness are due to manufacturing process.Besides, Type C cylinders possess less 90 o oriented layers than others.(a) (b) Figure 15: (a) Force vs. time experimental and simulation results for Type C cylinders.(b) Displacement vs. time simulation and experimental results.

Figure 16 :
Figure 16(b) shows the damaged area caused by the impactor.Simulation results for type C cylinder: (a) kinetic energy; (b) damaged area.

Table 3 :
Type A -experimental speed and kinetic energy.

Table 4 :
Type B experimental speed and kinetic energy

Table 5 :
Type C experimental speed and kinetic energy.