Quantified assessment of modelled impacts on asymmetric composite laminates

Abstract

Modern approaches to exploring the impact behaviour of composite components involve computer models such that the components are virtually tested. However, for these virtual tests to be trusted they require careful validation using experimental data. This study demonstrates a quantitative validation approach for confirming the effectiveness of a simplified impact model. A set of six asymmetric laminates was impacted at 30 J whilst their bottom surface was monitored using high-speed cameras. The camera images were processed to obtain full-field strain which was then used to validate the model by the orthogonal decomposition technique. The validation procedure has confirmed that the model provides reasonable predictions of impact behaviour and gives insights into how it could be improved. The overall approach demonstrates how virtual testing of composite components can be vigorously validated. Validated virtual testing can then provide a low-cost and rapid route to the certification of lightweight composite structures.

Keywords

Impact damage composite materials high-speed DIC finite element analysis quantitative validation

Introduction

The application of composite materials is vital for the reduction of emissions across a wide range of sectors. This is due to their use in lightweight vehicle structures, as well as being the principal structural material for the generation of clean electricity by wind. The reason for this stems from their microstructure, which typically combines high-strength fibres with a tough polymer matrix. However, this microstructure adds substantial complexity when considering the damage tolerance of composite structures, in particular when subjected to impact damage.¹ As such there is a need to characterise the behaviour of new material systems when developing composite structures. Historically this characterisation has been performed experimentally with great cost to manufacturers. With increasing levels of computational power, it is becoming possible to virtually test new materials to determine their suitability for applications.^2,3 These virtual tests often require high levels of complexity to simulate all the damage mechanisms involved.^4–6 The extent of this complexity tends to vary with the length scales at which the model is defined. The benefit of the virtual testing approach is that it allows for many more loading scenarios to be considered. It also reduces the amount of expensive experiments and thus is increasingly applied to more complex structures.⁷ However, for virtual test results to be used in place of experimental studies it is essential that the models that underpin the virtual testing are fully validated.

Validation is a key process to establish confidence in computational models and their outcomes, within the scope of their intended application.⁸ The process consists of comparison between model outcomes and physical measurements taken under equivalent conditions. This can be performed qualitatively, such as visually evaluating similarities or data trends,⁹ and quantitatively by employing statistical approaches, with latter approaches preferred to quantify the agreement between datasets. The quality and quantity of physical measurements can directly affect the outcomes of the validation process and consequently confidence in a model. Obtaining high quality physical measurements is important, but it is also necessary to efficiently utilise the data and extract as much information from as few experiments as possible. Simulations based on computational models produce complicated 3D maps, such as strain fields, which contain gigabytes of data. Advances in digital cameras have allowed similar 3D maps to be obtained from experiments with the aid of optical measurement systems, such as digital image correlation (DIC). These datasets can be presented as images and thus visually compared to validate model predictions. However, a direct quantitative comparison is challenging due to differences in scale, mesh and coordinate systems. It has been previously demonstrated that the quantitative comparison between predicted and measured strain fields can be performed by implementing orthogonal decomposition.¹⁰ This technique reduces the dimensionality of 2D and 3D datasets¹¹ into feature vectors, while preserving key strain features.

Models of composite impacts are often validated through the comparison of one-dimensional data such as force-time data from the impact machine.^12,13 More specific comparisons can be made through comparisons of final damage morphology, such as comparing the shape predicted by finite element analysis with X-ray computed tomography data.¹³ Other approaches focus on specific defects, such as delaminations predicted using cohesive layers in models with experimental data from ultrasonic C-scans.¹² More recently, methods of obtaining full-field stress data during impact events have become available through the use of high-speed DIC.¹⁴ High-speed DIC of composite impact has been used for a number of applications such as damage assessment¹⁵ as well as qualitative, side-by-side validation of models of damage propagation¹⁶ and quantitative validation of impact dynamics.¹⁷

Composite structures are typically laminate structures, where each lamina is orientated in a particular direction to impart desired strength characteristics to the composite. Whilst the laminae within the laminate can have highly varied orientations, in most applications they are arranged such that they are symmetric around the midplane of the laminate. This is to prevent complicated coupling effects, where in-plane forces result in out-of-plane deformation and vice-versa. There are however some situations where this effect is potentially beneficial, with particular attention given to it for aerodynamic and hydrodynamic applications. For example, turbine blades can be designed that twist to reduce loads at different rotation speeds,¹⁸ and lifting surfaces can deform to maintain an optimal angle-of-attack at different speeds.¹⁹ As these sorts of laminates are not widely applied there has been limited experimental study into the damage tolerance of asymmetric laminates and no quantitative validation of models of asymmetric laminates during impact.

In this study, finite element analysis is used to model asymmetric composite laminates as they are impacted to estimate deformation behaviour and damage creation. A novel technique is used to obtain the material properties for the model based on a small number of mechanical tests. The model is then validated through quantitative comparison with strain data captured using high-speed DIC of six specimens undergoing low-velocity impact. The paper is structured as follows, in the next section the experimental method for the impact experiment is described. Then the approach for obtaining the full elastic properties from a small number of experiments is introduced and used to provide the inputs for a finite element model. The results section explains the qualitative and quantitative validation of the model results, and the following section discusses the validation procedure as well as the strengths and limitations of the model. The paper is concluded in the final section.

Experimental method

Six asymmetric laminates with a stacking sequence of [0₂/90₂]₄ were produced using unidirectional prepreg (RP542-4, PRF Composite Materials, UK). This prepreg is a mid-temperature cure epoxy resin with 12K carbon-fibre tows, where the fibres have an areal weight of 300 g/m². The specimens were cured in a hot press operated with a pressure of 6 bar at a cure temperature of 120°C, according to the manufacturer’s recommendations. The pressure was first applied and then the specimens were heated to the cure temperature and held at this temperature for 60 minutes. After curing, the specimens were left to naturally cool to room temperature before pressure was removed. Finally, the laminates were cut down to 100 mm by 150 mm specimens using a wet diamond saw (Versatile Power Pro Bridge Saw, Vitrex, UK).

