Simulation analysis of deformation spring-back of free-form surface composite with cellular foam filling

Due to the complex thermo-chemical-mechanical coupling during the curing and forming process of composite materials, spring-back is commonly observed after forming, which makes the shape of parts after curing and forming deviate from the design shape, resulting in difficult or even impossible assembly. In order to reduce the influence of the spring-back deformation on the forming, the finite element analysis software was used to analyze the curing deformation of the free-form composite part with foam filling, and the numerical simulation and analysis of the curing deformation of the complex curved parts were realized. The simulation analysis involves warping deformation under different angles and radii, which breaks through the unitary shape of previous research objects. In addition, the effect of foam on skin is also considered in the simulation process, and it is found that the shrinkage of the foam will affect the spring-back deformation of parts. Through finite element software, the spring-back simulation of multi-surface composite material after curing was carried out, and the results were compared with those obtained by digital measurement, the comparative analysis shows the effectiveness of the curing deformation prediction model in this paper.


Introduction
Composite materials are widely used in aircraft due to their high specific strength, high specific modulus, and strong designability. However, the physical and chemical changes in the curing process of the composite will cause curing deformation, which will make the molded surface exceed the assembly tolerance and affect the service performance. The traditional mold compensation method is time-consuming and laborious. In recent years, the method of curing molding prediction by computer simulation has become a hot research topic [1].
The research on deformation prediction has developed from focusing only on the thermal shrinkage effect of the curing process to a more specific model which may include many deformation mechanisms and more advanced material models and analysis measures [2]. The research shows that the curing process of thermosetting composites usually goes through three stages: viscous state, rubber state, and glass state. To simulate the physical and chemical changes of these three stages, the coupling method of heat transfer-curing model, flow-compaction model, and stress-deformation model is commonly used [3]. For large and thin components, the flow compaction process can be ignored [4]. Curing deformation is mainly affected by anisotropic thermal expansion, chemical shrinkage, mold-part interaction and other factors. The equations involved include the thermochemical equation, stress-strain equation, and so on.
Based on neglecting the flow-compaction model, Çınar et al predict the curing deformation of composite components by sequential coupling method, studies the influence of thickness and radius on the deformation of composite components, and finds that the two-dimensional plane strain model is not enough to meet the complex deformation [5]. Ma Y R et al applied the material property change equation in the simulation model and predicted the curing deformation of curved components with three-dimensional solid modeling [6].
Min R et al established a curing model considering the thermal-chemical behavior, the thermal expansion and contraction behavior, the viscoelastic effect of the matrix resin, and the anisotropy of the composite. The simulation results were compared with the experiment data, and the reliability of the model was verified [7]. At present, the research on curing deformation mainly focuses on a flat panel, rib or rod, 'L', 'V', 'C', 'T', 'U', 'Z', circular, cylinder or tube, and so on [2], while the research on the free-form surface is less.
Free-form surface and foam filling structures are widely used in wing and other components. They are lightweight and economical [8], but difficult to manufacture. The free-form surface is difficult to be described by simple mathematical functions and has the characteristics of complex curvature and size changes, so it has higher requirements for the simulation process. Foam-filled structure made of light core and the composite shell is commonly used in aircraft wings. Foam filling can greatly increase the flexural and buckling resistance of the structure [9]. During the curing process, the foam can transfer the molding pressure evenly to the laying layer of the lower panel to improve the forming quality of the panel [10]. At the same time, the phenomenon of thermal deformation will also occur in foam, and its thermal expansion coefficient is different from that of composite materials. At present, there are few studies on the deformation simulation of free-form structures with foam filling.
Therefore, we established a finite element simulation model for the deformation of composite materials with foam bonding structure in the actual production of foam-filled free-form surface samples. Using the heat transfer-curing model and stress-deformation multi-module coupling method, we used ABAQUS and its secondary development subroutine to simulate the curing deformation. The model simulates the heat transfer in the curing process, the residual stress after curing, and the deformation spring-back after demolding. The reliability of the model is verified by the digital measurement experiment of the actual material.

