Sensitivity Analysis and Multi-Objective Optimization Strategy of the Curing Profile for Autoclave Processed Thick Composite Laminates

To mitigate the risk of manufacturing defects and improve the efficiency of the autoclave-processed thick composite component curing process, parameter sensitivity analysis and optimization of the curing profile were conducted using a finite element model, Sobol sensitivity analysis, and the multi-objective optimization method. The FE model based on the heat transfer and cure kinetics modules was developed by the user subroutine in ABAQUS and validated by experimental data. The effects of thickness, stacking sequence, and mold material on the maximum temperature (Tmax), temperature gradient (ΔT), and degree of curing (DoC) were discussed. Next, parameter sensitivity was tested to identify critical curing process parameters that have significant effects on Tmax, DoC, and curing time cycle (tcycle). A multi-objective optimization strategy was developed by combining the optimal Latin hypercube sampling, radial basis function (RBF), and non-dominated sorting genetic algorithm-II (NSGA-II) methods. The results showed that the established FE model could predict the temperature profile and DoC profile accurately. Tmax always occurred in the mid-point regardless of laminate thickness; the Tmax and ΔT increased non-linearly with the increasing laminate thickness; but the DoC was affected slightly by the laminate thickness. The stacking sequence has little influence on the Tmax, ΔT, and DoC of laminate. The mold material mainly affected the uniformity of the temperature field. The ΔT of aluminum mold was the highest, followed by copper mold and invar steel mold. Tmax and tcycle were mainly affected by the dwell temperature T2, and DoC was mainly affected by dwell time dt1 and dwell temperature T1. The multi-objective optimized curing profile could reduce the Tmax and tcycle by 2.2% and 16.1%, respectively, and maintain the maximum DoC at 0.91. This work provides guidance on the practical design of cure profiles for thick composite parts.


Introduction
Thermosetting composites have been widely used in aerospace, automotive, civil engineering, and other fields due to their light weights, high specific strength, fatigue resistance, and strong designability [1][2][3][4]. The composites are usually produced by an autoclave with a determined cure cycle profile in order to initiate and sustain an irreversible cross-linking of the resin [5][6][7][8]. With the composite increasingly used in complex structures, the thickness would exceed the applicable range of curing parameters recommended by Manufacturer Recommended Cure Cycles (MRCC) [9][10][11]. It will inevitably lead to thermal gradients, overshoot, and insufficient curing or energy consumption [12,13]. Therefore, it is of great importance to analyze the curing process comprehensively and to develop an optimization approach.
Optimization strategies based on finite element numerical simulation are used as a general method to drive the appropriate cure profile for thick laminate. Most of the research is limited to single-objective optimization [14,15], and some used multi-objective optimization to strike a balance between process costs and product quality. For example, Yuan et al. [16] and Gao et al. [17] established a multi-objective approach to optimize the curing process for thick composites based on a multi-field coupled model with a surrogate model. The t cycle (cure time duration), ∆T (maximum temperature gradient), and DoC (degree of curing) using the optimal cure profile have been reduced by about 30.9%, 45.76%, and 16.88%, respectively, in comparison with the MRCC cure cycle. Tang et al. [18] introduced a multi-objective optimization method based on finite element simulation to control the t cycle and cure-induced defects of C-shaped composites. Compared with the original profile, the t cycle was shortened by 19% to 14,686 s. Similarly, Li et al. [19] proposed a method to optimize the fiber-reinforced composite injection molding process by combining the combined Taguchi response surface methodology and the NSGA-II approach. It indicated that NSGA-II was an effective method to solve the multi-objective optimization problem for the quality optimization of fiber-reinforced composite injection molding. Dolcum et al. [20] developed a novel approach based on the finite element method and a multi-objective genetic algorithm to optimize the cure profile for thick thermoset composites. The results showed that compared to the original curing profile, the optimized one led to approximately a 56% reduction of the maximum difference DoC, a 71% decrease in the maximum difference in T max , and a 33% reduction in t cycle . Struzziero et al. [21] also developed a method combining a finite element solution with a genetic algorithm to optimize the curing process of thick components. The optimized cured profile indicated improvements of about 70% in overshoot and a reduction in process time of about 14 h. However, the studies mentioned above rarely analyzed the material and structure effects on the thick composite curing process. Meanwhile, the problem that to what degree of curing profile parameters affected the simulation results also remains to be addressed.
The present work aims to comprehensively analyze the curing process of thick composites based on finite element simulation, parameter sensitivity analysis, and multi-objective optimization. Firstly, a FE model based on a heat transfer module and a cure kinetics module was developed for a thick laminate and validated by experimental data. The effects of composite thickness, stacking sequence, and mold material on T max , ∆T, and DoC were comprehensively discussed. Then, parameter sensitivity was tested to identify the critical process parameters which have significant effects on the curing results. Finally, a multi-objective approach to optimize the curing process for thick composites by integrating the RBF and NSGA-II algorithms was established. A decision-making method was also used to select the final optimal solution from the Pareto optimal set.