The asymmetric specimens were then impacted using a drop-weight impact machine (CEAST 9340, Instron, Italy) equipped with a 16 mm hemispherical tup and a mass of 4.215 kg. The tup is the part of the falling weight that makes contact with the specimen during the impact. All six specimens were impacted at an energy level of 30 J, a similar energy to that imparted by tool drops. During the impact events force data was measured by the tup at a rate of 2 MHz. A high-speed DIC system (Q-450, Dantec Dynamics, Germany), operated with two high-speed cameras (Phantom v711, Vision Research, USA) was used to monitor changes in strain fields as specimens were being impacted. Both cameras were fitted with 100 mm focal length lenses (Makro Planar 2/100 T*, Carl Zeiss, Germany) and were set to record the impact events at a frame rate of 10 kHz. The camera recording was triggered by a light gate built into the test machine. A speckle pattern was a key requirement for the DIC system to track changes on laminates as they were impacted. The pattern was created using a white base coat (Matt White, PlastiKote, USA), followed by black speckles (Carbon Black, Liquitex, USA) to introduce a random texture for optical tracking. Istra4D (Dantec Dynamics, Germany) was then used to process these captured images to calculate surface displacement data with a facet size of 25 pixels and grid spacing of 10 pixels.

Ultrasonic characterisation was performed on all impacted specimens to obtain damage information using a pulse-echo ultrasound machine (Midas-NDT, UK). This involved an ultrasonic transducer emitting ultrasound at a frequency of 5.88 MHz, chosen as it was the closest selectable frequency to the ply resonance frequency of 6 MHz. Time-of-flight C-scans were captured to obtain accurate measurements such as the size, shape and location of damage formed in the specimens after impact. During the scans the ultrasound machine operated in a raster pattern with a line spacing of 0.2 mm and measurements captured every 0.1 mm along each line.

Elastic properties of the composite material were required for the model. To obtain these, a set of 25 mm wide coupons was produced to determine material properties in the axial and transverse directions. Five specimens with a [0]₄ layup and five with a [90]₈ layup were produced using the same curing conditions as described previously. Eight layers of prepreg were used for the 90° specimens as these were expected to fail at much lower stresses than the 0° specimens. The additional thickness meant that the failure load was well within the effective range of the test machine’s load cell. A similar speckle pattern to that used on the impact specimens was applied to these coupons. The coupons were then quasi-statically loaded to failure using a servo-hydraulic load machine (8501, Instron, USA) according to ASTM D638 whilst axial strain was measured using a stereoscopic DIC system (Q-400, Dantec Dynamics, Germany). This DIC system consisted of two cameras (Manta G-505B, Allied Vision Technologies, Germany) located a metre away from the specimens to measure three-dimensional displacements. Both cameras had a resolution of 2452 by 2056 pixels and were equipped with a pair of 25 mm lenses (FL-CC2514-2M, RICOH, Japan). All specimens were recorded using an aperture of f8.0 and illuminated with a LED light source (Dantec Dynamics, Germany), ensuring good light exposure during the experiments. The same facet size and grid spacing as used during the impact experiments was used to obtain surface strain on the specimens. These tests provided measurements of the ultimate strength and elastic moduli in the axial and transverse directions of the composite material. These properties were then further processed using a method described in the next section to obtain the full elastic properties of the composite.

Finite element analysis of impact events

Estimation of orthotropic material properties

To model the deformation and failure of composite specimens during impact it is necessary to have a large number of material properties. It is common practice in academic studies for such properties to be taken from multiple sources, often for very different material systems. As this study involves quantitative validation it was necessary to obtain more accurate estimations of material properties. A full experimental characterization of a new material system would be prohibitively expensive. Thus, a hybrid approach that uses experimentally measured properties in orthogonal directions was used to obtain estimates of the full orthotropic material properties. The composite material properties for the volumetric stress state of a unidirectional carbon-fibre reinforced polymer layer can be estimated using the rule of mixtures, however, this does not take account of fibre shielding, where the rigid fibres result in low stresses in parts of the matrix when loads are transverse to the fibres.²⁰ Numerous methods exist for adjusting the rule of mixtures to account for shielding, for this study, a method that accounts for fibre shielding by defining a second fibre volume fraction for the transverse direction was used.²¹ Using this method for orthotropic materials, the calculated dependencies of elastic characteristics for a fibre-reinforced layer can be calculated as follows:

E_{1} = ψ_{1} E_{f} + \frac{(1 - ψ_{1}) (1 + ψ_{3})}{1 - ψ_{3}} E_{m}

(1)

E_{2} = \frac{(1 + ψ_{1}) (1 + ψ_{3})}{(1 - ψ_{1}) (1 - ψ_{3}) (1 - {ν_{f}}^{2})} E_{m}

(2)

E_{3} = ψ_{3} E_{f} + \frac{(1 + ψ_{1}) (1 + ψ_{3})}{(1 - ψ_{1}) (1 - {ν_{f}}^{2})} E_{m}

(3)

ν_{12} = \frac{ν_{f} (1 + ψ_{3}) (1 + ψ_{1})}{ψ_{1} (1 - ψ_{3}) (1 - ψ_{1}) (1 - {ν_{f}}^{2})} \cdot \frac{E_{m}}{E_{f}}

(4)

ν_{23} = ν_{13} = ν_{f} ψ_{3} + (1 - ψ_{3}) ν_{m}

(5)

G_{12} = \frac{1 + ψ_{1}}{(1 - ψ_{1}) (1 + ψ_{3})} G_{m}

(6)

G_{23} = \frac{1 + ψ_{3}}{(1 - ψ_{3}) (1 - ψ_{1})} G_{m}

(7)

G_{13} = \frac{(1 + ψ_{1}) (1 + ψ_{3})}{(1 - ψ_{1}) (1 - ψ_{3})} G_{m}

(8)

where

E_{i}

are the moduli of elasticity in the three primal directions,

G_{i j}

are the shear moduli and

ν_{i j}

, are the Poisson’s ratio for the three primal planes.

ψ_{1}

ψ_{3}

are the relative volumetric contents of reinforcement in the direction of the fibres, and normal to them.

E_{m}

and

E_{f}

are the moduli of elasticity of the matrix and fibre respectively,

G_{m}

