Prediction of Damage in Non-Crimp Fabric Composites Subjected to Transverse Crushing: A Comparison of Two Constitutive Models

Abstract

Non-crimp fabrics (NCFs) are increasingly used in industry for manufacturing of composite structures due to a combination of high mechanical properties and excellent manufacturability. As with other composites, in-service damage can be a cause for severe reduction in load-carrying capacity of NCF-reinforced plastics. In this experimental and numerical study, two constitutive material models previously used only for damage prediction in unidirectional (UD) tape and woven fabric-reinforced materials (LS-DYNA’s *MAT_ENHANCED_COMPOSITE_DAMAGE—MAT54 and *MAT_LAMINATED_COMPOSITE_FABRIC—MAT58) were evaluated for simulating transverse crushing of composite parts processed from a non-crimp carbon fabric. For this purpose, UD NCF components of tubular shape were subjected to transverse crushing through a controlled indentation of a metallic cylinder, and results of the experiment were compared with numerical modeling. Considered verification metrics included the observed and the predicted patterns of interlaminar damage, the extent of delamination, as well as the ability of the models to replicate force-displacement response exhibited by the tested specimens.

1. Introduction

Non-crimp fabrics are increasingly used in industry for manufacturing of composite structures due to a combination of high mechanical properties and excellent manufacturability. Typical applications of non-crimp fabrics include load-bearing aerospace structures [1], automotive components [2], and wind energy applications [3].

NCFs consist of multiple straight and parallel yarn bands joined by polyester stitching (binder), as exemplified in Figure 1, for a unidirectional non-crimp carbon fiber fabric. On the mesoscale, such a structure usually produces distinct fiber bundles separated by resin-rich regions (unlike homogeneously distributed fibers in conventional UD tape composites). On the macrolevel, presence of the binder enhances manufacturability as compared to UD tape-based composites, while significantly lower crimp, as compared to woven fabric-based materials, has a positive effect on the mechanical properties of NCFs.

However, similarly to conventional unidirectional tape- and woven fabric-based materials (see, e.g., [4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19]), in-service damage, such as tool dropping or unintended collision, can be a cause for severe reduction in load-carrying capacity of NCF composites [20,21,22,23,24,25]. For example, Bibo et al. [24] demonstrated that resistance to propagation of damage and the compression after impact strength of NCF panels were not superior to their unidirectional prepreg tape equivalents. The latter reduced quickly with the increase of impact energy from the level of strength established for the intact panels.

Assessments of damage and residual load-bearing capacity of composite parts, regardless of the type of reinforcement, can be conducted using the integrated application of experimental destructive and non-destructive methods, which usually requires considerable financial and time investments. Understandably, the availability of a simulation approach suitable for this purpose would provide significant cost savings and accelerate the decision making when such assessments are required. While such processes are generally established for composites with conventional UD tape and fabric-based reinforcements, damage modeling in NCF structures using numerical simulation tools remains a challenge. Previous studies in this field considered detailed meso-scale simulations of non-crimp fabrics (see, e.g., [26,27,28]) and proposed new macroscale constitutive models for NCF materials [29,30,31,32]. However, such models are often not available to the general industry as part of commercial software, which limits their utilization at this point.

While development of novel macroscale constitutive relationships for NCFs is a promising approach in the long-term perspective, this study focuses on a short-term goal of evaluating existing material models available to the industry as part of commercial software and investigates their applicability to modeling of damage formation in structural elements made of NCF. In particular, two LS-DYNA material models were evaluated: *MAT_ENHANCED_COMPOSITE_DAMAGE (MAT 54) and *MAT_LAMINATED_COMPOSITE_FABRIC (MAT 58). Both constitutive models have been widely used for modeling of variety of in-service damage scenarios in UD tape and woven fabric composites, including damage induced by axial crushing (MAT 54: [5,12,33,34,35,36]; MAT 58: [12,37,38,39,40,41,42]) and transverse loading (MAT 54: [36,43,44]; MAT 58: [6,45,46,47]). The main takeaway from these studies is that both constitutive models can quite accurately predict the behavior of tape and woven fabric materials but require determination of their non-physical parameters (so-called stress limit factors, element erosion strain, etc.) through calibration with experimental data or, often, using the guidelines that can be deduced from the previous investigations.

