Numerical Simulation of Composite Structural Failure: A Review of Methods and Correlation with Destructive Tests

Early design of composite structures often employed linear elastic analysis up to first-ply failure using simple failure criteria. Classical criteria such as the maximum stress/strain criterion, the Tsai-Hill criterion, and the Tsai-Wu polynomial criterion were widely used to predict whether a given combination of stresses in a ply would cause failure [28]. Tsai-Wu theory (dating back to 1971) provides a quadratic interaction equation for anisotropic strength, and can predict initial failure load reasonably for many laminates. However, these criteria are essentially one-shot failure indicators – they do not simulate the gradual damage accumulation beyond the first crack. In reality, composite laminates often have reserve capacity after initial ply failure; they may redistribute load and continue to carry increasing load until final collapse. Linear criteria cannot capture this post-failure behavior.

To more accurately simulate composite structural failure, researchers developed progressive damage models (PDMs) that include material degradation after the onset of damage. A seminal approach was proposed by Hashin, who introduced physically-motivated failure criteria distinguishing fiber and matrix failures in tension or compression [29]. Hashin’s 1980 criteria define separate modes (fiber tension, fiber compression, matrix tension, matrix compression) with associated stress-based limits. For example, in one common formulation of Hashin’s theory:

Fiber tension failure: Fft=(σ11XT)2+(τ12S)2+(τ13S)2≥1. ${F_{ft}} = {\left( {{{{\sigma _{11}}} \over {{X_T}}}} \right)^2} + {\left( {{{{\tau _{12}}} \over S}} \right)^2} + {\left( {{{{\tau _{13}}} \over S}} \right)^2} \ge 1.$

Fiber compression failure: Ffc=| σ11 |XC≥1. ${F_{fc}} = {{\left| {{\sigma _{11}}} \right|} \over {{X_C}}} \ge 1.$

Matrix tension failure: Fmt=(σ22YT)2+(τ12S)2+(τ13S)2≥1. ${F_{mt}} = {\left( {{{{\sigma _{22}}} \over {{Y_T}}}} \right)^2} + {\left( {{{{\tau _{12}}} \over S}} \right)^2} + {\left( {{{{\tau _{13}}} \over S}} \right)^2} \ge 1.$

Matrix compression failure: Fmc=−σ22YC+(τ12S)2≥1, for σ22<0. ${F_{mc}} = {{ - {\sigma _{22}}} \over {{Y_C}}} + {\left( {{{{\tau _{12}}} \over S}} \right)^2} \ge 1,{\rm{\;for\;}}{\sigma _{22}} < 0.$

where

Here X_T, X_C are the longitudinal (fiber direction) tensile and compressive strengths, Y_T, Y_C are the transverse (matrix direction) strengths, and S is the in-plane shear strength. Failure occurs when the relevant condition is met (failure index ≥ 1). In a progressive analysis, once a failure criterion is met in an element (or an integration point of a shell/solid representing a ply), the stiffness of that element is reduced to simulate the loss of load-carrying capacity. Often a sudden drop to near-zero stiffness (complete “element death”) is used for fiber failure, whereas matrix cracking might be modeled by a more gradual stiffness reduction. This approach falls under continuum damage mechanics, where internal damage variables degrade the stress-strain response (e.g. E_eff =(1- D) E). Modern finite element codes allow user-defined material subroutines to implement such degradation models or provide built-in composite damage models. For example, Abaqus and ANSYS offer Hashin-type progressive failure models that automatically reduce ply stiffness after the criterion is exceeded. Yeh and colleagues extended such criteria (e.g., Yeh-Stratton) to refine predictions, and many variations (e.g., LaRC04/05 by NASA) have been proposed to improve accuracy for specific stress states.

An advanced alternative is the Puck inter-fiber failure theory, a physically based model focusing on matrix-dominated failure modes [30]. Puck’s criterion considers matrix cracking under combined normal and shear stresses and predicts the fracture plane orientation based on a Mohr-Coulomb-type concept. In practice, Puck’s criterion tends to give more accurate predictions for matrix compression failures (which often occur on inclined planes as a form of shear fracture). Implementation in finite element models is more involved due to the need to track the fracture angle, yet researchers have demonstrated improved correlation with experiments in compression and shear-loaded laminates.

For fiber-dominated failures (tensile or compressive fiber breakage), Puck’s and Hashin’s criteria are comparable, as both essentially reduce to a maximum strain or stress in the fiber direction. Comparative studies (including the World-Wide Failure Exercise) have found that most modern criteria can predict fiber failure with similar accuracy, but differ significantly in matrix failure predictions [31]. As an example, Puck’s failure criterion for a UD composite lamina under transverse compression can be expressed in a simplified form as:

