Safe Curing Limits of Thick Composite Shells

1. Introduction

The curing of thick fibre-reinforced polymer laminates (FRP) is challenging due to the strong exothermic nature of thermoset cross-linking together with the low thermal conductivity of polymer composites. When the heat generated by the curing reaction cannot be dissipated efficiently, large temperature gradients can develop within the laminate. These gradients may lead to a heterogeneous degree of cure, accumulation of residual stresses, matrix degradation, or even delamination. During cure, chemical shrinkage and modulus development further induce residual stresses and warpage [1], which in turn can initiate transverse cracking and delamination [2]. Thermal control is therefore essential for both dimensional stability and structural integrity.

Industrial cure cycles are commonly developed on the basis of experience with thin laminates, where heat removal is relatively effective and strong temperature gradients do not arise. However, as laminate thickness increases or when rapid cure cycles are employed, the balance between heat generation and heat dissipation may be disrupted. Under such conditions, manufacturers rely on numerical simulations to predict thermal behaviour during curing. Fully coupled thermal-kinetic models, typically solved using finite-difference, finite-volume, or finite-element methods, can reproduce thermal gradients and overshoot for a wide range of composite systems [3,4,5,6]. Although these models are accurate and versatile, they can take a lot of time to apply in practice, be sensitive to the uncertainty of the kinetic parameters, and offer limited physical insight into the stability boundary that separates safe processing from thermal runaway.

Analytical approaches can provide additional insight by revealing the underlying physics and yielding explicit relations between overheating, material parameters, and geometry. However, analytical predictions of thermal runaway in composite curing remain limited. Previous work established a perturbation-based runaway criterion for flat laminates [7]. Farjas et al. [8] obtained analytical critical-thickness conditions for planar geometries with convective boundaries. In parallel, process-orientated studies [6,9] have used these flat-plate criteria to define safe operating windows in advanced manufacturing strategies. Yet, no analytical framework currently generalises these criteria to curved shells or expresses the critical condition in a form that separates geometric/boundary effects from chemistry/temperature sensitivity.

In this work, we derive a stability criterion for thermal runaway during the cure of thick composite shells by reducing the full 3D thermo-kinetic problem to a locally one-dimensional through-thickness model evaluated pointwise over the shell surface. This reduction yields a closed-form critical condition for the Damköhler number in which geometry and boundary heat transfer enter only through a single stability factor, while chemistry and processing temperature enter through an explicit Arrhenius scaling. We compute the geometric stability factor from a nonlinear boundary-value problem and represent it by compact response surfaces for symmetric and asymmetric Robin boundary conditions, allowing rapid stability mapping over arbitrary shell geometries and Biot numbers. We validate the resulting criterion against high-precision transient simulations and a quasi-three-dimensional flux assessment, which confirms that transverse gradients dominate for industrially relevant curvatures and that the locally 1D approximation is conservative. Finally, we introduce an objective numerical classification of runaway onset based on temperature derivatives, providing a consistent way to locate the critical threshold in transient simulations and suggesting a route to real-time process monitoring.

The paper is organised as follows. Section 2 formulates the locally 1D shell model, introduces the perturbation framework, and derives the stability condition and its parameter separation. Section 3 validates the criterion against fully coupled transient simulations and examines the transient behaviour near criticality. Section 4 discusses the implications for process design and the limits of the underlying assumptions. Appendices A and B provide the derivation details and numerical procedures, respectively.

2. Theory and Numerical Methods

We formulate the governing thermo-kinetic problem for a general thick-walled shell in body-fitted coordinates and derive a thermal-runaway stability criterion. The formulation starts from the full three-dimensional energy balance with standard cure kinetics, which makes the geometric scope explicit: locally, any smooth shell is characterised by its two principal curvatures. We then introduce a controlled reduction to a locally one-dimensional through-thickness model, yielding a conservative stability problem that can be evaluated pointwise over the surface. Details of the perturbation derivation are given in Appendix A; numerical procedures are summarised here and documented in Appendix B.

2.1. Problem Formulation

The theory was developed for a general curved shell of half-thickness w whose mid-surface is described by two principal curvatures ${\kappa }_1$ and ${\kappa }_2$. A body-fitted coordinate system $(u,v,z)$ is defined with orthogonal principal-curvature coordinates $(u,v)$ on the mid-surface and a thickness coordinate z that measures the distance from the mid-surface along the surface normal (see Figure 1). Any point in the laminate can be expressed as

$$\mathbf{R}(u,v,z)=\mathbf{r}(u,v)+z\mathbf{n}(u,v),$$

where $z\in [-w,w]$.

Since the goal is to develop a theory that is valid for shells with arbitrary shape, the analysis will start from the most general formulation of the energy conservation equation combined with a standard autocatalytic cure kinetic model [10]:

$$\begin{array}{cc}\hfill \rho c_p{\displaystyle \frac{\partial T}{\partial t}}& =\nabla ·(\mathsf{\Lambda }\nabla T)+\rho q_c{\displaystyle \frac{\partial \alpha }{\partial t}},\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial \alpha }{\partial t}}& =Aexp\left(-{\displaystyle \frac{E}{RT}}\right){\alpha }^m{(1-\alpha )}^n,\hfill \end{array}$$

(2)

where T is the temperature, $\alpha$ is the degree of cure, $\mathsf{\Lambda }$ is the anisotropic thermal conductivity tensor and all other parameters have their standard meaning (see e.g. [7]).

The general problem has no known analytical solution, but it can be reduced to a one-dimensional equation along the thickness direction if in-plane conduction is neglected. In-plane conduction mainly redistributes heat laterally and therefore acts as a stabilising mechanism with respect to the risk of thermal runaway. Neglecting in-plane conduction thus yields a conservative (worst-case) approximation in which each surface point behaves as an independent through-thickness "thermal pillar". Crucially, because these thermal pillars are decoupled from their surroundings and depend strictly on the local principal curvatures, this single 1D formulation remains valid for every point on any smooth geometric shape. This decoupling allows the critical limit for thermal runaway to be precomputed over a broad range of curvatures and represented by compact response surfaces for rapid stability mapping on complex geometries.

Under this assumption, and using standard differential-geometry identities for orthogonal surface coordinates with principal curvatures ${\kappa }_1$ and ${\kappa }_2$, the energy equation (1) reduces to

$$\rho c_p{\displaystyle \frac{\partial T}{\partial t}}=\lambda \left[{\displaystyle \frac{{\partial }^{2}T}{\partial z^2}}+\left({\displaystyle \frac{{\kappa }_1}{1+{\kappa }_{1}z}}+{\displaystyle \frac{{\kappa }_2}{1+{\kappa }_{2}z}}\right){\displaystyle \frac{\partial T}{\partial z}}\right]+\rho q_c{\displaystyle \frac{\partial \alpha }{\partial t}},$$

(3)

where $\lambda$ is the transverse (through-thickness) conductivity associated with $\mathsf{\Lambda }$. To accommodate heated tooling, convection heating, or insulating boundaries within a unified mathematical framework, a Robin (convective) boundary condition is applied on both surfaces:

$$-\lambda {\displaystyle \frac{\partial T}{\partial n}}=h_{\mathrm{eff}}(T_{\mathrm{surf}}-T_{\infty }).$$

where n is the outward unit normal, $h_{\mathrm{eff}}$$\left[\mathrm{W}{\mathrm{m}}^{-2}{\mathrm{K}}^{-1}\right]$ is the heat transfer coefficient, $T_{\infty }$$\left[\mathrm{K}\right]$ is the reference temperature in Newton’s law of cooling [11] and $T_{\mathrm{surf}}$$\left[\mathrm{K}\right]$ is the temperature at the surface of the part.

To reveal the governing non-dimensional groups and reduce the total number of independent process parameters, all variables are scaled using the same procedure as in [7]:

• through-thickness coordinate scaled with the laminate half-thickness $w\to \tilde{x}=z/w$,
• time with the thermal-diffusion time from the boundary to the centreline $\tau =w^2\rho c_p/\lambda \to \tilde{t}=t/\tau$,
• temperature with the adiabatic temperature rise $\Delta T_{\mathrm{adi}}=q_c/c_p\to \tilde{T}=T/\Delta T_{\mathrm{adi}}$,
• and curvature with the thickness $w\to {\sigma }_i={\kappa }_{i}w$.

With these scalings, the governing equations take the dimensionless form:

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial \tilde{T}}{\partial \tilde{t}}}& ={\displaystyle \frac{{\partial }^2\tilde{T}}{\partial {\tilde{x}}^2}}+\left({\displaystyle \frac{{\sigma }_1}{1+{\sigma }_1\tilde{x}}}+{\displaystyle \frac{{\sigma }_2}{1+{\sigma }_2\tilde{x}}}\right){\displaystyle \frac{\partial \tilde{T}}{\partial \tilde{x}}}+{\displaystyle \frac{\partial \alpha }{\partial \tilde{t}}},\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial \alpha }{\partial \tilde{t}}}& =Daexp\left(-{\displaystyle \frac{A_r}{\tilde{T}}}\right){\alpha }^m{(1-\alpha )}^n,\hfill \end{array}$$