To achieve the declared goals, the following approach was used in this study. First, the NCF material—a unidirectional non-crimp carbon fabric—was characterized using mechanical testing. Second, representative structural components—hollow CFRP tubes—were manufactured and subjected to transverse loading, in which damage was induced in the structural components through a controlled indentation of a metallic cylinder. Third, the physical experiment was replicated in LS-DYNA with MAT 54 and MAT 58 representations of the NCF composite. Finally, results of the physical experiments and numerical simulations were compared, and conclusions were made regarding the applicability and the best practices for using these constitutive models for simulating the behavior of non-crimp fabric materials.

2. Physical Experiments

2.1. Materials and Manufacturing

A unidirectional intermediate modulus (IM) non-crimp carbon fabric pre-impregnated with epoxy resin was used in this study. The fabric has the areal density of 140 g/m², the resin content of 38 ± 3%, and is built from 12 K carbon fiber tows stitched together with polyester binder (see Figure 1). As per the manufacturer’s recommendation, the curing cycle involved heating the material to 132 °C with the 4 °C/min ramp rate, holding at this temperature for 4 h and cooling to the room temperature with the same rate. The curing was conducted using a desktop oven with a programmable logic controller (PLC).

In addition to the temperature, fabrication of composite parts from the pre-impregnated non-crimp fabric implied application of pressure. For this purpose, and to obtain composite parts of the desired shapes, two aluminum molds were designed and manufactured: one to process flat composite panels (used to cut specimens from for material characterization) and another for manufacturing of the tubular components. A steel mandrel consisting of two wedge-shaped sections was used with the latter mold to form the interior of composite tubes. Pressure to composite laminates in both molds was applied by adjusting torque on the closing screws. The molds along with processed composite parts are shown in Figure 2. In-plane dimensions of the processed flat composite panels were 100 × 300 mm. The tubular specimens had the length of 200 mm and the dimensions of the internal cross-section of 25.4 mm × 25.4 mm. Thickness of a single ply in the laminates post-curing was measured to be equal to 0.15 mm.

The flat panels of different thickness were subsequently used to cut test specimens—with parameters listed in

Table 1. Parameters of the processed tests specimens for NCF mechanical characterization.

#	Test Type	In-Plane Dimensions, mm × mm	Number of NCF Layers	Test Procedure
1	Longitudinal tensile test	250 × 15	8	ASTM D3039 [48]
2	Transverse tensile test	175 × 25	14	ASTM D3039 [48]
3	10°off-axis tensile test	250 × 15	8	Ref [49]
4	Double-cantilever beam (DCB) test (specimen with Teflon insert)	125 × 25	32	ASTM 5528b [50]
5	End-notched flexure (ENF) test (specimen with Teflon insert)	120 × 25	32	Ref [51]

