A machine learning assisted preliminary design methodology for bolted composite joints in large structures

Abstract

Damage initiation hotspots around features, such as bolts and ply drops, must be investigated during the preliminary design phase of large composite structures, such as composite airframes. A global-local modelling approach is commonly employed to perform this investigation, whereby a global low-fidelity model is used to drive high-fidelity local models around the features of interest. However, this methodology is slow, repetitive and expert-dependent. In this investigation, we address these issues by applying machine learning techniques to this global-local modelling framework and demonstrate the time-saving benefit when predicting damage initiation of bolted composite joints. Feature engineering of model inputs and outputs, and appropriate customisation of machine learning methods enables damage initiation prediction. Special consideration is given to the boundary conditions that must be varied to simulate the response of the bolted composite joints. Results show over three orders of magnitude time-saving benefit and satisfactory accuracy of the proposed methodology. This indicates its potential to be developed further into a rapid design and optimisation tool.

Nomenclature

${E_{11}}$

fibre direction Young’s modulus

${E_{22}}$ , ${E_{33}}$

transverse Young’s modulus

${G_{12}}$ , ${G_{13}}$

in-plane shear modulus

${G_{23}}$

out-of-plane shear modulus

${S_C}$

in-plane shear strength

${X_C}$

fibre direction compressive strength

${X_T}$

fibre direction tensile strength

${\upsilon _{12}}$ , ${\upsilon _{13}}$

in-plane Poisson’s ratio

${\upsilon _{23}}$

out-of-plane Poisson’s ratio

$\sigma $

stress

1.0 Introduction

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 [2].

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 [3], acceleration of concurrent multiscale FE2 analyses [4, 5], and to surrogate microscale behaviour via domain decomposition methods [6]. 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 [7]. In previous work, we have also applied ML to accelerate uncoupled global-local analyses [8, 9]. 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 [8].

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 [10]. 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 [10]. 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 [11]. 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 [12].