is the shear modulus of the matrix, and

ν_{m}

is the Poisson’s ratio of the matrix. The reinforcement coefficient

ψ_{1}

characterizes the relative volumetric content of fibers and thus is equivalent to the volume fraction of the composite. The value of

ψ_{3}

is determined using empirical dependencies and varies such the

ψ_{1} / ψ_{3}

is typically between 0.05 and 0.15.

The elasticity ratios for an orthotropic unidirectional reinforced layer in its symmetry axes considering the engineering constants in matrix form are given by:

ε = B σ

(9)

where,

B = [\begin{array}{l} \frac{1}{E_{1}} & - \frac{ν_{21}}{E_{2}} & - \frac{ν_{31}}{E_{3}} & 0 & 0 & 0 \\ \frac{1}{E_{2}} & - \frac{ν_{32}}{E_{3}} & 0 & 0 & 0 \\ \frac{1}{E_{3}} & 0 & 0 & 0 \\ \frac{1}{G_{23}} & 0 & 0 \\ s y m . & \frac{1}{G_{13}} & 0 \\ \frac{1}{G_{12}} \end{array}]

(10)

and

σ

and

ε

are the stress and strain tensors in Voigt notation. Elastic properties were determined numerically based on varying the values of

ψ_{1}

and

ψ_{1} / ψ_{3}

within the range [0.05, 0.15], such that the values of

E_{1}

and

E_{2}

corresponded to the values determined experimentally. For this it was necessary to specify the mechanical properties of the constituents of the composite. The properties of the constituents of composite materials do not vary as significantly from composite to composite and thus the values in Table 1 were used. The resulting orthotropic properties for the plies after numerical iteration are given in Table 2. The standard deviation of

E_{1}

and

E_{2}

were propagated through equations (1) to (10) to obtain the uncertainties on the elastic constants, also shown in Table 2.

Table 1.

Input data for determining the composite elastic constants.

$E_{f}$ , (GPa)³⁶	350
$ν_{f}$ , (−)³⁷	0.27
$G_{f}$ , (GPa)^a	138
$E_{m}$ , (GPa)³⁸	3.29
$ν_{m}$ , (−)³⁸	0.32
$G_{m}$ , (GPa)^a	1.25

^aShear moduli were calculated from elastic moduli and Poisson’s ratios.

Table 2.

Elastic properties for a composite layer with orientation of 0°. The numbers in parentheses are the uncertainties of the elastic constants. The values for E₁ and E₂ were experimentally measured, whilst all other values were estimated based on the method described in the finite element section.

E₁, (GPa)	129.7 (8.61)
E₂, (GPa)	7.763 (0.69)
E₃, (GPa)	12.38 (3.25)
ν₁₂, (−)	0.220 (0.0017)
ν₁₃, (−)	0.345 (0.0006)
ν₂₃, (−)	0.320 (0.0004)
G₃₂, (MPa)	2017 (65.6)
G₁₂, (MPa)	2605 (12.05)
G₃₁, (MPa)	1724 (26.30)

Modelling of impact event

The experiment was conducted in accordance with ASTM D7136/D7136 M and thus the specimen on the impact fixture was modelled in Ansys (2024 R1, Ansys Inc., Pennsylvania) with the following assumptions. The fixture base was assumed to be a steel plate with dimensions of 300 mm in width, 300 mm in length, and thickness of 20 mm, with a rectangular planar cut-out of 75 mm by 125 mm. In the experiment the specimen was held over the cut-out by four rubber feet, for simplicity these were neglected. The impact tup was modelled as a steel hemisphere.

The centroids of the test specimen, the cut-out in the fixture base, and the steel hemisphere were aligned along the direction of the impact direction. At the beginning of the model’s simulation, the laminate made contact with the fixture base through hinged movable supports along the perimeter of the contact line between the specimen and the rectangular cut-out of the support fixture. The contact was assumed to be frictionless. The fixture base was rigidly clamped along its outer edges. Delamination between the eight plies was simulated through the interface delamination failure model in Ansys.^22–24 The possibility of crack formation through each ply from the centre of the specimen along the reinforcement directions was simulated using debonding model interfaces based on separation distance.²⁵ The impact load on the system was assumed to arise kinematically through the corresponding displacement of the flat boundary of the hemisphere. A static simulation was performed such that the dynamic properties of the material did not need to be determined.

For meshing, the specimen was divided into a set of 3D finite elements consisting of 81,500 hexahedral SOLID186 elements for all bodies except the tup and 9706 tetrahedral SOLID187 elements for the tup hemisphere. The SWEEP method and path confirm method were employed, with element counts specified in the respective directions. Contact between pairs of opposing surfaces of the test specimen with the fixture base and the test specimen with the tup hemisphere were simulated using 4902 contact elements CONTAC174 and 4902 target elements TARGE170. Additionally, 26,600 interface elements INTER204 were defined between the unidirectional layers of the test specimen to model the regions susceptible to delamination using the cohesive zone modelling method built into the package. These delaminations were anticipated between layers reinforced in different directions and within the plane formed by the reinforcement direction of the layer and the normal to the opposing surface, passing through the centre of the corresponding opposing surface. The elastic properties of the specimen were taken from Table 2. The hemispherical tup and fixture were modelled as a generic steel, with an elastic modulus of 250 GPa and Poisson’s ratio of 0.33. In total, the mesh consisted of 459,562 nodes and 91,206 elements.

The final number of elements was determined by an iterative process of gradually refining the mesh, ensuring convergence of the results of the stress-strain state compared to the results of reactions of the tup. The computational mesh of the structural model has been created using sweep and tetrahedron methods, allowing for a complete decomposition of the model into finite elements. The sweep method predominates in the model (used for all the elements except the impact tup), resulting in a structured mesh. The placement of nodes in the structured mesh was clearly defined, predictable, and implicitly linked from a programming perspective. This enabled efficient programming through data storage and led to efficiency in execution time and utilisation of system resources.²⁶ Quadratic elements utilizing a nonlinear shape function have been employed. Displacements between nodes were interpolated using higher-order polynomials. Such elements are better suited for representing complex geometries and bending deformations.