Thermo-Chemical Description of the Curing Process
The FE model included the heat transfer module and the cure kinetics module. The evaluation of the DoC and temperature is calculated by the thermo-chemical-coupled heat transfer module as in [22]: in which ρ c , C c and k ii (i = x, y, z) are the density, specific heat capacity, and thermal conductivity of the composite, respectively; V f is the fiber volume fraction in the composite. ρ r is the density of resin and H R is the total quantity of heat released from the curing reaction of a unit mass of the resin. T and t are the current temperature and curing reaction time, respectively. dα dt indicated the instantaneous DoC, which can be computed in the incremental step in FE analysis using the cure kinetics model of resin.
The curing kinetics model of resin is mainly divided into two kinds: the phenomenological model and the mechanism model. The phenomenological model used a single reaction to replace the whole reaction process; the mechanism model is more inclined to analyze the kinetic mechanism of the reaction process. Most researchers use phenomeno-logical models to describe the curing reaction process of resin, and the expression of the curing kinetic model can be seen in [23,24]: in which K is the activation energy; m and n are the first and second exponential constants, respectively; A is the pre-exponential coefficient; C and R are the diffusion constant and gas constant, respectively. α C is the temperature-dependent DoC; α C0 is the constant at T = 0 K, and α CT is DoC increasing ratio with temperature. The density of the composite ρ c and specific heat capacity of the composite C c can be calculated by the rule of mixture as in [25]: in which ρ f and C f are the fiber density and specific heat capacity, respectively. Similarly, the thermal conductivity of the composite k ii (i = x, y, z) can also be computed according to the rule of mixture as in [25,26]: For unidirectional ply as the transversely isotropic material, in-plane thermal conductivity of composites perpendicular to the fiber direction (k yy ) and in the thickness direction (k zz ) are assumed to be equal and can be calculated as in [25,27]: This research used the composite materials in Ref. [25] to conduct the finite element simulation. The thermal physical properties for 8552 resin and AS4 fiber as well as the cure kinetic constants for 8552 resin are listed in Table 1. length. The 8-node linear heat transfer brick (DC3D8) was for composite laminate and mold. The prescribed temperature boundary condition was set to be equal to the autoclave air temperature (Figure 1a). used include HETVAL, USDFLD, DISP, and FILM. HETVAL can define the reaction heat inside the composite; USDFLD is used to describe the curing degree field of the composite curing process; DISP is used to determine the temperature boundary conditions of the composite curing process, and FILM can determine the convective heat transfer boundary conditions of the composite curing process (Figure 1b).
The curing cycle was followed as in Ref. [25]. It can be divided into four stages: (a) the initial temperature was set to 25°C; (b) the temperature was increased to 110 °C with a rate of 2.0 °C/min (r1); (c) the dwell temperature of 110°C (T1) was held for 1 h (dt1); (d) the temperature was increased again to 180°C with a rate of 2.0 °C/min (r2); (e) the dwell temperature of 180°C (T2) was held for 2 h (dt2) and then (f) decreased to 25 °C at a rate of −2.0 °C/min ( Figure 1c).

