Stochastic simulation of the inﬂuence of cure kinetics uncertainty on composites cure

A stochastic cure simulation methodology is developed and implemented to investigate the inﬂuence of cure kinetics uncertainty due to different initial resin state on the process of cure. The simulation addresses heat transfer effects and allows quantiﬁcation of uncertainty in temperature overshoot during the cure. Differential Scanning Calorimetry was used to characterise cure kinetics variability of a commercial epoxy resin used in aerospace applications. It was found that cure kinetics uncertainty is associated with variations in the initial degree of cure, activation energy and reaction order. A cure simulation model was coupled with conventional Monte Carlo and an implementation of the Probabilistic Collocation Method. Both simulation schemes are capable of capturing variability propagation, with the collocation method presenting beneﬁts in terms of computational cost against the Monte Carlo scheme with comparable accuracy. Simulation of the cure of a carbon ﬁbre–epoxy panel shows that cure kinetics uncertainty can cause considerable variability in the process outcome with a coefﬁcient of variation in temperature overshoot of about 30%.


Introduction
Uncertainty in composite manufacture has begun to receive attention in the last decade as a result of the adoption of advanced composites as the main material of choice for large aero-structures.Variability is present in all forms of pre-impregnated and dry textiles and is mainly associated with in-plane and out-of plane fibre misalignment generated during production, handling or storage [1][2][3].Experimental results have shown that in-plane fibre misalignment can be represented by a normally distributed variable [2,4] exhibiting strong spatial autocorrelation over the textile [2].This type of uncertainty can affect the forming process introducing considerable variability in defect generation [2].In addition, fibre misalignment as well as nesting effects can introduce significant scatter in permeability contributing to the formation of voids or dry patches [5][6][7][8].Similarly to fibre misalignment, macroscopic permeability values can be represented by a normally distributed random variable [5][6][7][8].However, simulation studies have indicated that permeability at the meso-scale (unit cell size) shows higher variations due to local inhomogeneities and cannot be represented by a normally distributed variable [9], a result that is also valid in the case of random mats [10].In addition to variations in fibre architecture, the cure process can be potentially influenced by environmental/boundary condition uncertainty as well as resin behaviour variability due to variation in handling and storage conditions.Naturally these effects can play a role in cure process defects such as severe temperature overshoots or under-cure and also introduce variability in residual stresses/shape distortion.
The issue of variability during the cure step has received limited attention so far in the literature.The effect of cure temperature variations and cure kinetics uncertainty on cure time has been investigated in a pure simulation study by coupling a cure kinetics model with a Latin Hypercube sampling scheme showing that cure temperature variations tend to dominate cure time variability [11].Furthermore, taking into account uncertainty in the optimisation of the cure time, has shown that optimal cure time increases with increasing variability [12].These results, which are based on simulation and hypothesised levels of uncertainty, gain significant practical importance when their conclusions are combined with experimental studies of uncertainty in cure kinetics.In addition to material behaviour and process parameters uncertainty, cure kinetics can show significant variations due to experimental characterisation and data reduction discrepancies [13].
This study aims at the quantification of cure kinetics uncertainty and investigation of its propagation through the manufacturing process.The variability can be due to variable resin handling/ storage conditions, variations in resin formulation and uncertainty introduced during resin characterisation.A series of experiments was carried out using Differential Scanning Calorimetry (DSC) to characterise cure kinetics uncertainty of a commercial epoxy resin used in aerospace applications.The variability in experimental behaviour was attributed to certain parameters of cure kinetics as expressed by a phenomenological model and the corresponding stochastic object was developed.
The resulting stochastic simulation problem was addressed by coupling a finite element analysis based cure simulation model with conventional Monte Carlo (MC) and the Probabilistic Collocation Method (PCM).The two stochastic simulation approaches were applied to the cure process of a thick carbon fibre-epoxy laminate to investigate the influence of cure kinetics uncertainty on variability in temperature overshoot.Temperature overshoot is defined as the maximum difference between the laminate temperature and mould temperature during the cure cycle.The two schemes were compared in terms of accuracy and efficiency.

Deterministic cure simulation model
A simulation of heat transfer effects occurring during the cure was implemented in the finite element analysis solver MSC.Marc.The model was three dimensional and transient.The materials considered were Hexcel G1157 [14] pseudo unidirectional carbon fibre reinforcement and Hexcel RTM6 [15] epoxy resin.The material properties depend on both temperature and degree of cure and the material sub-models of cure kinetics, specific heat capacity and thermal conductivity were implemented in user defined subroutines UCURE, USPCHT and ANKOND [16].