Due to the use of the tetrahedrons method, an unstructured tetrahedral mesh has been constructed on the impact tup. Employing the patch conforming method approach, this bottom-up approach involves building the mesh starting from the edges of faces and then moving into the volume. All faces and their boundaries adhere to and connect, creating a high-quality, clean mesh. The primary advantage of unstructured meshes lies in their ability to be generated on complex geometries, often absent in structured grids. The main criteria for mesh quality evaluation included orthogonal quality and skewness. It was recommended to use mesh models where the minimum orthogonal quality value was greater than 0.1 or the maximum skewness value was less than 0.95.²⁷

The overall quality of elements was assessed based on the ratio of volume to the sum of the squares of edge lengths for two-dimensional quadrilateral and triangular elements. For three-dimensional elements, it was the ratio of volume to the square root of the cube of the sum of the squares of edge lengths. In the model, the dimensional range of elements varied from 0.22842 to 1, where a value of 1 indicated an ideal cube or square, while a value of 0 indicated that the element had zero or negative volume. The aspect ratio for the sides of a triangle provides a comparison of the height and width of the triangle, ranging from 1 to infinity. For quadrilaterals, the aspect ratio was defined as the ratio of the long side to the short side of the quadrilateral. In the considered problem, this ranged from 1 to 8. The finalized mesh is shown in Figure 1.

Figure 1.

Computational mesh viewed from above (top), the side (middle), and a close up view of a cross-section through the tup in contact with the laminate (bottom).

Although the static structural problem has been selected, a time step size parameter was still used. The time-step size parameter delineates the intervals in physical time at which the model conducted steady-state simulations. Fully implicit methods exhibit first-order accuracy. Thus, a reduction in the time-step size proportionally linearly enhances the accuracy of the model. While models with larger time-step sizes are invariably stable and physically plausible, they may lack temporal precision in physical meaning. In addition, the number of inner iterations within each time-step for the steady-state simulations was resolved using iterative algorithms. At each temporal level, a specified number of iterations was necessary to achieve convergence. Consequently, the process spans from time 0 to a unified 1 s interval. The finite element model was subject to the following boundary and initial conditions at the hypothetical start time of the simulation. At time 0 s, the system was in an unstressed and undeformed state, with no delaminations present. Over time, the hemispherical tup moved purely along the inward normal direction to the upper face of the test specimen linearly and uniformly. By 1 s, the displacement was 6 mm. The simulation of flaw development has been conducted based on the cohesive zone material (CZM) failure model²⁸ which models interface delamination and progressive failure for two materials joined together. Fracture or delamination, occurring at the interface between phases, significantly impacts the toughness and flexibility of multi-phase materials, particularly those found in laminated composite structures. This method explicitly integrates fracture mechanisms by utilizing softening relationships between traction and separation, thereby introducing critical fracture energy. This concept represents the energy required to disrupt the interface surfaces. The interface surfaces of materials can be represented using a specialized set of interface elements or contact elements, with the CZM model employed to describe the constitutive behaviour of the interface. The CZM model encompasses a constitutive relation between the traction (

T

) acting on the interface and the corresponding interfacial separation (

δ

), which denotes the displacement jump across the interface. Index

n

was used for normal process indicators and

t

for tangential ones. Interlaminar delaminations have been simulated with a mixed-mode bilinear cohesive zone material model including bilinear cohesive law under the mixed-mode fracture. For this, the separation of material interfaces depended on both the normal and tangential components of displacement jumps, considering the difference in their contributions to the separation of material interfaces. The non-dimensional effective displacement jump

λ

for mixed-mode fracture was defined as²⁷:

λ = \sqrt{{(\frac{δ_{n}}{δ_{n}^{c}})}^{2} + β^{2} {(\frac{δ_{t}}{δ_{t}^{c}})}^{2}}

(11)

where the constants for bilinear interface delamination were taken from studies using similar materials^29,30 and are shown in Table 3.

Table 3.

Constants for bilinear interface delamination.

Maximum normal traction	$T_{n}^{\max}$ (MPa)	250
Normal displacement jump at completion of debonding	$δ_{n}^{c}$ (µm)	2.888
Maximum tangential traction	$T_{n}^{\max}$ (MPa)	95
Tangential displacement jump at completion of debonding	$δ_{t}^{c}$ (µm)	22.80
Ratio of displacements jumps	$α = \frac{δ_{n}^{*}}{δ_{n}^{c}}$ (−)	0.1111
Non-dimensional weighting parameter	β (−)	1 × 10⁻⁵

Cracks through the plies have been modelled using the bilinear behaviour contact approach. The bilinear cohesive zone material model for mode one debonding delineates a mode of interface surface separation where the normal separation predominates over tangential slip along the interface. The material behaviour model simulates linear elastic loading followed by linear softening. The maximal point at the end of linear elastic loading marks the attainment of maximum normal contact stress, initiating debonding, which concludes at the end of linear elastic loading as the normal contact stress diminishes to zero. Further separation then proceeds without any normal contact stress. The area beneath the curve represents the energy liberated during debonding, denoted as the critical fracture energy. The slope of the first line determines the contact gap at maximum normal contact stress. It thereby characterizes the reduction of normal contact stress with increasing contact gap, indicating the brittleness or ductility of the fracture. Following initiation of debonding, it was assumed to be cumulative, with any unloading and subsequent reloading occurring linearly along the second line at a gentler slope. The constants for separation-distance based debonding were taken from^29,30 and are shown in Table 4, no tangential slip was allowed when under compression. As the values in Table 4 were taken from literature, a simple sensitivity analysis was performed to identify how much they effect model results. Each parameter was varied by +/− 5% and the maximum strain on the surface recorded. The normalised displacement jump had the most effect on predicted strains, causing variations of around 100 µε, whilst the other parameters had smaller effects.

Table 4.

Mode 1 debonding Interface constants.

Maximum normal contact stress	$σ_{\max}$ (MPa)	57
Contact gap at the completion of debonding	$u_{n}^{c}$ (µm)	22.8
Power law exponent for mixed-mode debonding	(−)	2
Artificial damping coefficient	d_V (s)	5 × 10⁻⁴

Results

