A Novel Multi-Region, Multi-Phase, Multi-Component-Mixture Modeling Approach to Predicting the Thermodynamic Behaviour of Thermoplastic Composites during the Consolidation Process

In the processing of thermoplastic composites, great importance is attributed to the consolidation step, as it can significantly reduce the porosity of the semi-finished product and thus influence considerably the quality of the final component. This work presents an approach to modeling the thermodynamic behavior of composite materials during hot-press consolidation. For this purpose a multi-region, multi-phase and multi-component-mixture model was developed using the simulation toolbox OpenFOAM®. The sensitivity of the model was tested by varying the thermal parameters and mesh resolution, confirming its robustness. Validity of the model was confirmed by comparing simulation results to experimental data for (i) polycarbonate with 44% carbon fiber by volume and (ii) polypropylene with 45.3% glass fiber by volume. The simulation allows very precise estimation of when a particular temperature, such as the glass transition temperature or melting point, will be reached at the core of a composite. In relation to the total process time, maximum deviation of the simulation from the experimental data amounted to 2.84%. Therefore, the model is well suited for process optimization, it offers a basis for further model implementations and the creation of a digital twin.


Introduction
Achieving climate targets, such as a reduction in CO 2 emissions, requires efficient solutions in the transportation sector, first and foremost in the aviation and automotive industries. Lighter components with consistently good mechanical properties are used to reduce fuel and energy consumption and thus lower harmful emissions. Composite materials that consist of a polymer matrix with fiber reinforcement have been shown to fulfill these requirements [1][2][3][4].
Modern development of composites gained traction in the 1930s with the introduction of glass fibers by the Owens-Illinois glass company and with patents awarded to Carlton Ellis and Paul Schlack for polyester and epoxy resins, respectively. It was found that combining polymer resin with glass fiber resulted in materials with excellent mechanical properties at low weight; these materials were further improved during World War II, which led to their increased commercial availability and production in the USA [5].
Research into thermoplastic matrix systems began in the 1980s, but only became significant when thermoplastic composites were used in the construction of the Airbus A340-600 in 2002: Glass-fiber reinforced thermoplastics were employed to form most of the During tape laying, the individual tapes are stacked on top of each other. Since UD tapes are highly loadable only in the fiber direction, the fiber orientation is usually varied during laying. For aerospace applications, a layup of [0 • | ± 45 • |90 • ] S is prominent [17]. Tape laying can be fully automated: either pick-and-place, automated tape laying (ATL), or automated fiber placement (AFP). In pick-and-place tape laying, the individual pre-cut tapes are stacked on top of each other by a robot and welded locally, for instance, by hot-stamping or ultrasonic welding. Since the pick-and-place principle involves welding the tapes together locally, subsequent consolidation is required. The consolidation process can be carried out using hydraulic heating and cooling presses, where the layup is first heated under pressure in a heating press beyond the glass transition or melting point of the matrix material and then cooled in a cooling press, again while applying pressure, to a temperature below the glass transition or melting point of the matrix material.
In order to be formable the consolidated part must subsequently be heated to a temperature above the melting point for semi-crystalline polymers and above the glass transition temperature for amorphous polymers. Infrared or convection ovens are usually used in the preheating process [18], during which deconsolidation may occur depending on the properties of the product, such as the degree of porosity, moisture content, fiber network crosslinkage, degree of crystallization and matrix viscosity [19,20].
The forming process is often carried out by hydraulic presses, in which a semi-finished product is held between the mold halves and formed when they are brought together. This process can also be carried out in an injection molding machine, which allows simultaneous overmolding and thus functionalization.
As mentioned above, consolidation of sufficient quality can be achieved without a specific consolidation step if optimal process parameters are used for in-situ consolidation during the ATL or AFP process. However, a consolidation step by using a pressing operation before pre-heating and forming increases quality in terms of porosity and interlaminar shear strength (ILSS) [21].
As described in [22], ATL and AFP are mostly used during the production of thermplastic composite parts. Hence, lastest literature about the consolidation process and its modeling is related to ATL and AFP, whereas hardly any literature can be found on the consolidation process using a hot press. This work focuses exclusively on the consolidation by hot press. To predict the thermodynamic behavior of the tape stack during consolidation, we formulated and solved numerically a new mathematical model, using the open-source CFD (Computational Fluid Dynamics) Software OpenFOAM ® . To this end, we generated a new solver that represents a multi-region, multi-phase and multi-component-mixture flow of a compressible fluid under non-isothermal, transient conditions.