(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 [8, 9] 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 [14] to develop the machine learning models.

Figure 2. Previous work on a ML-assisted tool to predict stresses around a hole-in-plate feature [8].

First, we developed a sampling technique based on a modified Latin hypercube sampling scheme called Maximum Projection [15]. 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 [16] 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.0 Modelling methodology

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 [12].

Table 1. Material properties of composite material [12, 17]

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 [18]. 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 [19].

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 [11] 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 [8], 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.

3.0 Results

3.1 Improving the baseline ML model

The mean absolute error between volume averaged ${\sigma _{11}}$ ML predictions based on test data and the FEA simulation results averaged across three repeat tests was 9.84 + 2.52 MPa/mm³. Learning curves show the evolution of training and validation error over the training period and are given for the baseline test repeats in Fig. 7. This figure shows that the models consistently suffer overfitting, as the validation error exceeds the training error. This means that the models are overfitting to the training data only and thus not generalising to other bolted joints. The mean absolute validation error is 2.84 MPa/mm³ greater than the mean absolute training error.

Figure 7. Learning graphs for baseline ML model of all repeat tests (left) and best performing model test (right).

By using the trained ML model to make predictions for test data we can find further issues with the baseline model. Figure 8 shows the volume averaged stress per ply in the 11-direction for two different degrees of freedom as determined by FEA and as predicted by the baseline ML model. It is observed that the baseline model performs poorly at the transitions between the bottom and top laminate. To address this issue, we perform feature engineering to aid the ML network in deducing the difference between plies in the bottom laminate region, non-countersunk top laminate region and the countersunk top laminate region. We achieve this by creating an additional input feature labelled with the numbers 0,1 or 2 to classify to which region in the bolted joint each ply belongs.

Figure 8. Baseline ML model prediction vs FEA test data of ply-by-ply volume averaged stresses in 11-direction for two different degree of freedom loadings. Red region indicates error of concern.

The mean absolute error between ML predictions based on test data and the FEA simulation results averaged across three repeat tests for the feature engineered ML model was 3.82 + 0.07 MPa/mm³. This represents a 61% improvement on the baseline performance. The learning curves for this feature engineered ML model show much less overfitting, see Fig. 9. The mean absolute validation error is now 0.81 MPa/mm³ greater than the mean absolute training error. This represents a 71 % decrease in overfitting.

Figure 9. Learning graphs for model with additional feature engineering of all repeat tests (left) and best performing model test (right).

Figure 10 shows that the predictions made by this improved ML model follow much closer to the stresses determined from the FEA simulation. The ML model can determine the difference between behaviour in the top laminate and the bottom laminate much better, and this improvement has also led to the reduction of overfitting. However, the training error is still decreasing by the end of the training period epochs, whilst the validation error has converged. Therefore, there is still room to improve overfitting to improve ML model accuracy. To achieve this, we use neural network dropout, a common technique used to regularise networks and reduce overfitting. A small dropout of 0.1 is applied to the bi-directional LSTM layer.

Figure 10. Feature engineered ML model prediction vs FEA test data of ply-by-ply volume averaged stresses in 11-direction for two different degree of freedom loadings.

The mean absolute error between ML predictions based on test data and the FEA simulation results averaged across three repeat tests for the featured engineered ML model with additional dropout regularisation was 3.37 + 0.03 MPa/mm³ and the mean absolute validation error is now 0.35 MPa/mm³ greater than the mean absolute training error. The additional dropout has therefore resulted in a further 12 % error reduction and a 57 % reduction in overfitting.

However, when comparing the performance of the three ML models investigated thus far, as visible in Fig. 11, it is clear that the feature engineering delivers a much larger model improvement than the subsequent dropout regularisation. It is also observed that the lower the model error, the lower the standard deviation in model performance. Further model hypertuning may reduce the error further, however, as with the effect of dropout regularisation, this reduction may not result in drastic error reduction as overfitting is already close to zero.

Figure 11. Showing overall improvement of feature engineering and dropout.

3.2 Effect of amount of training data on ml model performance

A key parameter that may affect the model performance, that has not yet been varied, is the amount of training data. In the former investigations, the FEA-modelled behaviour of 67 bolted joints was used to train the ML models. However, by varying the k-fold partitions and the train-test split, we can investigate the performance of ML models on various amounts of training data. The feature engineered model without dropout was used for comparison. The results of this investigation, shown in Fig. 12, indicate two key points.

Figure 12. Effect of amount training data on model error.

Firstly, increasing the amount of training data is likely to improve performance further as errors haven’t converged. However, the likely improvement in ML model performance comes at a trade-off with the FE model simulation time required to generate this extra data. Secondly, the error with only 30 samples is still relatively low in the context of the range of volume averaged stresses in this study which reach up to 150 MPa/mm³. As standard deviation and overfitting increase with a lower amount of training data, there is scope for further model hypertuning to be used to reduce this error further. This means that only a low amount of training data (<50) may be necessary to achieve high prediction accuracy.

3.3 Testing ML model on loaded bolted joint

To further test the accuracy of our ML modelling we use the best developed model in this study, namely the feature engineered model with dropout regularisation, to predict the behaviour of a bolted joint under arbitrary boundary loading. We pre-load the bolts to 1kN. The FE model of the loaded bolted joint is shown in Fig. 13 and consists of a shear component and a tension component, to represent a mixed loading which may be typical in structural applications.

Figure 13. FE model with mixed boundary loading.

We use our ML model to predict the volume averaged stresses per ply due to each degree of freedom component of the given loading and superimpose the results to determine the total effect of the loading. As the Yamada Sun criterion can be used to evaluate bolt failure, and this criterion requires stresses in the 11 and 12 directions, we generate a machine learning model to predict each of these stress variants. The ply-by-ply variation of volume averaged stresses in the 11 and 12 directions are compared to the respective FEA model results in Fig. 14.

Figure 14. ML model predictions for key stress components as compared with FEA test data for bolted joint under given mixed boundary loading.

The stresses in the 11-direction are much higher and so are therefore the maximum errors in the machine learning model, as errors for a given unit loading scale with the applied loading. However, it is clear from Fig. 14, that for both stress components, the machine learning model successfully picks up the relative trends between stresses per ply. The magnitudes of the error may decrease upon further hypertuning of the machine learning model and the collection of further training data. However, it is important to note the time-saving benefit of this machine learning enhanced approach as, once we have fit the machine learning model, predictions take ∼200 ms using low computational resources (4CPU), whereas simulating each FE model takes 130 min using high computational resource (32CPU).

3.4 Further work on degrees of freedom

In this section, we observe issues in the application of the methodology for further representative loading conditions of bolted composite joints. We thereafter suggest and test a further degree of freedom which addresses this issue.

When testing the methodology to further loading conditions, it is found that tensile loading does not give the same magnitude of stress response as equivalent compressive loading. In a simple hole in plate case, this would be expected, however, bolt-hole interaction creates a difference between these two forms of loading. This can be addressed by doubling the number of degrees of freedom when generating training data, ensuring that a positive and negative displacement is applied to each of the x, y and z degrees of freedom at each boundary node. Furthermore, it is found that when the entire bottom and top laminate are displaced in relation to each other, not just partially as in Fig. 13, the stress responses of corresponding unit displacements do not linearly superimpose to give the same stress result, see Fig. 15.

Figure 15. Loading conditions (top) and ply-by-ply stress variation (bottom) demonstrating issue with linear superposition for unit relative movement between laminates.

This occurs as, for a given unit displacement on the boundary, the rest of the boundary is assumed to be fixed in position, and therefore tensile and compressive stresses are created by interaction with the bolt which remains in place. However, when the entire top laminate moves with respect to the bottom laminate, these tensile and compressive stresses are relieved as the bolt rotates. The contact stresses at the bolt-hole interface during this relative displacement are also not considered in the unit displacement-loading scenarios, and therefore they are not suitably considered during superposition.

However, linear superposition is typically considered for bolted joints within the bearing-bypass interaction. This interaction allows the superposition of stresses as a result of bearing loading, whereby force is reacted by the bolt, and bypass loading, whereby force is reacted by the laminate. Therefore, linear superposition should be attainable by modification of the degrees of freedom to account for bearing and bypass-related displacements.

Further work has demonstrated that this linear superposition can be achieved by adding an additional degree of freedom, via a relative unit displacement between bottom and top laminate. In this case, the displacements can be decomposed into a superposition of unit relative displacement between the top and bottom laminate, and unit displacements of boundary nodes for each of the bottom and top laminate, see Fig. 16. The relative displacement degree of freedom relates physically to bearing loading, whereby stresses arise due to bolt contact pressure. The displacement degrees of freedom relate physically to bypass loading, whereby laminate stresses arise regardless of bolt position.

Figure 16. Loading conditions (top) and ply-by-ply stress variation (bottom) demonstrating successful linear superposition with bearing and bypass-based displacements.

4.0 Conclusions

In this work, we have further developed a machine learning assisted framework for global-local stress analysis. This framework trains expert-gathered finite element analysis data of a given design feature, such as an open hole or bolted joint, in an offline training process. This offline training process can happen on high performance computing resources in a parallel manner to generate a database of feature-specific modelling data. The trained machine learning model can then be used much faster in situ to investigate stress and thereafter failure in local areas of interest in large structural models. This machine learning assisted framework is therefore useful to accelerate the virtual testing of composite structures during the design phase.

In this work, the proof of concept for this tool was extended, from application for open holes, to application for bolted composite joints through: consideration of the degrees of freedom for a bolted joint including pre-load; feature engineering to represent the through-thickness transitions between bottom laminate, top laminate and countersunk laminate section; and dropout regularisation to address model overfitting. These steps have enabled the tool to reach satisfactory accuracies for preliminary design.

Issues with respect to chosen degrees of freedom are addressed and future work aims to successfully incorporate this adjustment of degrees of freedom within the machine learning methodology. This work presents the stress distribution for a given characteristic distance. However, going forward, further work is required to determine the characteristic distance of a bolted joint given its laminate stacking sequence and bolt geometry parameters. This characteristic distance can then be used in conjunction with the developed ML framework to flag failure when damage has accumulated to this distance. This must be developed with respect to point stress or average stress criterion. The application of this tool must also be tested for other design features such as ply drops and ply cuts.

In this study, the offline training phase was performed using a high-performance computing environment in one week, and enabled online predictions which take 200 ms as opposed to 7800 s. This methodology is therefore of primary interest for repetitive design features, where the increase in offline training times may be accepted for a drastic reduction in online prediction times of the machine learning model. This tool is also especially useful for design optimisation studies, where potentially thousands of objective function predictions are required. Such studies can be better enabled using computationally inexpensive machine learning models as opposed to computationally expensive high-fidelity finite element models.