Cure kinetics model
The cure kinetics model used has been developed for the resin system of this study and is a combination of an nth order model and an autocatalytic model [17,18].
The cure reaction rate is expressed as follows: where a is the current degree of cure, k 1 , k 2 the reaction rate constants, and m, n 1 , n 2 the reaction orders.The reaction rate constants incorporated diffusion rate limitation terms [19]: where k i;c are chemical rate constants following an Arrhenius temperature dependence, and k d is a diffusion rate constant defined as [19]: here T is the instantaneous temperature, E i activation energies, A i pre-exponential factors, A D and E D the pre-exponential factor and activation energy of the diffusion process respectively, R the universal gas constant, b a constant, and f the equilibrium free volume expressed as [19]: where w and g are constants, whilst T g is the instantaneous glass transition temperature defined as [20]: where T go and T g1 denote the glass transition temperatures for the uncured and fully cured material, respectively and k is a fitting parameter [18].The parameters of the glass transition development model are reported in Table 1.

Specific heat capacity
The composite specific heat capacity is calculated using the rule of mixtures as follows: where w f is the weight fraction of the fibre, c pf the fibre specific heat capacity and c pr the specific heat capacity of the resin.The fibre specific heat capacity is a linear function of temperature and can be expressed as follows [21]: where A fc p is the slope and B f cp the intercept of the linear function [21].The values of these two coefficients for carbon fibre [21] are reported in Table 1.
The specific heat capacity of the resin depends on both temperature and degree of cure.The dependence on degree of cure is expressed using a dependence on the instantaneous glass transition temperature of the resin (Eq.( 6)) as follows: here A rcp and B rcp are constants expressing the linear dependence of the resin specific heat capacity on temperature for constant material state, while D rcp ; C rcp and s are constants referring to the strength, breadth and temperature shift of the transition around T g .The values of the coefficients used in Eq. ( 9) were estimated by fitting to experimental data produced by modulated scanning calorimetry [21] using a Genetic Algorithm and are reported in Table 1.

Thermal conductivity tensor
Each composite lamina is considered a transversely isotropic material.
The thermal conductivity in the longitudinal direction can be calculated as follows [22]: K lf and K r are the thermal conductivity of the fibre in the longitudinal direction and of the resin respectively.In the transverse direction the thermal conductivity can be computed as follows [22]: Table 1 Parameter values for glass transition temperature, specific heat and thermal conductivity sub-models [18,20,21,23].

Parameter
Value where K tf is the fibre conductivity in the transverse direction.The fibre conductivity in the longitudinal direction is as follows [23]: while the fibre conductivity in the transverse direction is [23]: The thermal conductivity of the resin can be expressed as [21]: The parameters used in Eqs. ( 10)-( 14) are summarised in Table 1.

Analysis of cure kinetics uncertainty
Investigation of cure kinetics uncertainty is based on a series of DSC tests.Here, the aim is to quantify cure kinetics uncertainty by fitting the experimental data with the cure kinetics model described in Section 2.1.

Experimental results
DSC tests were carried out at a heating rate of 1 °C/min from 80 °C to 240 °C after equilibration at the initial temperature.Samples from four different batches were tested.All batches were within their shelf life at the time of testing.Tests were duplicated within each batch.Integration of the heat flow versus time data using an iterative baseline [24] allows calculation of the degree of cure evolution and the corresponding reaction rate.Fig. 1 illustrates the results expressed as reaction rate versus temperature for the eight tests carried out.All curves have the same qualitative characteristics with a peak at intermediate temperatures and a shoulder towards the end of the reaction.The main peak of the reaction tends to be slightly asymmetric and ends with a plateau at a low reaction rate.This is followed by a drop to negligible reaction at very high temperatures.The repeatability within the same batch is very high with curves being almost identical.In particular, the absolute differences between curves of the same batch for values of degree of cure of 20%, 40%, 60% and 80% varied from 8EÀ06 to 2EÀ05 [1/s], which is more than one order lower than the values of cure reaction rate obtained by the experiments.This shows that the experimental and signal analysis method related to sample preparation, measurement quality and baseline decision induces negligible variations to the results.In contrast, significant variability can be observed in the comparison between batches, implying that variation in cure behaviour is attributed to batch to batch variability caused by variations in thermal history related to different resin handling/storage conditions.This variability is manifested in the maximum reaction rate, the position of the maximum, the temperature at which the reaction starts and the temperature of the high temperature shoulder.