Basic Equations
To predict the consolidation process of composite materials, a mathematical model was developed and solved numerically in OpenFOAM ® . Figure 3 shows a schematic of the model setup, which comprises a solid and a fluid region. The former includes the heating/cooling plates, which heat/cool and exert pressure to the composite, and the tools needed for transportation between the heating and cooling press. The fluid region consists of two phases: (i) the composite material and (ii) the air that surrounds it laterally. The composite material is represented by a homogeneous multi-component-mixture model that comprises polymer matrix and fiber fraction, and ignores the individual tape layers. Thus, a model of a compressible fluid flow under transient, non-isothermal conditions is obtained.
alpha.compositea Division of the computational grid into solid and fluid domains means that different assumptions are made for the respective regions and different equations are solved. For the solid domains, including the heating/cooling plates and the tools, the energy conservation equation is solved based on pure heat conduction.
Equation (1) describes the change in specific enthalpy h over time as a function of thermal diffusivity α th where λ is the heat conductivity, ρ the density, and c p the specific heat capacity at constant pressure. The heating and cooling plates and the tools are made of steel, for which the relevant values were obtained from the literature: λ: 54 W/(m*K), ρ: 7854 kg/m 3 , c p : 461 J/K. The transport of the phases in the fluid area (composite and air) is generally described by the conservation equations of mass, momentum, and energy in combination with the Volume of Fluid (VOF) model (see Equations (3), (4), (5) and (8)) where α is a dimensionless parameter that indicates, whether a cell contains composite (α = 1) or air (α = 0) (see Figure 3). The last term corrects for the smearing of the two immiscible phases (i.e., 0 < α < 1), with u r directed against the flow [23]. Mass, momentum and energy conservation (Equations (4), (5) and (8)) are solved for both phases composite and air in combination with Equation (3), considering the respective material parameters: with F S representing surface forces: where S vis is the viscous stress tensor, and p the pressure. Outer body forces F B , such as gravity, are ignored in this work. Further, to express the viscous stress tensor (Equation (7)), Newtonian behavior is assumed, from which follows that the dynamic viscosity η M of the mixture of fiber and matrix is independent of the shear rate. However, the viscosity is assumed to be a polynomial function of the temperature T fitted to experimentally obtained data and Equation (9), and is calculated from: where η is the dynamic viscosity and I the identity tensor. To include non-isothermal effects, the conservation of the inner energy e is considered by: The first three terms describe the change in inner energy e with time, convection of the inner energy e and compression heating. The remaining terms represent the change in mechanical energy K with time and convection of the mechanical energy K, respectively. The terms on the right-hand side describe the heat conduction, with α th,e f f being the effective thermal diffusivity. The optional source term S allows crystallization or melting energy to be included. However, since the focus was on polycarbonate, an amorphous matrix material, this effect was omitted in this work.
The material properties of the mixture of matrix and fibers, that is, density ρ M , thermal conductivity λ M , specific heat capacity c p,M , and dynamic viscosity η M , are calculated using the rule of mixture (Equation (9)) [24]: Here, indices m and f refer to the matrix and the fiber, respectively. Further, f v f is the fiber volume fraction with V as the volume: The density of the matrix ρ m was measured by a high capillary rheometer (HKR Rheograph 25, GÖTTFERT Werkstoff-Prüfmaschinen GmbH, Buchen, Germany). The thermal conductivities of the matrix λ m and the fibers λ f were taken from literature. The specific heat capacity of the matrix c p,m was determined by differential scanning calorimetry with a Mettler Toledo DSC 1 (Mettler Toledo Group, Columbus, OH, USA), and the dynamic viscosities of the matrix η m and the mixture η M were evaluated with an ANTON-Paar MCR 302 plate-plate rheometer (Anton Paar, Graz, Austria). The sources of the data are listed in Table 1. In order to study the effect of varying thermal conductivity and specific heat capacity on the simulation results, a sensitivity study was conducted (see Section 2.2.1). The multimixture model describes the composite as a homogeneous material and ignores any local irregularities, such as fiber accumulations. Furthermore, anisotropic thermal and flow behaviors, which heavily depend on the fiber orientation, were ignored.

Boundary Conditions
The boundary conditions have been chosen to the best of our knowlegde to match the reality as exactly as possible. To model the heat transfer between heating/cooling plates, tools and composite, a boundary condition for the heat conduction is used. To this end, a value fraction v f , is calculated at the interface of two regions (solid/solid or solid/fluid), and used to determine the wall temperature T w : The indices F and S refer to the fluid and solid regions, respectively, and d describes the height of a finite volume element at the corresponding boundary wall.
An impeded heat transfer due to surface roughness is considered by the degree of intimate contact D IC in Equation (11), which is based on the work presented in [25,26], and described in [27][28][29][30][31][32] and calculated by: Intimate contact is based on a simplified view of the surface roughness, which is assumed to be approximately rectangular and described by the initial geometrical values w 0 , b 0 and a 0 , the applied pressure P app and the temperature-dependent dynamic zeroviscosity η 0 (T). Figure 4 shows a schematic of the contact region. Assuming that there is no ideal contact between the heating/cooling plates and the tools, which limits heat conduction, a thin layer of air is considered in the boundary condition, which yields: where the index A refers to the air layer.
To represent the wall adhesion of the composite material to the tools, a noslip boundary condition is chosen at the interface to exclude velocity along the boundary face. A constant value is specified for the pressure at the boundary to the lower tool; for all other interfaces the pressure is calculated from the velocity.