Influence analysis of spring-back deformation
The curing process of composite components mainly includes heating, heat preservation, cooling, and demolding. In the heating stage, the curing reaction occurs, and the material is mainly viscoelastic; in the heat preservation and cooling stage, the resin gradually changes to glass state, and the composite component is linear elastic state, and the modulus is significantly higher than that of viscoelastic state. At present, it has been theorized that the residual stress that causes curing deformation mainly occurs at the linear elastic stage with high modulus of composite material, that is, during the cooling process of components [11].
In this paper, the factors such as raw material [12][13][14] (laminate and resin), process [13], curing parameters [13,14], mold material, mold structure, curing environment, heat transfer law of parts during heating and pressurization, heat release of resin internal curing reaction [15], the interaction between mold and part [16,17], and other factors will be comprehensively considered. The main factors are the anisotropy of the thermal expansion/contraction of the parts, the internal reaction heat during the resin curing process, the interaction between the mold and the part surface and the effect of foam on skin. Based on the analysis of previous work, we will research the deformation spring-back of composite components after curing based on the following thermochemical model, constitutive equation of foam, non-mechanical strain model and model of effect of foam on skin.

Thermochemical model
The temperature distribution inside the composite parts determines the degree of curing. In the curing process of composite materials, heat transfer from curing reaction of resin matrix should be considered [7,18]: Where: r c is the density of composite material; C c is the specific heat capacity of composite material; l , x l y and l z are the thermal conductivity of the composite materials in x, y and z directions; T is the temperature; t is the time; Q is the thermal generation rate, and its expression is: Where: r r is the resin density; V f is the volume fraction of fiber; H r is the total heat released by the curing reaction of unit mass resin; a is the degree of curing; a d dt / is the curing reaction rate, the expressions of different resins are different. Taking epoxy resin 3501-6 as an example, the curing reaction rate can be expressed as [12]: a a a a a a a = + --- Where: K i is the reaction rate constant of the autocatalytic model; A i is the frequency factor of the autocatalytic model; DE i is the activation energy of the autocatalytic model; R is the ideal gas constant.

Constitutive equation of foam
In recent years, with the large-scale development and utilization of polymer materials, the mechanical properties of many new polymer materials belong neither to the category of elastic theory, nor viscous theory. The viscoelastic mechanical behavior of foam materials under large deformation exhibits nonlinear characteristics. The deformation of the materials does not abide by the assumption of small deformation. Therefore, the coordinate of each particle within the material in the reference configuration and in the current configuration must be considered. It can be obtained by analysis that the nonlinear viscoelastic constitutive model independent of the coordinate system of foam materials under large deformation [19] is: Where s ij is Cauchy stress components; p is hydrostatic pressure; d ij is Kronecker delta; r is mass density; t is current time; F is deformation gradient; W is Helmholtz potential per unit mass; E is green strain tensor; T is the absolute temperature.

Non-mechanical strain model
During the curing process, the strain of composite is composed of thermal strain e th caused by temperature change, shrinkage strain e sh caused by curing shrinkage of matrix resin and effective strain e e caused by mechanical load. e th and e sh are called non mechanical strains. Accurate characterization of the changes of e th and e sh is one of the prerequisites for predicting the curing deformation of composites [3].
The thermal expansion coefficient a r of pure resin is a function of curing degree and temperature, When the temperature difference is DT, the thermal strain increment e D r th of the resin is: After predicting the thermal expansion coefficient of the composite, the thermal strain change of the composite in the incremental step can be obtained by the following formula: The curing shrinkage strains e sh 1 and e th 1 in the main direction of the material are also derived by the above method.