(5)

with dimensionless boundary conditions

$$\begin{array}{cc}\hfill {\left.{\displaystyle \frac{\partial \tilde{T}}{\partial \tilde{x}}}\right|}_{\tilde{x}=-1}& =+B{i}_{\mathrm{d}}\left({\tilde{T}}_{\mathrm{surf},\mathrm{d}}-{\tilde{T}}_{\infty }\right),\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill {\left.{\displaystyle \frac{\partial \tilde{T}}{\partial \tilde{x}}}\right|}_{\tilde{x}=+1}& =-B{i}_{\mathrm{u}}\left({\tilde{T}}_{\mathrm{surf},\mathrm{u}}-{\tilde{T}}_{\infty }\right),\hfill \end{array}$$

(7)

where the Damköhler number $Da$, the Arrhenius number $A_r$ and the Biot numbers are defined as:

$$Da={\displaystyle \frac{\rho c_{p}A{w}^2}{\lambda }},A_r={\displaystyle \frac{E}{R\Delta T_{\mathrm{adi}}}}B{i}_{\mathrm{d}},\mathrm{u}={\displaystyle \frac{{\left(h_{\mathrm{eff}}\right)}_{\mathrm{d}},\mathrm{u}w}{\lambda }}.$$

In spite of the simplification, the system in Eqs. (4) - (6) has no known general solution, but it can be solved with the approximate perturbation method discussed in the next section.

2.2. Perturbation Solution

From this point onwards, the tilde notation is dropped and all variables are dimensionless unless otherwise stated.

To make the coupled thermal–kinetic problem tractable, two simplifying steps are introduced (see Appendix A for details).

First, the kinetic term is frozen at its maximum value, replacing

$$f\left(\alpha \right)={\alpha }^m{(1-\alpha )}^n$$

with a constant, $C_{\alpha }=max\left(f\left(\alpha \right)\right)$ . This corresponds to the most heat-intensive stage of cure and therefore yields a conservative (worst-case) approximation of the runaway risk. With this approximation, the heat-generation term no longer depends on the evolving degree of cure, which decouples the energy equation from the cure kinetics.

Second, the analysis focuses on the steady-state temperature distribution that would arise if the laminate was exposed to its highest admissible background temperature $T_{\infty }$. If the part can reject reaction heat under these conditions, it will remain stable through the transient leading up to the stationary state. The temperature is therefore written as a small correction to this reference state,

$$T=T_{\infty }+\epsilon \vartheta$$

which leads to a static through-thickness problem for the perturbation temperature profile $\vartheta \left(x\right)$ (see Appendix A):

$${\vartheta }^{\prime \prime }+\left({\displaystyle \frac{{\sigma }_1}{1+{\sigma }_{1}x}}+{\displaystyle \frac{{\sigma }_2}{1+{\sigma }_{2}x}}\right){\vartheta }^{\prime }+\beta e^{\vartheta }=0.$$

(8)

Here the dimensionless parameter $\beta$

$$\beta =Da{C}_{\alpha }exp\left(-{\displaystyle \frac{A_r}{T_{\infty }}}\right){\displaystyle \frac{1}{\epsilon }}$$

(9)

measures the balance between heat generation and heat removal.

The corresponding boundary conditions involve only the Biot numbers and are independent of the reference temperature $T_{\infty }$, a direct consequence of expanding the temperature field around this reference state:

$${\vartheta }^{\prime }(-1)-B{i}_{\mathrm{d}}\vartheta (-1)=0,{\vartheta }^{\prime }\left(1\right)+B{i}_{\mathrm{u}}\vartheta \left(1\right)=0.$$

The perturbation problem admits a steady solution only if $\beta$ is below a certain threshold. When $\beta$ exceeds this limit, the steady-state temperature profile can no longer exist and the temperature becomes self-accelerating. Because the parameters governing the perturbation equation are the curvatures $({\sigma }_1,{\sigma }_2)$ and the Biot numbers $(B{i}_d,B{i}_u)$, the critical value

$${\beta }_{\mathrm{crit}}={\beta }_{\mathrm{crit}}({\sigma }_1,{\sigma }_2,B{i}_d,B{i}_u)$$

is determined entirely by geometry and surface heat transfer.

Once ${\beta }_{\mathrm{crit}}$ is known, the corresponding critical Damköhler number follows directly from Equation (9):

$$D{a}_{\mathrm{crit}}={\beta }_{\mathrm{crit}}{\displaystyle \frac{\epsilon }{C_{\alpha }}}exp\left({\displaystyle \frac{A_r}{T_{\infty }}}\right),$$

(10)

This expression separates the effects of geometry and heat transfer (contained in ${\beta }_{\mathrm{crit}}$) from those of chemistry and processing temperature. The critical Damköhler number can be recast as a criterion for the maximum safe half-thickness of the laminate for given process conditions:

$$w_{\max }=\sqrt{{\displaystyle \frac{\lambda }{\rho c_p}}{\displaystyle \frac{1}{A}}{\displaystyle \frac{1}{C_{\alpha }}}{\beta }_{\mathrm{crit}}{\displaystyle \frac{T_{\infty }^2}{A_r}}exp\left({\displaystyle \frac{A_r}{T_{\infty }}}\right)}.$$

(11)

The perturbation BVP in Equation (8) is still too complicated to solve analytically. However, a semi-analytical solution can be developed by solving the BVP numerically for a range of curvatures and Biot numbers and fitting an analytical expression to the numerical results. The determination of ${\beta }_{\mathrm{crit}}$ for each parameter set is done by an iterative search for the upper limit for valid solutions to the stationary problem (see Appendix B for details). The resulting map of ${\beta }_{\mathrm{crit}}$ versus curvatures and Biot numbers is used to fit analytical expressions. The resulting semi-analytical result can then be used as the basis for the universal thermal runaway criterion.

The perturbation solution also predicts the peak temperature rise at the stability limit:

$$\Delta T_{\mathrm{crit}}={\vartheta }_{max,\mathrm{crit}}\epsilon ,$$

where ${\vartheta }_{max,\mathrm{crit}}=ma{x}_x\left({\vartheta }_{\mathrm{crit}}\left(x\right)\right)$ is an $\mathcal{O}\left(1\right)$ shape factor that only depends on $({\sigma }_1,{\sigma }_2,B{i}_d,B{i}_u)$.

For a flat laminate with Dirichlet boundary conditions ${\vartheta }_{max,\mathrm{crit}}=1.186$ [7]. However, for other curvatures, it is difficult to determine ${\vartheta }_{\max },\mathrm{crit}$ accurately from the BVP because even a very small uncertainty in ${\beta }_{\mathrm{crit}}$ is strongly amplified in the peak-temperature estimate. This sensitivity makes the stepping method used to locate ${\beta }_{\mathrm{crit}}$ (see Appendix B) poorly suited for extracting $\Delta T_{\mathrm{crit}}$ directly from the perturbation equation.

Instead, the peak temperature at critical conditions is obtained from the full transient numerical model described in Appendix B, which exhibits much weaker sensitivity to small errors in the critical limit. The resulting values of $\Delta T_{\mathrm{crit}}$ are then fitted to the analytical scaling law above, ensuring consistency between the perturbation theory and the exact numerical solution.

In summary, the perturbation method offers an intuitive description of thermal-runaway onset: geometry and surface heat transfer determine the shape of the stability boundary, while chemistry influences the critical Damköhler number through a simple exponential dependence. The resulting stability criterion in Equation (10) is concise, physically transparent, and highly accurate, as confirmed by comparison with fully coupled numerical simulations in the Results section.

The discussion that follows examines the shape of the ${\beta }_{\mathrm{crit}}$ surface for two practically important heating configurations: (i) symmetric heating, where both surfaces have the same Biot number, and (ii) asymmetric heating, where one surface is strongly coupled to the surroundings and the other has a variable heat-transfer coefficient. These two cases cover the vast majority of industrial curing setups and illustrate the essential trends in how curvature and boundary conditions influence stability.

2.2.1. Symmetric Heating

Figure 2 shows the stability surface ${\beta }_{\mathrm{crit}}({\sigma }_1,{\sigma }_2)$ for the case of symmetric heating, here illustrated with $Bi=1000$. Changing the Biot number shifts the overall height of the surface but does not alter its geometric shape. What is immediately apparent is the symmetry about the diagonal ${\sigma }_1={\sigma }_2$. This reflects the exchange symmetry of the governing equations: interchanging the principal curvatures does not change the physical problem, and therefore ${\beta }_{\mathrm{crit}}({\sigma }_1,{\sigma }_2)={\beta }_{\mathrm{crit}}({\sigma }_2,{\sigma }_1)$.