Sensitivity Study
The specific heat capacity of the matrix material c p,m was measured by differential scanning calorimetry (Equipment: Mettler Toledo DSC 1), while the specific heat capacity of the fibers c p, f , thermal conductivity of the matrix λ m and thermal conductivity of the fibers λ f were taken from the literature. To determine the sensitivity of the simulation to these thermodynamic material parameters, a parametric study was carried out. For this purpose, a 10-layer UD tape layup consisting of polycarbonate with 44% carbon fiber by volume was heated from 60 • C to 250 • C within 60 s. For the specific heat capacities c p,m , c p, f and the thermal conductivities λ m , λ f we additionally used values that were 20% higher and lower than those taken from the literature. All values used are listed in Table 2.

Mesh Study
A mesh study was additionally carried out to assess the influence of cell size distribution on heat transfer. Again, a 10-layer UD tape layup consisting of polycarbonate with 44% carbon fiber by volume was heated from 60 • C to 250 • C within 60 s. Both a coarse and a fine mesh were investigated. The coarse mesh had a cell size of 0.35 mm in the thickness direction; which means that one cell consisted of two tape layers. For the fine mesh the cell size was halved, and thus each cell in the thickness direction represented one layer of tape, as illustrated in Figure 5. The coarse and fine meshes consisted of a total of 1.5 million and 2.2 million cells, respectively.

Validation Study
The consolidation process modeled in this work is carried out by a consolidation unit which includes a heating and a cooling press. Within the consolidation unit, the tape layup is transferred from one press to the other between two steel plates, which are moved within the machine fully automatically (see Figure 6). During the heating process, less pressure is usually applied to the tape layup compared to the cooling process, to avoid extensive squeeze flow. Typical temperature and pressure profiles for this process can be seen in Figure 7. To analyze the validity of the simulation approach, an experimental parameter study was carried out with various temperature and pressure profiles using various materials, and the results were compared to those of the simulation. The process parameters considered are listed in Table 3. Table 3. Validation study process parameters: heating temperature (T H ), cooling temperature (T C ), heating pressure (p H ), cooling pressure (p C ) and cycle times for heating (t H ) and cooling (t C ) of polycarbonate with carbon fibers (PC-CF) and polypropylene with glass fibers (PP-GF). For case 1, two previously consolidated sheets were consolidated with each other, while in cases 2-4, 18 individual tapes were consolidated into one sheet. The material used for cases 1-5 was a polycarbonate with 44% carbon fiber by volume (PC-CF). In case 5, a 10-layer tape layup of polypropylene with 45.3% glass fiber by volume (PP-GF) was used and only the heating process was investigated.

Material Layup T H T C p H p C t H t C
In order to record the temperature within the composite part, a thermocouple (TC Type K) was placed between the sheets in case 1. In cases 2 to 5, three individual thermocouples were placed between the individual tape layers. In cases 2-4, they were placed between layers 1 and 2, 9 and 10, and 17 and 18. In case 5, they were placed between layers 1 and 2, 5 and 6, and 9 and 10. This made it possible to obtain a complete picture of the temperature behavior at multiple points in the composite and thus allowed comprehensive comparison between experiment and simulation.

Model Sensitivity to Changes in Thermal Material Properties
In the sensitivity study, each thermal parameter (i.e., specific heat capacities c p,m , c p, f and thermal conductivities λ m , λ f (see Table 2)) was varied separately by ±20% while keeping all other thermal parameters at their original values from the literature. The results of this study are shown in Figure 8. Due to the large temperature difference between the heating plates and the composite material, the temperature increased sharply at the beginning and reached a plateau after about 60 s. For all thermodynamic parameters, the numerically calculated temperature curves hardly changed between the different settings. This suggests that, within a specific range, the thermal parameters of the individual components-that is, specific heat capacity c p,m , c p, f and thermal conductivity λ m , λ fhave little influence on the heat transfer.

Mesh Study
For the mesh study, a 10-layer composite material that consisted of polycarbonate with 44% carbon fiber by volume was heated from 60°C to 250°C. The meshes were defined as described in Section 2.2.2. The coarse mesh consisted of 1.5 million cells and the fine mesh of 2.2 million cells. For comparison, the temperature in the center of the composite part over time was investigated, as illustrated in Figure 9. The results of this study indicate that the mesh resolution has minimal influence on the temperature behavior. The largest difference between the temperature curves amounted to 6.6 ºC after 30 s of heating, which corresponds to a difference in temperature of about 3%. The influence of mesh resolution on other phenomena that occur during consolidation, which were ignored in this work, (e.g., degree of bonding and squeeze flow) remains to be investigated. To save computational time, the subsequent parameter study was conducted using the coarse mesh.

