*Article* **The Impact of Fiber Orientation on Structural Dynamics of Short-Fiber Reinforced, Thermoplastic Components—A Comparison of Simulative and Experimental Investigations**

**Alexander Kriwet 1,\* and Markus Stommel <sup>2</sup>**


**Abstract:** The quality of the fiber orientation of injection molding simulations and the transferred fiber orientation content, due to the process–structure coupling, influence the material modeling and thus the prediction of subsequently performed structural dynamics simulations of short-fiber reinforced, thermoplastic components. Existing investigations assume a reliable prediction of the fiber orientation in the injection molding simulation. The influence of the fiber orientation models and used boundary conditions of the process–structure coupling is mainly not investigated. In this research, the influence of the fiber orientation from injection molding simulations on the resulting structural dynamics simulation of short-fiber reinforced thermoplastic components is investigated. The Advani–Tucker Equation with phenomenological coefficient tensor is used in a 3- and 2.5-dimensional modeling approach for calculating the fiber orientation. The prediction quality of the simulative fiber orientations is evaluated in comparison to experiments. Depending on the material modeling and validation level, the prediction of the simulated fiber orientation differs in the range between 7.3 and 347.2% averaged deviation significantly. Furthermore, depending on the process–structure coupling and the number of layers over the thickness of the model, the Kullback–Leibner divergence differs in a range between 0.1 and 4.9%. In this context, more layers lead to higher fiber orientation content in the model and improved prediction of the structural dynamics simulation. This is significant for local and slightly for global structural dynamics phenomena regarding the mode shapes and frequency response behavior of simulative and experimental investigations. The investigations prove that the influence of the fiber orientation on the structural dynamics simulation is lower than the influence of the material modeling. With a relative average deviation of 2.8% in the frequency and 38.0% in the amplitude of the frequency response function, it can be proven that high deviations between experimental and simulative fiber orientations can lead to a sufficient prediction of the structural dynamics simulation.

**Keywords:** short-fiber reinforced thermoplastic components; injection molding simulation; fiber orientation; structural dynamics; material modeling

### **1. Introduction**

Short-fiber reinforced thermoplastics are an essential group of engineering materials in modern vehicle powertrains due to their significant lightweight potential. Under Noise-Vibration-Harshness (NVH) aspects, short-fiber reinforced plastics offer good vibration isolation and thus noise isolation behavior. This is due to the favorable stiffness and damping behavior. For efficient prediction of the stiffness and damping of short-fiber reinforced plastics, established material models are based on multi-stage homogenization methods. The principle of all methods is that the microstructure and thus the properties of the composite are described with mathematical–physical models. Thereby, the consideration of

**Citation:** Kriwet, A.; Stommel, M. The Impact of Fiber Orientation on Structural Dynamics of Short-Fiber Reinforced, Thermoplastic Components—A Comparison of Simulative and Experimental Investigations. *J. Compos. Sci.* **2022**, *6*, 106. https://doi.org/10.3390/ jcs6040106

Academic Editor: Stelios K. Georgantzinos

Received: 9 March 2022 Accepted: 30 March 2022 Published: 1 April 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

the fiber orientation via the process–structure coupling is fundamental. By conditioning the material models with direction-dependent properties, discrete material parameters are assigned to each element of the simulation depending on the local fiber orientation tensor. Therefore, the quality of the prediction of the fiber orientation in the injection molding simulation influences the prediction quality of the structural dynamics simulation. Further information corresponding to the homogenization methods and the process–structure coupling is provided in [1–8]. Glaser shows in [9] that the prediction quality of structural dynamics simulations of a short glass fiber reinforced thermoplastic intake pipe can be significantly increased by considering the fiber distribution and thus the anisotropic stiffness and damping. The main focus of this study is on the procedure and advantages of the process–structure coupling. However, information and boundaries to calculate the fiber orientation are not provided. Studies by Schmachtenberg et al. [10] additionally use advanced experimental methods to calibrate the simulation models. Arping [11] and Kremer [12] use a similar method to homogenize the properties of the fiber-reinforced plastic component and recalibrate the models of the structural dynamics simulation through reverse engineering. Thereby, material modeling is the focus of the research. The simulative fiber orientation is considered a fixed boundary condition. The disadvantage is that errors caused by an insufficient simulative fiber orientation are compensated by the reverse engineering of the material model. Influences on the results caused by the quality of the simulative fiber orientation are largely not considered. State-of-the-art extended approaches additionally pursue the consideration of the material properties depending on the boundary conditions, such as frequency, temperature or humidity. There exists a large number of publications dealing with the structural dynamics design of fiber-reinforced plastic components under NVH aspects [13–18]. However, all existing publications have in common that the reliable prediction of the fiber orientation in the injection molding simulation is assumed. The boundary conditions and fiber orientation models used by the injection molding simulations are mainly not explained. However, there can be significant differences in the prediction quality of the fiber orientation, depending on which fiber model is used and how the models are calibrated. Furthermore, as a result of the process–structure coupling, there is always a loss of information about the fiber orientation between the injection molding model and the finite element (FE) model, which is not sufficiently considered in existing publications [1].

The aim of this research is to investigate the influence of the fiber orientation from the injection molding simulation on the resulting structural dynamics FE simulation of short-fiber reinforced thermoplastic components. In the first section of this paper, the state of art and the applied method for calculating the fiber orientation using injection molding simulations are explained. Furthermore, the method used to investigate the influence of the fiber orientation on the prediction quality of the structural dynamics simulation is presented. In the last section, the results of the simulative and experimentally investigated fiber orientation, as well as the prediction quality of the corresponding structural dynamics investigations, are discussed.

