**1. Introduction**

Due to the ever-decreasing development times and the steady progress of digitization, numerical methods for predicting unknown target and design variables are gaining in importance. This also applies to the processing of plastics and the resulting part properties in general. Especially in the case of fiber-reinforced plastics, the strongly varying part properties directly depend on the fiber microstructure and must be taken into account during the design process [1–4]. In the following, the fiber microstructure is used as a generic term for fiber orientation, fiber length, and fiber concentration.

A reliable prediction of the process-induced fiber microstructure, especially for long fiber reinforced thermoplastics (LFT), still poses a considerable challenge today [5–7]. In contrast to short fiber reinforced thermoplastics (SFT), not only the fiber orientation must be taken into account in the prediction of the mechanical properties of LFTs, but effects such as fiber breakage and the migration of fibers during processing also play a significant role. In order to predict these properties correctly, suitable calculation models are required on the one hand, and correspondingly well-calibrated model parameters on the other hand.

However, the reliable prediction of fiber orientation in general, but especially for LFTs, is still the subject of current research and requires further improvement in the achievable prediction quality. Current commercial software environments that calculate the injection molding process do not allow a purposeful and fast optimization of the resulting fiber orientation distribution. Only a manual adaptation of model parameters with a subsequent recalculation of the entire flow field is possible. Due to the lack of automation and the necessary recalculations, the calibration of the model parameters becomes remarkably ine fficient and time-consuming. Therefore, a new calibration method was developed to automatically adjust the resulting fiber orientation distribution.

#### *1.1. State of the Art*

During the processing of discontinuous fiber reinforced plastics, a specific fiber microstructure (fiber orientation, fiber length and fiber concentration) is formed within the part cavity that depends on the material, the geometry of the molded part, and the process settings. This microstructure determines the parts main properties. With suitable accuracy of the predicted fiber microstructure the resulting anisotropic mechanical properties can be reliably calculated [8–12]. Thus, the part's behavior can also be determined in structural simulations.

#### 1.1.1. Fiber Microstructure

The resulting part properties of fiber reinforced materials are mainly determined by process induced fiber microstructure and its properties (Figure 1), such as:


**Figure 1.** Influence of fiber orientation, length, and concentration on the resulting mechanical properties (schematically) according to [16,26,29].

The increase in the respective fiber microstructure properties (orientation, length, and concentration) also leads to an increase in mechanical properties of the composite material, such as sti ffness, strength, and impact strength [23]. Only at high fiber concentrations can a decrease in strength and impact strength be observed. With regard to an entire fiber reinforced part, the local distribution of these properties along the flow path as well as along the cavity height (part thickness) must also be taken into account [31]. Even if in some cases a classification cannot be clearly distinguished, fiber reinforced materials are categorized into three classes according to their initial or resulting fiber length [32]: short fiber reinforced plastics (0.1–1 mm), long fiber reinforced plastics (1–50 mm), and continuous fiber reinforced plastics (>50 mm).

#### 1.1.2. Process Induced Microstructure

The properties of the fiber microstructure (orientation, length, and concentration), especially in the case of discontinuously reinforced plastics such as SFT or LFT, can usually only be influenced to a very limited extent by its constituent properties and compound constitution (initial fiber orientation, initial fiber length and initial fiber concentration). However, these properties are rather a result of the processing [13,26,33–37]. One of the most widely used processing methods is injection molding, in which even complex parts can be produced in large quantities [32,37].

The materials used are subjected to special stresses during the injection molding process. This can be seen in the orientation of the fibers as well as in the migration and fracture phenomena. During the filling of a part with molten plastic by injection molding, a flow field is formed inside the cavity that depends on the material, the part geometry, and the process settings. This varies not only along the flow direction, for example through cooling effects, but also along the part's thickness or the cavity height, as shown in Figure 2. Due to this velocity profile and the resulting shear stresses, characteristic layers are formed along the part thickness (Figure 3) in which the fibers


**Figure 2.** Velocity profile and the resulting fiber orientation distribution according to [51].

**Figure 3.** Locally varying fiber microstructure along the normalized part thickness (schematically).

These individual microstructure properties also affect each other. Due to the interaction of the fibers, the fiber orientation depends on the fiber length and concentration [5], which also applies to interchanged dependencies. The fiber length depends on the fiber orientation and concentration while the distribution of fibers depends inversely on the fiber orientation and length [36]. Due to wall adhesion and mold cooling, high shear rates are present in the skin or shear layers, and lower shear rates are found in the core layer in the middle of the cavity. A characteristic profile for each property is created inside a part, which can be subdivided into different layers [38,41,52] as shown in Figure 3.

According to the crystalline structure [34,53] or the fiber orientation distribution [46,54], the flow profile can be divided into at least a core layer and two shear or skin layers. This three-layer model is shown schematically in Figure 2 (right), and the individual microstructure properties are shown in Figure 3. The specific characteristic ratio of the layers is determined by the quotient of the individual layer thicknesses and is generally referred to as the skin core layer ratio.

