Thermal runaway criterion for thick polymer composites

An analytical solution has been developed for the curing of thick polymer composite laminates that shows how the temperature profile responds to arbitrary changes to the material properties and process parameters and that curing with slow reactions and a low exotherm temperature is impossible if the Damk ¨ ohler number is above a well-defined limit. The thermal runaway criterion can be recast as a criterion for the maximum allowable thickness of the laminate. The thermal runaway criterion was found to agree well with some results for thick laminates from the literature, but the peak temperature in the laminate was underpredicted for stable conditions. The model has a constant that can be adjusted to improve the peak temperature prediction, but more validation data is needed before the model can be optimized to simultaneously predict the peak temperature and thermal runaway with high accuracy.


Introduction
Thick thermosetting polymer composite laminates are used in many technical applications, e.g.pressure vessels, thick-walled pipes, high voltage electrotechnical equipment and marine applications.However, thick components are prone to overheating during cure and this can degrade the laminate properties [1,2].In practice there is an upper limit for the thickness of the laminate when the cure is done with a processing temperature according to the resin manufacturers recommendation [3,4].
In order to understand the behavior of thick parts during cure, detailed computer models can provide accurate information about the temperature in the laminate and the solidification of the matrix material [1,4,5].However, many simulations with different settings are needed to get a complete map of how the solution will change if the material or the processing conditions are changed.An analytical solution, on the other hand, will immediately provide the response to the change in any of the parameters of the problem.
The method of scaling analysis [4] is a helpful tool in theoretical analysis since it reduces the number of independent parameters of the problem.In this method the different terms in the governing equations are scaled with characteristic quantities (temperature T 0 , length L 0 and time t 0 ), which makes it possible to compare the magnitudes of the different terms.In the work by Secord et al. [4] the resulting dimensionless problem was solved numerically for many boundary temperatures and dimensionless problem parameters so that a detailed map over the peak temperature in the laminate could be constructed.From the results it was found that the behavior of thin and thick parts can be well defined and predicted by the following dimensionless groups: a modified Damköhler number defined as  f(α) max , an Arrhenius number defined as E/RT 0 , the dimensionless temperature ramp rate, the dimensionless boundary temperatures, and the exponents in the reaction rate model.The resulting map is quite complex but when the modified Damköhler number is below a certain value the peak temperature in the laminate will always be low.However, the investigation by Secord et al.
was only done for one specific resin-matrix combination and at one fiber volume fraction and many additional simulations would be needed to get a universally valid result.
One hurdle that makes it difficult to derive analytical solutions is the strong non-linear coupling between the reaction rate and the temperature in the laminate.Typically, the results in the literature are therefore derived for idealized conditions [3].However, in the present paper a more general solution is developed by using a perturbation method that is commonly used in fluid dynamics and combustion theory [6].
The present paper begins by developing a perturbation series solution for symmetric heating, followed by a result section where the analytical results are assessed against results from the literature.Finally, the results are discussed in relation to the results from the literature and the main conclusions from the study are given.The details of the analysis, which is quite lengthy, are given in the appendix.E-mail address: rikard.gebart@associated.ltu.se.