Six asymmetric laminates were impacted at 30 J. During these impacts the test machine recorded the force data from the moment the tup first contacted the specimen until after the tup left the specimen. This force data demonstrated that there was significant variation in the dynamics of the impact event, with specimens typically having a similar force-time profile but with variation in peak force and impact duration, see Figure 2. This variability in impact behaviour was also seen in the ultrasound scans performed after the impact event. Using these scans a typical impact damage morphology can be seen where delaminations are created at the interfaces within the laminate, increasing in size with depth, as shown in Figure 3. However, there was variation in the projected damage area and some delaminations near the top and bottom surfaces of the specimens appeared irregular in shape, an example of this irregularity is shown at the bottom of Figure 3.

Figure 2.

Force-time data captured by the impact machine showing variability in specimen behaviour.

Figure 3.

Ultrasound scans for two specimens impacted at 30 J, corresponding with specimen 1 (top) and specimen 6 (bottom) from Figure 2.

The simulation process was divided into steps according to a unified time scale from 0 to 1 s, where at 1 s the displacement applied by the simulated tup was 6 mm. Initial approaches for comparison between experimental measurements and model predictions focused on comparing damage visible on the rear surface of the specimens. Finding an adequate method for measuring the length of the crack on the back surface was a non-trivial task. A visual optical method of flaw detection was chosen. To avoid chemical damage to the specimen surfaces, the option of colour flaw detection was discounted. To verify the adequacy of the measurements, two techniques were selected for obtaining the crack lengths in the specimens: using a digital microscope and using a photo scanner. For the digital microscope, the length was determined by an aperture grid with an accuracy of ±15 µm, based on accumulated error of three measurements. Whilst the photo scanner had a resolution of 1200 dpi, equivalent to 47 px/mm. The dimensions measured using the photo scanner were verified by direct measurement with a calliper of the longer sides of the specimens with an error of ±5 µm, according to the reported precision of the instrument. The scanned digital images did not allow for the exact position of the crack tip to be determined, with the crack tip only visible if it happened to occur on a black speckle. Considering the aforementioned uncertainties, the overall uncertainty for the crack lengths amounts to ±15 µm. Across all the specimens, the mean crack length on the bottom surface was 2.03 mm with a standard deviation of 6.04 mm. An exemplar scanned specimen with the rear surface crack marked is shown in Figure 4.

Figure 4.

Bottom surface of a specimen after impact showing the speckle pattern and with the extent of the bottom surface crack indicated (top) and zoomed in view of crack tip (bottom).

The simulated crack length value of the specimen was determined at the moment the tup had been displaced by 6 mm in the model. This corresponds to the average value of the measured vertical displacements of the tup upon impacting all specimens. The crack opening as measured between the nearest nodes of the interface elements, whose mutual displacement exceeds the contact gap at the completion of debonding $u_{n}^{c}$ , given in Table 4. The crack length on the bottom surface of the modelled laminate was 72.86 mm. This is 1.79 standard deviations greater than the mean measured crack length.

Due to the variability encountered during the experiments, it was decided to base further model validation efforts on the strain field on the bottom surface of the specimens, which was measured using digital image correlation during the impact event. To study deformation behaviour throughout the duration of the impact event and validate the model, the interest lies in the model predictions at specific displacements of the impact tup hemisphere in the simulation at 2, 3, 4, and 5 mm. For the experiments the time these displacements were reached was estimated by twice integrating the force data measured by the load machine with respect to time. An example of the tup displacement is shown for one of the specimens in Figure 5, where the instances that were modelled are shown as square markers. Thus, the displacement of the hemisphere and the actual time elapsed since the contact of the tup with the laminate were obtained and the specific DIC strain fields at these times were selected for comparison with the strain predicted by the simulation. These specific instances are shown in Figure 5 as the square markers. For validation, strains from the simulation and experiment were converted to the planar Lagrangian principal strains. This measure of strain was used to prevent any influence on results due to the different coordinate systems of the simulation and the DIC system. A visual comparison between the simulation and one of the experiments at a displacement of 2 mm is shown in Figure 6.

Figure 5.

Displacement-time data captured by the impact machine for a specimen where positive displacement indicates the top surface of the specimen has been pushed down. The squares indicate the instances that were modelled and compared with the DIC data.

Figure 6.

First (left column) and second (right column) principal strains captured by the DIC system at a tup displacement of 2 mm (top row) and corresponding model predictions (bottom row).

Whilst there are differences between the experimental and modelled strain fields it can be difficult to determine how much of these differences are due to measurement uncertainty in the DIC data. As evident in Figure 6 there are also some visual differences due to the model data being plotted as contours with a finite number of levels resulting in blocks of colour whilst the DIC data is more continuous. Thus, to avoid subjective comparisons a quantitative validation approach was used based on orthogonal decomposition. The mathematics for performing decomposition are described in detail in Ref. 10 but are briefly summarised here. In this method, the strain fields, $I$ , were first decomposed into feature vectors, $f$ , by projecting them onto a set of discrete Chebyshev polynomials³¹:

f_{i} = \frac{1}{ρ_{i} \sqrt{M N}} \sum_{x, y = 1}^{M, N} I (x, y) T_{i} (x, y)

(12)

Where

f_{i}

is the

i

th coefficient in the feature vector,

T_{i}

is the

i

th Chebyshev kernel calculated from the recursive function given in Ref. 31,

I

is the strain field that is to be decomposed,

ρ_{i}

is a normalisation constant equal to the Euclidean norm of the

i

th Chebyshev kernel.

M

and

N

are the number of rows and columns of pixels in the strain field, and

x

and

y

are the pixel locations in the strain field. Decomposition results in the dimensionality of the data being reduced whilst retaining all the important information regarding the distribution of strain. When performing decomposition, it was necessary to choose the number of polynomials used during the process, with higher numbers meaning that when the data was reconstructed it was more similar to the original dataset. For this study the number of polynomials used was chosen by ensuring the representation error of the data, calculated by performing the inverse of equation (12), was less than the measurement uncertainty of the DIC system, found to be 459 µε in a previous study.¹⁵ This resulted in feature vectors containing on average 514 terms, compared to 6000 strain vectors in the experimental data and 91,206 elements from the simulation.