The influence of the velocity field on the resulting fiber orientation distribution is a complex interrelation. The flow field changes in response to many different influencing factors. Important parameters are, for example, the injection speed and the processing temperature, as well as the material's viscosity [14]. Figure 4 shows the correlation of the shear rate and injection speed [55] along the cavity height.

**Figure 4.** Influence of the injection speed on the velocity profile according to [55].

At high injection speeds (Figure 4b), the maximum shear rate shifts towards the mold wall, and the core layer width increases [39]. Whereas, for lower injection speeds (Figure 4a) the maximum of the shear rate profile shifts towards the cavity core, the range of the high shear rates generally increases, and the core layer width decreases. Depending on the shear and extensional flow inside the cavity, individual layers are formed with their associated fiber orientation distribution [19].

Furthermore, the viscosity of the melt changes with the type of matrix polymer and the reinforcing filler material (material, geometry, concentration, length, diameter, orientation, adhesion etc.) [42] as shown in Figure 5a. For example, a higher fiber concentration results in a higher viscosity, which leads to a different velocity and shear rate profile, as shown in Figure 5b,c [32,56]. In areas of low shear rates in particular, the viscosity of filled polymers differs strongly from that of unfilled polymers. Higher fiber concentrations mostly result in smaller skin layers due to their shear thinning behavior [19,39,57,58]. For this reason, different viscosity models are required to describe specific material behavior at low shear rates for reinforced or filled polymers (Figure 5b). For injection molding simulation, a frequently used model was developed by Herschel and Bulkley. A detailed description of the Herschel–Bulkley model is given in Section 1.1.7.

**Figure 5.** Influence of the viscosity to the velocity profile according to [59].

#### 1.1.3. Fiber Orientation Distribution of Short and Long Fiber Reinforced Thermoplastics

As mentioned before, the properties of the reinforcement itself and the fiber microstructure properties also a ffect each other. In this context, the resulting fiber length inside a part plays a crucial role in the resulting fiber orientation [1]. The fiber length in particular causes a varying layer ratio, which can be seen in the fiber orientation distribution of short and long fiber reinforced plastic parts [39,44]. Even if the fiber breakage and the resulting fiber length are also determined by its orientation inside the flow, the resulting fiber orientation is more dominated by the fiber length. As Figure 6 schematically demonstrates, the degree of orientation in the skin layer decreases and the core layer width increases with increased fiber length [39,44,47,54]. Furthermore, in some cases, the degree of fiber orientation perpendicular to the flow direction ( *A*22) in the core layer can also exceed the degree of fiber orientation parallel to the flow direction ( *A*11), as shown in Figure 6b, for a long fiber reinforced plastic. In that special case, the mechanical properties perpendicular to the flow can also exceed the mechanical properties in the direction of flow.

**Figure 6.** Characteristic fiber orientation distributions along the parts thickness in discontinuous fiber reinforced plastics.

#### 1.1.4. Tensorial Description of the Fiber Orientation

The fiber orientation is usually described as a 2nd order orientation tensor *A* = *aij* according to Advani and Tucker [60]. The tensor describes the orientation of all fibers in a defined volume or element where <sup>Ψ</sup>(*p*) is a probability density function, which is a statistical description of orientation states. The orientation tensor *A* is calculated as an integral over all possible fiber orientations from the distribution function of the fibers <sup>Ψ</sup>(*p*) and the dyadic product of the unit vectors *p*. The first entry *a*11 describes the orientation of the fibers in the flow direction and is represented in Figure 7a. The eigenvectors Λ*i* of the tensor indicate the principal direction of fiber orientation, while the eigenvalues *ei* indicate the degree of fiber orientation, i.e., the statistical proportion of fibers in the respective principal direction.

**Figure 7.** Tensorial Description of the Fiber Orientation.

The calculation of fiber orientation is based on a previously computed flow field. The velocity gradient tensor *L* and the shear rate .γ are required as input variables. The result is a tensor *A* that describes the orientation of the fibers in the melt or flow channel in three dimensions [60].

$$A = a\_{ij} = \begin{bmatrix} a\_{11} & a\_{11} & a\_{11} \\ a\_{11} & a\_{22} & a\_{11} \\ a\_{11} & a\_{11} & a\_{33} \end{bmatrix} \tag{1}$$

The diagonal entries *a*11, *a*22, and *a*33 describe the orientation in the three spatial directions.

#### 1.1.5. Experimental Determination of the Fiber Orientation

Several methods have been developed to analyze the fiber microstructure, such as ultrasound methods [61,62], optical methods like grinding pattern analysis [35,40,60,63–70], electron microscopy [14], Radiography [45,71,72] and X-ray computed tomography (CT) [73–78]. However, all methods are usually associated with high effort and are therefore used infrequently. To determine the spatial fiber orientation by means of cross section analysis, each fiber ellipsis, as shown in Figure 8a, has to be measured by its major and minor axis according to Figure 8b. With the trigonometric functions depicted in Figure 8c, the fiber orientation tensor components can be calculated for each fiber [40,60]. In addition, the reduced cutting probabilities of the inclined fibers can be considered by applying a weighting function according to [40,79].