A second symmetry appears along the opposite diagonal, where $({\sigma }_1,{\sigma }_2)\to (-{\sigma }_1,-{\sigma }_2)$. Under symmetric heating, flipping both curvatures simultaneously and changing the direction of the thickness coordinate x corresponds to reversing the shell orientation while leaving the heat-transfer conditions unchanged. Because the perturbation equation depends on the curvatures only through the combination $\left({\displaystyle \frac{{\sigma }_1}{1+{\sigma }_{1}x}}+{\displaystyle \frac{{\sigma }_2}{1+{\sigma }_{2}x}}\right)$ the stability limit is unaffected by this transformation. Together, these two symmetries explain the characteristic appearance of the stability surface in Figure 2.

These symmetries strongly constrain the functional form of any analytical expression for ${\beta }_{\mathrm{crit}}$. By rewriting the surface in terms of the curvature invariants $S={\sigma }_1+{\sigma }_2$ and $P={\sigma }_1{\sigma }_2$, both symmetry conditions are satisfied automatically. This makes it possible to approximate the entire stability surface using a compact polynomial in S and P, which closely reproduces the numerical data:

$${\beta }_{\mathrm{crit}}(S,P,Bi)=c_0\left(Bi\right)+a{S}^2+bP,$$

(12)

The coefficients $a=0.121$ and $b=0.427$ are obtained from a least-squares fit to the numerical data. While these parameters exhibit a weak $Bi$-dependency for low Biot numbers $(Bi<5)$, they are treated here as constants to maintain the analytical simplicity of the universal criterion. For typical laminate curvatures $(-0.4<{\sigma }_{1,2}<0.4)$, this approximation introduces an error of less than 1% for $Bi\ge 5$, and remains highly accurate even at $Bi=1$ with a maximum error of approximately 2.2%. Consequently, the curvature effects are decoupled from the surface heat transfer in the model, and only the offset term $c_0$ is required to vary with $Bi$:

$$c_0\left(Bi\right)={\beta }_{\mathrm{crit},\mathrm{Dirichlet}}\left({\displaystyle \frac{Bi}{Bi+K}}\right),$$

(13)

where $K=2.172$ is a fitting parameter. The reference value ${\beta }_{\mathrm{crit}},\mathrm{Dirichlet}=0.87846$ represents the analytical stability limit for a flat slab $({\sigma }_1={\sigma }_2=0)$ under Dirichlet boundary conditions $(Bi\to \infty )$ [7]. Equation (13) captures the dominant $Bi$-dependence with high fidelity, reflecting the rapid boundary-layer saturation observed as the process moves toward a diffusion-limited regime (see Figure 3).

Taken together, the symmetric-heating results provide a clear physical picture of how geometry and heat-transfer conditions influence the onset of thermal runaway. Increasing curvature consistently lowers the stability margin: as the shell becomes more curved, heat is removed less efficiently through the thickness, so runaway occurs under milder chemical conditions, embodied by the Arrhenius number and the processing temperature $T_{\infty }$. In contrast, the Biot number affects the overall sensitivity more directly and independently of curvature; a larger Biot number simply means better heat transfer at the surfaces, raising the entire stability surface uniformly. The predicted peak temperature at the stability limit shows only a weak dependence on curvature but is strongly governed by the small parameter $\epsilon$, and thus by the Arrhenius number and $T_{\infty }$.

Beyond these points, several additional conclusions might be of practical interest to engineers. First, the symmetric case highlights that geometric effects act primarily as a penalty, not as a source of interaction: synclastic (dome-like) and anticlastic (saddle-like) curvatures follow the same general trends, meaning the stability limit can be characterised using only the scalar invariants S and P. Second, the rapid saturation of $c_0\left(Bi\right)$ shows that improving surface heat transfer beyond a moderate Biot number yields diminishing returns, an important observation for tooling and process design where changes in convection or tool contact may be costly. Third, the smoothness and symmetry of the ${\beta }_{\mathrm{crit}}$ surface imply that manufacturing tolerances in curvature are unlikely to produce abrupt stability changes; instead, stability degrades gradually as curvature increases, which improves predictability in engineering assessments. Finally, the weak curvature dependence of $\Delta T_{\mathrm{crit}}$ means that temperature overshoot predictions can be made with reasonable confidence even before curvature-specific calibration, whereas safe-thickness predictions require using the full stability surface.

These conclusions together provide engineers with an intuitive and practically useful understanding of how geometry and heat-transfer efficiency shape the stability limits for thick composite shells under symmetric heating conditions.

2.2.2. Asymmetric Heating

For asymmetric heating conditions, where one surface is cooled more efficiently than the other $(B{i}_u=1000,B{i}_d=Bi<1000)$, the coordinate-flip symmetry of the system is broken. In the symmetric case, the stability limit was invariant under the simultaneous transformation of reversing the thickness coordinate $(x\to -x)$ and flipping the signs of the principal curvatures $(({\sigma }_1,{\sigma }_2)\to (-{\sigma }_1,-{\sigma }_2))$. In the asymmetric case, however, this transformation is no longer a symmetry of the problem because it swaps the positions of the well-cooled and poorly-cooled surfaces relative to the direction of curvature.

Consequently, the stability limit becomes sensitive to whether the laminate is curving "towards" or "away from" the side with inferior heat extraction. This manifests as a distinct "tilt" and distortion of the ${\beta }_{\mathrm{crit}}$ stability surface (see Figure 4), particularly at low Biot numbers (Bi<10). In this regime, the onset of thermal runaway depends strongly on the orientation of the shell relative to the thermal bottleneck at the less-cooled surface. As the cooling efficiency increases toward $Bi\approx 10$, this sensitivity to orientation diminishes, and the stability surface eventually recovers the symmetric characteristics of the $Bi\to \infty$ limit.

To capture this physical distortion mathematically, the semi-analytical fitting polynomial must be modified. Specifically, an additional linear term in the curvature sum S is required to account for the first-order sensitivity to curvature direction that appears when the cooling is unbalanced:

$${\beta }_{\mathrm{crit}}(S,P,Bi)=c_0\left(Bi\right)+L\left(Bi\right)S+a\left(Bi\right)S^2+b\left(Bi\right)P$$

(14)

In this asymmetric regime, all four coefficients $(c_0,L,a,b)$ exhibit a dependency on the Biot number. The linear coefficient $L\left(Bi\right)$ explicitly quantifies the "tilt" of the surface, while the remaining terms describe the evolution of the surface height and curvature sensitivity. These dependencies are accurately captured by the following fitting functions:

$$\begin{array}{cc}\hfill c_0\left(Bi\right)& ={\displaystyle \frac{c_{1}Bi+c_2{c}_3}{Bi+c_3}},\hfill \end{array}$$

(15a)

$$\begin{array}{cc}\hfill L\left(Bi\right)& =L_{1}exp\left(-{\displaystyle \frac{Bi}{L_2}}\right),\hfill \end{array}$$

(15b)

$$\begin{array}{cc}\hfill a\left(Bi\right)& =a_1+a_{2}exp\left(-{\displaystyle \frac{Bi}{a_3}}\right),\hfill \end{array}$$

(15c)

$$\begin{array}{cc}\hfill b\left(Bi\right)& =b_1-b_{2}exp\left(-{\displaystyle \frac{Bi}{b_3}}\right).\hfill \end{array}$$

(15d)

with the optimised coefficients listed in Table 1. These functions reproduce the ${\beta }_{\mathrm{crit}}$ surface for all fitted data with errors well below one percent (see Figure 5).

The asymmetric-heating results reveal several important differences compared to the symmetric case. The most prominent is that unequal boundary conditions introduce a directional bias into the heat-removal capability of the shell. This appears as a characteristic “tilt” of the stability surface, captured by the coefficient $L\left(Bi\right)$, which shifts the stability margin depending on whether the dominant curvature faces the hotter or the better-cooled surface. For engineers, this has a clear implication: orientation of the heating matters. A curved part with one side poorly cooled is inherently more vulnerable to runaway when its convex face is exposed to the weaker boundary condition.

Despite this asymmetry, many trends remain consistent. Mean curvature $(S/2)$ and Gaussian curvature $\left(P\right)$ still control the overall penalty through the coefficients $a\left(Bi\right)$ and $b\left(Bi\right)$, and these vary smoothly with $Bi$. Once the poorly cooled face reaches moderate heat-transfer efficiency, the influence of asymmetry decreases rapidly, mirroring the saturation behaviour discussed above. This means that in many practical setups, especially autoclave or oven curing where one surface is strongly coupled to the environment, only modest improvements to the weaker surface are needed to bring the system close to the more stable symmetric-heating case.

