A Theory-Guided Machine Learning and Molecular Dynamics Approach for Characterizing Fast-Curing Polyurethane Systems

1. Introduction

Polyurethane (PU) is a thermoset polymer formed by reactions between isocyanates with a segmented molecular structure consisting of soft segments and hard urethane linkages [1]. The chemical structure of polymers, together with their processing methods, ultimately governs the macroscopic properties of the material, thereby determining its suitability for downstream functional applications [2]. Owing to the versatility of their chemistry, PU can be synthesized from a wide range of structural units with varying molecular sizes and chemical functionalities. This flexibility enables the design and fabrication of a diverse family of PU materials with tailored chemical architectures, allowing their mechanical and thermal properties to be systematically tuned to meet specific application requirements [3]. PU is ubiquitous in structural, automotive and aerospace applications, including insulators, foams, and binders [1,4]. Pultrusion is a continuous manufacturing process widely used to produce fiber-reinforced polymer composites [5]. Among thermoset polymers, PU has emerged as a preferred matrix material due to its rapid curing kinetics and superior mechanical performance, making it highly suitable for pultruded applications in aviation, transportation, and construction [6,7]. Despite its technological importance, quantitatively tracking the transition from a low-viscosity liquid to a gel and ultimately solid during curing remains challenging, particularly for the glass transition behavior. Experimental characterization techniques such as differential scanning calorimetry (DSC) and thermomechanical analysis (TMA) are widely used to measure curing heat flow and thermomechanical responses such as thermal expansion, respectively [8,9,10]. However, for fast-curing PU, accurately capturing $T_g$ at low degrees of cure is difficult because rapid kinetics and overlapping thermal signals limit experimental resolution in the early curing regime. On the other hand, the $T_g$ is not governed by a fixed cure degree but rather by the attainment of a critical network density, ν_e, corresponding to the formation of an infinite molecular network [11]. Current industrial practice often assumes the $T_g$ and gelation as a fixed material property or directly correlates it with a single degree of cure (DoC), α, which oversimplifies the physical mechanisms governing network evolution and segment migration during curing [11,12]. This can lead to inaccurate predictions of $T_g$ and gelation onset, resulting in defects such as resin-rich zones, incomplete fiber wet-out, or premature solidification. A more fundamental, physics-based description of polymerization is needed to improve process control and reproducibility. Recent reviews of polyurethane (PU) cure-kinetics characterization emphasize that calorimetry/spectroscopy alone is often insufficient and that rheological monitoring is essential for capturing evolving physical properties (e.g., viscosity, gelation, pot life), particularly once diffusion and heterogeneity become important [13]. Related multi-technique characterization and modeling workflows have also been demonstrated for fast-curing PU adhesives to link conversion to thermomechanical transitions and processing windows [14].

Molecular dynamics (MD) simulation provides a molecular-level prediction of glass transition behavior from an atomistic structure and thermomechanical response. According to Han et al. [15], the $T_g$ of polyethylene, polystyrene, and polyisobutylene were captured by evaluating the thermal expansion of each polymer, and the $T_g$ was defined at the intersection point of lower temperatures and higher temperatures. In addition, the Enhanced Monte Carlo (EMC) package developed by Pieter J. in-’t Veld [16] has been used to generate equilibrated polymer configurations and structural descriptors prior to thermomechanical analysis. In fast-curing PU, these simulation and experimental datasets are often sparse and noisy in the low-DoC regime, motivating data-fusion approaches that enforce physical consistency. Recent advances in machine learning (ML) provide a powerful pathway to address these challenges. Probabilistic surrogate models, particularly Gaussian Process Regression (GPR), are well suited for learning smooth structure–property relationships from sparse and noisy datasets while simultaneously providing quantified predictive uncertainty [17,18]. When coupled with physics-informed priors, such as curing kinetics or thermodynamic constraints, GPR based surrogates can further guide MD simulations and experimental design, enabling more efficient exploration of the parameter space and targeted data acquisition. In this work, MD simulations with topological crosslinking and thermal expansion analysis were used to estimate crosslink density evolution and $T_g$ across degrees of cure [19,20,21,22]. These results were integrated with experimental characterization to construct a conversion-dependent transition model and an uncertainty-aware predictive $T_g$-(α) relationship using theory-guided machine learning. The resulting framework supports process-modeling efforts for fast curing PU systems incorporated in pultrusion. Similar MD-based determinations of $T_g$ and the volumetric coefficient of thermal expansion (CTE) from density/specific-volume–temperature trends have been reported for thermoset polymer nanocomposites [23]. Hybrid all-atom/coarse-grained MD strategies have further been used to connect microphase morphology and hydrogen bonding to macroscopic toughness in related segmented urea-based elastomers (polyureas), illustrating the value of multiscale simulation for urethane/urea chemistries [24].

