Analysis of composite/titanium alloy hybrid bolted structures based on improved LaRC05 criterion and backpropagation neural network

Abstract

Reliable prediction of ultimate load in composite/titanium bolted joints is hindered by intricate, interacting damage modes under multi-axial stress. To overcome the limitations of current high-fidelity yet computationally expensive models, we present a hybrid data-physics framework that couples an enhanced shear-driven three-dimensional LaRC05 failure criterion with a streamlined backpropagation (BP) neural network. This synergistic coupling represents a paradigm shift, enabling a leap from high-cost, high-fidelity simulation to high-fidelity, low-cost prediction. The progressive damage process is embedded in ABAQUS through a user-defined material subroutine, accurately capturing damage initiation and evolution without resorting to excessive mesh refinement. Bayesian optimisation and Shapley Additive Explanations (SHAP)-based feature selection yield a compact network architecture that retains only the most influential inputs, ensuring robust generalisation across a broad design space. Validation against experimental data demonstrates markedly improved accuracy and computational efficiency, enabling rapid evaluation during early-stage design. The resulting tool is readily deployable, offering practitioners a swift and dependable route for the safety assessment of lightweight hybrid structures.

Keywords

bolted joints composite materials failure prediction BP neural network ultimate load

Introduction

The research on hybrid structures of composite and metal bolted connections has made remarkable progress in recent years, but the reliability of the connection is still a core problem that needs to be solved in engineering practice. Bolted connections are widely used in various types of engineering due to their removable characteristics, but the anisotropic properties of composites often lead to complex multi-modal damage around the connection holes.^1–3

Researchers have studied metal-composite hybrid connections under static, fatigue, and pressure loading to understand damage mechanisms.⁴ Olivier et al.⁵ examined FRP-steel bolted joints under static and long-term loading, while Kabche et al.⁶ analysed hybrid joints under uniform pressure. Combining metal strength with composite specific strength, these connections are vital for lightweight, high-performance structures. Luo et al.⁷ studied hybrid connections under bending, and Caccese et al.⁸ investigated large panels under uniform pressure. Feldman⁹ pioneered removable stainless-steel bushings, effectively transferring extrusion loads to metal parts and increasing hole load capacity. Elhady et al.¹⁰ found fuselage cracks can propagate under fatigue, affecting safety. Camanho et al.¹¹ developed monolithic bonded aluminum alloy bushings, increasing model strength by 33% and achieving seamless bonding. Junhong Li et al.¹² proposed a novel titanium-to-steel joining method with superior loosening resistance. Akbarpour¹³ developed a co-curing embedded stainless steel foil technique, boosting load capacity by 60% but revealing foil interlayer peeling failure. Consequently, Liu Yuwei et al.¹⁴ improved the process with titanium alloy localised doping and precision drilling, enlarging the bonding region to improve load transfer and suppress sudden failure. However, a systematic framework for matching foil geometry parameters to solidification kinetics in doping is lacking; optimal design and detection systems are needed to improve reliability and safety.

ANN models, advantageous for repetitive tasks and complex nonlinear relationships, are increasingly used in composite mechanics for intelligent design and performance prediction. Kazi et al.¹⁵ built a multi-parameter prediction system for cotton fibre composites, including filler content, load–displacement, and crashworthiness optimisation. Liu’s team¹⁶ applied ANN to predict composite stress-strain response using empirical formula databases. For failure analysis, ANN offers limited efficiency gain with the fast Hashin criterion,¹⁷ but its predictive advantage is prominent with iterative Puck and LaRC05 criteria for damage index and fracture angle.¹⁸

The difference between the anisotropic properties of composites and the intrinsic structure of metals makes stress concentration and damage evolution mechanisms at the interface of the two materials complex and variable. For this reason, based on the development of damage prediction models and intelligent assessment techniques for physical mechanisms, this paper aims to develop an efficient and stable composite material ultimate load assessment method for the UMAT user material subroutine for failure prediction in the three-dimensional state, combining the shear-driven effect with the application of the improved LaRC05 criterion and the adaptive stiffness degradation model, and expanding the failure prediction of two-dimensional open-cell laminates to three-dimensional composites/titanium alloys with hybrid bolted joint structures. To develop an efficient and stable method to evaluate the ultimate load of composite materials, combine the prediction ability of neural network with the mechanical properties of materials, and construct an intelligent evaluation system to solve the computational bottleneck in the traditional analysis method.

Prediction model of backpropagation neural network based on TCT double-lap structure

LaRC05 failure initiation criterion

LaRC05 criterion is developed by NASA Langley Research Center on the basis of Puck failure theory, which is a composite material progressive failure criterion. This criterion is able to characterise the mechanical properties of thin-layer structures strengthened by interlayer constraints and has become an important theoretical tool for the refinement design of modern aerospace composite structures.

The traditional LaRC05 criterion defines the following major failure modes:

(1) Fibre Tension Failure:

This failure mode typically occurs in composite materials where the fibres may break or undergo tensile failure when the fibre portion of the material is subjected to stresses in excess of its tensile strength. This failure mode usually affects the load carrying capacity of the material, especially in parts of the aircraft structure that are subjected to large tensile forces.

f_{f t} = \frac{σ_{11}}{X_{T}} (σ_{11} \geq 0)

(1)

where $σ_{11}$ for the fibre main direction stress, $X_{T}$ for the longitudinal tensile strength, $f_{f t} \geq 1$ when the failure is judged.

(2) Fibre Compression Failure:

Fibre compression failure occurs when the composite material is subjected to compressive stresses, and the fibres may yield or fail in compression. This failure mode usually occurs in areas subjected to high compressive forces.

f_{k i n k} = f_{s p l i t} = {(\frac{τ_{12}^{m}}{S_{L}^{i s} - η_{L} σ_{2}^{m}})}^{2} + {(\frac{τ_{23}^{m}}{S_{T}^{i s} - η_{T} σ_{2}^{m}})}^{2} + {(\frac{σ {_{2}^{m}}_{+}}{Y_{T}^{i s}})}^{2} (σ_{11} < 0)

(2)

When $| σ_{1} | ⩽ | \frac{X_{C}}{2} |$ , the splitting mode failure of the fibre is predicted; when $| σ_{1} | > | \frac{X_{C}}{2} |$ , the kinking mode failure of the fibre is predicted. Where $f_{k i n k}$ is the fibre bending, $f_{s p l i t}$ is the splitting failure index. $X_{C}$ is the longitudinal compressive strength of the layup. $τ_{12}^{m}$ , $τ_{23}^{m}$ and $σ_{2}^{m}$ are the stresses on the kinking surface, which are obtained by converting the coordinate system of the initial direction of the fibre to the fibre deflection coordinate system. $S_{L}^{i s}$ is the longitudinal shear strength, $S_{T}^{i s}$ is the transverse shear strength, $Y_{T}^{i s}$ is the normal tensile strength at the fracture surface taking into account the in-situ effect. $η_{T}$ and $η_{L}$ are the slope and friction coefficient to characterise the effect of compression on failure.

There is an initial angular deviation in fibre orientation in fibre-reinforced composites, and under longitudinal compressive loading, this initial deviation leads to fibre rotation. Conventional fibre bending theory usually assumes that deformation occurs within the plywood face, and a large number of experimental studies have simulated this in-face deformation mechanism by constraining the out-of-face displacement of the specimen. However, recent studies have shown that the three-dimensional characteristics of fibre bending damage cannot be ignored. Based on the two-dimensional theoretical framework, Pinho et al.^18,19 developed a three-dimensional bending model analysis method. The failure mode has similarity with the mechanical mechanism of matrix damage, and the difference between the two lies in the selection of the reference plane for stress calculation. The matrix failure is calculated based on the fracture surface, while the fibre bending failure uses the bending surface as the basis for stress calculation. As shown Figure 1, the coordinate system 1–2–3 represents the material coordinate system, and $ψ$ is the angle between the bending plane and the 2-axis. Fibre bending occurs in the plane of the local coordinate system $1^{m} - 2^{m} - 3^{m}$ , which is obtained by rotating the angle $φ$ from the coordinate system $1^{ψ} - 2^{ψ} - 3^{ψ}$ .