**Figure 8.** Cross section analysis to experimentally determine fiber orientation.

#### 1.1.6. Injection Molding Simulation

For the simulation of the unsteady, non-Newtonian and non-isothermal injection molding process, the three-dimensional basic Equations for fluid mechanics are usually solved using the finite volume method (FVM). The basic equations for conservation of mass (2), conservation of linear momentum (3), and energy (4), are solved by FVM numerically [16,80].

conservation of mass:

$$\frac{\partial \rho}{\partial t} + \nabla \cdot \rho \,\overrightarrow{\boldsymbol{v}} = 0 \tag{2}$$

conservation of momentum:

$$
\rho \left( \frac{D \ \overrightarrow{v}}{Dt} \right) = \rho \,\overrightarrow{\text{g}} + \nabla \cdot \underline{\text{r}} - \nabla \cdot \text{p} \tag{3}
$$

conservation of energy:

$$\frac{\partial \rho \boldsymbol{T}}{\partial t} = -\nabla \cdot (\rho \,\overrightarrow{\boldsymbol{v}} \,\boldsymbol{T}) + \nabla \cdot \left(\frac{\Lambda}{c\_p} \cdot \nabla \boldsymbol{T}\right) + \underline{\boldsymbol{\tau}} \,:\, \boldsymbol{\nabla} \,\overrightarrow{\boldsymbol{v}}\tag{4}$$

Based on these three equations, flow velocities ∇ and their gradients <sup>∇</sup><sup>→</sup>*v* are calculated, which represent the major determining variables for the subsequent fiber orientation calculation. The velocity gradient tensor <sup>∇</sup><sup>→</sup>*v* describes the change of velocity in the three spatial directions of a

stationary flow and can be split into a symmetrical *D* and a rotational *W* component (Equation (5)). From the deformation velocity tensor *D*, the shear rate . γ can be derived (Equation (6)). However, for a detailed description of the FVM for injection molding simulation, refer to [80].

$$
\nabla \overrightarrow{\boldsymbol{v}} = \boldsymbol{L} = \frac{1}{2} (\nabla \overrightarrow{\boldsymbol{v}} + (\nabla \overrightarrow{\boldsymbol{v}})^T) + \frac{1}{2} (\nabla \overrightarrow{\boldsymbol{v}} - (\nabla \overrightarrow{\boldsymbol{v}})^T) = \underline{\boldsymbol{D}} + \underline{\boldsymbol{W}} \tag{5}
$$

$$\dot{\boldsymbol{\gamma}} = \begin{vmatrix} \underline{\boldsymbol{D}} \end{vmatrix} = \sqrt{\frac{1}{2} \underline{\boldsymbol{D}} : \underline{\boldsymbol{D}}^T} \tag{6}$$

## 1.1.7. Viscosity Models

To describe the flow behavior of plastic melts [81], several di fferent viscosity models were developed. The Herschel–Bulkley model is frequently used to accurately represent the viscosity of long fiber reinforced thermoplastics especially in the range of low shear rates. In this model, on the one hand, the yield stress viscosity at low shear rates is considered, which can be determined experimentally by rotational rheometry [82,83]. On the other hand, the range of high shear rates is precisely covered where so-called shear thinning [84] is present. This range can experimentally be determined by means of capillary rheometry. The model is composed of the Herschel–Bulkley yield stress (first term) and the Cross-WLF viscosity (second term), as shown in Equation (7) with power-law index *n*, temperature *T*, pressure *p*, and zero shear viscosity η0.

$$\eta\left(\dot{\mathbf{y}},T,P\right) = \frac{\tau\_y}{\dot{\mathbf{y}}} + \frac{\eta\_0(T,p)}{1 + \left(\frac{\eta\_0}{\tau^\ddagger}\dot{\mathbf{y}}\right)^{1-n}}\tag{7}$$

The model consists of seven parameters *N*, τ<sup>∗</sup> , *D*1, *D*2, *D*3, *A*1, and *A*2 in order to approximate the shear rate distribution of the materials viscosity.

#### 1.1.8. Fiber Orientation Models

To calculate fiber orientation distributions, the tensorial description developed by Advani and Tucker [60] is usually used. The corresponding calculation models are mostly based on the extension of the Je ffery model [85] made by Folgar and Tucker [86]. Current commercial injection molding solvers also use combined model extensions, such as the RPR (retarding principal rate) and ARD or iARD (improved anisotropy rotary di ffusion) models. The iARD-RPR model is a combination of Je ffery's hydrodynamic model (HD), the extension of the IRD model (isotropic rotary di ffusion) from Folgar–Tucker to the iARD model and the RPR model (Equation (12)), according to Tseng et al. [87–89]. With *Ci* (fiber–fiber interaction coe fficient), *Cm* (fiber–matrix interaction coe fficient), and α (slow-down parameter), a changing of the fiber orientation based on the existing flow field can be calculated. The influencing factors of the calculated fiber orientation are in this case the gradient tensor *L* and the shear rate . γ of the flow field.

