A Novel CAE Method for Compression Molding Simulation of Carbon Fiber-Reinforced Thermoplastic Composite Sheet Materials

: Its high-speciﬁc strength and stiffness with lower cost make discontinuous ﬁber-reinforced thermoplastic (FRT) materials an ideal choice for lightweight applications in the automotive industry. Compression molding is one of the preferred manufacturing processes for such materials as it offers the opportunity to maintain a longer ﬁber length and higher volume production. In the past, we have demonstrated that compression molding of FRT in bulk form can be simulated by treating melt ﬂow as a continuum using the conservation of mass and momentum equations. However, the compression molding of such materials in sheet form using a similar approach does not work well. The assumption of melt ﬂow as a continuum does not hold for such deformation processes. To address this challenge, we have developed a novel simulation approach. First, the draping of the sheet was simulated as a structural deformation using the explicit ﬁnite element approach. Next, the draped shape was compressed using ﬂuid mechanics equations. The proposed method was veriﬁed by building a physical part and comparing the predicted ﬁber orientation and warpage measurements performed on the physical parts. The developed method and tools are expected to help in expediting the development of FRT parts, which will help achieve lightweight targets in the automotive industry. and mold temperature on mold filling and curing were investigated using finite element analysis (FEA) analysis. In summary, all the previous studies on the simulation of sheet materials for the compression molding process focused on SMC, i.e., discontinuous thermoset materials, continuous fiber laminate structures, or woven fabric with thermoset or thermoplastic resins. The majority of research activities have concentrated only on the thermoforming simulation. Thermoset materials flow well due to low viscosity and therefore were preferred. However, thermoset materials take a longer time to cure and hence the cycle time is longer. Therefore, thermoset materials are preferred for low-volume parts. we in automotive applications that a higher thermoplastic materials that