Fmc=−σ22YC+(τ12S)2≥1, for σ22<0. {F_{mc}} = {{ - {\sigma _{22}}} \over {{Y_C}}} + {\left( {{{{\tau _{12}}} \over S}} \right)^2} \ge 1,{\rm{for}}{\sigma _{22}} < 0.

for σ₂₂ < 0. Here -σ₂₂/Y_C is a term that linearly accounts for increased matrix strength under compression (analogous to a frictional cohesion), and S_L is the longitudinal shear strength of the lamina. Failure is predicted when the combination of normal and shear stresses on the fracture plane reaches this condition. More complex formulas exist for other loading cases (including matrix tension and the calculation of the fracture plane angle), but the above gives a flavor of Puck’s approach.

Progressive damage models require more parameters than linear criteria, notably in defining stiffness reduction after failure initiation – i.e., the damage evolution law. Some models use a sudden drop to zero (brittle failure), while others apply gradual softening characterized by fracture energy to avoid mesh sensitivity. Zobeiry et al. [32] characterized strain-softening behavior in composites under tension and compression, providing data to calibrate damage evolution laws. The material shows a softening response where stress decreases with increasing strain after the peak, due to microcrack accumulation. This is often implemented via a damage variable D that increases from 0 (undamaged) to 1 (fully failed), reducing the effective stiffness as E_eff =(1- D) E. Damage evolution is calibrated by matching stress-strain curves from experiments (including the post-peak region) or by using fracture toughness values.

For fiber breakage, a brittle assumption (immediate D = 1) is usually sufficient. For matrix cracking, a gradual degradation better represents distributed damage. More elaborate continuum damage models (e.g., by Ladevèze) incorporate such concepts, though they are beyond the scope of this review. In summary, while linear criteria offer quick estimates and are useful in early design, progressive damage models enable more accurate simulation of damage development and ultimate failure – particularly when correlated with physical test results.

Delamination is one of the most critical failure mechanisms in composite laminates, often initiated under transverse loading, impact, or during fatigue. It results in the separation between plies and leads to a loss of structural integrity and stiffness. Unlike fiber and matrix failure, delamination occurs between layers and is governed by interlaminar stresses, particularly Mode I (opening), Mode II (sliding), and mixed-mode fracture. To model delamination, two main strategies are commonly used in finite element analysis: the Virtual Crack Closure Technique (VCCT) and Cohesive Zone Modeling (CZM). VCCT is based on linear elastic fracture mechanics and uses calculated energy release rates to simulate crack propagation. While effective in modeling the growth of pre-existing delaminations, VCCT requires predefined crack paths and is less suited to spontaneous crack initiation. This makes it more appropriate for problems where the crack front is known a priori or derived from tests. Cohesive Zone Modeling has become the preferred approach for simulating delamination initiation and growth. In this method, interface elements (cohesive elements) are placed between plies and assigned a traction-separation law. Initially, these elements behave elastically. Once the interfacial strength is exceeded, damage initiates, and as the separation increases, traction reduces, eventually reaching zero when the critical fracture energy is dissipated. CZM offers several advantages: it can handle arbitrary delamination paths, initiate cracks without a pre-existing flaw, and simulate mixed-mode behavior.

CZM’s application has been widespread in both static and dynamic simulations, including low-velocity impact (LVI) and compression-after-impact (CAI) tests. For example, Lemanski et al. [33] used cohesive elements to simulate delamination in specimens with fiber waviness defects under compressive loading. Their model captured the interaction between local buckling and delamination propagation, enabling accurate prediction of failure loads.

Similarly, Calvo et al. [34] demonstrated that digital image correlation (DIC) data could validate CZM predictions, with good agreement in delamination patterns and strain localization. Modern implementations of CZM often use bilinear traction-separation laws, with parameters calibrated from Double Cantilever Beam (DCB) and End-Notch Flexure (ENF) tests. The cohesive strength and fracture toughness in Mode I (G_Ic), Mode II (G_IIc), and the mixed-mode criterion (e.g., Benzeggagh-Kenane) are key inputs. When these parameters are accurately defined, CZM can reproduce delamination onset and growth with high fidelity. Challenges in cohesive modeling include computational cost (due to the added elements between each ply) and mesh sensitivity (requiring fine mesh near interfaces). Nonetheless, CZM has been successfully applied in commercial FEA codes such as Abaqus (COH3D8, COH2D4), LS-DYNA, and ANSYS. Researchers have used adaptive meshing and submodeling to manage computational demands.

