Research on Anisotropic Viscoelastic Constitutive Model of Compression Molding for CFRP

The carbon-fiber-reinforced polymer (CFRP) is a mainstream material for lightweight products from the end of the 20th century to the present day. Its compression molding process has obvious advantages in mass production. This paper attempts to establish the constitutive models of compression molding of the CFRP materials and study their mechanism. Based on anisotropic linear elastic mechanics, viscoelastic mechanics, and thermodynamics, as well as the Maxwell viscoelastic constitutive model, we first establish the constitutive model of thermorheologically simple CFRP materials (TSMs). Then, considering the influence of temperature on the initial stiffness and equilibrium stiffness, the concept of temperature stiffness coefficient is introduced, and the Cartier coordinate system is converted into a cylindrical coordinate system, thereby establishing the constitutive model of thermorheologically complex materials (TCMs) using the tensor form. Finally, by comparing to the structure of the Zocher model, the two constitutive models established in this study are verified. The research findings have important theoretical research significance for studying the compression molding mechanism of carbon fiber and further improving the quality of product molding.

the constitutive model of the CFRP molding process should be established. It is of great significance for reasonably predicting the residual stress development of the product molding process, optimizing the molding process, and improving product quality.
Scholars worldwide have conducted a lot of research based on the constitutive model of the CFRP materials and achieved rich results. The linear elastic constitutive equation combined with laminate theory is a commonly used method to predict the residual stress of composite materials. Stango et al. [8] assumed that the material properties during the curing process were constants and adopted the laminate theory to study the temperature-induced composite residual stress. Bogetti et al. [9] believed that the performance of the composite material was related to the curing degree and proposed an improved linear elastic constitutive model, namely the CHILE (α) model. Johnston [10] stated that the performance of the composite material during curing is related to the glass transition temperature Tg of the matrix, and they established the CHILE (t g ) model. In essence, both the CHILE (α) and CHILE (t g ) models are linear elastic constitutive models, which cannot reflect the characteristics of material stress relaxation and creep; the CFRP resin matrix exhibits obvious viscoelastic properties at high temperatures. Thus, the viscoelastic constitutive model can truly describe the mechanical properties of composite materials during curing. Zocher et al. [11] used the generalized Maxwell model and Prony technology to express the relaxation stiffness of composite materials, and they derived the anisotropic constitutive equation and incremental equation. Kim et al. [12,13] assumed that Poisson's ratio of the resin during the curing process was unchanged, and combined the time-temperature-curing degree equivalence principle to obtain the viscoelastic properties of related resins, which provided important references for understanding the viscoelastic behavior of the materials during curing. Zobeiry et al. [14] used a differential type to represent the viscoelastic constitutive equation and its increment. Abouhamzeh et al. [15] combined the Laplace transform and inverse transform and proposed a solution to predicting the viscoelastic behavior of anisotropic composite materials. Svanberg et al. [16] proposed a path-dependent constitutive model and studied the curing and post-curing deformation mechanism of L-shaped parts.
Based on the previous research, this paper establishes the CFRP anisotropic viscoelastic constitutive models of thermorheologically simple CFRP materials (TSMs) and thermorheologically complex materials (TCMs) by combining with anisotropic linear elastic mechanics, viscoelastic theory, and thermodynamics, and the Maxwell viscoelastic model. Considering the effect of material temperature on the stiffness and equilibrium stiffness of the carbon fiber molding process, the temperature-dependent thermoelastic stiffness coefficient was introduced, and the Cartesian coordinate system was converted into a cylindrical coordinate system, thereby establishing the CFRP anisotropic elastic constitutive equation and increment equations of the TCM materials. The establishment of the CFRP anisotropic viscoelastic constitutive models, especially for TCM materials, will lay a theoretical foundation for revealing the mechanism of carbon-fiber-reinforced polymers [17], and then provide theoretical guidance to improve the quality of products molding, to further speed up the industrialization process of CFRP.

Theoretical Basis of CFRP Compression Molding
Generally, if the performance of each point in the object is the same, it is a homogeneous material, otherwise, a heterogeneous material. In addition, materials that show the same properties in every direction at every point in the object are called isotropic materials, while those with different properties at every point in the object are anisotropic materials. For a CFRP material, the matrix is an isotropic material, and the reinforcing fiber is a transversely isotropic body, so the CFRP material is an anisotropic material. In the molding process of CFRP materials, the matrix material has undergone a transition process from the viscous fluid state to rubber state to glass state, etc. It is necessary to combine the elastic mechanics of the composite material and the viscoelastic mechanics in the research.