Introduction
Fiber-reinforced composites offer exciting new possibilities for weight reduction due to their excellent mechanical properties, economical fuel consumption, and lower carbon footprint. For typical automotive components where large volume, low cost, and fast cycle time are desired, chopped reinforcing fibers with thermoplastic resin are increasing in popularity [1][2][3]. Usually, such composite materials are made using injection molding, compression molding, or resin transfer molding processes [4,5]. Of these processes, compression molding offers the highest potential to maintain longer fibers and hence better mechanical properties [6][7][8]; therefore, compression molding is the preferred method by which to manufacture fiber-reinforced thermoplastic components in the automotive industry.
There are typically two different initial formats for the fiber-reinforced charge used in the compression molding process, i.e., (1) bulk charge, as seen in Figure 1, and (2) sheet charge, as per Figure 2. Typically, the bulk charge is made using a single/twin screw low shear plasticator; the desired size of the bulk coming out from the plasticator is collected and placed in the die for compression. Sheet charge manufacturing is more complex. Chopped fibers are dispersed to form mats or spread on the resin sheets and are then partially consolidated. The sheet charge is preferred because it starts with longer fibers and, since it covers most of the mold cavity surface to begin with, the material flow required to fill the mold is shorter; therefore, there is a better possibility of maintaining a longer fiber length.  As the use of fiber-reinforced polymer components in the industry is increasing, interest in using numerical simulation-often called computer-aided engineering (CAE)-to accelerate the development process is also increasing. There have been numerous research activities on the CAE simulation of fiber-reinforced polymeric resin materials in the past 30 years. Early on, Oswald [9,10] and Advani [11,12] developed the basic formulation for the compression molding simulation for fiber-reinforced polymeric materials. As a result of their work, commercial simulation software called Cadpress [13] was developed. Cadpress used simplified Barone-Caulk [14] flow geometry, which was suitable for sheet moulding compound (SMC) (discontinuous fiber + thermoset resin). Wang et al. [15] investigated compression molding simulation of a multilayer composite with continuous fibers and developed a method to estimate the draped shape and temperature distribution of the laminates. Kikuchi [16] studied the compression molding of woven-fabric thermoplastic composite laminates using finite element analysis, and developed an algorithm to optimize temperature distribution to minimize warpage. Three-dimensional thermo-viscoplastic analysis was conducted by Kim et al. [17], where the rheological characteristic of the SMC sheet material was modeled and the effects of dwelling time, and mold temperature on mold filling and curing were investigated using finite element analysis (FEA) analysis.
In summary, all the previous studies on the simulation of sheet materials for the compression molding process focused on SMC, i.e., discontinuous thermoset materials, continuous fiber laminate structures, or woven fabric with thermoset or thermoplastic resins. The majority of research activities have concentrated only on the thermoforming simulation. Thermoset materials flow well due to low viscosity and therefore were preferred. However, thermoset materials take a longer time to cure and hence the cycle time is longer. Therefore, thermoset materials are preferred for low-volume parts. Since we are interested in automotive applications that require a higher volume, thermoplastic materials that do not require curing and can be formed in a very short cycle time are desired. The is the preferred method by which to manufacture fiber-reinforced thermoplastic components in the automotive industry. There are typically two different initial formats for the fiber-reinforced charge used in the compression molding process, i.e., (1) bulk charge, as seen in Figure 1, and (2) sheet charge, as per Figure 2. Typically, the bulk charge is made using a single/twin screw low shear plasticator; the desired size of the bulk coming out from the plasticator is collected and placed in the die for compression. Sheet charge manufacturing is more complex. Chopped fibers are dispersed to form mats or spread on the resin sheets and are then partially consolidated. The sheet charge is preferred because it starts with longer fibers and, since it covers most of the mold cavity surface to begin with, the material flow required to fill the mold is shorter; therefore, there is a better possibility of maintaining a longer fiber length.  As the use of fiber-reinforced polymer components in the industry is increasing, interest in using numerical simulation-often called computer-aided engineering (CAE)-to accelerate the development process is also increasing. There have been numerous research activities on the CAE simulation of fiber-reinforced polymeric resin materials in the past 30 years. Early on, Oswald [9,10] and Advani [11,12] developed the basic formulation for the compression molding simulation for fiber-reinforced polymeric materials. As a result of their work, commercial simulation software called Cadpress [13] was developed. Cadpress used simplified  flow geometry, which was suitable for sheet moulding compound (SMC) (discontinuous fiber + thermoset resin). Wang et al. [15] investigated compression molding simulation of a multilayer composite with continuous fibers and developed a method to estimate the draped shape and temperature distribution of the laminates. Kikuchi [16] studied the compression molding of woven-fabric thermoplastic composite laminates using finite element analysis, and developed an algorithm to optimize temperature distribution to minimize warpage. Three-dimensional thermo-viscoplastic analysis was conducted by Kim et al. [17], where the rheological characteristic of the SMC sheet material was modeled and the effects of dwelling time, and mold temperature on mold filling and curing were investigated using finite element analysis (FEA) analysis.
In summary, all the previous studies on the simulation of sheet materials for the compression molding process focused on SMC, i.e., discontinuous thermoset materials, continuous fiber laminate structures, or woven fabric with thermoset or thermoplastic resins. The majority of research activities have concentrated only on the thermoforming simulation. Thermoset materials flow well due to low viscosity and therefore were preferred. However, thermoset materials take a longer time to cure and hence the cycle time is longer. Therefore, thermoset materials are preferred for low-volume parts. Since we are interested in automotive applications that require a higher volume, thermoplastic materials that do not require curing and can be formed in a very short cycle time are desired. The As the use of fiber-reinforced polymer components in the industry is increasing, interest in using numerical simulation-often called computer-aided engineering (CAE)-to accelerate the development process is also increasing. There have been numerous research activities on the CAE simulation of fiber-reinforced polymeric resin materials in the past 30 years. Early on, Oswald [9,10] and Advani [11,12] developed the basic formulation for the compression molding simulation for fiber-reinforced polymeric materials. As a result of their work, commercial simulation software called Cadpress [13] was developed. Cadpress used simplified  flow geometry, which was suitable for sheet moulding compound (SMC) (discontinuous fiber + thermoset resin). Wang et al. [15] investigated compression molding simulation of a multilayer composite with continuous fibers and developed a method to estimate the draped shape and temperature distribution of the laminates. Kikuchi [16] studied the compression molding of woven-fabric thermoplastic composite laminates using finite element analysis, and developed an algorithm to optimize temperature distribution to minimize warpage. Three-dimensional thermo-viscoplastic analysis was conducted by Kim et al. [17], where the rheological characteristic of the SMC sheet material was modeled and the effects of dwelling time, and mold temperature on mold filling and curing were investigated using finite element analysis (FEA) analysis.
In summary, all the previous studies on the simulation of sheet materials for the compression molding process focused on SMC, i.e., discontinuous thermoset materials, continuous fiber laminate structures, or woven fabric with thermoset or thermoplastic resins. The majority of research activities have concentrated only on the thermoforming simulation. Thermoset materials flow well due to low viscosity and therefore were preferred. However, thermoset materials take a longer time to cure and hence the cycle time is longer. Therefore, thermoset materials are preferred for low-volume parts. Since we are interested in automotive applications that require a higher volume, thermoplastic materials that do not require curing and can be formed in a very short cycle time are desired. The molten thermoplastic materials are highly viscous and therefore the flow geometry in filling the mold cavity is different when compared to thermoset materials. Consequently, the fiber orientation and fiber lengths are also different. A major challenge in CAE simulation of the discontinuous long fiber thermoplastic parts is the ability to simulate compression molding to predict the fiber orientation and warpage of the final part. For continuous fibers, the fiber orientation can be estimated by using the draping process simulation. Recently, CAE methods to simulate the compression molding process of fiber-reinforced thermoplastic composites using the bulk charge format have been developed and validated with physical tests by Song et al. [18]. However, the compression molding simulation approach developed for bulk materials does not work for the compression molding of sheet materials. This is because, when heated fiber-reinforced thermoplastic sheets are placed in the mold cavity, the sheets have very little stiffness and hence drape easily in the mold cavity under very little compressive force. This draping process involves significant rigid body motion of the sheets. The fluid mechanics-based approach used for compression molding of bulk material [19,20] is therefore not adequate to simulate the large rigid body motions that occur during the draping process.
In this paper, a novel CAE method was proposed to address this challenge. The compression molding of sheet materials was simulated in two steps. First, the draping of the heated sheet was simulated using an explicit finite element approach. Such an explicit finite element approach, used for metal sheet stamping simulation, has been proven to be effective for large structural deformation as well as rigid body motions. Next, the draped shape of the fiber-reinforced thermoplastic sheet was used as a prepreg in the compression molding. This is similar to the compression molding of a bulk charge with a given shape, which is an already developed method in our previous work [18]. Furthermore, the proposed approach was demonstrated on a complex 3D shape part. The compression molding simulation was carried out using the proposed two-step approach and the fiber orientation and warpage of the finished part were estimated. To verify the proposed approach, actual parts were made and the measured warpage and fiber orientation results were compared with the predictions. Finally, the effect of the draping distance for the compression molding process was discussed and the ideal draping stroke for the simulation of the given part were estimated. Material properties required to support the simulation, i.e., high temperature stiffness and strength for the draping simulation and material flows and other temperature properties, were measured experimentally.