Fig. 3.

Experimental photographs and finite element predictions of failed specimens with failure locations highlighted (adopted from [33]).

Fig. 4.

Comparison of experimental and numerical load-displacement analysis, showing delamination pattern (adopted from [34]).

In summary, CZM is a powerful and versatile tool for modeling interlaminar failure in composites. It provides predictive capabilities for delamination initiation and progression and is a crucial component in simulating impact events, CAI behavior, and bonded joint failures.

Simulating the failure of composite structures requires careful selection of finite element (FE) modeling parameters. These include element type, mesh density, time integration scheme, and how damage initiation and evolution are numerically handled. For structural-scale models, shell or layered shell elements (e.g., S4R in Abaqus) are often used due to their efficiency and ability to represent laminate layups. Shells provide good accuracy for in-plane stresses but are limited in representing through-thickness effects like delamination unless combined with cohesive elements or enrichment. For higher fidelity, solid elements (e.g., C3D8 or SC8R) are preferred, especially when delamination, local buckling, or impact responses are critical. However, solid modeling increases the element count significantly and requires thin elements through the thickness to capture stress gradients accurately.

Mesh density is another critical parameter. Coarse meshes may delay or suppress damage initiation, while fine meshes capture localized stresses but increase computational cost. Progressive damage models are particularly sensitive to mesh size due to strain localization. To mitigate mesh sensitivity, researchers apply fracture energy-based softening laws with regularization techniques. Cohesive elements used for delamination (e.g., COH3D8 or COH2D4 in Abaqus) must be placed between each ply. Their effectiveness depends on correct calibration of traction-separation laws and sufficient mesh resolution along the interface. Mesh refinement around high-gradient zones (e.g., bolt holes, notches) is often necessary.

Regarding time integration, both implicit and explicit solvers have their use cases. Implicit solvers (e.g., Abaqus/Standard) are suited for quasi-static problems but may struggle with convergence after damage initiates. Nonlinear behavior, softening, and snap-back responses are challenging to capture. Techniques like arc-length control, damping, or viscous regularization are employed to stabilize simulations. Explicit solvers (e.g., Abaqus/Explicit or LS-DYNA) are more robust in handling highly nonlinear, transient, or progressive failure problems, including crash and impact. They are increasingly used for quasi-static simulations by artificially reducing loading rates while keeping kinetic energy low. Explicit analysis avoids convergence issues and can track large damage progression, although it requires small time steps for stability.

Contact modeling is another key aspect. Accurate definition of contact between bolts and laminate (with friction, penalty stiffness, and potential bolt preload) influences stress redistribution and damage location. Models using surface-to-surface contact with friction coefficients calibrated from experiments offer better fidelity than simplified tied interactions.

Boundary conditions should realistically represent the test setup. Symmetry, loading, constraints, and preloads must match the physical test to ensure correlation. Global-local modeling is sometimes used to focus mesh density only in regions of interest (e.g., around holes or impact zones), while the rest of the structure uses coarser models.

Calibration of material degradation parameters (stiffness reduction, fracture energy, interfacial strengths) is usually done using coupon tests (tension, compression, shear, DCB, ENF). These properties are input into the PDM and CZM models and validated at sub-component level before being applied to full-scale simulations.

In summary, accurate FE modeling of composite failure requires a well-balanced approach between model detail, computational efficiency, and physical realism. Selection of solver, element type, contact modeling, and damage criteria – all directly affect correlation with experimental data. With appropriate choices and calibration, FE tools can simulate failure progression, mode interaction, and even post-failure behavior reliably.

Bolted joints in carbon/epoxy laminates remain one of the most critical and experimentally validated structural features in aerospace composite design. Numerous studies have focused on understanding and predicting their failure behavior under tensile loads using both destructive testing and numerical modeling approaches. Xiao and Ishikawa [35, 36] carried out a foundational two-part study analyzing single-bolt composite joints, where Part I detailed experimental failure modes (bearing, net-tension, shear-out), and Part II developed finite element models that closely replicated load-displacement responses and failure progression. Liang et al. [37] extended this by reinforcing bolted zones with local CFRP insert patches and demonstrated that progressive damage modeling via user-defined material subroutines (VUMAT in Abaqus) provided excellent correlation with experimental strength and failure locations.

