Prediction of shape distortions in thermosetting composite parts using neural network interfaced visco-elastic constitutive model

Abstract

The work aims to enhance the capabilities of a Finite Element tool, specifically related to a rheological thermo-chemo-viscoelastic constitutive model. This enhancement is intended to improve the tool’s ability to predict the distortions in composite parts caused by the polymerization of the thermoset composite matrix. These distortions occur due to internal residual stress generated by the inherent anisotropic properties of the thermoset composite material, including coefficients of thermal expansion and chemical shrinkage. The research work’s improvement is tied to the precise modelling of curing behaviour, which literature acknowledges as having a significant impact on manufacturing defects. In order to accommodate the influence of curing behaviour on various process variables—specifically, different thermal loading rates—a neural network model is implemented as an alternative to a standard diffusion cure-kinetics model. The neural network model is trained using Differential Scanning Calorimetry data and is integrated with the classical visco-elastic constitutive model to more accurately predict the progression of distinct thermoset resin states. This transition between cure states is assessed using two cure state variables: the degree of cure and the glass transition temperature. The enhanced predictions of state transitions lead to accurate assessments of internal residual stresses, especially when dealing with thick components subjected to thermal fluctuations. The anisotropic properties of thermoset composites, crucial for numerical analysis, are captured at various stages of cure. Ultimately, this methodology is employed to compare process-induced defects in the case study of the Z-shaped carbon/epoxy woven part, and the defects closely align with experimental measurements.

Research Highlights

• Enhancement of a Finite Element tool that utilizes rheological thermo-chemo-viscoelastic constitutive models.

• The cure state variables are observed to be dependent on thermal loading variables, even for cases below 3°C/min.

• The cure model parameters must be treated stochastically to account for uncertainties in thermal loading conditions.

• A non-parametric neural network model has been developed using characterization tests under various thermal loading conditions.

• The neural network-enhanced constitutive model outperforms the standard cure kinetics-interfaced constitutive model in accurately predicting distortions for the Z-shaped case study.

Introduction

Thermosetting composite materials have found extensive use in the aerospace industries due to their high strength-to-weight ratio characteristics. It has been observed that numerous iterative analyses, encompassing aspects like composite materials (resin and fibre) selection and optimization, mould characteristics, etc., are necessary to establish an efficient manufacturing process. This is primarily due to the defects encountered during the manufacturing stage.^1,2 During the manufacturing phase, composite parts acquire distorted shapes and internal residual stresses while being processed within the mould. These deformations are a consequence of external forces and thermal gradients, and will result in escalated costs during the final assembly of composite structures.^3,4 Moreover, studies^5–8 indicate that the developed internal residual stress fields have a significant impact on damage propagation, specifically delamination and transverse cracks. Therefore, it is imperative to develop novel approaches that can efficiently and accurately predict and validate internal residual stress fields in composite parts and structures.

Background - manufacturing distortions

The formation of manufacturing defects resulting from internal residual stress is attributed to various factors operating at different temperature levels.^4,9 In fact, during the polymerization reaction, also known as curing process, the thermosetting resin makes a transition from viscous to rubbery and finally to the glassy state, and the composite material experiences internal residual stresses and deformations related to the thermal and chemical behaviour of the resin. ¹⁰

• Thermally induced stress fields arise from the disparity in coefficients of thermal expansion (CTE) between the resin and the fibre, as discussed in the works.^2,11,12 When subjected to heating, the resin within the matrix expands due to its higher CTE, but this expansion is constrained by the fibres, resulting in compressive residual stresses. This phenomenon is numerically demonstrated in study. ¹³ Furthermore, the CTE of thermoset composites can exhibit variations at different stages of the manufacturing process. Hence, meticulous consideration of this property is essential to mitigate distortions in the final part.

• Another source of residual stresses arises from chemical shrinkage, which is linked to coefficients of chemical shrinkage (CCS) and involves the reduction of free volume within the resin during the curing process. The resin molecules are pulled inwards throughout the end linking reactions at higher temperatures, leading to transition between different states. Consequently, the available resin free volume decreases, exerting pressure on the fibres in the form of tensile residual stresses, which are uniformly distributed over the entire composite laminate. The strain associated with chemical shrinkage is assumed to follow a linear relationship with the cure state variable, as demonstrated in the study. ¹⁴

• Furthermore, deformations are constrained during curing by the contact between the material and the mould, resulting in the development of additional internal residual stresses.

Upon de-moulding, the stresses at the free surface are relieved, leading to the distortion of the manufactured composite parts. The residual stress field that persists within the component inevitably reduces its load-bearing capacity as it approaches the nominal strength of the resin and, consequently, contributes to the development of micro-cracks. It is also crucial to understand that other factors, such as lay-up configuration, mould material, composite part thickness, curing temperature etc., have a profound influence on the distorted shape.^15–19 At macro-scale, residual stresses in the laminates arises due to differences in thermal expansion coefficients, elastic moduli, and fibre orientations between the layers. ²⁰ Furthermore, the coefficient of friction between the mould and the composite part has been observed to have a significant impact on the final shape. ²¹ Additionally, studies^22–24 show that the moisture absorption, resin age and porosity can contribute to residual stresses, although, for simplicity, they have been neglected for the current study.

Existing methods to mitigate manufacturing defects

Several methods exist to optimize the manufacturing process with the goal of minimizing residual stresses or incorporating compensating features in the moulds used for composite part production. However, in the industry, trial-and-error approaches for predicting curing distortions are both costly and time-consuming. The introduction of analytical methods, relying on curing state variables and physical properties at different stages, has limitations and has been observed to be specific to particular applications.^25–27 Furthermore, such analytical approaches have yielded inconclusive results concerning the impact of thickness and shear stiffness during curing in the through-thickness direction. Aiming to overcome this limitation, Wisnom et al. ²⁸ introduced a concept that considers the influence of low shear stiffness in the rubbery state on process-induced defects. In a different study, Wucher et al. ²⁹ characterized distortion modes and optimized control points to determine an optimal mould configuration, yet this approach is specific to particular applications. Regrettably, the approaches discussed above lack reliability for complex geometries, such as L-shaped or Z-shaped parts, where flat sections tend to suppress shear deformation, as pointed out by Kappel et al. ¹⁸ Recently, modified analytical solutions have been developed based on the work of Wisnom et al. ²⁸ to address curved thermosetting parts, including cases like L-shaped parts with consideration for the flange effect. These solutions have demonstrated results comparable to Finite Element (FE) numerical analyses, as shown in studies.^30–33 However, the necessity to incorporate residual stress fields in addition to distortion predictions for further damage progression and failure analysis emphasizes the demand for constitutive models.

