1.1 Global-local modelling

Digital design and validation tools are key cross-cutting technologies to enable the faster developmental cycles required to help the aviation sector transition to a sustainable future [1]. One challenge is to reduce the time required to predict failure at the design stage. For large composite structures, such as composite airframes, it is computationally unfeasible to use high-fidelity finite element analysis (FEA) to model failure across the entire structure. Therefore, it is typical to use a global-local modelling framework whereby low-fidelity FEA can determine the nodal displacements at the macroscale structural level. These can then drive the high-fidelity mesoscale models localised in the region of the design features of interest, to determine the detailed 3D stresses required to accurately predict damage initiation using a relevant failure criterion [Reference Ostergaard, Ibbotson, Le Roux and Prior2].

As there may be thousands or even millions of such features in large composite structures, this global-local modelling framework can become very time-consuming and repetitive. A large amount of expertise is required to generate these FEA models and the failure data generated for a given feature is not re-used. These problems can be addressed using machine learning (ML) within this preliminary design framework. Well-conditioned machine learning models can make fast predictions of complex problems, given appropriate training data, and have demonstrated success in a variety of fields.

ML has been used to accelerate multiscale finite element analyses of composite structures in literature. For example, ML has been used for the prediction of constitutive relationships [Reference Logarzo, Capuano and Rimoli3], acceleration of concurrent multiscale FE2 analyses [Reference Le, Yvonnet and He4, Reference Xu, Yang, Yan, Huang, Giunta, Belouettar, Zahrouni, Ben Zineb and Hu5], and to surrogate microscale behaviour via domain decomposition methods [Reference Krokos, Bui Xuan, Bordas, Young and Kerfriden6]. However, there is limited literature applying machine learning to global-local modelling, which involves the macro- and mesoscales. A data-driven method has been developed to accelerate coupled global-local analyses [Reference Reille, Champaney, Daim, Tourbier, Hascoet, Gonzalez, Cueto, Duval and Chinesta7]. In previous work, we have also applied ML to accelerate uncoupled global-local analyses [Reference Azeem and Iannucci8, Reference Azeem and Pinho9]. A benefit of this latter methodology is that local stress distributions can be evaluated, as required for subsequent failure analysis under the global-local framework.

The concept for this machine learning assisted global-local FEA framework was demonstrated for an open hole-in-plate design feature. By using high-fidelity FEA data to train a machine learning model, we demonstrated accurate stress predictions with an order of magnitude improvement in prediction speed.

In this work, we further develop this framework for a more commonly used and complicated design feature: the countersunk bolted composite joint. The developments in this paper contribute towards a proposed workflow given in Fig. 1, which aims to accelerate the development cycle of large composite structures.

Figure 1.

Proposed workflow using a ML-assisted preliminary design tool [Reference Azeem and Iannucci8].

1.2 Bolt failure

Countersunk bolted joints are still extensively used in composite airframes despite developments in adhesive joints. Bolted joints can fail in multiple modes, such as bearing failure, net-tension failure, shear-out failure and bolt failure among others [Reference McCarthy and McCarthy10]. However, at the preliminary design stage, bearing failure is of particular interest as this is a progressive form of failure that often precedes other modes [Reference McCarthy and McCarthy10]. Typically, engineers may employ methods such as the point stress criterion or the average stress criterion whereby the stress at a point a given length away or the averaged stress across a given length away from the bolt hole must not exceed a threshold stress level [Reference Whitney and Nuismer11]. This given length is referred to as the ‘characteristic distance’. The Yamada Sun stress criterion given in Equation (1) has been successfully used as the failure criterion to determine this threshold stress level, and in this case, bearing failure is based on fibre failure [Reference Zhang, Liu, Zhao, Chen and Fei12].

