Method for residual strain modeling taking into account mold and distribution of heat transfer coefficients for thermoset composite material parts

The article presents a technique for process-induced residual strain modeling for thermoset composite material parts. The model takes into account the mechanical and thermal contact between the part and the mold. The technique is implemented in the ABAQUS software using user subroutines. Using the technique, it is possible to clarify the distribution of the heat transfer coefficient on the surface of the part and mold using the CFD method. Distribution of heat transfer coefficients are obtained in ANSYS CFX under the appropriate process conditions. The method is verified for the U-shaped sample. Also, the results of modeling the stringer-stiffened curved composite panel using the developed technique without taking into account the mold and heat transfer coefficient distribution are presented.


Introduction
Products made of thermoset composite materials are widespread in the aviation, shipbuilding, automotive, energy, and other industries. The main methods for obtaining such materials are liquid composite molding (LCM) and autoclave molding. As a result of the technological process, various technological defects may appear. To reduce costs at the stage of development a manufacturing process, numerical modeling is used. In general, the textile reinforcement geometry, draping, impregnation, curing, residual stresses, and strains are modeled. Simulation captures effects such as fiber shear; the formation of air voids; incomplete curing zones; excessive residual stresses that can lead to the destruction of the part during the manufacturing process; and unplanned shape distortion of the part [1][2][3]. The general goal of our research is to take into account all stages of the manufacturing process to determine the final properties of the product, including physical and mechanical properties, residual stresses, and shape distortion. This article is focused on the determination of thermal boundary conditions and subsequent modeling of process-induced residual strains.
Works [4][5][6] examines the effect of a number of design and process parameters on spring-in by experimental methods. Also, the spring-in angle is determined from analytical and semi-empirical relationships. These formulas include inplane and through-thickness thermal and isothermal chemical shrinkage, temperature change, and initial angle. In [7][8][9], curing simulation methods is presented, when interaction between part and mold is imitated with using constraints. In [10][11][12], Svanberg studied three different boundary conditions: freestanding for imitation of mold which is very soft compared to the composite; fully constrained part for rigid mold imitation; and with frictionless contact conditions between part and mold. Thermal contact was not taken into account in these models. As will be shown below, thermal contact in the model significantly influences the simulation results. In [13] a coupled thermo-mechanical problem is modeled taking into account thermal and mechanical contact. Thermal contact is assumed to be ideal. In [14][15][16], friction coefficient during cure changes for mechanical contact. The state of composite material changes significantly during cure, from viscous to rubbery and then to stiff solid. Therefore, contact changes from stick behavior to slip behavior. Various approaches and aspects of cure modeling and process-induced strains are presented in [17][18][19][20][21].
In these works, two approaches to modeling a laminated composite are used. In the first case, a layer-by-layer approach is used, when one ply corresponds to one solid element in thickness. Otherwise, the averaged properties of the laminate are used. The averaging techniques are described in more detail in [22,23].
In the works presented above, to simulate the heat transfer from the heated air of the oven/autoclave, convection boundary condition is used. For this, the sink temperature as a temperature cycle and constant typical heat transfer coefficient (HTC) on the open surface of the part and mold is set. In a real process, an uneven temperature and HTC distribution on the surface part occurs in an autoclave due to a complex air flow pattern. It is especially important to take this phenomenon into account to simulate for large-sized products of complex shape. For some of these thin-walled structures, increased requirements are imposed on shape accuracy (for example, antenna reflectors [24,25]) and the uneven distribution of HTC can significantly influence the accuracy of modeling residual strains. This fact is explained by influence of the HTC on the cure degree. Berdnikova et al. [26] studied uneven distribution temperature effect with using computational fluid dynamics (CFD) simulation. But material model in this research does not take into account the effect of resin chemical shrinkage. Airflow simulation and estimation of convection boundary conditions in autoclave are presented in [27][28][29].
In this research, the above approaches are combined and a new technique for modeling a laminated composite taking into account thermal and chemical shrinkage, as well as mechanical and thermal contact with the mold is presented. Two ways of setting the HTC are used. First, the HTC is set the same for all surfaces of the part and mold. Then the HTC are obtained from CFD analysis in ANSYS CFX. In the CFD analysis, for the part and the mold, the airflow in an autoclave is simulated. The possibility of simplifying the problem to a stationary process and transfer of HTC to a finite element software package ABAQUS are also investigated. The goal of this research is to propose an engineering approach for shape distortion prediction taking into account these phenomena.