Recent studies have also explored multi-bolt configurations. Yoon et al. [38] evaluated the effect of edge distance, bolt torque, and laminate thickness on bearing strength in single-lap joints and confirmed that proper modeling of clamping and contact can predict transitions between bearing and net-tension failure within 10-15% accuracy. Yang et al. [39] investigated multi-bolt composite-to-metal joints with countersunk fasteners, using both full-field strain measurements and nonlinear finite element simulations. Their results confirmed that numerical models could accurately capture the progressive load redistribution and sequential bolt failure. El-Sisi et al. [36, 37] presented experimental and numerical analyses of single- and double-lap CFRP joints with multiple fasteners, showing that woven laminates tend to absorb more strain energy before final failure, while unidirectional laminates fail more abruptly. Finally, Belardi et al. [42] developed a specialized Composite Bolted Joint Element (CBJE) to efficiently simulate multi-bolt assemblies in large structures. Their approach, validated against physical tests of 3-bolt specimens, reproduced load-sharing behavior and out-of-plane displacements while significantly reducing computational time compared to full-detail FE models.

Fig. 5.

omparison of experimental X-ray radiography and numerical damage progression (adopted from [36]).

Fig. 6.

Comparison of experimental and numerical results for progressive damage in a bolted joint (adopted from [37]).

Fig. 7.

Failure modes of composite bolted joints: (a–d) geometrical variables and schematic failure modes, (e) experimental bearing damage at varying e/D ratios, and (f) numerical strain fields under tension and compression (adopted from [38]).

Taken together, these studies illustrate that with appropriate calibration, modern finite element methods can not only match experimental strength values but also replicate detailed damage sequences in bolted composite joints. Table 1 summarizes a few representative results from the literature on single-lap bolted joint tests and simulations.

Shear tear-out and pull-through are critical failure modes in bolted composite joints, particularly where edge distances are insufficient or through-thickness stresses dominate. Both have been the subject of extensive experimental and numerical investigation. These mechanisms are of special concern in thin CFRP flanges subjected to high bearing loads near free edges, or under out-of-plane loading conditions. Shear tear-out occurs when material behind the bolt head fails in shear due to localized overstress. It often manifests as a block of laminate being pulled from the edge, with crack propagation along the 0° fiber direction and delamination planes. Pull-through, on the other hand, is typically observed when the bolt head is loaded through the thickness of the laminate, and the surrounding plies are crushed and separated in a conical or pyramidal pattern. Gupta et al. [43] analyzed pin-loaded composite plates, finding that a well-calibrated progressive damage model could reproduce the characteristic bearing failure progression observed in tests even in the absence of bolt clamp-up – essentially establishing a baseline case of pure bearing failure correlation. Cooper and Turvey [44] conducted early experiments showing how decreased edge-to-diameter ratios (e/D < 2) lead to shear tear-out rather than bearing failure.

Later work by Egan et al. [45] revealed that multi-bolt joint configurations could also experience localized pull-through in flanges with insufficient reinforcement or when subjected to prying forces. Advancements in numerical modeling, particularly finite element methods (FEM), have significantly enhanced predictive capabilities regarding shear tear-out. Camanho and Matthews introduced progressive damage models incorporating intralaminar and interlaminar failure criteria, successfully correlating numerical predictions with experimental findings [46]. Shipsha and Burman [47] performed detailed static and post-damage tests on composite structures with bolted joints, identifying interaction between shear-out and delamination. Their numerical simulations using cohesive interface elements accurately predicted tear-out initiation and propagation. Such simulations increasingly rely on cohesive elements to accurately represent delamination, crucially affecting tear-out predictions [48]. Catalanotti et al. [49] and Fu et al. [50] and Kelly and Hallström [51] focused on pull-through testing using quasi-static tests on CFRP panels with countersunk bolts. Their studies showed that progressive damage models, when calibrated with experimental strength and toughness values, can predict damage initiation near the bolt and through-thickness delamination leading to catastrophic pull-through.