The experiments were repeated six times and thus for each displacement interval six comparisons could be made. However, there are multiple formulations of strain that could be used during the validation process, with the first and second principal strains being the leading contenders as they fully define the strain state at every surface location, as shown in Figure 6. To choose which formulation to use during the validation process, it was decided to choose the quantity where the deviations between the simulation and experiment were greatest. The feature vector difference, defined as the Euclidean distance between the pairs of model and experiment feature vectors,³² was used to determine when the deviation was greatest. It was found that the deviation between the model and experiment increased as the displacement was increased and that it was consistently higher for the first principal strain, shown in Figure 7. Thus, quantitative validation was performed using the first principal strain only.

Figure 7.

Increase in feature vector difference due to greater uncertainties at higher tup displacements for both first (red) and second (blue) principal strains. The trends are indicated by a pair of linear lines-of-best-fit (LoBF).

Two approaches to validation were used in this study. The first provides a simple binary output regarding whether the deviation of the model predictions from the experiment can be explained by the measurement uncertainty of the DIC system. This approach is described in Ref. 10 and involves plotting the pairs of coefficients from the two feature vectors to be compared against each other, shown in Figure 8. If each coefficient in the simulation feature vector is the same as the corresponding coefficient in the experiment feature vector then it would lie on a straight line,

f_{s, i} = f_{e, i}

(13)

where

f_{s, i}

and

f_{e, i}

are the

i

th coefficients in the model and experiment feature vectors respectively. However, the uncertainties associated with the measurement and data processing mean this condition is impossible and thus a region called the expanded uncertainty bounds is defined,

f_{s, i} = f_{e, i} \pm 2 u_{\exp}

(14)

where

u_{\exp}

is the expanded uncertainty given by,

u_{\exp} = \sqrt{{u_{c a l}}^{2} + {u_{d e c o}}^{2}}

(15)

u_{c a l}

is the measurement uncertainty, and

u_{d e c o}

is the representation error for the experimental feature vector. Once the bounds represented by equation (13) were plotted it was found that the majority of coefficient pairs sit within the expanded uncertainty bounds but that typically one to three were outside the bounds, examples for four different displacements for the same specimen are shown in Figure 8. This means that the deviations cannot be explained through measurement uncertainty alone. Four validation diagrams out of a total of 24 indicated that the model compared to certain experimental datasets was valid given the measurement uncertainty of the DIC system. These were all for the 2 mm displacement simulations.

Figure 8.

Validation diagram showing how coefficients are more scattered within the expanded uncertainty (Exp. Unc.) bounds as tup displacement increases.

To quantify the deviations in relation to the expanded uncertainty a validation metric has previously been defined.³³ In this study the metric was implemented by first normalising the errors between the experiment and model coefficient pairs as,

w_{i} = 100 \cdot \frac{| f_{s, i} - f_{e, i} |}{\sum_{i = 0} | f_{s, i} - f_{e, i} |}

(16)

Next, the expanded uncertainty was similarly normalised to obtain a threshold value,

e_{t h} = 100 \cdot \frac{2 \cdot u_{\exp}}{\sum_{i = 0} | f_{s, i} - f_{e, i} |}

(17)

Finally, the normalised errors that are less than the threshold given by equation (16) were summed. This results in a scalar value that is 100% if every normalised error is less than the threshold indicating that the simulation is valid given the measurement uncertainty of the DIC system. For this study, the validation metric was calculated for each data pair, i.e. 24 times, and was found to have a median value of 88.9%. The interquartile range, indicating the spread of values of the metric, was 9.2%.

Discussion

Whilst the model was validated when considering some of the simulated instances of the impact, the differences for the majority of instances were greater than can be explained by the measurement uncertainty of the DIC system. This is to be expected as the complex microstructure of composite materials and the variability in impact behaviour means that it is challenging to accurately simulate what will happen to a specimen. The variability in the impact behaviour of composites is well reported and typically stems from the material microstructure as well as slight variations in impact angle and velocity.³⁴ Thus, the intention of this study was not to generate a highly complex accurate model, but to see how well a simplified model based on cohesive zone and brick elements can simulate typical damage creation in composites. Such a model will allow scaling up to layered composite thin-walled and medium-thickness shell structural elements in the aerospace industry, enabling the simulation of low-velocity collision problems of foreign object damage. Thus the approach used here would only be suited to simulating impact energies that result in barely visible impact damage. Higher velocity impacts would require reducing the interaction time and thus necessitate accounting for wave effects in the stress-strain state propagation through the layered shell. This would mean solving the transient problem or the explicit dynamic problem. Under these conditions, further research and the creation of appropriate non-classical rheonomic CZM models would be required. Overall, experimental data show that the simulation reliably estimates the extent of crack formation and delamination using the semi-inverse modelling method even with reduced interaction time, as material stiffness increases with strain rate at least in the initial deformation stage. The complexity of a model can be set at any level, with the most basic being simulating the whole laminate as a single layer of orthotropic bricks. Recent techniques such as extended finite element models allow damage propagation to be simulated at the fibre length scales however are computationally intensive.³⁵ The mechanisms that are simulated can also vary. In this study inertia and strain rate effects have been neglected as the impact velocity was just 3.773 ms⁻¹, well below the velocity at which dynamics affect damage propagation.³⁴ Specifying the required complexity of a model is often based on maximising it based on the computational resources available. The approach used to characterise the performance of the model in this study allows an optimal balance between model complexity and confidence in model results to be determined.

The validation metric demonstrates that all strain fields simulated by the model are close to the measured ones. This suggests that by making small modifications to the model it may be possible for it to be fully validated. There are a lot of potential areas where the model could be modified, such as reducing the number of assumptions. For example, it was assumed that transverse shear and compression stresses at the contact boundary are equal. These stresses are considered here because they are complementary to the internal force components in the tangential and normal linear displacements relative to the interface surface, respectively. The rigid layer contact was achieved through interface elements with CZM behaviour in a zero-thickness layer in the delaminated area. It was assumed that in some local areas of the layer, interface elements provide one-sided contact between the rigid layers based on a bilinear model: a linear trend between the linear displacements and complementary internal force stress factors unfolds and continues to decrease linearly to zero. An approach to improving the cohesive zone is through the possible replacement of the bilinear failure model with an exponential one.