Parameter Sensitivity Analysis
The typical curing profile of six parameters including the two heating rates (r1, r2), two dwell times (dt1, dt2), and two dwell temperatures (T1, T2) were extracted from the curing profile as the design variables.
Sobol sensitivity analysis is widely used in various fields in the industry to identify the most critical input variables and improve the accuracy and robustness of models. It was employed to quantify the impact of input variables on the output of a model or simulation and provided a measure of the relative importance of each input variable by decomposing the total variance in the model output into contributions [30,31]. Herein, the Sobol sensitivity analysis method was employed to investigate the effects of curing profile parameters on the results. The contribution of each parameter is derived based on variance decomposition. The system input-output function f(x) can be decomposed into a summary of increasing dimensions: in which f is the constant term, which is equal to the expectation value of the output. Each input variable x , i = 1,2, ⋯ n is randomly distributed in the range of [0, 1]. In this paper, the commercial finite element software ABAQUS and FORTRAN subroutines are used to simulate the curing process of composite materials. The subroutines used include HETVAL, USDFLD, DISP, and FILM. HETVAL can define the reaction heat inside the composite; USDFLD is used to describe the curing degree field of the composite curing process; DISP is used to determine the temperature boundary conditions of the composite curing process, and FILM can determine the convective heat transfer boundary conditions of the composite curing process (Figure 1b).
The curing cycle was followed as in Ref. [25]. It can be divided into four stages: (a) the initial temperature was set to 25 • C; (b) the temperature was increased to 110 • C with a rate of 2.0 • C/min (r 1 ); (c) the dwell temperature of 110 • C (T 1 ) was held for 1 h (dt 1 ); (d) the temperature was increased again to 180 • C with a rate of 2.0 • C/min (r 2 ); (e) the dwell temperature of 180 • C (T 2 ) was held for 2 h (dt 2 ) and then (f) decreased to 25 • C at a rate of −2.0 • C/min (Figure 1c).

Parameter Sensitivity Analysis
The typical curing profile of six parameters including the two heating rates (r 1 , r 2 ), two dwell times (dt 1 , dt 2 ), and two dwell temperatures (T 1 , T 2 ) were extracted from the curing profile as the design variables.
Sobol sensitivity analysis is widely used in various fields in the industry to identify the most critical input variables and improve the accuracy and robustness of models. It was employed to quantify the impact of input variables on the output of a model or simulation and provided a measure of the relative importance of each input variable by decomposing the total variance in the model output into contributions [30,31]. Herein, the Sobol sensitivity analysis method was employed to investigate the effects of curing profile parameters on the results. The contribution of each parameter is derived based on variance decomposition. The system input-output function f(x) can be decomposed into a summary of increasing dimensions: Polymers 2023, 15, 2437 5 of 12 in which f 0 is the constant term, which is equal to the expectation value of the output. Each input variable x i , i = 1, 2, · · · n is randomly distributed in the range of [0, 1]. The decomposed items in Equation (3) can be derived as the following functions: (8) in which dx ∼i is the integration of all variables. The total variance is defined as: The partial variances corresponding to the subset of parameters are defined as: Under the case that the input variables are mutually orthogonal, the variance decomposition can be derived as The Sobol sensitivity indices for a subset of parameters are defined as follows.

Multi-Objective Optimization
The composite laminates are typical an-isotropic materials. The heat transfer coefficients are different along the fiber direction and in-plane as well as along the out-of-plane vertical fiber direction. A large amount of cross-linking reaction heat cannot be eliminated in time. The temperature overshoot phenomenon would occur and result in an uneven temperature field inside the laminate and incomplete curing of the components. To optimize the curing profile of thick composite laminates, the first and second dwell temperatures T 1 and T 2 , the first and second dwell times dt 1 and dt 2 , and the first and second heating rates r 1 and r 2 were used as designed variables. To minimize the T max and the curing profile time, t cycle was the optimization object. Meanwhile, to ensure that the overall curing of the composite material components was complete after the curing stage was completed, the minimum value of DoC should be greater than 0.9: Find X = (dt 1 , dt 2 , T 1 , T 2 , r 1 , r 2 ) Min T max , t cycle S.T. DoC ≥ 0.9 (13) The steps of the multi-objective optimization design of the curing process were shown in Figure 2: Firstly, the Latin hypercube technique was used to generate random samples in the design space. Then the sample was re-organized to generate Python scripts for an efficient finite element simulation. The output values (T max , t cycle , DoC) from the FE model were saved for the surrogate model. Thirdly, the RBF was established as the surrogate model due to its applicability for higher-order nonlinear and multi-variable problems. Finally, NSGA-II multi-objective optimization was conducted. The radial basis function (RBF) maps the inputs to an output value based on its distance from a center point, which can be used to approximate complex functions and provide a good balance between model complexity and accuracy. Non-dominated Sorting Genetic Algorithm-II (NSGA-II), an extension of the original NSGA algorithm, is based on the concept of non-dominated sorting and is used to rank the solutions according to their dominance relationships. Polymers 2023, 15, x 6 of 12