Key to accurate modeling was the implementation of a combined intralaminar/interlaminar damage framework, capturing both matrix cracking and interfacial failure. Pearce et al. [52] and Perogamvros et al. [53] extended this to high-rate loading, using stacked-shell models and rate-sensitive cohesive laws to reflect dynamic pull-through behavior. Experimental validation confirmed that increased loading rate tends to raise pull-through resistance slightly due to strain rate sensitivity of the matrix. Zhang et al. [54] examined four-point bending of multi-bolt CFRP flanges, capturing partial pull-through and associated delamination in a C-beam geometry. The numerical results showed high correlation with DIC data and microscopy, reinforcing the value of combined experimental-numerical investigation. Montagne et al. [55] conducted comprehensive numerical and experimental analyses of single-lap bolted joints made from carbon-epoxy laminates, demonstrating strong agreement between numerical predictions and experimental results for stiffness and failure loads. However, differences were noted in accurately capturing transitions between shear-out and bearing failures based on geometric parameters. Similarly, Liu et al. [56] developed a progressive failure model utilizing continuum damage mechanics that closely matched experimental load-displacement curves for single-bolt composite joints. Their findings highlighted the model’s capability to accurately predict initiation and progression of damage, including matrix cracking and delamination, critical for joint strength predictions. Further detailed finite element investigations by Stocchi et al. [57] and Giannopoulos et al. [58] highlighted the critical role of bolt torque and preload effects on composite joint strength and damage evolution, reinforcing the necessity to incorporate realistic boundary conditions and fastener interactions in numerical simulations.

In summary, shear tear-out and pull-through failures can be modeled successfully using progressive failure analysis and cohesive zone modeling. Accurate prediction requires a realistic representation of bolt contact, edge geometry, laminate layup, and interlaminar properties. These modes are highly sensitive to bolt preload, edge margins, and localized laminate damage – all of which must be properly characterized in the model. (No separate table is provided for Section 3.2, as many tear-out and pull-through studies are also covered under bolted joint results in Table 1, and their accuracy falls in similar ranges.)

Open-hole tension (OHT) and bearing failure modes are classic scenarios used to assess the tensile performance and damage tolerance of composite laminates. These failure modes are commonly observed in mechanical joints and stress concentrators and are among the best-studied cases in the literature. In OHT specimens, a central hole introduces stress concentrations that initiate damage. The dominant failure mode depends on laminate layup, hole size, width-to-diameter ratio (w/D), and specimen geometry. Common failure mechanisms include net-tension (fracture across the ligament), shear-out (at low edge distances), and bearing (localized compression damage around the hole). Linear failure criteria such as Tsai-Wu and Hashin can predict the location of first-ply failure around the hole, but often underestimate ultimate strength. This is because after the first matrix cracks, load is redistributed to adjacent plies and damage can grow progressively. Progressive damage modeling (PDM), combined with cohesive elements at the hole interface, better captures this evolution. Yoon et al. [34] performed experiments and FE simulations on single-shear bolted joints in CFRP laminates, showing that for specimens with sufficient edge distance, bearing failure dominated and was gradual. For specimens with lower interlaminar strength, splitting and net-tension occurred instead. Their numerical model combined Hashin-based failure criteria with cohesive elements near the hole and predicted the observed failure mode transitions well.

Cameron et al. [59] compared thick-ply and thin-ply laminates in bearing tests and found that thin-ply laminates exhibited more gradual, less catastrophic failure due to improved transverse toughness and delayed splitting. Their simulations used continuum damage mechanics and matched experimental strength values and crack patterns within ~5–10%. Pin-loaded bearing-bypass tests further challenge modeling, as they combine bearing contact with global tension. Dogan et al. [60] reviewed multiple criteria and found that Puck’s theory better captured the inclined shear fracture paths than traditional maximum stress or Hashin models. Finite element simulations incorporating cohesive layers and frictional contact also improved correlation with tests.

In summary, OHT and bearing failures are well-suited to progressive damage modeling. When calibrated with experimental material data, numerical models can capture the sequence of damage events – matrix cracking, fiber kinking, delamination, and final rupture – and reproduce load-displacement curves, failure loads, and crack morphology with high accuracy. Accurate representation of contact, hole clearance, stacking sequence, and boundary conditions is essential to match experimental results.