$$A = A\_{\rm HD} + A\_{\rm iARD}(\mathbb{C}\_i, \mathbb{C}\_m) + A\_{\rm RPR}(\alpha) \tag{8}$$

This leaves three model parameters *Ci*, *Cm*, and α, which can be used to adjust or calibrate the fiber orientation distribution. The first part is defined as iARD (Equation (9)) [87].

.

.

.

$$\dot{\vec{A}}^{\dot{A}R\dot{R}D} = \dot{\gamma} \left[ 2D\_{\rm r} - 2tr(D\_{\rm r})A - 5D\_{\rm r}A - 5A \; D\_{\rm r} - 10A\_{4} \right] \tag{9}$$

.

In addition to the second order orientation tensor, the fourth order orientation tensor is required for the calculation, and is solved with invariant-based closure approximation (IBOF5) according to Chung and Kwon [90], which is a model that o ffers a good compromise of a precise approximation and computation e ffort.

.

The rotary diffusion tensor *Dr* (Equation (11)) is determined by the scalar *D*<sup>2</sup> (Equation (10)) of the square of the symmetrical part *D* from the velocity gradient tensor *L*.

$$\|\|D^2\|\| = \sqrt{\frac{1}{2}D^2 : D^2} \tag{10}$$

$$D\_r = \mathbb{C}\_I \Big( I - \mathbb{C}\_m \frac{D^2}{\| \| D^2 \| \|} \Big) \tag{11}$$

The second part describes the retarded principal rate (Equation (12))

$$
\dot{A}^{RPR} = -R \, \Lambda^{IOK} \, R^T \tag{12}
$$

with Λ*IOK* (Intrinsic Orientation Kinetics), as the substantial derivative of a certain diagonal tensor. In general, an eigenvalue calculation of both the orientation tensor *A* and its change . *A* is coupled here. This is realized as follows (Equations (13) and (14)):

$$
\Lambda\_{\vec{\mu}}^{\rm IOK} = a \dot{\Lambda}\_i,\text{ with }i = 1,2,3\tag{13}
$$

$$R = [\mathfrak{e}\_1, \mathfrak{e}\_2, \mathfrak{e}\_3] \tag{14}$$

where Λ*i* represents the eigenvalues of the orientation tensor and *R* the rotation matrix, which is composed of the eigenvectors of the change in the orientation tensor. This relationship is extended as follows (Equation (15))

$$A\_{ii}^{\text{IOK}} = a \left[ \dot{\Lambda}\_i - \beta \left( \dot{\Lambda}\_i^2 + 2 \dot{\Lambda}\_j \dot{\Lambda}\_k \right) \right], \text{ with i, j, k = 1, 2, 3} \tag{15}$$

with an additional parameter β, which is defined as a time constant.

In order to determine the different influences of the velocity field more precisely, the parameters *Ci* and α were calculated depending on the shear rate according to Tseng et al. [91]. Figure 9 shows that α was adjusted in the core layers, which damps this parameter at lower shear rates, while the opposite is true for *Ci*, which is damped at higher shear rates.

**Figure 9.** Shear rate-dependent adaption of the parameters *Ci* and α [91].

For this purpose, the parameters are changed as follows (Equations (16) and (17)) [91].

$$\frac{\mathbf{C}\_i(\mathbf{\dot{\nu}})}{\mathbf{C}\_{i0}} = \frac{1}{1 + \left(\frac{\mathbf{\dot{\nu}}}{\mathbf{\dot{\nu}}\_c}\right)^2} \tag{16}$$

$$\frac{\alpha\left(\dot{\boldsymbol{\nu}}\right)}{\alpha\_0} = 1 - \frac{1}{1 + \left(\frac{\dot{\boldsymbol{\nu}}}{\dot{\boldsymbol{\nu}}\_c}\right)^2} \tag{17}$$

The α0 and *Ci*0 input values are given by the user and then changed to a specific value for each element of the simulation model. With . γ*c*, a critical shear rate is defined, which has to be determined empirically. In commercial simulation tools like Moldex3D ® R17 this variable is implemented as a constant determined by experience of the developers.

#### 1.1.9. Calibration of Fiber Orientation Models

The model parameters play a crucial role in the use of fiber orientation models, because only with an appropriately calibrated model the fiber orientation can be reliably predicted [7,40,88,92–94]. The main challenge here is that the model parameters are not exclusively material-dependent, but may even depend on the process and part geometry [94,95]. Therefore, standard parameters can only be used rarely or only for a few specific cases. Especially for LFTs, a reliable prediction with current models is not possible when using standard parameters.

