Residual stresses in thick composites considering through-thickness temperature gradients during manufacturing

In the composite materials of wind turbine rotor blades, residual stresses form during manufacturing due to the different thermo-chemo-mechanical behaviors of the fibers and the matrix. These residual stresses decrease the ultimate and fatigue strength of the blade because the failure stress and cycles to failure are reached at lower levels of external loading. It is thus of utmost importance to consider residual stresses in the design stage to obtain reliable blade designs. This work deals with the prediction of residual stresses in a thick composite representative for a spar cap in a wind turbine rotor blade. A model is formulated for the thermo-chemo-mechanically coupled problem. An emphasis is given on the non-uniform and time-variable temperature distribution across the thickness of the composite. The different sources of residual stresses (e.g., chemical shrinkage, non-uniform distribution of different thermal expansion coefficients of fibers and matrix, exotherm, etc.) are analyzed, quantified, and compared. The impact of the time-dependent temperature distribution on the formation of residual stresses is highlighted by means of a representative numerical example.


Introduction
In the composite materials of wind turbine rotor blades, residual stresses are a phenomenon that is relevant to consider.Residual stresses not only decrease the ultimate and fatigue strength of the rotor blade [1][2][3], but can also lead or at least contribute to the formation of unwanted fiber waviness or undulations [4][5][6].However, the calculation of residual stresses is subject to intensive research.Some authors used a viscoelastic model with a group of Maxwell elements [7][8][9][10].Others assumed a linear relation between the stiffness of the matrix and the curing degree, and calculated the residual stress from the curing process [11][12][13].
Residual stresses are generated during the manufacturing process, and if it is to be analyzed, the manufacturing process must be considered and a numerical model needs to be built in order to simulate the composite behavior during the curing process [14].In this paper, the following mechanisms are considered: First, there is chemical shrinkage in the matrix due to the network formation on curing.Then, there is a non-uniform distribution of thermal expansion coefficients due to different thermo-mechanical characteristics of the fibers and the matrix.There is also a non-uniform and time-dependent temperature distribution, especially in thick laminates.This is caused by the exothermal curing reaction in the matrix and the imbalance of heat exchange 1293 (2023) 012038 IOP Publishing doi:10.1088/1757-899X/1293/1/012038 2 across the surfaces, i.e., there is heat transfer into the laminate on the tool side and heat transfer out of the composite on the opposite side.A thermo-chemo-mechanical strong coupling approach is employed to calculate the residual stresses in a thick UD glass/epoxy laminate that is representative for a spar cap of a wind turbine rotor blade.

Curing kinetics model
In the wind energy industry, resin transfer molding technology is commonly used for the manufacturing of rotor blades [15].Therein, dry fiber mats are placed in the mold and packed into a vacuum bag.Then a vacuum is drawn so that the resin, which was previously mixed with the hardener, flows in and wets the fiber mats.To decrease the viscosity and therefore increase the flowability of the resin and to initiate the curing process, the mold is heated, and the heat is transferred to the resin and fiber system.The temperature is slowly increased to accelerate curing and the temperature is fixed at the final curing temperature for a specified time period until the matrix is fully cured.After curing is finished, the heating system is switched off and the composite cools down to room temperature.To simulate this process, the chemical, thermal and mechanical are modeled in the following.

Chemical model
The chemical analysis covers the polymerization reaction in the matrix, i.e., the formation of molecular chain networks.These networks result in an increase of stiffness, i.e., the change from the liquid state to the solid state.This so-called curing process is measured by means of the curing degree denoted by β.The curing degree is a function of time defined by the expression Herein, Q(t) is the specific heat generation due to the exothermic reaction up to the current instance of time t and Q ex is the specific total exothermic heat generation.Hence, the curing degree is described by the ratio between the current amount of heat generation and the total heat generation during the exothermic curing reaction.It holds that 0 ≤ β ≤ 1.A value of β = 0 means that the curing reaction did not start yet while a value of β = 1 means that the curing reaction is already finished.The rate equation of the curing degree is defined by which uses a Kamal-Sourour reaction type, see [16,17].The coefficients κ 1 , κ 2 , m, and n can be calculated using the relationships with i = 1, 2 (3) , and C 4 are material-dependent parameters, R is the gas constant, and T is the absolute temperature.The parameter set used in this study are from [16] and can be seen in Table 1 in the appendix.

Thermal model
Heat is introduced by the tool heating system and the exothermic reaction.The heat conduction through the material requires a transient thermal analysis.Here, Fourier heat conduction analysis is used [18].In addition to the heating from the outside, the heat generated from the inside by the curing exotherm also needs to be considered.The governing equation for the transient thermal analysis is given by where i, j = 1, 2, 3 are the directions of the Cartesian coordinate system, x j is the position vector element of a point of interest in Cartesian coordinates, ρ is the mass density of the composite material, c p is the specific heat capacity, and k ij is the tensor element of thermal conductivity.