Another valuable insight is that the asymmetric case does not introduce any abrupt or unexpected behaviours: the stability surface remains smooth, monotonic, and predictable across the $({\sigma }_1,{\sigma }_2)$ domain. Thus, while asymmetric heating reduces the overall stability margin, it does so in a way that is systematic and quantifiable, enabling engineers to evaluate the effect of tooling contact, insulation, layup orientation, and part placement without resorting to full numerical simulations. In all, the asymmetric-heating analysis completes the picture by showing how unequal boundary conditions shift, but do not fundamentally change, the underlying geometric structure of the thermal-runaway limit.

2.2.3. General Boundary Conditions

While the symmetric and strongly asymmetric configurations discussed above capture the essential physical phenomena, a general heating scenario may involve two arbitrary, unequal Biot numbers. Conceptually, the stability surface for such a general case can be understood through a superposition of the two established behaviors. If one initially assumes a symmetric state based on the higher of the two Biot numbers, the overall elevation of the ${\beta }_{\mathrm{crit}}$ surface is correspondingly reduced. Subsequently reducing the remaining surface’s Biot number to its true, lower value breaks the symmetry, introducing the characteristic "tilt" and, at sufficiently low values, the flattening observed in the asymmetric case.

Rather than developing an exhaustive semi-analytical formulation for all possible $(B{i}_u,B{i}_d)$ combinations, a computational approach is recommended for general conditions. The software provided in the supplementary material (detailed in Appendix B) allows the user to input any arbitrary pair of Biot numbers to automatically generate a custom polynomial fit for ${\beta }_{\mathrm{crit}}$. Once this specific fit is generated, the thermal stability of the part can be straightforwardly evaluated using the theoretical framework $(D{a}_{\mathrm{crit}},\Delta T_{\max },\mathrm{and}{w}_{\max })$ established in the preceding sections.

2.3. Numerical Methods (Summary)

The perturbation boundary-value problem for ${\beta }_{\mathrm{crit}}$ is solved using MATLAB’s bvp4c with a continuation strategy in $\beta$ (details in Appendix B). The fully coupled transient 1D thermo-kinetic model used for validation of the perturbation model is discretized by a second-order cell-centered finite-volume method (method of lines) and integrated with MATLAB’s ode15s using a relative tolerance of ${10}^{-6}$ (Appendix B). A grid-convergence study with Richardson extrapolation indicates that the production grid yields relative errors below $2\times {10}^{-4}$ for the peak temperature. Thermal runaway in the numerical solution is classified objectively using the sign of the peak second time derivative of a virtual sensor temperature at the hottest node, and $D{a}_{\mathrm{crit}}$ is located using MATLAB’s fzero solver (Appendix B). The complete MATLAB (R2025a) source code used to generate the results is provided as supplementary material.

3. Results

This section validates the perturbation-based predictions against high-precision numerical solutions of the fully coupled thermal-kinetic problem (Eqs. 4 - 5). The standalone transient solver (Appendix B) resolves the through-thickness energy balance with Robin boundary conditions and autocatalytic kinetics without the simplifying assumptions used in the perturbation analysis. A grid-convergence study with Richardson extrapolation showed that, for $N=100$ control volumes, the relative error was below $2\times {10}^{-4}$. The MATLAB source code used to generate all results is provided as supplementary material.

An objective classification scheme (Appendix B) was used to distinguish sub-critical from super-critical solutions based on the temperature-time history at an internal sensor point. Coupled with a bracketed search in the Damköhler number, this made it possible to determine the critical Damköhler number $D{a}_{\mathrm{crit}}$ for any prescribed set of parameters. These “exact-solution” results are compared with the semi-analytical perturbation predictions, focusing on three aspects: (i) the geometric scaling through ${\beta }_{\mathrm{crit}}$, (ii) the Arrhenius/temperature scaling in the prefactor, and (iii) the predicted peak temperature at criticality.

The remainder of the section is organised as follows:

(i) Characteristics of the transient temperature: Detailed study of the transient behaviour near critical conditions and how it depends on the Arrhenius number;

(ii) Geometrical validation: comparison of the numerically determined $D{a}_{\mathrm{crit}}$ with the ${\beta }_{\mathrm{crit}}(S,P)$ surface and its polynomial fit;

(iii) Chemical/temperature validation: calibration of the kinetic factor $C_{\alpha }$ and assessment of the Arrhenius/processing-temperature scaling;

(iv) Generalisation across kinetics: extension to a wider range of kinetic exponents $(m,n)$ and evaluation of model robustness.

3.1. Characteristics of the Temperature History at Critical Conditions

When the Damköhler number exceeds its critical value, the reaction becomes self-accelerating. However, the qualitative nature of the resulting temperature rise is dependent on the Arrhenius number $A_r$, which dictates whether the runaway is "mild" or "explosive."

• Mild (gradual) runaway - low$A_r$: For Arrhenius numbers below a practical threshold, defined here by

  $$A_{r,\mathrm{crit}}\approx 16.67T_{\infty }(1+T_{\infty }),$$

  (16)

  the peak temperature increases gradually with $Da$, exhibiting a nearly linear dependence on $Da$ up to and beyond the stability limit (see $A_r=30$ in Figure 6). In this regime, crossing $D{a}_{\mathrm{crit}}$ does not trigger a catastrophic jump but rather a manageable increase in the exotherm.
• Ignition-like (explosive) runaway - high $A_r$: For Arrhenius numbers exceeding the threshold in Equation (16), the system behaviour changes character. Sub-critical cases remain stable, but once $D{a}_{\mathrm{crit}}$ is reached the system exhibits a long induction period followed by a rapid temperature rise toward the adiabatic limit (see $A_r=110$ in Figure 7), a behaviour closely resembling ignition in combustible mixtures [12].

At the critical limit ($Da/D{a}_{\mathrm{crit}}=1$), the numerical peak temperature rise $\Delta T_{\max }$ follows the scaling $\Delta T_{\max }\propto T_{\infty }^2/A_r$ predicted by the perturbation analysis. For the cases illustrated in Figure 6 ($A_r=30,40,60,100$), the numerical peak values at critical condition are approximately $0.12,0.08,0.06$ and $0.05$, respectively. These values are slightly lower than the absolute peaks predicted by the perturbation analysis, but the ratios of the magnitudes are consistent with the theory.

The higher peak temperatures in the perturbation model compared to the numerical simulations stem from the simplified kinetics used in the perturbation model, specifically the replacement of the kinetic function $f\left(\alpha \right)$ by a constant, which overestimates the reaction rate throughout the cure. By neglecting reactant consumption and the specific $(m,n)$-order kinetics, perturbation theory provides a conservative upper bound for the temperature rise. In contrast, the numerical model captures the reduction in reaction rate as conversion increases, leading to the slightly lower peak temperatures observed here.

The threshold in Eq. (16) represents a Semenov-type criterion [12] adapted to the current shell geometry and boundary conditions. Its approximate numerical value was determined by identifying the lowest $A_r$ at which the $\Delta T_{\max }$ vs. $Da$ response transitions from a continuous curve to a nearly discontinuous jump when exceeding $D{a}_{\mathrm{crit}}$. Although classical Semenov theory is derived for lumped systems, this empirically calibrated form successfully captures the observed bifurcation in behaviour, explaining why minor process variations in $Da$ can be benign at low $A_r$ but catastrophic at high $A_r$.

3.2. Validation of the Curvature-Chemistry Separation

The validation of the separation of curvature effects from the chemistry $D{a}_{\mathrm{crit}}\propto {\beta }_{\mathrm{crit}}({\sigma }_1,{\sigma }_2,B{i}_d,B{i}_u)$ was done with a fixed chemistry, $(A_r=32.4,T_{\infty }=2.335,m=0.313,n=1.66)$, taken from a reported ${120}^{\circ }\mathrm{C}$ epoxy system [4]. Anchoring the study to a documented material avoids unphysical parameter combinations. This choice does not bias the geometric results because ${\beta }_{\mathrm{crit}}$ is independent of chemistry; chemical effects enter only when converting to $D{a}_{\mathrm{crit}}$.

In what follows, the geometry-Biot dependence of the stability limit ${\beta }_{\mathrm{crit}}$ is taken directly from the polynomial fits in Eqs. (12) and (14) from Section 2. The resulting global stability limit $D{a}_{\mathrm{crit}}$ defined by Equation (10) was fitted to the numerical results for the different heating options using the kinetic scaling factor $C_{\alpha }$ as the only fitting parameter.