Most publications concerning the prediction of fiber orientation either use standard parameters [2] or simply adjust or fine tune them manually by trial and error [1,96]. Sometimes there is a so-called fit, but no algorithm for the parameter identification is described [8,88,92,97,98]. Independently of the geometry, objective fitting methods has been described for simple shear flows of short fiber [99] and also long fiber [100] reinforced plastics. Some of the necessary model parameters like the fiber-fiber interaction coe fficient *Ci* were extensively analyzed and associated approximation functions [95,101,102] were developed, respectively. The influence of the parameters of recent fiber orientation models [87,92,103,104] or viscosity [6,105] has been investigated in some publications, but basic correlations to the material or approximated functions to evaluate the model parameters are still missing.

Several methods are available to optimize the necessary parameters of the prediction models for the microstructure of fiber reinforced plastics. In a first step, however, the properties of the microstructure must be determined, or methods as described in Section 3.2 must be developed in order to compare certain properties with each other directly [106]. Scalar values, such as the spatial fiber concentration, are less problematic than the fiber orientation and length, which are subject to a certain distribution even in local areas. If these requirements are satisfied, the parameters of the respective prediction models can be optimized by means of a reverse engineering procedure, as shown in Figure 10. However, even though it is possible to automatically optimize the parameters of the fiber orientation model, only a few algorithms are described in the literature [94,106–109].

**Figure 10.** Reverse engineering of fiber microstructure properties at a specific position.

#### 1.1.10. Optimization Methods

The aim of this study is to predict the experimentally determined fiber orientation with minimum error variance by optimizing the parameters of the fiber orientation model ( *Ci*, *Cm*, and α for example for the iARD model, see Section 1.1.8). According to the previous Section 1.1.9, the most suitable calibration of the fiber orientation model should be identified. Since the intention of this work is not the invention of an improved optimization algorithm, the basics of optimization methods will only be described in principle. A detailed description and a comparison of di fferent approaches to optimization will be omitted and reference will be made to further literature [110,111].

In general, parameter optimization deals with the problem of finding a parameter vector *x*<sup>∗</sup> for which a given target function is minimized. Non-linear variation methods of the form *min f*(*x*) are solved. This function can be approximated by quadratic Taylor series approximation Equation (18).

$$\min f(\mathbf{x}) \approx h(\mathbf{x}) := f(\overline{\mathbf{x}}) + \nabla f(\overline{\mathbf{x}})^T (\mathbf{x} - \overline{\mathbf{x}}) + \frac{1}{2} (\mathbf{x} - \overline{\mathbf{x}})^t H(\overline{\mathbf{x}}) (\mathbf{x} - \overline{\mathbf{x}}) f(\mathbf{x}) \tag{18}$$

where ∇ *f*(*x*) represents the gradient of *f*(*x*), and *<sup>H</sup>*(*x*) represents the Hessian matrix of *f*(*x*). The gradient of *h*(*x*) is calculated as follows by finding the minimum of the gradient (Equation (19))

$$0 = \nabla f(\overline{\mathbf{x}}) + H(\overline{\mathbf{x}})(\mathbf{x} - \overline{\mathbf{x}}) \tag{19}$$

For the use of the Newton method for problems with constraints, an additional state variable α*C* is added to the original Equation, which leads to the following variation Equation (20).

$$\Psi(\mathbf{x} - \overline{\mathbf{x}}) = -a\_{\mathbb{C}} H(\overline{\mathbf{x}})^{-1} \nabla f(\overline{\mathbf{x}}) \tag{20}$$

In order to determine α*C*, various line search techniques or methods can be used to directly satisfy the decisive variable constraints. In addition, the Hessian matrix can be approximated by a diagonal matrix *D H* This simplifies the optimization to a quasi-Newton method [112,113].

#### *1.2. Objectives and Hypotheses*

According to the major challenge to reliably predict the behavior of discontinuously fiber reinforced parts produced by injection molding, a holistic approach to predict all considerable properties of the fiber microstructure is desirable. In a first step, the following hypotheses are drawn and investigated to improve the predicted fiber orientation. In the future, these investigations can be extended to the remaining microstructural properties such as fiber length and fiber concentration.

On the one hand, the prediction accuracy of the fiber orientation distribution can be significantly improved by using optimization methods. On the other hand, the required optimization time can be significantly reduced with a new calibration approach.

In order to optimize the parameters of the fiber orientation model to predict the experimentally determined fiber orientation with a minimum error variance, two di fferent calibration approaches are investigated. However, since there have been several further developments of the iARD fiber orientation model in the last few years; the most suitable version is initially identified. In the first calibration approach, a commercial simulation software is used to calculate the fiber orientation model, and only the parameters of the fiber orientation model are externally optimized in another tool to obtain the most accurate prediction for the fiber orientation distribution. As second approach, a new method is developed to calibrate the fiber orientation model. The aim of the new approach is to achieve the optimization goal as quickly and e fficiently as possible. Therefore, in contrast to the first approach, which requires a time-consuming recalculation of the entire flow field with each optimization step, the new approach is based on only one initial calculation of the flow field. The subsequent parameter optimization of the fiber orientation model is performed without any further calculation of the flow field. With this approach, a suitable calibration of the fiber orientation model should be achieved within a few minutes.