Effect of foam on skin
During the forming process, the effect of foam on the skin is mainly manifested in the expansion of the foam and the expansion pressure P of the skin. The expression of the expansion pressure P is [20]: The K v is the bulk modulus of the foam, and the a v is the coefficient of expansion of the foam. The T and T 0 are the current temperature and the initial temperature respectively. Due to the difference of thermal expansion coefficient between the outer skin composite and the foam, the interfacial effect will cause residual stress in the component when the temperature changes.
In the formula, E g is the elastic modulus in x-direction of the composite, T 0 is the initial temperature, T is the current temperature, a 1 is the thermal expansion coefficient of the x-direction of the composite, the calculation method is shown in equation (6), a s is the thermal expansion coefficient of the x-direction of the foam sandwich structure, the calculation method is as follows: Where E g is the elastic modulus in x-direction of the composite, a a is the thermal expansion coefficient of the resin, E a is the elastic modulus of the resin, A g is the cross-sectional area of the skin, A f is the cross-section area of the foam, A ai is the cross-sectional area of the interfacial resin, a f is the thermal expansion coefficient of the foam, a ai is the thermal expansion coefficient of the interfacial resin, and E ai is the elastic modulus of the interfacial resin.
Considering the thermo-chemical model and stress-strain relationship of the skin and the foam, the temperature and residual stress at each point can be obtained.
At the same time, the formula for thermal expansion of PMI foam is as follows [22]: Where, a L is the thermal expansion coefficient of the foam, DT is the temperature difference, e D is the thermal strain increment of the foam.
The thermal expansion coefficient a L of PMI foam is´ --3.7 10 C , 5 1 which is larger than that of CFRṔ  -5 10 C 7 and´ -3.53 10 C. 5 Meanwhile, in the heating stage of the curing process, the composites present the flowing state and the rubber state respectively. Due to the fluidity of resin in the flow phase and the low modulus of resin in the rubber phase, the contribution of curing shrinkage and thermal expansion to the residual stress of the composite is small. In the cooling stage, the composite material is glassy and the resin modulus is larger. In this stage, the thermal deformation mismatch between composite components and other parts has a great influence on the residual stress and deformation [3], so the influence of this stage should be considered when investigating the residual stress and deformation. In the high curing degree and large modulus stage of the cooling process of the composite part, the cooling effect of foam will also occur. The contraction of the foam will drive the solidified skin to shrink inward, causing the curvature radius and bending angle of the skin to become larger.

FEA numerical simulation
The research object of this paper is the free-form surface (it is equal to free-form modeling for CAE) composite part shown in figure 2, the spanwise length of the part is 800 mm, the chord length is 580 mm. The skin of the part is composed of composite materials laid continuously, the skin material is 3238A/CF3052, the thickness of the single layer is 0.31 mm, there are four layers in total, and the laying angle is [0°/0°/45°/0°]. This stacking sequence comes from the actual manufacturing part and is designed to meet the requirements of tension, shear, bending moment, and torque bearing of the whole component.
The skin is opened from the trailing edge. The inner skin is filled with D shaped beam and foam. The parts are glued together, and the area of each part is shown in figure 1. The section width of D beam is about 232 mm and the section height is about 66 mm. The foam filling material is isotropic material A52 foam, which is a kind of polymethacryl imide (PMI) material. The material is cross-linked closed cell foam [21]. It is often used in the core of aerospace composites. The chord length of foam is 295.315 mm, the vertical chord width is 66.281 mm, and the cross-section area is 0.012 m 2 . In this paper, the finite element simulation software ABAQUS is used, combined with the user subroutines written in FORTRAN language. In the simulation, the anisotropic elastic modulus, Poisson's ratio, shear modulus, thermal expansion coefficient, thermal conductivity coefficient, density, specific heat, thickness and direction of carbon fiber laminate, as well as the change of curing degree, the change of material properties with curing degree, and the internal heat generation of material curing should be considered.
The part and D beams are carbon fiber reinforced polymer (3238A/CF3052), and foams are made up of A52 foam. In the simulation, the elastic modulus, Poisson's ratio, thermal expansion coefficient, heat conduction coefficient and specific heat capacity of the composite are mainly set. The material of the external molds is invar steel. The elastic modulus, Poisson's ratio, density, specific heat capacity, thermal expansion coefficient, and thermal conductivity coefficient of the molds are defined in the simulation. The material properties in the simulation are shown in tables 1-3.

Settings of FEA parameter
For the settings of the FEA parameter, we comprehensively considered various factors in the curing process: anisotropic thermal expansion of composite materials, curing chemical shrinkage, etc To simplify the simulation calculation, the molds will not consider the heat loss, and keep the molds temperature change consistent with the molding process design temperature change curve.
The process condition of FEA simulation is hot pressing forming. The simulation mainly includes two processes: molding and spring-back. That is, first, the part and internal and external molds are set in the process conditions to complete the curing molding; then the external molds are removed to complete the process of spring-back. The molding process parameters are shown in figure 2. The initial pressure is added to the part after molds closing, and the pressure is set as (6±5MPa), and then the temperature is increased to 70±5°C at a rate less than 1.5°C min −1 . After 30±5 min of heat preservation, external pressure is added: the pressure is increased to (12±0.5MPa), and the temperature is further increased to 125±5°C and kept for 135±15 min, and then the pressure is relieved after cooling to 25°C. The spring-back process removes the external molds and sets the temperature of the predefined field to 25°C.
In the setting of the curing simulation environment, the pressure is applied to the contact surface of the upper mold and the part, and the mold and the part are compacted. The lower surface of the lower mold is completely fixed; the relationship between the external mold surface and the part surface is set as a surface-tosurface contact, and the relationship between the foam and the internal surface of the part is set as tie constraint.   Table 3. Invar steel material properties.