### **2. State of Art and Methods**

In the first part of this section, the state of art for fiber orientation calculation with injection molding simulations is discussed. Following the proposed method for the simulative and experimental characterization of the fiber orientation is presented. In the last part, the focus is on the method for the experimental characterization by modal analysis and the structural dynamics simulation model of the short-fiber reinforced plastic components. In this research, the proposed methods are implemented for two types of short-fiber reinforced thermoplastics. On the one hand, polyamide 6.6 (PA66) is a common engineering plastic, and on the other hand, polyphthalamide (PPA) is a high-performance plastic, each with 50 wt.% short glass fiber reinforcement (GF50) [19,20]. As a representative composite component for investigating the influence of fiber orientation on structural dynamics, the so-called engine bracket is investigated. In modern combustion engines, the engine bracket

is mounted on the crankcase and transmits the powertrain-induced operating vibrations into the vehicle structure. For further information about the functionality of the engine bracket, please refer to [15–18,21,22].

### *2.1. State of Art of Injection Molding Simulations for Fiber Orientation Calculation*

The prediction of the fiber orientation in injection molding simulations is based on the so-called continuity equation according to Fokker-Planck [23,24]. This is derived from the velocity field of the fluid and thus the hydrodynamic forces acting on the fiber. Extended by a diffusion term to describe the fiber interaction in the fluid field, the Folgar–Tucker equation [25] is formed:

$$\frac{d\psi}{dt} = -\nabla \cdot (\dot{p}\psi) + D\_i \nabla^2 \psi,\tag{1}$$

*ψ* corresponds to the fiber probability density function (PDF) at time *t*, ∇ is the gradient operator, . *p* is the change of the fiber orientation and *Di* is the fiber interaction coefficient. The first part of the equation can be summarized as a hydrodynamic term and the other as a diffusion term. Applied to injection molding simulations of composite components, solving Equation (1) is numerically extremely cost-intensive. For this reason, the PDF is mainly substituted by an evolution equation of the fiber orientation tensor **A**:

$$\frac{d\psi}{dt} \approx \frac{d\mathbf{A}}{dt} = \dot{\mathbf{A}} = \dot{\mathbf{A}}^h + \dot{\mathbf{A}}^d \,. \tag{2}$$

with

$$\mathbf{A} = \mathbf{A}\_N = \int\_0^{2\pi} \int\_0^{\pi} \mathbf{p}\_N \boldsymbol{\uppsi}(\Theta, \Phi) \sin \Theta d\Theta d\Phi. \tag{3}$$

Thus, **. A** *h* corresponds to the hydrodynamic part and **. A** *d* to the diffusion part of the Folgar–Tucker Equation. Equation (2) is numerically stable for calculating the fiber orientation of plastic components in injection molding simulations. On this basis, a large number of publications exist which follow different approaches for the formulation of the hydrodynamic and diffusive parts. Advani and Tucker formulate in [26] an extension of the Fokker–Planck Equation to consider the rate of change of the second-order fiber orientation tensor **A**, which is calibrated by a phenomenological approach via the parameter *CI*. Huynh shows in [27] that the prediction of the calculated fiber orientation can be further improved by a scalar reduction factor κ. Wang et al. [28] extend this reduction factor by a reformulation of the second-order fiber orientation tensor **A** through a decomposition of the eigenvalues λ and eigenvectors **e**. This forms the reduced strain closure (RSC) model, which is an established method for calculating fiber orientation. A comparable established approach is shown by Phelps and Tucker in [29], whereby the phenomenological parameter *CI* is replaced by a rotary diffusion tensor **C**. This forms the so-called anisotropic rotary diffusion (ARD) model and thus the basis of advanced models [29–32]. An overview of macroscopic fiber orientation models is discussed in [24].

Furthermore, the decomposition of the second-order fiber orientation tensor *<sup>d</sup>***A***<sup>N</sup> dt* is dependent on the fiber orientation tensor of the next higher-order A*N*:

$$\frac{d\mathbf{A}\_N}{dt} = f(\mathbf{A}\_N). \tag{4}$$

Established methods use the formulation of a closure approximation. This approximation of the higher-order fiber orientation tensor A*<sup>N</sup>* is based on mathematical assumptions [23,24,33]. In general, the increased information of the higher-order tensors is usually not fully captured by the closure formulations. Advani and Tucker show in [26] a linear summation of all products from the components of the second-order fiber orientation tensor **A***ij* and the Kronecker delta *δij* to approximate the higher-order tensor. By neglecting the linear terms, the quadratic closure is formed and by combining the linear and quadratic approach, the hybrid closure is formed [26]. Furthermore, the so-called natural

closure of Verleye and Dupret should be mentioned [34], where the fourth-order orientation tensor A*ijkl* is defined as a function of the second-order tensor *<sup>f</sup>* **A***ij* . This forms the basis of various advanced approximations methods [24,35]. Advanced models pursue an exact formulation of the fourth-order fiber orientation tensor [36,37], the use of neural networks [38,39], distribution reconstruction methods [23,33,40,41] or the reconstruction of even higher-order fiber orientation tensors [33,42]. However, to ensure numerically costeffective fiber orientation calculations in injection molding simulation, hybrid or natural closures are still widely used today [24]. An overview of existing closure formulations is discussed in [23,24,33].

### *2.2. Methodology for Simulative and Experimental Fiber Orientation Investigations*

According to the state of art, the Advani–Tucker Equation is used as the basis for calculating the fiber orientation in this research. The hydrodynamic part is defined as:

$$\dot{\mathbf{A}}^{\dot{h}} = (\mathbf{W} \cdot \mathbf{A} - \mathbf{A} \cdot \mathbf{W}) + \zeta (\mathbf{D} \cdot \mathbf{A} + \mathbf{A} \cdot \mathbf{D} - 2\mathbf{A} : \mathbf{D}),\tag{5}$$

and the diffusion part as

$$\dot{\mathbf{A}}^d = 2\dot{\gamma} \cdot \mathbf{C}^v \cdot (I - A),\tag{6}$$

with **<sup>W</sup>** as vorticity tensor, *<sup>ξ</sup>* as particle shape function, **<sup>D</sup>** as strain rate tensor, . *γ* as the magnitude of the strain rate tensor and *I* as a unit tensor. *C<sup>v</sup>* corresponds to a phenomenological coefficient tensor for describing the fiber interaction and is calibrated through experimentally determined fiber orientations. Thus, the diffusion part is equivalent to a rotary diffusion approach by the appropriate definition of *Cv*. To evaluate and calibrate the phenomenological parameters of the injection molding model, μCT investigations are carried out. For this purpose, material plates are created by varying the geometry or the process parameters, and μCT specimens are investigated at selected positions according to the fiber orientation. It is important to choose specimen positions that allow a representation of the microstructure and thus a reliable investigation of the fiber orientation. Under the usage of simple plate geometry, a small number of μCT specimens are taken in an evenly distributed way. When transferred to the injection molding simulation of short-fiber reinforced plastic components, significant differences in the prediction quality of the fiber orientation can occur [4,8,21]. As a result, extended methods include a recalibration of the default parameters of the simulation in comparison to μCT specimens taken from plastic components [8,43]. At the component level, μCT specimen positions are selected with significantly different geometric or material-specific characteristics in order to calibrate the simulation for robustness. The geometric aspects include, for example, different wall thicknesses, ribs or triple points. This allows a direct recalibration of the injection molding model. On the other hand, material-specific aspects are, for example, inor ejection positions or impact points of the plastic melt. These aspects allow a recalibration of the injection molding model in case of process-related errors or damage analyses. Based on the experimentally determined fiber orientations, the simulation model parameters are calibrated with numerical or mechanistic methods. The injection molding material model is sufficient if the deviation between simulated and experimentally determined fiber orientation is minimized for several parameter variations [8,43]. Figure 1 schematically shows the procedure for a multi-stage calibration of the calculated fiber orientation from injection molding simulations.

In this research, injection molding simulations are performed using two commercial software programs. On the one hand, 2.5-dimensional (2.5D) injection molding simulations were conducted with the software CADMOULD®from SIMCON. The simulations were conducted in the scope of contract simulations by PART Engineering respectively. Thereby, the hydrodynamic part of the Advani–Tucker Equation equals a natural closure (NC) formulation to approximate the fourth-order tensor A*ijkl*. Furthermore, the fiber interaction tensor *C<sup>v</sup>* = *CNC* in the diffusion part depends on the alignment factor in the middle layer

*αK*, in the boundary layer *α<sup>R</sup>* and on the rotational velocity factor *βR*. This corresponds to the phenomenological coefficients in CADMOULD®, which are calibrated in comparison to μCT results. The initial coefficient tensors *CNC*,*<sup>S</sup>* 0,*CM* for the respective material *S* originating from a database of the CADMOULD® software. These are defined as default values for both the PA66GF50 and the PPAGF50 material as *α<sup>K</sup>* = 0.92, *α<sup>R</sup>* = 0.92, *β<sup>R</sup>* = 0.15. To maximize the calculated information of the fiber orientation from the injection molding simulation, 10 elements are defined uniformly over each wall thickness. In comparison, 3-dimensional (3D) injection molding simulations were conducted with the software MOLDFLOW® from AUTODESK. The simulations were conducted in the scope of contract simulations by Daimler Trucks AG respectively. Here, an orthotropic closure (OC) is defined in the hydrodynamic part and the coefficient tensor *C<sup>v</sup>* = *COC* depends on the scalar interaction coefficient *CI* and the asymmetric coefficients of the rotational diffusion *D*1, *D*2, *D*3. This corresponds to the phenomenological coefficients in MOLDFLOW®. The initial coefficient tensors *COC*,*<sup>S</sup>* 0,*MF* originated from a database of the company BASF SE. These are defined as default values *CI* = 1.0, *D*<sup>1</sup> = 1.0, *D*<sup>2</sup> = 0.8, *D*<sup>3</sup> = 0.15. To maximize the fiber orientation information from the injection molding simulation with MOLDFLOW®, 14 elements are defined uniformly over each wall thickness.

**Figure 1.** Method for calibrating the calculated fiber orientation of injection molding simulations in comparison to experimental μCT investigations.

To evaluate the prediction quality of the simulative fiber orientation, experimental μCT investigations are evaluated in comparison at the plate and component level, as seen in Figure 2. At the plate level, the μCT specimens are extracted equally distributed from the plate geometry. The dimensions (length × width × thickness) of the plate geometry are (180 mm × 180 mm × 2 mm). In comparison, μCT specimens at the component level are extracted with a focus on geometric aspects of the engine bracket. As a result, a μCT specimen is extracted from a side rib and the upper end of the bracket. The approximate dimensions of the engine bracket are (215 mm × 135 mm × 135 mm) with a thickness of the surfaces and ribs between 2–5 mm. In the current case, the dimensions of the μCT specimens both at the plate and component level are approximate (5 mm × 5 mm × 2 mm). Figure 2 provides an overview of the geometry on the plate and component level and shows the specimen positions. For evaluation, A11 corresponds to the component of the second-order orientation tensor **A***ij* in the main flow direction. A22 corresponds to the component transverse to the main flow direction and A33 to the component in the thickness direction. Figure 2 schematically shows the evaluation directions of the components of the second-order orientation tensor **A***ij*. The μCT investigations are performed with the 3DnanoCT Phoenix Nanotom from GE Sensing Technologies GmbH. The fiber orientation is evaluated using VGSTUDIO MAX software from Volume Graphics GmbH. To compare the