Finite Element Analysis (FEA) coupled with mechanical constitutive models, on the other hand, has proven to be highly effective for in-depth analysis, capturing various mechanisms during the manufacturing process. Numerous mechanical constitutive models, such as elasticity and viscoelasticity, have been implemented over the years, as referenced in Ding et al. ³⁴ Although full viscoelastic models^13,35–38 accurately describe curing behaviour, they demand extensive material characterization, leading to high computational costs. ⁹ In contrast, the Cure Hardening Instantaneously Linear Elastic (CHILE) model introduces less complex incremental elastic relationships, with the elastic modulus changing as a function of the cure state variables, as elaborated in studies.^14,39–42 A classical instance of the CHILE model was presented by Svanberg and Holmberg, ⁴³ wherein a simplified viscoelastic model incorporating path dependence on cure state variables is implemented.

Although constitutive models have been established, the modelling of cure state evolution remains a significant challenge due to uncertainties inherent in the manufacturing process, as highlighted by Mesogitis et al. ⁴⁴ and has received limited attention. The cure evolution is influenced by factors such as resin thermal history, resin formulations and characterization uncertainties.^45–47 In recent times, stochastic numerical analyses based on machine learning approaches have been employed to address these uncertainties in manufacturing. For instance, Mesogitis et al. ⁴⁸ demonstrated the quantification of variability in fibre orientation and its impact on maximum residual stress. It was also noted that a substantial portion of cure state uncertainty can be captured by specific parameters in the cure model. ⁴⁹ Moreover, the modelling of cure kinetics for a thick thermoset part subjected to temperature gradients becomes more challenging when using a deterministic parametric cure-kinetics model. ⁵⁰ The majority of such optimization procedures are combined with numerical analyses to accurately estimate these parameters and model optimal curing cycles that result in fewer process-induced defects. This approach is evident in works.^51,52 Alternatively, these model parameters can also be directly estimated from the DSC characterization tests using regression methods.^53,54

However, the cure-kinetics model does not account for different process condition variables. Consequently, there is the need to address the model parameters stochastically to encompass various thermal loading conditions. Regrettably, the absence of prior knowledge about model parameters for different resins presents challenges for stochastic simulations. Moreover, if the resin is expected to undergo substantial changes in its chemical composition, more intricate cure kinetics formulations are essential.^55,56 These complex variations of the Arrhenius model introduce additional parameters, making the regression procedure for their determination laborious, especially with the available DSC characterization tests. Consequently, achieving an accurate representation of the cure kinetics of the thermoset material is crucial in minimizing variability in the manufacturing process outcome.

In this work, we propose interfacing a neural network model with the user material subroutine of the FE software ABAQUS to enhance the numerical tool for predicting process-induced distortions and internal residual stresses. The performance of the numerical analysis depends on the accurate modelling of the material behaviour during the curing process. To alleviate the burden associated with the choice of cure kinetics model and its optimized parameters, a neural network trained with Differential Scanning Calorimetry (DSC) experimental datasets is used to model the cure behaviour, which the viscoelastic constitutive model depends on. Section Cure kinetics and material models of the paper describes the cure kinetics models and the FE-based rheological thermo-chemo-viscoelastic constitutive model within the study. In Section Characterization tests and neural network model, we describe the accurate modelling of cure behaviour using a neural network model trained on the experimental dataset of epoxy carbon fibre characterization. Section Results and discussion presents the implementation on the case studies of a flat plate and Z-shaped thermoset parts cured in autoclave, specifically comparing the process-induced defects obtained by using two different constitutive models in simulation: one based on a standard diffusion cure kinetics model and the other on a neural network model. Finally, the conclusions and remarks from the study are summarized in Section Summary and conclusions.

Cure kinetics and material models

Cure-kinetics and associated state variables

The properties relevant to curing and the modelling of residual stresses include the degree of cure, the instantaneous glass transition temperature, the volumetric free strains, the CTEs, the CCSs, and the composite laminate properties. The degree of cure, denoted as X, represents the ratio of heat released at a given time, denoted as t, to the total amount of heat released during the curing process, q(t). The critical points in the cure evolution are X_gel (gelation) and X_vitr (vitrification) which define the transition between viscous and rubbery states and between rubbery and glassy states respectively. The temperature cycle takes the form of a one-hold isothermal dwell associated with temperature rates, denoted as r_i. The value of X ranges between 0 (indicating a fully uncured resin) and 1 (indicating a fully cured resin), and it is determined by,

$\mathrm{X}\left(\mathrm{t}\right)=\frac{\mathrm{H}\left(\mathrm{t}\right)}{{\mathrm{H}}_{\text{overall}}}\hspace{0.17em}\text{where}\hspace{0.17em}\mathrm{H}\left(\mathrm{t}\right)={\int }_0^{\mathrm{t}}\left(\frac{\partial \mathrm{q}}{\partial \mathrm{t}}\right)d\mathrm{t}$

(1)

where ∂q/∂t represents the heat generation rate, H(t) is heat distribution resulting from the polymerization reaction during the curing process, and H_overall is the total enthalpy of the completely cured resin. Several phenomenological models have been developed to describe the progression of cure and establish correlations with the experimental measures obtained from DSC. These models adhere to a generalized form of the Arrhenius-type equation. ⁵⁷ Specifically, the cure reaction rate, dX/dt is described as a function of the degree of cure, X and temperature, T as follows:

$\begin{array}{l}\frac{d\mathrm{X}}{d\mathrm{t}}=\left({\mathrm{K}}_1+{\mathrm{K}}_2{\mathrm{X}}^{\mathrm{m}}\right){\left(1-\mathrm{X}\right)}^{\mathrm{n}}\\ \text{where} {\mathrm{K}}_{\mathrm{i}}={\mathrm{A}}_{\mathrm{i}}{e}^{-\mathrm{\Delta }{\mathrm{E}}_{\mathrm{i}}/\text{RT}}\end{array}$

(2)

where K_i represents the rate constants, R is the molar gas constant, A_i denotes the exponential law coefficients and; m and n are the first and second exponents. Finally, ΔE_i represents the activation energies. Following the vitrification stage during curing, there is a notable reduction in molecular mobility due to the decrease in resin volume. Consequently, it becomes imperative to control the cure rate. The transition from a rubbery to a glassy state is known to occur when the process temperature aligns with the instantaneous glass transition temperature, denoted as T_g. This variable is dependent on X, and is defined by Di Benedetto’s formula, ⁵⁸

${\mathrm{T}}_{\mathrm{g}}=\frac{\lambda \mathrm{X}}{1-\left(1-\lambda \right)\mathrm{X}}\left({\mathrm{T}}_{\text{gv}}-{\mathrm{T}}_{\text{gg}}\right)+{\mathrm{T}}_{\mathrm{g}0}$

(3)

