Wenyuan Tang1, Yingjie Xu1,2, Xinyu Hui1, Wenchang Zhang1. 1. State IJR Center of Aerospace Design and Additive Manufacturing, Northwestern Polytechnical University, Xi'an 710072, China. 2. Shaanxi Engineering Laboratory of Aerospace Structure Design and Application, Northwestern Polytechnical University, Xi'an 710072, China.
Abstract
The contribution of this work is introducing a multi-objective optimization method based on finite element (FE) numerical simulation to simultaneously control the curing time and cure-induced defects of C-shaped composites during a curing cycle. Thermochemical and thermomechanical coupled analysis is performed and validated experimentally to understand the evolution details of temperature, degree of cure and curing deformation. Aiming to achieve the simultaneous control of manufacturing cost and composite quality, the curing profile is optimized by employing the critical factors including the total curing time, the maximum degree of cure difference, and the curing deformation. The optimization result shows that the designed curing profile can effectively reduce the curing time and guarantee the curing quality. The total curing time of the optimization is reduced by 19%. The verification experiment is also conducted to prove the accuracy and effectiveness of the proposed optimization method.
The contribution of this work is introducing a multi-objective optimization method based on finite element (FE) numerical simulation to simultaneously control the curing time and cure-induced defects of C-shaped composites during a curing cycle. Thermochemical and thermomechanical coupled analysis is performed and validated experimentally to understand the evolution details of temperature, degree of cure and curing deformation. Aiming to achieve the simultaneous control of manufacturing cost and composite quality, the curing profile is optimized by employing the critical factors including the total curing time, the maximum degree of cure difference, and the curing deformation. The optimization result shows that the designed curing profile can effectively reduce the curing time and guarantee the curing quality. The total curing time of the optimization is reduced by 19%. The verification experiment is also conducted to prove the accuracy and effectiveness of the proposed optimization method.
In recent years, carbon-fiber-reinforced thermoset polymer (CFRP) composite structures have been widely used in aeronautics [1] and aerospace [2], automotive [3], marine [4] and civil engineering [5] owing their high specific stiffness and strength. Autoclave solidification process is one of the most used methods to manufacture the CFRP composites, particularly the primary load-bearing structures in aircraft [6,7,8,9]. In an autoclave process, composites undergo a crosslinking reaction under the prescribed temperature and pressure. The quality of the composites strongly depends on the temperature and pressure curing profiles of the manufacturing process. Unreasonable curing profiles can result in some defects, such as temperature overshoots and non-uniform degree of cure within the composite parts, even the curing deformation [10]. Most of the industry manufacturers adopt the relative conservative curing profiles with long curing time to produce composites. Conventionally, the long curing time negatively affects the manufacturing productivity and leads to a very high energy cost.Therefore, it is necessary to achieve a more appropriate curing profile which reduces both the process induced defects and curing duration and then improves the high-performance of the composites with low costs. Better understanding of the physical state of composite during crosslinking reaction is essential for developing an optimization method of curing profile. Several experimental contributions have listed. Hui et al. [11] investigated the composite state evolution as temperature, pressure and viscosity during the composite solidification. Almazan-Lazaro et al. [12] observed the impregnation velocity of resin in VARI process by computer vision. Kahali Moghaddam et al. [13] studied the influence of pressure on curing process by embedding different types of pressure sensors in composite. Sorrentino et al. [14] monitored composite temperature evolution using the thermocouples inserted in prepreg sheets. Silva et al. [15] detected the temperature and fluidity of the polymer under compound extrusion manufacturing processes. However, experimental approach is inefficient and the sensor introduced would affect the curing process. Consequently, numerical simulation of the curing process is a relatively efficient and economical approach. Process modeling of polymerizing of thermoset polymer composites have improved from simple one-dimensional finite-difference models [16] to advanced three-dimensional finite element models [17]. Compared with the experimental method in Liu et al. [18], simulation approach could provide comprehensive information during the cure cycle more efficient, such as the trend of temperature [19], degree of cure [20], residual stress of cured part [21] and so on.Thus, the optimizations of curing profiles based on numerical simulation models of the curing process have become increasingly popular in the literature. For instance, Struzziero and Skordos [22] introduced a multi-objective optimization procedure that combined a finite element solution of process modeling with a genetic algorithm to mitigate the temperature overshoot and cure duration. Axenov et al. [23] proposed a Pareto-based multi-objective optimization method based on curing process simulation to capture optimal two-dwell curing profile for shell-like composite cure quality controlling in autoclave. Rubino et al. [24] studied a fuzzy logic and artificial neural network based computing method to predict the evolution of degree of cure and temperature during curing in autoclave, and optimize the cure profile of the composite. Oh and Lee [25] studied the curing of glass-epoxy composite laminates, and found the optimal value of the curing profile parameter to minimize the overshoot of the laminate. Pantelelis et al. [26] combined genetic algorithm and one-dimensional process model to optimize the optimal curing profile parameters of thermoset matrix composites, thereby shortening the curing duration under the constraints of maximum final degree of cure. Zhang et al. [27] developed a four dwell curing profile for solving the overshoot of degree of cure and temperature during curing process for thick CFRP composite laminates.However, the present existing work rarely introduces the cure-induced defects into the optimization problem. Motivated by this idea, this paper develops a methodology integrating numerical simulation of curing process with particle swarm optimization algorithm to optimize the curing profile for the C-shaped CFRP composite specimens. The simultaneous control of curing time and process-induced defects is achieved by minimizing the curing time with constraints of the typical defects including non-uniform degree of cure, insufficient cure and curing deformation. Finally, two groups of C-shaped CFRP specimens are produced by the optimal curing profile and the profile recommended by manufacturer. The test and comparison of the parts illustrate well the capability of the developed optimization methodology.
2. Materials and Methods
2.1. Curing Experiment of C-Shaped Specimens
The C-shaped AS4/3501-6 composite specimens were manufactured under the recommended curing profile (see in Figure 1a) and the deformation after curing process was measured to evaluate the composite quality. Three groups of specimens with layup , and were designed, where 0° direction was the length direction. The size of the each C-shaped specimen was 230 mm in length, 200 mm in web depth and 40 mm in flange width, as illustrated in Figure 1b. To quantify the deformation caused by curing process, the spring-back angle was introduced and recorded by the HandySCAN 300 laser scanner, as shown in Figure 2. The average values of the measured specimens with different layups were summarized in Table 1.
Figure 1
Experimental design of the C-shaped composite specimen: (a) Curing profile; (b) Geometric size.
Figure 2
The predicted spring-back angel of the composites with layup [0/90]4S and schematic of laser scanning experiment.
Table 1
The measurement value of the spring-back angle.
Layup
Experiment (∘)
0/904S
1.34
0/45/−45/902S
1.36
45/−454S
1.48
2.2. Numerical Simulation of Curing Process and Verification
To improve the accuracy of curing process simulation, many physics-chemical based views of composite curing involve numerous interacting systems and a variety of process-dependent parameters that must be accounted for. The thermo-chemical model is established to predict the evolution of degree of cure and temperature of the composites. Then the obtained temperature and degree of cure fields are employed into the thermo-mechanical model. The detailed modeling framework is presented as follow.
2.2.1. Thermo-Chemical Coupled Simulation
The evolution of the degree of cure and temperature during the C-shaped specimen curing process is calculated by a thermos-chemical coupled heat transfer simulation. The thermo-chemical coupled heat transfer equation is expressed as [28]:
where and are the density and specific heat of the composite, respectively. T is the current temperature and t is the curing reaction time. represent the thermal conductivities of the composite. is the resin density, is the fiber volume fraction, and is the total heat energy released from curing reaction of a unit mass of the resin. is the instantaneous curing rate of the resin and can be determined by the cure kinetics model.Considering the structural symmetry of the C-shaped composites, 1/2 model is constructed to solve the heat transfer analysis using finite element package ABAQUS. The finite element model meshed with three-dimensional 8-node linear heat transfer brick (DC3D8) elements. The temperature boundary conditions of the model are illustrated Figure 3, where is the mold temperature, and are the air temperature and pressure in the autoclave, respectively. The heating air is used to heat the composites, and the film condition is loaded on the composites surface. The curing profile shown in Figure 2 is imposed on the inner surface of the composites contacting with the mold.
Figure 3
Temperature boundary conditions.
The cure kinetics model of the resin, specific heat capacity and thermal conductivities of the composite are all implemented through the user subroutines UMATHT and USEFLD. The cure kinetics model for the 3501-6 resin is expressed as [29]:
where represent the curing rate constants and can be further defined by the Arrhenius equations:
where are the pre-exponential factors. are the reaction activation energies and R is the universal gas constant.The specific heat capacity of the composite C is calculated using the rule of mixture [30,31]:
where is the weight fraction of the fiber. and denote the specific heat capacity of fiber and resin, respectively.The longitudinal () and transverse direction ( and ) thermal conductivities of the composite are computed as [32]:
where and are the longitudinal and transverse thermal conductivities of the fiber. is the thermal conductivity of the resin. Considering the model in this study is macro scale and the size of the C-shaped specimens is small, the influence of Kapitza resistance [33,34] on the heat transfer of air-composite interface (and also at the fiber-resin interface) is ignored in the simulation process.The parameters values for Equations (1)–(7) are reported in Table 2. Note that the specific heat capacity and thermal conductivity of AS4 fiber has a linear dependence on temperature, while the specific heat capacity and thermal conductivity of the resin depends on both the degree of cure and temperature.
Table 2
Parameters values for the cure kinetics, specific heat capacity and thermal conductivities of resin [35].
Material
Property
Value
3501-6 resin
pr(kg/m3)
1272
Cr(J/(kg·K))
4184×(0.468+5.975×10−4×T−0.141a)
Kr(W/(K·m))
0.04184×(3.85+(0.035×T−0.141)a)
Hr(J/kg)
473,600
R(J/(mol·K))
8.3143
A1(min−1)
2.101×109
A2(min−1)
−2.014×109
A3(min−1)
1.96×105
ΔE1(J/mol)
8.07×104
ΔE2(J/mol)
7.78×104
ΔE3(J/mol)
5.66×104
AS4 fiber
pr(kg/m3)
1790
Cf(J/(kg·K))
750+2.05×T
Kfx(W/(K·m))
7.69+0.0156×T
Kfy(W/(K·m))
2.4+0.00507×T
2.2.2. Thermo-Mechanical Coupled Simulation
The calculated degree of cure and temperature are further used to predict the residual stress and strain in the cured composites. It has been demonstrated that the thermoset resin composite exhibits viscoelastic behavior during curing process [36]. The strain tensor of the composite in the curing process can be described as [37]:
where and represent the viscoelastic and thermal expansion strain, which are determined by the viscoelastic and thermal expansion properties of the composite. is the chemical shrinkage strain of the composite and can be calculated by
where and represent the longitudinal and transverse chemical shrinkage strains of the composite, respectively. and are the longitudinal and transverse chemical shrinkage coefficients of the composite.In this study, a linear viscoelastic constitutive equation [38] is used to model the composite:
where represents the stress tensor of composite. t and denote the current time and past time, respectively. is the stiffness matrix of composite and can be modeled by generalized Maxwell equation. The stiffness modeled by a parallel connection of n spring-dashpot elements, as shown in Figure 4, can be further approximated by the Prony series as:
where is the stiffness matrix of the fully relaxed (rubbery) composite and is the stiffness matrix of the unrelaxed (glassy) composite. and denote the weight factor and the relaxation time for the ith Maxwell element, respectively.
Figure 4
Illustration of the generalized Maxwell model.
It is assumed that the and of AS4/3501-6 composite are constants independent of temperature and degree of cure. The relation between and are described as:
where r is a constant between 0 and 1. For the AS4/3501-6 composite, the value of r is 0.14 [36]. The unrelaxed materials properties for AS4/3501-6 composite are listed in Table 3. The stiffness matrix coefficients of can be derived from the unrelaxed modulus and Poisson’s ratio.
Table 3
Unrelaxed properties of AS4/3501-6 composite [39].
Property
Value
Ex(GPa)
126
Ey(GPa)
8.3
Gxy(GPa)
4.1
Gyz(GPa)
2.8
vxy
0.3
vyz
0.5
βx(K−1)
0.5
βy(K−1)
35.3
δx
−0.017%
δy
−0.88%
The relaxation time is the function of the degree of cure and can be expressed as:
where is the relaxation time under reference degree of cure. and are formulated as [40]:The weight factor and the relaxation time at reference degree of cure for AS4/3501-6 are listed in Table 4.
Table 4
Relaxation times and weight factors at the reference degree of cure ( = 0.98) [40].
i
τi (min)
wi
i=1
2.922×101
0.0569
i=2
2.922×103
0.066
i=3
1.82×105
0.083
i=4
1.10×107
0.112
i=5
2.83×108
0.154
i=6
7.94×109
0.262
i=7
1.95×1011
0.184
i=8
3.32×1012
0.049
i=9
4.92×1014
0.025
The thermo-mechanical coupled simulation is conducted within ABAQUS. User subroutines UMAT and UEXPAN are used to model the constitutive behavior of the composites. A static implicit analysis is used to calculate the stresses and strains during the curing process. Then, another analysis step is conducted to calculate the curing deformation induced by the curing residual stresses. The displacement boundary condition in the static analysis is shown in Figure 5, where the axial movement of x and y and the rotation of the z axis are constrained.
Figure 5
The established finite element model.
2.2.3. Verification for Simulation
By using the thermal-chemical and thermal-mechanical analysis, the cure-induced deformations of the C-shaped composites are obtained. Figure 6 presents the deformation of the composite with layup . The predicted spring-back angle of the composite can be defined by:
where is the maximum displacement of the composite after curing process and L is the length of the composite flange.
Figure 6
The predicted curing deformation of the composites with layup .
To validate the accuracy of curing simulation, the calculated spring-back angle of the model is compared with the experiments, as summarized in Figure 7. The comparison shows good agreement with the composites of different layups, which highlights the prediction capacity of the proposed models.
Figure 7
The predicted curing deformation of the composites with layup .
In order to further evaluate the overall curing quality of the composite, the evolutions of the minimum degree of cure of the whole model are extracted. The minimum value exceeds 0.95 at the end of the curing process, which indicates that the composite part has completed cure reaction. In order to avoid residual stress caused by degree of cure distributed unevenly in the composites, the simultaneous difference of degree of cure is estimated from the center and bottom of the model, as shown in Figure 8. The first peak of appears after the first heating stage, since the first degree of cure rises rapidly. Then keep holding for a while and tends to reduce owing the thin-walled composites are easy to achieve isothermal. For the second heating stage and the early stage of second dwell stage, the curing reaction becomes severe, which induces the temperature difference between the middle and outer regions of the composite part. A transient relief appears when the degree of cure of surface is exceeded by center. With the end of the curing reaction, is reduced gradually and finally remains stable. The global maximum difference of degree of cure at the end of the curing process is 0.01.
Figure 8
Evolution history of maximum difference of degree of cure.
2.3. Optimization of Curing Profile
2.3.1. Optimization Problem Formulation
Based on the curing simulation result, an optimization problem is established for simultaneous control of curing time and process-induced quality of the composites. C-shaped AS4/3501-6 with layup is used as the optimization model. Two-dwell curing profile as shown in Figure 9 is optimized to shorten the curing time and the typical defects (non-uniform degree of cure and spring-back angel) with maintaining sufficient cure level of the composites. Six parameters including the two heating rate (, ), two dwell time (, ), and two dwell temperature (, ) are extracted from the curing profile as the design variables.
Figure 9
Two-dwell curing profile.
The optimization problem can be formulated mathematically as:
where is the total curing time. is the maximum difference of degree of cure during the curing process. is the minimum degree of cure in the composite specimen after cure and is spring-back angel of the composites after cure. To solve the optimization problem, the multi-objective particle swarm optimization (MOPSO) algorithm is introduced to yield a global optimum from the cure simulation model.
2.3.2. Particle Swarm Optimization Algorithm
The MOPSO algorithm can handle both linear and nonlinear optimization problems [41]. In MOPSO algorithm, every particle is designed to simulate the behavior of birds that fly in the design space. Their velocity and position not only could be affected by the best position that the particles flied, but also the position saved in archive. Under the effect of particle swarm and archive swarm, the MOPSO particle swarm algorithm retains diversity, besides demonstrates its search capabilities. According the good position, the particles update their position and fly velocity based on update formula as follows:
here and represent position coordinate and velocity coordinate of the ith individual particle, respectively, (the subscripts k and refer to the current and the next iterations, respectively). The first part of the formula used for updating the velocity describes the influence of the previous velocity on the current velocity. Inertia weight as proxies to quantification the influence, and is set to 0.875 in this study. The second part of the formula relies on the distance between the recent position and the best position of the particle that has flown so far. is the best previous position of the ith particle. is stochastic acceleration terms as proxies to the acceleration factor that attracts each particle towards the good position in its own experience. The third part depends on the distance between the recent position and the chosen particle from the archive. is a particle position taken from the archive. The index j is selected in the following way: After collecting the new position respecting the non-dominated solution, the adaptive gird the particle points in the archive can be wrapped by the hypercube, and a fitness value is assigned to estimate the density of the hypercube, which is inversely proportional to congestion. If the archive is full, the hypercube with the lowest fitness value will delete one particle to release storage for the new particle. After estimation, roulette-wheel selection is used to choose the hypercube based on the fitness value, and j is randomly selected within such hypercube. The size of archive is 20. is stochastic acceleration terms as proxies to the acceleration factor which attract each particles towards the chosen particle from the archive. In this work, = 2 and = 2 are set. and are two random constants from [0,1].The problem-specific constraints and the layout of the variable are included in most optimization problems. For the optimization in this study, the problem-specific constraints are the minimum final degree of cure. The layout of the variable is design bounds of the curing profile parameters, where the allowable ranges are listed in Table 5. The solution would be ignored, if the particle searches out of the layout of the variable, even the problem-specific constraint is satisfied. So making sure none of the particles reach the outside layout of the variable and verifying whether met the problem-specific constraint is crucial. In this paper, the particles reach outside range is handled by the harmony search algorithm-based method. The detailed descriptions of the harmony search algorithm-based method can be found in previous studies [42]. The MOPSO algorithm is integrated with the curing simulation model and the Pseudo code of optimization is shown in Algorithm 1. For the present optimization problem, a population of 20 particles is used, the maximum number of generations is set as 50, and the maximum velocity is 5.
Table 5
Range of the design variable.
Design Variable
Parameter
Range
Minimum Increment
First heating rate
a1 (K/min)
[0.5, 5.0]
0.1
Second heating rate
a2 (K/min)
[0.5, 5.0]
0.1
First dwell temperature
T1 (K)
[388, 428]
0.1
Second dwell temperature
T2 (K)
[448, 488]
0.1
First dwell time
t1 (s)
[600, 6000]
60
Second dwell time
t2 (s)
[600, 9000]
60
Cooling rate
a3 (K/min)
[0.5, 3.0]
0.1
3. Results and Discussion
The optimization process is executed using a twelve-core processor (2.8 GHz) and 50 GB memory. The optimization converges after 700 iterations and takes about 50 h. The multi-objective optimization results are illustrated in Figure 10 and Figure 11. As can be seen, the 3D contour plot of the optimal Pareto front suggests the three objectives ( difference of degree of cure, spring-back angel and total curing time) have the following relationship: If the differences of degree of cure and the spring-back angel are prioritized, the curing time will be prolonged. If focusing on the total curing time, differences of degree of cure and spring-back angel will be unsatisfactory. In order to achieve simultaneous control of curing time and process-induced defects, it is necessary to select the most reasonable curing profile among the number of Pareto fronts according to balancing the weight of total curing time, the maximum differences of degree of cure and spring-back angle. Therefore under the premise of ensuring the curing quality, the result of minimizing the curing time is selected.
Figure 10
Three-dimensional scatter plot of the Pareto front of the multi-objective optimization result.
Figure 11
Surface plot of the Pareto front of the multi-objective optimization result.
The optimized curing profile selected from Pareto front are reported in Table 6. For comparison, the parameters of the initial curing profile are also included. It is found that the optimized curing profile reduces to 14,686 s by about 19%. Figure 12 presents the optimized curing profile. It is observed that the first heating rate decreases from 2.5 K/min to 2.0 K/min and the first dwell temperature rises to 420.1 K. Properly increasing the dwell temperature helps to advance the curing reaction, and reducing the heating rate can effectively avoid the aggravation of the temperature difference caused by the curing exothermic. Besides, it is considered that the resin exhibits a viscous or highly elastic state in the early stage when the degree of cure is low. The premature peak of degree of cure difference weakens the cure-induced residual stress and further deformation. The first dwell time decreases from 3600 s to 720 s. To improve the cure efficiency, the first dwell time should not be long, because the temperature at this time is not sufficient to completely cure the composites. After the first peak, there has been an equilibrium period and the optimal curing temperature has not yet been reached. If it enters the second heating stage at this time, the curing cycle can be effectively shortened. After optimized, the second heating rate increase from 2.5 K/min to 4.0 K/min. In the second heating stage, higher heating rate push the temperature to reach the required temperature faster, which can accelerate curing reactive at heating stage. A higher heating rate causes the surface temperature of the composites to rise rapidly, which thereby reduces the temperature overshoot caused by internal curing exothermic and keeps the temperature and curing rate of the center and surface of the composites at a relatively same level. The second dwell temperature increase from 450 K to 475.2 K and the dwell time rises from 3600 s to 5760 s. The higher dwell temperature promotes the curing rate and the longer dwell time makes the reaction complete in the second heating stage to avoid large temperature and degree of cure difference. Although the spring-back deformation of the part can be reduced by reducing the cooling rate, the change from 2.5 K/min to 2.9 K/min has little effect on the deformation but can effectively reduce the curing time.
Table 6
Parameters of initial cure profile and optimized curing profile.
Parameters
Initial
Optimized
a1 (K/min)
2.5
2.0
a2 (K/min)
2.5
4.0
T1 (K)
389
420.1
T2 (K)
450
475.2
t1 (s)
3600
720
t2 (s)
3600
5760
a3 (K/min)
2.5
2.9
Sa (∘)
1.218
1.256
ttotal (s)
18,144
14,686
af
1
1
ad
0.01038
0.01006
Figure 12
Maximum difference of degree of cure evolution.
The self-organizing map (SOM) can provide visualization between design parameters and optimization objectives and inspire designers how design variables affect optimization objectives [43]. Figure 13 describes the relationships between the design variables and objectives with respect to the curing process using SOM. The label of each map represented by color indicates the weight value of the design variable of the optimization network, where the blue denotes the lower limit and the yellow denotes the higher limit. The relationship between the design parameters and objectives can be found by comparing the distribution law of the color map in a visualized approach.
Figure 13
Maximum difference of degree of cure evolution.
It is observed in Figure 13 the distribution of the first heating rate is similar to degree of cure difference , but opposite to total curing time , which suggests a higher make the increase and leads to a shorter . Based on the color map, the influences of other variables can be derived by the same way. For instance, the can be reduced by reducing the heating rate of the two stages and , increasing the dwell time of the first stage , or reducing the dwell temperature of the two stages and .
Experimental Verification
In order to verify the proposed multi-objective optimization method, the initial and optimized curing profiles shown in Table 6 are both implemented for the AS4/3501-6 C-shaped composites using autoclave system. The geometric size and layup property of the composites are remain consistent. Figure 14 illustrates the prepregs and the mold before entering the autoclave. After curing process, the average values of the spring-back angle are measured and compared with the numerical results, as shown in Table 7.
Figure 14
Prepregs and mold.
Table 7
Comparisons of measurement and simulation results of the spring-back angle.
Measurement (∘)
Simulation Result (∘)
Initial
1.218
1.225
Optimized
1.21
1.18
As far as the manufacturing cost and efficiency is concerned, the curing time of the optimized profile is shortened by 19% compared with the initial profile, which indicates the feasibility of the proposed optimization method.
4. Conclusions
A multi-objective optimization method is proposed to design the two-dwell curing profile of the CFRP composites to simultaneously control curing time and cure-induced defects. Through the thermo-chemical and thermo-mechanical coupled analysis, the evolution history of uneven cure and spring-back deformation is firstly evaluated and validated by the experiments. Then the curing profile which controls the composite manufacturing is optimized to improve the cure quality and costs through the multi-objective particle swarm algorithm. The verification experiment is also conducted to prove the accuracy and effectiveness of the optimization method. The result indicates that using the optimized curing profile achieves a better standard including high degree of cure, short curing time, and reliable quality controlled. Compared with the initial profile, the total curing time is shortened to 14,686 s, which is reduced by 19%. The self-organizing map further visualizes the relationship between the design parameters of the curing profile.
Figure 1
Figure 2
Figure 3
Figure 4
Figure 5
Figure 6
Figure 7
Figure 8
Figure 9
Figure 10
Figure 11
Figure 12
Figure 13
Figure 14