Sheet Material Compression Molding Issues
Detailed steps for compression molding of sheet materials are presented in Figure 3. There are four main steps.
Step 1: Sheet material is heated above its melt temperature at 270 • C.
Step 2: The heated sheet material is draped in the heated mold cavity, which includes a slight compression in the mold.
Step 3: Compressing the sheet in the mold and holding at high pressure until solidification.
Step 4: Solidified finished part is ejected from the mold and cooled.
From the literature review provided in the first section, it was observed that most commercial software can undertake draping analysis for a continuous unidirectional fiber-reinforced thermoset composite, mainly for wrinkling and thickness analysis. For the compression of discontinuous thermoplastic composite materials, a method to simulate the mold filling process, which can estimate fiber orientation and fiber length, is not available. This is mainly because fiber-reinforced thermoplastic sheets need to be draped onto the mold cavity and then compressed; a single computational tool to address both in the same step is not available. This is because of elastic-plastic structural behavior during draping and the flow behavior of resin melt during compression; the integration of structural and fluid flow simulation for a continuum in one physical system are difficult. Therefore, we propose a two-step approach, as shown in Figure 4. The entire compression molding of sheet material is divided into two steps. In the first step, the draping analysis is conducted using stamping capability using an explicit finite element analysis (LS-DYNA V971)Livermore Software Technology Corporation, Livermore, CA, USA). In the second step, the draped part is treated as a prepreg for the compression molding process using mold fill analysis (Moldex3D, R14) CoreTech System Co., Ltd., Hsinchu County, Taiwan). A translator to convert the deformed shape from LS-DYNA to the initial shape for compression molding was developed.
predictions. Finally, the effect of the draping distance for the compression molding process was discussed and the ideal draping stroke for the simulation of the given part were estimated. Material properties required to support the simulation, i.e., high temperature stiffness and strength for the draping simulation and material flows and other temperature properties, were measured experimentally.

Sheet Material Compression Molding Issues
Detailed steps for compression molding of sheet materials are presented in Figure 3. There are four main steps.
Step 1: Sheet material is heated above its melt temperature at 270 °C.
Step 2: The heated sheet material is draped in the heated mold cavity, which includes a slight compression in the mold.
Step 3: Compressing the sheet in the mold and holding at high pressure until solidification.
Step 4: Solidified finished part is ejected from the mold and cooled. From the literature review provided in the first section, it was observed that most commercial software can undertake draping analysis for a continuous unidirectional fiber-reinforced thermoset composite, mainly for wrinkling and thickness analysis. For the compression of discontinuous thermoplastic composite materials, a method to simulate the mold filling process, which can estimate fiber orientation and fiber length, is not available. This is mainly because fiber-reinforced thermoplastic sheets need to be draped onto the mold cavity and then compressed; a single computational tool to address both in the same step is not available. This is because of elastic-plastic structural behavior during draping and the flow behavior of resin melt during compression; the integration of structural and fluid flow simulation for a continuum in one physical system are difficult. Therefore, we propose a two-step approach, as shown in Figure 4. The entire compression molding of sheet material is divided into two steps. In the first step, the draping analysis is conducted using stamping capability using an explicit finite element analysis (LS-DYNA V971)Livermore Software Technology Corporation, Livermore, CA, USA). In the second step, the draped part is treated as a prepreg for the compression molding process using mold fill analysis (Moldex3D, R14) CoreTech System Co., Ltd., Hsinchu County, Taiwan). A translator to convert the deformed shape from LS-DYNA to the initial shape for compression molding was developed.

Sheet Material Properties Measured for Draping Analysis
For accurate draping analysis, temperature-dependent material properties such as elastic modulus and Poisson's ratio are needed for the sheets. The sheet materials are made from a water slurry process, with 35% carbon fiber. The initial resin was made of polyamide 6 (PA6) fibers. Figure  5a shows the procedures used to measure the mechanical properties of the sheet material. As shown, the two layers of fluffy sheet materials are stacked together under a hot press and consolidated to form the composite plate. The tensile bar samples are cut from the consolidated plate. Tests were performed following the ASTMD638-02 standard [21], and the samples were tested in a temperature chamber to obtain the mechanical properties at different temperatures. The modulus of the material was obtained from the tensile test. The Poisson's ratio was also measured as per the ASTMD638 standard [21] using the optical imaging system. Table 1 shows the measured mechanical properties. The increase of the Poisson's ratio in the end is due to the polymeric phase change at around 180 °C, and crystallization happens again. The modulus and Poisson's ratio results are the average of three repeated tests.

Sheet Material Properties Measured for Draping Analysis
For accurate draping analysis, temperature-dependent material properties such as elastic modulus and Poisson's ratio are needed for the sheets. The sheet materials are made from a water slurry process, with 35% carbon fiber. The initial resin was made of polyamide 6 (PA6) fibers. Figure 5a shows the procedures used to measure the mechanical properties of the sheet material. As shown, the two layers of fluffy sheet materials are stacked together under a hot press and consolidated to form the composite plate. The tensile bar samples are cut from the consolidated plate. Tests were performed following the ASTMD638-02 standard [21], and the samples were tested in a temperature chamber to obtain the mechanical properties at different temperatures. The modulus of the material was obtained from the tensile test. The Poisson's ratio was also measured as per the ASTMD638 standard [21] using the optical imaging system. Table 1 shows the measured mechanical properties. The increase of the Poisson's ratio in the end is due to the polymeric phase change at around 180 • C, and crystallization happens again. The modulus and Poisson's ratio results are the average of three repeated tests. performed following the ASTMD638-02 standard [21], and the samples were tested in a temperature chamber to obtain the mechanical properties at different temperatures. The modulus of the material was obtained from the tensile test. The Poisson's ratio was also measured as per the ASTMD638 standard [21] using the optical imaging system. Table 1 shows the measured mechanical properties. The increase of the Poisson's ratio in the end is due to the polymeric phase change at around 180 °C, and crystallization happens again. The modulus and Poisson's ratio results are the average of three repeated tests.