simulative and experimentally determined fiber orientations, a process–structure coupling to FE models of the geometry of the μCT specimens is performed. Regarding the injection molding simulations, a discretization of 10 elements through the wall thickness of the FE model is defined to ensure comparability. The process–structure coupling is performed with the software CONVERSE® from PART Engineering GmbH. Figure 3 shows the proposed approach for evaluating and recalibrating the simulative fiber orientation from injection molding simulations in comparison to the experimental μCT results.

**Figure 2.** Geometry of the plate (**a**) and component level (**b**) and corresponding μCT specimen positions (P1–P5, B1–B2). Evaluation directions of the fiber orientation (**c**).

**Figure 3.** Proposed approach for calibrating and comparing the calculated fiber orientation of injection molding simulations to experimental μCT investigations.

### *2.3. Methodology for Experimental and Simulative Structural Dynamics Investigations*

To investigate the structural dynamics of the short-fiber reinforced plastic components, modal analyses are carried out. Excitation of the system forces a structural response at the resonant frequency. This is characterized by transfer behavior. The response function represents the relationship between the out and input signals of the system. According to existing investigations, modal analyses are carried out with an elastic bearing of the plastic components [21]. The aim of the elastic bearing is that the impact effects of the surrounding system on the plastic component are minimized. Thus, the stiffness and damping of

the bearing are defined, that there is no superposition on the resonance frequencies and amplitude of the response function of the plastic component. By applying the plastic component with bearing in a climate chamber, the temperature and humidity can be set reproducible.

The basis for the simulation of the structural dynamics of short-fiber reinforced, thermoplastic components is the material model of stiffness, damping and viscoelasticity. According to existing investigations, the Arbitrary-Reconsidered-Double-Inclusion (ARDI) material model is applied [4]. This corresponds to a two-stage homogenization method of the stiffness and damping from the properties of the composite. In the first step, the effective stiffness and damping are calculated from the properties of the fibers, the matrix and the interphase properties between them. In the second step, a material database is generated by transforming the effective stiffness and damping tensors over a discrete number of fiber distribution functions (ODF). An assignment of the properties is based on the minimum error deviation between the discrete orientation tensors of the material database and the orientation tensor of each element of the composite component. The orientation tensors of the component are taken from corresponding injection molding simulations, Section 2.2. Next, the geometry of the plastic component is discretized into an FE model. Here, the number of elements of the discretized model controls the provided fiber orientation tensor significantly from the injection molding simulation [1]. As a result, to investigate the influence of the fiber orientation on the prediction quality of the structural dynamics simulation, the surface mesh remains identical, and the number of elements over the wall thickness is iteratively increased. In the last step, the boundary conditions of excitation and bearing are modeled. To avoid additional simulative impacts, a two-dimensional node-force excitation is applied. Furthermore, an unbound numerical bearing is assumed. Figure 4 provides an overview of the method for the experimental and simulative characterization of the structural dynamics of the short-fiber reinforced plastic components.

**Figure 4.** Method for the experimental and simulative characterization of the structural dynamics of short-fiber reinforced thermoplastic components.

Figure 5a shows the experimental setup for investigating the structural dynamics of the short-fiber reinforced, thermoplastic engine bracket by modal analyses in freesupported bearing. The engine bracket is suspended with 4 aramid ropes. Preliminary investigations confirmed that this leads to an optimum between stiffness and damping of the bearing method and thus minimizes a reaction on the engine bracket [44–46]. An electrodynamic shaker from BRÜEL & KJAER type 4809 [47] is used for excitation. This is equipped with a BRÜEL & KJAER type 8001 impedance measuring head and stinger to record the excitation force [48]. The structural dynamics response is recorded with a 3D laser vibrometer type PSV-500 from Polytec [49]. The measurement mesh consists of 145 scanning points with a focus on the qualitative evaluation of the mode shapes, as well as quantitative evaluation of the force reaction function (FRF), as shown in Figure 5b. In the

present case, the experimental FRF is defined as the ratio between the spatially averaged acceleration spectrum of all scanning points on the front surface of the engine bracket and the excitation force. The experimental modal analysis is performed at constant conditioning of 23 ◦C and 0% relative humidity. To simulate the structural dynamics of the engine bracket, the surface geometry is discretized into 137,128 elements in the first step, as shown in Figure 5c. Next, the number of FE elements is varied via the component thickness to investigate the resulting influence. The simulative FRF is defined as the ratio between the spatially averaged acceleration spectrum of all nodes on the front surface of the engine bracket and the excitation force, as shown in Figure 5d. Figure 6 shows the proposed approach for evaluating and comparing the experimental and structural dynamics of the short-fiber reinforced thermoplastic engine bracket.

**Figure 5.** (**a**) Setup for the experimental investigation of the structural dynamics of the engine bracket. (**b**) Experimental scanning point mesh. (**c**) Simulative surface mesh of the engine bracket geometry. (**d**) Evaluation nodes for simulative FRF.

**Figure 6.** Proposed approach for the experimental and simulative characterization and comparison of the structural dynamics of the short-fiber reinforced thermoplastic engine bracket.