Main relationships for thermal and mechanical behavior
Transient Fourier anisotropic heat conduction equation, with a heat generation term from the exothermic resin cure reaction Heat flux due to internal heat generation in a material In this research, Kamal-Sourour semi-empirical cure reaction formulation is used [30,31]. It can be determined by examining heat flow measurements (e.g., from differential scanning calorimetry (DSC)) where K(T) is the Arrhenius-type relation for the reaction temperature dependency; A is the pre-exponential factor; m and n are model constants; E a is the apparent activation energy; R is the universal gas constant; and T is the temperature. During cure, three states of the matrix are successively replaced-liquid, rubbery, and glassy. The process of curing the resin is associated with the appearance of a nonmechanical strains [22].
Nonmechanical strains represent the sum of strains from thermal expansion ε T ij and from chemical shrinkage ε C ij : Thermal strains are defined as: where the instantaneous coefficients of thermal expansion α ij depend on temperature and degree of cure as where α l ij , α r ij , and α g ij are linear coefficients of thermal expansion in the liquid, rubbery, and glassy states, respectively; X is the degree of cure; and X gel is the degree of cure at gelation. Chemical shrinkage is defined as where the instantaneous coefficients of chemical shrinkage β ij depend on temperature and degree of cure as, where β l ij , β r ij , and β g ij are the linear coefficients of chemical shrinkage in the liquid, rubbery, and glassy states, respectively.
The dependence of the glass transition temperature T g on the degree of cure α was defined using the Di Benedetto equation [32]: where T g0 and T g∞ are the glass transition temperatures of the uncured (X = 0), respectively, fully cured system (X = 1); and λ is a material constant. In [14,22], the main constitutive relations are given, which serve to describe the curing process of a composite material with a thermosetting matrix. These models describe the material in rubbery and glassy states.
In this paper, several defining relations of the elasticity model in simulating the process-induced stresses are used. When the modules are piecewise constant and are written as [22] where C R ijkl and C G ijkl components of the elastic module tensors in the rubbery and glassy states, respectively.
And when elastic module tensors change calculated by the Cure Hardening Instantaneous Linear Elastic (CHILE) model [14,33,34] ð12Þ where T * = T g (X) − T and T C1 and T C2 are material parameters.
To be able to simulate using the finite element method, the direct Euler scheme is used [22] where where ΔTи ΔX are are the temperature, respectively, degree of cure increment over the time step. According to [22], the coefficients for linear shrinkage in the liquid state α kl и β kl are taken as zero.

FE Implementation
А nonstationary simulation of heat transfer process according to a given curing temperature cycle in an implicit ABAQUS Standard solver is carried out with taking into account geometric nonlinearity. Relationships for thermal and mechanical behavior are implemented with using subroutines USDFLD, HETVAL, UEXPAN, and UMAT, written in the FORTRAN programming language.
With using USDFLD, the dependence of thermal conductivity and specific heat on temperature and cure degree is determined.
The user subroutine HETVAL defines the internal heat generation of the material as a result of the exothermic reaction as a function of the cure rate. The reaction rate expression is a differential equation and is solved by an improved Euler method.  With using UEXPAN, the orthotropic thermal expansion and chemical shrinkage coefficients are determined depending on current cure degree and temperature.
The UMAT was used to set the elastic properties of the material for the rubbery and glassy states of the matrix, and also to calculate the mechanical stresses.
Thus, the cure degree, the increment in thermal and shrinkage strains, and the update of the elasticity tensor are performed at each integration point at each time step, which allows taking into account the history of the process.
In the FE model, each ply represents a solid part. Each ply was represented by one solid element in thickness. The element type C3D8T is used to simulate a coupled thermomechanical problem. The laminate consists of 8 plies and has the lay-up [0/90/45/−45]s. The ply properties are presented in Table 2 [15].
The properties of the mold material were set as for aluminum.
At the «Initial» step, thermal and mechanical contact, temperature 20°C, and constraints on the bottom surface of the mold are set.
In the «Curing» step, the heat transfer coefficient on the free surface of the part and mold equal to 45 W/(m 2 • K −1 ) and the cure temperature cycle are set. On the outer surface of the part a pressure of 0.7 MPa is set.
The «Springback» step simulates the part demolding. The following boundary conditions are deactivated: heat transfer coefficient and external pressure; mechanical and thermal contact between part and mold. For the part, boundary conditions are set that allow it free to move.
The interaction between the part and the mold was defined using «Surface-to-surface contact».
In the «Normal behavior» mechanical contact property, by default, a «Hard» contact pressure-overclosure relationship is used for surface-based contact. When surfaces are in contact, any contact pressure can be transmitted between them. The surfaces separate if the contact pressure reduces to zero.
In the «Tangential behavior» mechanical contact property, Coulomb friction is used. The Coulomb friction model defines  Sidewalls --Opening critical shear stress (τ crit ) at which sliding of the surfaces starts as a fraction of the contact pressure (p) between the surfaces (τ crit =μp). The fraction (μ) is known as the coefficient of friction and is set to 0.25. Thermal contact was set using the conductance and clearance parameters.
The model has been validated based on research [11,12].