Quantification of cure kinetics uncertainty
The procedure for the estimation of the cure kinetics model parameters for the different tests is shown schematically in Fig. 2. The parameters in Eqs.1-5 were estimated using the hybrid Genetic Algorithm implemented in the Solver of Microsoft Excel [25].In addition to these parameters the procedure estimates the initial degree of cure which is involved in the kinetics model as the initial condition of the integration.The overall model for the kinetics parameters is obtained by carrying out an overall fitting based on published values for the resin system [17][18][19] in order to define relatively narrow ranges for the search.The mean value of the total heats of reaction computed after integration of each experimental curve was used to obtain an initial guess for the initial degree of cure.The determination of the cure model parameters and initial degree of cure for each sample then follows by fitting each experimental curve to the cure kinetics model.This procedure estimates the mean value l and standard deviation r of each parameter.This procedure can potentially be influenced by the non-linearity of the cure kinetics problem and the interdependencies of some of the fitting parameters making subsequent analysis steps phenomenological.Table 2 reports the values of the kinetics parameters estimated following this procedure.A sensitivity analysis was carried out by using the cure kinetics model described in Section 2.1 to determine which of the parameters should be considered as stochastic.Each parameter was varied by one positive standard deviation and one negative and the average relative absolute differences of predicted reaction rates from the cure kinetics model using the estimated parameter value was computed.The mean between these two values was used as an indication of the model sensitivity to the level of variability of each of the model parameters (Table 2).It should be noted that this relative difference corresponds to increments equal to one standard deviation around the mean, i.e. a 68% probability assuming a normal distribution of the variables.
The activation energy E 2 , the reaction order m and the initial degree of cure a o introduced the highest discrepancies, presenting a relative difference of 11%, 10% and 4%, respectively.Consequently, the main sources of uncertainty in the experimental results are the variability in the initial degree of cure a o , activation energy E 2 and reaction order m.
The initial degree of cure can influence the temperature of reaction onset introducing a shift to the cure reaction rate -temperature curve.The reaction order m and the activation energy E 2 mainly affect the height and the position of the main reaction rate peak.The nth order term parameters in Eq. ( 1), especially E 1 , influence the height of the main reaction rate peak as well; however, given the observed spread of values, they introduce relatively lower discrepancies compared to E 2 , m and a o (Table 2), and therefore can be considered as deterministic.

Statistical properties of cure kinetics
The basic statistical properties of the three stochastic variables are summarised in Table 2.The initial degree of cure shows the highest level of variation among the three stochastic variables.This can be attributed to thermal history variations between the different batches during their storage and transport.Due to the limited number of experimental data, the Kolmogorov-Smirnov goodness of fit test was carried out to investigate which statistical distribution is suitable to fit the stochastic variables.It was suggested that all variables can be represented by a normally distributed random variable with 99% level of confidence.The activation energy E 2 and the reaction order m are strongly correlated, while a o and m show moderate correlation (Table 3).

Stochastic simulation approach
Two stochastic simulation approaches were developed based on conventional Monte Carlo (MC) and the Probabilistic Collocation Method (PCM), respectively.The two stochastic simulation models were coupled with the finite element based cure simulation model.The Monte Carlo method presents a robust but usually computationally expensive solution.The main concept of PCM is to construct a response surface for every output parameter, as a function of uncertain parameters in the form of orthogonal polynomials, termed as the polynomial chaos [26].An uncertainty analysis is then carried out to quantify output parameters uncertainty using this surrogate representation of the response surface.Following this strategy, the stochastic problem is reduced to calculating the unknown polynomial chaos coefficients, which are computed using the probabilistic collocation approach [27,28].The collocation points are the roots of the next higher order orthogonal polynomial than the order of the response surface and are chosen so that the residuals between each response surface and the corresponding model output are zero, implying that the collocation points are selected from regions of high probability.Consequently, a system of linear equations is obtained to calculate the polynomial chaos coefficients for every output parameter [27,28].

Representation of input parameters
Cholesky decomposition was implemented in both stochastic simulation schemes to generate realisations of the stochastic variables.The Cholesky method decomposes the covariance matrix, C as follows [29]: here L is a lower triangular matrix and is the Cholesky root.The product of the Cholesky root with the vector Y of independent identically distributed standard normal variables is a vector V with the desired statistical properties of the stochastic variables, and is defined as follows [29]: Evaluation of Eq. ( 15) has relatively high computational cost; however, it needs to be executed only once.Realisations of the random vector Y are generated and transformed to realisations of  the vector V using Eq. ( 16).This step is computationally inexpensive and is executed as many times as the number of required realisations of the stochastic variables.