We first consider symmetric heating at a high Biot number ($B{i}_d=B{i}_u=1000$). Figure 8 shows the numerically computed stability surface $D{a}_{\mathrm{crit}}({\sigma }_1,{\sigma }_2)$. The surface exhibits the same symmetries as the ${\beta }_{\mathrm{crit}}({\sigma }_1,{\sigma }_2)$ surface from the perturbation analysis, namely invariance under exchange $({\sigma }_1,{\sigma }_2)\mapsto ({\sigma }_2,{\sigma }_1)$ and, under symmetric boundary conditions, invariance under the combined transformation $({\sigma }_1,{\sigma }_2,x)\mapsto (-{\sigma }_1,-{\sigma }_2,-x)$. The one-parameter fit yields $C_{\alpha }=0.213$, about 50% of $maxf\left(\alpha \right)$ for the chosen kinetics, indicating that the conservative choice $C_{\alpha }=max\left(f\left(\alpha \right)\right)$ in the Theory section is unnecessarily strict. The RMS relative error for the fit is 1.8% with a maximum deviation of 4.6% over the whole $({\sigma }_1,{\sigma }_2)$-plane.

For the moderate-convection case $B{i}_d=B{i}_u=5$, the same one-parameter fitting procedure yields an almost identical coefficient, $C_{\alpha }=0.208$, with an RMS relative error of 2.4% and a maximum deviation of 6.4%. The closeness of $C_{\alpha }$ to the high-Bi fit ($0.213$) confirms that the chemistry/process scaling is effectively independent of surface heat-transfer strength, while the Biot number acts primarily through the geometry–boundary factor ${\beta }_{\mathrm{crit}}$ (i.e., a uniform vertical shift of the surface without altering its shape). In other words, even away from the near-Dirichlet limit, the predicted separation of $\mathrm{geometry}+\mathrm{boundaries}$ from $\mathrm{chemistry}+\mathrm{temperature}$ holds: the semi-analytical ${\beta }_{\mathrm{crit}}(S,P)$ surface from the perturbation solution can be used as is, and a single $C_{\alpha }$ brings the model into close agreement with the full numerical results.

For asymmetric heating, the outer surface is taken nearly isothermal $(B{i}_u=1000)$ while the inner surface has finite thermal resistance $(B{i}_d\ge 0)$. The computed stability surfaces again exhibit the symmetries predicted by the perturbation analysis, including the expected tilt along the mean-curvature direction $S={\sigma }_1+{\sigma }_2$ for small $B{i}_d$. The semi-analytical surface ${\beta }_{\mathrm{crit}}(S,P;B{i}_d,B{i}_u)$ from the perturbation analysis is kept unchanged and only the single kinetic factor $C_{\alpha }$ is fitted for each $B{i}_d$. The fitted values lie in a narrow band $0.21\le C_{\alpha }\le 0.22$ over $(B{i}_u=1000,B{i}_d\le 100)$, while the errors decrease rapidly as $B{i}_d$ increases (Table 2). The nominally “worst” case is the perfectly insulated side $(B{i}_d=0)$, where the RMS relative error is 11.6% and the maximum deviation 55.9%. However, this seemingly large error is misleading as a diagnostic: for $B{i}_d=0$ the span of $D{a}_{\mathrm{crit}}$ over the curvature grid exceeds one decade, so uniform absolute errors translate into larger percent errors at the low end. Consistent with this, the parity plot for $B{i}_d=0$ is among the tightest, points cluster closely around the $y=x$ line across the full range (see Figure 9). Overall, the results confirm that the geometry–boundary contribution can be taken directly from the perturbation BVP, while the chemistry–temperature influence is captured by a single scalar $C_{\alpha }$ (here $\approx 0.215$) without re-tuning the geometric fit.

Across all symmetric and asymmetric cases, the best-fit $C_{\alpha }$ remains nearly constant ($\approx 0.21-0.22$) over the entire $({\sigma }_1,{\sigma }_2)$ grid and Biot-number range, reinforcing that geometry/boundary effects are captured by ${\beta }_{\mathrm{crit}}$ while chemistry is well represented by a single scalar $C_{\alpha }$.

3.3. Chemical/Temperature Calibration

Having confirmed that the surface curvature and boundary heating dependence of the stability limit is fully captured by the perturbation theory parameter ${\beta }_{\mathrm{crit}}(S,P,Bi)$, we now turn to the calibration of the chemical prefactor $C_{\alpha }$. This factor represents the effective exothermic strength of the reaction and appears multiplicatively in the closed form criterion for the critical Damköhler number. The objective of this section is to determine whether a single constant value of $C_{\alpha }$ is sufficient to represent the chemical contribution across a broad range of Arrhenius numbers and surface temperatures.

3.3.1. Calibration over $A_r$ and $T_{\infty }$

To isolate chemical effects from geometry and boundary conditions, all calibration simulations were performed for a fixed configuration: vanishing curvature $({\sigma }_1={\sigma }_2=0)$, symmetric heating with $Bi=1000$, and fixed kinetic exponents $(m=0.313,n=1.66)$, identical to those used in the geometrical validation. The validated separation between geometry–heating intensity and chemistry–temperature ensures that this restriction does not bias the calibration.

The calibration spanned a rectangular domain in the $(A_r,T_{\infty })$ plane, $25\le A_r\le 120,2\le T_{\infty }\le 4$, subject to two constraints: (i) $\epsilon \le 0.3$, ensuring the perturbation solution remains accurate, and (ii) exclusion of the region $A_r>16.67T_{\infty }(1+T_{\infty })$, where the system is prone to a sudden, catastrophic jump in temperature and severe overheating for even small excursions above $D{a}_{\mathrm{crit}}$. The resulting admissible domain is shown in Figure 10.

Within the calibration domain, the critical Damköhler number $D{a}_{\mathrm{crit}}$ was extracted from fully coupled transient simulations using the objective classification scheme of Appendix B. These numerical “truth” values were then compared with the perturbation-theory prediction. A least-squares fit over the entire $(A_r,T_{\infty })$ region gives $C_{\alpha }=0.251$, slightly larger than the value obtained in the geometrical study ($\approx 0.213$). The fitted constant represents approximately 59% of $max\left(f\right(\alpha \left)\right)$ for these kinetics and yields an average relative error of 20.5% in $D{a}_{\mathrm{crit}}$. Figure 11 shows that, despite the moderate absolute error, the shape of the isolines and the combined Arrhenius–temperature scaling are reproduced extremely well. A similar fitting procedure was applied to the predicted peak temperature at criticality, yielding an optimal temperature-scaling factor $C_T=1.313$, with an average relative error of approximately 9.5%. As illustrated in Figure 12, the critical-temperature isolines align closely with those of constant $\epsilon$ in Figure 10, consistent with the theoretical structure of the perturbation solution.

Overall, the comparison demonstrates that the semi-analytical perturbation framework captures the dominant chemical and thermal scaling across a broad parameter space. The fitted constants $C_{\alpha }$ and $C_T$ are effectively uniform over the calibration domain, and the remaining discrepancies are of the magnitude anticipated for a first-order perturbation approximation. Most importantly, the isoline alignment confirms that the perturbation solution retains the correct topology of the stability boundary, validating its use as a reduced-order predictive model for thermal runaway.

3.4. Generalisation Across Kinetics

Autocatalytic reaction models of the form $f\left(\alpha \right)={\alpha }^m{(1-\alpha )}^n$ are widely used to describe thermoset curing [10]. However, only a restricted subset of exponent pairs $(m,n)$ produces conversion-rate profiles that resemble experimentally observed behaviour. In particular, very small values of m or unusually large values of n cause the reaction rate to become highly skewed: the peak of $f\left(\alpha \right)$ shifts unrealistically close to either $\alpha =0$ or $\alpha =1$, the maximum value $f_{max}$ becomes excessively sharp, and the overall curve loses the broad, single-hump shape characteristic of epoxy and polyester systems. Such behaviours are inconsistent with DSC measurements reported in the literature and lead to runaway thresholds that do not represent industrial thermosets.

To ensure physically meaningful kinetics while still covering a broad range of autocatalytic behaviours, we therefore restrict attention to

$$0.3\le m\le 0.5,1.0\le n\le 1.7,$$

a domain consistent with reported kinetic fits for common resin families and free from pathological rate-curve shapes. Within this region, $f\left(\alpha \right)$ remains single-peaked, moderately skewed, and shows the characteristic early-stage acceleration followed by depletion-driven decay. Most commercial resin systems fall well within this window, and the associated rate curves differ primarily in width and skewness rather than in qualitative shape. Should values outside this region be required, a new calibration can easily be performed using the software provided in the supplementary material.

Dense Sweep and Global Response-Surface Fit

A full systematic sweep across the $(m,n)$ domain was carried out to obtain a smooth functional form for $C_{\alpha }$. Each $(m,n)$ pair was calibrated against a dense grid in $(A_r,T_{\infty })$, yielding a corresponding value of $C_{\alpha }$ that was used to define a damping factor that relates $C_{\alpha }$ to the maximum of the kinetic function:

$$K(m,n)={\displaystyle \frac{C_{\alpha }}{f_{max}(m,n)}},$$