Even though the new approach focuses on the calibration of the fiber orientation model, the method should be designed in a way that it can be easily transferred to other models to predict microstructural properties. For example, in the future, this approach should also be used to calibrate the models to predict fiber length and fiber concentration.

#### **2. Experimental Setup**

#### *2.1. Materials, Parts, and Processing*

The fiber orientation distribution was experimentally determined for injection molded plates according to ISO 294-5 with dimensions of 80 × 80 × 4 mm, as shown in Figure 11a.

**Figure 11.** Injection molded plates (**a**) and scheme of evaluation points along one side of the plate (**b**).

A long fiber reinforced polypropylene (PP-LGF) from TechnoCompound GmbH, Bad Sobernheim, Germany was chosen, because it is the most frequently used LFT material [114]. The granules are 10 mm long, rod-shaped pellets with fiber mass fractions of 20%, 40%, and 60% of the product type TechnoFiber. The plates were injection molded with two constant injection speeds of 30 cm<sup>3</sup>/s (low) and 100 cm<sup>3</sup>/s (high), a melt temperature of 230 ◦C and a mold temperature of 80 ◦C on an injection molding machine, Allrounder 520S 1600–400 from Arburg, Lossburg, Germany.

For the experimental determination of the fiber orientation distribution, small samples (10 mm × 10 mm) were taken at the beginning (A), middle (B), end (C) and aside (D) of the flow path according to Figure 11b. These samples were analyzed by inspecting the cross sections highlighted with red in Figure 11, using optical reflection microscopy as described in the following section.

#### *2.2. Experimental Determination of the Fiber Orientation Distribution*

With the image processing program FIJI (ImageJ), developed by Wayne Rasband, in combination with MATLAB® by MathWorks®, Natick, USA, spatial fiber orientation for the individual cross section were determined using the fiber orientation tensor *A* as explained in Section 1.1.5. First of all the high-resolution subsection images (1.920 × 1.200 px) of each cross section (ca. 150) were stitched together and converted to a binary image by means of different image processing algorithms (image stitching, background subtraction, local threshold, watershed) as shown in Figure 12a before and b after digital image processing. Based on these binary images a particle analysis is performed in ImageJ. By means of the second order central moment to fit the best ellipse, the major and minor axis of each particle is determined. Afterwards, the fiber orientation tensor of each ellipsis was calculated in a MATLAB® tool by using the equations given in Figure 8c, developed by [60], and corrected with the weighting function according to [40]. In order to analyze the fiber orientation distribution along the thickness with sufficient accuracy, the whole cross section was discretized along the z-direction into a structured grid with at least twenty individual fibers per cell. Subsequently, the fiber orientation tensors of all fibers within a discretized cell were averaged according to the cells marked with red in Figure 12.

(**a**) grayscale (**b**) binary 

**Figure 12.** Different states of image processing to determine the fiber orientation.

The fiber orientation distribution along the sample thickness (z-direction) of each individual sample is obtained from overall 15 grid points by the calculated cell averages. In Addition, two cross sections of each material and process setting were analyzed and averaged to eliminate measurement and processing deviations. For example, Figure 13a shows the fiber orientation results for two samples and the averaged curve.

**Figure 13.** Processing of the fiber orientation results for objective comparison (PP-LGF20 v30).

Even though two samples were averaged, an objective calibration of the fiber orientation model on the basis of only few measuring points is usually not reasonable. Furthermore, the experimentally determined fiber orientation shows a slightly off-center position of the core layer, which is caused by an uneven mold cooling as already observed in [101]. To obtain comparable results for experiment and simulation, the measuring points are symmetrized towards the center and smoothed by approximation. However, the characteristic properties of the curve progression regarding the fiber properties in the core and shear layer should be retained. For this purpose, the measuring points are represented by a fifth-order Taylor series approximation, as shown in Figure 13b for the orientation tensor entry *A*11. This method allows an objective comparison between simulation and experiment, independently of the number of measured data points. In addition, the curves are symmetrized without a significant change in their characteristic shape.

#### *2.3. Process Simulation*

The process simulation was performed with Moldex3D® of CoreTech System, Chupei City, Taiwan, which numerically solves the governing Equations (2)–(4). The calculated geometry is a complete model of the injection molded plate as shown in Figure 11a. The mesh was made by Rhinozeros 5 from Robert McNeel & Associates, Seattle, USA, and is shown in Figure 14a. The plates were meshed with a structured grid consisting of 21 hexahedral elements along the thickness of 4 mm, whereas tetrahedral elements were used for the runner system and sprue. A constant velocity according to the injection speed was set at the Inlet. The viscosities of the individual materials PP-LGF20, PP-LGF40, and PP-LGF60 were determined by means of a high-pressure capillary rheometer in the range of high shear rates, and a rotational rheometry for low shear rates. The experimental data was fitted by the Herschel–Bulkley viscosity model as shown in Figure 14b.