(1) \begin{align}Modified\ Yamada\ Sun\ Failure\ Index = \left\{\begin{array}{l@{\quad}l} \sqrt {{{\left( {{{{\sigma _{11}}} \over {{X_T}}}} \right)}^2} + {{\left( {{{{\sigma _{12}}} \over {{S_C}}}} \right)}^2}} & \left( {{\sigma _{11}} \gt 0} \right) \\[13pt] \sqrt {{{\left( {{{{\sigma _{11}}} \over {{X_C}}}} \right)}^2} + {{\left( {{{{\sigma _{12}}} \over {{S_C}}}} \right)}^2}} & \left( {{\sigma _{11}} \lt 0} \right) \end{array}\right.\end{align}

We therefore need to accurately predict the σ₁₁ and σ₁₂ stress components within a given characteristic distance from the bolt hole to be able to predict bearing failure. In this study, we use a volume-averaged stress criterion. Though there is a loss of information through volume averaging, namely the radial location of particular stresses, the average effect of stresses across the circumference is conserved. This, therefore, serves as a reduced output vector space to develop our baseline machine learning modelling strategies.

1.3 Previous work

The previous work [Reference Azeem and Iannucci8, Reference Azeem and Pinho9] on this tool had four key steps, see Fig. 2. These four steps were used to generate an ML surrogate of a high-fidelity open hole local model, to quickly predict ply-by-ply stress distributions given boundary conditions interrogated from a low-fidelity c-spar global model. We used Abaqus/2021 [13] for the finite element modelling and TensorFlow [Reference Abadi, Barham, Chen, Chen, Davis, Dean, Devin, Ghemawat, Irving, Isard, Kudlur, Levenberg, Monga, Moore, Murray, Steiner, Tucker, Vasudevan, Warden, Wicke, Yu and Brain14] to develop the machine learning models.

Figure 2.

Previous work on a ML-assisted tool to predict stresses around a hole-in-plate feature [Reference Azeem and Iannucci8].

First, we developed a sampling technique based on a modified Latin hypercube sampling scheme called Maximum Projection [Reference Joseph, Gul and Ba15]. This could generate a set of symmetric, balanced laminates which were well distributed with regard to the percentages of plies in each of the classic quad ply angles [0/45/-45/90]. Composite design rules [Reference Niu16] such as the 10% rule and the contiguity rule were considered when building laminates. The design feature parameter, i.e. the hole diameter, was also varied within this sampling scheme.

Secondly, degrees of freedom were placed on the four corners and four mid-side nodes of the laminates to reduce the dimensionality of the problem. This was achieved during the offline training process by applying these degrees of freedom to a coarse shell mesh that drove the high-fidelity local laminates in a submodelling procedure. In the online application, a bi-linear homogenisation scheme was developed to determine the magnitude of displacements in each degree of freedom given a global model simulation. This homogenisation scheme ensures that the work done on each edge is conserved, with a reduced number of nodal degrees of freedom on the boundary of the area of interest.

Thirdly, to generate the training data, a small displacement was applied in turn to each of the degrees of freedom along the boundary of the laminate. The underpinning idea is that under the principle of linear superposition, by predicting the stress response near the hole of the laminate due to each boundary perturbation, we could therefore understand the stress response of the laminate under an arbitrary boundary loading. The volume-averaged stress per ply at a given characteristic distance from the hole was determined using FEA.

The challenge of the machine learning model was therefore to predict the 3D volume averaged stresses per ply of an open hole laminate, with a given laminate stacking sequence and given hole diameter, under a given degree of freedom loading. So fourthly, a customised neural network machine learning model was used to achieve these predictions with high accuracy, as demonstrated for a case study of an open-hole feature in a composite C-spar. This neural network consists of a long-short term memory neural network with a bi-directional wrapper. Bi-directionality improved performance as the laminates in this study were symmetric.

1.4 Adaptations considered in this study

To adapt this methodology towards a countersunk bolted joint we must consider three additional steps. The method and results of each of these three adaptations are covered in this study.

Firstly, as we have two laminates bolted together, as well as the fact that we need to consider the effect of bolt pre-load, we re-consider which degrees of freedom we vary in our FEA simulation. Secondly, as we have multiple design feature geometry parameters, in this case, bolt radius and countersink depth, we need to engineer these features appropriately before passing them into the machine learning model. We compare two ways of representing these features in this study. Thirdly, we need to ensure the machine learning model results in high-accuracy predictions, and this requires further model hypertuning. For this, we investigate the use of neural network dropout, a common technique used to regularise networks and reduce overfitting.

2.1 Finite element analysis (FEA) modelling

Abaqus/2021 is used to model the bolted joints and run the general, static analyses. We model the bolt using aerospace grade 6AI-4V titanium (E = 110 GPa, ν = 0.29), and the nut with steel (E = 210 GPa, ν = 0.3). The material properties of the IM7/8552 carbon fibre composites used for the laminate are given in Table 1. The ply thickness is 0.125 mm. In-plane shear strength is derived using methods from literature [Reference Zhang, Liu, Zhao, Chen and Fei12].

Table 1.

Material properties of composite material [Reference Zhang, Liu, Zhao, Chen and Fei12, Reference Kevin O’Brien and Krueger17]

The nut and bolt are modelled as the same part but their respective sections are partitioned, see Fig. 3. For the bolted joint, we follow the modelling guidelines found in validated studies [Reference Egan, McCarthy, McCarthy and Frizzell18]. Linear, fully integrated solid elements (C3D8I) with fully incompatible modes are used, and one element is used per ply thickness. Extra circumferential partitions are created, which closely follow the profile of the countersunk bolt. Finite sliding surface-surface contact is employed, with the stiffer surfaces assigned the master role. The contact is modelled as frictionless, and we have zero bolt-hole clearance for this preliminary study, though the effect of these variables may be investigated in future work.

The length and width of the laminates are fixed to 32 mm. The model generation is scripted given the sample laminate stacking sequence and the bolt geometry parameters provided in Table 2. Seven samples did not result in converged finite element simulations. In total, 113 finite element simulations were performed in parallel, each taking ∼106 h to run. Note that the bolt countersink depth varies up to 0.7 times the laminate thickness, as suggested by bolt design guidelines [Reference Shyprykevich19].

Table 2.

Parameters varied in study and number of samples

Figure 3.

FEA modelling of bolted joint focused on meshing strategy.

In each simulation, for each step, we apply a nodal displacement of 0.1 mm on each degree of freedom, holding other degrees of freedom fixed. There are 24 degrees of freedom on each laminate, see Fig. 4, resulting in 48 steps. In a global-local context, these local boundary conditions would be extracted from a global shell model in the location of the local area of interest, where the bolted joint exists. The global shell model could be a large structural model of a c-spar or wing box, for example.

Figure 4.

Degrees of freedom used in study. Bottom and top laminate planes offset for diagrammatic purposes.

We apply an additional preliminary step to apply a 1kN pre-load on the bolted joint using Abaqus’ built-in bolt load tool. The effect of the pre-load can be scaled and superimposed for arbitrary pre-loads so linear superposition can be maintained, as long as care is taken to separate the effect of the pre-load from the effect of boundary loads on the stress field and as long as we are aware that only linear-elastic behaviour can be represented (i.e. no large displacements). For each step, the volume-averaged stress variants in the 11 and 12 directions are saved. The radial distance we averaged the volumetric stress per ply in this study is 0.5 mm from the hole boundary.

Note that we use linear-elastic analyses for our modelling, as according to the characteristic length method [Reference Whitney and Nuismer11] for predicting the failure of bolted joints. Therefore, non-linearities that are expected due to progressive damage are captured in the definition of the characteristic distance, and so are not considered in this study.

2.2 Machine learning (ML) modelling

Figure 5 shows an example of how the input data is represented for a given laminate. The maximum number of plies across bottom and top laminates in this study is 160. So, for a given laminate the input data consists of a (160,3) matrix to represent the laminate stacking sequence and bolt geometry parameters. Ply 0 represents the bottom ply. Rows beyond the number of plies for a given laminate are assigned a value of 0. In Fig. 5, for example, the first 28 rows have values for the 28-ply laminate to define the ply angle, hole radius and countersunk depth. However, the remaining 132 rows, not shown in the figure, would have a value of 0 in each of these columns. We investigated two methods of representing countersunk depth in this study, as seen in Fig. 5. Data structure A defines a repeating scalar value of the countersunk depth in millimeters. Data structure B defines an encoding of 0,1 and 2 for each ply to identify whether the ply belongs to the bottom laminate, non-countersunk region of the top laminate and countersunk region of the top laminate, respectively.

Figure 5.

Input data representations for an example laminate with [45/0/0/-45/90]s bottom laminate, [45/0/-45/90]_s top laminate and a 4 mm radius bolt with 0.5 mm countersunk depth.

For each laminate, the output data is a (160,49) matrix which contains the volume averages stress in a given direction, due to the 49 varied degrees of freedom. As before, volume average stresses for non-existent plies are given the value of 0. Given this input and output data, we can generate a baseline machine learning model to map from input to output. We generate train:validation:test sets using k-fold validation to monitor the performance of the machine learning model on held-out data. This nominally results in 67 training samples, 34 validation samples, and 12 test samples.

As the learning problem is sequence-sequence we use a neural network appropriate for sequential time-series data i.e. the Long Short-Term neural network (LSTM). We have shown previously that bi-directional LSTMs are especially accurate for predicting ply-by-ply stresses for symmetric layups [Reference Azeem and Iannucci8], as used in this study. The shape of our neural network is given in Fig. 6. Mean squared error is used as the loss metric but mean absolute error is monitored as well. Early stopping is used to terminate model training once the mean absolute validation error is no longer reducing at a significant rate (0.1/250 epochs).

Figure 6.

LSTM neural network architecture.