where $f_{max}={max}_{\alpha }f\left(\alpha \right)$. The sweep used a $5\times 5$ grid in $(m,n)$ and, for each point, a $21\times 10$ grid in $(T_{\infty },A_r)$, with the fitting procedure for $C_{\alpha }$ identical to that used for the chemical calibration.

The resulting optimised parameters $C_{\alpha }$ and $C_T$, their associated relative errors for $D{a}_{\mathrm{crit}}$ and $\Delta T_{\max }$, and the normalised damping factor

$$K={\displaystyle \frac{C_{\alpha }}{f_{max}}}$$

are summarised in Table 3. The consistency of $C_{\alpha }$ and $C_T$ across these widely separated kinetic exponents confirms that the calibration procedure is consistent across the parameter range.

The resulting $K(m,n)$ data were fitted to a low-dimensional polynomial response surface:

$$K(m,n)=p_{00}+p_{10}m+p_{01}n+p_{20}{m}^2+p_{11}mn,$$

(17)

with coefficients

$$\begin{array}{cccccc}\hfill p_{00}& =0.8266,\hfill & \hfill p_{10}& =-0.6373,\hfill & \hfill p_{01}& =-0.1630,\hfill \\ \hfill p_{20}& =2.5335,\hfill & \hfill p_{11}& =-0.0204,\hfill \end{array}$$

yielding ${\mathrm{R}}^2=0.9993$.

A similar fit was obtained for the temperature-scaling factor:

$$C_T(m,n)=p_{00}+p_{10}m+p_{01}n+p_{20}{m}^2+p_{11}mn+p_{02}{n}^2,$$

(18)

with

$$\begin{array}{cccccc}\hfill p_{00}& =1.8510,\hfill & \hfill p_{10}& =-1.2733,\hfill & \hfill p_{01}& =-0.0315,\hfill \\ \hfill p_{20}& =-2.9990,\hfill & \hfill p_{11}& =0.9710,\hfill & \hfill p_{02}& =-0.1082,\hfill \end{array}$$

and ${\mathrm{R}}^2=1.0000$.

Final Predictive Forms

For operational use, the effective rate constant is taken as

$$C_{\alpha }=K(m,n)f\left({\alpha }_{max}\right),$$

(19)

with

$${\alpha }_{max}={\displaystyle \frac{m}{m+n}},f\left({\alpha }_{max}\right)={\left({\displaystyle \frac{m}{m+n}}\right)}^m{\left({\displaystyle \frac{n}{m+n}}\right)}^n.$$

(20)

Within the calibration region

$$0.3\le m\le 0.5,1.0\le n\le 1.7,25\le A_r\le 120,2\le T_{\infty }\le 4,{\sigma }_{1,2}\in [-0.8,0.8],$$

subject to the constraint $\epsilon =T_{\infty }^2/A_r\le 0.3$ and excluding the folded S-curve regime, the damping factor remains bounded between $0.57$ and $0.97$. Outside this region, K should be taken as capped at $K=1.0$, corresponding to the limiting case $C_{\alpha }=f_{max}$, which is guaranteed to yield a critical limit on the safe side of the value obtained from the exact solution.

Substituting (19) into the analytical stability criterion (Eq. (10)) yields:

$$D{a}_{\mathrm{crit}}({\sigma }_1,{\sigma }_2,B{i}_d,B{i}_u,m,n,T_{\infty },A_r)={\displaystyle \frac{{\beta }_{\mathrm{crit}}({\sigma }_1,{\sigma }_2,B{i}_d,B{i}_u)}{K(m,n)f\left({\alpha }_{max}\right)}}{\displaystyle \frac{T_{\infty }^2}{A_r}}exp\left({\displaystyle \frac{A_r}{T_{\infty }}}\right),$$

(21)

where $K(m,n)$ and $f\left({\alpha }_{\max }\right)$ are defined in Eqs. 17 and 20, and ${\beta }_{\mathrm{crit}}({\sigma }_1,{\sigma }_2,B{i}_d,B{i}_u)$ is defined in Eqs. 12 and 14 for symmetric and asymmetric heating respectively.

Temperature-Rise Generalisation

An analogous generalisation applies to the maximum temperature rise $\Delta T_{max}$ at $Da=D{a}_{\mathrm{crit}}$. As shown in Table 3, the factor $C_T$ varies systematically with $(m,n)$, ranging from approximately $0.81$ to $1.36$. These variations are consistent with perturbation-theory expectations: changes in the skewness of $f\left(\alpha \right)$ alter the heat-release rate, but $C_T$ remains $\mathcal{O}\left(1\right)$ across all practical cases. The relative error in $\Delta T_{max}$ remained between $4.6\%$ and $14\%$ across all tested exponent pairs.

Fitting the numerical data to $C_T(m,n)$ provides a smooth functional dependence that reproduces numerical calibrations with high fidelity. Although the model could theoretically be refined by incorporating explicit semi-analytical corrections for local curvatures $({\sigma }_1,{\sigma }_2)$ and Biot numbers $(B{i}_d,B{i}_u)$, the current formulation already demonstrates excellent accuracy. Therefore, this simpler approach was retained to prioritize ease of practical implementation without compromising predictive reliability. The resulting expression for the peak temperature rise is

$$\Delta T_{max}({\sigma }_1,{\sigma }_2,B{i}_d,B{i}_u,m,n,T_{\infty },A_r)=C_T(m,n){\displaystyle \frac{T_{\infty }^2}{A_r}},$$

(22)

using the fit for $C_T(m,n)$ from Equation 18.

3.4.1. Verification of Global Accuracy

To rigorously verify the global accuracy of the fully assembled model, a final validation series was conducted using a fractional factorial experimental design [13]. The parameter space was bounded by the limits of extreme industrial part geometries $({\sigma }_{1,2}=\pm 0.4)$, spanning flat, strongly synclastic, and strongly anticlastic shells. Boundary conditions were restricted to either fully symmetric cooling $(B{i}_d=B{i}_u=1000)$ or heavily asymmetric cooling $(B{i}_d=10,B{i}_u=1000)$. The kinetic exponents were tested at the boundaries of the admissible domain $(m\in [0.3,0.5],n\in [1.0,1.7\left]\right)$.

An 8-run orthogonal test matrix was constructed to sample the corners of this multidimensional space without redundant geometric configurations. For each of the 8 cases, a direct numerical calibration sweep over the $(A_r,T_{\infty })$ plane was performed to extract the true values of $C_{\alpha }$ and $C_T$. These were then compared against the closed-form predictions of the global model. Because $D{a}_{\mathrm{crit}}\propto 1/C_{\alpha }$ and $\Delta T_{\max }\propto C_T$, the relative errors in these prefactors directly represent the model’s predictive error for the runaway threshold and exotherm severity. The results, summarized in Table 4, demonstrate that the global equations successfully capture the compounding effects of asymmetric curvature, heat transfer, and complex kinetics.

3.5. Assessment of In-Plane Thermal Fluxes and Model Assumptions

The theoretical framework presented in Section 2 assumes that through-thickness (transverse) heat conduction is the dominant mode, allowing in-plane conduction along the shell surface to be neglected. To validate this assumption, a quasi-three-dimensional (Quasi-3D) numerical analysis was performed on a representative saddle geometry. The computational domain is discretized via a 3D finite volume mesh, constructed by extruding a body-fitted orthogonal mid-plane mesh toward the inner and outer surfaces. Temperature profiles for the resulting through-thickness "thermal pillars" are initially computed using the 1D solver established in this study. Finally, the potential in-plane fluxes are integrated over the lateral faces of each thermal pillar and compared against the total transverse flux through the pillar’s end faces. This ratio provides a quantitative measure of the 1D assumption’s validity. Detailed algorithmic implementation and the discretization scheme are provided in Appendix B, with the corresponding MATLAB source code available as supplementary material.

A critical observation in the design of shell structures is that geometries that appear highly curved may still correspond to small dimensionless curvatures $({\sigma }_i={\kappa }_{i}w)$ because of the small thickness-to-radius ratio typical of composite laminates. While the perturbation theory in Section 2 remains mathematically valid for dimensionless curvatures up to $\left|\sigma \right|\approx 0.8$, most practical engineering applications, including the saddle geometry in Figure 13 for which $max\left|\sigma \right|\approx 0.1$, fall well within this range.

The results of the Quasi-3D analysis, summarized in Figure 13, demonstrate that even at the onset of thermal runaway, the "worst-case" scenario for thermal gradients, the net in-plane heat flux remains orders of magnitude smaller than the transverse flux. For the saddle geometry, the maximum flux ratio was found to be $0.0058$, corresponding to an in-plane flux of about $0.6\%$ of the transverse flux. This negligible ratio confirms that the high aspect ratio of the laminate and the magnitude of the transverse cooling effectively decouple the thermal behavior of adjacent surface elements. Consequently, the locally 1-D universal model can be used to predict safe process parameters for complex shell geometries with high confidence, as the error introduced by neglecting in-plane conduction is significantly smaller than typical uncertainties in material properties and kinetic parameters.