Compression testing of composite laminates is one of the most challenging to simulate due to complex nonlinear failure mechanisms such as fiber microbuckling, delamination, and kink-band formation. These issues are amplified when impact damage is present, which can initiate premature buckling and catastrophic failure. Compression After Impact (CAI) tests follow a standard sequence: a composite panel is impacted (usually at low velocity) and then compressed to failure. This scenario is representative of real-world threats such as tool drops, hail, or bird strikes followed by operational loading. The CAI test exposes weaknesses in interlaminar toughness, matrix cracking resistance, and ply waviness tolerance. Early simulations of CAI used simplified shell models and linear failure criteria, which could predict stiffness loss but not accurate failure loads. Modern CAI modeling incorporates 3D solid elements, progressive damage models, and cohesive elements to simulate delamination. Lu et al. [61] examined CAI behavior in carbon/PEKK thermoplastic composites and observed higher residual strength compared to conventional thermoset laminates, attributed to the thermoplastic’s superior ability to arrest delamination growth. Their progressive damage simulation, incorporating the tougher interlaminar properties, predicted the post-impact compressive strength within 5% of the experimental value, correctly reflecting the enhanced damage tolerance of the thermoplastic system. Lemanski et al. [33] used detailed 3D models with fiber waviness defects and cohesive interfaces to simulate CAI behavior. They successfully captured the interaction between matrix damage, local buckling, delamination growth, and final collapse. Their simulations showed strong agreement with DIC strain maps and test results. Bull et al. [62] introduced image-based modeling by importing CT scan data of real impacted panels to reproduce fiber misalignment and initial delaminations. Their model predicted the CAI strength within 5% of experiments and matched failure patterns in post-mortem analysis. Yang et al. [63] and Zhang et al. [64] modeled tubular and sandwich structures under impact and compression. Their results emphasized the importance of modeling skin-core debonding, stress wave interactions, and local crushing. Recent benchmark studies show that with appropriate calibration, numerical CAI models predict residual strength within 5–15% of test results. Guo et al. [65] performed a parametric study of CAI simulations, showing that factors like layup sequence and interlaminar toughness significantly affect predicted residual strength. They reported that only by adjusting input parameters such as Mode I/II fracture energies for each layup could their FE models consistently match the CAI strength and damage extent observed in tests – underlining the importance of proper calibration for CAI correlation Parameters such as interlaminar fracture toughness (G_Ic, G_IIc), matrix compression strength, and ply waviness distribution play dominant roles. Explicit solvers are widely adopted for CAI because they can capture sudden failure without convergence issues. Submodeling techniques also help reduce computational cost by refining the mesh only near the impact site.

In summary, CAI modeling has matured significantly. The use of cohesive elements, CT-informed defects, and calibrated damage models allows accurate prediction of failure modes and post-impact residual strength. CAI serves as a crucial benchmark for validating advanced FEM strategies. Interestingly, the use of thin-ply laminates has been shown to improve damage resistance and residual strength in impact scenarios. Sebaey and Mahdi [66] demonstrated that spreading plies into thinner layers delays matrix cracking and delamination, resulting in higher CAI strength experimentally. Cugnoni et al. [67] further confirmed that these thin-ply composites exhibit significantly smaller damage areas for a given impact energy. To accurately simulate such behavior, failure criteria and damage models may require modification to account for ply-thickness effects, as standard models calibrated on conventional plies might over-predict damage in thin-ply systems. This is an active research area, bridging manufacturing advances (thin-ply technology) with simulation needs. Table 3 provides a few representative results for CAI simulations:

Repairing damaged composite components is essential in extending the service life of aerospace structures. Composite repairs restore strength to regions affected by manufacturing defects, impact, or fatigue. Modeling repair effectiveness is crucial to quantify residual strength and support certification decisions. Flush bonded repairs, such as scarf and step-lap patches, are common in aerospace. These methods aim to restore load transfer while minimizing stress concentrations. Orsatelli et al. [23] reviewed modeling strategies for such repairs, showing that cohesive zone modeling combined with progressive damage in the parent laminate accurately predicts residual strength. Baker and Bitton [24] demonstrated that scarf angle strongly influences the load-carrying capacity: shallow angles (1:50 or 1.2°) spread load more evenly and reduce interfacial shear. Their finite element model predicted failure away from the patch interface when the scarf geometry was optimal. Feraboli et al.[25] simulated repairs with step-lap adhesive bonding using cohesive elements to model the adhesive failure. Their model reproduced the experimental load-displacement curve and captured both adhesive failure and parent laminate damage. Duong et al. [26-27] tested panels repaired using resin injection in regions with delaminations. Their finite element simulations included reduced stiffness and altered properties in the injected region, matching experimental failure loads within 10%. The simulations also revealed that drilled holes used for resin access caused minimal stress concentrations. In all these studies, digital image correlation (DIC) and ultrasonic NDT were used to validate model predictions. Numerical models matched strain fields, damage progression, and failure locations observed in tests.

Repairs are also compared in simulation: bonded patches provide better performance than bolted patches for most static load cases. Bonded joints distribute load more evenly and eliminate stress concentrations around bolt holes. Simulations of bolted patch repairs often predict early bearing failure or interfacial delamination.