Parameter Study
As described in Section 2.2.3 a parameter study was carried out to test our model's prediction accuracy against data from experiments using the consolidation unit shown in Figure 6. Each experiment was performed three times to record and minimize inconsistencies between runs, and the mean values were used to compare experiment with simulation.

Case 1
In case 1 the consolidation of two single polycarbonate sheets with 44% carbon fiber by volume was observed. Figure 10 shows the numerically and experimentally obtained temperature profiles, which are in good accordance. The dymanic temperature behavior during heating and cooling was reproduced very well by the simulation, but the target temperature of the heating or cooling processes was reached earlier than in the experiment. For process control (and, in a further step process optimization), it is important to know when a particular temperature is reached in the core of the composite to ensure that the material is completely molten during heating and completely solidified during cooling. Here, the glass transition temperature of polycarbonate (=145°C) was used. This temperature was reached after 10 s of heating in the experiment, while the simulation predicted a period of 11 s. During cooling, this temperature was reached 161.7 s into the experiment and after 165 s in the simulation. Based on the total process time of 250 s experiment and simulation differed by 0.4% in the heating process and by 1.32% in the cooling process.

Cases 2, 3 and 4
Cases 2, 3 and 4 are presented together because they used the same material, layup and cooling process, and differed only in the heating temperature set at the consolidation unit (see Table 3).
As mentioned in Section 2.2.3, the temperature was recorded and analyzed at three positions within the composite for these cases, where one position was located in the core and the other two at the outer layers. Figures 11-13 show the results of a comparison between experiment and simulation. It can be seen that the difference between experiment and simulation is comparable to that in case 1, which means that the dynamic behaviour during heating and cooling was predicted very well by the simulation.
For these cases we also determined when the glass transition temperature was reached in the core (between layer 9 and 10) of the composite material during heating and cooling, the results of which are summarized in Table 4.

Case 5
In order to assess the flexibility of the model in relation to the materials used, the temperature behavior of a 10-layer layup of polypropylene with 45.3% glass fiber by volume was investigated. As mentioned in Section 2.2.3 only heating from 50°C to 250°C was considered. The results of this experiment can be seen in Figure 14. In contrast to the cases shown above, the temperature curve in the outer layers (Figure 14a,b) shows irregularities around the melting temperature (173°C) which are not apparent in the simulation. A possible reason for the temperature fluctuations could be the melting behavior, which was not taken into account in the simulation and was not investigated further in this work.
Comparison between experiment and simulation in terms of the time when the melting temperature was reached shows a difference of 2.84 s. Relative to the total process time of 100 s this amounts to a deviation of the simulation from the experiment by 2.84%.
As already mentioned, the experimental results show a change in the slope of the curve in the outer layers (Figure 14a,b) around the melting temperature (T m = 173°C). Since the model does not consider phenomena related to crystallization, the simulation results cannot capture differences in temperature behavior between semi-crystalline and amorphous plastics.

Discussion and Conclusions
We have shown that our multi-region, multi-phase and multi-component-mixture approach is suitable for modeling the heat transfer during the consolidation process. It is robust within a range of ±20% of the thermal parameters taken from the literature and insensitive to variations in mesh resolution. Simulation results and experimental data were in very good agreement, especially in relation to the dynamic phase during heating and cooling. The simulation allows very precise estimation of when a particular temperature, such as the glass transition temperature or melting point, will be reached at the core of a composite. In relation to the total process time, deviation of the simulation from the experimental data were the following: 0.4% during heating and 1.32% during cooling for Case 1, 1.58% during heating and 0.55% during cooling for Case 2, 0.63% during heating and 0.98% during cooling for Case 3, and 0.65% during heating and 1.2% during cooling for Case 4.
Case 5 from our parameter study, using polypropylene with 45.3% glass fiber by volume, demonstrated that the model is also capable of predicting the temperature behavior of semi-crystalline matrix materials. Concerning the time when the melting temperature is reached within the core of the layup, the difference between simulation and experiment is 2.84% with respect to the total process time.
Our model, and in particular the accuracy of its predictions, enables process optimization in the form of reduction of cycle time and therefore improving the energy efficiency, and the building of digital twins. It will form the basis for implementing further models to investigate, for instance, degree of bonding, compression and squeeze-flow behavior and crystallinity. This gives the opportunity to consolidate thermoplastic composites in a more time-and energy-saving way.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.