—for NCF characterization. Digital image correlation (DIC) speckling was applied to all tensile test specimens (groups #1–3 in Table 1).

Two groups of tubular components (see Figure 2) featuring different layups—${\left[0_4,{90}_3\right]}_{\mathrm{S}}$ or ${\left[{90}_5,0_2\right]}_{\mathrm{S}}$ (angles measured from the longitudinal axis of the tube)—were manufactured for this study.

2.2. Material Characterization

All material characterization tests have been conducted in this study using the electromechanical load frame MTS Criterion C43.504, equipped with the 50 kN LPS-504 load cell. Mechanical V-shaped wedge grips were affixed to the test pieces on both ends (see Figure 3). Strain measurements during the tests have been conducted using the Advantage AVX Video Extensometer Manta G146B ASG camera and digital image correlation software (MTS AVX Digital Extensometer software for all tensile tests except 10° off-axis experiments, where test video was post-processed using the GOM Correlate software).

Stress–strain diagrams obtained for the NCF material during the mechanical characterization campaign are shown in Figure 4, while the measured mechanical properties are summarized in

Table 2. Mechanical properties of the unidirectional NCF material.

Property	Units	Mean Value
Longitudinal Young’s modulus, $E_1$	MPa	149,018
Transverse Young’s modulus, $E_2$	MPa	6071
Major in-plane Poisson’s ratio, ${\nu }_{12}$	-	0.32
In-plane shear modulus, $G_{12}$	MPa	4217
Longitudinal tensile strength, $X_t$	MPa	2060
Longitudinal compressive strength *, $X_c$	MPa	1539
Transverse tensile strength, $Y_t$	MPa	29.1
Transverse compressive strength *, $Y_c$	MPa	126.6
In-plane shear strength, $S_L$	MPa	44.5
Longitudinal tensile strain-at-failure, ${\epsilon }_{1f}$	%	1.37
Transverse tensile strain-at-failure, ${\epsilon }_{2f}$	%	0.40
In-plane shear strain-at-failure, ${\gamma }_{12f}$	%	1.71
Mode I strain energy release rate, GIc	$\mathrm{kJ}/{\mathrm{m}}^2$	0.66
Mode II strain energy release rate, GIIc	$\mathrm{kJ}/{\mathrm{m}}^2$	2.77
Shear stress at onset of non-linearity	MPa	30
Shear strain at onset of non-linearity	%	0.711

* Note: this parameter was not measured experimentally in this study and was estimated using the Lo–Chim micromechanical model for the longitudinal compressive strength [52] or the Chamis equation for the transverse compressive strength [53].

. Along with the measured parameters, Table 2 includes compressive properties of the material, which were evaluated in this study using the available micromechanical models [52,53] rather than determined experimentally.

2.3. Component Testing

Physical experiments involving transverse crushing of the NCF specimens of the tubular cross-section were carried out in this study. In the experiment, central loading of the tube by a cylindrical indenter introduced local deformation and damage of the NCF material (fiber breakage, matrix cracking, and delamination), as well as a reduction of cross-section in the middle of the tube. Crosshead displacements and forces were recorded to compare these test outputs with predictions of numerical simulations. Loading cylinder(s) were allowed to move 12 mm downward with the speed of 0.5 mm/min.

The experiments were conducted with composite tubes of two different layups: ${\left[0_4,{90}_3\right]}_{\mathrm{S}}$ and ${\left[{90}_5,0_2\right]}_{\mathrm{S}}$ (0 and 90—ply orientation, subscripts—number of plies in the corresponding direction, “S” stands for “symmetry”). Force–displacement diagrams resulted from these tests are also shown in Figure 5. After crushing, a 100 mm central segment of the tube with ${\left[0_4,{90}_3\right]}_{\mathrm{S}}$ layup was subjected to X-ray computed tomography to provide the measure of the extent of induced damage phase and assist with verification of numerical models. Some of the CT scan outputs are exemplified in Figure 6, illustrating the multiplicity of failure modes experienced by the specimen.

3. Numerical Model

3.1. General Description of the Numerical Model

A numerical model that replicates the tubular component crushing test described in the previous section was developed in LS-DYNA. Its major features are depicted in Figure 7. In the model, the composite layup was represented using stacked TSHELL elements (same formulation as conventional shells with Reissner–Mindlin kinematics, but shaped as solid elements), each representing one “physical” ply of the NCF material. A uniform element size of 0.7 mm was used. The support and the loading cylinder were represented as rigid bodies.

As both investigated material models—MAT54 and MAT58—are typically used in modeling of dynamic processes using LS-DYNA’s explicit solver, the same solver was employed in simulations with these constitutive models in this study. To reduce the computational time when modeling the quasi-static processes using an explicit solver, the loading speed was increased from 0.5 mm/min as in the experiments to 2 m/s. To further reduce simulation time and take advantage of the geometrical and loading symmetry, only quarter of the NCF tube was modeled and symmetry boundary conditions were applied. The developed numerical models incorporated one of the two layups used in the physical experiments: ${\left[0_4,{90}_3\right]}_{\mathrm{S}}$ or ${\left[{90}_5,0_2\right]}_{\mathrm{S}}$.

3.2. Modeling of Intra-Ply Damage

3.2.1. MAT54—*MAT_ENHANCED_COMPOSITE_DAMAGE

Material failure in this model is governed by the Chang–Chang stress criterion, given as (see, e.g., [54]):

$${\mathrm{e}}_{\mathrm{fT}}^2={\left(\frac{{\mathsf{\sigma }}_{11}}{{\mathrm{X}}_{\mathrm{T}}}\right)}^2+\mathsf{\beta }{\left(\frac{{\mathsf{\tau }}_{12}}{{\mathrm{S}}_{\mathrm{L}}}\right)}^2-1$$

(1)

$${\mathrm{e}}_{\mathrm{fC}}^2={\left(\frac{{\mathsf{\sigma }}_{11}}{{\mathrm{X}}_{\mathrm{C}}}\right)}^2-1$$

(2)

$${\mathrm{e}}_{\mathrm{mT}}^2={\left(\frac{{\mathsf{\sigma }}_{22}}{{\mathrm{Y}}_{\mathrm{T}}}\right)}^2+{\left(\frac{{\mathsf{\tau }}_{12}}{{\mathrm{S}}_{\mathrm{L}}}\right)}^2-1$$

(3)

$${\mathrm{e}}_{\mathrm{mC}}^2={\left(\frac{{\mathsf{\sigma }}_{22}}{2{\mathrm{S}}_{\mathrm{L}}}\right)}^2+\left[{\left(\frac{{\mathrm{Y}}_{\mathrm{C}}}{2{\mathrm{S}}_{\mathrm{L}}}\right)}^2-1\right].\left(\frac{{\mathsf{\sigma }}_{22}}{{\mathrm{Y}}_{\mathrm{C}}}\right)+{\left(\frac{{\mathsf{\tau }}_{12}}{{\mathrm{S}}_{\mathrm{L}}}\right)}^2-1$$

(4)

where ${\mathsf{\sigma }}_{\mathrm{ij}}$ are in-plane stresses in a ply and ${\mathrm{e}}_{\mathrm{fT}}$, ${\mathrm{e}}_{\mathrm{fC}}$, ${\mathrm{e}}_{\mathrm{mT}}$, and ${\mathrm{e}}_{\mathrm{mC}}$ are the failure indices for longitudinal tensile, longitudinal compressive, transverse tensile, and transverse compressive failure, correspondingly. Parameter β (BETA) is defined in

Table 3. The non-physical parameters for MAT 54.

Parameter	Meaning	Units	Value	Comment for the Chosen Initial Value
DFAIL	Maximum strains for directional straining at which element will be eroded.	mm/mm	–	Disabled to control elements’ erosion by timestep (TFAIL) and effective strain (EPS) only.
TFAIL	Element is deleted when its time step is smaller than the given value.	s	10−7	Element is deleted when current timestep is less 10−7s
EPS	Effective failure strain	mm/mm	0.55	Chosen as to be significantly higher than any directional strain at failure initiation, as used in [12]
SOFT	Softening reduction factor for material strength in crash front elements	–	0.90	Slight (10%) reduction of strength in crashfront elements is assumed
SOFT2	Optional transverse softening reduction factor	–	no input	Softening is assumed to be isotropic
PFL	Percentage of layers which must fail until crash front is initiated	–	100	Default value
BETA	Weighting factor for shear term in tensile fiber mode (see Equation (1))	–	0	No effect of shear stresses on fiber tensile failure (max stress criterion), which usually provides good agreement with experimental data.
SLIMT1	Factor to determine the minimum stress limit after stress maximum (fiber tension).	–	0.010	Small but non-zero residual strength is assumed after tensile failure to avoid numerical instabilities
SLIMC1	Factor to determine the minimum stress limit after stress maximum (fiber tension).	–	0.375	A post-calibration value obtained in [12]
SLIMT2	Factor to determine the minimum stress limit after stress maximum (matrix tension).	–	0.010	Small but non-zero residual strength is assumed after tensile failure to avoid numerical instabilities
SLIMC2	Factor to determine the minimum stress limit after stress maximum (matrix compression).	–	0.375	A post-calibration value obtained in [12]
SLIMS	Factor to determine the minimum stress limit after stress maximum (shear).	–	1.000	A recommended value [54]
FBRT	Reduction factor for fiber tensile strength after matrix compressive failure	–	0	A zero effect of transverse matrix cracking on fiber tensile strength is assumed
YCFAC	Reduction factor for compressive fiber strength Xc after matrix compressive failure	–	12	Assume no influence of the matrix failure on the longitudinal compressive strength

. In addition, the material model takes into account a decrease of longitudinal compressive strength of a ply (XC) in case of transverse matrix failure, caused by reduction of matrix efficiency in supporting fibers against microbuckling. This is represented by a reduction factor YCFAC, such that the compressive strength in the fiber direction after compressive matrix failure is reduced to ${\mathrm{X}}_{\mathrm{C}}=\mathrm{YCFAC}·{\mathrm{Y}}_{\mathrm{C}}$.

In addition to parameters that can be and were determined experimentally (see Table 1), MAT54 has a set of non-physical input parameters, which can be categorized into 3 groups:

(1)

(element deletion) parameters;

(2)

parameters controlling crashfront softening;

(3)

those characterizing material behavior after failure initiation.

Parameters of the first group control automatic removal of heavily distorted elements from the model upon satisfaction of a timestep- or strain-based criteria. In this study, in order to reduce the number of parameters requiring calibration, only two out of six available criteria for erosion were used: one based on a minimal timestep (TFAIL), and another based on effective strain of an element (EPS). In the latter criterion, the effective strain that triggers element deletion was set to 55%, i.e., a value that significantly exceeds any of the directional failure strains of the composite. This ensures that erosion never precedes the complete physical failure of elements.

Parameters of the second group are intended to represent a damage zone that may develop in the material ahead of the crashfront (i.e., in elements that are adjacent to those that will be deleted next in the simulation) and can also be viewed as a simple numerical measure to avoid global buckling during simulation. As elements are removed from the simulation due to erosion, the elements which share nodes with them may have their strength reduced by some factor given by the SOFT parameter.

Post-failure response of composite is mainly governed by an array of stress limit factors (SLIM_), which represent the amount of residual strength that composite retains after satisfaction of any of the failure criteria provided by Equations (1)–(4). For example, even a completely crushed composite usually can retain some resistance to compressive loading. All of the non-physical parameters of MAT54, together with the rationale used for their initial choice, are listed in Table 3.

3.2.2. MAT58–*MAT_LAMINATED_COMPOSITE_FABRIC

MAT058 or *MAT_LAMINATED_COMPOSITE_FABRIC—a continuum damage mechanics (CDM)-based model, which accounts for both pre- and post-peak softening of composite plies—was used to model the intra-ply behavior of the NCF-reinforced composite tubes. Details of this model’s implementation can be found in [55]. The following set of failure criteria was used in this study to represent failure of the NCF material:

$${\mathrm{e}}_{\mathrm{fT}}^2={\left(\frac{{\overline{\mathsf{\sigma }}}_{11}}{{\mathrm{X}}_{\mathrm{T}}}\right)}^2-1$$

(5)

$${\mathrm{e}}_{\mathrm{fC}}^2={\left(\frac{{\overline{\mathsf{\sigma }}}_{11}}{{\mathrm{X}}_{\mathrm{C}}}\right)}^2-1$$

(6)

$${\mathrm{e}}_{\mathrm{mT}}^2={\left(\frac{{\overline{\mathsf{\sigma }}}_{22}}{{\mathrm{Y}}_{\mathrm{T}}}\right)}^2+{\left(\frac{{\overline{\mathsf{\tau }}}_{12}}{{\mathrm{S}}_{\mathrm{L}}}\right)}^2-1$$

(7)

$${\mathrm{e}}_{\mathrm{mC}}^2={\left(\frac{{\overline{\mathsf{\sigma }}}_{22}}{{\mathrm{Y}}_{\mathrm{C}}}\right)}^2+{\left(\frac{{\overline{\mathsf{\tau }}}_{12}}{{\mathrm{S}}_{\mathrm{L}}}\right)}^2-1$$

(8)

The effective stresses (${\overline{\mathsf{\sigma }}}_{\mathrm{ij}}$) in the expressions (5)–(8) are related to the nominal stresses through the damage parameters ${\mathrm{d}}_{\mathrm{ij}}$, also known as area loss parameters, such that:

$$\left[\overline{\mathsf{\sigma }}\right]=\left[\begin{array}{c}{\overline{\mathsf{\sigma }}}_{11}\\ {\overline{\mathsf{\sigma }}}_{22}\\ {\overline{\mathsf{\tau }}}_{12}\end{array}\right]=\left[\begin{array}{ccc}\frac{1}{1-{\mathrm{d}}_{11}}& 0& 0\\ 0& \frac{1}{1-{\mathrm{d}}_{22}}& 0\\ 0& 0& \frac{1}{1-{\mathrm{d}}_{12}}\end{array}\right].\left[\begin{array}{c}{\mathsf{\sigma }}_{11}\\ {\mathsf{\sigma }}_{22}\\ {\mathsf{\tau }}_{12}\end{array}\right]$$

(9)

where damage evolution with straining is assumed as ${\mathrm{d}}_{\mathrm{ij}}=1-\exp \left[-\frac{1}{\mathrm{me}}.{\left(\frac{ϵ}{{ϵ}_{\mathrm{f}}}\right)}^{\mathrm{m}}\right]$ and m, $ϵ,$ and ${ϵ}_{\mathrm{f}}$ are the parameter controlling shape of the stress–strain response, strain and strain at maximum directional stress, correspondingly. Thus, the components of the constitutive tensor C(d) can be represented as functions of the damage parameters and the properties of undamaged layer:

$$\mathrm{C}\left(\mathrm{d}\right)=\frac{1}{\mathrm{D}}\left[\begin{array}{ccc}\left(1-{\mathrm{d}}_{11}\right){\mathrm{E}}_{11}& \left(1-{\mathrm{d}}_{11}\right)\left(1-{\mathrm{d}}_{22}\right){\mathsf{\nu }}_{21}{\mathrm{E}}_{22}& 0\\ \left(1-{\mathrm{d}}_{11}\right)\left(1-{\mathrm{d}}_{22}\right){\mathsf{\nu }}_{12}{\mathrm{E}}_{11}& \left(1-{\mathrm{d}}_{22}\right){\mathrm{E}}_{22}& 0\\ 0& 0& \mathrm{D}\left(1-{\mathrm{d}}_{12}\right){\mathrm{G}}_{12}\end{array}\right]$$

(10)

where $\mathrm{D}=1-\left(1-{\mathrm{d}}_{11}\right)\left(1-{\mathrm{d}}_{22}\right){\mathsf{\nu }}_{12}{\mathsf{\nu }}_{21}>0$.

Such CDM-based formulation provides smooth increase of damage and, upon failure initiation, prevents the immediate drop of stresses to the level determined by the SLIM_ factors, as it is the case with MAT 54, but rather makes it a gradual reduction to that level. Moreover, it should be noted, that the two damage parameters ${\mathrm{d}}_{11}$ and ${\mathrm{d}}_{22}$ assume different values for tension (${\mathrm{d}}_{11+}$ and ${\mathrm{d}}_{22+}$) and compression (${\mathrm{d}}_{11-}$ and ${\mathrm{d}}_{22-}$). Additional non-physical parameters associated with MAT 58, as well as rationale for their choice, are described in

Table 4. The non-physical parameters for MAT 58.

Parameter	Meaning	Units	Value	Comment for the Chosen Initial Value
TFAIL	Timestep for automatic element deletion	s	10−7	Element is deleted when current timestep is less 10−7 s
ERODS	Maximum effective strain for element failure. If lower than zero, element fails when the effective strain calculated from the full strain tensor exceeds ERODS	mm/mm	−0.55	Chosen as to be significantly higher than any directional strain at failure initiation
SOFT	Softening reduction factor for material strength in crash front elements	–	0.90	Slight (10%) reduction of strength in crashfront elements is assumed
SLIMT1	Factor to determine the minimum stress limit after stress maximum (fiber tension)	–	0.010	Small but non-zero residual strength is assumed after tensile failure to avoid numerical instabilities
SLIMC1	Factor to determine the minimum stress limit after stress maximum (fiber tension).	–	0.375	A post-calibration value obtained in [12]
SLIMT2	Factor to determine the minimum stress limit after stress maximum (matrix tension)	–	0.010	Small but non-zero residual strength is assumed after tensile failure to avoid numerical instabilities
SLIMC2	Factor to determine the minimum stress limit after stress maximum (matrix compression)	–	0.375	A post-calibration value obtained in [12]
SLIMS	Factor to determine the minimum stress limit after stress maximum (shear)	–	1.000	A recommended value [54]

.

3.3. Modeling of Delamination

The walls of the tube were modeled using 14 layers of stacked TSHELL elements (TSHELL ELFORM 1 in LS-DYNA), each representing one physical ply of the UD NCF material. To represent delamination, the contact interfaces between the 14 TSHELL layers were added using *CONTACT_AUTOMATIC_ONE_WAY_SURFACE_TO_SURFACE_TIEBREAK with OPTION 9 [56] This contact algorithm is equivalent to using zero-thickness cohesive zone elements and is based on the fracture model with bilinear traction-separation law, mixed mode delamination criterion, and damage formulation [54]. This delamination model’s parameters include normal and shear failure stresses at the interface between the adjacent layers (NFLS and SFLS parameters, respectively), mode I and mode II critical strain energy release rates (${\mathrm{G}}_{\mathrm{Ic}}$ and ${\mathrm{G}}_{\mathrm{IIc}}$), and the normal (CN) and tangential (CT; CT = CT2CN × CN, where CT2CN is a coefficient between 0 and 1) stiffness of the material in the interlaminar region. For each mode of fracture (normal or shear), the ultimate separation displacement (δ_ult) is calculated based on the provided values of initial stiffness, peak stress, and fracture toughness. After satisfaction of the failure criterion ($\mathsf{\delta }>{\mathsf{\delta }}_{\mathrm{ult}}$), the failed interface (master segment–slave node pair) can only resist compressive force. The parameters of the delamination model used in this study, as well as the rationale for their choice, are summarized in

Table 5. Input data for *CONTACT_AUTOMATIC_ONE_WAY_SURFACE_TO_SURFACE_TIEBREAK.

Property	Unit	Value	Rationale
NFLS	MPa	75.00	The ultimate strength of bulk epoxy resin multiplied by a scaling factor of 0.95 to account for the mesh dependency (see the recommendation provided in [57] for meshes with element sizes of 0.7 mm).
SFLS	MPa	43.30	Assumed as $\mathrm{SFLS}=\mathrm{NFLS}\sqrt{3}$(von Mises criterion)
G_Ic	$\mathrm{kJ}/{\mathrm{m}}^2$	0.66	Measured experimentally, see Table 1
G_IIc	$\mathrm{kJ}/{\mathrm{m}}^2$	2.77	Measured experimentally, see Table 1
CN	MPa/mm	200,000	CN = Eepoxy/δRRR, where Eepoxy is the Young’s modulus of epoxy matrix (~3650 MPa) and δRRR is the thickness of the interlaminar resin-rich region (typically within 0.01 and 0.10 mm). Thus, the lower and upper bounds for CN correspond to 36,500 MPa/mm and 365,000 MPa/mm, accordingly. This averages to 200,000 MPa/mm as an estimate for the CN parameter.
CT2CN	-	0.37	CT2CN = CT/CN = Gepoxy/Eepoxy = 1/2 × (1 + νepoxy), where Gepoxy and νepoxy are the shear modulus and the Poisson’s ratio (~0.35) of epoxy resin, correspondingly.

.

4. Results and Discussion

Metrics to evaluate the numerical predictions specific to the crushing phase included comparison of the damage in the outer layer of the specimen (visual vs. predicted), extent of delamination (CT-scan vs. models), and the force-displacement diagrams (recorded in experiments vs. produced by LS-DYNA).

Figure 8 and Figure 9 show visual damage induced due to crushing to the outer layer of the [0₄, 90₃]_S specimen, as well as predictions of it obtained using MAT54 (Figure 10) and MAT58 (Figure 11). As can be seen in Figure 8, the outer layer exhibited fiber breakage (crack propagates under the crushing cylinder perpendicular to the fibers in the outer layer) and multiple matrix cracks (propagate along fibers), which are especially noticeable near the top and bottom surfaces of the specimen in the central region and also seen on the sides.

The predicted damage is depicted in terms of the history variables of the material models, which included longitudinal tensile and compressive and transverse tensile and compressive failure (see Equations (1)–(4)) for MAT54, as well as longitudinal, transverse, and shear damage for MAT58. It should be noted that MAT54 does not have a separate history variable for in-plane shear, while MAT58 does not distinguish between tensile and compressive modes (only cumulative longitudinal or transverse damage). As MAT54 predicts only instant failure (rather than gradual damage), any colors other than grey (not failed material) and red (failed material) appear in Figure 8 due to nodal averaging employed during post-processing.

As can be deduced from the figures, both material models predict fiber breakage observed experimentally in the central region of the specimen’s outer ply. The shapes of the predicted crack (true for both MAT54 and MAT58) is more simplistic (straight path) than that observed experimentally, which can be attributed to the quarter-symmetrical nature of the models. With MAT54, by the end of the crushing phase it is represented by the row of eroded elements, while erosion had not yet happened at that point in MAT58 model, and the crack can be identified by the elements fully damaged in the fiber direction (MAT58 history variable 1). Additional fiber damage that can be seen in Figure 9 (blue area) can be attributed to the initial accumulation of defects during the pre-peak softening phase in MAT58.

Both models seem to be able to replicate the qualitative pattern of matrix cracking observed experimentally and described above (history variables 3 and 4 in Figure 8, and 2 and 3 in Figure 9). This includes predicted high density of matrix cracks in the upper and lower areas of the central region of the specimen, as well as multiplicity of scattered matrix cracks at the sides (MAT54 history variable 3 and MAT58 history variable 2).

The model of the crushed tube with [0₄, 90₃]_S layup obtained using X-ray computed tomography was processed using the Porosity Analysis module of myVGL software to visualize the region of the tube affected by delamination. In particular, the “Equivalent diameter” feature of the module was used, which indicates the diameter of a sphere that has the same volume as the detected defect in the material. In Figure 10a, delamination is represented by pink color, which highlights defects with the equivalent diameter of 3 mm and larger.

In modeling, delamination was visualized as contact gap between stacked TSHELL elements connected using *CONTACT_AUTOMATIC_ONE_WAY_SURFACE_TO_SURFACE_TIEBREAK (OPTION 9), as described in Section 3.3. In total there are 13 ply interfaces in each model on which delamination (contact gap) can be displayed. To make the modeling results comparable with the X-ray output, the contact gap images from ply interfaces were made half-transparent and placed on top of each other. Comparison of the regions of the crushed specimen affected by delamination—as measured experimentally and predicted—is shown in Figure 10. It should be noted that unlike the equivalent diameter feature in the CT-scan processing software, contact gap in LS-DYNA is a unitless measure of delamination. Pink color in Figure 10b,c represents completely debonded regions between plies (contact gap = 1). As can be deduced from the figure, all three images show similar extent of delamination. Conducted measurements indicated that the area affected by delamination constitutes 26.6% of the image (a) in Figure 10, 26.9% of the image (b), and 25.8% of the image (c).

The force–displacement diagrams corresponding to the crushing phase of the specimens with [0₄, 90₃]_S and [90₅, 0₂]_S layups are shown in Figure 11, comparing the experimental results and the modeling predictions. As can be deduced from the figure, both constitutive models equally adequate predict the force–displacement response of the specimen during this phase. This is true for both considered composite layups.

5. Conclusions

This numerical and experimental study was conducted to investigate the applicability of two constitutive material models previously used only for UD tape and woven fabric-reinforced materials to modeling of damage in composite parts processed from non-crimp carbon fabric. Despite the differences in their formulations (strength criteria-based MAT54 vs. CDM-based MAT58), both material models were found to be robust and able to predict damage induced by transverse crushing in UD NCF parts reasonably well. This included the representation of qualitative patterns of interlaminar damage and the extent of intralaminar damage (delamination), as well as the prediction of force-displacement response exhibited by the tested tubular components. Future work must be focused on investigating the applicability of MAT54 and MAT58 for prediction of residual load-bearing capacity of damaged NCF parts under transverse loading, as well as consideration of other loading scenarios, such as axial crushing.