Challenges in modeling repairs include accurately representing the interface condition (e.g., residual adhesive thickness, surface preparation), predicting long-term durability (moisture, fatigue), and simulating damage initiation within the patch. Advanced simulations now consider durability effects – for example, fatigue simulations degrade interfacial strength over cycles. Environmental simulations adjust matrix stiffness and fracture toughness based on temperature and moisture exposure.

In summary, modeling composite repairs with FEM and CZM enables accurate assessment of strength recovery. Proper calibration and validation ensure that analytical substantiation can be used instead of extensive retesting. This is critical in aerospace certification, especially when defects are found during manufacturing or operation.

Blier et al. [68] carried out tests on hybrid bonded–bolted joints under tension and reported that the hybrid configuration improved joint strength and altered the failure mode compared to purely bolted joints. They developed an FE model including both cohesive elements for the adhesive and contact/friction for the bolt, which successfully predicted the joint’s load–displacement response and the sequence of adhesive vs. bolt failure observed in experiments. Mehrabian et al. [69] extended simulation studies to optimize hybrid joint configurations. They used FEM calibrated against experimental data to determine how varying the overlap length and bolt torque affected failure load. The simulation outcomes, verified by their experiments, helped identify an optimal balance where the adhesive carries load up to a certain point before bolts engage – achieving higher strength than either bonding or bolting alone.

Table 4 summarizes two representative outcomes in repair modeling.

Despite major progress in simulation methods, modeling composite failure remains inherently complex. One major issue is the natural variability of composite materials. Differences in fiber alignment, matrix content, voids, and ply thickness lead to unpredictable local behavior. This makes it challenging for numerical models to achieve exact correspondence with experimental results unless all such variations are explicitly accounted for – which is rarely feasible.

Mesh sensitivity is another persistent challenge. Progressive damage models (PDM), especially those involving element deletion or strain-softening, can yield mesh-dependent results. If the mesh is too coarse, damage initiation may be delayed or smeared; if too fine, numerical instability or unrealistic localization may occur. Mesh regularization or energy-based softening laws are often used to mitigate these issues.

Material property calibration is a third limitation. Most advanced models require not only elastic and strength parameters but also damage evolution data – fracture energies in multiple modes, degradation rates, and interaction laws. These are often difficult to measure and vary significantly between manufacturing batches or processing methods.

Inconsistent boundary conditions between tests and models are another source of error. Test fixtures may introduce unintended constraints, friction, or misalignment, which are not captured in idealized models. As a result, even well-tuned models can produce accurate results in some tests and poor ones in others. Furthermore, models often struggle to predict failure mode transitions. A laminate might fail in net-tension, shear-out, or bearing depending on subtle geometric or material differences, and many models are biased toward one outcome unless explicitly tuned.

In structural-scale models, error accumulation becomes a concern. Effects such as local wrinkling, fiber waviness, or ply gaps – which may be negligible in small coupons – can shift the failure load or path in full-scale parts.

Despite these limitations, the best-calibrated models typically achieve prediction accuracy within 5–15% for strength and can qualitatively reproduce failure modes, making them useful tools for design and certification. Continued refinement and validation against physical tests remain essential.

The selection of a failure criterion significantly affects the outcome of composite simulations. The most common criteria include Hashin, Puck, and NASA’s LaRC04/05 models. Hashin’s criterion (1980) is widely implemented in commercial FEA software and distinguishes four failure modes: fiber tension, fiber compression, matrix tension, and matrix compression. It is easy to apply and produces reasonable results in fiber-dominated failures. However, it lacks physical modeling of crack angles or matrix-dominated fracture evolution, limiting accuracy in delamination or inter-fiber failure.

Puck’s criterion, in contrast, introduces a detailed physical model for matrix failure, including the computation of the fracture plane angle based on shear and transverse stress components. It provides better predictions in matrix-dominated or off-axis loading scenarios and captures complex damage paths more accurately. However, Puck’s formulation is more complex, requires more input parameters (e.g., inter-fiber shear strength, fracture angle coefficients), and can increase computation time.

The LaRC04 and LaRC05 models, developed at NASA, refine matrix and fiber kinking modeling and include more robust criteria for compressive failure. These models have shown excellent correlation in open-hole compression (OHC), CAI, and buckling cases. However, their implementation is not yet standardized in most commercial solvers and often requires user-defined subroutines.