2. Materials and Methods

2.1. Material

In this study, ELASTOCOAT® 74850R resin and ELASTOCOAT® 74850T isocyanate were supplied by BASF North America. These two components were used as received without further purification and mixed according to the manufacturer-recommended formulation. ELASTOCOAT® 74850R serves as the polyol component, while ELASTOCOAT® 74850T acts as the reactive isocyanate, together forming a fast-curing polyurethane system suitable for curing-kinetics investigations.

2.2.1. Curing Kinetics Evolution

To evaluate curing kinetics model for PU, dynamic DSC measurements were applied with DSC 250 from TA. Polyol and isocyanate were stored in an insulated container with freezer packs to maintain a temperature around 5 °C. Plastic film and zip bags were used to prevent moisture from the environment. Since the PU system is known to react even at subzero conditions (initiating as low as −25 °C), these precautions slowed down the reaction kinetics and ensured that minimal curing occurred prior to DSC measurement. Thus, the recorded exothermic heat accurately reflects the intended curing process. Non-isothermal DSC conducted at multiple heating rates, combined with isoconversional and model-fitting approaches, is a common strategy for identifying apparent activation energies and multi-step mechanisms in fast-curing PU systems [25].

To investigate the dynamic curing kinetics of the PU system, the samples were subjected to a non-isothermal heating program over the temperature range of -50–250 °C at a constant heating rate of 3, 5, 10, and 20 °C/min.

2.2.2. Gelation

A TA Instruments Discovery HR30 rheometer was used to measure the viscosity evolution of the PU system using dynamic oscillatory tests. To correlate viscosity with the DoC, experiments were conducted at 50 °C under various angular frequencies. The evolution of curing behavior was monitored to identify the gelation point. All tests were performed using 20 mm diameter parallel plates with a gap of 300 microns. Rheometry is particularly valuable for PU systems because it directly captures physical-property evolution during cure; for heterogeneous reactions, these physical changes can dominate and are not fully accessible through thermal or spectroscopic techniques alone [13].

2.2.3. Glass Transition Temperatures

To capture glass transition temperature, the following isothermal & dynamic scanning DSC method was used: Partial curing of the PU samples was achieved through isothermal holding at 100 °C for durations of 1, 3, 5, 15, and 30 min, respectively. Cooling down to -60 °C holding 5 mins and ramping to 250 °C with 10 °C/min. Under the same measure groups and testing conditions, the exothermic behavior during the cooling steps was neglected. $T_g$was determined by DSC from the second heating cycle when a distinguishable glass transition was observed. At low degrees of cure (below 0.5), $T_g$could not be reliably resolved due to weak or overlapping thermal transitions.

2.2.4. Molecular Dynamic Simulation capturing T_g

Due to the limited understanding of glass transition behavior below the 0.5 DoC, accurately predicting $T_g$ evolution in the early curing regime remains challenging. MD simulations can be employed to help bridge this gap by providing molecular-level insight into network formation and thermomechanical behavior at low DoC.

A PU MD model was set up in LAMMPS [26] with a 50.00 x 50.00 x 50.00 ų cubic simulation box and a density of 1120 kg/m³. The optimized potentials for liquid simulation all atom (OPLS-AA) force field was used to define the interaction between atoms for the crosslinking simulation. Time integration was performed using the velocity Verlet algorithm with a time step of 0.25 fs. Long-range electrostatic interactions were computed using the particle–particle particle–mesh (PPPM) method, with a cutoff distance of 1.2 nm [22]. In addition, “Condensed-phase Optimized Molecular Potential for Atomistic Simulation Studies” (COMPASS) and “Polymer consistent force field” (PCFF) were applied to investigate the volume expansion capturing $T_g$ [27,28].