Validation of the Finite Element Model
The mid-point temperature profile and DoC profile from numerical simulation and experiments in Ref. [25] were compared in Figure 3, and the values were listed in Table 2. It can be seen that the simulation results agreed well with the experimental results. The values of dwell temperature T1 from the temperature profile, finite element simulation, and experiment were 110 °C, 110 °C, and 104 °C. The relative error between the experimental result and the numerical result was 5.8%. The values of dwell temperature T2 from the temperature profile, finite element simulation, and experiment were 180 °C, 181 °C, and 180 °C. The relative error between experimental and numerical results was 0.6%. The maximum DoC for experimental and numerical results was 0.93 and 0.91. The relative error was 2.2%. It can be concluded that the established finite element model was accurate and reliable for the following parameter sensitivity analysis and multi-objective optimization.

Validation of the Finite Element Model
The mid-point temperature profile and DoC profile from numerical simulation and experiments in Ref. [25] were compared in Figure 3, and the values were listed in Table 2. It can be seen that the simulation results agreed well with the experimental results. The values of dwell temperature T 1 from the temperature profile, finite element simulation, and experiment were 110 • C, 110 • C, and 104 • C. The relative error between the experimental result and the numerical result was 5.8%. The values of dwell temperature T 2 from the temperature profile, finite element simulation, and experiment were 180 • C, 181 • C, and 180 • C. The relative error between experimental and numerical results was 0.6%. The maximum DoC for experimental and numerical results was 0.93 and 0.91. The relative error was 2.2%. It can be concluded that the established finite element model was accurate and reliable for the following parameter sensitivity analysis and multi-objective optimization.

Validation of the Finite Element Model
The mid-point temperature profile and DoC profile from numerical simulation and experiments in Ref. [25] were compared in Figure 3, and the values were listed in Table 2. It can be seen that the simulation results agreed well with the experimental results. The values of dwell temperature T1 from the temperature profile, finite element simulation, and experiment were 110 °C, 110 °C, and 104 °C. The relative error between the experimental result and the numerical result was 5.8%. The values of dwell temperature T2 from the temperature profile, finite element simulation, and experiment were 180 °C, 181 °C, and 180 °C. The relative error between experimental and numerical results was 0.6%. The maximum DoC for experimental and numerical results was 0.93 and 0.91. The relative error was 2.2%. It can be concluded that the established finite element model was accurate and reliable for the following parameter sensitivity analysis and multi-objective optimization.

Parameter Sensitivity Analysis
The 4 mm, 6 mm, 8 mm, 10 mm, and 12 mm thick laminates with [90 n /0 n ] s stacking sequence were used to evaluate the effect of thickness on the T max , ∆T (temperature gradient at the mid-point and surface-point), and DoC after the curing process. The cross-section temperature distributions of laminates with different thicknesses are shown in Figure 4. It can be seen that the maximum temperature always occurred in the mid-point regardless of the laminate thickness.

Parameter Sensitivity Analysis
The 4 mm, 6 mm, 8 mm, 10 mm, and 12 mm thick laminates with [90n/0n]s stacking sequence were used to evaluate the effect of thickness on the Tmax, ΔT (temperature gradient at the mid-point and surface-point), and DoC after the curing process. The cross-section temperature distributions of laminates with different thicknesses are shown in Figure  4. It can be seen that the maximum temperature always occurred in the mid-point regardless of the laminate thickness.   The cross-ply, uni-direction, quasi-static, and angle ply is the most commonly used stacking sequence in the industry, and a helicoidal bio-inspired stacking sequence with excellent impact resistance has drawn more attention in recent days [32][33][34][35]