Theoretical analysis
The aim of the analysis is to determine the necessary processing conditions to avoid thermal runaway and to achieve a low exothermal peak temperature.Under ideal conditions the heat release from the reactions will be balanced by heat conduction to the mold walls and the resulting temperature increase will be small.However, in thermal runaway the heat release will be several times higher than the heat conduction.This will lead to a feedback loop between increasing temperature leading to increasing heat release leading to increasing temperature and so on until all reactants have been consumed or the heat release is slowed down by decreasing mobility of the reactants.The resulting peak temperature in the laminate will be high but it will always be lower than the adiabatic temperature rise since a significant amount of heat will have time to be conducted to the walls before thermal runaway occurs.In the analysis it is assumed that the laminate is flat, and that the in-plane temperature gradients are negligible, but the solution is valid for curved laminates if the radius of curvature is large compared to the thickness of the laminate thickness [5].
Fig. 1 shows a schematic of the geometry and boundary conditions with symmetric heating, which is often used in resin transfer molding and is approximately correct for autoclave and oven curing.
The analysis follows the same principles as the analysis by Carslaw & Jaeger [7] of the steady temperature in a slab with heat generation [8].The analysis is lengthy and to improve the readability of the paper only an outline of the solution will be presented in the main part of the paper while the details are given in the appendix.
The laminate material with a reactive matrix and a non-reactive and heat conducting fiber phase was approximated as a homogeneous reactive material with equivalent thermal properties that can be estimated from rule-of-mixture [8][9][10] formulas or measured directly using experimental methods.From experiments it is known that the throughthickness thermal conductivity of polymer composites varies with temperature [11].However, the variation is modest, e.g. for a carbon fiber epoxy prepreg laminate the variation is only about 2 % per 10 • C in the range 25 -165 • C. Hence, in the following analysis, where the goal is to find solutions with a limited temperature variation, the thermal conductivity was approximated as constant.
The reactions in the laminate were assumed to be governed by an overall global reaction: The corresponding total heat of reaction q c ( J kg ) and mass reaction rate w for the composite are given by: where V f is the fiber volume fraction, q r is the total heat of reaction for the resin and ρ r is the resin density.The degree of cure α is defined in the customary way [2][3][4][12][13][14] as: where q is the instantaneous heat release from the reactions in the material.Substitution of eqs.( 2) and (3) in the energy equation for the laminate yields: where the coordinate x is defined in Fig. 1 and the thermal properties λ, ρ and c p are the equivalent properties for the composite mentioned above.
In order to close the problem, a mathematical model for dα dt must be defined.There are many proposals in the literature [14], but most of these models have the same general form: The common factor Ae − E RT is an Arrhenius type expression that describes how the reaction rate depends on the absolute temperature (K).For simplicity and to avoid misunderstandings, all temperatures in the rest of the paper are therefore absolute temperatures.The main difference between the models lies in how the dimensionless function f(α) depends on the degree of cure α [15].This function has the same qualitative behavior for all models, with a low value at low degree of cure, a maximum at some intermediate degree of cure and a vanishing value as the degree of cure approaches an upper limit.The upper limit for α is either unity or a slightly lower value, depending on the choice of model.Quantitatively the main difference between the models lies in the behavior at low and high degrees of cure when the reaction rate is low, but for the main part of the curing process when the reaction rate is high, they all agree closely if the model parameters have been properly fitted to the same experimental data.
The complete problem described by eqs.( 5) and ( 6) is very hard to solve analytically since both equations must be solved simultaneously and because of the non-linearity.However, the two equations can be decoupled by approximating the function f(α) as a constant C α : This is a very conservative approximation if the value of C α is set to the maximum of f(α).However, a slightly smaller value should still give a conservative prediction since the reaction rate is over-estimated both at high and low degree of cure.The approximation in Eq. ( 7) can also be motivated by the experimental observation that the sensitivity of the reaction rate to typical changes in the degree of cure is much smaller than the sensitivity to typical changes in temperature.
In order to determine suitable values of C α , the behavior of three different resins from the literature [5,12,15] were compared for a model where the model parameters A, E, n and m were determined by a leastsquare fit of Eq. ( 8) to DSC data.The resulting parameters are shown in Table 1.
The behavior for one of the resins in Table 1 is shown in Fig. 2, where it can be seen that dα dt changes by one order of magnitude when the temperature changes with 30 • K, but with less than a factor of 2 when α changes from 0.03 to 1.
For the model in Eq. ( 8) the constant C α must obey: Fig. 1.Schematic of the geometry and the coordinate system used for the analysis.
where the upper limit has been derived using calculus.The optimum choice of C α , in the sense that it will minimize the difference between the approximate and the exact solution, is still unknown but it can be determined from experimental or numerical data.As a first guess, C α was set to the average of f(α) over the interval [0, 1]: but it is likely that the optimum value is closer to the upper limit of f(α) in Eq. ( 9).The integral in Eq. ( 10) can be easily evaluated since it is equal to the Euler Beta function B(m + 1, n + 1), which is tabulated in standard mathematical handbooks and is available as a subroutine in numerical libraries.The value of C α in Eq. ( 10) is approximately 55 % of the upper limit value in Eq. ( 9) for typical values of m and n (see Table 1 and Fig. 2b where the horizontal line indicates the average value).Substitution of Eq. ( 7) in Eq. ( 6) yields the uncoupled energy equation: which can be made dimensionless, using a characteristic length scale L c and a characteristic time scale for heat conduction t char = ρcpL 2

C
λ .The length scale will for convenience be chosen to L/2 for the symmetric heating case.After normalizing the temperature with the adiabatic temperature rise ΔT adi = q c /c p the dimensionless energy conservation equation becomes: where the tilde sign indicates non-dimensional variables.The dimensionless parameters in Eq. ( 12) are the Damköhler number: which can be interpreted as the ratio of the heat conduction time scale to the reaction time scale, and the Arrhenius number:

R. Gebart
Ta = E R • ΔT adi (14) which can be interpreted as the ratio of the activation temperature E/R to the maximum temperature change that can occur if all reactants are consumed.
The boundary conditions to Eq. ( 12) are The initial condition for the temperature is T = Tu everywhere.The dimensionless formulation in Eq. ( 12) with its boundary conditions in Eq. (15) shows that the only parameters that can influence the behavior of the non-dimensional solution are the Damköhler number, the Arrhenius number and the dimensionless boundary temperature.
In the analysis the goal is to find under which conditions it is possible to obtain a quasi-steady balance between heat release from chemical reactions in the laminate and heat conduction to the boundaries.In the detailed analysis (see appendix) it is shown that the least stable case will be the steady-state solution when the laminate has had time to be heated through the thickness by heat transfer from the boundaries.

Criterion for thermal runaway; symmetric case
With symmetric heating, the solution must be symmetric, and the temperature gradient must be zero at the centerline.Hence, the solution from the mold wall to the centerline also describes the case with one heated and one insulated side (zero gradient) of the mold, which would be a good approximation for a sandwich laminate with a foam core.
The solution was derived as a perturbation series in a small parameter ∊ : where T0 = Ts is the steady state solution without chemical reactions ("frozen reactions") and ϑ is the first order correction, representing the response to heat release from the reactions.The small parameter ∊ was defined as: which greatly simplified the perturbation solution (see appendix for details).It can be interpreted as a typical amplitude of the (dimensionless) exothermal overshoot ( T − Ts ) that can be expected during cure, if the solution exists.It is important for the analysis that ∊ is small, and to estimate its magnitude it is helpful to rewrite it in dimensional terms as ∊ = T 2 s /ΔT adi T a , where T a = E/R and ΔT adi is the adiabatic temperature rise.The discussion can be restricted to resins that are used in structural polymer composites and composites with a high demand on dimensional stability.Those resins are designed to have a limited reactivity at the processing temperature to avoid excessive temperatures that can give rise to material degradation or warping of the finished part.This means that the Arrhenius factor in Eq. ( 6), which is the rate determining factor, must be relatively small and this is only possible if the ratio T a /T s has a large value.The remaining factor in the definition of ∊ is T s /ΔT adi , where ΔT adi depends on the total heat of reaction and the specific heat of the composite.The ratio is significantly smaller than 10 for most composites, which means that ∊ will be small if T a /T s > 10.For the resins in the result section, with the manufacturer's recommended processing temperature, T a /T s = 17.5 − 28.2 and ∊ ≈ 0.13 (see Table 3).For other commercial resins it is likely that ∊ will have a similar small value but it must be checked on a case-by-case basis.

Da e
− Ta Ts ∊ .The boundary conditions to Eq. ( 18) are: Eq. ( 18) can be integrated twice, using the first boundary condition, to yield: where the constant ϑ 0 represents the unknown value of ϑ at x = 0.The unknown constant ϑ 0 can be determined from the second boundary condition in Eq. ( 19), which gives the following transcendental equation for ϑ 0 : This equation has two possible solutions for 0 < β < 0.88 while for β > 0.88 there are no solutions.This means that for β > 0.88 the assumed steady-state balance between heat release and heat conduction is impossible.Hence, a criterion for the initiation of runaway reactions with symmetric heating can be defined as: Equation ( 22) can be recast as a runaway criterion for the Damköhler number by using the definition of β: When the thickness is below the critical thickness in Eq. ( 24) one can expect curing with low exotherm temperatures and above this thickness thermal runaway with higher temperatures will occur.The complete dimensional temperature profile at the critical Damköhler number can be obtained from Eq. ( 20) with ϑ 0 = 1.186 and β = 0.88: The shape of this profile is shown in Fig. 3.The corresponding temperature overshoot at the centerline is: Since ∊ is a small number one may conclude that the peak temperature overshoot in a stably cured laminate, i.e. with Da ≤ Da crit , will be a small fraction of the adiabatic temperature rise.The dimensional time to reach complete cure at the centerline can be estimated by substituting the result from Eq. ( 26) and Eq. ( 7) in Eq. ( 6): The righthand side of Eq. ( 27) is a constant which means that it can be integrated to: where the last term is an approximation of the time it takes to heat the laminate and reach stationary conditions and where it has been assumed that the initial value of α is negligible and the final value is α = 1.