### **3. Results**

### *3.1. Simulative and Experimental Fiber Orientation Investigations*

Figure 7 shows the experimental fiber orientation tensors (CT) in comparison to the simulative fiber orientation tensors with CADMOULD® (CM) and MOLDFLOW® (MF) over the normalized position for the initial coefficient tensors *CNC*,1 0,*CM* and *<sup>C</sup>OC*,1 0,*MF* for the

material PA66GF50 and for an evaluation on plate level. For quantitative comparison of the calculated fiber orientation with CM and MF to the CT results, the relative deviation of the averaged fiber orientation tensors is shown in Table 1.

**Figure 7.** Experimental (CT) and simulative (CM, MF) components of the fiber orientation tensor **A***ij* plotted over normalized position for PA66GF50 at plate level (P1–P5).

**Table 1.** Relative deviation of the averaged orientation tensors from experiment *Aij*,*CT* and simulation *Aij*,*CM*, *Aij*,*MF* of PA66GF50 on plate level (P1–P5).


Figure 7 shows that the boundary layer equals the main part of the fiber orientation tensor of the μCT results. Furthermore, the majority of the μCT results have a pronounced middle layer. P1 and P4 even show a change in the main fiber orientation *A*11,*CT*, whereby the fibers in the middle layer are oriented transverse to the main flow direction. The comparison between *A*11,*CT*, *A*22,*CT* and *A*11,*CM*, *A*22,*CM* shows, in general, a good prediction of the fiber orientation by the simulation with CM approach, as shown in Table 1. The comparison between *A*11,*CT*, *A*22,*CT* and *A*11,*MF*, *A*22,*MF* shows a good prediction, too. However, both the CM and MF results show, in general, that the development of the middle layer is not sufficiently represented. On the one hand, according to the evaluation method, elements are missing in the middle of the FE model. This results in a loss of information in the process–structure coupling since only a weighted average of the calculated fiber orientation tensor from the injection molding simulation is transferred [1]. On the other hand, the Advani—Tucker Equation with the rotary diffusion approach only allows a limited representation of the fiber orientation in the middle layer [8,34,38]. Furthermore, Figure 7 shows that the component *A*33,*MF* has a high deviation compared to *A*33,*CT*. This difference can be explained by the method used to calculate the fiber orientation of the highly filled material coupled with the 3D modeling. As a result, the fiber orientation tensor is not only influenced by the phenomenological coefficient tensors, but also by the three-dimensional flow field and thus by the hydrodynamic part of the Advani–Tucker Equation, see [8]. Figure 8 shows the experimental (CT) and simulative (CM, MF) fiber orientation tensors for the coefficient tensors *CNC*,1 0,*CM*, *<sup>C</sup>OC*,1 0,*MF* for the material PA66GF50 at the component level. Table 2 shows the corresponding deviations of the averaged fiber orientation tensors.

**Figure 8.** Experimental (CT) and simulative (CM, MF) components of the fiber orientation tensor **A***ij* plotted over normalized position for PA66GF50 at component level (B1–B2).



Figure 8 shows the CT results of B1, where a pronounced boundary layer with a change in the main fiber orientation occurs. In comparison, no middle layer can be identified. Thus, a significantly different trend of the fiber orientation compared to the plate level can be determined at specimen position B1. Figure 8 also shows that the fiber orientation tensors of the CT result B2 are comparable to those in the plate level. Following this allows a robust check of the injection molding simulations according to the calculated

fiber orientation. Figure 8 shows differences in the prediction quality of the CM and MF results compared to the CT results of specimen B1. Thereby, the CM results show a good prediction of the boundary layer, while the prediction of the MF results is not sufficient. An inverse prediction quality can be identified for the middle layer. On the other hand, the CM and MF results show a good prediction of the CT results at specimen B2. Thus, it can be assumed that good predictions of the fiber orientation can be achieved at the component level for similar flow conditions compared to the plate geometry. The quantitative comparison between *A*11,*CT*, *A*22,*CT* and *A*11,*CM*, *A*22,*CM* and *A*11,*MF*, *A*22,*MF* shows a sufficient prediction of the fiber orientation, Table 2. Furthermore, it can be shown on the component level that *A*33,*MF* has a high deviation compared to *A*33,*CT*. Finally, it should be mentioned that both at the plate and component level, the prediction quality of the CM and MF results could not be further improved by a parameter study of *CNC*,1 *i*,*CM* and *COC*,1 *<sup>i</sup>*,*MF*.

Next, the robustness of the injection molding simulations is checked in the presence of a material variation. Figure 9 shows the CT, CM and MF fiber orientation tensors for the coefficient tensors *CNC*,2 0,*CM* and *<sup>C</sup>OC*,2 0,*MF* for the material PPAGF50 at plate level. The CT results show a small boundary and significant middle layer. Thereby, the CM and MF results show a significant deviation compared to the CT results. It is assumed that this deviation is due to the loss of information of the evaluation method and the process–structure coupling. Further, the coefficient tensors *CNC*,2 0,*CM* and *<sup>C</sup>OC*,2 0,*MF* are insufficient. The evaluation of the relative deviation of *A*11,*CT*, *A*22,*CT*, *A*33,*CT* to *A*11,*CM*, *A*22,*CM*, *A*33,*CM* and *A*11,*MF*, *A*22,*MF*, *A*33,*MF* also show an insufficient prediction of the fiber orientation, as shown in Table 3.

**Figure 9.** Experimental (CT) and simulative (CM, MF) components of the fiber orientation tensor **A***ij* plotted over normalized position for PPAGF50 at plate level (P1–P5).