Parameter Sensitivity Analysis
The 4 mm, 6 mm, 8 mm, 10 mm, and 12 mm thick laminates with [90n/0n]s stacking sequence were used to evaluate the effect of thickness on the Tmax, ΔT (temperature gradient at the mid-point and surface-point), and DoC after the curing process. The cross-section temperature distributions of laminates with different thicknesses are shown in Figure  4. It can be seen that the maximum temperature always occurred in the mid-point regardless of the laminate thickness.   The cross-ply, uni-direction, quasi-static, and angle ply is the most commonly used stacking sequence in the industry, and a helicoidal bio-inspired stacking sequence with excellent impact resistance has drawn more attention in recent days [32][33][34][35]  The cross-ply, uni-direction, quasi-static, and angle ply is the most commonly used stacking sequence in the industry, and a helicoidal bio-inspired stacking sequence with excellent impact resistance has drawn more attention in recent days [32][33][34][35]. Thus, the laminates with [90 8  in Figure 6a showed that the T max , ∆T, and DoC were affected slightly by the laminate stacking sequence. listed in Table 3. The results are shown in Figure 6b. It was indicated that the mold material mainly affected the uniformity of the temperature field. The ΔT value of aluminum mold was the highest, followed by copper mold and invar steel mold.

Parameter Sensitivity Analysis
Parameter sensitivity analysis requires a large number of samples to support accurate variance analysis and provide reliable sensitivity information. Herein, RBF models were trained to serve as the surrogates of time-consuming FE simulation models. Figure 7 showed the scatter plot of finite element prediction and RBF surrogate model prediction. The red lines are the perfect fitting function. The more scatter points around the perfect fitting function, the higher the accuracy of the surrogate model. The coefficient of determination (R 2 ) was used to evaluate the accuracy of the surrogate model. When R 2 equals 1, it indicated that all the scatters were located on the perfect-fitting line and that the surrogate model has the same prediction value as the FE model. As shown by the scatter plots in Figure 7, the surrogate models exhibit good consistency with the FE model in predicting Tmax, Doc, and tcycle since the scatters cluster around the perfect-fitting lines and the R 2 square value is above 0.9. These surrogate models will be used to support sensitivity analysis to identify the critical parameters that affect the simulation results of the curing process significantly.  Invar steel, aluminum mold, and cooper mold were used to investigate the mold material on the curing result. The thermal parameter comparison of the three materials is listed in Table 3. The results are shown in Figure 6b. It was indicated that the mold material mainly affected the uniformity of the temperature field. The ∆T value of aluminum mold was the highest, followed by copper mold and invar steel mold.

Parameter Sensitivity Analysis
Parameter sensitivity analysis requires a large number of samples to support accurate variance analysis and provide reliable sensitivity information. Herein, RBF models were trained to serve as the surrogates of time-consuming FE simulation models. Figure 7 showed the scatter plot of finite element prediction and RBF surrogate model prediction. The red lines are the perfect fitting function. The more scatter points around the perfect fitting function, the higher the accuracy of the surrogate model. The coefficient of determination (R 2 ) was used to evaluate the accuracy of the surrogate model. When R 2 equals 1, it indicated that all the scatters were located on the perfect-fitting line and that the surrogate model has the same prediction value as the FE model. As shown by the scatter plots in Figure 7, the surrogate models exhibit good consistency with the FE model in predicting T max , DoC, and t cycle since the scatters cluster around the perfect-fitting lines and the R 2 square value is above 0.9. These surrogate models will be used to support sensitivity analysis to identify the critical parameters that affect the simulation results of the curing process significantly. Figure 8 provides the sensitivity analysis of T max , DoC, and t cycle . It was indicated that T max and t cycle were affected by the six parameters simultaneously as all the correlation factors were all higher than 0.29; meanwhile, the dwell temperature T 2 was the most significant process parameter, followed by dwell time dt 1 and temperature increasing ratio r 2 . DoC was mainly affected by dwell time dt 1 , temperature increasing ratio r 2 , and dwell temperature T 1 , and the influences of temperature increasing ratio r 1 and dwell temperature T 2 were not obvious as their correlation factors were less than 0.1. rogate model has the same prediction value as the FE model. As shown by the scatter plots in Figure 7, the surrogate models exhibit good consistency with the FE model in predicting Tmax, Doc, and tcycle since the scatters cluster around the perfect-fitting lines and the R 2 square value is above 0.9. These surrogate models will be used to support sensitivity analysis to identify the critical parameters that affect the simulation results of the curing process significantly.   Figure 8 provides the sensitivity analysis of Tmax, DoC, and tcycle. It was indicated that Tmax and tcycle were affected by the six parameters simultaneously as all the correlation factors were all higher than 0.29; meanwhile, the dwell temperature T2 was the most significant process parameter, followed by dwell time dt1 and temperature increasing ratio r2. DoC was mainly affected by dwell time dt1, temperature increasing ratio r2, and dwell temperature T1, and the influences of temperature increasing ratio r1 and dwell temperature T2 were not obvious as their correlation factors were less than 0.1.

Multi-Objective Optimization
The optimized Pareto solution was shown in Figure 9. The two objectives, Tmax and tcycle, were opposite to each other, which means that a decrease in tcycle will inevitably lead to an increase in Tmax. Although Pareto frontiers can provide designers with many optimization solutions, designers can obtain desired solutions from different perspectives and obtain appropriate solutions. However, there are also designers who want to find an optimal and most satisfactory solution (Knee point) in the Pareto frontier, making each objective function as optimal as possible. In order to obtain the optimal solution from the Pareto front, the minimum distance method is used to select the Pareto front to obtain a Knee point. As shown in Figure 9a, the distance (Utopia point) of the ideal optimal solution from the Knee point is the smallest.

Multi-Objective Optimization
The optimized Pareto solution was shown in Figure 9. The two objectives, T max and t cycle , were opposite to each other, which means that a decrease in t cycle will inevitably lead to an increase in T max . Although Pareto frontiers can provide designers with many optimization solutions, designers can obtain desired solutions from different perspectives and obtain appropriate solutions. However, there are also designers who want to find an optimal and most satisfactory solution (Knee point) in the Pareto frontier, making each objective function as optimal as possible. In order to obtain the optimal solution from the Pareto front, the minimum distance method is used to select the Pareto front to obtain a Knee point. As shown in Figure 9a, the distance (Utopia point) of the ideal optimal solution from the Knee point is the smallest.
Polymers 2023, 15, x 9 of 12 Figure 8 provides the sensitivity analysis of Tmax, DoC, and tcycle. It was indicated that Tmax and tcycle were affected by the six parameters simultaneously as all the correlation factors were all higher than 0.29; meanwhile, the dwell temperature T2 was the most significant process parameter, followed by dwell time dt1 and temperature increasing ratio r2. DoC was mainly affected by dwell time dt1, temperature increasing ratio r2, and dwell temperature T1, and the influences of temperature increasing ratio r1 and dwell temperature T2 were not obvious as their correlation factors were less than 0.1.

Multi-Objective Optimization
The optimized Pareto solution was shown in Figure 9. The two objectives, Tmax and tcycle, were opposite to each other, which means that a decrease in tcycle will inevitably lead to an increase in Tmax. Although Pareto frontiers can provide designers with many optimization solutions, designers can obtain desired solutions from different perspectives and obtain appropriate solutions. However, there are also designers who want to find an optimal and most satisfactory solution (Knee point) in the Pareto frontier, making each objective function as optimal as possible. In order to obtain the optimal solution from the Pareto front, the minimum distance method is used to select the Pareto front to obtain a Knee point. As shown in Figure 9a, the distance (Utopia point) of the ideal optimal solution from the Knee point is the smallest. The Tmax, tcycle, and DoC from RBF were 178.9 °C, 0.91, and 17,667 s. Compared with the corresponding optimized values of 177.1 °C, 0.91, and 16,862 s from finite element simulation, the relative errors were 1%, 0, and 4.8%, which can be acceptable and prove the reliability of the surrogate model. The histories of temperature and DoC at the mid-