Parameters Value
Density r m /t/mm3´-8.1 10 9 Thermal conductivity /mW/(mm·K) 11 Specific heat C m /mJ/(t·K)5. 15 10 8 Elasticity (Young's modulus/kPa-Poisson's ratio) 1250000.22 Expansion a m /°C-1´-1.6 10 6 In the setting of the spring-back environment, the surface relationship between foam and part is still tie constraint. Different from the spring-back simulation model in the existing literature, the foam is still filled inside the parts after forming studied in this paper, and foam is cemented with the inner surface of the part. Therefore, the influence of cellular foam on the spring-back of the part should be considered in the spring-back process.
Although ABAQUS has powerful functions, the material library of ABAQUS software cannot analyze the resin curing reaction. Therefore, the secondary development function of ABAQUS software needs to be used [7]. The user subroutines of USDFLD, DISP, HETVAL, and UEXPAN are used in the numerical simulation of curing deformation of composite materials. The user subroutines of USDFLD, DISP, and HETVAL are applied in the heat conduction curing analysis module, and the UEXPAN user subroutine is applied in the stress displacement module [23].
USDFLD is used to define the initial curing degree, which is used to define the field variables as functions related to time or any other variables in the output variable list. The defined variables can be called by other subroutines to realize the data transfer between subroutines; HETVAL is mainly used to define the internal heat generation in the heat transfer analysis, and calculate the heat flow caused by the heat generation of the curing reaction of the part, to obtain the temperature and curing degree of each incremental step; DISP is used to apply temperature boundary conditions; UEXPAN represents the thermal expansion coefficient of the composite and

Analysis of numerical simulation results
FEA numerical simulation was used to complete the whole process, including curing and spring-back. The parts experienced three stages: pre-curing, post-curing (pre-spring-back), and spring-back. Figures 4 and 5 show the stress and strain nephogram of the part after curing. It can be seen from the stress nephogram that the maximum stress is 88.83 MPa, mainly concentrated at the position of the wing surface corresponding to the foam, while the stress of the corresponding part of the D beam is smaller.
The stress distribution of the part in each region is relatively uniform. The displacement nephogram shows that the maximum displacement is 1.027 mm. Displacements are distributed uniformly on the whole, in which the displacement at the spanwise edge of the part is larger than that in the middle.
The main reason for the difference between the two regions in the stress nephogram is that the D beam has stronger load-bearing capacity and the foam is softer. Therefore, during the molding process, the position of the wing surface corresponding to the foam deforms seriously due to insufficient supporting force. Besides, due to the special properties of the material, the uneven heat release of the resin in the curing process makes the part heated unevenly, which is also one of the reasons for the difference of stress and deformation at different positions of the part.  According to the results of the stress and strain nephogram, it can be found that a large amount of residual stress remains in the part after the heating and pressure curing process, especially in the position of the wing surface corresponding to the foam. After the part is cured and cooled and the mold pressure is relieved, the internal residual stress of the part will gradually release, which will cause spring-back and make the part deviate from the design shape.
The deformation and spring-back process of the part is carried out on a stable wooden frame, and the surface of the wooden frame is paved with a cotton protective layer for protection. Under the supporting force and selfweight, the parts are deformed and rebounded.
Therefore, according to the actual working conditions, we simulate the spring-back process in the finite element software. The setting shelf is shown in figure 6, in which the contact part between the frame and the wing surface is made of cotton material; the lower part of the frame is made of wood material, and the lower bottom surface of the frame is fixed. Figures 7 and 8 respectively show the displacement nephogram of the part and the shelf after the spring-back process of the part. It can be seen that due to the soft material of the protective layer, the position of the leading edge of the part at the beginning is higher, and the position of the center of gravity of the part is closer to the leading edge, the displacement of the leading edge and the contact position of the part is larger, and the maximum displacement is 9.314 mm. At the same time, the displacement of the part has a certain angle of torsion along the spanwise and chord directions, which may be related to the initial placement position of the part and the torsion angle of the wing surface itself.   From the analysis of the strain and stress nephogram after curing and spring-back, it can be concluded that the part deformed after curing due to the combined action of material shrinkage and internal stress. With the removal of external mold, the internal stress release of foam and part further increased the spring-back of the part.
According to the simulation results, the whole FEA model can be used to simulate the heat transfer, residual stress evolution, and curing deformation after demolding. Finally, to verify the effectiveness of the simulation results, a comparison analysis is made between the simulation results and the measured numerical model of the actual hot-pressed part.

