A computational process modeling framework, informed by accurate material characterization, is presented for virtual manufacturing of wind energy thermoset composites . Process modeling simulations of composite microstructures are carried out to predict the in -situ matrix property evolution and performance- altering residual stress generation. To achieve this, a comprehensive material characterization effort is carried out . Consequently, a novel material property dataset for a widely -used wind energy thermoset system is generated as a function of the temperature and degree of cure. Informed by these material propert ies, the ability of the process model to reliably estimate manufacturing -induced residual stresses is highlighted. For a prescribed cure cycle, the in -situ elastic modulus evolution , chemical and thermal strains, and random fiber distribution are shown to significantly influence the residual stress generation. The results also show that a full process modeling analysis that includes the complete cure cycle (instead of the standard approach of just considering the post -processing cool -down cycle) is necessar y to accurately predict manufacturing -induced residual stresses. Keywords: A. Thermosetting resin B. Thermo mechanical B. Residual/internal stress E. Cure 3 1. Introduction Polymer matrix composites (PMCs) are finding widespread applications as primary load-bearing components in the wind energy industry due to their desirable properties and high strength-to -weight ratios [1]. However, manufacturing such comp osites remains a challenging task. The chemo -rheological and thermo -mechanical properties of the constituent matrix evolve as cur ing progresses. Several multi -physical phenomena , such as chemical shrinkage, exothermic heat generation , and mismatch between constituent material properties, further influence the in -situ matrix response and result in performance-altering residual stress build -up, deformation , and damage during manufacturing [2]–[9]. The extent to which these residual stresses alter the composite performance is not fully understood . Process modeling simulations informe d by accurate and comprehensive material characterization can fill this knowledge gap . The resulting process model can be used to optimize composite manufacturing to improve the performance of composite parts for wind energy applications . Process modeling leverages robust computational analysis tools to enable accurate prediction of the composite microscale response when subjected to various thermo-mechanical loads [2]–[8], [11]– [17]. Numerous computational studies, based on aerospace-grade composites , have reported phenomenological and constitutive relations to predic t the evolution of matrix cure [18], estimate the process-induced residual stress generation [19]–[28], and evaluate its effect on the processed composite performance [3]–[6], [8], [11], [12], [29] . However, d ue to the lack of a complete cure- and temperature -dependent material property data set, such studies typically rely on room temperature data or employ deterministic matrix properties for their analysis. Consequently, the prediction accuracy of the numerical methods presented in these studies strongly depends on the quality of the input material propert y dataset [30]. In-depth knowledge of the cure- and temperature- dependent material behavior is critical to fully understand the development of residual stresses within composite materials and accurately quantify their effect on composite part performance [10]. Few studies have reported comprehensive material 4 characterization methodologies to accurately and reliably determine the matrix curing kinetics and the development of thermo-mechanical properties as a function of the degree of cure and temperatu re [10], [30]–[36]. However, t hese studies are limited to the characterization of aerospace-grade resin systems. Similar material property dataset s are extremely scarce for commercial systems used in wind energy applications . Consequently, the implementation of robust process modeling techniques to optimize wind blade manufacturing proves challenging and thus remain s unexplored. Therefore, the need for a comprehensive and reliable material property dataset is paramount to facilitate process modeling of wind energy composites for their optimized manufacturing and improved performance and durability. In this study, a robust computational process modeling framework, informed by comprehensive material characterization , is presented to predict residual stress generation in composite microstructures during manufacturing. Using established characterization procedures typical to aerospace-grade materials, this work presents a novel material property dataset fo r a widely -used wind energy material system as a function of the degree of cure and temperature. For ease of implementation within the process modeling framework , resin -specific material models are developed that establish the dependence of the resin therm o-mechanical properties on the degree of cure and temperature. Subsequently, finite element (FE) -based process modeling analyse s are carried out to predict residual stress generation in composite microstructures subjected to an arbitrary temperature profil e. A separate linear -elastic thermal cooldown analysis is performed using constant room temperature material properties. The end-of-cure residual stresses from the process model and the thermal cooldown are compared . The sensitivity of residual stress generation to the material property evolution during manufacturing is thus evaluated. The manuscript is organized as follows: the material characterization methodology is outline d and the relevant results are di scussed in Section 2, the computa tional mode l generation and process modeling analysis of composite microstructures are detailed in Section 3, and the m ain conclusions are elaborated in Section 4. 5 2. Material and Experiment al Methods 2.1. Material and Manufacturing The present research focuses on a low viscosity thermosetting epoxy system developed by Hexion Inc. for wind energy applications. The system, commercially sold as EPIKOTETM Resin MGSTM RIMR 135 with the EPIKURETM Curing Agent MSGTM RIMH 1366 ( henceforth referred to as RIM R135/H1366), is a two -part formulation consisting of a base epoxy and an amine -based hardener. This system has a low processing temperature and long gel time, making it an ideal choice for infusion processes popularly used in wind turbine blade manufacturing. For all the characterization techniques discussed in this section , the resin mixture was prepared at room temperature by thoroughly mixing the RIM R135/H1366 in a stoichiometric ratio of 100:30 parts by weight for two minutes. The mixture was then degassed for a minimum of twenty minutes. Depending on the characterization technique, the uncured mixture was either directly tested for specific properties or injected into molds to prepare the test specimens. Test specimens were cured at a manufacturer -recommended temperature of 80 °C. 2.2. Cure kinetics The resin cure kinetics were measured with a D ifferential Scanning Calorimeter (DSC Discovery, TA Instruments). For these tests, a small amount of uncured mixture ( 6 ± 2.5 mg) was sealed in a hermetic aluminum pan. Before each DSC scan, the pan was equilibrated to - 10°C for two minutes to avoid premat ure curing. Non-isothermal measurements were performed to measure the total exothermic heat of the reaction 𝐻𝐻T. During these tests, specimens were heated at a constant temperature rate of 10 °C/min over a temperature range of -10 °C to 300° C. The total heat of reaction was computed by integrating the area under the normalized heat flow versus temperature plot shown in Figure 1a. 6 The relationship between the heat release rate d 𝐻𝐻/d𝑡𝑡 and the cure time 𝑡𝑡 for a given isothermal cure temperature 𝑇𝑇 was established with isothermal DSC measurements. For these tests, uncured resin specimens were rapidly heated from the equilibration tem perature ( 𝑇𝑇eq= -10°C) to 50 °C, 60 °C, and 70° C (typical processing conditions based on manufacturer recommendations) and cured for 720 minutes ( the time at which full cure is expected). The normalized heat flow versus time plots for the three temperatures, shown in Figure 1b, were used to compute the instantaneous d𝐻𝐻/d𝑡𝑡 and the partial heat of the reaction 𝐻𝐻 pr up to 𝑡𝑡. The rate of reaction d 𝜙𝜙/d𝑡𝑡 and the degree of cure 𝜙𝜙 were determined using d𝜙𝜙 d𝑡𝑡=1 𝐻𝐻Td𝐻𝐻 d𝑡𝑡 (1) 𝜙𝜙=1 𝐻𝐻T�d𝐻𝐻 d𝑡𝑡d𝑡𝑡𝑡𝑡 0=𝐻𝐻pr 𝐻𝐻T (2) For each value of 𝑇𝑇, d𝜙𝜙/d𝑡𝑡 was expressed as a function of 𝜙𝜙 using a phenomenological autocatalytic semi -empirical kinetic model developed by Kamal et al. [18], d𝜙𝜙 d𝑡𝑡=(𝐾𝐾1+𝐾𝐾2𝜙𝜙𝑚𝑚)(1−𝜙𝜙)𝑛𝑛 (3) 𝐾𝐾i=𝐴𝐴iexp�−Δ𝐸𝐸i 𝑅𝑅𝑇𝑇∗� i=1,2 (4) where 𝐾𝐾1 and 𝐾𝐾2 are Arrhenius rate functions , 𝑚𝑚 and 𝑛𝑛 dimensionless and real modeling parameters, 𝐴𝐴1 and 𝐴𝐴2 are pre-exponential factors , Δ𝐸𝐸1 and Δ𝐸𝐸2 are activations energies, 𝑇𝑇∗ is the absolute temperature in Kelvin, and 𝑅𝑅 is the gas constant. For each isothermal cure temperature, d 𝜙𝜙/d𝑡𝑡 and 𝜙𝜙 were fitted to Equation (3) using least-squares non-linear regression. Exponents 𝑚𝑚 and 𝑛𝑛 (constants for all isothermal measurements) were determined through trial -and-error in order to obtain the best fits. This yielded the 7 relevant rate functions 𝐾𝐾i for each isothermal cure temperature. Subsequently, 𝐴𝐴i and Δ𝐸𝐸i were obtained from a linear fitting of the ln(𝐾𝐾i) versus 1/𝑅𝑅𝑇𝑇∗. These kinet ic constants are summarized in Table 1. The derived model predictions were superimposed onto the experimental d𝜙𝜙/d𝑡𝑡 versus 𝜙𝜙 curves computed using Equations (1) and ( 2) in Figure 2a. A strong match was observed between the experiments and model predictions. To gauge the model’s accuracy, p artially cured specimens were tested under dynamic conditions to determine their degrees of cure. The kinetic model was very accurate in predicting the final cure state of these specimens, as shown in Figure 2b. Figure 1: Normalized heat flow versus (a) temperature plot from a non-isothermal DSC measurement and (b) time plot from three independent isothermal DSC measurements. Table 1: Summary of the kinetic constants determined from isothermal DSC measurements of RIM R135/H1366 epoxy system. Property Value Unit Exponent 𝑚𝑚 0.4 [-] 𝑛𝑛 1.5 [-] Rate constant 𝐴𝐴1 3.6 × 109 [sec-1] 𝐴𝐴2 0.01246 [sec-1] Activation energy Δ𝐸𝐸1 85.3 [kJ/mol] Δ𝐸𝐸2 11.1 [kJ/mol] 8 Figure 2: Comparison of experimental data with the predictions from t he autocatalytic cure kinetic model. The p lot of (a) the rate of reaction versus degree of cure computed from isothermal DSC measurements and (b) degree of cure and temperature versus time obtained from non -isothermal DSC measurements. 2.3. Glass Transition Temperature The resin glass transition temperature 𝑇𝑇g was determined using dynamic DSC. During the test, the normalized heat flow was measured as a function of the test temperature. A discontinuity in the slope of the normalized heat flow was regarded as the glass transition temperature. For these tests, uncured resin specimens were sealed in hermetic pans and cured under isothermal conditions at 80° C for several time intervals to attain intermediate cure levels. Following the isothermal cure, the specimens were rapidly cooled to -70° C to stop further curing. The specimens were subsequently subjected to a dynamic temperature sweep from -70° C to 300° C at a temperature rate of 10 °C/min. Figure 3a shows the normalized heat flow versus temperature plot for partially cured specimen s. The mid -point of the step -change in the normalized heat flow versus temperature plot (marked with arrows) was regarded as 𝑇𝑇g𝜙𝜙. The residual heat of the reaction 𝐻𝐻rs following the glass transition phase was used to determine the degree of cure using 9 𝜙𝜙=1−𝐻𝐻rs 𝐻𝐻T (5) The evolution of the glass transition temperature with the degree of cure is presented in Figure 3b. This evolution of 𝑇𝑇g𝜙𝜙 with 𝜙𝜙 was modeled with the DiBenedetto equation [37], 𝑇𝑇g𝜙𝜙−𝑇𝑇g0 𝑇𝑇g1−𝑇𝑇g0=𝜆𝜆𝜙𝜙 1−(1−𝜆𝜆)𝜙𝜙 (6) where 𝑇𝑇g0= -53.2° C and 𝑇𝑇g1= 85.6°C are the glass transition temperatures of the uncured and fully cured resin, 𝑇𝑇g𝜙𝜙 is the glass transition temperature for an intermediate degree of cure 𝜙𝜙 , and 𝜆𝜆 is a dimensionless fitting parameter used to obtain the best fit. For this material system, 𝜆𝜆=0.2373 . This model accurately predict ed the evolution of the glass transition temperature with an 𝑅𝑅2-value of 0.999. The glass transition temperature s obtained by two other methods, t hermo -mechanical analysis (TMA) and dynamic mechanical analysis ( DMA, discussed below) , were also compared to the predicted temperature as shown in Figure 3b. The DiBenedetto equation predictions of the 𝑇𝑇g𝜙𝜙 agreed well with the measured values from the TMA and DMA methods. The 𝑇𝑇g values between the testing methods were expected to differ margin ally since the glass transition passage does not affect the thermodynamic and mechanical properties at the same point, but rather occurs over a range of temperatures [34], [38]. 2.1. Gelation The time to gelation 𝑡𝑡 gel was measured with a rotational rheometer (ARES -G2, TA Instruments). The uncured resin mixture, prepared at room temperature, was injected between the parallel plates with an initial gap of ℎ0= 1.5 ± 0.05 mm. The mixture was then rapidly heated to 80 °C. To avoid temperature -induced volumetric changes in the test specimen, the tests were carried out under isothermal conditions at 80 °C. 10 Figure 3: Plot of (a) normalized heat flow versus temperature from dyn amic DSC measurements of partially cured specimens and (b) glass transition temperature versus degree of cure computed as the mid -point of the step -change , model fitting obtained from the DiBenedetto equation. For these tests, the rheometer was set in gap control mode, which mean t the instrument maintain ed a constant gap of ℎ0 for the duration of the test. As the resin cured , it was subjected to an oscillating shear strain 𝛾𝛾 with a maximum value of 200% and an oscillating frequency of 1 Hz. The storage modulus 𝐺𝐺′ and the loss modulus 𝐺𝐺′′ are components of the complex shear modulus 𝐺𝐺, which relates the applied 𝛾𝛾 to the measured shear stress 𝜏𝜏. Their evolution, as a function of the 𝑡𝑡, for one isothermally -cured test specimen is shown in Figure 4. The crossover between 𝐺𝐺′ and 𝐺𝐺′′, which indicates that the resin has sufficiently solidified to transfer stresses, was used to determine the resin gel point. This crossover occurred approximately at 𝑡𝑡gel= 38 minutes of isothermal curing at 80° C. Considering the time -temperature history of the rheology test, the degree of cure corresponding to 𝑡𝑡gel was calculated fr om the kinetic model given by Equation ( 3). The average gel point from the three tests 𝜙𝜙gel= 0.78. Literature studies on other RTM resin systems have reported similar gel point values [31], [34], [39]– [42]. 11 Figure 4: Plot showing the evolution of the storage and the loss moduli as a function of cure time. The crossover between storage and loss moduli is defined as the gel point. 2.2. Chemical Shrink age The post -gelation volumetric chemical shrinkage 𝜖𝜖sh𝜙𝜙 was measured using the same rheom eter setup as the gel point measurement. The test procedure for the shrinkage measurement was divided into two steps. The first step was identical to the gel point measurement where the instrument operated in gap control mode. During this step, the liquid resin was allowed to cure under zero normal force 𝐹𝐹N when 𝑡𝑡<𝑡𝑡gel. The complex viscosity measured during this step was used to switch from gap control mode to force control mode when 𝑡𝑡=𝑡𝑡gel. During the second step, the instrument applie d a constant 𝐹𝐹N= 0.1 N to maintain contact between the shrinking specimen and the plates. Additionally, a maximum torque of 100 mN was applied at a 30 Hz frequency. The steady decrease in the gap distance (ℎ) resulting from chemical shrinkage in the specimen, was monitor ed by the instrument. The decrease in ℎ relative to the 𝑡𝑡 is presented in Figure 5. When 𝑡𝑡<𝑡𝑡gel, ℎ remain ed constant with a value of ℎ 0 indicating a gap control mode operation. When 𝑡𝑡 =𝑡𝑡gel, the instrument switch ed to a force 12 control mode and applie d a 𝐹𝐹N= 0.1 N. This caused a sudden decrease in ℎ , as seen in Figure 5. Upon stabilizing, ℎ steadily decreased as the specimen volume reduced due to chemical crosslink formation. The specimens were allowed to cure for an additional 120 minutes after gelation. By then, the specimens were expected to be fully cured and additional change in ℎ was found to be insignificant. Neglecting the decrease in ℎ due to the application of 𝐹𝐹N, the post-gelation volumetric chemical shrinkage 𝜖𝜖sh was computed using 𝜖𝜖sh=�1+1 3 �ℎ−ℎ0 ℎ0��3 −1 (7) This relationship assumes the in -plane strains in the specimen are zero and the material is incompr essible (𝜈𝜈= 0.5) [40]. The chemical shrinkage as a function of the cure time for three isothermally cured test specimens is shown in Figure 6a. The average post -gelation volumetric shrinkage for a fully cured specimen 𝜖𝜖sh1= 2.1 ± 0.07% . This value was within the typical range of shrinkage values observed in the literature (2 – 7%) for epoxy systems [34], [43], [44] . Based on the time-temperature history of the test, the degree of cure was computed using Equation (3). The chemical shrinkage as a function of the degree of cure is presented in Figure 6b. 𝜖𝜖sh𝜙𝜙 was modeled as a linear function of 𝜙𝜙 𝜖𝜖sh𝜙𝜙=�0, 𝜙𝜙≤𝜙𝜙gel 𝛽𝛽(𝜙𝜙− 𝜙𝜙gel),𝜙𝜙>𝜙𝜙𝑔𝑔𝑔𝑔𝑔𝑔 (8) where 𝛽𝛽 is the shrinkage coefficient determined using the least square linear regression of chemical shrinkage versus time curves. For this material system, 𝛽𝛽 =0.111. 2.1. Coefficient of Thermal Expansion The post -gelation coefficient of thermal expansion (CTE) 𝛼𝛼i𝜙𝜙 of the RIM R135/H1366 epoxy was measured using a TMA (Q400, TA Instruments). The TMA employs a highly sensitive expansion probe to measure the through-thickness dimension changes Δℎ in the specimen relative to the temperature. 13 Figure 5: Plot showing the variation in the gap between the parallel plates (circular, square, and triangular symbols in blu e) during various phases of the test and the evolution of the degree of cure (dot -dash curve in orange). Figure 6: Evolution of the measured volumetric chemical shrinkage under isothermal conditions as a function of (a) time and (b) degree of cure. The test specimens were prepared in an open mold. The d egassed mixture was injected into the mold, which was then transferred into a convection oven that was programmed to run the manufactu rer- recommended cure cycle. Specimens were allowed to cure in the oven until the target degrees of cure were 14 achieved. The specimens were removed from the oven and rapidly cooled to arrest further curing. The specimens were then cut into 3- 4 mm cubes and polished to achieve top and bottom sur faces that were flat and parallel for the probe to rest on [45]. The specimen s were placed in the TMA furnace which was equilibrated at 20°C and heated at a constant rate of 5 °C/min to 130 °C. A n 𝐹𝐹N= 0.01 N was applied for the entire duration of the test to maintain contact with the specimen. As the specimen was heated, the probe measured the through-thickness dimension change Δℎ in the specimen as a function of the temperature 𝑇𝑇 . The CTE is defined as the relative change in the specimen dimension Δℎ/ℎ0 with respect to the change in the temperature Δ𝑇𝑇 [34], [45] . This value was obtained as the slope of the relative dimension change versus the temperature plot generated by the TMA. Figure 7a shows one such plot for a fully cured specimen. A distinct linear region between 40 °C and 80° C (region I) was observed , followed by a transition region between 80 °C and 95 °C (gray area), where the slope changed . After this transition region, a second linear region was observed between 95 °C and 11 0°C (region II). The transition region represent ed the material’s physical transition from a glass y state to a rubbery state as it was heated, and the inflection point was defined a s the 𝑇𝑇g. The slopes of the plot before and after the transition region corresponded to the pre - and post -transition CTE. The specimen register ed a 𝑇𝑇g1= 91.2° C. The pre -transition CTE was 𝛼𝛼I1= 7.39 × 10-5/°C while the post -transition CTE was 𝛼𝛼II1= 17.1 × 10-5/°C. The variation in the post-gelation CTE is illustrated as a function of 𝜙𝜙 in Figure 7b. For 𝑇𝑇<𝑇𝑇g𝜙𝜙, 𝛼𝛼I𝜙𝜙= (7.92±0.538 )×10−5/°C was obtained. This value did not change with 𝜙𝜙 . By contrast, above 𝑇𝑇 g, the material registered a CTE value ranging from 𝛼𝛼II𝜙𝜙=(18.95±0.86)×10−5/°C for 𝜙𝜙=0.9 to 𝛼𝛼II𝜙𝜙= (16.86±0.21)×10−5/°C for 𝜙𝜙=1. A 11% drop was observed in the material CTE with an increasing degree of cure which was attributed to crosslink network formation, which restricts molecular mobility as cure progresses. The variation in the post-gelation CTE with the degree of cure and glass transition temperature was defined using 15 𝛼𝛼i𝜙𝜙=� 7.93×10−5, 𝑇𝑇≤𝑇𝑇g𝜙𝜙 �𝑎𝑎−𝑏𝑏�𝜙𝜙− 𝜙𝜙gel��×10−5, 𝑇𝑇>𝑇𝑇g𝜙𝜙 (9) where 𝑎𝑎 and 𝑏𝑏 are dimensionless fitting parameters obtained using linear regression. For this material system, 𝑎𝑎= 21.21 and 𝑏𝑏= 18.25. Figure 7: (a) Relative dimension change as a function of temperature for a fully cured specimen and (b) variation in the coefficient of thermal expansion with the degree of cure before and after the glass transition region. 2.2. Elastic Modulus The variation of the post -gelation elastic modulus 𝐸𝐸T𝜙𝜙 with the processing temperature and the degree of cure was characterized with a DMA (Q800, TA Instruments). The test specimens were prepared in a closed mold. T he degassed mixtures of resin were injected into the mold, which was then transferred into a convection oven that was programmed to run the manufacturer -recommen ded cure cycle. Specimens were allowed to cure in the oven until target degrees of cure were achieved. The specimens were removed from 16 the oven and rapidly cooled to arrest further curing. The specimens were then cut and polished to conform to the ASTM D70 28-07 standard test specimen geometry [46]. During the test, fully cured specimens were mounted in a single cantilever beam configuration inside the DMA chamber. The specimens were first tested isothermally at 20 °C for 15 minutes by applying an oscillating displacement of 10 𝜇𝜇 m at a 1 Hz frequency. Following this, a temperature sweep was performed, where the temperature was linearly increased to 120 °C at a rate of 3 °C/min. For the entire duration of the test, the instrument measure d the storage modulus 𝐸𝐸 ′ and the loss modulus 𝐸𝐸′′. The complex elastic modulus 𝐸𝐸 was computed using 𝐸𝐸=𝐸𝐸′+𝑖𝑖𝐸𝐸′′ → 𝐸𝐸=√𝐸𝐸′2+𝐸𝐸′′2 (10) The variation in the complex elastic modulus with the temperature for a fully cured specimen (𝜙𝜙=1) is shown in Figure 8. For the initial 20 minutes of the test, the temperature was held constant, and the specimen register ed an 𝐸𝐸201= 2482 MPa. As the temperature gradually increased , the material elastic modulus beg an to drop. A sharp drop in the complex elastic modulus was observed close to 85 °C, where the glass transition passage was expected. Past the glass transition window, the material register ed a significantly low modulus. Figure 8: Variation in the matrix elastic modulus with time , temperature, and degree of cure. 17 Partially cured specimens were tested under similar test conditions and the evolution of the ir elastic modul i as a function of the degree of cure and temperature is also presented in Figure 8. It is clear from this plot that all specimens exhibit ed an elastic mo dulus close to 2400 MPa at 20 °C. Subsequently, the elastic modulus beg an to drop as the temperature approached the specimen 𝑇𝑇g𝜙𝜙 corresponding to its cure level. This trend suggest ed that the post -gelation elastic modulus of the vitrified resin ( 𝑇𝑇<𝑇𝑇g𝜙𝜙) had a low sensitivity to the degree of cure. However, as the resin temperature approache d its glass transition temperature, which is highly dependent on the degree of cure (see Figure 3b), its elastic modulus drop ped significantly . The post -transition elastic modulus, when the resin is in the rubbery state ( 𝑇𝑇≥𝑇𝑇g𝜙𝜙), was highly sensitive to not only the temperature but also the degree of cure (Figure 8). Therefore, the effect of temperature and the degree of cure on the elastic modulus should be carefully accounted for during process modeling analysis. To model th is complex dependence of the elastic modulus on the degree of cure and temperature, an approach consistin g of a dimensionless, normalized independent parameter and regression was introduced. The dimensionless parameter 𝜂𝜂 is defined as 𝜂𝜂≡1−𝑇𝑇∗ 𝑇𝑇𝑟𝑟∗ (11) where 𝑇𝑇𝑟𝑟∗ is the reference temperature in degrees Kelvin , which herein was assigned as the highest temperature for which data was available. The parameter 𝜂𝜂 was thus defined to be independent of units and provide d a modulus function with direct proportionality to the independent variable , that is, as 𝜂𝜂 increases, so does the modulus . Here, according to Figure 8, 𝑇𝑇𝑟𝑟∗ = 393.15 K (or 120° C) was utilized. The E versus 𝜂𝜂 curve was fitted to a sigmoidal function (Boltzman n equation) for all degrees of cure: 18 𝐸𝐸(MPa )=𝐴𝐴2(𝜙𝜙)+𝐴𝐴1(𝜙𝜙)−𝐴𝐴2(𝜙𝜙) 1+exp�𝜂𝜂−𝜂𝜂0(𝜙𝜙) d𝜂𝜂(𝜙𝜙)� (12) The data fitting was perfo rmed using the Origin Pro [47] analysis tool and t he phenomenological parameters 𝜂𝜂0,d𝜂𝜂,𝐴𝐴1 and 𝐴𝐴2 were all determined as a function of degree of cure : 𝜂𝜂0= −0.860𝜙𝜙+0.955 d𝜂𝜂= 0.070𝜙𝜙−0.041 𝐴𝐴1(𝜙𝜙)= −10−13exp(34.55𝜙𝜙 ) 𝐴𝐴2(𝜙𝜙)= 1471 .3𝜙𝜙+801 .39 (13) This set of equations can be used to co mpute 𝐸𝐸 at arbitrary cure levels and temperatures. 3. Computational Micromechanics Computational studies of composite microstructures traditionally assume that the thermo -mechanical response of composite material can be determined by analyzing a R epeating Unit Cell (RUC). A unit cell can be termed as a representative volume element (RVE) if the cell is large enough to capture the intended physical mechanisms that influence the composite response without any boundary effects [8]. A convergence analysis must be carried out to determine the appropriate size of the R VE for a particular microstructure. The results of a size effect study performed previously [8] were used for the current study . The micromechanical and process modeling details are discussed in the following sections. 3.1. RVE Generation In this study, RVEs comprising a random dispersion of glass fibers ( fiber diameter, 𝑑𝑑 f=14 𝜇𝜇m) in a RIM R135/H1366 epoxy matrix were generated with a random RVE generator that was developed in- house with MATLAB. Here, 40 fiber RVEs with a fiber volume fraction 𝑣𝑣f= 0.55 were used. Perfect bonding 19 was assumed between the fiber and the matrix. Five distinct renditions of the 40 fiber RVE were analyzed , as shown in Figure 9. The fibers were modeled as isotropic linear -elastic solids, the thermo -mechanical properties of which are summarized i n Table 2. The matrix was modeled as an isotropic material . The material -specific models described above were used to define the matrix behavior during curing. The process modeling and thermal cooldown analysis were carried out with the commercial FE code Abaqus/STANDARD, supplemented with user -written UMATHT and UMAT [48] subroutines . Each RV E was meshed with C3D8T elements (eight -node, brick elements with temperature degrees of freedom) [48]. Flat Boundary Conditions (FBC s) were chosen over P eriodic Boundary Conditions (PBCs) due to their efficient and low-cost implementation [8]. Figure 9: Five distinct realizations of 40 fiber composite microstructu res gene rated by the random RVE generator . Table 2: Constituent thermo -mechanical properties of E -glass fiber. Property Value Unit Density 𝜌𝜌f 2,550 [kg/m3] Elastic Modulus 𝐸𝐸11f=𝐸𝐸22f=𝐸𝐸33f 73,000 [MPa] Poisson’s ratio 𝜈𝜈12f=𝜈𝜈13f=𝜈𝜈23f 0.22 [-] Shear Modulus 𝐺𝐺12f=𝐺𝐺13f=G23f 30,000 [MPa] Coefficient of Thermal Expansion 𝛼𝛼11f=𝛼𝛼22f=α33f 5 × 10-6 [K-1] Thermal Conductivity 𝑘𝑘f 1.2 [W/m -K] Specific Heat 𝑐𝑐pf 800 [J/kg -K] 20 3.2. Process Modeling The process modeling replicated the complete manufacturing of the composites by accounting for (a) kinetic progression of the cure which was quantified by 𝜙𝜙 and (b) the evolution of the in-situ thermo - mechanical matrix properties which led to residual stress generation. The progression of the cure f or a prescribed temperature profile was defined by Equation (3) and the kinetic constants summarized in Table 1. A coupled temperature-displacement analysis was carried out in Abaqus/STANDARD with user - subroutine s UMATHT and UMAT [48]. The kinetic model was solved simultaneously with the 3 - dimensional Fourier heat transfer model to predict the progression of cure and the temperature distribution as a result of the exothermic heat of reaction. Figure 10 shows the prescribed cure cycle for virtual manufacturing of the composite microstructure along with the computed degree of cure. The evolution of the in-situ matrix properties and residual stress generation was modeled with a UMAT user subroutine [48]. For a given cure state, the instantaneous material properties were computed with the help of the material models described above. Residual stress generation was then modeled with an instantaneous linear -elastic constitutive model [4], [8] , 𝜎𝜎i(𝑡𝑡)=�𝐶𝐶ij(𝑡𝑡)��𝜖𝜖jtotal(𝑡𝑡)−�𝜖𝜖jth(𝑡𝑡)+𝜖𝜖jsh(𝑡𝑡)�𝛿𝛿j� where , �𝛿𝛿j=1 if j=1,2,3 𝛿𝛿j=0 if j>3 (14) where i and j are Voigt notation indices; 𝜖𝜖jtotal(𝑡𝑡), 𝜖𝜖jth(𝑡𝑡) and 𝜖𝜖jsh(𝑡𝑡) are the total, thermal, and chemical shrinkage strains, respectively; 𝐶𝐶ij(𝑡𝑡) is the stiffness matrix as a function of the time of cure; 𝜎𝜎i(𝑡𝑡) is the accumulated residual stress governed by the development of the in -situ matrix elastic modulus and the chemical and thermal strains experienced by the matrix. 21 Figure 10: (a) Evolution of the degree of cure as a function of the applied cure cycle. The evolution of the in-situ matrix elastic modulus in a representative RVE with the degree of cure and temperature is shown in Figure 11a. Figures 11 (b)-(d) show the development of residual stresses with the time of cure at three distinct points in the same RVE resulting from the elastic modulus evolution and volumetric strains experienced by the matrix . In Figure 11a, the pre -gelation epoxy matrix (𝜙𝜙≤𝜙𝜙gel) exhibits a very low elastic modulus. Therefore, no significant residual stress generation is observed in Figures 11 (b)-(d) during the pre -gelation phase, even though the matrix experienced thermal and shrinkage strains during the initial heating and isotherma l hold stages of the cure cycle, respectively. However, the formation of sufficient crosslink networks led to a rapid increase in the material elastic modulus post- gelation , as seen in Figure 11a. A high post -gelation elastic modulus combined with chemical shrinkage strain s during the isothermal cure stage and thermal strains during the cooldown stage of the cure cycle result ed in a maximum residual stress generation of 20.28 MPa at point 3 in the RVE , as seen in Figure 11(d). The representative RVE registered volume -average residual stress of 7.6 MPa. 22 Figure 11: (a) Evolution of the matrix elastic modulus as a function of the prescribed cure profile and degree of cure for process modeling analysis and development of residual stresses (max. principal) at (b) point 1, (c) point 2, and (d) point 3 within an RVE for process modeling and thermal cooldown analysis. The residual stresses were also sensitive to the local fiber packing. Figure s 11 (b)-(d) show a variation in the magnitude of the residual stresses at three distinct points within a representative RVE . Owing to the random fiber distribution, the end -of-cure residual stresses at three points within the RVE varie d between 12.95 MPa and 20.28 MPa . To further illustrate the influence of random fiber packing on residual stress generation, the contour plots of the end-of-cure maximum principal stress for each of the 40-fiber RVE s are presented in Figure s 12 (a)-(e). The figure s demonstrate a considerable difference in the distribution of 23 the residual stresses o wing to the unique fiber packing topology in each RVE . Between the RVEs, there was a maximum and volume -average end -of-cure residual stress of 20.45 ± 0.25 MPa and 7.56 ± 0.28 MPa , respectively . Figure 12: Contour plots of the end -of-cure residual stresses (maximum principal) obtained from (a) -(e) process modeling analysis and (f) -(j) thermal cooldown analysis. 3.3. Thermal Cooldown Analysis Simple thermal cooldown analyses are typically employed to estimate the end-of-cure residual stresses in composites when extensive material property characterization and computational modeling are not feasible. These analyses only consider the final cooling stage of the primary curing cycle during which the matrix material is expected to be fully cured and is gradually cooled from the processing temperature to room temperature . Consequently, t he progression of cure, matrix elastic property evolution, and chemical shrinkage strains are ignored in such analyses, and only differential thermal volumetric changes experienced by the constituent fiber and matrix are thus assumed to contribute to the total residual stress generation. A comparison of the end-of -cure residual stresses from the full process model described her ein and a simple thermal cooldown analysis can quantify the contributing sources that influence residual stress generation . 24 To perform a thermal cooldown analysis in FE, constant room temperature matrix properties corresponding to the full y cured state were allotted to the RVEs shown in Figure 9 with a prescribed linear cooldown from 80°C to 20°C. The residual stress evolution during the thermal cooldown analysis at the same three points in the representative RVE is overlayed on the process modeling results in Figures 11 (b)- (d). As seen in the figure, the residual stresses increased linearly at each point as the RVE temperature decreased from 80°C to 20°C. Stress values as high as 23.84 MPa were predicted in the representative RVE shown in Figure 11(d). Compared to the process modeling results at the same location (𝜎𝜎rsmax= 20.28 MPa), the thermal cooldown analysis overesti mated the maximum residual stresses by 17.5%. Furthermore, the volume -average residual stress over the entire representative RVE from the thermal cooldown analysis was 9.02 MPa. Compared to the process model (𝜎𝜎rsavg= 7.6 MPa) , the cooldown analysis overestimated the volume -average stress by 18.6% as well. This overestimation is further evident in Figure 12 which compares the end -of-cure residual stress contour plots from the simple thermal cooldown analysis ( Figure 12 (f)-(j)) to those from the full process model ( Figure 12 (a)-(e)). The thermal cooldown analyses predicted maximum and volume -average residual stresses of 25.2 ± 0.6 MPa and 9.03 ± 0.3 MPa, respectively. These values were 23.2 % and 19.4 % higher than the maximum and volume -average stresses reported by the full process model, respectively. From these results, it is clear that simple thermal cooldown analyses for this material fail to accurately predict the process-induced residual stress. In addition to ignoring the mechanical property evolution and chemical shrinkage during curing, the higher initial CTE , which gradually decreases as the cure progresses, may cause significant volumet ric expansion during the initial heating stage of the cure cycle, which was also not accounted for. Thus, a combination of l ow elastic modulus, high initial volumetric expansion, and chemical shrinkage result in non -negligible residual stress generation. T herefore, reliable predictions of residual stress development are not possible with the simple thermal cooldown analysis approach. With the full process modeling approach, the evolution of cure - and temperature-dependent material properties, as 25 well as thermal and chemical shrinkage strain predictions, are possible for any prescribed cure cycle. This implies that the process modeling analysis developed herein can be employed to tailor the prescribed cure cycles to mitigate residual stress generation or to achieve a desired internal stress state for enhanced composite part performance. 4. Conclusions This wor k expands the premise of composite process modeling to wind energy applications. A robust and reliable computational process modeling framework was presented to facilitate virtual manufacturing of composite microstructures and accurate residual stress prediction during their manufacturing. Using established characterization procedures typical to aerospace-grade materials, this work present ed a novel material property dataset for a wind energy material system as a function of the cure and temperature. Informed by the se accurate ly measured material properties, t he ability of p rocess model ing tools to reliably estimate the manufacturing -induced local residual stress distribution within composite microstructures was highlighted. Process modeling techniques such as those presented in this research can greatly benefit the efforts to reduce the generation of residual stresses or to achieve the desir ed internal stress state for enhanced composite part performanc e. Optimization studies can be performed with process models to assess the influence of each processing parameter , such as the heating rate, hold temperature, and cooling rate on the final stress state. This in turn will provide optimization opportunities and facilitate the establishment of rapid, reliable, and cost -effective manufacturing procedures. 26 Acknowledgments This material is based upon the work supported partially by the National Science Foundation under grant number IIP – 1362022 (Collaborative Research: IUCRC for Wind Energy, Science, Technology and Research) and the Air Force Office of Scientific Research and the National Science Foundation under grant number IIP – 1826232. Any opinions, findings and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation. The authors would li ke to thank NASA for their support of this research under grant number s 80NSSC19K1246 and NNX17AJ32G. The authors would like to acknowledge the WindSTAR Industrial Advisory Board (IAB) members Steve Nolet, Amir Salimi (TPI Composites Inc.) for their techni cal contributions and Paul Ubrich, Nathan Bruno, Mirna Robles (Hexion Inc.) for their technical insight to the project and providing the resin material .