Sheet Material Properties Measured for Compression Molding Analysis
To accurately simulate the compression molding process, actual material properties such as pressure-volume-temperature (PVT), viscosity, and thermal properties are also needed for this process [22]. The sheet material properties were measured by following the corresponding ASTM standards [23][24][25][26]. After the physical measurement, all the measurement results were converted to the moldex3D material format. These measured composite properties were then imported into the compression molding software Moldex3D's material database and were applied to the process simulation for the part we designed, as Figure 6 shows.  To accurately simulate the compression molding process, actual material properties such as pressure-volume-temperature (PVT), viscosity, and thermal properties are also needed for this process [22]. The sheet material properties were measured by following the corresponding ASTM standards [23][24][25][26]. After the physical measurement, all the measurement results were converted to the moldex3D material format. These measured composite properties were then imported into the compression molding software Moldex3D's material database and were applied to the process simulation for the part we designed, as Figure 6 shows.

Overall CAE Simulation Concept and Steps
Based on the potential needs for the light weighting of automotive components, a three-cavity tool was designed and the sample geometry is presented in Figure 7.

Overall CAE Simulation Concept and Steps
Based on the potential needs for the light weighting of automotive components, a three-cavity tool was designed and the sample geometry is presented in Figure 7. After the initial design of the part, the question was how to use the CAE software to predict the forming process from sheet materials. Another question of interest was how to predict the warpage and microstructural details such as fiber orientation, fiber length, etc. of the part post-compression molding. As proposed in Section 2, the integrated simulation process combining LS-DYNA and Moldex3D was used for this part. Figure 8 shows the overall CAE simulation workflow for the part. The first step was the draping analysis in LS-DYNA. After the draping analysis, the draped part was transferred to the Moldex3D for compression molding process analysis. The details of each simulation process will be discussed separately in the following sections.

Draping Analysis with LS-DYNA
For the first step of the simulation, LS-DYNA was used to predict the draping behavior of the sheet materials. A thermoforming module from LSTC was used in the simulation. The initial sheet material was modeled using solid elements, as is shown in Figure 9. In the model, the total number of solid elements in the model was 15,000, the total number of shell elements in the model was 25,896, and the total number of nodes was 49,280. After the initial design of the part, the question was how to use the CAE software to predict the forming process from sheet materials. Another question of interest was how to predict the warpage and microstructural details such as fiber orientation, fiber length, etc. of the part post-compression molding. As proposed in Section 2, the integrated simulation process combining LS-DYNA and Moldex3D was used for this part. Figure 8 shows the overall CAE simulation workflow for the part. The first step was the draping analysis in LS-DYNA. After the draping analysis, the draped part was transferred to the Moldex3D for compression molding process analysis. The details of each simulation process will be discussed separately in the following sections.

Overall CAE Simulation Concept and Steps
Based on the potential needs for the light weighting of automotive components, a three-cavity tool was designed and the sample geometry is presented in Figure 7. After the initial design of the part, the question was how to use the CAE software to predict the forming process from sheet materials. Another question of interest was how to predict the warpage and microstructural details such as fiber orientation, fiber length, etc. of the part post-compression molding. As proposed in Section 2, the integrated simulation process combining LS-DYNA and Moldex3D was used for this part. Figure 8 shows the overall CAE simulation workflow for the part. The first step was the draping analysis in LS-DYNA. After the draping analysis, the draped part was transferred to the Moldex3D for compression molding process analysis. The details of each simulation process will be discussed separately in the following sections.

Draping Analysis with LS-DYNA
For the first step of the simulation, LS-DYNA was used to predict the draping behavior of the sheet materials. A thermoforming module from LSTC was used in the simulation. The initial sheet material was modeled using solid elements, as is shown in Figure 9. In the model, the total number of solid elements in the model was 15,000, the total number of shell elements in the model was 25,896, and the total number of nodes was 49,280.

Draping Analysis with LS-DYNA
For the first step of the simulation, LS-DYNA was used to predict the draping behavior of the sheet materials. A thermoforming module from LSTC was used in the simulation. The initial sheet material was modeled using solid elements, as is shown in Figure 9. In the model, the total number of solid elements in the model was 15,000, the total number of shell elements in the model was 25,896, and the total number of nodes was 49,280.
The basic model information and boundary condition are shown in Figure 9. During the thermoforming process, thermal properties of the sheet material and the tool were defined, and the mechanical properties of the sheets were also defined. Due to the draping/tearing of the sheet material at the corners, the adaptive remeshing technique was used in the simulation to reduce the instability or any numerical convergence issues [27,28]. Figure 10 shows the displacement contour of the draped part. Then, the shape of the draped part was exported to Moldex3D for the next simulation step. The basic model information and boundary condition are shown in Figure 9. During the thermoforming process, thermal properties of the sheet material and the tool were defined, and the mechanical properties of the sheets were also defined. Due to the draping/tearing of the sheet material at the corners, the adaptive remeshing technique was used in the simulation to reduce the instability or any numerical convergence issues [27,28]. Figure 10 shows the displacement contour of the draped part. Then, the shape of the draped part was exported to Moldex3D for the next simulation step.