A model of 400 monomers was used in the simulation box, corresponding to a polyol-to-isocyanate weight ratio of 46:54 as provided by BASF. The number of monomers was determined based on their molecular weights, and the resulting simulation configuration is shown in Figure 1. To enable crosslinking between the target reactive functional groups and to capture the $T_g$ of the polymer system, four sequential simulation steps were implemented in LAMMPS, such as constant Number of particles, Pressure, and Temperature (NPT) densification, crosslinking, creating chemical bonds, and T_g scanning. During NPT densification, the environment was equilibrated at 300 K under 1atm pressure. After densified the simulation box, the crosslinking in this work was introduced using a topology algorithm without relying on chemical reactions. The hydroxyl oxygen atoms in the polyol and the double-bonded carbon atoms in the isocyanate were identified as active atoms. The crosslinks were formed when the distance between active atoms equal to or less than the cutoff radius of 3.0 Å [22] and the schematic of this behavior is presented by Figure 1; The chemical bonds between atoms were introduced based on the parameter data file generated from the EMC package; The $T_g$scanning was performed after the chemical bonds mapping. First, the crosslinked PU system was loaded along with force field parameters, and the initial structure undergoes energy minimization to eliminate unreasonable contacts. Subsequently, isothermal stabilization is performed under NPT conditions at high temperature to fully homogenize the system and erase initial configuration history. Based on this, the external pressure is held constant while the system temperature is linearly decreased from $T_{high}$ to $T_{low}$ under NPT conditions to obtain the density–temperature response during cooling. A brief isothermal NPT stabilization at the low-temperature end minimizes nonequilibrium effects. Finally, under identical pressure conditions, the temperature is linearly raised from $T_{low}$back to T_high to derive the density vs temperature curve during heating. The $T_g$was determined by segmented linear fitting of the density vs temperature curves from the cooling and heating phases, with the intersection point of the two segments’ slope changes serving as the $T_g$ value.

2.3. Physics-Informed Prediction

To predict the glass transition behavior of the PU system based on experimental DSC data and MD simulation results, GPR was employed as a bridging framework to fit the target T_g curve [17,18]. In this approach, T_g (α) data obtained from experiments and MD simulations were used as training inputs, guided by the physically motivated DiBenedetto equation as a prior constraint.

The DiBenedetto relationship is given by Eq. 1[29]:

$$T_g\left(\alpha \right)=T_{g0}+\left(T_{g\infty }-T_{g0}\right){\displaystyle \frac{\lambda \alpha }{1-\left(1-\lambda \right)\alpha }}$$

(1)

The outputs of this model include the DiBenedetto parameter (λ), a continuous T_g – α relationship, and the associated data uncertainty. To enable smooth and continuous material property prediction, a Matérn kernel was selected [30]. A constant term was introduced to provide additional flexibility, and the kernel was combined with a noise term to account for random noise arising from experimental measurements and MD fluctuations [31]. The Matérn kernel function can be expressed as in Eq. 2 and the final kernel function presented by Eq. 3:

$$k\left(x_i,x_j\right)={\displaystyle \frac{1}{\mathsf{\Gamma }\left(v\right)2^{v-1}}}\left({\displaystyle \frac{\sqrt{2v}}{l}}d\left(x_i,x_j\right)\right)X{K}_v\left({\displaystyle \frac{\sqrt{2v}}{l}}d\left(x_i,x_j\right)\right)$$

(2)

where d(x_i, x_j) is the Euclidean distance, K_v is a modified Bessel function, ε is noise and Γ is the gamma function; $v$ is the smoothness parameter controlling the differentiability of the Matérn kernel; $l$ is length scale parameter.

$$Kernel=C·k\left(x_i,x_j\right)+\epsilon \left({\sigma }_n^2{\delta }_{ij}\right)$$

(3)