Comparative studies (e.g., WWFE-II [31]) reveal that all three models can match experimental data under specific conditions but differ significantly in their ability to generalize. Hashin tends to underpredict damage in matrix-dominated cases. Puck and LaRC handle these cases better but are harder to implement. For high-fidelity predictions, LaRC05 appears to provide the best balance of physical accuracy and general applicability. Hinton et al. [31] compared 19 different 3D failure criteria against a common set of experimental benchmarks – showing that while many criteria could be calibrated to specific failure cases, their predictive accuracy varied widely outside those conditions.

In practice, modelers often select a criterion based on available test data, implementation effort, and expected failure mode. For early-stage design, Hashin may suffice. For final certification or complex loadings, Puck or LaRC-based models offer better accuracy, provided proper material calibration is available. Details on these and other criteria, along with earlier comparative results from the first World-Wide Failure Exercise, can be found in the reference volume by Hinton et al. [31]. For instance, Dogan et al. [60] directly compared Hashin and Puck criteria in simulations of pin-loaded composite plates. They found that Puck’s formulation (with its angle-dependent matrix failure model) predicted the onset of damage and ultimate load more accurately, matching the test failure mode (matrix cracking followed by net-tension) better than Hashin’s criterion, which tended to underestimate the extent of matrix damage. Table 5 compares these criteria qualitatively.

Composite structures in service are exposed to cyclic loads, high strain rates, and environmental factors such as temperature and moisture – all of which influence failure. Fatigue modeling in composites remains an evolving area. Zhou et al. [70] recently demonstrated a progressive fatigue damage model on CFRP bolted joints, showing that by degrading stiffness each cycle based on an S–N curve, their simulation could estimate the number of cycles to failure within a factor-of-two of the test results. However, they also noted significant scatter in experimental fatigue lives, implying that while the model captured average behavior, conservative safety factors are still required. Many progressive damage models apply a stiffness or strength degradation law per load cycle, or use residual strength drop after a predefined number of cycles. Paris law-based delamination growth models or strain-based S–N curves are used, though experimental calibration is limited.

Strain rate sensitivity affects matrix-dominated properties such as interlaminar shear and compressive strength. Dynamic impact simulations (e.g., ballistic or fan blade-out events) show that higher strain rates slightly increase strength and delay damage. For example, Giannaros et al. [71] simulated an engine fan-blade fragment impact on a containment structure and found that incorporating strain-rate-dependent material properties was essential. Their high-velocity impact model, which included rate-sensitive composite strength, successfully reproduced the experimentally measured damage extent and absorbed energy, whereas a rate-insensitive model mispredicted both the damage extent and the absorbed energy. Accurate modeling requires rate-dependent stress-strain and fracture data, which are scarce and test-intensive.

Temperature and humidity reduce matrix modulus, glass transition temperature, and fracture toughness. Hot-wet conditions, particularly in thermoset composites, lower delamination resistance and promote microcracking. Some models scale damage evolution laws based on temperature or introduce moisture diffusion-coupled degradation. Elamin et al. [72] examined foam-core sandwich composites in arctic conditions (down to –60°C). Their tests showed reduced ductility and a shift in failure mode (more brittle fiber fracture), which their thermal-dependent material model was able to replicate. This study provides a validation point that even at extremely low temperatures, appropriately modified material properties in the FE model can match the observed loss of composite toughness.

Most commercial FEA solvers do not yet support full environmental interaction models for composites. Therefore, degradation is often included indirectly by modifying material parameters. High-fidelity analysis of service aging still requires combined modeling–testing workflows.

Despite progress, several challenges remain in predictive composite failure simulation:

• Multiscale modeling: Full-field simulations at the ply or sub-ply level are computationally expensive. Efficient global-local approaches, homogenization, or reduced-order models are needed.

• Probabilistic and defect-tolerant analysis: Real structures contain manufacturing defects. Probabilistic methods, stochastic FEM, and defect-informed modeling (e.g., CT-informed models) are under development.

• Automation and material databases: Many failure models require extensive calibration. Integrated databases with standardized material cards for multiple solvers (e.g., Abaqus, LS-DYNA, Ansys) could speed up adoption.

• Digital twins and SHM integration: Real-time structural health monitoring (SHM) data could be linked to validated FEM models to create digital twins that assess residual strength and support maintenance decisions.

• Certification by analysis: Regulatory acceptance of simulation as a partial substitute for destructive testing is increasing. However, clear validation paths and uncertainty quantification remain necessary. Continued efforts in benchmarking (e.g., CAI round robins, WWFE exercises), sharing of test data, and development of user-friendly simulation frameworks will be key to overcoming these challenges.