**1. Introduction**

#### *1.1. Motivation and Materials*

Fiber reinforced composites have gained importance in the automotive industry due to their potential for lightweight applications. During compression molding of discontinuous fiber composites (DFC) the long fibers undergo low shear stresses and can flow in the melt without significant fiber breakage that is common during plastification and mold filling in the injection molding process [1,2]. The shorter the required flow length to fill the mold, the lower is the possibility of fiber attrition [3]. Since typical sheet molding compounds (SMC) charge coverages lead to comparably short flow paths, the reinforcing fibers maintain their length [3,4] and compression molded DFCs show more advanced mechanical properties with high mass-specific stiffness and strength than compared to injection molded long fiber reinforced polymers [3,5].

The low tooling costs (compared to steel processing) and fast cycle times make compression molding a cost-efficient method to manufacture large DFC parts in a one-shot high volume production process, enabling DFC parts to replace automotive metal components for mass-reduction purposes [6,7]. As this material class and manufacturing method allow a high freedom of design [5], part integration (e.g., fasteners or inserts) [7], and can be used to mold complex three-dimensional (3D) shaped structural and non-structural components with corrugations, ribs and domes, it is widely used among the automotive industry and a key aspect for the endeavor to reduce the vehicle weight as much as possible [5,8]. Among this type of DFC, carbon fiber sheet molding compounds (CF-SMC) have been extensively used for interior and exterior, structural and non-structural composite applications in the automotive and aerospace industry [7,9–12].

The focus of this work is on a DFC consisting of transversely chopped unidirectional (UD) carbon fiber prepreg tows, so called 'strands', 'chips' or 'platelets', which are randomly distributed into a mat. SMCs with those randomly oriented strands (ROS) show a high degree of heterogeneity (variability in intra- and inter-part structure on the meso- and macro-scale) ye<sup>t</sup> seek to reach quasi-isotropic mechanical properties [13–15]. This material class is also called prepreg platelet molding compound (PPMC) [4,16–27] or tow-based discontinuous composite (TBDC) [28–33]. The material system is characterized by a high fiber volume fraction with good impregnation, comparable to continuous fiber layups, and therefore higher mechanical properties are reachable than with traditional SMCs [18]. High performance CF-SMCs, such as the epoxy-based material used in this study (see Section 2.1), are further characterized by a high delamination resistance, near quasi-isotropic in-plane stiffness, high out-of-plane strength and stiffness, and low notch sensitivity [13]. Moreover, state of the art resin systems enable very short curing times, leading to an 84% shorter molding time and an overall process time reduction of 44% for a one-piece inner monocoque compared to the same part produced in a resin transfer molding (RTM) process [34].

However, due to the part design together with the high fiber volume fraction and fiber length in the CF-SMC material, complex anisotropic material flow conditions occur [21], which induce a characteristic microstructure in compression molded components [5,31]. This process-induced microstructure is mainly characterized by locally varying fiber orientations and fiber concentrations. Besides defects such as voids (air entrapments), swirls or resin rich pockets in between strands [35–40], the flow-induced strand alignments have the biggest impact on the mechanical performance of DFC parts [21,31,41]. Therefore, obtaining realistic 3D information of the local strand orientations in a CF-SMC part is key for a better understanding of its related mechanical behavior and a sophisticated part design. 3D representations of the morphology of entire CF-SMC parts help engineers to grasp the compression molding process related fiber alignment and the gained fiber orientation information can be mapped to a structural simulation mesh and used as a digital twin in an integrative (coupled) Finite Element Analysis (FEA) [42,43]. This microstructure information can be acquired by precise and reliable process simulations or laborious micro-computed tomography (μCT) measurements. CT scans that allow fiber orientation measurements of complete parts may be used to compare 'real' with numerically predicted fiber orientations coming from molding simulations in order to evaluate the quality of the simulation results. However, especially for parts made of carbon fiber reinforced plastics, this is so far a costly challenge due to the low contrast between the carbon fibers and the polymer matrix system in CT scans.

This paper shows the application of a direct fiber simulation (DFS) tool featuring a novel strand generation approach to simulate the compression molding of ROS-based CF-SMC parts. The strand generation feature is a pre-processing step to initialize multi-bundle UD strands in the initial numerical material charge. In order to validate the filling and fiber orientation prediction accuracy on real parts, automotive-related demonstrator parts with a complex ribbed structure are molded with a high-performance ROS-based CF-SMC material. As the accuracy of the numerical flow prediction is crucial for fiber orientation predictions, it is compared with the real part filling behavior observed in short shot experiments. Subsequently, simulated fiber orientations are compared with true process-induced strand orientations inside the molded parts, which are determined by CT scan analyses. The micro-, or more appropriately, mesostructures of three entire parts are non-destructively analyzed using low-resolution CT scans. The determined fiber orientations of all three CT scanned parts are averaged to eliminate local mesostructure di fferences between the parts and to gain a representative image of the general fiber orientations in the ribbed structure. In order to eliminate further influencing factors when comparing measured and predicted fiber orientations, both the process simulation and the CT scan analyses are conducted with exactly the same tetra mesh, enabling a one-to-one comparison. The size of the molded, simulated and CT scanned CF-SMC parts and the conducted detailed analyses surpass previously published studies dealing with this material class. The objective of this paper is to provide a good understanding of the process-induced mesostructure in complex ribbed DFC parts and to show the application of a commercially available DFS tool in a DFC part development process to avoid costly trial and error molding experiments and part design loops.

#### *1.2. Compression Molding Simulations*

Compression molded DFCs show distinct long fiber e ffects such as complex fiber orientations and fiber matrix separation (FMS) due to fiber entanglements and accumulations leading to fiber content irregularities and consequently to significant changes in mechanical properties [44,45]. As this process-induced microstructure has a strong e ffect on the part performance, molding simulation results such as fiber orientations and fiber content distributions are nowadays mapped to FEA meshes and used in integrative (coupled) simulations to account for the created anisotropy and to obtain realistic mechanical part behavior predictions. Therefore, accurate fiber configuration predictions are of major importance for the quality of subsequent structural simulations. Furthermore, high-quality virtual process chains help to avoid cumbersome trial-and-error molding experiments and design loops, which decreases the associated costs for the trials and mold changes and also accelerates the part development process. In conclusion, accurate numerical predictions are crucial to keep costs as low as possible and enable e fficient material use and safe lightweight part design.

Process simulation tools are widely used in the automotive industry to predict the flow front advancement during molding and to indicate adverse e ffects like short shots or knit lines [6]. These tools are further used for the calculation of changing fiber configurations during molding of fiber reinforced plastics. Process simulations of CF-SMC parts can be done with traditional numerical tools that use phenomenological models or with more elaborate methods like smoothed particle hydrodynamics (SPH) or direct fiber simulations (DFS).

#### 1.2.1. Statistical Fiber Orientation Models

Common commercially available molding simulation tools show good results predicting fiber orientation and fiber length changes during processing of short fiber reinforced plastics using standard empirical calculation models designed for short fibers [6,46–48]. Based on Je ffery's hydrodynamic model, the fibers are modeled as ellipsoidal shaped rigid bodies rotating in a viscous flow [47]. Leveraging the work on short fiber orientation behavior in concentrated suspensions by Folgar and Tucker [48], the change of fiber orientations during processing can be described by an evolution equation introduced by Advani and Tucker [49] stated in Equation (1):

$$\frac{\mathrm{D}a\_{ij}}{\mathrm{D}t} = -\frac{1}{2}(\omega\_{\mathrm{ik}}a\_{\mathrm{j}j} - a\_{\mathrm{ik}}\omega\_{\mathrm{kj}}) + \frac{1}{2}\lambda \Big(\dot{\boldsymbol{\gamma}}\_{\mathrm{ik}}a\_{\mathrm{k}j} + a\_{\mathrm{ik}}\dot{\boldsymbol{\gamma}}\_{\mathrm{kj}} - 2\dot{\boldsymbol{\gamma}}\_{\mathrm{k}}a\_{ij\mathrm{k}l}\Big) + 2\mathrm{C}\_{\mathrm{I}}\dot{\boldsymbol{\gamma}}\Big(\delta\_{ij} - 3a\_{ij}\Big). \tag{1}$$