Analysis of measured numerical model and FEA simulation results
Three models were used in this research. (1) Design model: the model designed by the design department according to the actual use requirements. As the input of the simulation, the model shows the surface before curing deformation. It is used as a base model for comparison with other models, as shown in figure 9. (2) Simulation model: a finite element analysis model derived from the CAE post-processing program, as shown in figure 10. By comparing the simulation model with the design model, that is, the undeformed model, the simulation deformation can be obtained. (3) Scanning model: the surface model obtained by using a 3D scanner to scan the actual component surface and then through point cloud fitting, as shown in figure 11. Compared with the design model, the deformation of the component in the actual manufacturing process can be obtained.

Measured numerical model
The measured numerical model of the part is realized by the photogrammetric system combined with the scanning method of blue light 3D scanner. The dense point cloud of the target surface is obtained first, and then the surface is formed by post-processing. In the process of acquiring the integral point cloud of the part, the photogrammetric system mainly provides the splicing field for the scanner, imports the splicing field data into the scanner software, and then uses the scanner to obtain the point cloud data in each splicing field, and then splices through the mark points of the splicing field. The flow chart of modeling is shown in figure 12.
The V-STARS/S8 photogrammetric measurement system is used, and the 3D scanner is COMET L3D-5M blue light 3D scanner. After the integral point cloud is measured, Geomagic Design and Geomagic Control software are applied for post-processing. The specific process is shown in figure 13. The operation of the point  cloud and polygon phase is completed by Geomagic Design, and the error comparative analysis with the numerical model is completed by Geomagic Control.