where C is a constant; The term $\epsilon \left({\sigma }_n^2{\delta }_{ij}\right)$represents the observation noise variance added to the diagonal of the covariance matrix; ${\delta }_{ij}$is the Kronecker delta, equal to 1 when $i=j$ and 0 otherwise.

As this study relies on hybrid data sources with a limited number of data points, cross validation is critical to ensure robust and reliable predictions. Therefore, k-fold cross validation [32] was applied to evaluate model performance and mitigate overfitting. Within this framework, the relationship between the observed data and the latent target function can be expressed as:

$$Target:T_g\left(\alpha \right)=T_{g0}^{DiB}\left(\alpha ;\lambda \right)+\delta \left(\alpha \right)$$

(4)

$$T_{g,i}^{obs}\left({\alpha }_i\right)=T_g\left({\alpha }_i\right)+{\epsilon }_i{\epsilon }_i~N(0,{\sigma }_i^2)$$

(5)

$${\sigma }_i^2=\left\{\begin{array}{c}{\sigma }_{Exp}^2,\\ {\sigma }_{MD}^2,\end{array}\begin{array}{c}i\in Experiment\\ i\in MD\end{array}\right.$$

(6)

where ${\sigma }_{\mathrm{E}\mathrm{x}\mathrm{p}}^2$and ${\sigma }_{\mathrm{M}\mathrm{D}}^2$represent the noise variances associated with experimental and MD data points, respectively; $\delta \left(\alpha \right)$ is residual correction .

A constant experimental noise level σ_Exp = 0.4 was assumed, while the MD noise level σ_MD was automatically selected by minimizing a predefined objective function. The fitting performance was evaluated using the root mean square error (RMSE), defined as

$$RMS{E}_{CV}\left(\lambda ,{\sigma }_{MD}\right)={\left[{\displaystyle \frac{1}{N}}\sum _{k=1}^K\sum _{i\in D_{val}^{\left(k\right)}}{\left(T_{g,i}^{obs}-T_g^{\left(k\right)}\left({\alpha }_i\right)\right)}^2\right]}^{{\displaystyle \frac{1}{2}}}$$

(7)

where D_val^(k) denotes the validation subset in the k-th fold during the k-fold cross validation procedure; $K$is the number of cross-validation folds; N is the total number of validation data points across all folds, $T_{g,i}^{\mathrm{o}\mathrm{b}\mathrm{s}}$ is the observed glass transition temperature of the $i$-th data point, and $T_g^{\left(k\right)}\left({\alpha }_i\right)$ is the predicted glass transition temperature evaluated using the model trained in the $k$-th fold at degree of cure ${\alpha }_i$.

In addition, a curvature regularization term $J_{curv}$was introduced to discourage spurious oscillations in the predicted $T_g$–DoC relationship. Since the true evolution of $T_g$with curing was expected to be smooth and physically monotonic, model selection was performed by minimizing the cross-validated RMSE together with the curvature regularization term and the target function is described in Eq. 8 [33,34,35].

$$S\left(\lambda ,{\sigma }_{MD}\right)=\arg \min \left[RMS{E}_{CV}\left(\lambda ,{\sigma }_{MD}\right)+J_{curv}\right]$$

(8)

3. Results

3.1. Curing Kinetics Evolution

To capture the first glass transition temperature, $T_{g,0}$, a PU sample was tested from -80 to 250 °C and back at 5 °C /min. The initial glass transition temperature, $T_{g,0}$ was -56.32 °C and the final glass transition temperature, $T_{g,\mathrm{\infty }}$ was 89.12 °C. Figure 2 presents the normalized heat-flow curves of the PU system. As shown in Figure 2 (b), the dynamic curing processes under three different heating rates were recorded and analyzed using DSC. The corresponding curing characteristic parameters are summarized in Table 1, including the onset curing temperature (T_onset), peak temperature ($T_p$), endset curing temperature (T_endset), and the overall curing temperature range of the PU system. The total reaction enthalpy ($∆H_{tot}$) was determined by integrating the area under the exothermic peak in Figure 2 (a) according to Eq. 9 and the average value was around 138.10 J/g.

$$∆H_{tot}={\int }_{onset}^{endset}\dot{q}\left(t\right)dt$$

(9)

where $\Delta H_{\mathrm{tot}}$ is the total reaction enthalpy, $\dot{q}\left(t\right)$is the DSC heat flow rate as a function of time, and $T_{\mathrm{onset}}$ and $T_{\mathrm{end}}$denote the onset and endset temperatures of the curing reaction, respectively.

To evaluate the conversion of the PU system, degree of cure, $\alpha$ and curing rate, ${\displaystyle \frac{d\alpha }{dt}}$ were integrated based on Eq. 10 to 12 [36]. The conversion and curing-rate curves were calculated and are plotted in Figure 3 (a) and Figure 3 (b), respectively, for different heating rates.

$$\alpha ={\displaystyle \frac{∆H}{∆H_{tot}}}$$

(10)

$$∆H={\int }_{t_0}^t\dot{q}\left(t\right)dt$$

(11)

$${\displaystyle \frac{d\alpha }{dt}}={\displaystyle \frac{\dot{q}\left(t\right)}{∆H_{tot}}}$$

(12)

The curing rate can be evaluated in kinetic analysis and the kinetic process can be described as [37]

$${\displaystyle \frac{d\alpha }{dt}}={\displaystyle \frac{\dot{q}\left(t\right)}{∆H_{tot}}}=K\left(T\right)f\left(\alpha \right)$$

(13)

where $\dot{q}$ is heat flow during DSC test; K(T) as a temperature-dependent reaction rate constant is described by Eq. 14, f($\alpha$) a dependent kinetic model function which can be calculated using the logarithmic form of the kinetics below

$$K\left(T\right)=Aexp(-{\displaystyle \frac{E_a}{RT}})$$

(14)

$$ln{\displaystyle \frac{d\alpha }{dt}}=\mathit{ln}\left[Af\left(\alpha \right)\right]-{\displaystyle \frac{E_a}{RT}}$$

(15)

where $A$ is the pre-exponential factor; $R$ is the universal gas constant.

The activation energy ($E_a$), a key parameter in the curing kinetics model, was determined from the slope of the corresponding linearized isoconversional equation. To improve numerical stability and accuracy, the experimental data were smoothed and normalized by introducing the auxiliary functions $y\left(\alpha \right)$and $z\left(\alpha \right)$ [38,39]. Furthermore, the temperature integral is approximated using the fourth-order rational expression, $\pi \left(x\right)$ proposed by Senum and Yang [40].

$${y\left(\alpha \right)=\left({\displaystyle \frac{d\alpha }{dt}}\right)e}^x$$

(16)

$$z\left(\alpha \right)=\pi \left(x\right)\left({\displaystyle \frac{d\alpha }{dt}}\right){\displaystyle \frac{T}{\beta }}$$

(17)

where $\beta$ is DSC heating rate.

The pre-exponential factor, A can be determined by Eq. 18 [37,41]

$$A=\left({\displaystyle \frac{\beta x_p}{Tf\mathrm{\text{'}}\left({\alpha }_p\right)}}\right)exp{x}_p$$

(18)

where $f^{\mathrm{\text{'}}}\left({\alpha }_p\right)$ is the derivative of the kinetic model with respect to the degree of cure, ${\alpha }_p$ is the degree of cure corresponding to the maximum reaction rate on the DSC curve, and the subscript $p$ denotes quantities evaluated at the peak of the DSC exotherm; $x_p$ denotes the reduced activation energy, defined as $x_p=E_a/\left(R{T}_p\right)$, with $T_p$corresponding to the peak temperature of the DSC curve.

Figure 4. Variations of (a) y(α) and (b) z(α) functions vs conversion for PU system.

Figure 4. Variations of (a) y(α) and (b) z(α) functions vs conversion for PU system.

The maximum values of ${\alpha }_p$,${\alpha }_m$, and ${{\alpha }_p}_{\mathrm{\infty }}$ obtained respectively from the curing rate curve, y(α), and z(α) functions, are summarized in Table 3. According to [38], the autocatalytic model with two reacting factors m, n is applicable when ${{\alpha }_p}_{\mathrm{\infty }}$ < 0.632. To account for both autocatalytic and non-autocatalytic reaction pathways, the Kamal–Sourour (KS) kinetic model was employed [42]. The KS model can be expressed as follows:

$${\displaystyle \frac{d\alpha }{dt}}=[k_1+k_2{\alpha }^m]{\left(1-\alpha \right)}^n$$

(19)

where $k_1$and $k_2$represent the rate constants associated with the chemical-controlled and autocatalytic reactions, respectively; $m$denotes the autocatalytic reaction order, and $n$ represents the overall reaction order with respect to the unreacted fraction.

Table 2. Critical Curing Parameters of PU Determined from DSC at Heating Rates of 3, 5, 10, and 20 °C/min.

Heating rate	αp	αM	αp∞
3 °C/min	0.325	0.201	0.421
5 °C/min	0.381	0.200	0.448
10 °C/min	0.407	0.200	0.450
20 °C/min	0.389	0.198	0.450

Table 2. Critical Curing Parameters of PU Determined from DSC at Heating Rates of 3, 5, 10, and 20 °C/min.

Heating rate	αp	αM	αp∞
3 °C/min	0.325	0.201	0.421
5 °C/min	0.381	0.200	0.448
10 °C/min	0.407	0.200	0.450
20 °C/min	0.389	0.198	0.450

Based on the global optimizer, the KS kinetic parameters were determined through robust nonlinear least-squares fitting using a soft L1 loss function [43]. This approach mitigates the impact of experimental noise in the DSC data, particularly at the extremities of the curves, resulting in stable and physically meaningful parameter estimates. Due to the rapid curing behavior of the PU system, the dataset collected at the highest heating rate was excluded from the fitting. To ensure computational stability and avoid negative values, the parameters were log-transformed. The robust loss optimization function was employed to minimize the residuals, with the objective function described using the soft L1 loss function, as defined below:

$$r_i=w_i[{\alpha }_{pred}\left(A_1,A_2,m,n,E_a\right)-{\alpha }_{exp}]$$

(20)

$$\mathit{min}F=\sum _i\rho \left({\displaystyle \frac{r_i}{\sigma }}\right)$$

(21)

where $r_i$ denotes the residual of the $i$-th data point, ${\sigma }_i$ is the corresponding noise scale, and $\rho$ is a robust loss function introduced to reduce the influence of outliers.

The fitting curves are shown in Figure 5 and compared with the experimental data. The fitted DoC curves closely matched the experimental results, with $R^2$ values listed in Table 3, all exceeding 0.99. Overall, the $R^2$ values for the curing rate were consistently greater than 0.95, indicating high accuracy.

Table 3. Correlation (R²) Between Degree of Cure and Curing Rate at Heating Rates of 3, 5, and 10 °C/min.

Heating rate	R2
	DOC	Curing rate
3 °C/min	0.997	0.953
5 °C/min	0.999	0.972
10 °C/min	0.999	0.981

Table 3. Correlation (R²) Between Degree of Cure and Curing Rate at Heating Rates of 3, 5, and 10 °C/min.

Heating rate	R2
	DOC	Curing rate
3 °C/min	0.997	0.953
5 °C/min	0.999	0.972
10 °C/min	0.999	0.981

Based on the KS equation and fitted curves, the critical values of KS kinetics curing model are listed in Table 4.

3.2. Gelation

Time sweep oscillation tests were conducted at 50 °C with angular frequencies ($\omega$) of 1, 3, 10, 30, and 100 rad/s. Gel time can be estimated by measuring low-amplitude oscillatory shear and tracking changes in storage modulus ($G\text{'}$), loss modulus ($G″$), and complex viscosity (${\eta }^{\mathrm{*}}$). To identify the gelation point for the PU system, three methods were considered: the crossover point between storage and loss modulus, the frequency-independent intersection of the dissipation factor $\left(\mathit{tan}\left(\delta \right)={\displaystyle \frac{G^{\mathrm{\text{'}}\mathrm{\text{'}}}}{G^{\mathrm{\text{'}}}}}\right)$, and the moment when viscosity becomes practically infinite [44]. For a typical thermoset, viscous behavior dominates in the early stages of curing ($G″$ > $G\text{'}$). As the degree of crosslinking increases, there is often a rapid rise in elasticity, which begins to dominate the system ($G\text{'}$ > $G″$). The intersection, $G\text{'}$ = $G″$, is taken as the apparent gel point [45,46]. In some cases, no clear intersection point is observed within the measurement time window. This may be attributed to the characteristic timescale $\left(\tau ={\displaystyle \frac{1}{\omega }}\right)$, which determines how long the material has to respond within each oscillation cycle (i.e. a shorter timescale means the material is “probed” more frequently). At low angular frequencies, the large timescale enables polymer chain relaxation, allowing the loss modulus to dominate longer into the cure. In such cases, gelation can be assessed using criteria based on the dissipation factor or complex viscosity. Of the five angular frequencies tested, crossover between the storage and loss modulus occurred only during the 100 rad/s trial, at approximately 200 s. This is not taken to be the gel point – rather, it is attributed to the very short deformation timescale ($\tau =0.01$s) resulting in mildly elastic-dominant behavior as the polymer network rapidly began forming during the early stages of curing.

During the first 400-500 seconds of curing, the low magnitudes and close proximities of the storage and loss moduli in Figure 6 (a), coupled with experimental noise, led to significant instability in the dissipation factor for all trials. When computing $\mathit{tan}\left(\delta \right)$ for the curves in Figure 6 (b), $G\text{'}$ and $G″$ were first smoothed using a 3^rd order Savitzky-Golay filter with a frame length of 21. Winter-Chambon theory states that, at the gel point, the storage and loss moduli grow proportionally as a function of frequency [47]. Therefore, the gel point can be identified by determining when the dissipation factor becomes independent of angular frequency. In Figure 6 (b), the most notable intersection occurs between 600-625 s, for which four of the five curves intersect. The 10 rad/s trial intersects most of the other curves slightly earlier, between 475-550 s. Thus, the $\mathit{tan}\left(\delta \right)$ crossover method indicates a gelation window of 475-625 s, during which the true gel point is likely to exist.

The final method is based on the premise that fully crosslinked networks exhibit a complex viscosity which diverges toward infinity. This method often relies on visual cues to identify the gel point, such as spikes or plateaus across various frequencies. No obvious features are present in Figure 6 (c) to indicate the gel point via this method – in linear or log scale. When mapping the selected points to the DoC curve under 50 °C in Figure 6 (d), the gelation point occurred around DoC = 0.53.

3.3. Glass Transition Temperatures

Three glass transition temperatures $T_{g0}$, $T_{g,gel}$, $T_{g\mathrm{\infty }}$were measured for uncured resin, at gelation point, and cured resin, respectively. To determine $T_g$as a function of conversion $\left(\alpha \right)$for a fast-curing PU system, an integrated approach was employed [10]. The procedure combines isothermal curing within a DSC, rapid quenching, and subsequent dynamic scanning. The conversion $\left(\alpha \right)$ was calculated based on the residual enthalpy measured during the DSC analysis. All time references were synchronized to the DSC cell clock, effectively eliminating uncertainties related to sample preparation. The isothermal DSC test was conducted at 100 °C for 30 min, and the corresponding conversion curve is shown in Figure 7 (a). The PU system reached full cure at approximately 30 min under isothermal conditions. The time–temperature–transformation (TTT) diagram was constructed based on the Di Benedetto equation. The fitted glass transition temperature as a function of time is shown in Figure 7 (b), with the fitted constant $\lambda$of approximately 0.31. Similar calorimetry-based $T_g\left(\alpha \right)$ characterization and construction of TTT-type processing diagrams have been reported for fast-curing polyurethane adhesives to support model calibration and process-window definition [14].

3.4. Molecular Dynamic Simulation capturing T_g

To compare the MD model with experiment, the unreactive PU system was heated to 100 °C as previously described. The DoC vs step time curve is plotted by Figure 8 (a). Crosslinking was rapid in the beginning, followed by a gradual slowing until it reached 0 due to the limit of the reactive functional groups. The curing rate from MD simulation was normalized and compared to the curing rate with the experimental curve, which is shown in Figure 8 (b). The curing rates were mapped to each DoC point, the experimental curve shows a higher reaction rate happened in the beginning of the reaction compared to the MD curve, which shows the isocyanate–hydroxyl reaction is highly reactive, with reaction sites achieving instantaneous and complete contact [48,49]. By comparing the curing-rate trends of the two groups, it was observed that chemical reactions dominate the cross-linking process during curing. To evaluate the $T_g$ of the MD model, 10 DoC points were selected: 0.2, 0.3, 0.4, 0.65, 0.7, 0.75, 0.8, 0.85, 0.95 and 0.97. $T_g$was determined from the heating-induced specific volume temperature curves by identifying the intersection of the linear fits to the glassy and rubbery regimes and 5 curves of different DoC models are presented by Figure 8 (c). As the degree of cure increases, the specific volume of the polymer systematically decreases, reflecting the progressive formation of a denser crosslinked network and the reduction of free volume.

The resulting $T_g$– $\alpha$ relationship is shown in Figure 8 (d). The extracted $T_g$data were fitted using the DiBenedetto equation, with the initial and fully cured $T_g$values fixed. The fitted curve shows excellent agreement with the MD data, yielding a high coefficient of determination ($R^2=0.97$), indicating that the DiBenedetto model effectively captures the evolution of $T_g$with DoC.

A GPR-based predictive framework was implemented to robustly model the DiBenedetto relationship, while explicitly accounting for experimental and MD noise through additive noise terms.

3.5. Gaussian Process Regression Prediction

To fit the predictive model, in total 16 data points were combined from experimental and MD data, a cross validation was performed by k-fold cross validation method (k = 5), and the fitted curve is plotted as the blue dash line in Figure 9. The noise-based weighting controls both the data influence on the mean prediction and the width of the confidence interval. The best global fitted λ = 0.275 was selected by cross validation with the lowest RMSE in the defined region. This model captures the globally monotonic increasing behavior of T_g with DoC and strictly satisfies the predefined physical endpoints at DoC = 0 and 1.

By superimposing GP residuals onto the DiBenedetto baseline, the proposed model (solid orange curve) effectively corrects these local deviations (systematic and random noises) while preserving the physically meaningful trends imposed by the baseline and hard boundary constraints. The GPR correction primarily acted on regions where experimental and MD data persistently deviate from the baseline between 0.3 and 0.6, thereby enhancing data consistency across the entire DoC range without introducing non-physical oscillations. Notably, the model avoids overfitting sparse regions of the dataset, as evidenced by the final curve's smoothness and the absence of abnormal curvature near the endpoints. The wider 95 % confidence interval bands observed at low and intermediate degrees of cure arise from sparse experimental measurements and the relatively higher noise assigned to the MD data, where the automatically selected noise levels were 0.4 for MD and 0.4 for experiments. The predicted curve captures the theory-guided DiBenedetto trend well, and the model incorporates a reliable treatment of measurement and simulation noise.

4. Conclusions

In this study, the curing behavior of the PU system was characterized using thermal and rheological measurements conducted by DSC and a rheometer. The apparent E_a determined using the isoconversional method was 46.12 kJ/mol. A KS curing kinetics model was established to describe the curing evolution. In addition, the gelation point was identified based on combined DSC and rheological measurements, corresponding to a DoC of approximately 0.53.

Due to the rapid curing behavior of the PU system, the T_g was determined using a combination of experimental measurements and MD simulations. At the molecular scale, MD simulations were performed using LAMMPS, from which the corresponding T_g values were extracted at different DoC. GPR was employed as a bridging framework to integrate experimental DSC data and MD results, with crossvalidation used to ensure predictive robustness. A topological crosslinking approach was implemented in LAMMPS under different temperature conditions, from which a global DiBenedetto parameter of λ = 0.28 was determined. GPR enabled reliable forecasting the trend of T_g – α combined MD and experimental results with CI regions indicate DoC ranges dominated by MD data and limited experimental constraints. Furthermore, GP residuals effectively capture systematic, physics-driven biases in the T_g relationship, biases that cannot be fully described by the DiBenedetto model alone. In the T_g – α relationship, systematic, physics-driven biases exist that cannot be fully described by the DiBenedetto model alone.