Stiffened composite curved panel FE model
The

Computational fluid dynamics model
The method for HTC determination HTC is the same for LCM and autoclave molding methods. Depending on the molding method, the geometrical dimensions of the oven/autoclave, temperature and boundary conditions will differ. For a reliable description of the processes in the autoclave, it is necessary to simulate the complete geometry, as well as the unsteady process of heating the air flow that blows the part and the mold.
The need to model the complete geometry is due to the fact that the flow pattern inside the autoclave can be quite complex with the formation of eddies. The unsteadiness of the process is caused by a change in temperature according to the curing cycle. Also, as a result of temperature changes, the air thermos-physical properties change, and, accordingly, the HTC.
Curing of a thermoset composite material in an autoclave is a nonstationary process with a given temperature regime. In order to determine the boundary and initial conditions for setting, stationary CFD-modeling tasks used the results of the preliminary coupled temperature-displacement calculation in the commercial software package ABAQUS CAE.
A stationary model of the fluid dynamics and the heat transfer in autoclave has been set up in the commercial software package ANSYS CFX discretizing the incompressible Reynolds Averaged Navier-Stokes equations and total energy transport equation. The k-ε turbulence model has been chosen as one of the most common models for technical applications.
In order to reduce the scale of the case and also CPU time, some parts of the autoclave in which the air flow is not of interest are excluded from consideration (Fig. 3).
The discretization has been done manually using an unstructured mesh approach. Three different variants of discretization with different element sizes were used in order The description of the meshes is given in Table 3. All meshes feature a high mesh density at the plate of the autoclave and the U-shaped part to provide an accurate representation of the boundary layer.
Model and boundary conditions are shown in Fig. 4. Autoclave panel is considered to be adiabatic; boundary conditions on the sidewalls are defined as «opening». The inlet velocity was assumed to be constant over the inlet plane. The direction of the velocity vector was assumed to be perpendicular to the inlet plane. The description of boundary conditions is given in Table 4.
Boundary conditions are based on the results obtained in [28].

U-shaped part with constant heat transfer coefficient
First, simulation was carried out with a constant HTC (45 W/(m −2 K −1 )) on the surface of the part and mold. Figures 5 and 6 show the temperature distribution of the part and mold without thermal contact at the moment of the greatest local overheating of the part for constant air temperature of the autoclave-110°C and 180°C, respectively. This temperature rise above the cycle temperature is due to the exothermic heat generated by curing. Figure 7 shows the temperature distribution at the moment of the greatest local overheating of the part, taking into account the thermal contact between the part and the mold. Compared to the model without thermal contact, it can be seen, that the temperature rise in this model was 0.46°C more, than the cycle temperature. This is due to the fact that the heat supplied to the part and released as a result of the cure reaction is partially transferred to the mold. Figure 8 shows the distribution of the degree of cure after gelation of the resin. In the model with thermal contact at the same time, the degree of cure is lower.
Analysis of the results in Figs. 5-8 shows that, without taking into account thermal contact, the most heated zones in the part arise closer to the mold. In the same zones, the degree of cure is higher. Taking into account the thermal contact, the most heated zones and the higher the degree of cure arise at the radius of the part on the outer surfaces. Fig. 7 The maximum temperature of the part and mold for constant air temperature (180°C) in the model with thermal contact Fig. 8 The degree of cure after gelation of the resin in the model without thermal contact (left) and with thermal contact (right) at the time 7240 s  Fig. 11 Effect of thermal contact on U-shaped part temperature (2 mm) Fig. 12 Effect of thermal contact on U-shaped part temperature (6 mm) Figure 9 shows the change in cure degree, temperature of U-shaped part, and change of temperature in autoclave over time with thermal and mechanical contact between U-shaped part and mold.
The graphs show that the temperature of the part differs (~20°C) from the temperature of the autoclave cycle at the initial and final stages. This effect arises due to taking into account the thermal contact between U-shaped part and mold which allows heat exchange through the contact surface. On the initial stage and on the final stage of the thermal cycle mold and U-shaped part, temperature values practically coincide, but heat transfer still occurs between them.
It can be seen on Fig. 10, where negative heat flux-is a heat flux into the body, positive heat flux-is a heat flux from the body. Figure 11 shows a temperature graph for a 2-mm-thick part. With thermal contact, temperature differs little from temperature of cycle in autoclave at the regions with constant temperature framed on Fig. 11. Without thermal contact between mold and U-shaped part, temperature rises become more noticeable (~2°С). In the same time with thermal contact temperature, rises are less (~0.5°C). The increase over the autoclave temperature in these areas is explained by the occurrence of exothermic chemical reactions in the part.
The thickness of the part affects the amount of temperature rise in framed areas: with increasing thickness of the part the maximum temperature gets higher. It can be seen by comparing Figs. 11 and 12. As shown in Fig. 13 with an increase in the HTC, the value of the rise of part temperature in the thermal cycle sections with a constant temperature decreases: with HTC~8 W/(m −2 K −1 ) part temperature rise over autoclave temperature cycle is 30°C; with HTC~45 W/(m −2 K −1 ) part temperature rise over autoclave temperature cycle is~5°C. This effect can be explained by the fact that at the moment of intense exothermic heat generation the temperature of the part becomes higher than the air temperature in the autoclave, and with a higher HTC, the air flow cools the part better. Figure 14 shows the final shape of the part and the change from the nominal dimensions after demolding at the «Springback» step. Figure 15 shows the maximum contact pressures acting on the part from the mold side during the curing process.