Mechanical constitutive equations
As the curing reaction proceeds, the stiffness of the epoxy resin increases.At the beginning of the curing process, i.e., β = 0, the stiffness of the matrix is denoted by E m0 .When the curing is completed, i.e., β = 1, the stiffness of the matrix is denoted by E m1 .The initial stiffness of the matrix is supposed to be zero (liquid resin), but as an exact zero is numerically critical due to a singularity issue, it is assumed to be sufficiently small, i.e., E m0 = εE m1 , with ε = 10 −3 .A linear relationship between the elastic moduli and the curing degree is further assumed in this paper.Hence, the stiffness of the matrix at an arbitrary instance of time t denoted by E m (t) is proportional to the curing degree β (which is a function of time) and it holds (5)

Multi-scale multi-physics model for residual stress calculation
For the multi-physical problem that we need to solve, the following calculation steps are employed.First, the transient thermal analysis is executed to obtain the through-thickness temperature distribution as a function of time.Then, the temperature distribution in each time step is fed into the chemical analysis to calculate the curing degree.The curing degree is then put into the thermo-mechanical and chemo-mechanical couplings to calculate the stresses in the fibers and the matrix.To reduce the computational cost, the concept of a Representative Volume Element (RVE) is followed here, see e.g.[19][20][21][22].It is assumed that the major processes related to heat exchange take place in through-thickness direction, and the heat flux in the laminate plane is considered negligible.Hence, we reduce the full-scale model of Figure 1a) to 1D RVE model in throughthickness direction for the thermal analysis, which is named RVE T , see Figure 1b).It is a mesoscopic model because the laminate is treated in a smeared manner, i.e., the fibers are not resolved.The RVE T model is evaluated in a discrete number of evaluation points (points 1-5 in Figure 1) that are distributed in an equidistant manner across the laminate thickness.After the through-thickness temperature distribution in each time step is used to calculate the curing degree and with it the matrix stiffness, a 1D RVE model in fiber direction is used to calculate the stresses in the fibers and the matrix, respectively, see Figure 1c).This microscopic model that explicitly resolves the fibers is named RVE σ .

Chemical residual stress
During the curing process, the epoxy resin undergoes a chemical shrinkage of approximately 5% [18].If a rigid connection between the fibers and the matrix was assumed, see Figure 2a), the resultant strain in longitudinal direction in the 1D model would constant.Consequently, the fibers would be heavily compressed.However, the fibers cannot withstand such large compression without micro-buckling.Hence, there must be a strong shear deformation in the close vicinity of the fibers.This area is called interphase layer in the following and has a very small thickness, see Figure 2b).The initial stiffness of the interphase layer is very small and increases as the curing process progresses.It is assumed that the stiffness of the interphase layer denoted by E i (t) is proportional to the curing degree, but smaller than the matrix stiffness, because the 3D network formation in the interphase layer is disturbed by the fiber surface.It is assumed that

Thermal residual stress
For quick manufacturing processes in composite manufacturing, the curing temperature needs to be relatively high.After curing, the temperature must be cooled down to room temperature.There is a big difference in coefficient of thermal expansion between the fibers and the matrix.Due to higher thermal shrinkage in the matrix during cool-down, compression strains and stresses form in the fibers while tensile strains and stresses form in the matrix.For the thermo-mechanical cool-down analysis, one needs to know the temperature distribution at the end of the curing process.At this instance of time, the resin will be cooled down to room temperature.The respective temperature distribution is an outcome of the transient thermo-chemo-mechanical curing analysis and is an initial value condition for the calculation of the thermal residual stress.
It will be seen in the simulation results that the temperature distribution at the beginning of the cool-down process is nonlinear in the through-thickness direction, resulting in nonlinear distributions of thermo-mechanically triggered residual stresses from cool-down.

Numerical example
The numerical example is a 5 cm thick uni-directional composite laminate, which is representative for a spar-cap section of a wind turbine blade.The material parameters used for the analysis are given in the appendix, see Table 1, Table 2, and Table 3.There is a heating source from the bottom (considered the tool surface), which gives a heat flux of 400 J/m²/s.After heating for 10 hours, the heat source will be removed (tool heating system switched off).
On the top of the laminate there is a heat exchange with the surrounding air and the heat transfer coefficient is 8 J/m²/K.

Results and discussion
The multi-field coupled simulation described above was implemented in Abaqus by means of a user subroutine.The simulations are carried out without considering the exothermal peak as a first step.The exothermal peak is then included to allow for a thorough evaluation of its relevance for residual stress calculation.