where T_gg and T_gv represent the transition temperatures at the gel and vitrification points, respectively, T_g0 stands for the T_g at the initial time, and λ is a material-dependent parameter used to fit the equation to experimental data. To accurately model the curing behaviour, dX/dt has been individually solved for each state, as outlined in Sun et al. ⁵⁹ Nevertheless, this method demands the incorporation of a larger number of constants to model the cure kinetics. Alternatively, a more efficient cure kinetics model that incorporates this phenomenon through diffusion control has been introduced by Cole et al. ⁶⁰ The expression is a modified variant of equation (2):

$\begin{array}{l}\frac{d\mathrm{X}}{d\mathrm{t}}=\frac{{\text{KX}}^{\mathrm{m}}{\left(1-\mathrm{X}\right)}^{\mathrm{n}}}{1+e^{\mathrm{c}\left(\mathrm{X}-{\mathrm{X}}_{\mathrm{c}}\right)}}\\ \text{where}\hspace{0.17em}{\mathrm{X}}_{\mathrm{c}}={\iota }_0-\left({\iota }_t*\mathrm{T}\right)\end{array}$

(4)

where c represents the parameter controlling the breadth of the transition of the reaction from chemical to diffusion control, ι₀ stands for the X at the initial time and ι _t is a constant accounting for the increase of critical X with temperature. When the concentration of X reaches a certain critical value, diffusion control takes over the cure kinetics. The constants utilized in equations (3)–(4) are subject to optimization when employed for the curing process under various thermal loading conditions.

Additionally, the temperature across the thickness of the part significantly deviates from the control temperature set in the oven due to several factors. The part lags behind the control temperature during the heat-up process, influenced by X, thermal properties, and convective thermal resistance. The heat transfer in the thickness direction, z, can be described as follows:

$\zeta {\mathrm{c}}_{\mathrm{p}}\frac{\partial \mathrm{T}}{\partial \mathrm{t}}=\frac{\partial }{\partial \mathrm{z}}\left[\mathrm{L}\frac{\partial \mathrm{T}}{\partial \mathrm{z}}\right]+{\zeta }_{\mathrm{m}}{\mathrm{H}}_{\text{overall}}\frac{d\mathrm{X}}{d\mathrm{t}}\left(1-{\mathrm{V}}_{\mathrm{f}}\right)$

(5)

where ζ stands for laminate density, ζ_m corresponds to the resin matrix density, c_p signifies the specific heat capacity and L denotes thermal conductivity in the thickness direction. Meanwhile, V_f stands for the fibre volume fraction. Accurate temperature predictions for the part can be achieved through such thermal analysis, which incorporates thermo-chemical effects and process condition specifications.^37,39,61 In thick laminates, the temperature gradients can be more pronounced due to the lower thermal conductivity of the material, which can result in larger temperature differentials across the material thickness. ⁶² These gradients can result in regions of the laminate that are over-cured, under-cured, or have incomplete curing, which can lead to localized stresses and reduced mechanical properties.

Viscoelastic material model

The state-of-the-art solution for predicting the distorted shapes of composites after manufacturing relies on FE tools interfaced with thermo-viscoelastic constitutive models. The prevalent material behaviours employed to analyse defects encompass the CHILE model and the non-linear viscoelastic model. In the simplified case of the CHILE model, as proposed by Adolf and Martin ⁶³ and further developed by Svanberg and Holmberg, ⁴³ a path dependence on state variables, specifically volumetric free strains (ɛ^vf), X and T_g, is incorporated. The volumetric free strains, resulting from thermal expansion (ɛ^th) and chemical shrinkage (ɛ^ch), are defined as follows:

${\epsilon }^{\text{vf}}={\epsilon }^{\text{th}}+{\epsilon }^{\text{ch}}$

(6)

${\epsilon }^{\text{th}}={\int }_0^{\mathrm{t}}\alpha \left(\mathrm{X},\mathrm{T}\right)\frac{\partial \mathrm{T}}{\partial \mathrm{t}}d\mathrm{t} \text{and} {\epsilon }^{\text{ch}}={\int }_0^{\mathrm{t}}\beta \left(\mathrm{X},\mathrm{T}\right)\frac{\partial \mathrm{T}}{\partial \mathrm{t}}d\mathrm{t}$

(7)