The basis of Anisotropic Linear Elastic Mechanics
To study the mechanical properties of composite materials and solve the stress and strain of a continuous elastomer under external load, we need to balance the equations, geometric equations, constitutive equations, etc.
In the space coordinate system xyz, the stress state of any point in the elastic body under any load is represented by the normal stress components σ x , σ y , and σ z , and the shear stress components τ xy , τ yz , and τ zx . According to the theorem of conjugate shearing stress, there is τ xy = τ yx , τ yz = τ zy , and τ zx = τ xz . F x , F y , and F z respectively represent the external force components received by the elastic body in the x, y, and z directions. Ignoring the volume force, the equilibrium relationship of any point in the elastic body along the coordinate axes x, y, and z is: Similarly, in the space coordinate system xyz, the strain state of any point on the elastic body can be represented by the normal strain components ε x , ε y , and ε z and the shear strain components γ xy , γ yz , and γ zx at that point. If u, v, and w represent the displacement components in the three directions of x, y, and z, then the geometric equation of any point in the elastic body along the x, y, and z directions is: According to Equations (1) and (2), in the space coordinate system xyz, there are 15 unknown functions, namely 6 stress components, 6 strain components, and 3 displacement components. Their relationship is shown as: Among them, C mn (m, n = 1, 2, 3, 4, 5, 6) is the stiffness coefficient. The matrix of the stiffness coefficient is symmetric, and only 21 stiffness coefficients are independent. Equation (3) establishes the relationship between stress and strain, which is called the generalized Hooke's law or elastic constitutive equation.
For the matrix of a CFRP material, the number of independent stiffness coefficients is 2, and its stress-strain relationship is The reinforced fiber has 5 independent stiffness coefficients. Assuming that the xoy coordinate plane is isotropic, its stress-strain relationship is

The Basis of Anisotropic Viscoelastic Mechanics
The mechanical properties of viscoelastic materials such as shear modulus, loss modulus, and loss factor, etc. are usually related to ambient temperature, vibration frequency, strain amplitude, etc. Therefore, the constitutive relationship of viscoelastic materials is very complicated.
When describing the viscoelastic behavior of materials, the generalized Maxwell model is generally used, which consists of multiple Maxwell elements in parallel with a spring, as shown in Figure The reinforced fiber has 5 independent stiffness coefficients. Assuming that the xoy coordinate plane is isotropic, its stress-strain relationship is ( )

The Basis of Anisotropic Viscoelastic Mechanics
The mechanical properties of viscoelastic materials such as shear modulus, loss modulus, and loss factor, etc. are usually related to ambient temperature, vibration frequency, strain amplitude, etc. Therefore, the constitutive relationship of viscoelastic materials is very complicated.
When describing the viscoelastic behavior of materials, the generalized Maxwell model is generally used, which consists of multiple Maxwell elements in parallel with a spring, as shown in Figure 1: Each Maxwell element contains a spring with elastic modulus E and a viscous strain damper with strain viscosity η . In the generalized Maxwell model, the modulus is used to measure the performance of the spring, and the relaxation time is used to express the performance of the damper. The constitutive equation of anisotropic viscoelastic materials [18][19][20][21] is: where ij C is the stiffness matrix, eff j ε is the effective strain tensor, and i σ , α , and τ are the stress tensor, curing degree, and virtual time integral variable of the tensor, respectively. The effective strain tensor can be expressed by the total strain tensor j ε and the non-mechanical strain tensor tc j ε , namely: Each Maxwell element contains a spring with elastic modulus E and a viscous strain damper with strain viscosity η. In the generalized Maxwell model, the modulus is used to measure the performance of the spring, and the relaxation time is used to express the performance of the damper. The constitutive equation of anisotropic viscoelastic materials [18][19][20][21] is: where C ij is the stiffness matrix, ε e f f j is the effective strain tensor, and σ i , α, and τ are the stress tensor, curing degree, and virtual time integral variable of the tensor, respectively.
The effective strain tensor can be expressed by the total strain tensor ε j and the non-mechanical strain tensor ε tc j , namely: ε If the material exhibits simple thermo-rheological properties when cured, the above equation can be changed to: where ξ is the time and is expressed as: Materials 2020, 13, 2277 5 of 13 Among them, α T is the conversion equation, related to temperature and curing degree; s is the time integral variable.
The stiffness matrix C ij is represented by the Prony series: where C ∞ ij is the stiffness matrix after the material is completely relaxed, that is, the equilibrium stiffness; C m ij is the stiffness matrix of each m branch in the Prony series; τ m is the discrete relaxation time of the m-th branch.