Results
A literature search for suitable validation experiments found two independent studies for the symmetric case.The first assessment case is a carbon fiber-epoxy laminate with a fiber volume fraction of 60 % [5].In this paper, several cure cycles with one to three constant temperature stages (called temperature dwells) are investigated but only the one dwell case is relevant for the present study.
The second assessment case is a glass fiber-epoxy laminate with 40 % fiber volume fraction from a study by Parthasarathy et al. [15].In this work, the cure cycle is first optimized with numerical simulations to avoid high exotherm temperatures followed by experiments to validate that the optimized cure cycle gives the desired result.Experiments and simulations were performed for four different laminate thicknesses but only the two thickest that gave a significant exotherm are of interest for the present study.
The cure kinetic model parameters that were used for the analysis are listed in Table 1 and were taken from the original papers [5,15].The corresponding adiabatic temperature rise is listed in Table 2 together with the effective thermal properties for the two assessment cases.Table 1 and 2 also contains model parameters for a third material, Fiberite 976, to give additional information about the variation that can be expected in the model parameters.Notice that the kinetic model parameters for the two assessment cases in Table 1 are quite different, with a significantly higher activation energy and pre-exponential constant in Parthasarathy et al., but that the exponents m (0.3 and 0.5) and n (1.7 and 1.88) are relatively close to each other in both cases.Moreover, the heat of reaction and adiabatic temperature rise in Li et al. [5] is close to twice that in Parthasarathy et al. [15] while the thermal conductivity is about three times larger.
The procedure for using the present theory is illustrated in Table 3 where the various parameters in the solution are computed step by step.First, the Damköhler number, the Arrhenius number and the perturbation parameter ∊ are computed from the known input parameters.In this step a check should be made that the perturbation parameter is small.Secondly, the critical Damköhler number and the maximum allowable thickness are computed from the dimensionless mold temperature Ts and the Arrhenius number Ta .Finally, the input parameters are compared to the criticality conditions Da crit and L max to check the margin to thermal runaway.
In the study by Li et al. [5] the laminate thickness was 40 mm, and the boundary temperature was set to 131 • C, which resulted in a numerically computed peak temperature at the laminate centerline of 200 • C, i.e. an exothermal peak of about 69 • C.This exothermal Fig. 3. Normalized temperature excess (ϑ) from the centerline at x = 0 to the wall at x = 1, for critical conditions with β = 0.88 and ϑ 0 = 1.186 (see eq.( 20)).

Table 2
Equivalent homogeneous laminate properties for three prepreg laminates.The data for MXB7701/7781 and EPON 862/W was taken from Parthasarathy et al. [15] and Li et al. [5] respectively.The thermal conductivity for Fiberite 976 was assumed to be the same as that in Bard et al. [10] while the rest of the properties were estimated with a fit to DSC-data from Dusi et al. [12].temperature rise is about 34 % of the adiabatic temperature rise, which indicates that the cure process is close to thermal runaway.The present theory confirms this assumption since the real laminate thickness of 40 mm is very close to the maximum allowable thickness of 40.1 mm (see Table 3).Hence, the present approximate theory is in excellent agreement with the detailed numerical simulations by Li et al.
Using the present theory, it is also possible to determine the maximum mold temperature that can be used without thermal runaway, and it is only 0.5 K higher than the value used by Li et al. [5].This shows, as expected, that the safety margin to runaway conditions is very small when the optimized boundary temperature from Li et al. [5] is used.The centerline peak temperature in the laminate, calculated with the symmetrical solution in Eq. (25), was 30 K above the boundary temperatures while the simulation result in Li et al. is about two times higher than this (see Table 3).The reason for this discrepancy most likely comes from the approximation to the kinetic model that was done in the present paper, which underestimates the maximum reaction rate.
In the study by Parthasarathy et al. [15] the two thickest laminates, 14.4 mm and 21.6 mm thick, were found to exhibit significant exotherms when they were cured with the manufacturers recommended cure cycle at 120 • C. For the 14.4 mm laminate the experimental peak temperature at the center of the laminate was 25.5 K higher than the boundary temperature and for the thicker laminate it was more than twice this value [15].This can be compared with the adiabatic temperature rise, which is 103 K based on the data in Parthasarathy et al., which shows that thermal runaway must have occurred in the thicker laminate while the thinner laminate must be close to the limit of runaway.This was also confirmed by the present theory, which showed that the Damköhler number for the thinner laminate is 25 % below the critical value predicted by the present theory, while the thicker laminate is 69 % above the critical value in Eq. (23) (see Table 3).The theoretical maximum thickness that can be cured without thermal runaway (see Eq. ( 24)) at the same cure temperature (393 K) is 16.6 mm, which gives a thickness margin of 15 % for the thinner laminate but is exceeded by 30 % for the thicker laminate.Turning the attention to the thicker laminate, it is interesting to use the theoretical results to investigate how much the boundary temperature must be lowered to avoid thermal runaway.After a few hand iterations it was found that a relatively modest decrease of the mold temperature from 393 K to 385 K would prevent thermal runaway in the thicker laminate.Based on this result, it seems likely that it would be possible to achieve complete cure by an initial pre-curing of the thicker laminate at 385 K until a high degree of cure is achieved, followed by a post-cure at the manufacturers recommended cure temperature of 393 K to reach complete cure.