Compression Molding Analysis with Moldex3D
Based on the cavity design from Figure 6, a compression molding 3D model was built in Moldex3D, as shown in Figure 11. The element type used in the model was tetra and hex mesh; the total number of elements was 5,944,826 and the total number of nodes was 2,942,618. There were 11 elements through thickness. The green surface on the top is the compression surface, the pink area in the middle is the compression zone, and the gray area on the bottom is the cavity. Figure 12 shows the draped part as the charge layover the cavity. The red area on the top is assigned as the charge; the yellow area on the bottom is assigned as the cavity.  The basic model information and boundary condition are shown in Figure 9. During the thermoforming process, thermal properties of the sheet material and the tool were defined, and the mechanical properties of the sheets were also defined. Due to the draping/tearing of the sheet material at the corners, the adaptive remeshing technique was used in the simulation to reduce the instability or any numerical convergence issues [27,28]. Figure 10 shows the displacement contour of the draped part. Then, the shape of the draped part was exported to Moldex3D for the next simulation step.

Compression Molding Analysis with Moldex3D
Based on the cavity design from Figure 6, a compression molding 3D model was built in Moldex3D, as shown in Figure 11. The element type used in the model was tetra and hex mesh; the total number of elements was 5,944,826 and the total number of nodes was 2,942,618. There were 11 elements through thickness. The green surface on the top is the compression surface, the pink area in the middle is the compression zone, and the gray area on the bottom is the cavity. Figure 12 shows the draped part as the charge layover the cavity. The red area on the top is assigned as the charge; the yellow area on the bottom is assigned as the cavity.

Compression Molding Analysis with Moldex3D
Based on the cavity design from Figure 6, a compression molding 3D model was built in Moldex3D, as shown in Figure 11. The element type used in the model was tetra and hex mesh; the total number of elements was 5,944,826 and the total number of nodes was 2,942,618. There were 11 elements through thickness. The green surface on the top is the compression surface, the pink area in the middle is the compression zone, and the gray area on the bottom is the cavity. Figure 12 shows the draped part as the charge layover the cavity. The red area on the top is assigned as the charge; the yellow area on the bottom is assigned as the cavity. The basic model information and boundary condition are shown in Figure 9. During the thermoforming process, thermal properties of the sheet material and the tool were defined, and the mechanical properties of the sheets were also defined. Due to the draping/tearing of the sheet material at the corners, the adaptive remeshing technique was used in the simulation to reduce the instability or any numerical convergence issues [27,28]. Figure 10 shows the displacement contour of the draped part. Then, the shape of the draped part was exported to Moldex3D for the next simulation step.

Compression Molding Analysis with Moldex3D
Based on the cavity design from Figure 6, a compression molding 3D model was built in Moldex3D, as shown in Figure 11. The element type used in the model was tetra and hex mesh; the total number of elements was 5,944,826 and the total number of nodes was 2,942,618. There were 11 elements through thickness. The green surface on the top is the compression surface, the pink area in the middle is the compression zone, and the gray area on the bottom is the cavity. Figure 12 shows the draped part as the charge layover the cavity. The red area on the top is assigned as the charge; the yellow area on the bottom is assigned as the cavity.   The temperature information and deformed shape were taken into account in the compression molding simulation in Moldex3D. The governing equations of the fluid mechanics that describe the transient and non-isothermal compression flow motion are as follows: where ρ is the density; u is the velocity vector; t is the time; σ is the total stress tensor; g is the acceleration vector of gravity; p is the pressure; η is the viscosity; Cp is the specific heat; T is the temperature; k is the thermal conductivity; and γ  is the shear rate. Tait's pVT model is to express a thermodynamic state relationship, where the volume of the material is a function of the temperature/pressure. The power law model is used to describe complex viscosity behaviors, including the dramatic viscosity change in a lower shear rate range. Details of the numerical implementation are available elsewhere [29] and so are not repeated here. After the part is ejected from the mold, a free thermal shrinkage occurs due to the temperature difference. The equilibrium equation representing the stress is: The relationship between stress and strain is: where σ is the stress; C is a fourth-order tensor and a function of the relaxation modulus E; αCLTE is the coefficient of the linear thermal expansion tensor; ε is the strain tensor, which is transferred by the volume shrinkage during molding; and U is the displacement vector, respectively. In order to consider the viscoelastic behavior after molding, a master function for the relaxation modulus is described using the generalized maxwell model, which is a function of time and temperature with temperature shift factor T a described below: where C1, C2 are constant; and Tref is the reference temperature. The compression molding process simulations were used in Moldex3D, version R14 [13]. The iARD-RPR (the improved anisotropic rotary diffusion model combined with the retarding principal rate model) was used in the simulation with Ci = 0.01. The other parameters used in the iARD-RPR model were Cm = 0.005, Ci = 0.1, and the fiber matrix interaction alpha factor was 0.7. The processing The temperature information and deformed shape were taken into account in the compression molding simulation in Moldex3D. The governing equations of the fluid mechanics that describe the transient and non-isothermal compression flow motion are as follows: where ρ is the density; u is the velocity vector; t is the time; σ is the total stress tensor; g is the acceleration vector of gravity; p is the pressure; η is the viscosity; C p is the specific heat; T is the temperature; k is the thermal conductivity; and . γ is the shear rate. Tait's pVT model is to express a thermodynamic state relationship, where the volume of the material is a function of the temperature/pressure. The power law model is used to describe complex viscosity behaviors, including the dramatic viscosity change in a lower shear rate range. Details of the numerical implementation are available elsewhere [29] and so are not repeated here. After the part is ejected from the mold, a free thermal shrinkage occurs due to the temperature difference. The equilibrium equation representing the stress is: The relationship between stress and strain is: where σ is the stress; C is a fourth-order tensor and a function of the relaxation modulus E; α CLTE is the coefficient of the linear thermal expansion tensor; ε is the strain tensor, which is transferred by the volume shrinkage during molding; and U is the displacement vector, respectively. In order to consider the viscoelastic behavior after molding, a master function for the relaxation modulus is described using the generalized maxwell model, which is a function of time and temperature with temperature shift factor a T described below: where C 1 , C 2 are constant; and T ref is the reference temperature.
The compression molding process simulations were used in Moldex3D, version R14 [13]. The iARD-RPR (the improved anisotropic rotary diffusion model combined with the retarding principal rate model) was used in the simulation with Ci = 0.01. The other parameters used in the iARD-RPR model were Cm = 0.005, Ci = 0.1, and the fiber matrix interaction alpha factor was 0.7. The processing conditions in CAE are the same as the one used to make the physical part. The material properties used in the simulation are listed in Table 2 and Figure 6.