$\text{where} \alpha | \beta =\{\begin{array}{l}{\alpha }^{\mathrm{v}} | {\beta }^{\mathrm{v}}\forall \mathrm{X}<{\mathrm{X}}_{\text{gel}} \text{and} \mathrm{T}\ge {\mathrm{T}}_{\mathrm{g}}\left(\mathrm{X}\right)\\ {\alpha }^{\mathrm{r}} | {\beta }^{\mathrm{r}}\forall \mathrm{X}\ge {\mathrm{X}}_{\text{gel}} \text{and} \mathrm{T}\ge {\mathrm{T}}_{\mathrm{g}}\left(\mathrm{X}\right)\\ {\alpha }^{\mathrm{g}} | {\beta }^{\mathrm{g}}\forall \mathrm{T}<{\mathrm{T}}_{\mathrm{g}}\left(\mathrm{X}\right)\end{array}$

where αⁱ and βⁱ are the CTEs and CCSs in the different states during curing process. Upon gelation, it is observed that the effect of the reduction of free volume within the resin in the form of chemical shrinkage, coupled with matrix expansion in the form of thermal expansion, occurs. Furthermore, upon reaching the glassy state, i.e., when reaching the end of the cross-linking reaction, the effect of chemical shrinkage is reduced, coupled with no thermal effect due to the isothermal temperature hold. The thermal effect in the glassy state is only triggered upon reaching the cooling phase of the cure cycle.

Furthermore, the transversely isotropic mechanical characteristics, such as the elastic modulus (E), shear modulus (G) and Poisson’s ratio (μ) also vary in different states. Consequently, multiple characterization tests, including DSC, Thermal Mechanical Analysis (TMA), and Dynamic Mechanical Analysis (DMA), are necessary for the implementation of constitutive models.^64,65

A model describing the CHILE behaviour or the non-linear viscoelastic behaviour of partially cured resin takes into consideration the process conditions and cure state variables. The instantaneous relaxation modulus, corresponding to Maxwell’s element with a spring and dashpot system. The viscoelastic model of the partially cured resin, depends on the reduced time associated with the spring and dashpot, denoted by ψ and ψ′, respectively. These reduced times are dependent on the X and T_g. The thermo-viscoelasticity for anisotropic and thermo-chemo-rheological materials is expressed in integral form, where it is a function of the instantaneous relaxation modulus, C. The equation is given as:

$\sigma =\mathbf{C}\epsilon +{\int }_0^{\mathrm{t}}\delta \mathbf{C}\left(\psi -{\psi }^{\prime }\right)\frac{\partial \left(\epsilon -{\epsilon }^{\text{vf}}\right)}{\partial \tau }d\tau$

(8)

where C represents the linear elastic Hooke tensor, δC is the independent relaxation function of N Maxwell elements, C _∞ is the fully relaxed modulus of uncured resin, C _n is the spring constant and ρⁱ are the relaxation times, with i denoting the Maxwell element. The dependence of ψ and ψ′ on cure state variables is established using a shift factor, a_T, in the constitutive model.

$\psi ={\int }_0^{\mathrm{t}}\frac{1}{{\mathrm{a}}_{\mathrm{T}}}d{\mathrm{t}}^{\prime };\hspace{1em}{\psi }^{\prime }={\int }_0^{\mathrm{t}}\frac{1}{{\mathrm{a}}_{\mathrm{T}}}d{\tau }^{\prime }$

(9)

The generalized Maxwell’s model with N elements, provides the relaxation Hooke’s tensor as follows,

$\mathbf{C}=\{\begin{array}{l}0,\mathrm{X}<{\mathrm{X}}_{\text{gel}}\\ {\mathbf{C}}_{\infty }+{\displaystyle \sum }_{i=1}^{\mathrm{N}}{\mathbf{C}}_{\mathrm{i}}\left(e^{\left(\frac{-\mathrm{t}}{{\rho }^{\mathrm{i}}}\right)}\right),\mathrm{X}\ge {\mathrm{X}}_{\text{gel}}\end{array}$

(10)

The shift factor, a_T is assumed to be a function that approaches 0 in the rubbery state and infinity in the glassy state. A variable γ is introduced to define the functional approximation of a_T, which is given by:

${\mathrm{a}}_{\mathrm{T}}=\underset{\gamma \to 0}{\lim }\{\begin{array}{l}\gamma ,\mathrm{T}\ge {\mathrm{T}}_{\mathrm{g}}\left(\mathrm{X}\right)\\ \frac{1}{\gamma },\mathrm{T}<{\mathrm{T}}_{\mathrm{g}}\left(\mathrm{X}\right)\end{array}$

(11)

This assumption governs the behaviour of the Hooke tensor: in the rubbery state, C^r is denoted as the fully relaxed Hooke tensor, C _∞ while in the glassy state, C^g is represented as $\left({\mathbf{C}}_{\infty }+{\sum }_{\mathrm{i}=1}^{\mathrm{N}}{\mathbf{C}}_{\mathrm{i}}\right)$ . The stress increment, Δσ corresponding to incremental time step is approximated by,

$\mathrm{\Delta }\sigma =\left({\mathbf{C}}_{\infty }+\frac{1}{\mathrm{\Delta }\psi }{\displaystyle \sum }_{i=1}^{\mathrm{N}}{\rho }^{\mathrm{i}}{\mathbf{C}}_{\mathrm{i}}\left(e^{\left(\frac{-\mathrm{\Delta }\psi }{{\rho }^{\mathrm{i}}}\right)}\right)\right)\mathrm{\Delta }\left(\epsilon -{\epsilon }^{\text{vf}}\right)+\mathrm{\Delta }{\sigma }^{\mathrm{R}}$

(12)

The term, Δσ^R is given by,

$\mathrm{\Delta }{\sigma }^{\mathrm{R}}=-{\displaystyle \sum }_{i=1}^{\mathrm{N}}\left(1-e^{\left(\frac{-\mathrm{\Delta }\psi }{{\rho }^{\mathrm{i}}}\right)}\right){\mathbf{S}}^{\mathrm{i}}\left(\mathrm{t}\right)$

(13)

where S represents the internal residual stress field, initially defined as zero, defined by a recursive relation,

${\mathbf{S}}^{\mathrm{i}}\left(\mathrm{t}+\mathrm{\Delta }\mathrm{t}\right)=\left(1-e^{\left(\frac{-\mathrm{\Delta }\psi }{{\rho }^{\mathrm{i}}}\right)}\right){\mathbf{S}}^{\mathrm{i}}\left(\mathrm{t}\right)+{\rho }^{\mathrm{i}}\frac{\mathrm{\Delta }\left(\epsilon -{\epsilon }^{\text{vf}}\right)}{\mathrm{\Delta }\psi }{\mathbf{C}}_{\mathrm{i}}\left(e^{\left(\frac{-\mathrm{\Delta }\psi }{{\rho }^{\mathrm{i}}}\right)}\right)$

(14)

This is incorporated into the FE-tool following the approach outlined by Zocher et al. ⁶⁶ The reduced time is given by, Δψ = Δt/a_T. Inserting equation (11) and upon simplifying yields,

$\mathrm{\Delta }{\sigma }^{\mathrm{R}}=\{\begin{array}{l}-{\mathbf{S}}^{\mathrm{j}}\left(\mathrm{t}\right),\forall \mathrm{T}\ge {\mathrm{T}}_{\mathrm{g}}\left(\mathrm{X}\right)\\ 0,\forall \mathrm{T}<{\mathrm{T}}_{\mathrm{g}}\left(\mathrm{X}\right)\end{array}$

(15)

${\mathbf{S}}^{\mathrm{j}}\left(\mathrm{t}+\mathrm{\Delta }\mathrm{t}\right)=\{\begin{array}{l}0,\forall \hspace{0.17em}\mathrm{T}\ge {\mathrm{T}}_{\mathrm{g}}\left(\mathrm{X}\right)\\ {\mathbf{S}}^{\mathrm{j}}\left(\mathrm{t}\right)+{\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{N}}{\rho }^{\mathrm{i}}{\mathbf{C}}_{\mathrm{i}}\mathrm{\Delta }\left(\epsilon -{\epsilon }^{\text{vf}}\right),\forall \hspace{0.17em}\mathrm{T}<{\mathrm{T}}_{\mathrm{g}}\left(\mathrm{X}\right)\end{array}$

(16)

where S^j(t) is defined by ${\sum }_{\mathrm{i}=1}^{\mathrm{N}}{\rho }^{\mathrm{i}}{\mathbf{S}}^{\mathrm{i}\left(\mathrm{t}\right)}$ . This results in stress increment with dependency on cure state variables corresponding to a time step, Δt is given by,

$\mathrm{\Delta }\sigma =\{\begin{array}{l}{\mathbf{C}}^r\mathrm{\Delta }\left(\epsilon -{\epsilon }^{\text{vf}}\right)-{\mathbf{S}}^{\mathrm{j}}\left(\mathrm{t}\right),\forall \mathrm{T}\ge {\mathrm{T}}_{\mathrm{g}}\left(\mathrm{X}\right)\\ {\mathbf{C}}^g\mathrm{\Delta }\left(\epsilon -{\epsilon }^{\text{vf}}\right),\forall \mathrm{T}<{\mathrm{T}}_{\mathrm{g}}\left(\mathrm{X}\right)\end{array}$

(17)

The internal stress field, S^j evolves based on the loading history and is calculated as,

${\mathbf{S}}^{\mathrm{j}}\left(\mathrm{t}+\mathrm{\Delta }\mathrm{t}\right)=\{\begin{array}{l}0,\forall \mathrm{T}\ge {\mathrm{T}}_{\mathrm{g}}\left(\mathrm{X}\right)\\ {\mathbf{S}}^{\mathrm{j}}\left(\mathrm{t}\right)+\left({\mathbf{C}}^{\mathrm{g}}-{\mathbf{C}}^{\mathrm{r}}\right)\mathrm{\Delta }\left(\epsilon -{\epsilon }^{\text{vf}}\right),\forall \mathrm{T}<{\mathrm{T}}_{\mathrm{g}}\left(\mathrm{X}\right)\end{array}$

(18)

In the rubbery and viscous states, the material is assumed to be fully relaxed, resulting in S^j being assigned a value of 0. In contrast, during the glassy state, the S^j is progressively accumulated. This internal residual stress is ultimately released upon the completion of the curing process, resulting in distorted shapes. The subroutine HETVAL is used to perform Fourier heat conduction analysis (refer to equation (5)), which incorporates thermo-chemical effects due to thermal loading specifications. This yields the actual part temperature, which is then used as input for the CHILE model to predict post-manufacturing PIDs.

Characterization tests and neural network model

The material of interest is a unidirectional carbon/epoxy prepreg, specifically AS4/8552. The DSC test is conducted to determine the rate of polymerization, denoted as X as a function of time, t for a given temperature curing cycle. This test measures the disparity in heat flow between capsules containing the resin sample under analysis and a reference resin sample. A sample of 5 milligrams, is subjected to a curing cycle, and the variables corresponding to exothermic and endothermic phenomena, as outlined in equation (1), are recorded. The former signifies the heat released by the sample during the curing process, H_t, whereas the latter represents the total heat absorbed by the sample, H_overall.

Three DSC experiments, each with two repetitions corresponding to temperature cycles at three different temperature rates, r₁ of 1.50°C/min, 0.55°C/min and 0.50°C/min, with isothermal dwells (r₂ = 0.00°C/min) at 180°C, 175°C and 185°C, respectively, have been conducted, as shown in Figure 1. The temperature for all the DSC tests are gradually reduced to room temperature at a rate of r₃ = 1.50°C/min. The cure reaction rates, dX/dT and glass transition temperatures, T_g corresponding to the three DSC experiments are observed in Figure 2. In the initial stages of the curing process, the reaction rate is relatively slow due to the presence of unreacted functional groups within the resin. As curing progresses, the number of these unreacted functional groups decreases, leading to an increase in the reaction rate. Eventually, the reaction rate slows down again as the resin approaches complete cure, indicated by the T_g. This evolution of dX/dT results in the sigmoid-shaped graph for X, as demonstrated in Figure 1.

Figure 1.

Experimental DSC datasets with the cure state X at different thermal loading conditions.

Figure 2.

DSC outputs, the cure reaction rate (dX/dT) and the glass transition temperature (T_g) at different thermal loading conditions.

If the model parameters from studies^11,39 are employed, they do not accurately correlate with the observed dX/dT and T_g at different thermal loading conditions. Moreover, the relative difference between measured and predicted cure state variables are between 15%–20%, resulting in incorrect cure state transitions. Consequently, there is a need for stochastic cure model parameters in numerical analyses to accurately quantify the uncertainties associated with thermal loading conditions, which make the numerical analyses computationally expensive.

Due to the inaccuracies encountered in handling thermal loading uncertainties with deterministic parameters using the diffusion cure-kinetics model for the case of Hexcel epoxy 8552 resin, a simpler non-parametric model that relies on the information derived from DSC characterization tests is needed. In this work, such model is introduced in the form of neural networks.^67–69 As a supervised machine learning model, it enables capturing the non-linear relationships between input and output variables without assuming a specific functional form, unlike the Arrhenius model. This makes it an efficient tool for predicting the cure state variables of polymeric resin based on the data extracted from available DSC characterization tests.

The neural network structure consists of an input and output layers, and 2 hidden layers, each with 13 hidden neurons, determined as a compromise between model accuracy and complexity. The ANN is designed to predict the dX/dT and T_g through its output nodes, processing the X, the cure temperature, T and the thermal loading conditions, r_{i = 1,2,3} as inputs. This allows for taking varying rates and maximum curing temperatures into account.

The training algorithm is executed using the MATLAB Neural Network Toolbox, specifically enabling the cross-validation with early stopping as a means toward network regularization. Specifically, each of the three DSC tests consists of 2000 samples of which 1500, 300 and 200 samples are used for training, validating and testing the neural network model respectively. Despite establishing an optimized network architecture, the possibility of failing to achieve practical solutions persists. This is attributed to the network initialization commencing with distinct initial conditions and samplings. As a preventive measure, a total of 200 neural networks are trained utilizing diverse initial conditions through the Levenberg-Marquardt algorithm.

The neural network model with the lowest RMSE is chosen to predict the evolution of cure-state variables for the test case of a DSC experiment involving a heating ramp, r₁ of 1.00°C/min, as referenced in Van Ee and Poursartip ⁷⁰ and Shahkarami et al. ⁷¹ The comparison between the predictions of dX/dT between the neural network and diffusion cure-kinetics models is shown in Figure 3. It is observed that the peak and post summit regime of evolution of dX/dT from the neural network model correlates closer to experimental measures. While the parametric model captures the non-linear cure evolution, a trained neural network based model is observed to be the better choice due to it’s capability of handling non-linear relationships accurately at different thermal loading conditions.

Figure 3.

Comparison of the prediction of dX/dT with neural network and cure-kinetics models for r₁ = 1°C/min^70,71.

Results and discussion

The numerical analysis entails chemo-thermo-mechanical FE calculations performed by the commercial software ABAQUS standard/explicit, which is integrated with a FORTRAN user material subroutine. The model aligns with the approach of Svanberg and Holmberg, ⁴³ implementing a simplified linear viscoelastic constitutive model that incorporates superposition between cure state variables and curing conditions, as detailed in Section Cure kinetics and material models. The flowchart illustrating the procedure implemented is depicted in the Pseudocode 1.

This simplified constitutive model is implemented under the following assumptions:

• The CTE remain independent of X.

• The cure-kinetics model for the numerical analysis is taken from Ersoy et al. ¹¹

• The thermosetting material is modelled as orthotropic in the FE based constitutive model.

• The thermo-material properties with a fibre volume fraction, V_f =60% are detailed in Tables 1 and 2 as outlined in studies.^47,54 The effective properties that are missing are determined through micro-mechanical models, referring to works.^72–74 These models rely on semi-analytical and/or analytical expressions derived from the properties of the resin and fibres.

• It is assumed that there is sufficient stress build-up raise in the viscous state.

• The fibre resin matrix is challenging to model and is considered to be homogeneous at the ply-scale level.

• The mechanical properties of the ply laminates are maintained as constants across all states, and their continuous evolution in cooling phase is correlated using slopes, as illustrated in Figure 4.

Table 1.

Unidirectional prepreg AS4/8552 properties at different cure states.

Properties	Rubbery state	Glassy state
Longitudinal modulus, E X (MPa)	137665.50	137770.70
Transitional slope for E X in cooling phase (MPa/°C)		3.43
Transversal modulus, E Y (MPa)	165	7070.70
Transitional slope for E Y in cooling phase (MPa/°C)		13.51
Through thickness modulus, E Z (MPa)	165	7070.70
Transitional slope for E Z in cooling phase (MPa/°C)		13.51
In plane shear modulus, G XY (MPa)	44.30	3404.40
Transitional slope for G XY in cooling phase (MPa/°C)		12.50
Out of plane shear modulus, G XZ (MPa)	44.30	3404.40
Transitional slope for G XZ in cooling phase (MPa/°C)		12.50
Out of plane shear modulus, G YZ (MPa)	41.60	2594.10
Transitional slope for G YZ in cooling phase (MPa/°C)		5.77
In plane Poisson’s ratio, μ YX (−)	4.00 × 10−4	1.89 × 10−2
Out of plane Poisson’s ratio, μ ZX (−)	4.00 × 10−4	1.89 × 10−2
Out of plane Poisson’s ratio, μ YZ (−)	0.98	0.45
CTE in longitudinal direction, α X (1/°C)	−8.78 × 10−7	−7.98 × 10−8
CTE in transverse direction, α Y (1/°C)	4.06 × 10−5	3.05 × 10−5
CTE in through thickness direction, α Z (1/°C)	4.06 × 10−5	3.05 × 10−5
CCS in longitudinal direction, β X (−)	1.07 × 10−4	6.97 × 10−3
CCS in transverse direction, β Y (−)	−1.82 × 10−2	−1.80 × 10−2
CCS in through thickness direction, β Z (−)	−1.82 × 10−2	−1.80 × 10−2

Table 2.

Woven prepreg AS4/8552 properties at different cure states.

Properties	Rubbery state	Glassy state
Longitudinal modulus, E X (MPa)	66190	68000
Transitional slope for E X in cooling phase (MPa/°C)		8.87
Transversal modulus, E Y (MPa)	66190	68000
Transitional slope for E Y in cooling phase (MPa/°C)		8.87
Through thickness modulus, E Z (MPa)	165	10,000
Transitional slope for E Z in cooling phase (MPa/°C)		13.47
In plane shear modulus, G XY (MPa)	44.30	5000
Transitional slope for G XY in cooling phase (MPa/°C)		12.02
Out of plane shear modulus, G XZ (MPa)	42.90	4500
Transitional slope for G XZ in cooling phase (MPa/°C)		8.07
Out of plane shear modulus, G YZ (MPa)	42.90	4500
Transitional slope for G YZ in cooling phase (MPa/°C)		8.07
In plane Poisson’s ratio, μ YX (−)	1.00 × 10−3	0.22
Out of plane Poisson’s ratio, μ ZX (−)	2.00 × 10−3	0.07
Out of plane Poisson’s ratio, μ YZ (−)	0.83	0.49
CTE in longitudinal direction, α X (1/°C)	−7.11 × 10−7	2.50 × 10−6
CTE in transverse direction, α Y (1/°C)	−7.11 × 10−7	2.50 × 10−6
CTE in through thickness direction, α Z (1/°C)	3.50 × 10−4	5.88 × 10−5
CCS in longitudinal direction, β X (−)	9.31 × 10−4	5.29 × 10−2
CCS in transverse direction, β Y (−)	9.31 × 10−4	5.29 × 10−2
CCS in through thickness direction, β Z (−)	−3.45 × 10−2	−2.53 × 10−2

Figure 4.

Adaption of mechanical properties in cooling phase (e.g. in-plane modulus with transitional slope of 8.873 MPa/°C) for numerical modelling of woven 8552/AS4 prepegs.

Numerical analysis and discussion on flat plate thermoset part

A flat plate measuring 150 mm × 150 mm x 2.16 mm is considered, consisting of 8 plies of unidirectional prepreg AS4/8552 arranged as [(0)₄/(90)₄]. This configuration represents an exaggerated distorted case of an asymmetrical stacking sequence. The FE model is depicted in Figure 5. The composite part being modelled is manufactured on an aluminium mould with a CTE, α_mould of 24 × 10⁻⁶ (1/°C). Contact between the part and the mould surfaces is established with a static friction coefficient of 0.2 and a shear stress limit of 0.5 MPa. ²¹ The relevant thermal properties for solving Fourier’s heat conduction, such as ζ, ζ_m, c_p, and L, are given as 1.59 × 10⁻⁶ kg/mm³, 1.30 × 10⁻⁶ kg/mm³, 1289 J/Kg°C, and 1.17 × 10⁻³ W/mm°C, respectively.^39,61,76

Figure 5.

Illustration of the FE model for a flat plate and its associated boundary conditions.

The 3D model of the part and mould includes 28,800 and 3972 hexahedral C3D8 elements, both with linear geometric order, respectively. The cure cycles involve an initial temperature increase at rate of r₁ = 1.00°C/min, reaching up to 180°C. This is followed by an isothermal hold (r₂ = 0.00°C/min) for 120 minutes. Finally, the temperature is gradually reduced to room temperature at a rate of r₃ = 1.50°C/min.

Unidirectional plies exhibit significant variations in CTEs, CCSs and stiffness between different directions. As a result, the temperature evolution within the cure cycle leads to varying levels of volumetric free strains and residual stresses within laminates containing fibres oriented in different directions. The distortion shape after de-moulding, considering the adopted asymmetrical configuration, is shown in Figure 6. This shape aligns well with the findings presented in the work by Svanberg. ² The substantial differences in thermo-chemo and thermo-mechanical properties between the 0°and 90°plies contribute to the pronounced distortion, as observed in Cho et al. ⁷⁵

Figure 6.

Distorted saddle shape prediction comparison between cure-kinetics-based and neural network-based constitutive models upon de-moulding (deformation scaling factor: 2.0).

To assess the sensitivity of both models, we employed cure cycles by adjusting the heating ramp, r₁, with both one positive standard deviation and one negative standard deviation at 25%. The resulting values of maximum displacements serve as an indicator of the model’s sensitivity to the variability of r₁.

The comparison of displacements in the thickness direction with FE numerical analysis based on variale thermal loading, r₁ are shown in Table 3. The observations made from Table 3 reveal a significant average relative percentile difference of 37.63% in the final displacements along the thickness direction. As observed earlier in Figure 3, the neural network model predicts accurately the molecular mobility of the resin reaction. Consequentially, this associates the plies with higher volumetric free strains, ɛ^vf, leading to a relatively higher level of internal residual field.

Table 3.

Maximum displacement comparison at different thermal loading conditions, r₁.

Thermal loading condition (r1)	Neural network interfaced model	Cure-kinetics interfaced model
1.25°C/min	8.324	5.667
1.00°C/min	8.283	5.652
0.75°C/min	8.199	5.644

Meanwhile, when employing the diffusion cure-kinetics model, the dX/dT is notably lower than anticipated, leading to a lesser internal residual field, which, upon de-moulding, results in less displacement. Furthermore, the standard deviation in the maximum displacement predicted by the neural network-interfaced constitutive model among the three cases of r₁ is 5.45 times greater than that predicted by the diffusion cure-kinetics-interfaced constitutive model. This is because the ɛ^vf is highly influenced by the resin molecular mobility, with an increase or decrease in r₁. This demonstrates the significant influence of the thermal loading conditions on final displacements.

The example of the flat plate analysis is extended with 32 plies of prepreg arranged as [(0)₁₆/(90)₁₆] with a total thickness of 8.64 mm. The numerical analysis, based on neural network and diffusion-interfaced cure kinetics models, has been conducted using modified cure cycle, with central plies experiencing temperature overshoots. This is depicted in Figure 7. The maximum temperature overshoot occurs between in the central plies due to the heat released during the cure reaction, leading to an increase in temperature compared to the boundary temperature. Meanwhile, the regions near the top and bottom surfaces lose some of the heat through heat exchange, following more or less the control temperature. When dealing with thicker laminates, it becomes crucial to comprehend the spatial variation in ɛ^vf resulting from variations in X.

Figure 7.

Temperature distribution in central plies across thick laminates estimated with Fourier heat conduction analysis.

The observations from Figure 8 reveal a substantial percentile difference, ranging from 43% to 52%, in localized residual stress along the thickness direction. This disparity can be attributed to the thermal lag, which results in a delay in reaching the X_gel, coupled with a temperature overshoot that accelerates the attainment of the X_vitr. Consequently, the central plies associated with higher ɛ^vf cure faster when compared to the outer plies, with modified thermal loading condition. When employing the diffusion cure-kinetics model, the dX/dT and and T_g, assumed to be independent of the thermal loading condition, leads to a relatively slower cure of central plies. As a result, this underestimates the constraint imposed by fully cured central plies on the outer plies. With examples of thicker parts, it is observed that potential errors associated with thermal gradients in simple flat plates are effectively circumvented by employing the non-parametric cure model.

Figure 8.

Comparison of residual stress field upon de-moulding in thickness direction, S _Z (deformation scaling factor: 2.0).

Case study: Numerical analysis and discussion on Z-shaped thermoset part

The methods described in Section Cure kinetics and material models and implemented in Section Characterization tests and neural network model on an example of a flat plate are now implemented on a case study consisting of a Z-shaped part, as shown in Figure 9. The part is composed of 15 plies with a stacking sequence, [±45°/0 − 90°]_4s, and a length of 150 mm along the Z-axis. Materials are the woven prepregs AS4/8552, produced by Hexcel ⁷⁶ ; is under consideration for manufacturing. The part is created using a mould with a negligible CTE, α_mould = 0. The numerical analysis involves a cure cycle featuring a heating ramp with r₁ = 1.50°C/min, gradually raising the temperature from room temperature to 180°C. Subsequently, the component is held at an isothermal dwell of 180°C (r₂ = 0.00°C/min) for a duration of 120 minutes before gradually cooling back to room temperature at a rate, r₃ = 1.50°C/min. Throughout the entire cure cycle, a pressure of 7 bars is consistently applied to the part.

Figure 9.

Z-shaped part and its associated geometry.

The numerical analysis for the case of Z-shaped part involves additional assumptions:

• An Invar mould with a negligible CTE is considered, thereby being modelled as a rigid body within the FE model.

• A frictionless contact condition is assumed to exist between the part and the tool.

• The applied temperature cycle is assumed to be uniform across the entire part, thus rendering the temperature gradient across the thickness negligible.

The numerical analysis considers an autoclave thermal loading condition, with a pressure cycle applied to the exposed surface of the Z-shaped part. To create the 3D geometry, extrusion and meshing are performed in the direction of the defined shell geometry. Before extrusion, the materials, thicknesses and orientations of the plies are established, and draping modelling is executed with the assistance of the Simulayt Composite Modeller within ABAQUS. The FE numerical model is illustrated in Figure 10. The 3D model of the part includes 19,665 hexahedral C3D8 elements with a linear geometric order, along with assigned local orientations. The mould is considered as a rigid component, consisting of 2296 quadrilateral R3D4 elements with linear geometric order. The purpose of the imposed boundary condition in the numerical model is to closely resemble the experimental reality. In the actual scenario, the pressure holds the part against the mould, effectively preventing the curved sections in the Z-shaped part from closing freely. Finally, during de-moulding, the pressure and contact with the mould are suppressed, resulting in the component distortions.

Figure 10.

Illustration of the FE model of a Z-shape part subject to autoclave thermal loading conditions.

The evolution of volumetric free strains on an element of the front flange in the Z-shaped part during the curing process is depicted in Figure 11. The phenomenon of resin shrinkage dominates before the onset of vitrification. Upon vitrification, thermal effects come into play during the cooling phase of the cure cycle. The neural network model’s T_g reaches the process temperature and transitions to the glassy state much earlier compared to the diffusion cure kinetics model. The relative difference in strain between the X and Y directions with the neural network based simulation (as compared with the standard cure-kinetics model) upon de-moulding is 15.41% due to the earlier onset of the vitrification phenomenon. This is primarily because of the higher dX/dT predicted by the neural network model, enabling rapid resin molecular mobility. The volumetric free strains are integrated within the constitutive model to predict the the internal residual stress, S. According to equation (18), during the rubbery state, the constitutive model assumes the material to be fully relaxed. Conversely, in the glassy state, the S are stored incrementally, and it is finally released upon de-moulding.

Figure 11.

Comparison of volumetric free strain evolution between cure-kinetics-based and neural network-based constitutive models.

Figures 12–14 show a comparison of internal residual stresses, S, in the X, Y, and Z directions obtained by simulation with standard diffusion cure-kinetics and neural network interfaced constitutive models. The stress states before the cooling phase, after the cooling phase, and upon de-moulding are analysed. The differences before the cooling phase range between 42% and 44% in the primary directions. Following the earlier onset of vitrification, volumetric free strains are coupled with a significant increase in mechanical properties, resulting in the development of tensile residual stresses on the flanges and web. As a result, a relatively higher amount of residual stresses is observed in the proposed model before the cooling phase. By the end of the isothermal hold, the material has fully vitrified into the glassy state, and the evolution of the cure has ceased, stabilizing the material properties.

Figure 12.

Comparison of residual stress prediction in X-direction, S _X (deformation scaling factor: 5.0).

Figure 13.

Comparison of residual stress prediction in Y-direction, S _Y (deformation scaling factor: 5.0).

Figure 14.

Comparison of residual stress prediction in Z-direction, S _Z (deformation scaling factor: 5.0).

On the contrary, the standard diffusion cure-kinetics model is associated with relatively low volumetric free strains and a delayed onset of vitrification. This results in reduced tensile residual stresses. During the cooling phase, although cure shrinkage is negligible, there is a strong influence of the thermo-elastic effect on homogeneous mechanical properties. This effect leads to an increase in tensile residual stresses. Upon removing the boundary conditions in the de-moulding phase, the majority of the induced tensile residual stress field is released, resulting in distortions.

Validation of the predicted internal residual stress field is accomplished through a comparison with distortions. In the case study, the distortions in the form of spring-in angles were measured using the point clouds acquired with the Nikon MMDx100 laser scanner, followed by post-processing as described by Parmentier et al. ⁵⁴ These points are employed to establish a map of deviations between the initial reference and distorted geometry. The experimental measures of spring-in angles exhibit a standard deviation of 0.2°.

Comparison of average value of spring-in angles, A₁ and A₂ post de-moulding between the viscoelastic models interfaced with the standard diffusion cure-kinetics model and the neural network model is depicted in Figure 15. An example of post-processed spring-in predictions with the neural network-based constitutive model along various points in the curved sections is shown to validate the measured spring-in angles. Regarding the average spring-in angle predictions, the proposed numerical analysis provides a prediction with a relative error of 2.31% and 10.20% for angles A₁ and A₂, respectively. Meanwhile, the diffusion cure-kinetics interfaced model provides a prediction with a relative error of 29.90% and 35.20% for angles A₁ and A₂, respectively. While the proposed approach yields accurate predictions, a disparity in angle A₂ arises. This occurs when dealing with curved parts, as the pressure from the vacuum bag pushes the resin out of the radius of curvature region, resulting in a fibre-rich area. Conversely, when the vacuum bag does not perfectly follow the radius of curvature region, resin could accumulate between the bag and prepreg, leading to a resin-rich region. Hence, relatively more internal residual stress are retrieved during the curing process in the A₂ region, resulting in a relatively larger spring-in angle.

Figure 15.

(a) Spring-in angles measures at different sections using neural network based constitutive model and (b) Comparison of average spring-in measures between different models.

This issue is addressed by modelling a thin layup of hexahedral C3D8 elements, with a thickness of approximately 0.02 mm, corresponding to the isotropic properties of the 8552 resin in the resin-rich region. The comparison of the average values of spring-in angles, A₁ and A₂, post de-moulding with the resin-rich region is illustrated in Figure 16. The proposed numerical analysis, with the mitigation technique, provides a closer prediction with a relative error of 5.10% for angle A₂. Upon mitigation, the region around angle A₂, as observed in Figure 16, experiences more chemical shrinkage compared to the region near angle A₁, leading to an increase in localized internal residual stress. The localized residual stress, S _X is approximately 10 times when compared with the former case without modelled resin-rich region in Figure 12.

Figure 16.

(a) Comparison of average spring-in measures between different models, with resin and fibre-rich regions modelled and (b) Residual stress distribution in X-direction, S _X upon de-moulding with resin and fibre-rich regions modelled using neural network based constitutive model.

A comparison of average computational costs and memory requirements with the Abaqus implicit solver between both models, using 6 CPU cores per job in a High-Performance Computing environment during numerical analyses, is presented in Table 4. There is a slight increase in the costs associated with the numerical run using the proposed approach, primarily due to the non-linearity of the activation functions used in the trained neural network model. Despite this slight increase in computational aspect and training time, the model can be effectively generalized to accurately predict the localization of internal stress fields and process-induced defects.

Table 4.

Comparison of average computation cost and memory requirements with Abaqus implicit solver.

Model	Total memory (Megabytes)	Total CPU time (Hours)
Cure-kinetics interfaced constitutive model	992.272	01: 53: 04
Neural network interfaced constitutive model	1087.676	02: 03: 18

Summary and conclusions

Process-induced defects and residual stresses are inevitable during the manufacturing process due to the mismatch of free volumetric free strains resulting from thermal expansion and chemical shrinkage. The mechanism behind development of defects occurs at different mechanical levels. When assessing the residual stress field, the primary interest lies in the evaluation of distortions and performance-related attributes such as micro-cracks and reduced mechanical properties. Such defects could potentially be mitigated by adopting an optimal temperature cycle, which necessitates numerous iterative numerical analyses before the manufacturing process.

A linear viscoelastic constitutive model, based on the path dependence of cure state variables, has been implemented for numerical analysis on thermosetting composite parts subjected to autoclave conditions. Moreover, the numerical modelling of residual stress relies on the homogenized mechanical and thermo-chemo properties of the resin and fibre, depending on the curing behaviour. Therefore, accurately capturing the transitions between the viscous, rubbery, and glassy states is crucial. These cure behaviour and state transitions, which depend on cure state variables, are defined using cure-kinetics and DiBenedetto’s models, involving numerous constants. The modelling of such a cure-kinetics model requires extensive knowledge of various Arrhenius models and prior information about model parameters. However, such parametric cure models have been observed to be unable to accurately capture the non-linear relationships between the cure state variables and the thermal loading condition variables. These complex relationships are captured in this study by using a non-parametric neural network model that, however, requires extensive training. The study provides insights into the influence of thermal loading rates (below 3°C/min) on the evolution of cure state parameters for carbon/epoxy prepreg AS4/8552. The neural network model is integrated with the constitutive model to enable accurate predictions of process-induced distortions in the case study of a Z-shaped thermoset part. Several significant conclusions and remarks from the study are as follows:

• It has been observed that the cure state variables are dependent on process conditions, specifically the thermal loading rates.

• To capture the non-linear relationship between cure kinetics and process condition variables, a non-parametric neural network model is implemented using DSC characterization tests conducted under various thermal loading conditions. Such a model has been found to be particularly useful when dealing with thick thermoset parts that are subject to temperature gradients.

• When tested with new data for the carbon/epoxy prepreg AS4/8552, the neural network model provides more accurate state transitions between the viscous, rubbery, and glassy states compared to the conventional diffusion cure-kinetics and DiBenedetto’s models under below 3°C/min.

• For the case study of the Z-shaped thermosetting part, the proposed approach predicts process-induced distortions in the form of spring-in angles that correlate more closely with experimental measures obtained from laser scans.

By training and implementing the neural network model, the expenses tied to comprehensive characterizations under diverse boundary conditions for thermosetting composite materials can be circumvented, allowing the attainment of desired intermittent temperature rates and maximum curing temperatures. The trained networks can subsequently be integrated with the constitutive models to forecast process-induced defects. This approach furnishes an initial estimation of defects and streamlines the optimization of temperature profiles to mitigate risks and enhance manufacturing quality.