**Table 3.** Relative deviation of the averaged orientation tensors from experiment *Aij*,*CT* and simulation *Aij*,*CM*, *Aij*,*MF* of PPAGF50 on plate level (P1–P5).


Figure 10 shows the CT, CM and MF fiber orientation tensors for the coefficient tensors *CNC*,2 0,*CM* and *<sup>C</sup>OC*,2 0,*MF* for the material PPAGF50 at the component level. The CT results of B1 show a significant boundary layer with a change of the main fiber orientation and a small developed middle layer. The CT results of B2 show a pronounced boundary and middle layer comparable to the PA66GF50. The comparison between CT and CM results shows with the difference in the relative deviation of Δ *Aij*,*CT* <sup>→</sup> *Aij*,*CM* = 15.0% an improved prediction of the fiber orientation in comparison to the investigations at the plate level, Tables 3 and 4. In this context, with a difference of Δ *Aij*,*CT* <sup>→</sup> *Aij*,*MF* <sup>=</sup> <sup>−</sup>25.5% the comparison between CT and MF results shows an increased deviation compared to the plate level. Thereby, it can be shown that the boundary and middle layer of specimen B1 are sufficiently predicted by the CM results but not by the MF results, as shown in Figure 10. However, the CM and MF results show a comparable prediction of the fiber orientation of specimen B2. The quantitative comparison between *A*11,*CT*, *A*22,*CT*, *A*33,*CT* and *A*11,*CM*, *A*22,*CM*, *A*33,*CM* shows a sufficient prediction of the fiber orientation, Table 4. In this context, the comparison between *A*11,*CT*, *A*22,*CT*, *A*33,*CT* and *A*11,*MF*, *A*22,*MF*, *A*33,*MF* shows an increased deviation due to the insufficient prediction quality of *A*33,*MF*.

**Figure 10.** Experimental (CT) and simulative (CM, MF) components of the fiber orientation tensor **A***ij* plotted over normalized position for PPAGF50 at component level (B1–B2).

A parametrical study of the coefficient tensor *CNC*,2 *<sup>i</sup>*,*CM* shows, that the prediction quality of the fiber orientation can be improved at the plate level. The phenomenological coefficients of the optimized coefficient tensor *CNC*,2 *<sup>i</sup>*,*CM* are redefined as *α<sup>K</sup>* = 0.11, *α<sup>R</sup>* = 0.9, *β<sup>R</sup>* = 0.05 Figure 11 shows the CT and CM fiber orientation tensors for the optimized coefficient tensor *CNC*,2 *<sup>i</sup>*,*CM* of the PPAGF50 material at plate level. Figure 11 shows that the prediction quality of the CM results of the boundary and middle layer is sufficiently improved. Furthermore, the comparison between *A*11,*CT*, *A*22,*CT* to *A*11,*CM*, *A*22,*CM* confirms that the prediction quality of the fiber orientation can be optimized, as shown in Table 5.


**Table 4.** Relative deviation provided in [23,50]. Figure averaged orientation tensors from experiment *Aij*,*CT* and simulation *Aij*,*CM*, *Aij*,*MF* of PPAGF50 on component level (B1–B2).

**Figure 11.** Experimental (CT) and simulative (CM) components of the fiber orientation tensor **A***ij* plotted over normalized position for PPAGF50 at plate level (P1–P5) in optimized condition *CNC*,2 *<sup>i</sup>*,*CM*.

**Table 5.** Relative deviation of the averaged orientation tensors from experiment *Aij*,*CT* and simulation *Aij*,*CM* of PPAGF50 on plate level (P1–P5) in optimized condition *<sup>C</sup>NC*,2 *<sup>i</sup>*,*CM*.


Next, the robustness of the CM simulation for PPAGF50 with optimized coefficient tensor *CNC*,2 *<sup>i</sup>*,*CM* is tested on the component level, as shown in Figure 12. Specimen B1 shows higher deviations between CT and CM results in the expression of the middle layer. Furthermore, the comparison between the CT and CM results at specimen B2 shows significantly higher deviations. This is due to the significant adjustment of the alignment factor *α<sup>K</sup>* which leads to a higher deviation between CT and CM results, as shown in Table 6. In summary, it can be shown that an optimized coefficient tensor at the plate level does not necessarily constitute an improved prediction quality of the fiber orientation at the component level. Thus, limits of the calculation method of the fiber orientation and the modeling in injection molding simulation are reached. The default values of the initial coefficient tensors *CNC*,*<sup>S</sup>* 0,*CM* and *<sup>C</sup>OC*,*<sup>S</sup>* 0,*MF* are therefore considered as balanced optimum for

the calculation of the fiber orientation at plate and component level with an averaged relative deviation of *Aij*,*CT* <sup>→</sup> *Aij*,*CMPA*<sup>66</sup> <sup>=</sup> 19.2%, *Aij*,*CT* <sup>→</sup> *Aij*,*MFPA*<sup>66</sup> <sup>=</sup> 67.1% and *Aij*,*CT* <sup>→</sup> *Aij*,*CMPPA* <sup>=</sup> 24.1%, *Aij*,*CT* <sup>→</sup> *Aij*,*MFPPA* <sup>=</sup> 116.6%. This is an expected result since the suppliers of the material databases provide an extensive optimization of the phenomenological fiber interaction coefficient tensors, see [8,43]. If a significantly increased prediction quality of the simulative fiber orientation is necessary, case-specific optimizations of the component in injection molding simulations are required.