Discussion
The solution that was developed in this paper was partially validated against literature data, with promising result.However, for a full validation there is a need for more numerical or experimental data near critical conditions and for several different materials.In an extended validation study, the model constant C α should be optimized to give an accurate prediction of both the critical Damköhler number and the peak temperature in the laminate.
The theoretical solution in this paper should be seen as a complement, not a replacement, to detailed computer models since the computer models solve the full equations without approximations and can simulate much more complicated cure cycles [5].They can also be coupled to process control systems that can adjust the processing temperature in real time [15].However, the present theoretical solution has the advantage that it gives an excellent overview of the dependence on the process and material parameters that is lacking in a computer model.In principle, it is possible to obtain the same information from a computer model, but it would require a significant effort.As an example, a map over the dimensionless peak temperature in a symmetrically heated laminate versus the three dimensionless numbers that governs its magnitude (Da, Ta and Ts ) would require 1000 independent simulations if each parameter is resolved with ten steps.
A key simplification in the present analysis was the approximation of the function f(α) in the reaction rate equation in Eq. ( 6) with a constant value C α .Without this simplification it would not have been possible to use the perturbation method that led to the criterion for thermal runaway.The good agreement between the theoretical runaway criterion and the assessment cases discussed above gives support to the assumption that the fast exponential temperature response of the chemical reaction rate is more important for thermal runaway than the slow response to changes in the degree of cure.

Conclusions
A theoretical solution, with symmetric heating, has been developed for the temperature in thermoset composite laminates during cure.The analysis showed that three dimensionless parameters, the Damköhler number, the Arrhenius number, and the dimensionless boundary temperature, govern the behavior of the solution.
An important step in the solution was that the function f(α) in the reaction rate equation (Eq.( 6)) was approximated as a constant C α , which decoupled the energy equation and the degree of cure equation.The constant was as a first guess set to the average of f(α).However, it is likely that the accuracy of the solution can be improved by optimizing C α against detailed numerical solutions of the full equations.
The solution was derived as a perturbation series in a small parameter ∊ that was found to have a natural definition that simplified the analysis.The analysis showed that a quasi-stationary balance between heat release and heat conduction in the laminate is only possible if the Damköhler number is below a critical value (see Eq. ( 23) that depends on the Arrhenius number and the dimensionless wall temperature.However, when the Damköhler number exceeds the critical value, thermal runaway with a high peak temperature will occur.
The thermal runaway criterion based on the Damköhler number can, by using the definition of the Damköhler number, be recast as a criterion for the maximum laminate thickness that can be cured without thermal runaway (see Eq. ( 24)).The maximum laminate thickness will be different for different material combinations, because of differing Table 3 Process parameters and stability analysis results for symmetric heating using thermal properties and kinetic model parameters from Table 1 and Table 2.The thickness of the laminate L, the cure temperature T s , and the peak temperature in the laminate ΔT peak are taken from Parthasarathy et al. [15] for MXB7701/7781 and from Li et al. [5] for EPON 862/W.β crit e ϑ0 )1 2 , which gives ϑ 0 ≈ 1.186.The resulting solution for ϑ at critical conditions becomes:

Fig. 2 .
Fig. 2. Graphical illustration of the cure model dα dt = Ae − E RT α m (1 − α) n using the model parameters from Parthasarathy et al .[15] (see Table 1): a) Surface plot of reaction rate versus temperature and degree of cure; b) Reaction rate versus degree of cure at T = 370 K.The horizontal line represents the average value of the curve; c) Reaction rate versus temperature at α=0.2.
critical Damköhler number only depends on the dimensionless boundary temperature Ts and the Arrhenius number Ta .The critical Damköhler number can be used to define a maximum thickness criterion:

Table 1
[12]tic model parameters for three types of epoxy resin determined by fitting to DSC data.The parameters for MXB7701/7781 and EPON 862/W were taken from Parthasarathy et al.[15]and Li et al.[5]respectively, while the parameters for Fiberite 976 were determined by fitting to DSC-data from Dusi et al.[12].