Comparative analysis of results
After the part surface is generated, import it into Geomagic Control software, and then import the design numerical model to be compared. The design numerical model is taken as the reference data, and the measured part surface is taken as the measurement data for the best fitting alignment.
After meeting the required fitting degree, run 3D comparison to check the error of each point on the measurement surface in the normal direction compared with the simulation results. The comparison results are shown in figure 14 below. The green part indicates that the error is less than 1 mm; the yellow part indicates that the error of this part is between 1 mm and 1.5 mm, and it is the positive error of the designed numerical model protruding scanning surface. From the error distribution, it can be seen that most of the measured surface has high forming accuracy (the error is less than 1 mm). The larger curvature side has an error of 1 to 1.5 mm at the trailing edge, and the smaller curvature side has a large range of positive error near the leading edge. In this range, the numerical model is outside the scanning surface, which indicates that there is a certain spring-back in this part when forming the part. A small piece of blue in the upper right corner of the comparison results of the smaller curvature part, that is, the part of the scanning surface outside the numerical model, is considered to be caused by the too thick developer or accidental manufacturing error due to the small area.
To obtain the comparison error between the simulation results and the design numerical model, the simulated cured surface is taken as the measurement data in Geomagic Control software, and the design numerical model is used as the reference data. The Geomagic Control software also carried out the optimal fitting alignment and 3D comparison operation, and the comparison results as shown in figure 15 and the error histogram as shown in figure 16 were obtained.
It can be seen from the results that, similar to the comparison results between the scanning model and the design model, the comparison results between the simulation model and the design model also present a large range of positive error of 1 to 1.5 mm on the trailing edge part of the side with larger curvature and a small range of positive error in the leading edge of the side with smaller curvature. The error of other parts is less than 1 mm. It can be seen from the histogram that the average error is 0.0791 mm, and most of the errors are within 1 mm.
Comparing the simulation results of part deformation spring-back with the part design model, it is found that the local surface of the part simulation model is located in the inner part of the part design model, and the error area is basically the same, both in the smaller curvature side near the leading edge and the larger curvature side near the trailing edge. That is to say, the simulation results show that after the forming and demolding of the part, the whole part has an inward spring-back deformation. In the simulation results, the error areas of the larger curvature side near the trailing edge and the smaller curvature side near the leading edge are smaller than  the measurement results, which indicates that the spring-back degree of the simulation results is smaller than the actual situation.
It can be seen from the above analysis that there are some differences between the scanning results and the simulation results in terms of error areas and error distribution. To know the differences between the scanning results and simulation results more intuitively, the simulation results are used as reference data, and the scanning results are taken as the measurement data. The error distribution is checked by 3D comparison. The results are shown in figure 17.
From the comparison results, the measured results and simulation results in most areas of the spring-back of difference are not big, while a small part of the area appears blue negative errors, indicating that the scanning results are outside the simulation results, which indicates that the actual degree of inward spring-back is less than that given in the simulation, which may be related to the more idealized boundary conditions given in the simulation, such as Tie constraints. To sum up, the spring-back results obtained by the simulation are basically consistent with those of the actual manufacturing part. The spring-back mainly occurs in the area where the smaller curvature side near the leading edge and the larger curvature side near the trailing edge, both spring-back to the inner part.
Referring to the relevant literature and data, it is found that the simulation results of inward deformation spring-back of the part are contrary to the existing simulation spring-back results of composite materials, as shown in figure 18. In the research on spring-back simulation of composite parts with L-shape [12], Ω shape [24], and V-shape [24], it is found that the composite parts tend to spring back to the direction of smaller curvature radius and bending angle after forming. However, the simulation results of inward spring-back are contrary to it, and the inward spring-back increases the curvature radius of the part surface.
By comparing the simulation model in this paper with the simulation model in the existing research results in terms of molds, parts materials, boundary conditions, loads, and other relevant conditions, it is found that  there is no cellular foam structure inside the part in the existing research results, and when the composite part spring-backs, the part does not interact with other molds. The main reason for the spring-back of part is that the interaction between foam and part leads to spring-back of the part after forming. To prove the correctness of the hypothesis, the adhesive contact between foam and part is modified to normal contact, and the comparison between the simulation results and the design model is shown in figure 19.
According to the contrast results, it can be found that when the foam is not bonded to the part in the simulation, the parts will not be affected by the inward contraction during the cooling of the cellular foam. Therefore, when the residual stress is released, the surface of the part will spring-back outward. It is proved that the main reason for the spring-back is the interaction between the foam and the part. When the part is released from the mold to cool down, the residual stress inside the part is released, resulting in the spring-back trend of the part. At the same time, the foam shrinks completely when it is cooled down, and the part spring-back inwards eventually.
The differences between simulation results and scanning results are also affected by the fact that the setting of simulation conditions has not been fully restored to the actual situation, the errors generated by the mold and process in the actual manufacturing process, and the errors caused by the surface alignment of Geomagic Control software. Therefore, the simulation results can accurately reflect the spring-back of the part.

Conclusion
In this paper, ABAQUS finite element simulation software combined with multiple user subroutines is used to simulate the deformation and spring-back of thermo-chemical coupling curing deformation of free-form surface composite part with foam filling, and the reliability of the model is verified by digital measurement and comparative analysis technology. The analysis explored the effect of foam filling on the curing deformation of the complex curved surface with different angles and radii and revealed the deformation and distribution of the part after curing and spring-back. The results show that the part has different degrees of deformation after curing and spring-back, and the main cause of deformation is stress release. The interaction between the part and the foam will affect the spring-back trend of the part. Due to the adhesion between the part and the foam, the foam will shrink with the release of internal stress of the part and the foam, which will impede the outward deformation of the part. Moreover, due to the large deformation of the foam, the part will also deform inward.
The method can be applied to the study of curing deformation of the thin-walled shell part with foam core structure. It provides a good analysis model for deformation prediction and provides an important reference value for deformation control of parts forming. This method can be further extended to the study of deformation analysis of other types of core filling structures and other types of parts. The simulation method of foam-skin interaction can be further refined in future research.