Stochastic cure simulation
The cure of a 30 mm thick carbon fibre-epoxy laminate fabricated by infusion was modelled by coupling the two stochastic simulation schemes with the finite element based cure simulation model.The applied cure profile is illustrated in Fig. 3.The lay-up sequence of the laminate was [0°/90°/90°/0°] 25 .The initial temperature was set at 15 °C and was applied to all the nodes of the model.A prescribed temperature boundary condition defined by the cure profile was applied to the nodes in contact with the mould, whereas natural air convection with a surface heat transfer coefficient of 5 W/(m 2 K) was applied on the surface in contact with the vacuum bag.In the collocation method implementation a second order response surface was constructed to represent the output parameters and a modified regression-based collocation approach was implemented to improve accuracy.In this case, the number of collocation points used is larger than the number of the unknown polynomial chaos coefficients, implying that the effect of each collocation point is reduced.In particular, the number of the unknown coefficients for a three dimensional second order polynomial chaos is 10 [27] and 21 collocation points were used in this study, leading to 21 deterministic model runs.This was followed by Monte Carlo simulation using the surrogate model.

Stochastic cure simulation results
As shown in Fig. 3, the temperature at the centre of the laminate is initially lower than the prescribed temperature, whereas a temperature overshoot occurs at the beginning of the second dwell due to high exothermic heat.This is followed by a decrease until the end of the cycle.Similarly, the degree of cure at the centre of the laminate is lower than the nominal value at the initial stage of the process, whilst it increases abruptly at the time of temperature overshoot.This is attributed to the large out of plane temperature gradients, in this case 1.5 °C/mm, due to the low thermal conductivity of the material in the through the thickness direction.
Satisfactory convergence behaviour was obtained for the first and second statistical moments of temperature overshoot after 2000 Monte Carlo iterations for both stochastic simulation schemes (Fig. 4).A very good agreement is achieved between MC and the collocation method.The results suggest that temperature overshoot presents a coefficient of variation of 31%.Examination of the probability distribution of temperature overshoot shown in Fig. 4 indicates that temperature overshoot can be considered a normally distributed random variable.Fig. 5 illustrates the results regarding for the time of temperature overshoot.Similarly to temperature overshoot the simulation converges after 2000 Monte Carlo iterations, with the two stochastic models presenting satisfactory agreement.The time at which temperature overshoot occurs exhibits a coefficient of variation of 1.7%, and can be represented by a normally distributed variable, as illustrated in Fig. 5.The stochastic cure simulation results are summarised in Table 4.
Both MC and PCM have the capability of capturing variability propagation, with the MC presenting a computationally expensive and rich solution and the PCM offering an efficient approximation (e.g. for the given case, the computational cost of the PCM is 1.5% of that of the MC), with comparable accuracy.The results presented here show that cure kinetics uncertainty can introduce significant variability in temperature overshoot with considerable cost implications in industrial practice, especially in the case of ultra-thick complex geometry components where exothermic effects can be severe.Consequently, variability effects should be incorporated in cure profile optimisation to minimise the probability of resin degradation due to temperature overshoot variations while using a process design as efficient as possible in terms of duration and energy consumption.

Concluding remarks
The methodologies developed in this work allow the quantification of the influence of variability on the cure process outcome.The experimental results show that the cure behaviour of high specification thermosets involves uncertainty related to variations in resin thermal history due to different resin transport, handling and storage conditions, variations in resin formulation and characterisation uncertainty, which in turn can introduce significant variability to the process outcome.It is found that majority of uncertainty can be represented by variations in the initial degree of cure, activation energy and reaction order.The stochastic simulation results suggest that temperature overshoot can present a coefficient of variation of about 30%, with potential implications in the amount of scrap during the manufacturing process of composite materials.Furthermore, the time of temperature overshoot has a coefficient of variation of 1.7%.
The modelling approaches presented in this work can be extended to investigate the effect of other sources of uncertainty such as heat convection coefficient and thermal boundary conditions variations as well as fibre misalignment on both heat transfer and process stress effects taking place during the cure.The outcome of this work will allow incorporation of variability in process design/ optimisation to address robustness -performance trade-offs.In addition, development and incorporation of stochastic simulation schemes in existing commercial process simulation tools will lead to efficient and robust process designs.

Fig. 1 .Fig. 2 .
Fig. 1.Evolution of reaction rate as a function of temperature during dynamic cure at 1 °C/min.Letters denote the different batches of resin and numbers different samples within the same batch.

Fig. 3 .
Fig.3.Prescribed temperature boundary condition, evolution of laminate temperature at centre of laminate, nominal degree of cure, degree of cure at centre of laminate.Deterministic model results.

Fig. 4 .Fig. 5 .
Fig. 4. Convergence of statistics of mean (a) and standard deviation (b) of temperature overshoot (c) probability distribution of temperature overshoot.

Table 2
Estimated kinetics parameters, standard deviation of kinetics parameters, coefficient of variation of kinetics parameters, sensitivity analysis results.

Table 3
Correlation matrix of uncertain parameters.

Table 4
Stochastic cure simulation results; temperature overshoot, time of temperature overshoot.