Regarding crack formation in the quasi-homogeneous layer, it should be noted that the influence of tangential stresses through the application of the mixed debonding interface is not obvious. Although it requires further research, the deformation state shows the minor influence of such stresses on average and thin (in terms of shell class) layer thicknesses. The expected impact may concern local stability losses of a limited edge area due to the explanation of the cracked face surface.

For obtaining the elastic properties for the simulation, an approach based on experimental measurements and laminate theory was able to obtain good material properties. However, no such approach exists for material properties associated with damage propagation. Determining and experimentally refining the constants of the CZM and debonding interface models is thus needed for a more accurate simulation. These vary in each experimental sample, and the individual determination of such constants should significantly refine the simulation results, especially with a developed deformation pattern. Changes to these constants by small amounts could yield approximately similar stress-strain states and crack lengths, allowing for an upper-bound estimate of the maximum strength of laminate constructions.

The work described in this study is the confluence of two highly complex efforts. The high-speed measurements capturing the strains on the rear surface of a specimen during impact are non-trivial to obtain due to the expense of the measurement systems involved and test speeds. Similarly, due to the complex mechanics involved under impact of composites, there are a lot of potential different ways of modelling the impact event. By utilising these experimental measurements, it is easier to trial different approaches to modelling composite impact and determine the best approach whilst minimising the computation power required. Once such a modelling approach is found and validated, it is then possible for the model to be modified to more closely represent structural details which are too challenging to experimentally impact.

Conclusions

Ensuring impact tolerance of composite materials is essential for the application of laminates in aerospace and clean energy generation structures. However, exploring the damage behaviour of composite laminates under impact is a highly challenging procedure requiring substantial effort from the modeller as well as computational resources. In this study quantitative validation has been used to confirm that simplified modelling schemes can predict the damage behaviour of composite laminates during the early stages of impact events. A set of six asymmetric carbon fibre reinforced laminate specimens were manufactured and then impacted whilst their bottom surface was monitored using a high-speed digital image correlation system. The same laminates were modelled using finite element analysis to yield predictions of the first and second principal strain fields during the impact. Orthogonal decomposition was used to quantitatively confirm the validity of the model predictions during the early stages of the impacts, with the differences between model and experiment increasing as the extent of damage created in the specimens increased. By confirming the extent of validity for the simulations it is possible to identify the simulation parameters that need refinement. This makes it simpler to achieve the right balance between model fidelity and computation time, helping to reduce development costs for new composite structures.

Footnotes

Acknowledgments

This work is the result of a visit by Dmytro Zhyhylii and Oleksandr Pozovnyi to the University of Liverpool in 2023 as part of a Universities UK International’s (UUKi) UK-Ukraine R&I twinning grant funded by Research England. We are grateful to Professor Igor Potapov of the University of Liverpool for his work in securing the grant.

Author Contributions

Dmytro Zhyhylii: Conceptualization, Investigation, Writing - Original Draft. Oleksandr Pozovnyi: Conceptualization, Investigation, Writing - Original Draft. Eric Y.H. Chai: Investigation. Ksenija Dvurecenska: Conceptualization, Supervision, Writing - Review & Editing. William J.R. Christian: Conceptualization, Supervision, Validation, Writing - Review & Editing.

Declaration of Conflicting Interest

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: Visit by Dmytro Zhyhylii and Oleksandr Pozovnyi to the University of Liverpool was funded by Universities UK International’s (UUKi) UK-Ukraine R&I twinning grant funded by Research England. Eric Y.H. Chai received a scholarship from the University of Liverpool and National Tsing Hua University dual PhD programme.

ORCID iD

William JR Christian

Data availability statement

Data available on request.*

References

Christian

WJR

DiazDelaO

Patterson

. Strain-based damage assessment for accurate residual strength prediction of impacted composite laminates. Compos Struct 2017; 184: 1215–1223. DOI: 10.1016/j.compstruct.2017.10.022.

Hallett

Green

Jiang

W-G

, et al. The open hole tensile test: a challenge for virtual testing of composites. Int J Fract 2009; 158: 169–181. DOI: 10.1007/s10704-009-9333-8.

Okereke

Akpoyomare

Bingley

. Virtual testing of advanced composites, cellular materials and biomaterials: a review. Compos B Eng 2014; 60: 637–662. DOI: 10.1016/j.compositesb.2014.01.007.

Bruno

Esposito

Papa

, et al. Design, manufacturing, and numerical characterization of hybrid fiber reinforced polymer under dynamic loads. J Mater Eng Perform 2023; 32(9): 3905–3913. DOI: 10.1007/s11665-023-07926-1.

Falcó

Lopes

Sommer

, et al. Experimental analysis and simulation of low-velocity impact damage of composite laminates. Compos Struct 2022; 287: 115278. DOI: 10.1016/j.compstruct.2022.115278.

Lee

, et al. Damage characterization and modelling of FRP laminated composites subjected to external edge-on impact. Compos Part B – Eng 2024; 287: 111879. DOI: 10.1016/j.compositesb.2024.111879.

Guerrero

Sasikumar

Llobet

, et al. Experimental and virtual testing of a composite-aluminium aircraft wingbox under thermal loading. Aero Sci Technol 2023; 138: 108329. DOI: 10.1016/j.ast.2023.108329.

ASME V&V 10-2019 . Guide for verification & validation in computational solid mechanics. New York: American Society of Mechanical Engineers.

Papa

Bruno

Napolitano

, et al. Numerical and mechanical analysis of laminated PA11/twill basalt composites with enhanced flame behavior. Compos Struct 2024; 338: 118123. DOI: 10.1016/j.compstruct.2024.118123.

10.

Sebastian

Hack

Patterson

. An approach to the validation of computational solid mechanics models for strain analysis. J Strain Anal Eng Des 2013; 48: 36–47. DOI: 10.1177/0309324712453409.

11.

Amjad

Christian

WJR

Dvurecenska

, et al. Quantitative comparisons of volumetric datasets from experiments and computational models. IEEE Access 2023; 11: 123401–123417. DOI: 10.1109/ACCESS.2023.3328326.