This evolution equation for the second order orientation tensor *aij* is often used in contemporary molding software providing a compact expression for the fiber orientation state at reasonable calculation speed [2]. It includes the vorticity tensor <sup>ω</sup>*ij*, the rate of deformation . γ*ij*, the scalar magnitude of the rate of deformation tensor .γ, and the identity tensor δ*ij* (δ*ij* = I). λ denotes the particle shape parameter (for long fibers λ → 1), which is related to the fiber aspect ratio. The fiber interaction term 2*C*I .<sup>γ</sup><sup>δ</sup>*ij* − <sup>3</sup>*aij* results from the isotropic rotary diffusion (IRD) model by Folgar and Tucker and includes the empirical fiber interaction coefficient *C*I [48,50]. In order to determine the fourth order tensor *aijkl* in the evolution equation, a closure approximation is needed and various are proposed in the literature to achieve good simulation accuracy at reasonable calculation times [51–57].

In addition to this fundamental equation for isotropic rotary diffusion and fiber dynamics [48,49], some enhancements by successor models that are implemented in conventional software improved the ability to take fiber–fiber interactions into account by implementing various anisotropic rotary diffusion (ARD) models [50,58–61]. The observed orientation delay between the measured fiber orientation and the theoretical orientation evolution is considered through strain reduction models [62–64]. Excluded volume effects can be considered to describe the dependency of the fiber interaction and the volume fraction of long fiber filled materials [65]. Recently the flow and fiber orientation prediction was coupled using anisotropic viscosity models to analyze the behavior of concentrated fiber suspensions where the fiber orientation state affects the viscosity [66–71].

These statistical fiber orientation models strongly depend on empirical coefficients, which are cumbersome to determine experimentally [72,73]. Furthermore, the model simplifications and boundary condition assumptions made in these models are not valid for flexible fibers significantly longer than structural features of the part (scale separation) [21,45,73]. In addition, phenomenological models do not consider segments of one fiber being in different flow conditions (at different positions inside the material flow) at the same time, which is highly probable for long fibers [45]. Consequently, velocity and strain rate distributions along the fiber axes cannot be taken into account and the prediction of fiber bending is not feasible [2,74]. Since interactions of long bendable fibers lead to fiber entanglements and accumulations influencing the fiber movement, the capability to model fiber bending is a crucial element [15].

In conclusion, the used empirical models implemented in conventional molding simulation tools, such as Moldflow (Autodesk, Inc., San Rafael, CA, USA) or Moldex3D (CoreTech System Co., Ltd., Zhubei City, Taiwan), are designed for rigid particles with low aspect ratios and are therefore not suited for highly filled long fiber materials. As expected and seen in the literature, even with modified evolution equations, flow-induced long fiber specific effects like fiber bending, fiber–fiber interactions, and changes in fiber content distribution (especially FMS) cannot be accurately predicted by traditional short fiber models [2,6,45,46].

#### 1.2.2. Stochastic Particle-Based Simulation Method (Purdue University)