Curing kinetics
The first simulation was executed without consideration of exotherm.The results are shown in Figure 3 (temperature on the left, curing degree on the right).The temperature curves are very smooth in all points observed in thickness direction and have a parabolic shape.After heating for 10 hours and removing the heat source, the temperatures decrease.The curing degree plotted against time has a sigmoidal shape.This means that the curing reaction is slow at the beginning, then quickly accelerates (high slope after a short time) and slows down again shortly before the curing process is finished.It can also be observed that the higher the temperature, the quicker the curing process takes place, which is in line with practical experience in composite manufacturing.
As a next step, the exotherm is included in the simulation.The simulation results are shown in Figure 4.As the curing process proceeds, the temperature of the material increases gradually.
There is an additional increase after about 3 h associated with the exothermic peak.Because of the higher temperature, the curing process accelerates and is finished after approximately 9 h in point 5 (at the top surface).In this point, the curing process was not finished after more than 12 h when the exotherm was not considered, cf. Figure 3.This highlights the importance to include the exotherm, as the temperature distribution and curing degree are more accurate in this case, and with it the residual stress prediction as well.
It can be seen in Figure 4 that the curing process starts after approximately 2 h, although the temperature is already quite high earlier.This is because of the chosen material properties that were adopted from [16].In wind energy, a quicker curing process would be advantageous to decrease the tool occupation time and increase the blade output frequency.

Residual stresses
As the chemical shrinkage progresses, the epoxy resin compresses the fiber through the interphase layer.The resulting fiber stress is shown on the left of Figure 5.As the laminate is internally in equilibrium, the norm of the matrix stress is the same, but has the opposite sign.The increase of the fiber stress with time is consistent with the increase of the curing degree, which is shown on the right of Figure 4.    Due to the uneven heat exchange with the tool at the bottom and the surrounding air at the top, the temperature distribution is non-linear throughout the manufacturing process, not only during curing, but also during the cool-down procedure.Due to the lower temperature at the top, the time needed for curing increases from the bottom to the top, respectively.After curing is finished, the tool heating is switched off and the composite cools down.Due to different coefficients of thermal expansion in the fibers and the cured matrix (the matrix will experience larger thermal shrinkage than the fibers), the fibers experience significant compressive residual stresses on cool-down, see on the right of Figure 5.
The (compressive) residual stresses in the fibers are non-uniform across the thickness.This is a consequence of non-uniform through-thickness temperature distributions during manufacturing and different time periods for the completion of curing processes along the thickness.The through-thickness gradient in residual stress results in shear stresses between different layers (which are not considered yet), and, consequently, may trigger undulations and fiber wavinessor may at least amplify originally existing undulations and fiber waviness.A detailed evaluation of the coupling between residual stresses and formation of undulations is work in progress and will be investigated in future publications.

Conclusions
In this paper, a multi-physics and multi-scale approach was proposed for the calculation of residual stresses in thick composite laminates that form during the manufacturing processes.One essential part of the model is the inclusion of a transient thermal analysis revealing the required boundary conditions for the evaluation of the curing degree inside the composite structure, which in turn enters the mechanical constitutive equations.A fully chemo-thermomechanically coupled procedure was proposed comprising exothermal peak, chemical shrinkage, thermal shrinkage, and curing state-dependent constitutive equations for the matrix.
The material properties used in this paper were partially extracted from literature, and partially based on assumptions, so that physically meaningful results were obtained.Especially the thermal behavior is difficult to measure but will be calibrated against measurements recorded during the manufacturing of a spar-cap of a real wind turbine blade. is work in progress and will be available soon.
Due to an imbalance in heating (heat entry on the bottom tool side and heat exchange with the surrounding air on the upper surface) and the presence of an exothermic reaction, the temperature profile in thickness direction is a function of time.The temperature profile is non-linear and shows clear temperature gradients in thickness direction that in turn changes with time.Consequently, the curing process is non-uniform across the thickness and takes place quicker at locations with higher temperatures.Due to the non-linear temperature profile, the chemical shrinkage on curing and the thermal shrinkage on cool-down is non-uniform and non-linear across the thickness as well, resulting in a nonlinear distribution of microscopic residual stresses.These are compressive in nature in the fibers and tensile in the matrix.Due to compressive residual stresses in the fibers and a vertical gradient of residual stresses, residual shear stresses are introduced that may contribute to the formation or amplification of undulations.This effect will be investigated in upcoming studies.

Appendix
Table 1: Curing kinetics-related material parameters used in this paper (adopted from [16]).

Figure 1 :
Figure 1: Visualization of the different model scales: a) Full-scale model; b) RVE T : 1D mesoscopic representative volume element in through-thickness direction for the thermal analysis; c) RVE σ : 1D microscopic representative volume element in fiber direction for the mechanical analysis.

Figure 2 :
Figure 2: Two methods to calculate the chemical residual stress in the microscopic RVE σ model: a) Constant strain in longitudinal direction; b) Strain difference between fibers and matrix with interphase (shear) layer connection.

Figure 3 :
Figure 3: Simulation without consideration of the exothermic reaction: Temperature (left) and curing degree (right) plotted against time.

Figure 4 :
Figure 4: Simulation with consideration of the exothermic reaction: Temperature (left) and curing degree (right) plotted against time.

Figure 5 :
Figure 5: Residual stresses in the fibers.Fiber stresses due to chemical shrinkage (left) and after cool-down (right).

Table 2 :
Material parameters used for the analysis of thermal residual stresses.

Table 3 :
Material parameters used for the analysis of the mechanical behavior of the laminate.