{\begin{matrix} σ_{22}^{ψ} = σ_{22} \cos^{2} ψ + σ_{33} \sin^{2} ψ + 2 τ_{23} \sin ψ \cos ψ \\ τ_{12}^{ψ} = τ_{12} \cos ψ + τ_{31} \sin ψ \\ τ_{23}^{ψ} = (σ_{33} - σ_{22}) \sin ψ \cos ψ + τ_{23} (\cos^{2} ψ - \sin^{2} ψ) \\ τ_{31}^{ψ} = τ_{31} \cos ψ - τ_{12} \sin ψ \end{matrix}

(3)

Since the bending bands are located in the same bending plane, $τ_{23}^{ψ} = 0$ . The angle of the bending plane $ψ$ is determined by the equation:

\tan 2 ψ = \frac{2 τ_{23}}{σ_{22} - σ_{33}}

(4)

The stresses are then converted to the fibre deflection coordinate system according to equation (5):

{\begin{matrix} σ_{22}^{m} = σ_{11} \sin^{2} φ + σ_{22} \cos^{2} φ - 2 τ_{12}^{ψ} \sin φ \cos φ \\ τ_{12}^{m} = (σ_{22}^{ψ} - σ_{11}) \sin φ \cos φ + τ_{12}^{ψ} (\cos^{2} φ - \sin^{2} φ) \\ τ_{23}^{m} = τ_{23}^{ψ} \cos φ - τ_{31}^{ψ} \sin φ \end{matrix}

(5)

where m coordinate system represents the fibre deflection coordinate system, $σ_{22}^{m}$ , $τ_{12}^{m}$ , $τ_{23}^{m}$ for the fibre deflection coordinate system under the stress component.

(3) Matrix Failure

Matrix failure refers to the destruction of the matrix of a composite material when it is subjected to shear stress. This is usually manifested as a shear fracture in the bonding region of the material, which leads to the overall structural instability of the material. Tensile and compressive failure of the matrix can be given by an expression:

f_{m} = {(\frac{τ_{T}}{S_{T}^{i s} - η_{T} σ_{n}})}^{2} + {(\frac{τ_{L}}{S_{L}^{i s} - η_{L} σ_{n}})}^{2} + {(\frac{{⟨ σ_{n} ⟩}_{+}}{Y_{T}^{i s}})}^{2}

(6)

where $τ_{T}$ is the transverse shear stress, $τ_{L}$ is the transverse normal stress, and $σ_{n}$ is the positive stress acting on the fracture surface. When $f_{m} \geq 1$ , it means that matrix damage occurs.

(4) The improved LaRC05 criterion considers the influence of shear drive on the failure judgement based on the traditional criterion, and the criteria are as follows:

f_{s h e a r - d r i v e n} = | \frac{σ_{11}}{X_{C}} | + {| \frac{τ_{12}}{S_{L}^{i s}} |}^{n}_{\begin{matrix}  \end{matrix} | σ_{11} | > 0.8 X_{^{C}} : f i b r e k i n k i n g}^{N o t e s : \begin{matrix} | σ_{11} | \leq 0.8 X_{^{C}} : f i b r e s p l i t; \end{matrix}}

(7)

f_{f c} = max {f_{k i n k}, f_{s h e a r - driven}}

(8)

where $f_{s h e a r - d r i v e n}$ is the fibre shear damage failure judgement. In addition, to limit the range of pure fibre compression strengths, the second term in equation (7) is used to consider the effect of shear on initial failure. The introduction of the shear-driven term and the parameter *n* follows the physical rationale developed by Pinho et al.¹⁸ for capturing the effect of shear stress on fibre compression failure, and its numerical implementation is adapted from the work of Schirmaier et al.²⁰ for efficient fracture angle search. n represents the power of the contribution of shear failure initiation, n ≥ 2. As n decreases, the predicted pure fibre compressive strength tends to Xc. In the absence of experimental data, the model is adaptable as it is able to degrade to a more conservative form of fibre compression damage prediction based on the predicted initial fibre deflection angle.

Meanwhile, in order to improve the computational efficiency, the one-dimensional extended golden section method is selected to search for the fracture angle α.²⁰ This method is based on the golden section search method, which utilises three points on the search sub-interval to fit a parabola $P (α_{1}^{'}, f_{m} (α_{1}^{'})), P (α_{2}^{'}, f_{m} (α_{2}^{'})), P (α_{3}^{'}, f_{m} (α_{3}^{'}))$ , fm is available at the extreme points of the parabola:

α \approx α_{2}^{'} - \frac{1}{2} \frac{{(α_{2}^{'} - α_{1}^{'})}^{2} (f_{m} (α_{2}^{'}) - f_{m} (α_{3}^{'})) -}{(α_{2}^{'} - α_{1}^{'}) (f_{m} (α_{2}^{'}) - f_{m} (α_{3}^{'})) -} \frac{{(α_{2}^{'} - α_{3}^{'})}^{2} (f_{m} (α_{2}^{'}) - f_{m} (α_{1}^{'}))}{(α_{2}^{'} - α_{3}^{'}) (f_{m} (α_{2}^{'}) - f_{m} (α_{1}^{'}))}

(9)

where α is the fracture angle under any stress state, p(α) is approximated with fm, and the algorithm terminates the iteration when the difference between the extreme value points of the two search sub-intervals reaches the preset accuracy threshold. At this point, when the plywood substrate undergoes tensile or compressive failure, the fracture angle on the corresponding potential fracture surface will remain unchanged.

Under complex stress states, the fracture at material point failure $α$ is not necessarily equal to 53° or 0°. The angle that maximises the matrix failure indicator fm is searched in the range [0°, 180°] for any stress state is the angle of the fracture surface at the current stress state. In Figure 2, the golden section search method assumes that two points $α_{3}$ and $α_{4}$ are found in the search interval [ $α_{1}$ , $α_{2}$ ] according to the equation such that $\frac{b}{a} = φ$ , and $\frac{b}{a} = \frac{a}{c}$ , where $φ = \frac{1 + \sqrt{5}}{2}$ is the ratio.

{\begin{matrix} r = \frac{\sqrt{5} - 1}{2} \\ α_{3} = α_{1} + (1 - r) (α_{2} - α_{1}) \\ α_{4} = α_{1} + r (α_{2} - α_{1}) \end{matrix}

(10)

If $f_{m c} (α_{4}) \geq f_{m c} (α_{3})$ , then the maximum value occurs in $[α_{3}, α_{2}]$ , and the search interval is narrowed to $[α_{3}, α_{2}]$ . If $[α_{1}, α_{4}]$ , it means that the largest value appears in $[α_{1}, α_{4}]$ , then the search interval is reduced to $[α_{1}, α_{4}]$ , and so on iteratively reduce the search space until the search result meets the accuracy requirements.

The implementation of golden section algorithm needs to meet the requirement of single peak characteristic of the objective function. The distribution law of matrix failure index with azimuthal angle under four typical stress conditions is shown, in which the criterion for determining the fracture surface is: the azimuthal angle that makes the failure index reach the extreme value is selected as the fracture angle.

Adaptive material property degradation method

Under different failure modes, the failure mechanism and driving strain of damage evolution are different, so the corresponding material property degradation schemes are also different. In this paper, the damage evolution method from the literature is applied.²¹ The definitions of driving strains and damage factors for each failure mode in the improved LaRC05 failure initiation criterion are given below.

(1) Matrix damage factor dmat

Matrix failure is influenced by the stress component on the fracture surface, and the associated matrix damage evolution is also driven by the strain on the fracture surface. Therefore, the driving stress σmat is defined in the matrix failure mode:

σ_{m a t} = \sqrt{{⟨ σ_{n} ⟩}^{2} + τ_{m a t}^{2}}

(11)

The damage factors for matrix compression and matrix tension can be expressed uniformly as:

d_{m a t} = \frac{ε_{f}^{e, m a t} (ε_{max}^{e, m a t} - ε_{0}^{e, m a t})}{ε_{max}^{e, m a t} (ε_{f}^{e, m a t} - ε_{0}^{e, m a t})}

(12)

dmat= dmt in tension; dmat= dmc in compression.

In the formula:

ε_{f}^{e, m a t} = \frac{2 G_{C}^{m a t}}{σ_{0}^{m a t} l *}

(13)

where $ε_{f}^{e, m a t}$ is the driving strain when the matrix fails completely; $l *$ is the characteristic length of the unit, which is obtained directly in the finite element software.

G_{C}^{m a t} = G_{I C}^{m a t} {(\frac{σ_{n}^{0}}{σ_{0}^{m a t}})}^{2} + G_{I I C}^{m a t} {(\frac{τ_{T}^{0}}{σ_{0}^{m a t}})}^{2} + G_{I I C}^{m a t} {(\frac{τ_{L}^{0}}{σ_{0}^{m a t}})}^{2}

(14)

where $σ_{n}^{0}$ , $τ_{T}^{0}$ , $τ_{T}^{0}$ are the stress components on the fracture surface at the onset of failure; $G_{I C}^{m a t}$ , $G_{I I C}^{m a t}$ are the fracture toughnesses of the $I$ and $II$ types; $G_{I C}$ , $G_{II C}$ can be used instead of $G_{I C}^{m a t}$ and $G_{II C}^{m a t}$ ; $G_{C}^{m a t}$ is the energy release at the time of complete failure of the matrix.

ε_{max}^{e, m a t} = max_{t * \leq t} {ε_{m a t}^{e}}

(15)

$t *$ , t is the current moment; $ε_{max}^{e, m a t}$ is the maximum driving strain corresponding to the matrix damage.

σ_{0}^{m a t} = σ_{m a t} |_{f_{m} = 1}

(16)

$σ_{0}^{m a t}$ is the failure onset driving stress.

ε_{0}^{e, m a t} = ε_{m a t}^{e} |_{f_{m} = 1}

(17)

$ε_{0}^{e, m a t}$ is the onset driving strain corresponding to matrix damage; $ε_{m a t}^{e}$ is the driving strain in matrix failure mode.

ε_{m a t}^{e} = \frac{σ_{n}}{σ_{n}} ε_{n}^{e} \sin ω + γ_{m a t}^{e} \cos ω

(18)

$γ_{m a t}^{e}$ is the shear driving strain in the direction of the shear driving stress at the fracture surface.

γ_{m a t}^{e} = | γ_{T}^{e} \cos λ + γ_{L}^{e} \sin λ |

(19)

The angle between $τ_{T}$ and $τ$ is defined as $λ$ and the angle between $τ_{m a t}$ and $σ_{m a t}$ is defined as $ω$ at the action surface with the fracture angle $α$ . The elastic strain components $ε_{n}^{e}$ , $γ_{T}^{e}$ , $γ_{L}^{e}$ are expressed as follows:

ε_{n}^{e} = ε_{22}^{e} \cos^{2} α + ε_{33}^{e} \sin^{2} α + \frac{1}{2} γ_{23}^{e} \sin 2 α

(20)

γ_{T}^{e} = \sin 2 α (- ε_{22}^{e} + ε_{33}^{e}) + γ_{23}^{e} \cos 2 α

(21)

γ_{L}^{e} = γ_{12}^{e} \cos α + γ_{13}^{e} \sin α

(22)

$σ_{n}$ , $τ_{T}$ , $τ_{L}$ – for the stress component on the fracture surface, obtained by the stress transformation in the material coordinate system:

{\begin{matrix} σ_{n} = σ_{22} \cos^{2} α + σ_{33} \sin^{2} α + τ_{23} \sin 2 α \\ τ_{T} = \sin α \cos α (- σ_{22} + σ_{33}) + τ_{23} \cos 2 α \\ τ_{L} = τ_{12} \cos α + τ_{13} \sin α \end{matrix}

(23)

(2) Fibre tensile damage factor $d_{f t}$

d_{f t} = \frac{ε_{f}^{e, f t} (ε_{max}^{e, f t} - ε_{0}^{e, f t})}{ε_{max}^{e, f t} (ε_{f}^{e, f t} - ε_{0}^{e, f t})}, ε_{f}^{e, f t} = \frac{2 G_{C}^{f t}}{σ_{0}^{f t} l *}

(24)

$ε_{f}^{e, f t}$ – for the fibre tensile complete failure of the driving strain; $G_{C}^{f t}$ – for the fibre tensile fracture toughness, the material properties table has given the value; $ε_{\max}^{e, f t}$ – for the fibre tensile maximum driving strain.

ε_{max}^{e, f t} = max_{t * \leq t} {ε_{11}^{e}}

(25)

$ε_{0}^{e, f t}$ is the fibre tensile failure index; $f_{f t} = 1$ is the corresponding starting driving strain at failure:

ε_{0}^{e, f t} = ε_{11}^{e} |_{f_{f t} = 1}

(26)

$σ_{f t}$ is the corresponding failure onset driving stress.

σ_{0}^{f t} = σ_{11} |_{f_{f t} = 1}

(27)

The expressions for driving stress $σ_{f t}$ and driving strain $ε_{f t}^{e}$ are as follows:

σ_{f t} = σ_{11}

(28)

ε_{f t}^{e} = ε_{11}^{e}

(29)

(3) Fibre splitting damage factor dsplit

d_{s p i l t} = \frac{ε_{f}^{e, s p l i t} (ε_{max}^{e, s p l i t} - ε_{0}^{e, s p l i t})}{ε_{max}^{e, s p l i t} (ε_{f}^{e, s p l i t} - ε_{0}^{e, s p l i t})}

(30)

In the formula:

ε_{f}^{e, k i n k} = \frac{2 G_{C}^{f c}}{σ_{0}^{k i n k} l *}

(31)

$ε_{f}^{e, k i n k}$ – for the fibre compression complete failure of the driving strain; $ε_{max}^{e, k i n k}$ – for the maximum driving strain in the loading course.

ε_{max}^{e, k i n k} = max_{t * \leq t} {ε_{k i n k}^{e} |}

(32)

$σ_{0}^{k i n k}$ – for the fibre compression bending failure index; $f_{k i n k} = 1$ when the corresponding failure of the starting driving stress:

σ_{0}^{k i n k} = σ_{k i n k} |_{f_{k i n k} = 1}

(33)

$ε_{0}^{e, k i n k}$ – is the corresponding starting driving strain.

Because the shear force in the bending surface causes the accumulation of fibre deflection, as well as the longitudinal compression damage in the fibre, which ultimately leads to fibre bending damage, so the driving stress $σ_{k i n k}$ is defined as:

σ_{k i n k} = \sqrt{{(σ_{11}^{m})}^{2} + {(τ_{12}^{m})}^{2}}

(34)

The driving strain $ε_{k i n k}^{e}$ is:

ε_{k i n k}^{e} = \frac{| σ_{11}^{m} ε_{11}^{e, m} | + | τ_{12}^{m} γ_{12}^{e, m} |}{σ_{k i n k}}

(35)

$σ_{11}^{m}$ Therefore, the driving stress $γ_{12}^{e, m}$ is defined as: driving strain $ε_{11}^{e, m}$ is: $τ_{12}^{m}$ for the stress component in the fibre deflection coordinate system; for the elastic strain component in the fibre deflection coordinate system.

Therefore, the strain is first converted to the $ψ$ coordinate system where the bending plane is located:

ε_{22}^{e, ψ} = ε_{22}^{e} \cos^{2} ψ + ε_{33}^{e} \sin^{2} ψ + γ_{23}^{e} \sin ψ \cos ψ

(36)

γ_{12}^{e, ψ} = γ_{12}^{e} \cos ψ + γ_{13}^{e} \sin ψ

(37)

γ_{13}^{e, ψ} = γ_{13}^{e} \cos ψ - γ_{12}^{e} \sin ψ

(38)

γ_{23}^{e, ψ} = 0

(39)

Then the strain is converted to the deflected fibre in the m coordinate system:

ε_{11}^{e, m} = ε_{11}^{e} \cos^{2} θ + ε_{22}^{e, ψ} \sin^{2} θ + γ_{12}^{e, ψ} \sin θ \cos θ

(40)

ε_{22}^{e, m} = ε_{11}^{e} \sin^{2} θ + ε_{22}^{e, ψ} \cos^{2} θ - γ_{12}^{e, ψ} \sin θ \cos θ

(41)

γ_{12}^{e, m} = \sin 2 θ (- ε_{11}^{e} + ε_{22}^{e, ψ}) + γ_{12}^{e, ψ} (\cos^{2} θ - \sin^{2} θ)

(42)

Fibre compression bending damage factor $d_{k i n k}$ :

d_{k i n k} = \frac{ε_{f}^{e, k i n k} (ε_{max}^{e, k i n k} - ε_{0}^{e, k i n k})}{ε_{max}^{e, k i n k} (ε_{f}^{e, k i n k} - ε_{0}^{e, k i n k})}

(43)

Similar to the fibre compression bending failure mode, the driving stresses $σ_{s p l i t}$ and driving strains $ε_{s p i l t}^{e}$ for the splitting failure mode are given directly below:

σ_{s p l i t} = \sqrt{{⟨ σ_{22}^{m} ⟩}^{2} + {(τ_{12}^{m})}^{2}}

(44)

ε_{s p l i t}^{e} = \frac{| ⟨ σ_{22}^{m} ⟩ ε_{22}^{e, m} | + | τ_{12}^{m} γ_{12}^{e, m} |}{σ_{s p l i t}}

(45)

(4) Shear-driven damage factor dsd

d_{s d} = \frac{ε_{f}^{e, s d} (ε_{\max}^{e, s d} - ε_{0}^{e, s d})}{ε_{\max}^{e, s d} (ε_{f}^{e, s d} - ε_{0}^{e, s d})}

(46)

where the influence of the normal strain component on the acting surface on the matrix under pressure on the crack closure effect is analysed by introducing the $⟨ σ_{n} ⟩$ term.

τ_{m a t} = \sqrt{τ_{T}^{2} + τ_{L}^{2}}

(47)

$τ_{m a t}$ is the shear driving stress; $σ_{n}$ , $τ_{T}$ and $τ_{L}$ are the stress components on the fracture surface.

ε_{f}^{e, s d} = \frac{2 G_{C}^{f t}}{σ_{0}^{e, s d} l *}

(48)

$ε_{f}^{e, s d}$ is the driving strain at complete destruction; $G_{C}^{f t}$ is the fibre compression fracture toughness.

Similarly, $σ_{0}^{e, s d}$ is the driving stress at failure initiation corresponding to $f_{s h e a r - d r i v e n} = 1$ , and the driving strain at failure initiation corresponding to $ε_{0}^{e, s d}$ . The driving strain at failure initiation corresponding to is the driving strain at failure initiation:

σ_{0}^{e, s d} = \frac{σ^{e, s d}}{f_{s h e a r - d r i v e n}}

(49)

ε_{dr}^{sd, 0} = \frac{ε^{e, s d}}{f_{s h e a r - d r i v e n}}

(50)

ε_{\max}^{e, s d} = max_{t * \leq t} {ε^{e, s d}}

(51)

$σ^{e, s d}$ , $ε^{e, s d}$ is the elastic strain component acting on the fracture surface; $ε_{\max}^{e, s d}$ is the maximum value of the driving strain reached during the loading history.

The above parameters can be obtained in the following three ways: (1) directly from the material property tables, e.g., $G_{I C}$ and $G_{II C}$ , $σ_{22}$ ; (2) dynamically stored and called as state variables in the UMAT subroutine; (3) derived from other known parameters.

Based on the theoretical foundation of Kachanov²² and Matzen miller,²³ the extension of the dimension of damage response analysis is realised by extending the classical damage model to three-dimensional stress states. The model establishes a mathematical correlation between the real stress and the nominal stress by introducing the damage tensor, and adopts three damage variables, $d_{f}$ , $d_{m}$ and $d_{s}$ , to characterise the evolution of fibre orientation damage, matrix phase damage and shear damage, respectively. The weakening effect of fibre bending on the shear properties was determined by $f_{k i n k}$ as shown in equation (52);

\begin{aligned} {\begin{matrix} d_{f} = d_{f t} + d_{s d} - d_{f t} d_{s d}; f_{f c} = f_{s h e a r - d r i v e n} \geq 1 > f_{k i n k} \\ d_{f} = d_{f t} + d_{k i n k} - d_{f t} d_{k i n k}; f_{f c} = f_{k i n k} \geq 1 > f_{s h e a r - d r i v e n} \end{matrix} \\ {\begin{matrix} d_{m} = 1 - (1 - d_{m t}) (1 - d_{m c}); \\ d_{s} = 1 - (1 - d_{f}) (1 - d_{m}) \end{matrix} \end{aligned}

(52)

The damage model adopts the orthogonal anisotropic principal structure,²⁴ which effectively reduces the mesh sensitivity of numerical calculation by introducing the characteristic length to normalise the critical energy in the damage evolution process. After the material failure, different damage evolution models are used to degrade the corresponding material properties, which can effectively characterise the softening behaviour of the material during the damage accumulation process and accurately describe the degradation characteristics of different mechanical properties.

UMAT user material subroutine implementation

Considering the real degradation mode of the material and the fact that the inter-fibre fracture criterion needs to consider the stress on the non-fracture surface, this paper writes a UMAT subroutine based on the improved LaRC05 failure initiation criterion. The subroutine defines the material model by constructing the constitutive equations rather than directly inputting specific material parameters. The chosen constitutive equations are orthogonal anisotropic models commonly found in carbon fibre-reinforced composite laminates. Based on the stress hazard factor calculated from the inter-fibre failure criterion, an adaptive material degradation model is used in this paper to simulate the failure process. The specific implementation flow of the UMAT subroutine is shown in Figure 3, which first defines the system and the custom parameters, and subsequently establishes the linear constitutive equations. For fibre-reinforced composites, orthogonal anisotropic eigen equation models were used.

Finite element model establishment

The TCT double lap structure is shown in Figure 4, the T700/2510 composite plywood used,²⁵ and its mechanical properties are shown in Table 1.

Table 1.

Composite material properties of T700/2510.²⁶

Material properties	Value	Material properties value	value
E₁₁ (GPa)	127	X_T (MPa)	2200
E₂₂ (GPa)	8.41	X_C (MPa)	1470
E₃₃ (GPa)	8.41	Y_T (MPa)	49
G₁₂ (GPa)	4.21	Y_C (MPa)	199
G₁₃ (GPa)	4.21	S (MPa)	150
G₂₃ (GPa)	3.40	G_f (N/mm)	133
ν ₁₂	0.31	G_kink (N/mm)	40
ν ₁₃	0.31	GIC (N/mm)	0.23
ν ₂₃	0.3	GIIC (N/mm)	0.78

The layup sequence of the laminate is [45/0/−45/90]6S. Where 6S denotes a symmetric laminate formed in 6 repetitions containing a total of 48 layers (4 layers × 6 repetitions × 2 symmetric). The top and bottom laps and fasteners were manufactured from TC21 titanium alloy, with a modulus of elasticity E2 of 110 GPa and a strength N2 of 1100 MPa. The thickness of the top and bottom laps was 3.0 mm. The thickness of the composite laminate was 3.12 mm, and the nominal thickness of the laminates was 0.065 mm. The hole diameters were 6.0 mm, and the distance from the edges was 20 mm. The widths of both the laminates and the top and bottom laps were 38.1 mm. The width of the plywood and the top and bottom laps are 38.1mm and the length is 135 mm.

This paper combines an improved LaRC05 criterion and adaptive material degradation theory to establish a damage analysis model for thin-walled composite laminates. Based on material constitutive relationships, a FORTRAN UMAT subroutine is developed and embedded in ABAQUS/Standard to simulate the quasi-static tensile response of a TCT double lap structure. Components (titanium plates, laminates, fasteners) are discretised with C3D8R elements. Mesh density is coarser away from bolt holes and finer around them, incorporating feature length treatment to reduce mesh sensitivity. Zero-thickness COH3D8 cohesive elements between plies model delamination: damage initiation follows a quadratic stress criterion, degradation uses linear stiffness reduction, and propagation is controlled by the BK mixed-mode criterion. Boundary conditions: a 6 kN bolt preload is applied; the left titanium plate is fully fixed; a 6 mm x-direction displacement load is applied to the right end of the intermediate laminate (see Figure 4). Seven contact pairs are defined, using the stiffer surface as master. Contact friction coefficients were assigned based on values commonly reported and validated in the literature for similar material interfaces in composite-metal joints^25,27: titanium plate/laminate = 0.2; bolt shank/laminate & bolt shank/titanium = 0.1; nut and washer/laminate = 0.2. Composite properties are in Table 1.

BP neural network architecture

As shown in Figure 5, the backpropagation (BP) neural network in this study adopts a multi-layer cascade topology, comprising an input layer, multiple hidden layers for feature transformation, and an output layer for ultimate load prediction. The information flows unidirectionally; the input layer receives and pre-processes the feature data, which is then passed through the hidden layers via weighted summation and nonlinear activation functions to generate the final prediction.

The specific architecture, including the number of hidden layers and neurons, the types of activation functions (e.g., Tansig, Logsig, ReLU), and the training algorithm, were not predetermined. Instead, they were systematically determined through a Bayesian optimisation process, detailed in section ‘Optimisation of BP neural network model architecture and hyperparameters’. This data-driven approach ensures the network configuration is optimally tuned for the specific problem of predicting the ultimate load of hybrid bolted structures.

In batch training, the total error is defined as:

E_{t o t a l} = \sqrt{\frac{1}{2} \sum_{P = 1}^{P} \sum_{k = 1}^{n} {({\hat{z}}_{k}^{P} - z_{k}^{P})}^{2}}

(53)

where P is the number of training samples.

The batch training method aims to reduce the global error and follows the principle of ‘collectivism’ to ensure that the total error changes in the direction of decreasing. When the number of samples is large, the convergence speed of batch training is significantly faster than that of single-sample training. The specific process is shown in Figure 6.

Data set preparation

The value ranges for these input features, as listed in Table 2, were defined to cover a broad and realistic design space relevant to aerospace composite structures. The geometric parameter ranges (L, W, h, D) were selected based on common dimensions and design guidelines for bolted joints in airframe components. The material property ranges were determined by considering the nominal properties of the T700/2510 system²⁶ alongside typical material scatter and the properties of other common aerospace-grade carbon/epoxy systems. This approach ensures that the trained neural network model possesses good generalisation capability for practical engineering applications beyond the specific baseline configuration.

Table 2.

Input features and value ranges of BP neural network model.

	Input features and value ranges of BP neural network model	Value ranges
Geometric parameters	Length of composite plywood board $l$	130 ∼ 160 mm
	Composite plywood board width $w$	35 ∼ 40 mm
	Composite plywood board thickness $h$	2 ∼ 5 mm
	Hole diameter $r$	4∼8
Material Properties	Longitudinal modulus of elasticity: E₁₁	120 ∼ 140 GPa
	Transverse modulus of elasticity: E₂ = E₃₃	7 ∼ 10 GPa
	In-plane shear modulus: G₂ = G₃	3 ∼ 4 GPa
	Outside shear modulus: G₂₃	4 ∼ 5 GPa
	In-plane Poisson’s ratio: V₁₂ = V₁₃	0.30 ∼ 0.35
	Outer Poisson’s ratio: V₂₃	0.32 ∼ 0.37
	Longitudinal tensile strength: X_T	1700 ∼ 1800 MPa
	Longitudinal compressive strength: X_C	600 ∼ 800 MPa
	Transverse tensile strength: Y_T	30 ∼ 40 MPa
	Transverse compression strength: Y_C	130 ∼ 150 MPa
	Shear strength: S₂ = S₂₃ = S₁₃	95 ∼ 105 MPa
	Longitudinal shear friction coefficient: η_l	0.076 ∼ 0.084
	Transverse shear friction coefficient: η_r	0.27 ∼ 0.31

Latin hypercube sampling (LHS) was used to generate evenly distributed training samples across the multi-dimensional parameter space. To address significant magnitude differences between input features that could bias the BP neural network and affect convergence/prediction accuracy, maximum-value normalisation was applied. The normalised dataset was randomly split: 80% for training (with 20% of this reserved as a validation set to improve generalisation) and 20% for testing. Bayesian optimisation was employed to systematically tune the BP network’s structure and hyperparameters, first constructing a surrogate model to estimate objective function performance. Model prediction ability was evaluated using the mean squared error (MSE) and the coefficient of determination (R²). The MSE measures the average squared difference between predicted and actual values during training, while R² (range 0–1) indicates prediction quality (higher values = better performance). The definitions are as follows:

MSE = \frac{1}{n} \sum_{i = 1}^{n} {(y_{i} - {\hat{y}}_{i})}^{2}

(54)

R^{2} = 1 - \frac{\sum_{i = 1}^{n} {(y_{i} - {\hat{y}}_{i})}^{2}}{\sum_{i = 1}^{n} {(y_{i} - \bar{y})}^{2}}

(55)

where $y_{i}$ is the actual value, ${\hat{y}}_{i}$ is the predicted value, $\bar{y}$ is the mean value of the actual value, n is the number of samples.

Analysis of simulation results

Validation of the progressive damage model

This section evaluates the improved model by comparing its predictions (using the enhanced LaRC05 criterion) with the traditional LaRC05 criterion and experimental data.²⁵ Figure 7 shows the improved LaRC05 criterion predicts the initial load–displacement response more accurately than the traditional model against three experimental datasets (Exp1-3). All load histories exhibit a three-stage evolution²⁷: initial loading (0–0.8 mm displacement), damage initiation/extension, and final failure. In the initial stage (load 0–15 kN), a clear linear relationship exists, with high consistency between experiments and simulations. Both criteria prediction curves coincide here and match experimental data well, accurately simulating the material’s elastic behaviour (constant stiffness, no damage).

When displacement reaches 0.8–1.0 mm, the curve slopes decrease, marking the damage initiation stage. Within 1.0–2.0 mm, nonlinear behaviour emerges as load growth slows. The improved LaRC05 criterion shows a smooth transition aligning with experimental trends, while the traditional criterion exhibits an abrupt change at ∼1.0 mm followed by steep linear growth, peaking at an unrealistic 25 kN (10% above experiments) at 1.8 mm. This discrepancy stems from the traditional criterion neglecting shear stress effects on initial damage. In the final failure stage (2.0–3.0 mm), the improved criterion closely matches the experiments’ smooth decline, whereas the traditional criterion shows unstable fluctuations. Analysis of Figure 7 confirms the improved criterion’s superiority: it predicts damage onset more accurately, enables stable elastic-to-damage transition (due to the new shear-driven initiation criterion), and captures evolution/ultimate failure reliably. The Modified LaRC05 criterion demonstrates superior accuracy over the Traditional LaRC05: Inflection point load: both predicted 26.05 kN (experimental 23.43 kN), with Modified reducing relative error from 11.18% to 5.93%; Inflection displacement: Modified predicted 0.89 mm (3.26% error) vs. Traditional's 0.96 mm (4.33% error) at experimental 0.92 mm; Ultimate load: Modified achieved 26.99 kN (1.60% error) vs. Traditional's 26.58 kN (3.10% error) against experimental 27.43 kN, yielding error reductions of 5.25% (inflection load), 1.07% (displacement), and 1.50% (ultimate load).It quantifies this: the improved criterion reduces relative errors versus the traditional criterion for inflection point load (5.93% vs. 11.18%), inflection point displacement (3.26% vs. 4.33%), and ultimate load (1.60% vs. 3.10%), demonstrating enhanced prediction reliability for structural safety assessment.

The improved LaRC05 criterion is more refined in capturing the damage at the inflection point, and its displacement–load curves in the transition region are smoother and more consistent with the experimental data. The model also introduces a degradation mechanism in the early stage of fibre compression, which enhances the ability to describe the nonlinear behaviour of the material. From the results of the computational accuracy comparison, a larger denominator parameter $X_{C}$ reduces the contribution from the fibre direction and makes the prediction more dependent on the shear-driven effect in the early stage. A smaller parameter $S_{L}^{is}$ enhances the influence of shear stress, and at larger n boosts the shear effect sensitivity to obtain earlier failure prediction. At small fibre deflection angles, the contribution of shear stress $τ_{12}$ is amplified and dominates the failure process. The improved criterion significantly improves the ability to recognise the failure value at the critical point and effectively avoids the phenomenon of sudden load drop, thus enhancing the overall stability of the prediction process.

As shown in Figures 8 and 9, the prediction results of fibre damage (SDV61) and matrix damage (SDV62) for laminates at the inflection point are demonstrated for different layup angles of the conventional and improved LaRC05 guidelines. The cloud plots in the figure are represented in red–blue order, with red indicating the region where the damage occurs and blue indicating that the material remains intact. From the figure, it shows that the improved LaRC05 criterion predicts a larger region of fibre damage and a more coherent damage distribution in the ±45° orientation compared to the conventional criterion. In contrast, the damage prediction of the conventional criterion is mainly in the localised region of stress concentration. In the 0° and 90° layers, fibre breakage is mostly concentrated above and below the hole edge or on the left and right sides. In this part, the predictions of the two guidelines are basically the same, especially on the 0° layer, and both accurately identify the location of fibre damage on the upper side of the hole, which is subjected to the highest longitudinal tensile stress. However, the advantage of the improved criterion is that it has a smoother transition to the damage boundary, which better reproduces the ‘intermediate state’ of fibre breakage. The difference in matrix damage prediction between the two is more pronounced in the ±45° layer. The traditional LaRC05 criterion predicts a smaller range of damage. The improved criterion identifies a larger and more continuous region of matrix damage. Most of the damage extended radially outward from the hole edge, matching the shear yielding and cracking paths observed in the experiments. In the 0° layup, the improved LaRC05 criterion points out matrix damage on the left and right sides of the hole edge, and it is also sensitive to capture matrix cracking due to transverse stresses. The improved LaRC05 criterion also predicts a wider range of matrix damage in 90° layers. The improved LaRC05 criterion is closer to the actual physical process in damage initiation simulation, which improves the accuracy of prediction. Therefore, it provides theoretical support for the safety assessment of composite structures.

The prediction accuracy of the modified LaRC05 criterion is significantly improved, but there is still an error of up to 5.93% in the prediction of the ultimate load and the load at the turning points. This error comes from the sensitivity of the bolt-to-hole clearance, which is usually idealised to zero in FEA. This affects the realism of the simulation to some extent. However, the overall results show that the improved model has higher accuracy in predicting the failure of connectors and can reasonably reflect the trend of their failure.

Analysis of plywood damage evolution

The failure modes of composite thin plywood include fibre breaks, fibre kinks, matrix cracks, and combinations of these failure modes. Since the severe failure region is concentrated near the intermediate holes, this section focuses on analysing the failure mechanism of the intermediate plywood of the TCT double-lap structure under different degrees. The experimental and finite element simulation results, as shown in the comparison of Figure 10. The key features of the composite laminates at the final failure stage are revealed. The hole perimeter region presents a complex combination of multiple failure modes, and this combination of failure modes is the result of the material being subjected to complex stress states in the high stress concentration region.

Stress concentration triggers the damage initiation near the hole and subsequently undergoes an evolutionary process. In the experimental image (left panel), the failure initiation phase starts with the formation of localised matrix cracks at the hole perimeter. With increasing load, the cracks connect and expand into a crack network in the region of the hole edge perpendicular to the direction of the main load, the fibres undergo kinking and local buckling, and this kinking phenomenon is manifested as a material buildup in the experiment. This evolutionary trend is accurately captured by the red damage region in the finite element simulation (right), where the damage shows a gradual expansion from the hole edge outwards. In the experiment, the hole edge region in the opposite direction of the tensile load was locally deformed by the bolt rod compression and formed a ‘fold buildup’. Finite element simulation predictions model this process. The red area shows the location of plastic deformation and damage, which is consistent with the experimentally observed trend.

The displacement–load curves (Figure 9) are labelled with three characteristic state points A, B and C. Figure 11 shows the corresponding experimental observations and numerical simulation results. Where SDV61 is fibre damage and SDV62 is matrix damage. The figure shows the failure characteristics of the material: at point A, the initial damage starts to appear at the edge of the hole; the load is increased to point B, the damage region expands and the material starts to deform; at point C, the structure exhibits material damage. The damage in the peripheral region of the hole continues to accumulate and leads to elongation of the hole shape, which is consistent with the experimentally detected failure process.

Figure 12 shows the comparison between the experimental scanning detection results and the numerical simulation results of the hole perimeter at point A. The contact pressure between the bolt and the plywood at point A resulted in the damage distribution in the surface region, which was manifested by the matrix rupture and the formation of shear cracks. The observation results show that the middle region of the hole remains relatively intact and does not show obvious damage characteristics.

Figure 13 shows the experimental observations and numerical simulation results of the hole circumference semi-truncated surface at point B state. The hole circumference region of the material in this state exhibits a characteristic damage pattern. Fibre kinking develops in parallel with matrix fracture, forming a characteristic wedge-shaped distribution. The model does not fail completely at this stage, and the hole structure still maintains basic integrity, which confirms that the joint still has the load-bearing capacity under the current load. The direction of fibre fracture is consistent with the expansion path of the wedge-shaped crack, and both loading directions form about 60°. Since the thin-layer configuration effectively inhibited the interlaminar stress diffusion, the thin-layer composites in this study did not show the delamination damage pattern of thick plate compression. Therefore, the damage suppression capability of the thin pre-preg enables the laminate to continue to withstand more damage accumulation before final failure.

Figure 14 shows the finite element prediction knot at point C. Figure 15 shows the damage distribution comparison of different layup angles in this state, and the comparison of the two presents the damage characteristics of the laminated plywood bearing surface at point C. Obvious material fragmentation and folds have been formed on the bearing surface, and the damage area is obviously extended along the hole circumference in the opposite direction of loading. It is driven by the localised stress concentration generated by the contact pressure between the bolt and the hole wall. When comparing the different ply orientations, the damage zone is more extensive in the 0°/90° ply than in the ±45° ply because of the mutual shear-tension coupling effect exhibited by the anisotropic composites under non-spindle loading conditions. The numerical analysis results show that the damage has been extended to the whole hole circumference region.

Therefore, the material damage process is affected by the fibre-laying angle. The Modified LaRC05 criterion successfully modelled the damage patterns of different ply layers under complex stress states.

Optimisation of BP neural network model architecture and hyperparameters

The architecture and hyperparameter configuration of the BP neural network model have an important impact on its prediction performance. However, there is currently no clear theoretical basis to guide the setting of these parameters, so parameter optimisation becomes a key step to improve the model performance. In this paper, a Bayesian optimisation method is adopted, with the MSE of the test data as the optimisation objective. The optimisation variables include the network structure and multiple hyperparameters. The optimisation parameters and their value ranges are defined as follows: number of neurons in hidden layer 1 (20∼30), difference in neuron count between hidden layers 2 and 1 (0∼6), difference between hidden layers 3 and 2 (0∼6), activation functions for hidden layers (Tansig, Logsig, ReLU), training functions (Trainlm, Trainbr, Trainbfg, Trainscg), and learning rate (0.0001∼0.2).

According to the related research in the literature,²⁸ the number of neurons in the first hidden layer should usually be more than that in the input layer, which is set between 20 and 30 in this paper. In order to improve the performance of the model, this paper adopts the pyramid structure design, i.e., the number of neurons in each hidden layer decreases layer by layer. The difference in the number of neurons between the first and second hidden layers is set as an optimisation variable, and a difference of 0 means that the second hidden layer is not enabled. Similarly, the third hidden layer is configured using this strategy, and its neuron difference with the second hidden layer is also used as an optimisation variable, and if the value is 0, it means that the third hidden layer does not exist. In terms of hyper-parameter settings, the activation function (e.g., ReLU, Sigmoid, and Tanh), the type of loss function (e.g., MSE), and the learning rate of each hidden layer are included in this paper for optimisation. The value range of the learning rate is set from 0.001 to 0.2, and a smaller learning rate helps to improve the stability and generalisation ability of the model.²⁹ The dataset used for model training contains a total of 600 samples.

The Bayesian optimisation process of the ANN model is shown in Figure 16. The y-axis is the minimum MSE value and x-axis is the number of iterations. the smaller the MSE, the better the prediction performance. It shows that as the number of iterations increases, the minimum MSE value of the model decreases and eventually stabilises, indicating that the Bayesian optimisation has found the optimal configuration.

The final optimised ANN model (ANN_MIU1) achieves a minimum MSE of 2.145 × 10⁻³ with the following configuration: hidden layer 1 has 23 neurons (Tansig activation), layer 2 has 19 neurons (Logsig), and layer 3 has 15 neurons (Logsig). The model uses Trainlm training function with a learning rate of 0.00214. Specifically, the ANN model reaches a minimum MSE value of 2.145E-3 at the 27th iteration, indicating that the model has a high prediction accuracy in this configuration. The following is the detailed configuration of the optimisation results:

Optimisation of dataset size for BP neural network models

In addition to the neural network architecture and hyperparameters, the dataset size also affects the prediction ability of the neural network model. Too small a dataset can lead to insufficient learning of the neural network, while too large a dataset not only significantly increases the computational cost, but also may cause the neural network to learn redundant information due to the duplication of samples or overlapping information. Therefore, it is necessary to improve the learning efficiency while ensuring the prediction ability of neural network efficiency. To this end, dataset scaling cases: Case 1 (600 samples), Case 2 (900), Case 3 (1200), Case 4 (1500), Case 5 (1800), this section designs in five sets of experiments, each with a different size of experimental dataset, and all of them use the LHS technique to obtain representative sample points in the input feature parameter space.

In order to eliminate the effect of data division randomness on the experimental results, a five-fold cross-validation scheme is implemented for each experimental group. Five-fold cross-validation randomly divides the dataset into five mutually exclusive subsets, and through five iterations of training, four subsets are selected as training data each time, and one subset is retained as the validation set. As can be seen from Figure 17, when the dataset size increases, the MSE mean (marked by red circles) based on the five-fold cross-validation shows a trend of rapid decay followed by a slight rebound. The minimum MSE value is 1.896 × 10^–3 at a dataset size of 1200. The blue diamond markers are the MSE values for the 5-fold cross-validation and the black boxes are the standard deviations of the MSE values for the 5-fold cross-validation. Also at a dataset size of 1200, the standard deviation reaches a minimum value of 3.647 × 10⁻⁴. This shows that the neural network not only has high prediction accuracy but also the neural network has strong stability.

Importance assessment of input features based on SHAP values

When constructing the neural network model, all the parameters in Table 2 were used as input features. It has been shown in the literature that too many input features do not necessarily improve the prediction ability of the model, and may even reduce the generalisation ability of the neural network due to redundant input features.³⁰ Therefore, this section introduces an interpretable analysis method based on game theory–Shapley Additive Explanations (SHAP) value theory, and constructs a feature importance evaluation system by quantifying the incremental marginal effect of each input feature on the prediction results. Figure 18 shows the result of SHAP value of input features, the horizontal axis is the input features, and the vertical axis is the evaluation absolute SHAP value.

From Figure 18, it shows that the distribution of the average absolute SHAP values of each input feature presents obvious differences, which indicates that different input features have different importance to the neural network model. $X_{T}$ (Transverse tensile strength), $X_{C}$ (Transverse compressive strength), and $E_{11}$ (Longitudinal modulus of elasticity) have the highest SHAP values, indicating that these mechanical parameters are highly correlated with the predicted results. $Y_{T}$ (Longitudinal tensile strength), $Y_{C}$ (Longitudinal compressive strength) have the next highest SHAP values, indicating that they are of medium importance to the neural network. The rest of the parameters have small SHAP values indicating that these input features have less impact on the model.

Improvement and evaluation of BP neural network model

In order to improve the prediction efficiency and generalisation ability of the BP neural network model, this paper systematically constructs eight groups of comparative experiments based on the feature screening criteria of SHAP value analysis, as shown in Table 3.The core difference of the experiments of each group is reflected in the optimisation strategy of the input feature space:

Table 3.

Eight groups with different input feature combinations.

Case I	X_T, X_C, E₁₁, Y_T, Y_C, S₁₂, E₂₂, h, l, W, γ, η_L, G₁₂, G₂₃, ν₁₂, ν₂₃, ηT
Case II	X_T, X_C, E₁₁, Y_T, Y_C, S₁₂, E₂₂, h, l, W, γ, η_L, G₁₂, G₂₃, ν₁₂, ν₂₃
Case III	X_T, X_C, E₁₁, Y_T, Y_C, S₁₂, E₂₂, h, l, W, γ, η_L, G₁₂, G₂₃, ν₁₂
Case IV	X_T, X_C, E₁₁, Y_T, Y_C, S₁₂, E₂₂, h, l, W, γ, η_L, G₁₂, G₂₃
Case V	X_T, X_C, E₁₁, Y_T, Y_C, S₁₂, E₂₂, h, l, W, γ, η_L, G₁₂
Case VI	X_T, X_C, E₁₁, Y_T, Y_C, S₁₂, E₂₂, h, l, W, γ, η_L
Case VII	X_T, X_C, E₁₁, Y_T, Y_C, S₁₂, E₂₂, h, l, w, γ
Case VIII	X_T, X_C, E₁₁, Y_T, Y_C, S₁₂, E₂₂, h, l, w

Table 3 shows the predicted results of the eight scenarios compared to the actual values. In order to visualise the difference in prediction performance of each case more, the corresponding MSE and coefficient of determination of the case are labelled in the lower right corner of each figure $R^{2}$ . From the Figure 19, Case VI has the best prediction performance with the smallest MSE and the largest $R^{2}$ values, which are 1.100E-4 and 0.996, respectively. It indicates that Case VI has obvious advantages in prediction accuracy. Comprehensive analysis shows that the improved BP neural network model is successfully obtained by removing the input features that have less influence on the model performance.

Finally, Table 4 summarises the results of the improved BP neural network model. These results further demonstrate that Case VI has superior performance in prediction. Meanwhile, the relevant architecture and hyperparameter configuration of Case VI. The optimised BP neural network features hidden layers with 14 (Tansig), 10 (Tansig), and 7 (Logsig) neurons respectively, using the Trainbr training function with a learning rate of 0.00894.

Table 4.

MES and $R^{2}$ values for 8 sets of cases.

	MSE	R ²
Case I	0.00202	0.94214
MSE	0.00119	0.966
Case III	7.1691E-4	0.97948
Case IV	5.9382E-4	0.983
Case V	2.5366E-4	0.99274
Case VI	1.1003E-4	0.99685
Case VII	2.0948E-4	0.994
Case VIII	1.3086E-4	0.99525

Conclusion

In this paper, a three-dimensional progressive damage model of composite/titanium alloy hybrid bolted joint structure was constructed based on the improved LaRC05 criterion and combined with the adaptive material degradation method. The numerical simulation of the model was realised through the secondary development of UMAT subroutine in conjunction with ABAQUS software, and the corresponding data set was established. Meanwhile, a fast prediction tool for ultimate load was developed based on neural network, which accomplished the efficient evaluation of metal hybrid connection structures.

The three-dimensional progressive damage model based on the improved LaRC05 criterion is programmed by the UMAT user material subroutine. By introducing in-situ effect, shear-driven effect and orthogonal anisotropy principal structure, the problem of insufficient accuracy of damage prediction by the traditional model for complex stresses is solved and the influence of the thickness of the lay-up on the mechanical properties of the plywood is analysed. A three-dimensional asymptotic damage prediction model is constructed by combining the damage evolution theory and a more accurate fracture angle search method.

Numerical simulation of damage prediction for TCT metal hybrid double lap structure based on the developed 3D asymptotic damage prediction model. Comparison with the conventional criterion and experimental results reveals the effect of shear driving on the damage evolution prediction. It is shown that the two guidelines have similar prediction effects in the initial loading stage, but after entering the nonlinear region, the improved guideline shows obvious advantages. The criterion can more accurately reflect the actual damage evolution process and avoid the problems of reduced abruptness and prediction lag of the traditional method. The improved method improves the inflection point load prediction accuracy by 5.25% and reduces the limit load prediction error to 1.6%, confirming its superior performance in predicting the failure behaviour under simulated complex stress states.

The BP neural network ultimate load prediction model is established by LHS, applying Bayesian optimisation automatic search and using five-fold cross-validation. The feature importance analysis based on SHAP values screened out the input features with explanatory power, simplified the model structure and improved the prediction accuracy. The study shows that the prediction accuracy and efficiency can be significantly improved by reasonably controlling the dataset size, selecting the input features and optimising the network configuration. It provides efficient and reliable technical support for the rapid evaluation of composite material ultimate load.

Footnotes

Acknowledgements

We express our sincere gratitude to these funding agencies for their support.

ORCID iD

Litai Sun

Funding

The authors disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This research was supported by the National Natural Science Foundation of China (Grant No. 51975449), the Shaanxi Key Research and Development Program (Grant No. 2024GX-YBXM-206) and the Xi'an ‘Scientists + Engineers’ Team Construction Project (Grant No. 24KGDW0026).

Declaration of conflicting interests

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

References

. Mechanical properties of CFRP plywood bolted joints with localized reinforcement. Xi'an: Xi'an University of Technology, 2024.

Guangzheng

Feng

Zhuoyu

, et al. Research progress on mechanical properties of composite bolted structures. New Technol New Process 2022; 9: 1–5.

Meijuan

Libin

. Progress in design and analysis of reinforcement in composite bolted joints. J Compos Mater 2023; 40: 3771–3784.

Hongchao

Yunhe

Zhenshan

, et al. Experimental study on fatigue performance of composite glued bolted hybrid nodes. Ind Build 2015; 45: 31–36.

Olivier

Csillag

Christoforidou

, et al. Feasibility of bolted connectors in hybrid FRP-steel structures. Constr Build Mater 2023; 383: 131100.

Kabche

Caccese

Berube

, et al. Structural response of a hybrid composite-to-metal bolted connection under uniform pressure loading. Compos Struct 2007; 78: 207–221.

Luo

Yang

Xuan

, et al. Experimental investigation of double-lap bolted composite joints with different outer plates subjected to bending loads. Polym Polym Compos 2018; 26: 408–419.

Kabche

Caccese

Berube

, et al. Structural response of a hybrid composite-to-metal bolted connection under uniform pressure loading. Compos Struct 2007; 78: 207–221.

Thite

Ghazvini

Wallace

, et al. Interfacial adsorption controls particle formation in antibody formulations subjected to extensional flows and hydrodynamic shear. J Pharm Sci 2023; 112: 2766–2777.

10.

Elhady

Sallam

. Crack sensitivity of bolted metallic and polymeric joints. Eng Fract Mech 2015; 147: 55–71.

11.

Kolesnikov

Herbeck

Fink

. CFRP/Titanium hybrid material for improving composite bolted joints. Compos Struct 2008; 83: 368–380.

12.

J-H

. Application of advanced materials in offshore drilling platform design. China Mach 2024; 2: 65–69.

13.

Akbarpour

Hallström

. Enhancing the performance of bolted joints in composites by use of patched steel or titanium inserts. Compos Struct 2021; 275: 114464.

14.

Yuwei

Xia

Jiyun

, et al. Research progress on friction and wear properties of titanium alloys and wear reduction methods. Surf Technol 2024; 53: 1–21.

15.

Kazi

Eljack

Mahdi

. Data-Driven modeling to predict the load vs. displacement curves of targeted composite materials for industry 4.0 and smart manufacturing. manufacturing. Compos Struct 2021; 258: 113207.

16.

Liu

Tian

Tao

, et al. A review of artificial neural networks in the constitutive modeling of composite materials. Compos B Eng 2021; 224: 109152.

17.

Wang

Liu

Zhou

, et al. Highly efficient simulation of composites by determining failure initiation and fracture angle with artificial neural networks. Compos Struct 2023; 307: 116644.

18.

Pinho

Iannucci

Robinson

. Physically-based failure models and criteria for laminated fibre-reinforced composites with emphasis on fibre kinking: part I: development. Compos Part A Appl Sci Manuf 2006; 37: 63–73.

19.

Pinho

Davila

Camanho

, et al. Failure models and criteria for FRP under in-plane or three-dimensional stress states including shear non-linearity . Nasa Technical Memorandum, 2005.

20.

Schirmaier

Weiland

Kärger

, et al. A new efficient and reliable algorithm to determine the fracture angle for puck's 3D matrix failure criterion for UD composites. Compos Sci Technol 2014; 100: 19–25.

21.

Biao

Yazhi

Ruixiang

. Failure analysis of laminates based on physical mechanisms. J Compos 2013; 30: 158–162.

22.

Sevostianov

Kachanov

. Explicit cross-property correlations for anisotropic two-phase composite materials. J Mech Phys Solids 2002; 50: 253–282.

23.

Matzenmiller

Gerlach

. Parameter identification of elastic interphase properties in fiber composites. Compos B Eng 2005; 37: 117–126.

24.

Hoffarth

Rajan

Goldberg

, et al. Implementation and validation of a three-dimensional plasticity-based deformation model for orthotropic composites. Compos Part A Appl Sci Manuf 2016; 91: 336–350.

25.

Cao

Zuo

Zhao

, et al. Experimental investigation on bearing behavior and failure mechanism of double-lap thin-ply composite bolted joints. Compos Struct 2021; 261: 113565.

26.

Yuejie

. Mechanical behavior and failure mechanism of carbon/epoxy composite bolted joints. Xi'an: Northwestern Polytechnical University, 2020.

27.

Yang

Bao

Liu

, et al. Progressive failure analysis of composite/aluminum riveted joints subjected to pull-through loading. Chin J Mech Eng 2023; 36: 129–138.

28.

Hui

Hongyu

Man

, et al. Analysis of research hotspots and theme evolution in machine learning based on scientometrics. J Shenyang Univ (Soc Sci Ed) 2021; 23: 19–25.

29.

Queiroz

Bueno

VFDS

Santos

HBD

, et al. Surface tension prediction of pure organic components: an artificial neural network approach. Fuel 2025; 380: 133096.

30.

Zhao

Shi

, et al. Developing a machine learning model for heat pipes considering different input features. Int J Therm Sci 2025; 208: 109398.