To further assess the robustness of this approximation, a sensitivity analysis was performed by varying the inner-wall cooling conditions while maintaining high outer cooling $(B{i}_u=1000)$. In the extreme case of a perfectly insulated inner boundary $(B{i}_d=0)$, the maximum in-plane flux ratio was $0.0032$. As the inner cooling is increased, the ratio exhibits a non-monotonic trend: it increases to a peak of $0.022$ at $B{i}_d=1$, then decreases and stabilizes at $0.0056$ for $B{i}_d\ge 100$. This convergence toward the symmetric cooling value $\left(0.0058\right)$ is consistent with the finding in Section 2 that thermal asymmetry becomes negligible when $Bi$ exceeds 100. Because the omission of in-plane conduction effectively removes a potential cooling pathway, the 1-D model serves as a conservative "worst-case" predictor; it will exhibit thermal runaway slightly earlier than a full 3-D system, therefore providing an inherent safety margin for industrial process design.

4. Discussion

The Power of Dimensional Reduction and Universal Scaling. Although modern computational power enables direct numerical simulation (DNS) of thermal runaway in specific geometries, such point solutions provide limited insight into the governing physics. By reducing the problem dimensionally, the high-dimensional space of physical variables (thermal properties, kinetics, geometry, and process conditions) collapses into a small set of governing parameters $(Da,A_r,Bi,\sigma )$. This reveals scaling laws that are otherwise obscured. In particular, stability does not depend on curvature or thickness alone, but on their dimensionless product $\sigma =\kappa w$. The resulting framework defines a universal “safe-by-design” envelope through $D{a}_{\mathrm{crit}}$, independent of discretisation details or solver settings.

Validity of the Locally 1-D Approximation. A central question is whether neglecting in-plane conduction introduces a significant error. Quasi-3D analysis (Section 3.5) shows that for dimensionless curvatures and aspect ratios typical of composite laminates, transverse gradients dominate the thermal response. Even under strong curvature and asymmetric boundary conditions (e.g., $B{i}_d=0$), the in-plane flux remains on the order of a few percent of the total heat dissipation. This confirms that the locally 1-D formulation captures the dominant heat-transfer pathways and is not merely a mathematical simplification.

Generality across arbitrarily shaped geometries. The applicability to arbitrary shell geometries follows from the fact that any smooth surface is locally defined by its dimensionless principal curvatures, ${\sigma }_1$ and ${\sigma }_2$. By expressing the problem in terms of these local invariants, the model decouples the stability criterion from global geometry. Therefore, complex structures can be analysed as collections of independent “thermal pillars”, each governed by local curvature. In practice, even visually pronounced curvature often corresponds to moderate dimensionless values $(\parallel \sigma \parallel <0.5)$, where the analytical model is most accurate.

Robustness under Asymmetric Boundary Conditions. Perfectly symmetric cooling is rare in practice, which makes asymmetric boundary conditions essential to consider. The results show that stability is strongly influenced by the better-cooled surface: even modest increases in the corresponding Biot number significantly raise the runaway threshold. This insight has direct practical implications. For example, improved convection on the accessible surface can compensate for limited tool-side heat extraction, an effect that can be quantified directly using the analytical model without resorting to full 3D simulations.

Practical Implementation for Process Design. The framework is readily applicable in engineering practice. Local principal curvatures can be extracted from standard CAD tools, after which the semi-analytical expressions are evaluated pointwise over the surface. This identifies the least stable location through the minimum $D{a}_{\mathrm{crit}}$, automatically capturing the reduced stability of anticlastic regions relative to synclastic ones. Combined with estimates of Biot numbers and processing temperature, the model enables rapid assessment of safe operating conditions. If excessive peak temperature is predicted, mitigation can be achieved by increasing surface heat transfer or introducing intermediate dwell steps to reduce the accumulated reaction heat.

Consistency Between Dynamic Detection and Perturbation Theory. The second-derivative criterion aligns closely with the perturbation-based stability limit. In the perturbation framework, runaway occurs when no steady temperature profile can balance heat generation and conduction. Dynamically, this same transition is detected when ${\ddot{T}}_s=0$, marking the onset of self-acceleration. The agreement between these independent approaches demonstrates that the transient acceleration signal reflects the same underlying stability limit as the analytical theory. This provides strong support for using the second derivative as an early-warning indicator in practical applications.

Real-Time Monitoring and Active Process Control. The second-derivative criterion used to classify the numerical solutions also provides a practical basis for real-time process monitoring. During stable curing, heat removal dominates and the temperature response is decelerating $({\ddot{T}}_s<0)$. If ${\ddot{T}}_s$ starts increasing and approaches zero from the negative side, the system reaches a critical balance where heat generation begins to exceed cooling capacity. This transition is detectable before large temperature rises occur. Monitoring the curvature of the temperature signal therefore provides an early-warning indicator and enables timely intervention, for example by reducing the tool temperature or introducing a temporary dwell. In this way, the tooling can act as an active control element rather than a passive thermal boundary.

5. Conclusions

In this work, a universal stability framework for thermal runaway during the cure of composite shells has been developed and validated. The main conclusions are:

• A unified analytical criterion. A general analytical stability criterion has been derived for arbitrary shell geometries and convective (Robin) boundary conditions. The resulting closed-form expression for the critical Damköhler number, $D{a}_{\mathrm{crit}}$, separates the influence of geometry and heat transfer from that of chemistry and processing temperature. This reduces the full 3D stability problem to a single non-dimensional margin $Da/D{a}_{\mathrm{crit}}$, providing a universal “safe-by-design” framework.
• Validation of the locally 1-D approximation. The central modelling assumption, that through-thickness conduction dominates over in-plane transport, was validated using a quasi-3D flux analysis. Even for curved shells and asymmetric boundary conditions, in-plane heat fluxes remain small (typically $\ll 1\%$ of the transverse flux). This confirms that the locally 1-D formulation accurately captures the dominant heat-transfer mechanisms and provides a conservative estimate of stability.
• Geometric scaling of stability. The analysis shows that geometric effects enter the stability problem through the dimensionless curvature $\sigma =\kappa w$, rather than curvature or thickness separately. This scaling enables direct transfer of results across different materials and part sizes, and allows complex shells to be treated as collections of independent local “thermal pillars” governed by principal curvatures.
• Arrhenius-dependent runaway behaviour. The nature of thermal runaway is governed by the Arrhenius number $A_r$. Low $A_r$ leads to gradual, mild transitions, while high $A_r$ produces ignition-like behaviour with long induction periods followed by rapid temperature rise. The analytical framework captures the stability limit in both regimes, providing guidance on whether process deviations are likely to be manageable or catastrophic.
• Role of boundary conditions and tooling. Boundary heat transfer has a dominant influence on stability, particularly under asymmetric conditions. The results show that the better-cooled surface controls the stability margin, implying that improvements in accessible surfaces (e.g. enhanced convection) can compensate for poor cooling elsewhere. This highlights the dual role of tooling as both a heat source and a heat sink.
• Practical applicability and implementation. The semi-analytical model is computationally inexpensive and can be applied directly using curvature data from CAD models and standard heat-transfer estimates. It enables rapid identification of critical locations, prediction of safe thickness, and evaluation of process modifications such as dwell steps or improved cooling strategies.

Together, these results provide a physically transparent and practically applicable framework for predicting and preventing thermal runaway in the manufacture of thick composite structures.

Appendix A.1. Controlled Approximations

To transform the transient coupled system into a tractable stability problem, three approximations are introduced:

1.

Decoupling of kinetics

The autocatalytic factor $f\left(\alpha \right)$ is replaced by its maximum value $C_{\alpha }=f\left({\alpha }_{\max }\right)$, where ${\alpha }_{\max }=m/(m+n)$. This produces a conservative, worst-case estimate of heat generation.

2.

Small temperature-rise expansion

The temperature field is expanded as a perturbation series about the zeroth-order solution $T_0$:

$$T(x,t)=T_0(x,t)+\epsilon \vartheta (x,t)+\mathcal{O}\left({\epsilon }^2\right)$$

(A1)

where $\vartheta$ is the first-order perturbation function. The dimensionless small parameter $\epsilon$ will be selected later to simplify the resulting perturbation equation.

3.

Steady-state limit

Although the perturbation approach allows the analysis of the fully transient problem, the stability of the problem is determined by the behaviour of the solution after the initial heating phase. In the absence of internal heating, the boundary-driven temperature reaches the limit $T_{\infty }$ after approximately one dimensionless time unit since the scaling of time is done with the characteristic heat conduction time. If a balance between heat release and cooling by the boundaries can be maintained at this thermal maximum, it must also be valid during the preceding transient heating. The perturbation problem is therefore formulated at the steady-state limit where $\partial T/\partial t=0$. The resulting perturbation ansatz therefore becomes:

$$T\left(x\right)=T_{\infty }+\epsilon \vartheta \left(x\right)+\mathcal{O}\left({\epsilon }^2\right),$$

where the zeroth order solution $T_0(x,t)$ has been replaced by its constant stationary solution $T_{\infty }$.

Appendix A.2. The First-Order Perturbation Equation

The Arrhenius exponential in the energy equation creates a difficulty, but it can be eliminated by chosing the small parameter $\epsilon =T_{\infty }^2/A_r$, which after Taylor expansion of the argument yields:

$$exp\left(-{\displaystyle \frac{A_r}{T_{\infty }+\epsilon \vartheta }}\right)\approx exp\left(-{\displaystyle \frac{A_r}{T_{\infty }}}\right)exp\left(\vartheta \right),$$

Substituting the perturbation series into the decoupled energy equation and collecting all $\mathcal{O}\left(\epsilon \right)$ terms yields the first-order ODE:

$${\displaystyle \frac{d^2\vartheta }{d{x}^2}}+\left({\displaystyle \frac{{\sigma }_1}{1+{\sigma }_{1}x}}+{\displaystyle \frac{{\sigma }_2}{1+{\sigma }_{2}x}}\right){\displaystyle \frac{d\vartheta }{dx}}+\beta exp\left(\vartheta \right)=0,$$

(A2)

where the dimensionless stability parameter is:

$$\beta =Da{C}_{\alpha }exp\left(-{\displaystyle \frac{A_r}{T_{\infty }}}\right){\displaystyle \frac{1}{\epsilon }}$$

(A3)

Appendix A.3. Boundary Conditions

The dimensionless Robin conditions applied in the main text lead, at first order, to:

$$\begin{array}{cccc}\hfill {\displaystyle \frac{d\vartheta }{d\tilde{x}}}-B{i}_{\mathrm{d}}\vartheta & =0,\hfill & \hfill & \tilde{x}=-1,\hfill \end{array}$$

(A4)

$$\begin{array}{cccc}\hfill {\displaystyle \frac{d\vartheta }{d\tilde{x}}}+B{i}_{\mathrm{u}}\vartheta & =0,\hfill & \hfill & \tilde{x}=+1.\hfill \end{array}$$

(A5)

Equations (A2)–(A5) form a nonlinear boundary-value problem for $\vartheta \left(x\right)$ with $\beta$ as a parameter.

Appendix A.4. Definition of the Critical Stability Limit

A stationary solution exists only if the heat-generation term $\beta e^{\vartheta }$ can be balanced by conduction in the curved shell. The critical stability limit is therefore defined as the value ${\beta }_{\mathrm{crit}}$ where the solution to Equations (A2) - (A5) disappears. Because geometry and boundary conditions appear only in the differential operator and BCs, ${\beta }_{\mathrm{crit}}={\beta }_{\mathrm{crit}}({\sigma }_1,{\sigma }_2,B{i}_d,B{i}_u)$.

Appendix A.5. Peak Temperature Rise at Criticality

The first-order perturbation solution predicts a peak temperature rise at criticality of

$$\Delta T_{\mathrm{crit}}={\vartheta }_{max},\mathrm{crit}\epsilon ={\vartheta }_{max},\mathrm{crit}{\displaystyle \frac{T_{\infty }^2}{A_r}},$$

(A6)

where ${\vartheta }_{max},\mathrm{crit}={max}_x\vartheta \left(x\right){|}_{\beta ={\beta }_{\mathrm{crit}}}$.

Appendix B.1. The Perturbation Boundary-Value Problem

The geometric stability limit, ${\beta }_{\mathrm{crit}}({\sigma }_1,{\sigma }_2,B{i}_{\mathrm{d}},B{i}_{\mathrm{u}})$, is determined by identifying the threshold value of $\beta$ beyond which the nonlinear boundary-value problem (Equation A2) ceases to possess a stationary solution. This is implemented using MATLAB’s bvp4c solver, which utilizes a fourth-order collocation method based on Lobatto quadrature with adaptive mesh refinement.

To ensure robust convergence, a sequential continuation strategy is employed. For every new geometric configuration, the search initializes with a “cold start” at $\beta =0$ using an analytical parabolic profile. The parameter $\beta$ is then incrementally increased ($\Delta \beta ={10}^{-4}$), utilizing the converged solution $\vartheta \left(x\right)$ from the previous step as the initial guess for the next. The critical limit ${\beta }_{\mathrm{crit}}$ is defined as the last value of $\beta$ for which the solver successfully converges, guaranteeing an accuracy of $\le \Delta \beta$.

Appendix B.2. The Fully Coupled Transient Model

The fully coupled 1D transient model is solved using the Method of Lines [14]. The spatial derivatives are discretized using a second-order, cell-centered Finite Volume Method (FVM) [15]. To maintain global second-order accuracy and prevent order reduction at the boundaries, the Robin boundary fluxes are evaluated using temperatures extrapolated to the physical walls ($T_{\mathrm{wall},\mathrm{d}}$ and $T_{\mathrm{wall},\mathrm{u}}$) via a second-order one-sided stencil. The resulting system of stiff ordinary differential equations is integrated in time using MATLAB’s ode15s solver with a relative tolerance of ${10}^{-6}$.

A grid convergence study using Richardson extrapolation [15] on three nested grids ($N=100,200,400$) confirmed that the spatial discretization achieves a super-convergent order of $p\approx 3.0$ to $3.6$ for the peak temperature. The estimated spatial discretization error for the production grid of $N=101$ nodes is verified to be below $2\times {10}^{-4}$.

Appendix B.3. Quasi-3D Validation Framework

To validate the locally 1D assumption, a Quasi-3D analysis was performed by mapping the shell surface onto an orthogonal grid of local principal curvatures. The 1D transient solver is executed independently for each resulting volumetric "thermal pillar” (frustum) to generate a 3D transient temperature field, $T(u,v,x,t)$.

At the moment of maximum exotherm ($t=t_{\mathrm{peak}}$), the lateral heat leakage through the East, West, North, and South faces of every layer i is calculated using a central finite difference between neighboring thermal pillars. For example, the flux through the East face of a thermal pillar at layer i is determined by the temperature gradient between the current pillar and its neighbor in the u-direction. The total lateral flux magnitude for a given thermal pillar is summed over all N layers:

$$F_{\mathrm{lateral},\mathrm{total}}={\displaystyle \sum _{i=1}^N}\left|F_{u,i}^{\mathrm{East}}+F_{u,i}^{\mathrm{West}}+F_{v,i}^{\mathrm{North}}+F_{v,i}^{\mathrm{South}}\right|$$

(A7)

The validity of the 1D approximation is quantified by the ratio of this lateral leakage to the primary transverse heat exchange at the boundaries:

$$\mathrm{Ratio}={\displaystyle \frac{F_{\mathrm{lateral},\mathrm{total}}}{|F_{x,\mathrm{u}}|+|F_{x,\mathrm{d}}|}}$$

(A8)

A ratio $\ll 1$ confirms that transverse conduction strictly dominates the local energy balance. Detailed implementation of the geometric mapping and the face-flux calculations can be found in the computeLateralFlux.m function within the provided supplementary software.

Appendix B.4. Thermal Runaway Detection and Critical Limit Search

The transition to thermal runaway is identified by evaluating the temperature history, $T_{\mathrm{s}}\left(t\right)$, at a fixed virtual sensor located at the spatial node reaching the highest absolute temperature. The stability is classified by analyzing the second derivative with respect to time (${\ddot{T}}_{\mathrm{s}}$).

To isolate the chemical exotherm from the initial boundary heating, the signal is truncated for $t<0.2$. Furthermore, the analysis window is strictly bounded at $t=t_{\mathrm{peak}}$ to prevent the algorithm from locking onto the asymptotic cooling phase. Within this window, a sub-critical cure is characterized by a strictly negative acceleration peak (${\ddot{T}}_{\mathrm{s},\mathrm{peak}}<0$), indicating sufficient cooling. A super-critical runaway is identified by a positive acceleration peak (${\ddot{T}}_{\mathrm{s},\mathrm{peak}}>0$), marking the onset of self-acceleration.

The critical Damköhler number, $D{a}_{\mathrm{crit}}$, is identified as the root of this classification function (${\ddot{T}}_{\mathrm{s},\mathrm{peak}}=0$). This is located using a robust search algorithm that first expands to establish a strict bracket $[D{a}_{\min },D{a}_{\max }]$ spanning the sign change, followed by execution of MATLAB’s fzero routine within this bracket to resolve $D{a}_{\mathrm{crit}}$ to a relative tolerance of ${10}^{-5}$.