U-shaped part with heat transfer coefficient from CFD simulation
The above results were obtained under the condition of constant value of the heat transfer coefficient on the surface of the part.
To simplify the formulation of the CFD calculation, the cycle area with stationary value of temperature and zero heat flux is chosen. These sections are highlighted in Fig. 16 with a red frame.
In accordance with the above methodology and the CFD model, an analysis of the flow around a part and mold in an autoclave was carried out (Fig. 17).
According to the simulation results shown in Fig. 18, it can be seen that the HTC depends little on the surface temperature of the part. The values of the HTC averaged over the part surface are shown in Table 5.
The validation of the CFD simulation of the HTC in an autoclave was carried out according to the results given in [28].  Figure 19 shows the temperature distribution fields on the part surface at different times, corresponding to the thermal cycle areas with a constant air flow temperature (110, 180, 20°С). It can be seen from the figures that the temperature distribution over the surface of the part is nonuniform, but the difference between the maximum and minimum surface temperatures is in the range from 0.01 to 0.3°C. Figure 19 shows that for a given part, the influence of nonuniformity of the HTC does not significantly affect the temperature field distribution. The same conclusion is valid for the cure rate and degree of cure. Figure 20 shows that the degree of cure and the cure rate for cases of uneven and uniform distribution of the HTC over the part surface practically coincide.
However, for larger parts with greater thickness, the effect of nonuniformity of the HTC over its surface on the curing rate and the cure degree can be significant. Figure 21 shows the results of a stringer panel simulation in according to described methodology with constant typical HTC without account into mold. Maximum displacement from nominal geometry due to process-induced residual strains is predicted.

Conclusions
A new technique has been developed that can be used by engineers to assess residual technological deformations. The technique makes it possible to estimate process-induced residual strains with mold, depending on the requirements for modeling accuracy and available computing resources.  The presence of the mold significantly affects the simulation results. Mold through thermal expansion and mechanical contact affects the final shape of the part. Also, the mold has a significant impact on the modeling of heat conduction process between the part and the mold due to thermal contact.
The accuracy of solving the heat problem depends on the specified HTC on the surface of the part and mold. When specifying the typical HTC, using the literature data, one may encounter a large scatter of this value. The HTC itself on the surface of the part in the oven/autoclave will depend on the part geometry, temperature cycle, air properties, and boundary conditions. Therefore, it is important for parts requiring increased geometric accuracy to carry out CFD simulation to refine HTC.
In a future study, it is planned to validate the models based on the experiment results and to clarify recommendations for the method application. The influence of the size and type of time step (fixed/automatic) on the results should be investigated. Availability of data and material Not applicable.
Code availability Not applicable.

Declarations
Ethical approval Not applicable.

Consent to participate Not applicable.
Consent to publish Not applicable.

Competing interests The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article ' s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article ' s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.