A three-dimensional viscoelastic analysis of thermoplastic resin matrix composite laminates during hot stamping

In this paper, a three-dimensional thermo-viscoelastic constitutive equation is used to simulate and analyze the hot stamping process of thermoplastic resin matrix composites. Dynamic mechanical analyzer (DMA) was used to measure the relaxation behavior of resin matrix at different temperatures. Then, the relaxation times and the weight factors of the resin matrix were obtained by regression analysis of the data of the resin matrix relaxation behavior using genetic intelligent algorithm. According to the classical Maxwell model with N Maxwell elements and representative volume element (RVE), the integral constitutive model was established, and the corresponding integral finite element program was incorporated into commercial software ABAQUS with an UMAT subroutine.Thermal expansion was also taken into account in the present model. The constitutive model and its finite element program in this paper were verified by experiments with a self-designed hot stamping die for thermoplastic composites.In addition, the simulation results were compared with the actual failure areas to further verify the applicability of the constitutive model adopted in this paper. The results indicate that the ball edge of the hemispherical part is easy to be damaged.


Introduction
Polymer matrix composites are widely used in aerospace, shipbuilding, civil infrastructure and other highperformance fields because of their high specific strength and specific stiffness [1]. Continuous fiber reinforced thermoplastic resin matrix composites have lots of advantages like high toughness, high damage tolerance, easy forming and processing, recyclable and so on. The application potential of continuous fiber reinforced thermoplastic resin matrix composites has attracted extensive attention of researchers [2]. Thermoplastic fiber reinforced composites are composed of different thermoplastic resins (PPS [3], PEEK [4], PP [5], etc.) and different fibers (glass fiber [6], carbon fiber [7] and PBO fiber [8], etc.). The methods of composite molding include autoclave molding, resin transfer molding, vacuum-assisted molding, die forming,and etc. Despite the methods mentioned above,hot pressing is more suitable for continuous fiber reinforced composites for its advantages of short production cycle, suitable for mass production, adaptability to materials with complex shapes, and the ability to complete parts forming and curing in one step [9,10]. Schematic diagram of hotpressing molding of composite material is shown in figure 1. The prepreg is preheated, hot pressed, cooled, and demolded to obtain a workpiece.
The forming defects such as fiber breakage, delamination and wrinkling may appear in the process of continuous fiber reinforced resin matrix composites deep drawing. De Luca et al [11] studied the forming quality of the parts on the forming speed. The stress of the composite material processing and the residual stress afterwards have an important influence on the failure and springback deformation of the composite material. At present, researchers mostly use experimental methods to study the defects of composite hot stamping forming. However, it is usually expensive and time-consuming in this method. Through numerical simulation, the economic cost of research can be reduced and at the same time, and a deeper understanding of the stamping process can be obtained. Because of the advantages of numerical simulation, the establishment and use of a constitutive model that meets the characteristics of thermoplastic composite materials is definitely required.
The presence of the resin matrix makes the composite material gain the stress relaxation and creep properties [12]. The hot stamping of composite materials is mostly carried out at a temperature close to the melting point of the resin [13], so the viscoelasticity of the composite material is more obvious during the hot stamping process. The current constitutive models used in the molding process of composite materials mainly include linear elastic models [14], path-dependent models [15] and viscoelastic models [1,6,16]. The performance of composite materials is related to multiple factors such as temperature and fiber percentage,etc. The linear elastic model regards the performance of composite materials as constants, which obviously has great limitations. Bogetti et al [17] and Johnston et al [18] improved the linear elastic model by changing characteristics of the mechanical properties of thermoset composites with curing degree and temperature, and proposed different linear elastic models of CHELE. However, the nature of the two modified CHELE linear elastic models is still that the linear elastic constitutive model cannot describe the stress relaxation and creep properties of composite materials [19]. Studies on the viscoelastic properties of composite materials in this field can be traced back to the analysis of viscoelastic stresses of composite materials conducted by Schapery [20]. Subsequently, the thermosetting resin matrix composites were further studied under the condition of autoclave. Researchers have carried out a series of studies on thermal expansion, solidification shrinkage, equivalence principle and other factors in forming [6,[21][22][23][24][25][26][27], integral and differential forms of viscoelastic finite element, and their accuracy and efficiency [6]. At present, the research objects of viscoelasticity of composite materials are mostly thermosetting composites [1,6,[22][23][24][25][26], while researches on viscoelasticity of thermoplastic composites are rare and mostly focused on the normal temperature level [28,29].
When obtaining the performance parameters of composite materials, the advantages of using micromechanics to calculate the performance parameters of composite materials according to the performance parameters of the matrix resin and reinforcing fibers are obvious. Firstly, the performance of the fiber basically does not change during the molding process of the composite material, so the work can be focused on the measurement of the matrix resin. The resin is an isotropic substance, and its performance is easier to determine than anisotropic composite materials. Secondly, the performance of composite materials will change with temperature, fiber content, etc. The use of mesomechanical expressions is more conducive to numerical simulation. Simple models based on mesomechanics include Reuse and Voigt models, and complex Periodic Microstructure Model(PMM) [30] and Self-Consistent Field Micromechanics model(SCFM) [31] with periodic distribution. In addition, many numerical models also consider the interface phase to study the influence of the interface on the mechanical properties of composite materials [32].
In this paper, considering the factors such as thermal expansion and stress relaxation during the processing of composite materials, the generalized Maxwell model and SCFM model as well as the 3D integral thermoviscoelastic constitutive equation are adopted to establish the 3D integral program and apply it to the numerical simulation of thermal stamping of thermoplastic composite materials.

Theoretical modeling and performance research
The properties of thermoplastic continuous fiber reinforced composites are determined by the resin matrix and its reinforced fibers, and the mechanical properties of the composites also vary with temperature, fiber content percentage and other parameters. At present, the macro and micro research of composite materials provides a method for researchers to calculate the properties of composites in different states [31].

Thermodynamic model
The governing equation of the heat transfer model is as follows: Where ρ is the density of the material; C P is the specific heat capacity of the material; k x , k y , k z are the thermal conductivity in the x, y, and z axis, respectively; T is the temperature of the material at time t;  Q is the rate at which heat is produced inside the material.
Given that the thermal parameters of matrix resin and reinforced fiber are known, we can calculate the thermal parameters of composites by rule of mixture of composite. The mesomechanics theory is based on the following three assumptions: (1) The fibers are evenly distributed in the matrix; (2) The fibers and the matrix maintain an ideal bonding state, that is, the two are bonded together chemically or physically. The surface is in direct contact, and will not detach or slide before it is destroyed. (3) The total volume of bubbles and pores in the composite material is small enough to be ignored. In the law of composite material mixing, the weighted average method is employed to predict the various properties of fiber-reinforced composite materials [32].  Where C, k, CTE are specific heat capacity, thermal conductivity and thermal expansion coefficient respectively; V f  is fiber volume fraction; E, v are elastic modulus and Poisson's ratio respectively. The subscript p stands for composite material, the subscript r represents resin material, the subscript f denotes fiber material, the subscripts 1, 2, 3 stand for directions, and 1 refers to the direction parallel to the fiber. Micromechanics constructs the effective relationship of composite material as a whole by means of representative volume element(RVE). The volume representative element commonly are used to predict the properties of composite materials by micromechanics including A, B and C, as shown in figure 2. In this paper, SCFM model is adopted and the expression is shown in equation (3) [31]. Where, E, G, v are respectively elastic modulus, shear modulus and Poisson's ratio, subscripts r, f represent matrix resin and reinforced fiber, subscripts 1, 2, 3 stand for directions, and 1 denotes the direction parallel to the fiber. The properties of grade E glass fiber are shown in table 1 [33], and the properties of resin matrix PP are shown in table 2.
The relationship between the shear modulus and the elastic modulus of the matrix resin can be obtained by equation (4), where the superscript 0 represents the unrelaxed state, and the superscript ∞ represents the completely relaxed state. The viscoelastic properties of composite materials are mainly determined by the resin matrix, and the relaxation times of the composite material can be taken as the relaxation times of the resin matrix [1,29,34].

Stress relaxation behavior of resin matrix
The continuous glass fiber reinforced polypropylene composite, produced by Zhejiang Shenggang New Materials Co., Ltd., is chosen as the experimental material in this paper and the matrix resin is modified PP provided by it too. The polymeric matrix of polypropylene resin is a crystallized high polymer whose molecular aggregation state is affected by temperature and external force, thus showing creep and stress relaxation characteristics. In order to accurately describe the development of mechanical properties, time dependence and temperature dependence of material behavior should be considered. The stress relaxation test of resin matrix Parameters Value was carried out using dynamic mechanical analyzer (DMA). The main curve and offset function of resin matrix were obtained by the superposition principle of time and temperature. Small beam-shaped test specimens (80 mm × 9.7 mm × 3.4 mm) were manufactured to perform relaxation experiments in EPLEXOR1500N model DMA instrument. Figure 3 shows the experimental apparatus for tensile stress relaxation tests. Tensile stress relaxation tests were carried out at different temperatures. After the samples were buffed, the temperature was raised to the setting temperature and then kept for 30 min, and then the stress relaxation tests at the setting temperature were performed. Figure 4 shows the stress relaxation data of the resin matrix obtained at different setting temperatures. Figure 4(b) shows the experimental data plotted on a logarithmic time axis.
After obtaining the stress relaxation data of the resin matrix at various temperatures, the principle of timetemperature equivalence is needed to construct its master curve, Where ξ is the reduction time; a T is the conversion factor. The reference temperature of the master curve is selected as 30°C. The model describing the relaxation behavior of the resin matrix includes a generalized Maxwell model containing N Maxwell units [1].
Where E ∞ is the fully relaxed modulus (equilibrium modulus); E w is the stiffness of the element related to the stress relaxation time. If N increases infinitely, the integral form of the model is obtained, In theory, the relaxation time spectrum can be obtained from the data of relaxation time distribution in frequency domain measured by dynamic mechanical instrument using Laplace or Fourier transform. However,

Parameters
Value 168.46 the form of relaxation time spectrum function is very complex, so its analytical form is seldom used. Most researchers adopt the Alfrey approximation form like: Equation (8) shows that the relaxation time spectrum can be directly determined by the slope of the master curve.
In thermodynamics, the stress relaxation curve of thermorheological simple materials(TSMs) can usually be described by the power law or use discrete exponential series to describe, Where b is the material constant; E ∞ is the fully relaxation modulus; E u is the unrelaxed modulus; ξ is the reduction time; Ww is the weight factor; τ w is the discrete relaxation time. In the previous studies, a single b constant can not accurately describe the stress relaxation behavior of the matrix resin. In addition, the exponential series model provides a convenient viscoelastic solution technology, such as time superposition integral calculation of recursive formula. Therefore, a discrete exponential series model is used to describe the stress relaxation behavior of resin matrix in this paper.
In order to obtain model parameters,genetic algorithm was used to regress experimental data in this paper. Genetic algorithm is an algorithm that simulates the biological evolution mechanism of nature, that is, it follows the rule of survival of the fittest. Genetic algorithm has good global search ability and fast convergence speed, and is suitable for solving complex nonlinear data fitting problems [35]. In this paper, the population size in the genetic algorithm is 500, the value of the population generation is 10 000, the mating probability is 0.9, the mutation probability is 0.09. The fitness function set in the article is is the i-th value in the real (measured) data, and f model (i) is the i-th value in the model forecast series. Figure 5(a) shows the maximum and average values of the fitness function in the population in different generations. Figure 5(b) shows the comparison between the regression curve obtained by genetic algorithm and the stress relaxation experiment data. The reference temperature of the master curve is 30°C.
The coefficient of determination of the regression equation R 2 is 0.9998. The model parameters obtained by genetic algorithm are shown in table 3.
The same mechanical state of a material can be expressed at a higher temperature and in a shorter period of time, or at a lower temperature for a longer period of time. This relationship is the principle of temperature equivalence, Where a T is the conversion factor. The WLF equation is generally used to describe the relationship with the temperature T from the glass transition temperature T g of the matrix resin to about 100°C above it. The empirical form and its parameters are given out of this range. Figure 6 is a comparison between the curve obtained by using equation (12) and the experimental data. The relationship between a T and temperature T and its parameters obtained by WLF equation and empirical form are shown in equation (12).

Thermo-viscoelastic constitutive law for composites forming
The linear viscoelastic model can be expressed in integral and differential forms. At present, both integral and differential forms are used in resin matrix composites. In this paper, the classical integral linear viscoelastic model is adopted. The expression of integral linear viscoelasticity is shown in equation (13) [1,6,36].
where C ij is the material stiffness tensor; t and t′ refer to the current and the past time, respectively; σ i is the stress vector; e j eff is the strain vector and can be expressed as: Figure 5. Genetic algorithm processing results of the relaxation behavior data of resin matrix. where ε j and e j eff represent the total strain component and the nonmechanical strain vector, respectively; CTE eq j is the equivalent thermal expansion coefficient (CTE) and can be expressed as: where β j is the coefficient of hygrothermal expansion; CTE j is the thermal expansion coefficient of the composite material and c is the moisture content, i.e. the mass of the absorbed water per unit mass of material.
In view of the thermo-rheological property of the polymer material, both the relaxed function and the constitutive relationship can be simplified to decouple the temperature effect using the time-temperature superposition principle. When there is no strain history at t 0, equation (13) can be expressed as the following form: The viscoelastic properties of the material are expressed by the generalized Maxwell model as shown in figure 7, where the unrelaxed form of C ij u m of the m-th Maxwell element is [6]: The expression of C ij u is [6]:  The relaxation times of the composite material can be considered to be approximately the same as the relaxation times of the resin matrix. So the form of C ij is [36], According to equations (16) and (19), we can get the following equation, On the other hand, the relaxed functions ( ) of the composite at reduced time x +D t t can be written as: the state of stress of the composite at reduced time x +D t t can be expressed as: According to equations (20) and (21), we can get the following equation, In order to obtain the integral I 1 , the following formula is defined: Then S i at reduced time x +D t t can be obtained: In the case of small time intervals, there is such an approximate relationship [36], is determined from the previous time step. I 3 can be determined as: Then I 1 can be recursively determined as: Using the same method as I 1 , I 2 can be integrated in the closed form: where: Substituting equations (32) and (33) into equation (24), the stress increment during the time interval Δt can be expressed as: In order to verify the accuracy of the subroutine, a self-designed hot stamping die and thermoplastic composite material were used to conduct composite hot stamping experiments and the experimental results were compared with the numerical simulation results.
The mold device is shown in figure 8. The mold device monitors the heating process of the mold through the temperature control equipment, and the mold can be heated through the heating tube. The power of a single heating tube is 100W, and the efficiency of adding heat pipe in the numerical simulation process is set as 0.75. There are 5, 8 and 16 heating tubes for punch, blank holder and die respectively. In this paper, the die material is 45 steel with a density of 7.85 g cm −3 . Its thermal parameters, specific heat capacity C and thermal conductivity coefficient k are shown in table 4.
Temperature sensors are installed in the punch, blank holder, and die to detect the mold temperature. The blank in the experiment is a 3-layer composite prepreg (single layer thickness of 0.3 mm) laid up in accordance with [0/90/0]. After the mold is set to the preheating temperature and heated, the temperature of the mold is monitored and adjusted by the temperature sensor. The workpiece deformation process is shown in figure 9. NT11 represents the temperature (°C) in figure 9.
After the die reaches the preheating temperature, the heating process ends and the punch presses the composite material at a certain stamping speed. The blank holder pressure was 0.15 Mpa. The friction coefficient between the composite material and the die parts is set as 0.2 [37]. The stamping process of composite material is shown in figure 10. Figure 11 shows the relationship between stamping stroke and drawing force at different punching speeds V. The figure also shows the contrast between the shapes of the simulated pieces and the actual pieces when the punch stroke is 20 mm.
After cutting, embedding and polishing, the pieces with different drawing depths were observed and photographed under an optical microscope, and finally the pictures were spliced to obtain the pieces'   microscopic diagrams. Figure 12 shows the comparison between the cross-section shapes of the real and the simulated parts (V=50 mm Min −1 ) at the preheating temperature T=140°C. Figure 13 shows the stress distribution of the part with preheating temperature of 140°C (V=50 mm Min −1 , H=20 mm). The X direction is the fiber direction of the first layer and the third layer of prepreg, and the Y direction is the fiber direction of the second layer of preimpregnating material. Path 1 and path 2 in the figure are the paths where parts may fail. It can be known from the simulation results that the stress of the part is distributed symmetrically to path 1 and path 2. The maximum stress or stress component is also distributed on path 1 or path 2. The shape of the failure area of the workpiece is also symmetrical with respect to path 1 or path2.

The influence of forming factors on parts
The preheating temperature of the composite during hot stamping will influence the drawing force and the stress distribution of the composite. The stress distribution in hot stamping process is closely related to the  fracture, wrinkle and other defects of composites. In this paper, the validated three-dimensional thermoviscoelastic integral constitutive model is used to analyze the effect of the preheating temperature on the stamping of composite materials under different preheating temperatures. Figure 14(a), (b) and (c) show the stress distribution of the parts (V=50 mm Min −1 , H=20 mm) with preheating temperature of 110°C, 125°C and 140°C respectively. Figure 15 shows the change of the longitudinal stress σ 1 along path 2 and the change of transverse stress σ 2 along path 1 of the parts (V=50 mm Min −1 , H=20 mm) under different preheating temperatures. It can be seen from the figure that the maximum stress of the parts occurs near the waist. The increase of preheating temperature can reduce the stress value of the corresponding position of the workpiece to a great extent.
As shown in figure 16, the distribution of the longitudinal stress σ 1 of the workpiece varies with the depth when the preheating temperature of the workpiece is 140°C. In this paper, two points on the top and the waist of the workpiece that are most likely to fail are analyzed.
With the increase of the drawing depth H of the workpiece, the failure area will appear and expand. Figure

Conclusions
The development of stress relaxation behavior of matrix resin PP was presented in this paper. The DMA was used to measure the stress relaxation behavior of the thermoplastic resin matrix composite material matrix at different temperatures. Master curve and shift function were obtained by time-temperature superposition, and the genetic intelligence algorithm was used to perform regression processing on the data to obtain the relaxation time and displacement factor of the resin matrix.
In this paper, a generalized Maxwell model with N Maxwell elements was used to describe the viscoelastic complexity of the composite stamping process. In addition, according to the mechanical parameters of resin matrix and reinforced fiber, SCFM representative volume element model is used to obtain the mechanical parameters of composites in unrelaxed and completely relaxed states. A three dimensional (3D) incremental viscoelastic constitutive equation is established and implemented into finite element software ABAQUS to predict the deformation and stresses distribution of the parts during hot stamping. The correctness of threedimensional viscoelastic and finite element analysis was verified by using self-designed hot stamping die. The influence of the preheating temperature on the stress distribution of the parts was studied. The preheating temperature has a significant influence on the stress evolution and distribution. The fracture criterion of composite material forming will be further studied to quantitatively study the forming limit of the material.