Carbon Fiber Sheet Material Manufacturing
The first step in the production of fiber-reinforced thermoplastic composites is the combining of the fiber reinforcement with the matrix resin to form the sheet material. One technique that has been demonstrated to produce a high-quality self-supporting preform is the wet-lay paper making process [30][31][32]. In this process, the reinforcing fibers and thermoplastic polymer fibers are dispersed in water along with bonding agents. This creates a slurry that is cast onto a moving forming wire and then dewatered, leaving a nonwoven mat of thermoplastic and reinforcing fibers. The carbon fibers, which are made from recycled carbon fiber composites, were supplied by Carbon Conversion (Lake City, SC, USA). The carbon fiber mats have a random fiber orientation, and an average fiber length of about 12 mm, and diameter about 6-7 µm. The resin was PA6 fibers with an average length of about 10 mm, and diameter of 20 µm. Each sheet with a thickness of about 0.2 mm was stacked to form a 2 mm nominal thickness sheet. The mats were fluffy and had in-plane random fiber orientation. The size of the mats we used in this study were 355 mm × 254 mm. Figure 13 shows the steps used in the water slurry process to form the sheet materials. conditions in CAE are the same as the one used to make the physical part. The material properties used in the simulation are listed in Table 2 and Figure 6.

Carbon Fiber Sheet Material Manufacturing
The first step in the production of fiber-reinforced thermoplastic composites is the combining of the fiber reinforcement with the matrix resin to form the sheet material. One technique that has been demonstrated to produce a high-quality self-supporting preform is the wet-lay paper making process [30][31][32]. In this process, the reinforcing fibers and thermoplastic polymer fibers are dispersed in water along with bonding agents. This creates a slurry that is cast onto a moving forming wire and then dewatered, leaving a nonwoven mat of thermoplastic and reinforcing fibers. The carbon fibers, which are made from recycled carbon fiber composites, were supplied by Carbon Conversion (Lake City, SC, USA). The carbon fiber mats have a random fiber orientation, and an average fiber length of about 12 mm, and diameter about 6-7 µm. The resin was PA6 fibers with an average length of about 10 mm, and diameter of 20 µm. Each sheet with a thickness of about 0.2 mm was stacked to form a 2 mm nominal thickness sheet. The mats were fluffy and had in-plane random fiber orientation. The size of the mats we used in this study were 355 mm × 254 mm. Figure 13 shows the steps used in the water slurry process to form the sheet materials.

Actual Part Manufacturing
Three cavity parts were built using carbon fiber (CF) (35%) + PA6 sheet material. The basic procedure for the compression molding of the part is shown in Figure 14. After the sheet material was prepared from the water slurry process, it was moved to the infrared oven, and heated to 343 °C before it was moved to the press.