A particle-based flow simulation method with a stochastic approach was recently developed at the Purdue University (West Lafayette, IN, USA) to investigate ROS materials with flow-induced fiber alignment [4,17,21,25,75,76]. In contrast to conventional flow simulations typically executed in Eulerian frameworks [25], a smoothed particle hydrodynamics (SPH) method with a variable user-defined material subroutine (VUMAT) is applied in Abaqus/Explicit achieving a Lagrangian solution that allows large deformations required by mold filling [17,21,22,25,76]. The 3D flow molding simulation method is based on a stochastically generated random planar initial strand orientation distribution [17]. To predict the complex 3D strand deformations and the resulting orientations in the final part, it calculates the strand orientation evolution based on the assumption of affine motion (equivalent to Jeffery's equation) coupled with an anisotropic viscosity model [17,21,22,25,76].

The initially introduced variability of the charge orientation state is generated with a pseudo-random number generator (strand generation scheme) [25]. Coded in Python, the strand generator randomly defines in-plane strand centroid locations with random strand orientation angles as long as all SPH particles in the layer are assigned to one strand set [22,25]. Accordingly, the strands are modeled as a set of interconnected SPH particles ('groupings' [4,17]) forming a rectangular strand shape with one particle through the thickness [25]. This leads to a low resolution of the orientation distribution through the thickness of the flow simulation [17]. However, the particle size cannot be set to be equal to the strand thickness, since the computation time would be unreasonably high [17]. This is due to the fact that the particle size equally reduces in all three spatial directions so that a huge amount of particles would be needed to describe one strand.

During the flow simulation the initial charge is extended and the strand sets deform and progressively align along their flow direction tending to disaggregate [17]. Since an anisotropic (transversely isotropic) viscosity model is used, strands that are initially oriented in flow direction tend to translate, while those oriented in cross-flow direction widen [17]. Both e ffects are qualitatively consistent with physical observations [17,31]. Moreover, strand–strand interfaces [17] and interactions are only implicitly taken into account through the anisotropic viscosity model [21].

#### 1.2.3. Direct Fiber Simulations (DFS)

In order to overcome the deficiencies of the described phenomenological models, many single fiber based models were devolved. Single fiber simulations, also called direct fiber simulations (DFS), connect a number of particles or rigid rods (beams) to model long flexible fibers inside the polymer flow [2,45,74,77–84]. Because of the fiber segmentation DFSs can consider the strain rate distribution along the fiber axis and are therefore suited to predict long fiber behavior like bending [74]. These models use the lubrication theory combined with small flexible inextensible threads [85], stretchable, bendable, and twistable chains of bonded spheres [78] or rigid spheres connected by ball and socket joints [86–89]. Contact formulations describe the fiber–fiber interactions [90,91] and fiber–fluid interactions are modeled by hydrodynamic drag forces [92–94].

Two types of DFS techniques can be distinguished—the velocity-based method and the mechanistic model [45]. Kuhn et al. show notable amelioration in prediction precision when using DFS tools compared to commercial phenomenological simulation tools predicting long fiber behavior in a small rib structure [45]. However, modeling a large amount of connected particles leads to long computation times if applied to practical/industrial parts, whereas the duration can be decreased by using connected rigid rods to model long fibers with a velocity-based method [74]. The same observations are made by Kuhn et al. who describe decreased preparation e fforts and calculation durations when using a commercial velocity-based non-interaction DFS instead of a non-commercial, research-focused mechanistic DFS [45]. However, due to the enormous amount of calculated fibers necessary, and the therefore very numerically expensive task of single fiber simulations, this method is often restricted to predict the resulting fiber configuration in small volumes [45].

#### Mechanistic Model (University of Wisconsin-Madison)

Kuhn et al. use a bead chain model for a DFS of a single rib of a molded component using a reduced amount of fibers [6,45,46]. The presented study uses a particle level simulation approach based on a mechanistic model to simulate fiber bending, fiber interactions, and especially FMS during compression molding of long fiber reinforced plastics [6,45,46]. The mechanistic model simulation approach begins with a traditional mold filling calculation with a standard tool like Moldex3D [46]. The calculated flow field is then extracted and used as basis for the mechanistic model simulation. After a stack of fibers is randomly inserted into the cavity volume, the movement of all single fibers is determined by interaction with the flow field. In the mechanistic model each fiber is individually modeled as a chain of rigid rods connected by ball-and-socket joints [6,46,95,96]. The fiber discretization numerically enables fiber bending at the joints (nodes) between the rods, which is important for the prediction of long fiber behavior in component areas that are smaller than the fiber length [46,97]. Furthermore, the inside of the rods is modeled as a chain of balls to include fiber interactions and hydrodynamic e ffects [6]. Accordingly, the fibers are simulated as a collection of equidistant nodes where the connecting rods experience elastic deformation, hydrodynamic drag forces and excluded volume e ffects in the polymer flow field [46,97]. The contact and hydrodynamic forces are used in force and momentum balance equations to calculate the rotations and movements of each modeled

fiber. Additionally, complex interactions with the enclosing melt flow, other fibers and mold walls are considered [45,46,96,98–100]. The excluded volume forces act as repulsive forces to inhibit fibers to interpenetrate or overlap one another or with the mold wall and are used to model fiber–fiber as well as fiber–mold interactions [46,97,98].

Due to the high fiber sti ffness and the more significant impact of fiber bending, extensional fiber deformation is neglected [46]. In order to decrease the tremendous calculation times required, fiber–fiber and fiber–mold interactions are only computed once every 50 integrations [46]. However, due to the still limited amount of calculable fibers or rather the overall low simulation speed, the mechanistic model is currently limited to small part volumes [45]. In sum, the mechanistic model calculates the motion, bending, and rotation of single fibers considering all their complex interactions with the fluid, the mold walls, and other fibers. Although the mechanistic model describes only a one-way coupling from the fluid flow to the fibers, the direct fiber simulation approach results in more accurate fiber orientation and fiber content distribution predictions when compared to common process simulation software [6,45].

#### Direct Bundle Simulation (DBS)

For a DFS of a small compression molded cross-rib-shaped SMC part in LS-DYNA (Livermore Software Technology Corporation (LSTC), Livermore, CA, USA) Hayashi et al. use beam elements constrained in highly deformable solid elements representing the matrix [101]. In order to enable large deformations of the matrix elements, the 3D adaptive Element-Free-Galerkin (EFG) method is applied. Recently, this tool was enhanced to handle ROS composites. The fibers within one strand are modeled by multiple connected elastic beam elements in a row [102]. Since the real amount of fibers in ROS composites and even within a strand is very large, the calculated number of fibers per strand was drastically reduced, while the beam thickness was increased to ensure maintaining the nominal fiber volume fraction of the material [102]. This simplification method reduces the numerical calculation time by representing thousands of fibers as fiber bundles. Considering microstructure studies reported in the literature [4,73,103,104], it can be seen that most of the fiber bundles or strands widen and flatten, but do not disperse. Only the highly sheared strands are more likely to fan out. This observation justifies the simplifying assumption to represent thousands of carbon fibers within a strand as one bundle in order to shorten the calculation time and to enable full component scale simulations [73].

Recently, Meyer et al. [73] presented a DBS method that reproduces Je ffery's equation for single fiber bundles in shear flow, which are described as a chain of truss elements. Since the bundles are represented as one-dimensional instances, they can move independently from the matrix material, flow and interact through contact forces and with hydrodynamic drag forces of the surrounding flow. Since bundle–wall and bundle–bundle interactions are considered, there is no need for empirical interaction parameters, unlike in the commonly used statistical descriptors of fiber orientation changes (Folgar–Tucker-based models). Since this DBS approach additionally considers two-way coupling it allows for a more accurate calculation of fiber volume fraction distributions, knit lines and FMS at a component level with reasonable computation times [73].

#### 1.2.3.3. 3D TIMON CompositePRESS

#### Fiber Orientation Simulation in 3D TIMON CompositePRESS

A commercially available compression molding simulation tool with a similar DFS approach is 3D TIMON CompositePRESS (Toray Engineering D Solutions Co., Ltd., Otsu, Shiga, Japan) [ ¯ 105]. Here, each fiber within the initial fiber cluster is modeled using a flexible, non-elastic chain of rods (rigid bodies, constant rod length) connected by hinge nodes [2], as shown in Figure 1. The velocity-based DFS tool begins with the filling simulation determining the flow field, assuming homogeneous isotropic resin properties (isotropic viscosity) and using a one-way coupling from the fluid flow to the fibers. Subsequently, it calculates the motion of each single fiber based on the flow velocities [2,45]. The position and curvature of each single fiber are changed according to the strain rate distribution along the fiber length [74]. Thereby, the nodes function as universal joints allowing the rods to rotate in three independent axes and thus enabling the calculation of movement, bending, and rotation of each single long fiber [2,106]. In this velocity-based direct fiber simulation approach the fiber node position *ui* is a function of the enclosing fluid velocity *v*fluid at every time step *n* [2]. The melt velocity *v*fluid is interpolated from the node positions in the velocity field [2]. Since the motion *ui* of node *i* directly follows the surrounding melt flow velocity, the temporary new node position for the next time step *n* + 1 can be calculated by Equation (2) [45,107]:

$$
\mu\_i^{n+1} = \mu\_i^n + \sigma\_{\text{fluid},i} \prescript{n}{}{\text{dt.}}\tag{2}
$$

**Figure 1.** Schematic illustration of a flexible long fiber modeled with 6 hinge nodes and 5 rods in 3D TIMON shown in its initial straight form and bent after moved by fluid forces and applied rod length adjustments (node relocations).

This leads to an unrealistic and impermissible stretching of the fibers, which is numerically corrected by rod length adjustments (cf. Figure 1) or node relocations, respectively. At each time step the change of the rod length *L* is reviewed to a maximum elongation threshold ε (Equation (3)):

$$\max \left| 1 - \frac{L\_{ij}^{\;n+1}}{L\_{ij}^{\;n}} \right| < \varepsilon. \tag{3}$$

If the rod length change exceeds this threshold, node *i* is relocated considering the neighboring nodes and rod lengths by Equation (4):

$$
\mathfrak{u}\_{i}^{n+1} = \mathfrak{u}\_{i}^{n} + \mathfrak{d}\mathfrak{u}\_{i\star} \tag{4}
$$

with

$$\mathbf{d}u\_{i} = \left(\mathbf{u}\_{i}^{n+1} - \mathbf{u}\_{i}^{n}\right) + \frac{1}{2} \Big(1 - \frac{L\_{i\bar{j}}^{n+1}}{L\_{i\bar{j}}^{n}}\Big) \left(\mathbf{u}\_{\bar{j}}^{n+1} - \mathbf{u}\_{i}^{n+1}\right) + \frac{1}{2} \Big(1 - \frac{L\_{i\bar{k}}^{n+1}}{L\_{i\bar{k}}^{n}}\Big) \left(\mathbf{u}\_{k}^{n+1} - \mathbf{u}\_{i}^{n+1}\right),\tag{5}$$

where *Lij* is the rod length between node *i* and node *j* and *Lik* is the rod length between node *i* and node *k*. Equation (5) is iteratively solved until the rod length changes are smaller than ε and the node positions are updated by time integration based on an explicit Euler scheme [2].

**Figure 2.** Schematic illustration of a tetra element in the morphing method with one node at the top, 3 nodes at the bottom, and 20 calculation points through the thickness for Light 3D analyses in 3D TIMON CompositePRESS.

Outperforming conventional fiber orientation models, the velocity-based DFS method applied in 3D TIMON CompositePRESS enables the prediction of local fiber orientations and fiber content distributions in compression molded SMC parts where long fibers (typically 25–50 mm) are used [2]. Moreover, fiber breakage can be simulated through detachment of rod–rod connections [2]. However, due to the fact that fiber attrition is of minor importance during SMC compressing molding, this effect is typically neglected. Fiber–fiber interactions are not considered in the current R6.0 version of 3D TIMON 10 to keep the calculation time reasonable [2,73,105]. Furthermore, the software does not account for anisotropic viscosity and two-way coupling [73]. Compared to phenomenological models though, the long fiber behavior is superiorly predicted, despite being based on non-colliding, velocity-following nodes [45]. However, it is reported that fiber content distribution simulation results have some insufficiencies when compared with experiments. Kuhn et al. conclude that these discrepancies are induced by extensive fiber interactions causing FMS in reality, which cannot be accurately predicted by either phenomenological (statistical) models based on Morris and Boulay or the velocity-based DFS by 3D TIMON [45].

Due to the current status of the presented long fiber orientation models, Kuhn et al. recommend to use velocity-based DFSs like 3D TIMON for the fiber orientation simulation of large SMC components [45]. However, for smaller volumes, where FMS is pronounced and of interest, they sugges<sup>t</sup> to use the mechanistic model to be able to predict fiber agglomeration more accurately [45]. Kuhn et al. further propose to develop and implement a phenomenological or simplified model for FMS in common computational tools like Moldflow, Moldex3D or 3D TIMON [45].

#### Flow Simulation Simplifications in 3D TIMON CompositePRESS

In order to simulate the compression molding process, the flow analysis in 3D TIMON is based on the standard equations of traditional fluid dynamics: the continuity equation, the momentum equation, and the energy equation. These governing equations for the filling phase are given in the following. Theoretical background and a more detailed deduction of the general and special forms of the equations are given in Appendix A.

**Continuity equation:**

$$\frac{\partial \rho}{\partial t} + \frac{\partial}{\partial x}(\rho u\_\mathbf{x}) + \frac{\partial}{\partial y}(\rho u\_\mathbf{y}) + \frac{\partial}{\partial z}(\rho u\_\mathbf{z}) = 0,\tag{6}$$

where ρ is the density, *t* is time and *ui* stands for the fluid velocity vector in direction x, y, and z, respectively. Assuming incompressibility and a steady state, neglecting inertia and gravity, and considering the fluid motion in z-direction, the governing equations for non-Newtonian polymer melt flow can be simplified to a volume continuity equation:

$$
\frac{
\partial \mathfrak{u}\_{\mathbf{x}}
}{
\partial \mathbf{x}
} + \frac{
\partial \mathfrak{u}\_{\mathbf{y}}
}{
\partial y
} + \frac{
\partial \mathfrak{u}\_{\mathbf{z}}
}{
\partial z
} = 0.
\tag{7}
$$

**Momentum equation** in x-direction:

$$\rho \left( \frac{\partial \mathbf{u\_{x}}}{\partial t} + \mathbf{u\_{x}} \frac{\partial \mathbf{u\_{x}}}{\partial \mathbf{x}} + \mathbf{u\_{y}} \frac{\partial \mathbf{u\_{x}}}{\partial y} + \mathbf{u\_{x}} \frac{\partial \mathbf{u\_{x}}}{\partial z} \right) = -\frac{\partial p}{\partial \mathbf{x}} + \mu \left( \frac{\partial^{2} \mathbf{u\_{x}}}{\partial \mathbf{x}^{2}} + \frac{\partial^{2} \mathbf{u\_{x}}}{\partial y^{2}} + \frac{\partial^{2} \mathbf{u\_{x}}}{\partial z^{2}} \right) + \rho \mathbf{g\_{x}}.\tag{8}$$

where *p*: pressure, μ: viscosity of Newtonian fluids, *gi*: body force acting on the continuum, for example, gravity. This equation can also be expressed in terms of deviatoric stress τ and is then commonly called the Cauchy momentum equation:

$$\rho \left( \frac{\partial u\_{\mathbf{x}}}{\partial t} + u\_{\mathbf{x}} \frac{\partial u\_{\mathbf{x}}}{\partial \mathbf{x}} + u\_{\mathbf{y}} \frac{\partial u\_{\mathbf{x}}}{\partial \mathbf{y}} + u\_{\mathbf{z}} \frac{\partial u\_{\mathbf{x}}}{\partial \mathbf{z}} \right) = -\frac{\partial p}{\partial \mathbf{x}} + \left( \frac{\partial \tau\_{\mathbf{x}\mathbf{x}}}{\partial \mathbf{x}} + \frac{\partial \tau\_{\mathbf{y}\mathbf{x}}}{\partial \mathbf{y}} + \frac{\partial \tau\_{\mathbf{z}\mathbf{x}}}{\partial \mathbf{z}} \right) + \rho g\_{\mathbf{x}}.\tag{9}$$

**Energy equation** for a fluid with constant properties is given by:

$$
\rho c\_p \frac{\text{DT}}{\text{Dt}} = k \left( \frac{\partial^2 T}{\partial \mathbf{x}^2} + \frac{\partial^2 T}{\partial y^2} + \frac{\partial^2 T}{\partial z^2} \right) + \dot{Q}\_{\text{viscous heating}} + \dot{Q}\_{\text{\textquotedblleft}} \tag{10}
$$

.

where ρ: density, *cp*: specific heat, *T*: temperature, *t*: time, *k*: thermal conductivity, *Qviscous heating*: viscous dissipation and . *Q*: arbitrary heat source (i.e., exothermic reaction).

For polymer flows, the fundamental governing equations given above are non-linear, non-unique, complex, do not have a general solution, and are thus difficult to solve [108]. To gain analytical solutions and to reduce the required CPU time the balance equations must be simplified. 3D TIMON CompositePRESS avoids solving the Navier–Stokes equations to realize a 3D flow analysis of non-Newtonian polymer melts, and instead uses the traditional Hele-Shaw approximation for a 2.5D flow simulation (called Light 3D) conducted on a one layer tetra mesh when using the morphing method [109]. The Light 3D solver virtually divides each mesh element (tetra or hexa) into 20 calculation points across the thickness (between the nodes on the surfaces) in order to predict the flow front advancement realistically [109] (see Figure 2). This procedure allows the calculation of temperatures, velocities, shear rates, and viscosity changes through the thickness by using a single layer mesh [109]. The Light 3D method reduces the number of necessary elements for accurate flow predictions for large and thin parts drastically [109].

The Hele-Shaw simplification, which is commonly used in injection molding simulations, was also applied by Folgar and Tucker to solve the compression molding process [110,111]. By assuming a flow through a narrow cavity gap *h* between lower and upper platen an order-of-magnitude analysis shows that the flow rate in thickness direction is small compared to the x- and y-directions. Furthermore, it shows that the in-plane velocity gradient can be neglected [109]. Therefore, in a thin layer flow field the shearing stresses ∂τzx ∂*z* across the narrow cavity gap *h* are dominant [112]. These assumptions allow simplifying the volume continuity and the momentum equation to:

$$\frac{\partial u\_{\mathbf{x}}}{\partial \mathbf{x}} + \frac{\partial u\_{\mathbf{y}}}{\partial y} = 0,\tag{11}$$

$$
\frac{\partial p}{\partial \mathbf{x}} = \frac{\partial \tau\_{\mathbf{x}\mathbf{x}}}{\partial \mathbf{z}} \,, \tag{12}
$$

$$\frac{\partial p}{\partial y} = \frac{\partial \tau\_{\text{xy}}}{\partial z}. \tag{13}$$

The 2.5D Hele-Shaw simplification method applied in 3D TIMON CompositePRESS aims at solving the momentum equation with the continuity equation. Therefore, the shear stresses must be derived from the three unknown variables in a 2D flow: pressure *p* and the flow velocities *u*x, and *<sup>u</sup>*y. As the flow field through the part thickness is neglected in the Hele-Shaw model, the velocity across the thickness *h* is integrated from the lower platen to the upper platen to compute the gap-wise average velocities *u*x and *<sup>u</sup>*y. When the continuity equation is integrated over the thickness it takes the form:

$$\frac{\partial}{\partial x} \int\_0^h u\_\mathbf{x} \mathrm{d}z + \frac{\partial}{\partial y} \int\_0^h u\_\mathbf{y} \mathrm{d}z = 0. \tag{14}$$

This integrated continuity equation further reduces to:

$$
\frac{
\partial
}{
\partial \mathbf{x}
}(h\overline{\boldsymbol{u}}\_{\mathbf{x}}) + \frac{
\partial
}{
\partial y
}(h\overline{\boldsymbol{u}}\_{\mathbf{y}}) = 0.
\tag{15}
$$

In order to solve the continuity equation (Equation (11)) 3D TIMON CompositePRESS uses an assumption of potential viscous flow to express the flow velocity components *u*x and *<sup>u</sup>*y, which can be written as:

$$\overline{u}\_{\chi} = -\frac{S}{h} \Big( \frac{\partial p}{\partial \chi} \Big) \tag{16}$$

$$
\overline{u}\_{\rm y} = -\frac{\mathcal{S}}{h} \left( \frac{\partial p}{\partial y} \right) \tag{17}
$$

where *h* is the gap height between the mold halves and *S* is the viscosity-dependent flow conductance for 2.5D flow analysis of thin parts, defined by:

$$S = \int\_0^\hbar \frac{\left(z - \lambda\right)^2}{\eta} d\mathbf{z},\tag{18}$$

where λ is the local value of *z* at which the shear stresses are zero. Since most problems are symmetric, λ = *h*/2.

Equations (16) and (17) show that the flow rate in each direction is proportional to its pressure gradient. Substituting the unknown flow velocities *u*x and *<sup>u</sup>*y in Equation (11) with the gap-wise average velocities from Equations (16) and (17) reduces the number of unknowns at each node in the calculation from three to one (pressure only) in Equation (19), which is the classic Hele-Shaw model:

$$
\frac{
\partial
}{
\partial \mathbf{x}
} \left( \mathbf{S} \frac{
\partial p
}{
\partial \mathbf{x}
} \right) + \frac{
\partial
}{
\partial y
} \left( \mathbf{S} \frac{
\partial p
}{
\partial y
} \right) = 0.
\tag{19}
$$

However, for the compression molding process the *z*-velocity component cannot be fully neglected and therefore the mold closing speed . *h* is included as an extra term in the continuity equation:

$$
\frac{
\partial
}{
\partial \mathbf{x}
} (h\overline{\boldsymbol{u}}\_{\mathbf{x}}) + \frac{
\partial
}{
\partial y
} (h\overline{\boldsymbol{u}}\_{\mathbf{y}}) + \dot{\boldsymbol{h}} = \mathbf{0},
\tag{20}
$$

which, when implementing Equations (16) and (17), leads to the Hele-Shaw model for compression molding:

$$
\frac{
\partial
}{
\partial \mathbf{x}
} \left( \mathbf{S} \frac{
\partial p
}{
\partial \mathbf{x}
} \right) + \frac{
\partial
}{
\partial y
} \left( \mathbf{S} \frac{
\partial p
}{
\partial y
} \right) - \dot{h} = 0.
\tag{21}
$$

Since the flow conductance *S* depends on the temperature-dependent viscosity η, the temperature field must be calculated. In order to calculate the temperature distribution in the thermoset molding compound, 3D TIMON CompositePRESS assumes the heat conduction in thickness direction is

dominant [105] and therefore simplifies the energy equation (Equation (10) and Equation (A10) in Appendix A) to the energy equation for thermoset materials in 2.5D flow, written as

$$\rho c\_p \left( \frac{\partial T}{\partial t} + u\_\mathbf{x} \frac{\partial T}{\partial \mathbf{x}} + u\_\mathbf{y} \frac{\partial T}{\partial \mathbf{y}} + u\_\mathbf{z} \frac{\partial T}{\partial \mathbf{z}} \right) = k \frac{\partial^2 T}{\partial z^2} + \tau\_\mathbf{yz} \left( \frac{\partial u\_\mathbf{y}}{\partial z} \right) + \tau\_\mathbf{xc} \left( \frac{\partial u\_\mathbf{x}}{\partial z} \right) + \dot{Q}\_\prime \tag{22}$$

where ρ: density, *cp*: specific heat, *T*: temperature, *ui*: fluid velocity vector in direction x, y, and z, respectively, *k*: thermal conductivity, <sup>τ</sup>*ij*: deviatoric stress, and . *Q*: heat generation due to the exothermic reaction in the thermoset curing process.

Using the simplification <sup>τ</sup>yz = η<sup>∂</sup>*<sup>u</sup>*y ∂*z* and τzx = η ∂*<sup>u</sup>*x ∂*z* (assumption of simple shear flow) yields to the final energy equation used in Light 3D heat transfer analyses in 3D TIMON CompositePRESS:

$$
\rho c\_p \left( \frac{\partial T}{\partial t} + \mu\_\mathbf{x} \frac{\partial T}{\partial \mathbf{x}} + \mu\_\mathbf{y} \frac{\partial T}{\partial y} + \mu\_\mathbf{z} \frac{\partial T}{\partial z} \right) = k \frac{\partial^2 T}{\partial z^2} + \eta \left( \frac{\partial u\_\mathbf{y}}{\partial z} \right)^2 + \eta \left( \frac{\partial u\_\mathbf{x}}{\partial z} \right)^2 + Q \nu \frac{\mathbf{d}\alpha}{\mathbf{d}t} \tag{23}
$$

where *Q*0 denotes the total amount of generated heat in the exothermic curing reaction and dαd*t* is the curing reaction rate.

In conclusion, the presented simplification method of determining the flow conductance facilitates the pressure field calculation. Subsequently, the temperature field is calculated and the flow front advancement during the 2.5D filling simulation can be predicted. However, in more complex parts with thicker ribs the narrow gap assumption of this simplified approach has its limits. Furthermore, the model does not account for the plug flow with a slip boundary condition at the mold surface, so that a proper choice of the flow conductance is important. For a more detailed explanation of the balance equations and their transformations, the simplifications made, and general background knowledge, the reader is referred to Appendix A and [1,108,110,111,113–116]. For more information about the flow analysis in 3D TIMON CompositePRESS the reader is referred to [105,109].

#### Viscosity Models in 3D TIMON CompositePRESS

During the compression molding process the thermoset resin is continuously heated up by the hot mold walls and the self-generating curing reaction heat leading to a viscosity drop. Simultaneously, due to the temperature increase, the curing reaction accelerates causing a conversion rate and viscosity increase. In order to describe these dependencies between viscosity, shear rate, temperature, and the curing reaction rate of the polymer melt, 3D TIMON combines three models: the Arrhenius-like temperature-dependent Andrade model (Equation (24)) [117], the curing reaction rate-dependent Castro–Macosko model (Equation (25)) [118,119], and the shear rate-dependent Cross model (Equation (26)) [120].

**Andrade model:**

$$
\eta\_0(T) = A \exp\left(\frac{B}{T}\right). \tag{24}
$$

with η0: zero shear viscosity (initial resin viscosity before curing), *A*: empirical model constant (fitted), *B*: material specific constant in Kelvin, *T*: melt temperature in Kelvin.

The **Castro–Macosko model** in Equation (25) is a widely used model in molding simulation tools, which describes the viscosity of thermoset materials as a function of temperature *T* and degree of cure α that can be expressed as follows:

$$
\eta\_1(T, \alpha) = \eta\_0(T) \left( \frac{\alpha\_{\rm gel}}{\alpha\_{\rm gel} - \alpha} \right)^{(D+E\alpha)}.\tag{25}
$$

Here, η1 is the viscosity of the resin at a given degree of cure (conversion) α, η0 is the initial resin viscosity before curing, <sup>α</sup>ge<sup>l</sup> is the degree of cure at the gel point, and *D* and *E* are empirical model constants that fit the experimental data.

#### **Modified Cross model:**

$$\eta(T,\alpha,\dot{\gamma}) = \frac{\eta\_1(T,\alpha)}{1 + \left(\frac{\eta\_1(T,\alpha)\dot{\gamma}}{\pi^\*}\right)^{(1-n)^\*}},\tag{26}$$

with .γ: shear rate, τ<sup>∗</sup>: critical shear stress (stress level at which the viscosity is during the transition between zero shear region (Newtonian plateau) and the shear thinning (power law) region of the viscosity curve), *n*: power law index. This model considers the effect of shear rate and temperature on the viscosity and describes the shear thinning behavior by the power law index *n*.

The combination of all three viscosity models extends the Castro–Macosko model with a power-law type shear rate dependence and leads to a modified **Cross–Castro–Macosko model**:

$$\eta = \frac{A \exp\left(\frac{B}{T}\right) \left(\frac{a\_{\text{gel}}}{a\_{\text{gel}} - a}\right)^{(D+E\alpha)}}{1 + \left\{ \frac{\left| A \exp\left(\frac{B}{T}\right) \left(\frac{a\_{\text{gel}}}{a\_{\text{gel}} - a}\right)^{(D+E\alpha)}\right|}{\tau^\*} \right\}^{(1-n)'} } \tag{27}$$

which is the finally used viscosity model in 3D TIMON CompositePRESS [105]. This model considers the dependence of the thermoset resin viscosity on temperature, degree of cure (conversion), and shear rate.

#### Curing Model in 3D TIMON CompositePRESS

The widely used semi-empirical Kamal model [121], which is a further development of a model proposed earlier by Kamal and Sourour [122], is applied in 3D TIMON CompositePRESS in order to describe the temperature-driven, *n*th-order autocatalytic reaction of thermoset resins. Equation (28) expresses the curing reaction velocity as a product of a temperature term following the Arrhenius temperature dependence and a function of the degree of cure α. Typically, the resin curing behavior is measured by differential scanning calorimetry (DSC) analyses and the six unknown constants in Equation (28) can be fitted to the DSC data.

**Kamal model:**

$$\frac{d\alpha}{dt} = \left[A\_1 \exp\left(-\frac{E\_1}{RT}\right) + A\_2 \exp\left(-\frac{E\_2}{RT}\right) \alpha^m\right] \left(1 - \alpha\right)^n,\tag{28}$$

where dαd*t* stands for the curing reaction rate, *A*1, and *A*2 for pre-exponential factors, *E*1 and *E*2 for activation energies, *R* for the universal gas constant and *m* and *n* are reaction orders.

During the flow simulation this curing reaction model calculates the reaction rate, which is then used in the Cross–Castro–Macosko model to obtain the current viscosity of the flowing and simultaneously curing material. This coupling allows for a viscosity model that considers the effect of the curing kinetics on the conversion rate and the corresponding viscosity.

#### 1.2.4. Computed Tomography of ROS Materials

In order to collect experimental 3D information about fiber orientation distributions of heterogeneous fiber reinforced materials, X-ray computed tomography (CT) is widely used in industry as a non-destructive measuring method due to the fact that it is easy to prepare samples and only requires a difference in the density-dependent linear X-ray attenuation coefficients of the matrix and the reinforcement [123]. The morphology of inhomogeneous materials like CF-SMCs can be investigated three-dimensionally by micro-CT (μ-CT), which is a high resolution X-ray CT method, allowing an in-depth material characterization [123].

To gain microstructural information with this experimental technique the sample is exposed to X-rays while incrementally rotating on a platform between the X-ray source and the detector (cf. Figure 11a). For each angle increment the sample is penetrated with radiation for a specific exposure time. Denser parts of the sample (fibers) absorb more radiation and therefore appear brighter in CT images than less dense material (resin). Since the X-ray radiation is proportionally attenuated as a function of the material's density, the detector can record a shadow projection of the sample revealing its inner structure [124]. The sample is fractionally rotated and irradiated until radiographic projections of a full 180◦ sample rotation are recorded. This acquired set of angular projections is su fficient to be reconstructed into a 3D model consisting of a large number of parallel micro-slice images by applying specific mathematical algorithms [124]. In order to improve the scan quality, a 360◦ rotation is usually carried out and several recordings per increment can be taken and averaged. The achievable scan quality (resolution and image sharpness) further depends on the spot/focal size of the X-ray tube, geometrical magnification, detector quality, vibrations during the image recording and the chosen combination of scan parameters [125].

In the reconstructed volumetric representation of the internal sample structure the scanned material is defined and visualized by its grayscale values (high density = high grayscale value). Therefore, composites show two maxima in the grayscale value histogram—one for the fibers and one for the matrix material. In order to analyze the fiber constitution, especially the fiber orientation, a threshold between both maxima is set by the user to define the surface of the fibers or strands, respectively. However, as this is a user-dependent manual act, analyzed fiber volume content distributions should only be used normalized.

In order to ensure a proper component design regarding mechanical requirements and for quality assurance of manufactured composite parts, the industry's endeavor is to determine the material's microstructure for large areas or ideally for entire parts [43]. However, due to the lack of contrast between polymer resins and carbon fibers (similar linear X-ray attenuation coe fficients) in X-rayed composites, μ-CT scans of entire CF-SMC components are so far limited in size when receiving useful data for fiber orientation analyses is required. Normally, attaining fiber orientations by CT scan data analysis of fiber reinforced polymer parts requires a finer scan resolution than the fiber's diameter to distinguish between individual fibers, which by implication limits the scan volume size [126]. Obtaining CT data for a larger 3D part is therefore always a trade-o ff between scan volume size, or part size, respectively, achievable resolution (voxel size) and the required scanning time. In conclusion, with increasing sample size, the achievable scan resolution decreases, which makes it very time-consuming and costly or even impossible to resolve fiber-scale details in CT scans of entire composite parts [126].

Useful μ-CT scans for CF-SMC parts, however, just need a su fficient resolution so that it is clearly distinguishable between strands and resin by grayscale value di fferences related to local relative fiber volume fraction variations enabling a fiber orientation analysis with a common CT scan analysis software (e.g., the commercially available software package VGSTUDIO MAX 3.3., Volume Graphics GmbH, Heidelberg, Germany). The analysis algorithm within VGSTUDIO MAX 3.3 (VG hereafter) is indeed intended to be used for the orientation analysis of discretely visible fibers [42,127–130]. Yet with correctly set parameters and with local relative density gradients between resin and fibers the image analysis principles are suited to be used for scans with mesoscale resolutions, where no single fibers, only coarser structures like strands and fiber bundles, respectively, are visible [43,123]. Since microscopy shows that during compression molding the fibers within one strand mostly flow and orient together, deforming ye<sup>t</sup> remaining as an intact strand with locally highly aligned fibers, the density gradients of CT scanned ROS-based materials are su fficient for a determination of local average strand orientations even at a coarser scan resolution [4,42,43,123]. In a CT scan of a ROS-based material the smallest density gradient is present in fiber direction, intermediate density di fferences are visible transverse to the strand orientation, and the highest density gradient occurs normal to the strand plane and at strand–strand boundaries (matrix rich strand interfaces), respectively [4,20,21,42,43,123]. Due to these intra- and inter-strand density changes the determination of strand orientations from mesoscale resolution CT scan data is feasible [4,20,21,43].

Denos et al. use a CT scan with a mesoscale resolution of 53 μm (voxel edge length) to determine the heterogeneous internal microstructure of a prepreg carbon/polyetherketoneketone (PEKK) UD strand-based T-bracket part with maximum dimensions of 65 × 65 × 45 mm [20,42,43]. However, at that resolution and with the applied CT scan settings it is not possible to distinguish between single carbon fibers ( ∅ ≈ 5–7 μm) and the matrix or to discern strand boundaries (~100 μm thick) in thickness direction [42]. Although Denos' CT scan configuration is not able to represent distinct strand boundaries, it is still possible to receive a mean local fiber orientation due to su fficiently pronounced relative density gradients [20,42].

Another common method to achieve bigger CT scan volumes at a reasonable resolution is to merge several scan volumes generating a digital twin of the scanned part and its microstructure. Sommer et al., Kravchenko, and Denos et al. merged 8 partial scans of a prepreg carbon/PEKK UD strand-based tensile test specimen with a size of 30 × 30 × 5 mm each, using a scan resolution of 15 μm [17,123,131]. At that resolution the CT scan quality is high enough to discern between strands and suitable for precise fiber orientation analyses. An analysis mesh size of 0.7 × 0.7 × 0.1 mm is used to determine a single orientation tensor from each of the measured orientation vectors by a grayscale analysis [123,131]. The finer analysis resolution better resolves the thin strands in thickness direction and enables gathering more detailed information about the local strand orientation changes.

The spectrum generated by the X-ray source significantly depends on the elemental composition of the used target material. Tungsten (W) is widely used as target material for microfocus X-ray sources. However, depending on the applied X-ray voltage and the absorption behavior of the sample material, alternative target materials might deliver beneficial spectrum characteristics that can improve CT measurement quality concerning the separation capability of fiber and matrix for fiber orientation analysis [132]. A higher μ-CT scan quality, by means of a higher contrast between fiber and matrix, also enables to scan bigger composite part volumes, fulfilling industry demands, where the fiber orientation within a whole component is of interest.

#### **2. Material and Methods**

## *2.1. Compression Molding*

The high-performance CF-SMC used in this study is Hexcel's HexMC ® (Hexcel Corporation, Stamford, CT, USA). This DFC material is designed for compression molding of complex 3D shaped parts in a heated metal tool. HexMC consists of prepreg carbon/epoxy UD tapes that are longitudinally slit and transversely chopped into strands and then randomly distributed into a mat [14]. Those ROS have nominal in-plane dimensions of 50 × 8 mm and a thickness of approximately 0.15 mm (Figure 3). The strands contain high strength carbon fibers impregnated by fast-curing Hexcel HexPly ® M77 epoxy resin [14,133]. The carbon fiber content amounts to 62% by weight in the raw material, corresponding to 57% fiber volume content in a molded part and giving a material density of 1.55 g/cm<sup>3</sup> (Table 1) [14,15].

For the compression molding trial a research tool with a ribbed hat profile tool insert designed for compression molding of CF-SMCs is used (Figure 4a). The tool is made to mold a complex ribbed structure with ribs of varying heights and non-symmetrically alternating wall thicknesses with an attached plate (Figure 4b). The molded part has outer dimensions of 450 × 450 mm and a nominal wall thickness of 2 mm. The ribbed hat profile spans an area of 150 × 450 mm and has a hat height of 52 mm.

The ribbed hat profile parts analyzed in this study were manufactured with a mold coverage of approximately 80% using a 1000 ton compression molding machine (Die ffenbacher DCL-S 1000, Die ffenbacher GmbH, Eppingen, Germany) at a temperature of 140 ◦C, a pressure of 200 bar, a closing speed of 5 mm/s, and a closing time of 480 s (Table 2). The mold coverage would be high for standard SMCs, ye<sup>t</sup> for ROS-based high-performance CF-SMCs it is rather low giving the material a su fficiently long flow path to create a flow-induced fiber mesostructure especially in the hat end brim section, which is wanted in this study.

Material flow simulations in 3D TIMON CompositePRESS were used to develop an SMC charge pattern (Figure 5), which allows for a 100% filling of the complex ribbed structure without defects, such as FMS. Predicted pressure distributions, shear rates, flow velocities and flow front advancements are analyzed for various charge patterns with the goal to reach uniform filling of the cavity. Due to the slightly asymmetrical design of the ribbed structure with very thin (0.8 mm) and very thick ribs (8.2 mm) of different rib heights, the cavity volume in the hat profile area changes locally. Therefore, the finally identified charge pattern used in this study consists of several smaller single charge packages optimized in size and position in order to achieve a balanced raw material distribution in the cavity and to facilitate optimal rib filling. Furthermore, the edges of some charge packages are placed underneath some rib entries to ease the flow of the UD strands into the ribs. All smaller charge packages are stacked on one 440 × 440 mm base layer to ease the charge placement in the tool.

In order to reach a lower viscosity for optimal flow behavior, the material charge was pre-heated for 17 s. This pre-heating time was found in a previous flowability study for HexMC plates. For the pre-heating procedure the charge is placed in the bottom mold half and the upper mold half is fast lowered to the pre-heating position leaving a small gap between the cavity wall and the material surface. Heat conduction from the lower mold half and heat radiation from the upper mold half decrease the material viscosity enabling better flowability when the mold is closed after the pre-heating phase.

**Figure 3.** (**a**) HexMC raw material roll; (**b**) HexMC mesostructure with randomly in-plane oriented prepreg carbon/epoxy unidirectional (UD) strands.



1 layer (220 x 50 mm)

1 layer (170 x 50 mm)

**Figure 4.** (**a**) CAD images of the upper and lower mold half of a research tool with a ribbed hat profile tool insert designed for compression molding of carbon fiber sheet molding compounds (CF-SMCs); (**b**) CAD geometry of the ribbed hat profile part with varying rib heights and thicknesses.

 molding processing conditions.

**Table 2.**

> 5

layers (50 x 50 mm)

layers (110 x 40 mm)

5

Compression


2 layers (330 x 200 mm) 1 base layer (440 x 440 mm) 2 layers (330 x 40 mm)1 layer (220 x 50 mm) 5 layers (50 x 50 mm) 1 layer (170 x 50 mm) 5 layers (110 x 40 mm)

**Figure 5.** HexMC charge pattern for ribbed structure part; (**a**) schematic illustration with dimensions and number of layers per charge package (for visualization purposes the charge package are slightly separated in thickness direction); (**b**) prepared HexMC charge on a metal preform.

(**a**) (**b**)

#### *2.2. Short Shots*

For accurate fiber orientation predictions with DFS tools it is crucial to obtain as precise filling predictions as possible. In order to receive a first understanding about the simulation quality, short-shot experiments were conducted with a shimmed mold for a comparison of the true material flow behavior and the numerical flow front advancement inside the mold cavity. For the shimming technique, two different shims are used as spacers between both mold halves preventing a complete mold closure but allowing the material to undergo the real pressure applied in the unshimmed mold. The used spacers have dimensions so that the press contacts the shims at intermediate positions leaving a gap between the mold halves of 8.65 and 4.00 mm, respectively. The schemes in Figure 6 illustrate the shimming technique. For the short shot trials the material charge pattern and the molding parameters remain the same as for the normal press experiments (cf. Figure 5 and Table 2).

**Figure 6.** Schematic illustration of the shimming technique for short shot experiments in a cross-sectional view; (**a**) initial pre-heating position of the upper mold half; (**b**) final position of the upper mold half touching the 8.65 mm shims and showing the short shot of the pressed HexMC charge.

#### *2.3. Compression Molding Simulation with 3D TIMON CompositePRESS*

## 2.3.1. Filling Simulation

The 3D TIMON software module CompositePRESS is used to simulate the press process of the ribbed structure. It calculates the cavity filling and the resulting fiber orientations for the thermoset sheet-shaped HexMC charge. In 3D TIMON CompositePRESS there are two analysis methods available to simulate the compression molding process—the 'Euler method' and the 'morphing method' [105]. The Euler method uses a multi-layer voxel mesh to discretize the part by several elements over the thickness, whereas the morphing method utilizes a one-layer tetra mesh to represent the final part geometry (mold-closed shape). Advantages of the Euler method are the ease of 3D expressions for the charge definition and its flow as well as the possibility to discretize the part by several elements over the thickness. The accuracy of the movement, rotation and bending prediction of each single fiber improves when the number of elements across the thickness increases. However, this ultimately leads to a huge number of elements for more complex and/or larger components. Furthermore, due to the

cubic-shaped voxel elements, a 3D-shaped component cannot be meshed without terraced steps in the contour.

In order to calculate the material flow during the compression molding process with the morphing method, the mesh is initially morphed in negative pressing direction, which means that the tetra elements are artificially stretched in z-direction to represent the open mold condition (cf. Figure 7). During the mold closure the morphed elements are then pressed back into their initial form (mold-closed shape) enabling the material flow calculation. The biggest advantage of the morphing method is that a tetra mesh can be used. Tetra elements allow easy meshing of complex geometries without terraced steps. This enables a high quality mesh consisting of a much smaller number of elements than a comparable voxel mesh. Therefore, the morphing method allows for faster analyses than the Euler method.

**Figure 7.** Mold filling simulation in 3D TIMON CompositePRESS using the morphing method; (**a**) cross sectional view of the morphed tetra mesh in the initial open-mold position; (**b**) pressed back tetra mesh close to end of compression showing the proceeding flow front in the hat brim.

Considering the size and complexity of the ribbed part, the morphing method is the method of choice in this study since a high-quality tetra mesh without terraced steps can be used with a reasonable calculation time. As the applied Light 3D solver virtually divides each tetra element into 20 calculation points across the thickness, realistic flow predictions are possible even with a one-layer mesh. For the discretization of the entire part 141,385 tetra elements are used. Element sizes are varied to reduce the total amount of elements. Simple geometry areas of less interest, like the plate area, are meshed coarser with a maximal element edge length of 14 mm, whereas the hat area is discretized finer with a minimal element edge length of 0.4 mm in some rib fillets (cf. Figure 7) giving precise prediction results at a still reasonable calculation time. In total, 203 output steps are defined in order to gain a finely resolved filling analysis. Each charge package used in the physical moldings trials is also defined in 3D TIMON CompositePRESS by selecting all elements in the respective charge area and assigning the real charge thickness to those elements. Stacked charge packages are considered as one combined volume so that the actual thickness of each charge stack can be taken into account by the 20 virtual calculation points across the thickness of the stretched mesh (open-mold position). This enables a realistic modeling of the total charge volume.

Besides the chosen curing and viscosity models, good filling simulation quality depends on the material properties deposited in the material card. The interaction between the fibers and the fluid flow (two-way coupling) and the resulting anisotropic viscosity is not considered. Typical polymer characterization techniques, such as viscosity measurements by means of rotational viscometers, differential scanning calorimetry (DSC), and dielectric analysis (DEA) measurements are used in order to calibrate the HexMC material parameters for the applied Cross–Castro–Macosko viscosity model (Table 3) and to find fitting Kamal parameters (Table 4) to model the curing behavior appropriately.

Figure 8 shows three DSC curing curves at relevant molding temperatures of 150 ◦C, 140 ◦C, and 130 ◦C, one exemplary DEA measurement at 145 ◦C, and the corresponding fitted Kamal model

curves. Additionally, three Kamal model curves at 155 ◦C, 135 ◦C, and 125 ◦C are plotted in gray to show the modeled curing behavior at slightly varied temperatures. The generated Kamal model curves fit well to the experimentally observed curing behavior of HexMC.


**Table 3.** Cross–Castro–Macosko viscosity model constants.


**Table 4.** Kamal curing model constants.

**Figure 8.** Isothermal differential scanning calorimetry (DSC) measurements of curing HexMC samples at different temperatures and fitted Kamal model curing curves using the identified parameters given in Table 4.

#### 2.3.2. Direct Fiber Simulation (DFS) with Strand Generation Feature

A novel feature recently implemented in 3D TIMON CompositePRESS enables the simulation of ROS materials. In the pre-processing the new strand generation feature generates randomly oriented multi-bundle UD strands mimicking the real initial material charge. In contrast to the DBS of Meyer et al. [73], which uses one chain of truss elements per fiber bundle, 3D TIMON CompositePRESS represents strands by a grouping of several aligned numerical fibers (cf. Figure 9). Four fibers frame the outer strand dimensions and a user-defined number of additional inner fibers can be added inside the strand volume to represent a realistic amount of fiber bundles within a strand. Each fiber is

still divided into a user-defined number of divisions (rods) giving the fibers the possibility to bend naturally in the flow. If it is assumed that a single fiber represents one fiber bundle within a strand, the simulation approach can be considered to be representative and gives the opportunity to mimic the strand geometry and character properly. Nonetheless, so far all fiber bundles defining one strand are numerically handled separately without a cohesive force between them.

**Figure 9.** HexMC carbon fiber UD strands and a schematic illustration of their numerical counterpart modeled in 3D TIMON CompositePRESS with 4 boundary fibers defining the outer strand dimensions and 11 added inner fibers (here each fiber consists of 20 nodes and 19 rods, indicated on just one boundary fiber) (note: for visualization purposes length specifications are not true to scale).

More fiber divisions allow a more realistic bending behavior, especially in situations where fibers flow into small features such as ribs. Regarding the calculation effort, the number of divisions as well as the number of additional fibers should be carefully chosen. However, the more fibers are calculated, the better the simulation accuracy and the better the fiber volume fraction representation. The numerical strands have dimensions of 50 × 8 × 0.15 mm, equal to the real HexMC UD strands. For this study each strand is represented by 4 boundary fibers and 11 inner fibers, all divided into 19 rods connected by 20 nodes (cf. Figure 9). Using 15 fibers to represent the designated strand volume is chosen as a compromise between the ability to visually analyze the strands' movements and deformations during the DFS and the needed calculation time.

For the DFS of the ribbed structure a random in-plane strand orientation distribution in the 3D-shaped initial charge is generated (Figure 10a). The parameters for the strand generation algorithm are set so that one strand is generated per mesh element and that the minimum distance between initial fibers is zero, which means that fibers can be generated close to each other. The strand generation approach does not model the strands as enclosed volumes, so that strands intersect each other and compaction is not considered. In Figure 10b the numerically defined individual charge packages are displayed in different colors. Both the macroscopic shape of the generated charge and the distribution of strands are comparable to the physical HexMC. Since the strands are cut at each edge of a charge package, it is possible that incomplete and/or shorter virtual strands are generated—just like in reality. Even in the almost vertically oriented cavity zones the numerical representation of the initial charge pattern is reasonably realistic when compared to the true charge configuration (cf. Figure 5b). As only straight fibers are generated, a negligible amount of fibers is cut in order to fit into the mesh in tight radii. Initially, the generated strands are not bent in order to follow the curvature of the charge in those regions.

**Figure 10.** (**a**) Numerically modeled initial charge pattern for the ribbed hat profile part consisting of in-plane randomly oriented UD strands (5 UD strands are highlighted); (**b**) bottom view of the numerical charge definition where each charge package is colored differently for better discernibility.

#### *2.4. Micro-CT Measurements*

Three entire ribbed hat profile parts were μCT scanned with a CT-AlphaDuo device (Procon X-Ray GmbH, Sarstedt, Germany) operated by the Fraunhofer WKI, Hannover, Germany. The system is equipped with a 240 kV microfocus X-ray tube XWT-240-TCHE Plus (X-RAY WorX GmbH, Garbsen, Germany) with a high energy diamond/tungsten transmission target offering a JIMA (Japan Inspection Instruments Manufacturers' Association) test resolution of 0.9 μm. The X-rays are recorded by a 249 × 302 mm PaxScan® 2530DX detector (Varian Medical Systems, Inc., Salt Lake City, UT, USA) with a resolution of 1792 × 2176 pixels. In order to maximize the voxel resolution at the given sample diameter, the detector panel width was virtually extended by measuring field extension (MFE) (in situ horizontal movement of the detector panel). The focus object distance (FOD) and the focus detector distance (FDD) were 642 mm and 1500 mm, respectively. For best possible scan results, the parts' plates were removed from the hat section with a water jet cutter. This improves the CT scan quality since the penetration length is shorter and thereby induced scan defects (e.g., the resulting image noise) are kept low. Each ribbed part was mounted on a rotating table between X-ray source and detector (see Figure 11a). The whole scan setup is built on an air-damped granite base in order to minimize vibration-induced scan artifacts. For the measurements the X-ray tube was operated at a voltage of 160 kV and a current of 250 μA. The sample was rotated in 2400 rotation steps (0.15◦ angular increments) for a full 360◦ rotation with an integration time of 6 × 200 ms per angular increment. This was repeated for both detector positions. The measuring time for each partial scan accumulated to approximately 100 min.

All samples were scanned at a mesoscale resolution of 59.6 μm (voxel edge length) corresponding to approximately 16.8 pixels/mm. Each voxel contains a single 16-bit relative density-defined grayscale value between 0 and 65,535. With the applied scan resolution, which is approximately 10× the carbon fiber diameter, it is possible to clearly see in-plane strand boundaries. However, the strands are only resolved with ~2.5 voxels in thickness direction (minimum fiber bundle dimension), which is not enough for clear boundary detections in that strand dimension. Further, a 1 mm thick aluminum filter is used in front of the X-ray target in order to reduce the effect of beam hardening. As a result of different penetration lengths of low and high energy X-ray photons beam hardening causes grayscale value gradients inside the sample, which can especially affect grayscale value sensitive analyses such as the fiber orientation analysis. By using a filter, low-energy X-ray photons are suppressed and thereby the X-ray spectrum is shifted to higher energies. In previous experiments this specific CT scan configuration has proven to be conducive for reasonable fiber orientation analyses of larger CF-SMC components. Reconstruction of the 3D volumes was performed with the software 'Offline CT' (Fraunhofer Institute for Integrated Circuits IIS, Erlangen, Germany). Table 5 summarizes the applied μCT scan parameters.

**Figure 11.** (**a**) Exemplary CT measurement setup (note: here, focus object distance (FOD) and the focus detector distance (FDD) are not the ones used for the scans); (**b**) schematic illustration of merged CT scans with overlap areas.


**Table 5.** Micro-CT measurement parameters.

To capture the entire ribbed structure at the necessary resolution, seven partial scans were performed and afterward virtually stitched together using the CT scan visualization and processing software VGSTUDIO MAX 3.3. As the individual CT scans slightly vary in their grayscale value distribution, for each 3D data set the average grayscale values of background and material were determined from the grayscale value histograms. Based on these values the volumes were imported into VG applying the 'histogram calibration' function to create a homogeneous grayscale value profile over all scans. Each individual scan captures a sample section of about 150 mm in width, 52 mm in height, and 90 mm in axial direction (Figure 11b) leading to an overlap of at least 15 mm between adjacent scans that eases the manual alignment and merging of the imported partial volumes. The alignment and merging procedure was done with utmost care so that any negative effect on the final 3D data set and the subsequent CT scan analysis is ruled out. Regions of interest (ROIs) were created defining the area of two adjacent volumes to be merged. The filter 'merge volumes' in mode 'mean' was applied to combine two adjacent volumes based on these ROIs. By repeating this procedure, one continuous 3D volume of the entire sample was composed of the partial scans.

The surface determination was conducted on the merged volume so that the scanned parts' thicknesses match with the average true thicknesses of the molded parts measured with an electronic external measuring gauge. The final object volume is then used as ROI where the background (surrounding air) is removed. This reduces the 3D volume data to be analyzed leading to faster fiber orientation calculations. The final volume of each stitched ribbed part has outer dimensions of 150 × 450 × 52 mm and an analyzed volume of approximately 350,000 mm<sup>3</sup> captured by almost 1.78 billion voxels.

All three scanned ribbed structures are analyzed regarding their second order fiber orientation tensors (hereafter FOTs) in several sample areas using the fiber composite material analysis (FCMA) tool implemented in VG. The same tetra mesh that is used for the 3D TIMON CompositePRESS process simulation is imported to VG and is utilized as an integration mesh for the CT scan analysis. Therefore, each 3D volume is registered to the tetra mesh's coordinate system, so that the scanned ribbed structures are perfectly superimposed with the tetra mesh. With this procedure a one-to-one comparison with the process simulation results is possible since the same tetra elements used in the process simulation can be utilized for the fiber orientation analyses in VG. For each element of this analysis mesh the FOTs, eigenvectors, and eigenvalues are determined and evaluated.