The Basis of Thermodynamics
During the molding process of the composite material, the matrix undergoes a curing reaction to release heat, which is regarded as an internal heat source and controls the entire heat transfer process together with the applied temperature.
According to the Fourier heat conduction equation, the heat transfer process can be expressed as: where ρ, C p , k x , k y , and k z respectively represent the density, specific heat, and anisotropic thermal conductivity of the composite material; T is the temperature at time t; Q is the heat generated by the curing reaction of the resin, which can be calculated as In Equation (13), ρ m is the resin matrix density, V f is the fiber volume fraction, H t is the total heat released by the resin during complete curing, and dα dt is the instant curing rate of the resin. The curing degree can be expressed as The boundary conditions of the heat conduction model are: where T s and T m are the surface temperature and heating temperature of the composite material, respectively; K e f f and h e f f are the equivalent thermal conductivity and equivalent convection heat transfer coefficients of the composite material surface, respectively.

CFRP Constitutive Model of the TSMs
The viscoelastic constitutive modulus can reflect the creep and relaxation characteristics of the material. In this paper, the generalized Maxwell viscoelastic constitutive equation was used to reflect the viscoelastic characteristics of the simple thermal-rheological composite material.

The Constitutive Model Based on Viscoelastic Theory
For a single Maxwell element, the strain of the element is a combination of spring and viscous element strain, which is expressed as Using the relaxation modulus to express the relaxation characteristics of the material, where τ = η/E, i.e., the relaxation time.
For the generalized Maxwell model, the stress relaxation modulus is given as E m and τ m are the elastic modulus and relaxation time of the m-th branch, respectively. It can be seen from Figure 1 that the N + 1-th unit has only one spring, the relaxation modulus of the viscous element is negligible, and its elastic modulus is the modulus E ∞ after complete relaxation. The initial modulus of the generalized Maxwell model is E 0 , and it can then be expressed as where W m is the weight coefficient, which is expressed as The stress of the m-th branch in the generalized Maxwell model at time t is expressed as where σ t−∆t m and σ t m , ε t−∆t e f f , and ε t e f f represent the stress tensor and effective strain tensor before and after the time increment ∆t, respectively, t t m and C m , respectively represent the relaxation time and stiffness matrix of the m-th branch of the Maxwell model at time t. The effective stress tensor can be expressed as Among them, ε s represents the total stress tensor, K cte and K ccs represent the effective thermal expansion coefficient and chemical shrinkage coefficient, respectively, ∆t is the temperature increment, and ∆α is the curing degree increment. Combining the stresses of all branches of the Maxwell model and the balance spring, we obtain the total stress at time t: Based on the above, the stress increment at time t can be expressed as where C ij represents the stiffness matrix of the spring, and the stiffness matrix of the m-th branch can be expressed as C 0 ij is the initial stiffness matrix and C ∞ ij is the stiffness matrix after complete relaxation. Generally, C ∞ ij = rC 0 ij ; r is a constant, and different material systems have different values of r. The conversion equation can be expressed as Among them, α 1 and α 2 are constants, and different matrices have different values.

CFRP Viscoelastic Constitutive Equation and Incremental Equation of the TSMs
When the CFRP materials are TSMs, the material's equilibrium stiffness and initial stiffness are unrelated with the degree of curing. According to Equations (11), (19), and (27), the relaxation stiffness of the CFRPs can be expressed as: For anisotropic materials, the internal stress is expressed as: [13,22,23] Considering the effects of temperature T and curing degree on relaxation stiffness, the above equation can be expressed as:F Considering the stress increment ∆σ ∆t i within the time increment ∆t, there is When the time increment is small, the curing degree is considered to be approximately constant, that is, Therefore, the stress increment can be expressed as: [24] ∆σ where S i m is a historical state variable; the initial value is 0, which can be expressed as:

Thermoelastic Expressions of the TCMs
To simulate the high-temperature viscoelastic behavior of composite materials, it is generally assumed that the composite material is a simple thermal-rheological material, that is, the equilibrium stiffness and initial stiffness of the material are related to temperature and curing degree. Many experiments have shown that the initial stiffness and equilibrium stiffness are related to the material temperature [25]. In this paper, the viscoelastic materials whose equilibrium stiffness and initial stiffness change with temperature are called thermorheologically complex materials.
In the curing process of composite materials, the molding temperature and the curing degree of the matrix change with time, and the elastic modulus of the material is related to the temperature and the curing degree. On the one hand, it affects the conversion factor and relaxation time, similar to the performance of the viscous elements; on the other hand, it also affects C 0 ij and C ∞ ij , which is similar to the elastic properties.
When constructing the viscoelastic model of the material in this paper, the vertical movement coefficient representing the correlation between C ∞ ij and temperature was added to be the thermal viscosity stiffness coefficient. Through the experimental method and data processing method provided by Kim [12,13], the exponential function was used to express the thermoelastic stiffness coefficient of an epoxy resin base, and linearly fit it, as shown in Figure 2 below.