**Figure 12.** Experimental (CT) and simulative (CM) components of the fiber orientation tensor **A***ij* over normalized position for PPAGF50 at component level (B1–B2) in optimized condition *CNC*,2 *<sup>i</sup>*,*CM*.

**Table 6.** Relative deviation of the averaged orientation tensors from experiment *Aij*,*CT* and simulation *Aij*,*CM* of PPAGF50 on component level (B1–B2) in optimized condition *<sup>C</sup>NC*,2 *<sup>i</sup>*,*CM*.


Finally, the influence of the process–structure coupling is investigated. Thereby, the Kullback–Leibner divergence between the experimental and simulative orientation distribution is evaluated. The Kullback–Leiber divergence represents the similarity between probability distributions. Thus, low divergence equals a high similarity and vice versa. For simulative reconstruction of the orientation distribution, the maximum entropy method (MEM) is used. Further information on the Kullback–Leiber divergence and the MEM is provided in [23,50]. Figure 13 shows Kullback–Leibner divergences D*KL* between CT, CM and MF fiber orientation distribution for PA66GF50 and PPAGF50 at the component level (B1–B2) plotted over the layers of the FE model. Figure 13 shows that the CM results at 10 layers have the smallest deviation compared to the CT results. In comparison, the deviation of the MF results at 10 layers is significantly increased due to a high deviation between *A*33,*CT* and *A*33,*CM*. Furthermore, this influence is more pronounced for the PPAGF50. By reducing the number of layers in the process–structure coupling, the transferred fiber orientation content is reduced. This leads to an increased deviation between experimental and simulated fiber orientation distribution. In general, from 10 layers to 1 layer, there is an increased deviation due to a loss of fiber orientation content originating from the injection molding simulation. This is due to the process–structure coupling. By reducing the number of FE layers, a weighted average of the simulative fiber orientation tensors from the injection molding simulation is transferred. This averaging method significantly influences the deviation regardless of the material, the specimen position and the type of injection molding simulation.

**Figure 13.** Kullback–Leibner divergences D*KL* of the fiber orientation tensors of the simulation *Aij*,*CM*, *Aij*,*MF* compared to the experiment *Aij*,*CT* of the PA66GF50 (black) and PPAGF50 (grey) on component level (B1–B2) and plotted over the number of layers of the FE-model.

### *3.2. Experimental and Simulative Structural Dynamics Investigations*

Figure 14 shows the experimental FRF of the PA66GF50 and PPAGF50 engine bracket. The first relevant mode shape in the frequency range between 400–500 Hz corresponds to a global torsional mode, as shown in Table 7. Due to the higher stiffness of the PPAGF50, the resonance frequency of the torsional mode is shifted by 6.0% towards higher frequencies compared to the PA66GF50. Furthermore, the amplitude is increased by 5.7% due to the lower damping of the PPAGF50. A comparison via the Modal Assurance Criterion (MAC) shows comparability of 90.2% between the experiments. Above 1900 Hz, the second relevant mode shape is a global bending mode of the engine bracket, as shown in Table 7. Comparable to the first mode shape, the higher stiffness of the PPAGF50 leads to a shift of the resonance frequency by 5.5% towards higher frequencies. However, compared to the previous trend, the amplitude of the PPAGF50 is 61.3% lower than the PA66GF50. It can be assumed that the frequency-dependent structural damping is more pronounced than the material damping. This assumption is reliable by the fact that the MAC analysis shows comparability of 82.6% only. A similar comparison can be shown for the third relevant mode shape, a local surface vibration, as shown in Table 7.

Comparable to the second mode shape, the PPAGF50 shows a stiffer behavior and higher structural damping. Thus, the resonance frequency is shifted by 5.1% towards higher frequencies and the amplitude is 40.3% lower compared to the PA66GF50. Furthermore, the MAC analysis shows comparability of only 69.9%. This results in significantly different mode shapes. Pronounced mixed, local mode shapes characterize frequency ranges above 2300 Hz. A definite evaluation and interpretation of the resonance peaks are no longer ensured. Following, the three identified resonances according to Table 7 are used for further comparison with the simulation. Considering the deviations between the experimental results, it is assumed that not only a superposition of material and component properties occurs, but also the overall structure influences the results, see Section 2.3, Figure 5 and [21,46,51]. Thus, the experimental structural dynamics investigations represent a trend and not an exact behavior.

**Figure 14.** Experimental FRF of the modal analysis of the PA66GF50 and PPAGF50 engine bracket.

**Table 7.** Experimental mode shapes with resonance frequency, amplitude and MAC value of the modal analysis of the PA66GF50 and PPAGF50 engine bracket.

Concerning the proposed approach of Section 2.3 and the results of Section 3.1, more than six elements over the thickness do not lead to an increase in the fiber information content in the FE model, as shown in Figure 13. Furthermore, the engine bracket in the simulation is considered with a discretization between 1 to 6 element layers over the component thickness in this results section. The evaluation of the FRF of the simulations in the relevant frequency range allows the identification of the mode shapes, as shown in Table 8. Comparable to the experiment, the relevant global torsion and bending mode of the engine bracket can be identified in the simulation. Both in the simulation with fiber orientation from CM and MF results, a high MAC value can be verified for all layers. Furthermore, local surface mode shape can also be identified. However, with an average MAC value of 74.0%, the simulation shows noticeable differences compared to the experiment. Nevertheless, a comparison of all performed simulations with numerical fiber orientation from CM and MF results is sufficient. The evaluation of the simulations of the PPAGF50 engine bracket provides comparable mode shapes and MAC values for all layers. Corresponding to the experiments, this is an expected result, as the stiffness, damping and viscoelasticity of the PPAGF50 cause a shift in the resonance frequency and amplitude. Thus, the mode shape is only slightly affected. Furthermore, the mode shapes and MAC values of the simulation of the PPAGF50 engine bracket are not displayed.