**Figure 14.** Geometry and material setup of the simulation model.

The results of the viscosity determination show the expected rise of the viscosity with higher fiber content, especially for relatively low shear rates. The respective model parameters for the prediction of the fiber orientation model are determined by different iterative optimization methods described in the following sections, which compare the simulation results to the experimental measured data. The pvT-model for the thermodynamic state relations according to [115] was used to describe the density of the fluid as a function of temperature and pressure. All other necessary boundary conditions, such as injection speed as well as melt and mold temperatures, were defined according to the injection molding processing conditions described in Section 2.1.

In order to investigate the influence of the different versions of the iARD fiber orientation model as explained in Section 1.1.8 on the resulting fiber orientation distribution the following software versions of Moldex3D® with the different fiber orientation model versions


Because fiber length is significantly reduced within the plastification unit, and this is not part of the simulation, the number averaged fiber length (see Table 1) was experimentally measured as described in [36] at the sprue and specified as a boundary condition in the respective process simulations. As expected, the results show that, with an increasing fiber content, the probability of fiber–fiber interaction increases significantly, and fiber breakage occur more frequently.


**Table 1.** Experimentally determined number average fiber length at the sprue.

#### *2.4. Adjustment of the Fiber Orientation Model within Process Simulation*

#### 2.4.1. Comparison of Simulation and Experimental Data

In order to objectively compare the predicted fiber orientation distribution along the thickness *z* with the experimentally determined data, the least square error *S* is calculated by Equation (21) at all positions of a sample according to Figure 11b.

$$S = \sum\_{z=0}^{d} \left( A\_{\exp}(z) - A\_{\text{sim}}(z) \right)^2 \tag{21}$$

To reduce the optimization problem, only the three orientation tensor components *A*11, *A*22, and *A*33 were added as least square functions. Furthermore, the minimum square error functions can be weighted differently if the adaptation to a certain orientation direction takes priority.

#### 2.4.2. Influence of the Fiber Orientation Model Parameters

Another considerable factor is the accurate calibration of the model parameters *Ci*, *Cm*, and α, which can influence the core-layer thickness (mainly α) as well as the degree of orientation in the skin layers (mainly *Ci* and *Cm*). In comparison of different fiber orientation model versions implemented in Moldex3D® as described in Section 2.3 the influence of the three parameters changes significantly. As shown in Figure 15a change of *Ci* and *Cm* within the basic iARD model mainly result in a different degree of orientation inside the skin layers, whereas α leads to a change in the core layer thickness. With more recent enhancements of the orientation model iARD with shear rate dependency in R16 and also with fiber coupled viscosity in R17, the parameter α shows hardly any change in the core layer thickness Figure 15b,c. This also leads to a very small transition zone between the core and skin layers.

**Figure 15.** Influence of fiber orientation model parameters on different software versions.

Even though the basic iARD fiber orientation model implemented in Moldex3D® version R13 is considerably older than the other versions, only this version allows to adjust the core layer thickness by means of the corresponding fiber orientation model parameters. According to this study, the basic iARD model has the most individual configuration parameters to represent the experimentally determined fiber orientation as accurately as possible.

#### *2.5. Calibration of Fiber Orientation Model*

Based on the calibration and optimization methods explained in Sections 1.1.9 and 1.1.10, a detailed description of the methods used in this article is given below.

#### 2.5.1. Direct Optimization

In order to directly optimize the results of the fiber orientation model iARD, the respective model parameters *Ci*, *Cm* and α have to be adjusted until the defined optimization objective is reached. However, the optimization algorithms require input variables as well as output variables, which can be described by the functions or directly adjusted. Therefore, a specific routine is needed, which changes the fiber orientation model parameters inside the simulation environment, starts the simulation run with the defined parameters, waits until the simulation is successfully calculated, reads out the defined results (i.e., fiber orientation tensor component *A*11 along the thickness at a specified location of the part), and transfers these results as input parameters to the optimization algorithm, and compares the results with a defined optimization objective. Finally, depending on the optimization objective and the associated accuracy, the optimization procedure must be terminated or the routine must be restarted. In this work, the described optimization approach was implemented in a MATLAB® script, which directly accesses the input data of Moldex3D®, and thus creates new runs and triggers their calculation.

#### 2.5.2. A New Calibration Method

In a new method, the fiber orientation calculation and its optimization were performed with a MATLAB® script according to the process scheme shown in Figure 16.

**Figure 16.** Process scheme of the fiber orientation optimization tool.

To obtain the necessary flow field for the injection molding process, the results (velocity gradient tensor *L*) of a single initial injection molding simulation in Moldex3D® were exported and used in MATLAB® to calculate the corresponding fiber orientation model by means of Equations (5)–(15) (basic iARD model). Based on the initial calculation of the flow field, the whole parameter optimization of the fiber orientation model was carried out within a few minutes.