Conversion of Coordinate Systems
For a better description in the study of the CFRP mechanical properties, it is necessary to establish a local coordinate system for the geometric elements or convert the global coordinate system. In this paper, the Cartesian coordinate axis was converted into a cylindrical coordinate system for reducing the amount of calculation.

Conversion of Coordinate Systems
For a better description in the study of the CFRP mechanical properties, it is necessary to establish a local coordinate system for the geometric elements or convert the global coordinate system. In this paper, the Cartesian coordinate axis was converted into a cylindrical coordinate system for reducing the amount of calculation.
In the Cartesian coordinate system (x,y,z), the unit base vectors of the coordinate axes are i , j , and k , respectively, and the unit vectors converted into the cylindrical coordinate system are e r , e ϕ , and e z , respectively. It is shown as Then, the conversion formula in the Cartesian coordinate system and cylindrical coordinate system is: The stress vector in the cylindrical coordinate system can be expressed by the stress vector σ i in the Cartesian coordinate system, namely: where ββ T = 1.

CFRP Viscoelastic Constitutive Equation and Enhancement Equation of the TCMs
As the anisotropic viscoelastic constitutive model of TCMs contains many factors, the tensor was used in this paper to express the three-dimensional viscoelastic constitutive model and incremental equations of TCMs.
For the one-dimensional Maxwell element, the stress relaxation under unit constant strain is equal to the relaxation modulus of the model. In view of the thermoelastic stiffness coefficient, the relaxation modulus of the element at time t is: where d is the thermoelastic stiffness coefficient, d = E 0 (T)/E 0 (T(0)). E 0 (T) is the spring stiffness and E 0 (T(0)) is the spring stiffness at t = 0. For the three-dimensional generalized Maxwell model, the Cartier coordinate system was converted into a cylindrical coordinate system to express the relaxation modulus. At this time, it appeared to be isotropic on the xoy plane. When calculating, the three-dimensional generalized model is equivalent to the two-dimensional generalized model, which can be compared to the relaxation modulus representation of a one-dimensional Maxwell element, no sum on ϕ,k and l Then, the thermoelastic stiffness coefficient is d T(t) ϕkl = C ϕkl (T)/C ϕkl (T(0)). According to Equation (32), the stress value at time t is: According to Equations (44) and (45), the stress σ t ϕ at the reduced time ξ can be rewritten as: The relaxation modulus C t ϕkl is expressed as: Referring to (46), the stress σ t ϕ at the reduced time ξ t is expressed as: In the formula, the reduced time ξ t is discretized for the sum of the reduced time ξ t−∆t at time t and the reduced time increment ∆ξ t at time ∆t, that is: The stress increment equation ∆σ t ϕ at time t is expressed as: Among them, H1 ϕ can be written in recursive form, At very small time increments, the temperature T and the curing degree are extremely small and negligible, and the conversion factor is constant, so it can be expressed as: In addition, Then, Use the same method to deal with H2 ϕ and H3 ϕ . Finally, the stress increment ∆σ t ϕ at time t is expressed as: In addition,

Verification of the Constitutive Models
To verify the constitutive models of both TSM and TCM CFRP materials, this paper selects the molded flat structure of a CFRP material based on Kim's research [12,13]. Four layers were molded in the layering angle of [0/90/90/0]. Compared to the viscoelastic constitutive model proposed by Zocher [11], the curing degree change of the center point (0, 0, 0) of the laminate during the curing process is shown in Figure 3 below:  It can be seen from Figure 3 that the three curves are basically coincident, verifying the constitutive models given in this paper. It can be seen from Figure 3 that the three curves are basically coincident, verifying the constitutive models given in this paper.

Conclusions
This paper studies the constitutive models of the CFRP during the compression molding process. The main content and conclusions are as follows: • On the basis of the constitutive model of CFRP molding and viscoelastic theory, viscoelastic mechanics, and thermodynamics, the CFRP viscoelastic constitutive equation and incremental equation of the TSMs were established, which lays a foundation for explaining the solidification and relaxation behavior of the molding process.

•
Considering the effects of temperature on the initial stiffness and equilibrium stiffness, the thermoelastic stiffness coefficient was introduced, and the Cartesian coordinate system was converted to the cylindrical coordinate system. Thus, the viscoelastic constitutive equation and incremental equation of the TCMs were established, and the accuracy and application range of the constitutive model were improved while the calculation scale was reduced.