Actual Part Manufacturing
Three cavity parts were built using carbon fiber (CF) (35%) + PA6 sheet material. The basic procedure for the compression molding of the part is shown in Figure 14. After the sheet material was prepared from the water slurry process, it was moved to the infrared oven, and heated to 343 • C before it was moved to the press.  Table 3 shows the processing conditions of the manufacturing process. After the sample is demolded and cleaned, it is ready for warpage and fiber orientation measurement. Figure 15 shows the warpage measurement method used for the three-cavity tool. The FARO-Arm (FARO Technologies, Inc., Lake Mary, FL, USA) [33,34], which is a state-of-the-art non-contact scanning device for warpage measurement, was used to collect the point cloud of the final part.  The scanned point cloud was then superimposed with the original computer-aided design (CAD) data of the part, as shown. From the superimposed contour, the warpage information was obtained. Computed tomography (CT) scanning and volume graphics (VGSTUDIO MAX, Volume Graphics GmbH, Heidelberg, Germany) analysis [13,35] were used for measuring the fiber orientation. For fiber orientation tensor analysis, the part was cut from the molded part at locations A, B, C, etc., as shown in Figure 16. The sample size used was 10 mm × 10 mm × 2 mm, and the resolution used for the CT scanning was 1 µm.   Table 3 shows the processing conditions of the manufacturing process. After the sample is demolded and cleaned, it is ready for warpage and fiber orientation measurement. Figure 15 shows the warpage measurement method used for the three-cavity tool. The FARO-Arm (FARO Technologies, Inc., Lake Mary, FL, USA) [33,34], which is a state-of-the-art non-contact scanning device for warpage measurement, was used to collect the point cloud of the final part.
The scanned point cloud was then superimposed with the original computer-aided design (CAD) data of the part, as shown. From the superimposed contour, the warpage information was obtained. Computed tomography (CT) scanning and volume graphics (VGSTUDIO MAX, Volume Graphics GmbH, Heidelberg, Germany) analysis [13,35] were used for measuring the fiber orientation. For fiber orientation tensor analysis, the part was cut from the molded part at locations A, B, C, etc., as shown in Figure 16. The sample size used was 10 mm × 10 mm × 2 mm, and the resolution used for the CT scanning was 1 µm.  Figure 14. Three-cavity tool and the part made using sheet material. Table 3 shows the processing conditions of the manufacturing process. After the sample is demolded and cleaned, it is ready for warpage and fiber orientation measurement. Figure 15 shows the warpage measurement method used for the three-cavity tool. The FARO-Arm (FARO Technologies, Inc., Lake Mary, FL, USA) [33,34], which is a state-of-the-art non-contact scanning device for warpage measurement, was used to collect the point cloud of the final part.  The scanned point cloud was then superimposed with the original computer-aided design (CAD) data of the part, as shown. From the superimposed contour, the warpage information was obtained. Computed tomography (CT) scanning and volume graphics (VGSTUDIO MAX, Volume Graphics GmbH, Heidelberg, Germany) analysis [13,35] were used for measuring the fiber orientation. For fiber orientation tensor analysis, the part was cut from the molded part at locations A, B, C, etc., as shown in Figure 16. The sample size used was 10 mm × 10 mm × 2 mm, and the resolution used for the CT scanning was 1 µm.   Table 3 shows the processing conditions of the manufacturing process. After the sample is demolded and cleaned, it is ready for warpage and fiber orientation measurement. Figure 15 shows the warpage measurement method used for the three-cavity tool. The FARO-Arm (FARO Technologies, Inc., Lake Mary, FL, USA) [33,34], which is a state-of-the-art non-contact scanning device for warpage measurement, was used to collect the point cloud of the final part.  The scanned point cloud was then superimposed with the original computer-aided design (CAD) data of the part, as shown. From the superimposed contour, the warpage information was obtained. Computed tomography (CT) scanning and volume graphics (VGSTUDIO MAX, Volume Graphics GmbH, Heidelberg, Germany) analysis [13,35] were used for measuring the fiber orientation. For fiber orientation tensor analysis, the part was cut from the molded part at locations A, B, C, etc., as shown in Figure 16. The sample size used was 10 mm × 10 mm × 2 mm, and the resolution used for the CT scanning was 1 µm.   Figure 17 shows the warpage measurement location from both the measurement and CAE predication. The location of the point was the same for both cases except that the notation used for one was B and the other one was A. After the FARO-ARM laser scanning, the scanned data were superimposed with the original CAD design and the post-processing analysis was conducted using Polywork viewer (2016, Innovmetric Software Inc., Québec, QC, Canada) [36]. At location B, the green line was the warped condition, while the gray line was the initial geometry. Similarly, at location A, the gray line was the original geometry, and the colored contour was the warped shape. The warpage comparison between the CAE prediction at location A and measurement at location B is shown in Figure 17a. The effect of mold temperature to the warpage is shown in Figure 17b. It can be seen that with an increase in the mold temperature, the warpage increased in proportion. The overall error between the simulation and measurement was noted to be around 8%.  Figure 17 shows the warpage measurement location from both the measurement and CAE predication. The location of the point was the same for both cases except that the notation used for one was B and the other one was A. After the FARO-ARM laser scanning, the scanned data were superimposed with the original CAD design and the post-processing analysis was conducted using Polywork viewer (2016, Innovmetric Software Inc., Québec, QC, Canada) [36]. At location B, the green line was the warped condition, while the gray line was the initial geometry. Similarly, at location A, the gray line was the original geometry, and the colored contour was the warped shape. The warpage comparison between the CAE prediction at location A and measurement at location B is shown in Figure 17a. The effect of mold temperature to the warpage is shown in Figure 17b. It can be seen that with an increase in the mold temperature, the warpage increased in proportion. The overall error between the simulation and measurement was noted to be around 8%.

Fiber Orientation Comparison
The key to successful simulation is the ability to simulate the physics of the event. For the fiberreinforced analysis, the fiber orientation plays a key role in structural properties of the part and hence the warpage. Figure 18 shows the fiber orientation measurement results for location A as an overview. It is clearly shown that all the tensors added up to 1.

Fiber Orientation Comparison
The key to successful simulation is the ability to simulate the physics of the event. For the fiber-reinforced analysis, the fiber orientation plays a key role in structural properties of the part and hence the warpage. Figure 18 shows the fiber orientation measurement results for location A as an overview. It is clearly shown that all the tensors added up to 1.  Figure 17 shows the warpage measurement location from both the measurement and CAE predication. The location of the point was the same for both cases except that the notation used for one was B and the other one was A. After the FARO-ARM laser scanning, the scanned data were superimposed with the original CAD design and the post-processing analysis was conducted using Polywork viewer (2016, Innovmetric Software Inc., Québec, QC, Canada) [36]. At location B, the green line was the warped condition, while the gray line was the initial geometry. Similarly, at location A, the gray line was the original geometry, and the colored contour was the warped shape. The warpage comparison between the CAE prediction at location A and measurement at location B is shown in Figure 17a. The effect of mold temperature to the warpage is shown in Figure 17b. It can be seen that with an increase in the mold temperature, the warpage increased in proportion. The overall error between the simulation and measurement was noted to be around 8%.