The main intention of this method was not the development a new optimization algorithm, but the removal of the time-consuming recalculation of the entire flow field for each optimization step (in contrast to the direct optimization described in Section 2.5.1). By using the flow field of an initial injection molding simulation, the computation is reduced to the calculation of the fiber orientation model.

The goal of the optimization was to minimize the objective function, consisting of the error sum of squares between simulative calculated and experimental determined fiber orientations at pre-defined positions. A gradient-based algorithm according to Section 1.1.10 was used for this purpose. This means that at the current point of the optimization, the first derivative for the direction of the optimization step and second derivatives for the length of the iteration step for each variable are used.

Because this case is a multi-parameter optimization, a system of Equations (22) is required to solve a variation step. This is as follows for the three parameters of the fiber orientation model.

$$
\begin{bmatrix} \mathbb{C}\_{i} \\ \mathbb{C}\_{m} \\ \mathbb{C} \end{bmatrix} = \begin{bmatrix} \partial^{2} \mathcal{S} \begin{pmatrix} \mathbb{C}\_{i} \end{pmatrix} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \partial^{2} \mathcal{S} \begin{pmatrix} \mathbb{C}\_{m} \end{pmatrix} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \partial^{2} \mathcal{S} \begin{pmatrix} \mathbb{C} \end{pmatrix} \end{bmatrix}^{-1} \begin{bmatrix} \partial \mathcal{S}(\mathbb{C}\_{i}) \\ \partial \mathcal{S}(\mathbb{C}\_{m}) \\ \partial \mathcal{S}(\infty) \end{bmatrix} \tag{22}
$$

In this form, only the main diagonals of the matrix are filled, which changes the optimization to a quasi-Newtonian method. In this form, fewer derivatives are needed to calculate the variation vector, which makes the calculation much more timesaving. Furthermore, the addition of scaling is also simplified. The Hessian matrix is positive definite, and thus the correct direction of the variation step is determined. Additionally, all three parameters are assigned to limits.


For this reason, an additional line search method is introduced to guarantee that these limits are observed. To prevent the optimization from getting stuck at these limits, a damping factor <sup>∝</sup>*D* is introduced in Equation (23) in addition to the line search factor <sup>∝</sup>*k*.

$$\begin{bmatrix} \mathbb{C}\_{i}^{\*} \\ \mathbb{C}\_{m}^{\*} \\ \alpha^{\*} \end{bmatrix} = \begin{bmatrix} \mathbb{C}\_{i} \\ \mathbb{C}\_{m} \\ \alpha \end{bmatrix} \otimes\_{k} \begin{array}{c} \mathbb{c}\_{k} \\ \mathbb{c}\_{k} \end{array}, \text{with } \mathbb{C}\_{i}^{\*} = \mathbb{C}\_{i, i-1} + \infty\_{k} \begin{array}{c} \mathbb{C}\_{i} \geq 0.1 \text{ and } \begin{array}{c} \infty\_{k} = \frac{0.1 - \mathbb{C}\_{i, i-1}}{\mathbb{C}\_{i}} \\ \mathbb{C}\_{i} \end{array} \tag{23}$$

From all restrictions, the value <sup>∝</sup>*k* is determined. Furthermore, a line search factor is used to prevent oscillation. For this purpose, it is verified whether the target function value increases contrary to the expectation during the optimization. In this case, the optimization step is shortened greatly with <sup>∝</sup>*k* to allow only a small deterioration. For this purpose, the factor <sup>∝</sup>*k*= 0.05 is specified if *Si* > *Si*−1. Because the optimization problem can be dominated by the restrictions, the following is suggested for the damping factor:


This is to prevent the optimizer from getting stuck at the given variable boundaries. This problem is further improved by the introduction of scaling for the Hessian matrix (Equation (24)).

$$
\begin{bmatrix} \mathbb{C}\_{i} \\ \mathbb{C}\_{m} \\ \alpha \end{bmatrix} = \begin{bmatrix} f\_{1} \ \partial^{2} \mathcal{S} \ (\mathbb{C}\_{i}) & 0 & 0 \\ 0 & f\_{2} \ \partial^{2} \mathcal{S} \ (\mathbb{C}\_{m}) & 0 \\ 0 & 0 & f\_{3} \ \partial^{2} \mathcal{S} \ (\alpha) \end{bmatrix}^{-1} \begin{bmatrix} \delta \mathcal{S} (\mathbb{C}\_{i}) \\ \delta \mathcal{S} (\mathbb{C}\_{m}) \\ \delta \mathcal{S} (\alpha) \end{bmatrix}
\\
f\_{2} = f\_{3} = \left(0.1 - \mathbb{C}\_{i}^{\*}\right)^{-\mathrm{m}} \tag{24}
$$

With the damping factor <sup>∝</sup>*D*, the barrier is approached more slowly, and the other two parameters are weighted more heavily.