**Table 8.** Mode shapes of the structural dynamics simulation and MAC values in comparison to the experiment of the PA66GF50 engine bracket.

Figure 15 shows the experimental FRF of the PA66GF50 engine bracket compared to the structural dynamics simulation with fiber orientation tensor of the CM and MF results. Figure 15 shows for the global mode shapes (torsion, bending) the highest deviation of the FRF from the simulation with 1 layer compared to the experiment. This applies to the simulations with CM and MF fiber orientation, as shown in Table 9. This shows a good correlation to the deviations with a process–structure coupling of the fiber orientation, as shown in Figure 13. This trend is continued for a simulation with 2 and 3 layers. Thereby, increasing the number of layers generates an increased fiber orientation content and thus reduces the deviation between experimental and simulated FRF. Contrary, increasing the number of layers to 4 and 5 using the CM results shows a slight increase in the deviation between experimental and numerical FRF. It can be assumed that increases in the deviations are due to the weighted averaging of the orientation tensor in the process–structure coupling, see [1]. Further, this influences the material modeling and the structural dynamics simulation, as shown in Figure 6. However, by increasing the number of layers to 6, a reduction in the deviation between experimental and numerical FRF can be shown, using both the CM and MF results.

**Table 9.** Relative deviation of frequency and amplitude between simulations and experiment of the PA66GF50 engine bracket for identified mode shapes with fiber orientation from CM and MF results.


**Figure 15.** Experimental (EXP) and simulative (CM, MF) FRF of the structural dynamics investigations of the PA66GF50 engine bracket. Simulations are plotted as a function of used FE-layer (L).

A comparative evaluation of the local mode shape (surface) shows an inverse trend of the deviation when increasing the number of layers from 1 to 6, as shown in Figure 15 and Table 9. By increasing the number of layers, the boundary and middle layer of the fiber orientation are more pronounced in the process–structure coupling. Weighted over the component thickness, this leads to a reduction of the stiffness in the material model. As a result, the entire FRF shifts towards lower frequencies. The influence of the material modeling applies to both the local and the global structural dynamics phenomena and shows a good correlation to existing material modeling studies, see [4,7,15,21].

Figure 16 shows the experimental FRF of the PPAGF50 engine bracket compared to the simulation with fiber orientation tensor of CM and MF. Figure 16 shows that the simulation with CM results in the global torsion mode that increasing the layers from 1 to 6 leads to a shift in the resonance frequency. Thereby, increasing the layers improves the mapping of the fiber orientation, which positively affects the stiffness and damping of the material model. Furthermore, it can be shown that the simulations with 6 layers allow the best prediction of the structural dynamics, as shown in Table 10. However, this is not the case for the simulations with MF results. Thereby, the increase of layers leads to an increase in the deviation between experiment and simulation. Again, the fiber orientation of the MF results affects the stiffness and damping of the material model. The evaluation of the global bending mode and local surface mode shows a comparable trend for the CM and MF results. Increasing the layers reduces the deviation between simulated and experimental FRF, as shown in Figure 16 and Table 10. However, side resonances are formed in the simulation, which cannot be identified in the experiment. The influence of a high local deviation between experimental and simulated fiber orientation becomes significant and leads to different local stiffness and damping in the material model.

**Figure 16.** Experimental (EXP) and simulative (CM, MF) FRF of the structural dynamics investigations of the PPAGF50 engine bracket. Simulations are plotted as a function of used FE-layer (L).



Concerning existing works of literature, the results of this research significantly show the importance of reproducibly describing the boundary conditions of the injection molding simulation, the process–structure coupling and the material modeling. With the used framework of the integrative simulation, relative averaged deviations of 2.8% in the frequency

and 38.0% in the amplitude of the FRF can be proven in the structural dynamics simulation, as shown in Tables 9 and 10. Furthermore, it can be shown that high deviations between experimental and simulative fiber orientation tensors can lead to a sufficient prediction of the structural dynamics simulation. This can be verified for global and local structural dynamics phenomena.

### **4. Conclusions**

This contribution investigates the influence of fiber orientation from injection molding simulation on the structural dynamics simulation of the short-fiber reinforced thermoplastic components, e.g., engine brackets made of PA66GF50 and PPAGF50. The Advani–Tucker Equation with rotary diffusion approach and phenomenological fiber interaction coefficient tensor was used to calculate the numerical fiber orientation. The comparison of the experimental and simulative fiber orientations shows a sufficient prediction of the PA66GF50. Optimizations of the interaction coefficients on the plate and component level constitute no further improvements. On the other hand, optimizations of the interaction coefficients of the PPAGF50 on the plate level lead to an improvement and, on the component level, to an insufficient prediction of the simulative fiber orientation compared to the experiments. Simulations of the structural dynamics of the PA66GF50 and PPAGF50 engine bracket show that increasing the layers over the component thickness of the model leads to an improved prediction quality compared to the experiment. This can be shown for global and local structural dynamics phenomena. However, from a simulative point of view, the difference in the FRF with different layers is significantly smaller than the loss of information of the fiber orientation due to a process–structure coupling. In summary, the following conclusions can be derived:


**Author Contributions:** Methodology, A.K. and M.S.; software, A.K.; investigation, A.K.; validation, A.K. and M.S.; writing—original draft preparation, A.K.; writing—review and editing, M.S.; visualization, A.K.; supervision, M.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data is not publicly available.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**