Fiber Orientation Comparison
The key to successful simulation is the ability to simulate the physics of the event. For the fiberreinforced analysis, the fiber orientation plays a key role in structural properties of the part and hence the warpage. Figure 18 shows the fiber orientation measurement results for location A as an overview. It is clearly shown that all the tensors added up to 1.  We also compared the fiber orientation at selected locations such as A, B, and C (see Figure 19). A11 is usually the direction of the fiber along the X direction, as shown in Figure 19. The comparison of the fiber orientation tensor A11 between the prediction and measurement at these locations is also shown. We also compared the fiber orientation at selected locations such as A, B, and C (see Figure 19). A11 is usually the direction of the fiber along the X direction, as shown in Figure 19. The comparison of the fiber orientation tensor A11 between the prediction and measurement at these locations is also shown. According to the Department of Energy (DOE) findings, the discrepancy between the simulation and measurement was within the 15% accuracy criterion [37]. Comparison of the fiber orientation results at other locations such as D, E, and F also showed trends similar to A, B, and C. The comparison results are not presented here for brevity. As a reference, Table 4 shows the analytic computed Young's modulus E11 based on fiber orientation from the predicted vs. measured at locations A, B, and C. Such an analysis was based on the Mori-Tanaka Mean Field homogenization scheme [38]. It is clear that the E11 value showed less than a 10% difference between these two.  Figure 20 shows the fiber length distribution comparison. The fiber lengths of the parts were measured using the FASEP method [39]. In the compression molding process simulation, the fiber breakage calculation is used to obtain the fiber length distribution. The fiber breakage algorithm was According to the Department of Energy (DOE) findings, the discrepancy between the simulation and measurement was within the 15% accuracy criterion [37]. Comparison of the fiber orientation results at other locations such as D, E, and F also showed trends similar to A, B, and C. The comparison results are not presented here for brevity. As a reference, Table 4 shows the analytic computed Young's modulus E11 based on fiber orientation from the predicted vs. measured at locations A, B, and C. Such an analysis was based on the Mori-Tanaka Mean Field homogenization scheme [38]. It is clear that the E11 value showed less than a 10% difference between these two.  Figure 20 shows the fiber length distribution comparison. The fiber lengths of the parts were measured using the FASEP method [39]. In the compression molding process simulation, the fiber breakage calculation is used to obtain the fiber length distribution. The fiber breakage algorithm was developed from the Coretech company (Hsinchu County, Taiwan) and implemented in Moldex3D software, so the algorithm is available anywhere [40]. The comparison shows that the fiber length of the final part is within the range of 3-4 mm. developed from the Coretech company (Hsinchu County, Taiwan) and implemented in Moldex3D software, so the algorithm is available anywhere [40]. The comparison shows that the fiber length of (a) (b)

Stroke Distance Effect
There is one issue/concern in the compression molding of sheet materials using the proposed approach, i.e., when should the draping stop and how does a partially draped part affect the warpage? In production, the compression molding process occurs continually at the established cycle times, typically 2-3 min depending on the part size. To address this issue, the LS-DYNA analysis was conducted considering three different boundary conditions. As shown in Figure 21, the draping distance can be 100% of the gap between the top punch and the bottom die, or a fraction, i.e., 60%, 80%, etc. The draped part is then transferred to moldex3D for compression molding simulation.

Stroke Distance Effect
There is one issue/concern in the compression molding of sheet materials using the proposed approach, i.e., when should the draping stop and how does a partially draped part affect the warpage? In production, the compression molding process occurs continually at the established cycle times, typically 2-3 min depending on the part size. To address this issue, the LS-DYNA analysis was conducted considering three different boundary conditions. As shown in Figure 21, the draping distance can be 100% of the gap between the top punch and the bottom die, or a fraction, i.e., 60%, 80%, etc. The draped part is then transferred to moldex3D for compression molding simulation.
It is shown that the draping distance has an effect on the prepreg shape; therefore, the final warpage of the part is different for each case. From the comparison of warpage measurement results, it was seen that 80% of the stroke distance was best matched with the actual measurement results. When the mold temperature was 140 • C, similar CAE simulations were performed. It is shown that, for 80% of the draping distance, the warpage had matching results between CAE and measurement. approach, i.e., when should the draping stop and how does a partially draped part affect the warpage? In production, the compression molding process occurs continually at the established cycle times, typically 2-3 min depending on the part size. To address this issue, the LS-DYNA analysis was conducted considering three different boundary conditions. As shown in Figure 21, the draping distance can be 100% of the gap between the top punch and the bottom die, or a fraction, i.e., 60%, 80%, etc. The draped part is then transferred to moldex3D for compression molding simulation.

Conclusions
A novel two-step CAE method to simulate the compression molding of chopped carbon fiber-reinforced thermoplastic sheet material was developed. First, the draping simulation was carried out to model the draping of a heated carbon thermoplastic sheet on the mold cavity with very little pressure. Next, the draped part was transferred as a prepreg for a further compression molding process. The two-step approach helps with simulating the complex process whereby a solid elastic-plastic material is thermoformed and then compressed to flow like a highly viscous fluid. The effect of the draping distance to the final warpage was discussed and the ideal draping stroke for the given component was demonstrated. The predicted warpage, fiber length, and fiber orientation for the parts made from sheet material using the proposed approach were compared with the measurements from the physical parts. The comparison showed an acceptable match between the simulation and the actual measurements. In this study, it was assumed that the fibers in the sheet form were of in-plane random orientation, and, after draping, they remained as in-plane randomly oriented fibers. The research findings based on these models will benefit lightweight composite parts/components development for the automotive industry.