12.

Chen

Cui

, et al. Low-velocity impact analysis and multi-objective optimization of hybrid carbon/basalt fibre reinforced composite laminate. Compos Struct 2024; 343: 118305. DOI: 10.1016/j.compstruct.2024.118305.

13.

Wang

Wei

, et al. Low-velocity impact responses and failure of sandwich structure with carbon fiber composite honeycomb cores. Int J Impact Eng 2024; 192: 105034. DOI: 10.1016/j.ijimpeng.2024.105034.

14.

Flores

Mollenhauer

Runatunga

, et al. High-speed 3D digital image correlation of low-velocity impacts on composite plates. Compos B Eng 2017; 131: 153–164. DOI: 10.1016/j.compositesb.2017.07.078.

15.

Chai

EYH

Wang

Christian

WJR

. Experimental study on the effect of waviness defects on composite material impact dynamics. Compos B Eng 2024; 283: 111638. DOI: 10.1016/j.compositesb.2024.111638.

16.

Zhu

Song

, et al. Ballistic performance and energy transfer of ultrahigh molecular weight polyethylene laminate. Mater Lett 2024; 367: 136540. DOI: 10.1016/j.matlet.2024.136540.

17.

Burguete

Lampeas

Mottershead

, et al. Analysis of displacement fields from a high-speed impact using shape descriptors. J Strain Anal Eng Des 2014; 49: 212–223. DOI: 10.1177/0309324713498074.

18.

Bagherpour

Manolas

, et al. Modeling of material bend-twist coupling on wind turbine blades. Compos Struct 2018; 193: 237–246. DOI: 10.1016/j.compstruct.2018.03.071.

19.

Liao

Martins

JRRA

Young

. 3-D high-fidelity hydrostructural optimization of cavitation-free composite lifting surfaces. Compos Struct 2021; 268: 113937. DOI: 10.1016/j.compstruct.2021.113937.

20.

Hull

Clyne

. An introduction to composite materials. 2nd ed. Cambridge: Cambridge University Press, 1996.

21.

Vereshchaka

. Упругие постоянные и обобщенные xарактеристики жесткости материала слоистой структуры [Elastic constant and generalized stiffness characteristics of the material of a layered structure]. Bulletin of the National Technical University “KhPI”. Dyn. Strength Mach. 2008; 36: 29–40. https://repository.kpi.kharkov.ua/handle/KhPI-Press/33900.

22.

Zhaohui

Xiaodong

Jun

, et al. Study on delamination problems of composite hull structures under slamming loads. Polym Polym Compos 2011; 19: 433–438. DOI: 10.1177/0967391111019004-528.

23.

Jadhav

. Modeling of mode 1 delamination behavior in laminated composite structures for aerospace applications. Mater Today Proc 2023; 82: 61–67. DOI: 10.1016/j.matpr.2022.11.341.

24.

Jadhav

Yella

Azmeera

, et al. Mode II- delamination modeling in multi-directional laminates using cohesive zone methodology under static and dynamic loading. Mater Today Proc. 2023. DOI: 10.1016/j.matpr.2023.01.227, In press.

25.

Kubair

Spearing

. Cohesive zone model for direct silicon wafer bonding. J Phys D Appl Phys 2007; 40: 3070. DOI: 10.1088/0022-3727/40/10/010.

26.

Paluszny

Matthäi

Hohmeyer

. Hybrid finite element–finite volume discretization of complex geologic structures and a new simulation workflow demonstrated on fractured rocks. Geofluids 2007; 7: 186–208. DOI: 10.1111/j.1468-8123.2007.00180.x.

27.

Ansys . ANSYS help, 2024 R1. 15 ed. ANSYS, Inc, 2023.

28.

Alfano

Crisfield

. Finite element interface models for the delamination analysis of laminated composites: mechanical and computational issues. Int J Numer Methods Eng 2001; 50: 1701–1736. DOI: 10.1002/nme.93.

29.

Gao

Giannopoulos

, et al. Progressive damage numerical modelling and simulation of aircraft composite bolted joints bearing response. Materials 2020; 13: 5606. DOI: 10.3390/ma13245606.

30.

Benzeggagh

Kenane

. Measurement of mixed-mode delamination fracture toughness of unidirectional glass/epoxy composites with mixed-mode bending apparatus. Compos Sci Technol 1996; 56: 439–449. DOI: 10.1016/0266-3538(96)00005-X.

31.

Mukundan

Ong

Lee

. Image analysis by Tchebichef moments. IEEE Trans Image Process 2001; 10: 1357–1364. DOI: 10.1109/83.941859.

32.

Patterson

Wang

, et al. Prediction of residual strains due to in-plane fibre waviness in defective carbon-fibre reinforced polymers using ultrasound data. J Nondestruct Eval 2022; 42: 2. DOI: 10.1007/s10921-022-00910-0.

33.

Dvurecenska

Graham

Patelli

, et al. A probabilistic metric for the validation of computational models. R Soc Open Sci 2018; 5: 180687. DOI: 10.1098/rsos.180687.

34.

Schoeppner

Abrate

. Delamination threshold loads for low velocity impact on composite laminates. Compos Part A – Appl S 2000; 31(9): 903–915. DOI: 10.1016/S1359-835X(00)00061-0.

35.

Higuchi

Okabe

Nagashima

. Numerical simulation of progressive damage and failure in composite laminates using XFEM/CZM coupled approach. Compos Part A – Appl S 2017; 95: 197–207. DOI: 10.1016/j.compositesa.2016.12.026.

36.

Mirdehghan

. 1 - fibrous polymeric composites. In: Latifi

(ed). Engineered polymeric fibrous materials. Sawston: Woodhead Publishing, 2021, pp. 1–58.

37.

Krucinska

Stypka

. Direct measurement of the axial Poisson's ratio of single carbon fibres. Compos Sci Technol 1991; 41: 1–12. DOI: 10.1016/0266-3538(91)90049-U.

38.

MatWeb. Overview of materials for epoxy adhesive [dataset]. https://www.matweb.com/search/datasheettext.aspx?matguid=c1ec1ad603c74f628578663aaf44f261 (2024, accessed 28 September 2024).