# Open-Source Framework for Modeling the Evolution of Fiber Orientation

Bruno Ramoa, Ricardo Costa, Francisco Chinesta, João Miguel Nóbrega

► To cite this version:

Bruno Ramoa, Ricardo Costa, Francisco Chinesta, João Miguel Nóbrega. Open-Source Framework for Modeling the Evolution of Fiber Orientation. OpenFOAM journal, 2025, 5, pp.17-37. 10.51560/ofj.v5.131. hal-04949845

HAL Id: hal-04949845

https://hal.science/hal-04949845v1

Submitted on 20 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# An Open-Source Framework for Modeling the Evolution of Fiber Orientation

B. Ramoa<fnref n="1" />,*, R. Costa<fnref n="1" />●, F. Chinesta<fnref n="2" />, and J.M. Nóbrega<fnref n="1" />●

**DOI:** https://doi.org/10.51560/ofj.v5.131

**Results with version(s):** OpenFOAM®v2212

**Repository:** https://github.com/Computational-Rheology/paper_OFJ_FiberOrientation/tree/main

**Abstract.** Modeling fiber orientation plays a crucial role in predicting the behavior of fiber reinforced thermoplastic materials. The equation that governs the evolution of the fiber orientation is hyperbolic in nature and requires handling of fourth-order tensors, which are currently unavailable in OpenFOAM®. The current work explores the possibility of using OpenFOAM® and open-source symbolic computation for modeling the evolution of fiber orientation. For this purpose, a *functionObject* was programmed to work as a plug-in for any OpenFOAM® incompressible flow solver and compute the evolution of fiber orientation in a decoupled manner. Several fiber orientation models and closure relations available in the literature were implemented in the tool, which were verified by comparing their predictions with independent results for a single material point obtained by numerically integrating the associated governing equations.

## 1. Introduction

Thermoplastic materials are tied to our daily lives, encompassing a wide range of applications, from household products to advanced technical appliances. To meet the permanent demand for enhanced products, fiber-reinforced thermoplastic materials (FRTM) are often employed to achieve specific properties (property weighted by the density) that cannot be obtained with conventional materials. Industries such as automotive heavily rely on this technology to achieve lightweight designs, essential for optimizing fuel efficiency [1].

The final properties of thermoplastic parts depend on the microstructure developed during processing. For FRTM, the effects of fiber length, concentration, and orientation play a significant role in the performance of the final part. Particularly, fiber orientation has been extensively studied throughout the last decades [2–6], leading to the development of several phenomenological models [7–13] based on the seminal work of Jeffery [14] to quantify the state of orientation in FRTM parts. During the product development phase, knowledge of this final state of orientation is crucial for structural analysis and design, as it directly affects the parts’ mechanical performance [15–17].

In the context of OpenFOAM® only a few works have been dedicated to fiber orientation modeling, where only the Folgar-Tucker (FT) [7] and reduce strain closure (RSC) [8] models in tensor-based formulations have been reported [18–21]. Given the open-source character and outreach of OpenFOAM®, this work presents a novel strategy for modeling the fiber orientation evolution in incompressible fluid flows using the OpenFOAM® framework. The developed plug-in can be utilized for injection molding, compression molding, and any other solver that involves fluid flow. The objective is to incorporate the capability of predicting fiber orientation as an extension, without making modifications to an existing solver. The approach involves creating a *functionObject* to solve the three-dimensional orientation evolution equation with state-of-the-art models and closure approximations in a decoupled manner, using a one-way interaction, where the velocity field affects the fiber orientation evolution but the computed orientation will not influence the relation between the velocity gradient and the stress tensor (τ), i.e., the constitutive model, required to calculate the velocity field. Nevertheless, the strategy employed is

<footnote n="1">¹Institute for Polymers and Composites, University of Minho, Campus de Azurém, 4800-058 Guimarães, Portugal
*Email address: bruno.ramoa@dep.uminho.pt*</footnote>

<footnote n="2">²ESI Group Chair, PIMM Institute at ENSAM ParisTech, 155 Boulevard de l'Hôpital, F-75013 Paris, France</footnote>

<footnote>* Corresponding author
Received: 21 December 2023, Accepted: 21 June 2024, Published: 20 January 2025</footnote>
flexible enough for the interested user to extend the developed classes to account for the effect of fiber orientation into the stress tensor calculation [22].

The remaining paper is organized as follows: Sections 2 and 3 introduce, respectively, the fiber orientation modeling approaches and the state-of-the-art models together with the closure approximations used in commercial applications. In Section 4, the proposed approach to implement fiber orientation modeling within the OpenFOAM® framework is presented. Sections 6 and 7 showcase and discuss case studies involving simple-shear and the flow in a center-gated disk, tested to validate the numerical implementations. Finally, Section 8 summarizes the main conclusions drawn from this work.

## 2. Fiber orientation modeling with macroscopic descriptors

The employed formulation considers rigid, axi-symmetric, straight fibers, such as glass fibers typically used to manufacture FRTM [9], whose orientation characterization can be carried out at different scales. At the micro-scale, the orientation of individual fibers can be used as a complete descriptor. This can be done through a pair of $\theta$ and $\phi$ angles or by a unit vector $\underline{p}$ that is aligned with the fiber axis, as illustrated in Fig. 1. Both descriptions are related by:

$$ \underline{p} = \begin{bmatrix} p_1 \\ p_2 \\ p_3 \end{bmatrix} = \begin{bmatrix} \sin \theta \cos \phi \\ \sin \theta \sin \phi \\ \cos \theta \end{bmatrix}. \qquad (1) $$

<figure><img src="image_2.png" /><figcaption>Figure 1. Schematic of the orientation of a single fiber.</figcaption></figure>

The motion of a single, rigid, ellipsoidal fiber in a Newtonian fluid was formulated by Jeffery, assuming that the: (i) velocity field was not influenced by the fibers; (ii) inertial and body forces were negligible; (iii) the fluid velocity varied linearly with position and iv) the fluid where the fiber is suspended is incompressible [14, 23–25]. Under these conditions and in the absence of additional forcings, the equation for describing the rate-of-change of the fiber orientation is given by [25]:

$$ \dot{\underline{p}} = \underline{W} \cdot \underline{p} + \xi (\underline{D} \cdot \underline{p} - \underline{D} : \underline{p} \underline{p} \underline{p}), \qquad (2) $$

or, equivalently, in tensor notation as [25]:

$$ \dot{p}_i = W_{ij}p_j + \xi (D_{ij}p_j - D_{jk}p_k p_j p_i), \qquad (3) $$

where:

* $\dot{\underline{p}}$ is time rate-of-change of $\underline{p}$,

* $\underline{W}$ is the vorticity tensor, given by: $\underline{W} = \frac{1}{2}(\underline{L} - \underline{L}^T)$,

* $\underline{D}$ is the rate-of-deformation tensor, given by: $\underline{D} = \frac{1}{2}(\underline{L} + \underline{L}^T)$,

* $\underline{L}$ is the velocity gradient tensor, given by $L_{ij} = \frac{\partial v_i}{\partial x_j}$,

* $v_i$ are the components of the velocity vector, $v$,

* $\xi$ is the particle shape factor, given by $\xi = \frac{a_r^2 - 1}{a_r^2 + 1}$, and
* $a_r$ is the fiber aspect-ratio

Equation (3) reflects the effect of the flow on the fiber motion. The first term on the right-hand side accounts for the influence of rigid-body rotation and the second (between parenthesis) for the effect of fluid deformation on fiber orientation [25]. The term $D_{jk}p_k p_j p_i$ will result in a vector that is parallel to $p_i$ and scaled by $D_{jk}p_k p_j$. This subtracts from $D_{ij}p_j$ the component that is parallel to $p_i$, thus ensuring that $\dot{p}_i$ is always perpendicular to $p_i$ and that the magnitude of $p_i$ remains constant [25]. It is also known that Eqn. (3) is kinematically reversible. Meaning that the flow field might affect the fiber position and orientation, but if the flow is reversed the fiber will return to its initial position and orientation [26].

In FRTM, it is common to have several thousand fibers per cubic millimeter, consequently, keeping track of each orientation vector is a massive and mostly unnecessary task [25]. Thus, at the mesoscopic scale, statistical methods are used to describe the state of orientation for collections of fibers, through the use of a probability density function $\psi(p, t)$, also known in literature as an orientation distribution function (ODF) [27]. The evolution of the ODF is governed by a Fokker-Planck equation [9,27], reading:

$$ \frac{D\psi}{Dt} = -\nabla_s \cdot (\psi\dot{p}), \qquad (4) $$

where:

* $\frac{D\psi}{Dt}$ is the material derivative of $\psi$ in canonical space,

* $\nabla_s$ is the gradient operator evaluated at the unit sphere's surface, and

* $\psi\dot{p}$ represents the probability flux.

Equation (4) has several properties, and the reader is referred to references [25, 28, 29] for a more complete description. Since $\underline{p}$ is a unit vector, all possible combinations of $(\theta, \phi)$ or $\underline{p}$ are located on the surface of the unit sphere (conformational space). In this work, only the normalization condition of the ODF is reported. It states that the integral of $\psi$ over the conformational space is unitary [25], mathematically:

$$ \begin{gathered} \oint \psi(\underline{p}) d\underline{p} = 1 \\ \text{or} \\ \int_{\phi=0}^{2\pi} \int_{\theta=0}^{\pi} \psi(\theta, \phi) \sin(\theta) d\theta d\phi = 1. \end{gathered} \qquad (5) $$

Collisions between fibers will result in a diffusion-like behavior [26]. In 1984, Folgar and Tucker [26] proposed a new model based on Jeffery's work to account for fiber interactions by adding a phenomenological diffusive flux $q$, which reads [9]:

$$ \frac{D\psi}{Dt} = -\nabla_s \cdot (\psi\dot{p} + q), \qquad (6) $$

with:

$$ q = -D_r \nabla_s \psi, \qquad (7) $$

where $D_r$ is a rotary diffusivity. The authors proposed a scalar rotary diffusivity, which has since been known in the literature as the isotropic rotary diffusion (IRD) model and is given by:

$$ D_r = C_I \dot{\gamma}, \qquad (8) $$

where:

* $C_I$ is a phenomenological parameter modeling the randomizing effect of fiber-fiber interactions, and

* $\dot{\gamma}$ is the magnitude of the rate-of-deformation tensor, given by $\dot{\gamma} = \sqrt{2\underline{D}:\underline{D}}$.

The ODF offers keen insights as an orientation descriptor, however, its usage for computer-assisted engineering (CAE) in three-dimensional simulations is computationally expensive. To minimize this difficulty, a macroscopic orientation descriptor was proposed by Advani and Tucker [7] in 1987, where the information contained in the ODF was condensed into a tensor-based formulation, by operating on the moments of the ODF. The most commonly used macro-descriptor is defined as [9,27]:
$$ \underline{\underline{A}} = \oint \underline{p} \underline{p} \psi(\underline{p}) d\underline{p} $$
or
$$ A_{ij} = \int_{\phi=0}^{2\pi} \int_{\theta=0}^{\pi} p_i p_j \psi(\theta, \phi) \sin(\theta) d\theta d\phi. \quad (9) $$

This second-order tensor is symmetric and has a unitary trace. The tensor-based formulation of the model proposed by Folgar and Tucker (FT), Eqn. (6), reads [9]:

$$ \dot{\underline{A}} = \dot{\underline{A}}^{\text{H}} + \dot{\underline{A}}^{\text{IRD}}, \quad (10) $$

with:

$$ \dot{\underline{A}}^{\text{H}} = \underline{W} \cdot \underline{A} - \underline{A} \cdot \underline{W} + \xi (\underline{D} \cdot \underline{A} + \underline{A} \cdot \underline{D} - 2\underline{A} : \underline{D}), \quad (11) $$

$$ \dot{\underline{A}}^{\text{IRD}} = 2C_I \dot{\gamma} (\underline{I} - 3\underline{A}), \quad (12) $$

where:

* $\dot{\underline{A}}$ is the material derivative of $\underline{A}$,

* $\dot{\underline{A}}^{\text{H}}$ is the hydrodynamic component of motion that results from Jeffery's model – the first term of the right-hand side of Eqn. (6),

* $\dot{\underline{A}}^{\text{IRD}}$ is the isotropic rotary diffusion – the second term of the right-hand side of Eqn. (6),

* $\underline{A}$ is a fourth-order orientation tensor, and

* $\underline{I}$ is the second-order identity tensor.

Although computationally more efficient than the ODF, the use of orientation tensors have some drawbacks:

* The reconstruction of the orientation distribution function from $\underline{A}$ is ambiguous, as $\underline{A}$ might not represent a unique orientation state [28]. Figure 2 illustrates different fiber configurations that result in the same orientation tensor.

* In the derivation of the transport equation for $\underline{A}$, the next even-order tensor appears. To close the system, a closure relationship is needed, and several authors have proposed solutions to write $A_{ijkl} = f(A_{ij})$. In Section 3, the most relevant closure relations used in CAE will be presented.

* When an experimentally determined $C_I$ is unavailable, some authors use semi-empirical expressions [21, 30] based on the product of the fiber volume fraction and the fiber aspect-ratio.

Nevertheless, this descriptor retains the generality of the ODF and provides a compact, computationally efficient and convenient representation of the fiber orientation state, being the state-of-the-art descriptor used in CAE applications.

Over the last four decades several modifications have been proposed to the FT model. Different studies have shown that the orientation kinetics predicted by the FT model evolve much faster than the ones observed experimentally [9]. Accordingly, aiming to achieve a better agreement with the transient experimental observations, Wang et al. [8] proposed a new objective model (independent of the considered coordinate system) designated reduced-strain closure (RSC). This model decomposes $\dot{\underline{A}}$ into rate equations for both the eigenvalues ($\lambda_i$) and eigenvectors ($e_i$). The approach consists in slowing down the evolution for the eigenvalues by a phenomenological scalar parameter ($\kappa \le 1$), while leaving the rate equations for the eigenvectors unchanged [25]. The RSC model reads:

$$ \dot{\underline{A}} = \dot{\underline{A}}^{\text{RSC}} + \kappa \dot{\underline{A}}^{\text{IRD}}, \quad (13) $$

with:
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Illustration of fiber orientation states and the corresponding fiber orientation tensor, $\underline{\underline{A}}$.</figcaption></figure>

$$ \dot{\underline{A}}^{\text{RSC}} = \underline{\underline{W}} \cdot \underline{\underline{A}} - \underline{\underline{A}} \cdot \underline{\underline{W}} + \xi (\underline{\underline{D}} \cdot \underline{\underline{A}} + \underline{\underline{A}} \cdot \underline{\underline{D}} - 2[\underline{\underline{A}} + (1-\kappa)(\underline{\underline{L}} - \underline{\underline{M}} : \underline{\underline{A}})] : \underline{\underline{D}}), \quad (14) $$

$$ \underline{\underline{A}} = \sum_{i=1}^{3} \lambda_i (\underline{\underline{e}}_i \underline{\underline{e}}_i), \qquad (15) $$

$$ \underline{\underline{L}} = \sum_{i=1}^{3} \lambda_i (\underline{\underline{e}}_i \underline{\underline{e}}_i \underline{\underline{e}}_i \underline{\underline{e}}_i), \qquad (16) $$

$$ \underline{\underline{M}} = \sum_{i=1}^{3} (\underline{\underline{e}}_i \underline{\underline{e}}_i \underline{\underline{e}}_i \underline{\underline{e}}_i). \qquad (17) $$

The FT and RSC-IRD models work reasonably well for short FRTM [8]. However, when dealing with long FRTM, their accuracy reduces, as longer fibers tend to align less in the flow direction [9]. To better match experimental data, several authors [9–12] proposed the use of a tensorial rotary diffusion term, to make it direction-dependent, $\underline{\underline{C}}$. The general framework, known as anisotropic rotary diffusion (ARD), was proposed by Phelps and Tucker (PT) [9], and reads:

$$ \dot{\underline{A}} = \dot{\underline{A}}^{\text{H}} + \dot{\underline{A}}^{\text{ARD}}, \qquad (18) $$

with:

$$ \dot{\underline{A}}^{\text{ARD}} = \dot{\gamma} [2\underline{\underline{C}} - 2\text{tr}(\underline{\underline{C}})\underline{\underline{A}} - 5(\underline{\underline{C}} \cdot \underline{\underline{A}} + \underline{\underline{A}} \cdot \underline{\underline{C}}) + 10\underline{\underline{A}} : \underline{\underline{C}}]. \qquad (19) $$

Phelps and Tucker [9] proposed $\underline{\underline{C}} = f(\underline{\underline{A}}, \underline{\underline{D}})$ following a second-order polynomial as:

$$ \underline{\underline{C}} = b_1 \underline{\underline{I}} + b_2 \underline{\underline{A}} + b_3 \underline{\underline{A}}^2 + \frac{b_4}{\dot{\gamma}} \underline{\underline{D}} + \frac{b_5}{\dot{\gamma}^2} \underline{\underline{D}}^2, \qquad (20) $$

where the dimensionless parameters, $b_i$, $i = 1, \dots, 5$, are obtained by fitting experimental data. However, it is known that certain values of $b_i$ lead to unstable and non-physical solutions [25]. Additionally, this model simplifies to the IRD counterpart if $b_1$ is equal to $C_I$, and $b_{2...5}$ to zero. To encompass the need of
matching transient results, the authors linked the ARD framework with the RSC methodology, creating the ARD-RSC model, which is presently available in the proprietary software AutoDesk MoldFlow® [31].

Due to numerical stability problems and fitting dependency requirements of the rotary diffusion tensor proposed by Phelps and Tucker, Tseng and co-workers [10, 32] proposed an alternative model, the improved ARD (iARD), which the authors combined with a new formulation for slowing down the orientation kinetics, known as the retarding principal rate (RPR). Although with a different approach than the one used in the RSC, in essence, the RPR also aims at slowing down the evolution of the eigenvalues without modifying the eigenvectors rate-of-change. Their fiber orientation model reads:

$$ \underline{\dot{A}} = \underline{\dot{A}}^{\text{H}} + \underline{\dot{A}}^{\text{ARD}} + \underline{\dot{A}}^{\text{RPR}}, \qquad (21) $$

with:

$$ \underline{C}^{\text{iARD}} = C_{I}^{\text{iARD}} (\underline{I} - C_{M} \underline{\tilde{L}}), \qquad (22) $$

$$ \underline{\tilde{L}} = \frac{\underline{D}^2}{\|\underline{D}^2\|}, \qquad (23) $$

$$ \underline{\dot{A}}^{\text{RPR}} = -\underline{R} \cdot \underline{\dot{A}}^{\text{IOK}} \cdot \underline{R}^{\text{T}}, \qquad (24) $$

$$ \underline{\dot{A}}^{\text{IOK}}_{ii} = \alpha [\dot{\lambda}_i - \beta (\dot{\lambda}_i^2 + 2\dot{\lambda}_j\dot{\lambda}_k)], \quad i, j, k \in \{1, 2, 3\}, \qquad (25) $$

where:

* $C_M$ is a scalar parameter representing the fiber-matrix interaction,

* $\|\underline{D}^2\|$ is the magnitude of $\underline{D}^2$, given by $\|\underline{D}^2\| = \sqrt{\frac{1}{2}\underline{D}^2 : \underline{D}^2}$,

* $\underline{\dot{A}}^{\text{IOK}}$ is the material derivative of a diagonal tensor calculated under the intrinsic orientation kinetics (IOK) assumption [32],

* $\lambda_i$ is the $i^{th}$ eigenvalue of $\underline{A}$, ordered as: $\lambda_1 \ge \lambda_2 \ge \lambda_3$,

* $\underline{R}$ is an orthogonal rotation matrix whose columns are the eigenvectors of $\underline{A}$ ($\underline{R} = [\underline{e}_1, \underline{e}_2, \underline{e}_3]$),

* $\underline{R}^T$ is the transpose of $\underline{R}$, and

* $\alpha$ and $\beta$ are phenomenological scaling factors used to slow-down the eigenvalues evolution.

The proposed rotary diffusion tensor can be viewed as a special case of Eqn. (20) where only the first and last terms are used [25]. The iARD-RPR model reduces to the FT model by specifying $\alpha = 0$ and $C_M = 0$. Additionally, the RSC and RPR models are identical by setting $\alpha = 1 - \kappa$ and $\beta = 0$ [10]. This model is currently available in the proprietary injection molding simulation software Moldex3D® [33].

Further alternatives to model $\underline{C}$ were made in 2018 by Tseng et al. [11] and Bakharev et al. [12], who proposed a similar fiber orientation model, where the rotary diffusion tensor was considered to be coaxial with $\underline{A}$. Tseng et al. [11] proposed the principal ARD (pARD), in which the rotary diffusion tensor is defined as [27]:

$$ \underline{C}^{\text{pARD}} = C_{I}^{\text{pARD}} (\underline{e}_1 \underline{e}_1 + \Omega \underline{e}_2 \underline{e}_2 + (1 - \Omega) \underline{e}_3 \underline{e}_3), \qquad (26) $$

where $0.5 \le \Omega \le 1$ is a fitted scalar parameter. Bakharev et al. [12] proposed the Moldflow rotational diffusion model (MRD), where $\underline{C}$ is defined as [27]:

$$ \underline{C}^{\text{MRD}} = C_{I}^{\text{MRD}} (D_1 \underline{e}_1 \underline{e}_1 + D_2 \underline{e}_2 \underline{e}_2 + D_3 \underline{e}_3 \underline{e}_3), \qquad (27) $$

where $D_1, D_2$ and $D_3$ are non-negative scalar quantities.

Both $\underline{C}^{\text{pARD}}$ and $\underline{C}^{\text{MRD}}$ have a similar formulation. However, $\underline{C}^{\text{pARD}}$ is less flexible and cannot be reduced to the FT model whilst the $\underline{C}^{\text{MRD}}$ formulation recovers the FT model by setting all parameters equal to unity. The report of Bakharev et al. [12] provides suitable choices for the model parameters to achieve the best correlation with experimental data.

Other fiber orientation models were proposed in the literature, such as Wang's two constant model [25], Kugler et al. [13] flow dependent model, MRD/pARD with RSC formulation for slowing fiber orientation kinetics from Kugler et al. [34] and Kech et al. [35], and even data-driven approaches [36, 37]. The interested reader is referred to [25, 27, 29] for a comprehensive literature review of the most relevant associated topics.
### 3. Closure models

As pointed out in Section 2, one of the drawbacks of using the tensor-based approach as a macro-descriptor of fiber orientation is the need to calculate the next even-ordered orientation tensor, which results from the derivation of the evolution equation of the lower order tensor. Breuer et al. [28] identified in literature 17 closures to describe $A_{ijkl}$ and 4 to identify a sixth-order tensor. However, in commercial applications, only the hybrid, the invariant-based optimal fitting (IBOF), or an orthotropic fitted closure are typically used [33, 38].

The hybrid closure was referred in the review of Advani and Tucker [39] as being a blend between a quadratic and a linear closures. For a 3D case, the model reads:

$$A_{ijkl}^{\text{hyb}} = f A_{ijkl}^{\text{quad}} + (1 - f) A_{ijkl}^{\text{lin}}, \quad (28)$$

where:

* $f$ is the blending factor, which is usually calculated as $f = 1 - 27 \det(A_{ij})$ [39],

* $A_{ijkl}^{\text{quad}} = A_{ij}A_{kl}$, and

* $A_{ijkl}^{\text{lin}} = -\frac{1}{35}(\delta_{ij}\delta_{kl} + \delta_{ik}\delta_{jl} + \delta_{il}\delta_{jk}) + \frac{1}{7}(A_{ij}\delta_{kl} + A_{ik}\delta_{jl} + A_{il}\delta_{jk} + A_{jl}\delta_{ik} + A_{jk}\delta_{il} + A_{kl}\delta_{ij})$, with $\delta_{ij}$ being the second-order identity tensor.

The quadratic closure component is known to be exact when the fibers are perfectly aligned in a single direction, while the linear closure component is known to be exact for isotropic orientation [39]. This blended closure is used due to its small computational cost, when compared to other alternatives. However, it is known to over-predict flow-induced fiber orientation [40–42].

The eigenvalue-based version of the natural closure (ORE) [43] is an orthotropic eigenvalue-based fitted closure. The second-order tensor is decomposed into its eigenvalues and corresponding eigenvectors, with $\lambda_i$ sorted in descending order, $\lambda_1 \ge \lambda_2 \ge \lambda_3$. Due to the normalization condition, Eqn. (5), the summation of the eigenvalues of $\hat{A}$ should equal unity, consequently, only two eigenvalues are independent, typically $\lambda_1$ and $\lambda_2$.

Consider $\hat{A}$ to be a fourth-order tensor defined in the principal axis of $\hat{A}$. Due to the orthotropic character, this tensor can be written with contracted notation having 12 independent parameters as [40]:

$$\hat{A} = \begin{bmatrix} \hat{A}_{11} & \hat{A}_{12} & \hat{A}_{13} & 0 & 0 & 0 \\ \hat{A}_{21} & \hat{A}_{22} & \hat{A}_{23} & 0 & 0 & 0 \\ \hat{A}_{31} & \hat{A}_{32} & \hat{A}_{33} & 0 & 0 & 0 \\ 0 & 0 & 0 & \hat{A}_{44} & 0 & 0 \\ 0 & 0 & 0 & 0 & \hat{A}_{55} & 0 \\ 0 & 0 & 0 & 0 & 0 & \hat{A}_{66} \end{bmatrix}. \quad (29)$$

Assuming a complete symmetry of the fourth-order tensor, there are six independent variables as:

$$\begin{align} \hat{A}_{23} &= \hat{A}_{32} = \hat{A}_{44}, \nonumber \\ \hat{A}_{13} &= \hat{A}_{32} = \hat{A}_{55}, \nonumber \\ \hat{A}_{12} &= \hat{A}_{21} = \hat{A}_{66}. \tag{30} \end{align}$$

Moreover, taking into account that the higher order tensor, $A_{ijkl}$, includes information from the lower order tensor, and since $A_{ijkl} = A_{ij}$ [25, 29], the number of independent components can be further reduced to 3, by:

$$\begin{align} \hat{A}_{11} + \hat{A}_{12} + \hat{A}_{13} &= \lambda_1, \nonumber \\ \hat{A}_{21} + \hat{A}_{22} + \hat{A}_{23} &= \lambda_2, \nonumber \\ \hat{A}_{31} + \hat{A}_{32} + \hat{A}_{33} &= \lambda_3, \tag{31} \end{align}$$

and

$$\begin{align} \hat{A}_{11} &= f_1(\lambda_1, \lambda_2), \nonumber \\ \hat{A}_{22} &= f_2(\lambda_1, \lambda_2), \nonumber \\ \hat{A}_{33} &= f_3(\lambda_1, \lambda_2), \tag{32} \end{align}$$

where $f_1$, $f_2$, and $f_3$ are functions of the independent eigenvalues of $\hat{A}$. For the ORE closure, these are defined as:
$$f_i(\lambda_1, \lambda_2) = C_{i1} + C_{i2}\lambda_1 + C_{i3}\lambda_2 + C_{i4}\lambda_1\lambda_2 + C_{i5}\lambda_1^2 + C_{i6}\lambda_2^2 \\ + C_{i7}\lambda_1^2\lambda_2 + C_{i8}\lambda_1\lambda_2^2 + C_{i9}\lambda_1^3 + C_{i10}\lambda_2^3 + C_{i11}\lambda_1^2\lambda_2^2 \quad i \in \{1, 2, 3\}. \quad (33) \\ + C_{i12}\lambda_1^3\lambda_2 + C_{i13}\lambda_1\lambda_2^3 + C_{i14}\lambda_1^4 + C_{i15}\lambda_2^4,$$

The coefficients $C_{i1}, \dots, C_{i15}$ are reported in Verweyst's PhD thesis [43], but are also available in [44].

Considering the normalization condition, the entries for $\hat{A}_{44}$, $\hat{A}_{55}$ and $\hat{A}_{66}$ can be computed as [43]:

$$ \begin{aligned} \hat{A}_{44} &= \frac{1}{2} (1 - 2\lambda_1 + \hat{A}_{11} - \hat{A}_{22} - \hat{A}_{33}), \\ \hat{A}_{55} &= \frac{1}{2} (1 - 2\lambda_2 - \hat{A}_{11} + \hat{A}_{22} - \hat{A}_{33}), \\ \hat{A}_{66} &= \frac{1}{2} (-1 + 2\lambda_1 + 2\lambda_2 - \hat{A}_{11} - \hat{A}_{22} + \hat{A}_{33}). \end{aligned} \quad (34) $$

Once the components of $\hat{A}$ have been calculated, the tensor is rotated to the case study reference frame by:

$$A_{ijkl} = R_{im}R_{jn}R_{ko}R_{lp}\hat{A}_{mnop}. \quad (35)$$

Another alternative is the IBOF closure approximation, proposed by Chung and Kown [41]. This closure does not require the eigendecomposition of $\underline{A}$ and, instead, resorts to the invariants of the second-order orientation tensor. The IBOF closure reads [41]:

$$A_{ijkl} = \beta_1 S(\delta_{ij}\delta_{kl}) + \beta_2 S(\delta_{ij}A_{kl}) + \beta_3 S(A_{ij}A_{kl}) + \\ \beta_4 S(\delta_{ij}A_{km}A_{ml}) + \beta_5 S(A_{ij}A_{km}A_{ml}) + \beta_6 S(A_{im}A_{mj}A_{kn}A_{nl}), \quad (36)$$

where the function S gives the fully symmetric part of its fourth-order tensor argument:

$$S(T_{ijkl}) = \frac{1}{24} [ T_{ijkl} + T_{ijlk} + T_{ikjl} + T_{iklj} + T_{iljk} + T_{ilkj} \\ + T_{jikl} + T_{jilk} + T_{jkil} + T_{jkli} + T_{jlik} + T_{jlki} \\ + T_{kijl} + T_{kilj} + T_{kjil} + T_{kjli} + T_{klij} + T_{klji} \\ + T_{lijk} + T_{likj} + T_{ljik} + T_{ljki} + T_{lkij} + T_{lkji}]. \quad (37)$$

The six $\beta_i$ coefficients are functions of the second ($I_2$) and third ($I_3$) invariants of $\underline{A}$, given by:

$$I_2 = \frac{1}{2} [\text{tr}(\underline{A})^2 - \text{tr}(\underline{A}^2)], \quad (38)$$

$$I_3 = \det(\underline{A}). \quad (39)$$

Additionally, only three coefficients are independent, $\beta_3$, $\beta_4$ and $\beta_6$, whose values are given by:

$$\beta_i = a_{i1} + a_{i2}I_2 + a_{i3}I_2^2 + a_{i4}I_3 + a_{i5}I_3^2 + a_{i6}I_2I_3 + a_{i7}I_2^2I_3 \\ + a_{i8}I_2I_3^2 + a_{i9}I_2^3 + a_{i10}I_3^3 + a_{i11}I_2^2I_3 + a_{i12}I_2^3I_3 \\ + a_{i13}I_2I_3^3 + a_{i14}I_2^4 + a_{i15}I_3^4 + a_{i16}I_2^4I_3 + a_{i17}I_2^3I_3^2 \\ + a_{i18}I_2^2I_3^3 + a_{i19}I_2I_3^4 + a_{i20}I_2^5 + a_{i21}I_3^5, \quad i \in \{3, 4, 6\}. \quad (40)$$

where the fitting coefficients $a_{i1}, \dots, a_{i21}$ are provided in the Appendix of [41]. The remaining coefficients are formulated as [41]:

$$\beta_1 = \frac{3}{5} \left[ -\frac{1}{7} + \frac{1}{5}\beta_3 \left( \frac{1}{7} + \frac{4}{7}I_2 + \frac{8}{3}I_3 \right) - \beta_4 \left( \frac{1}{5} - \frac{8}{15}I_2 - \frac{14}{15}I_3 \right) \right. \\ \left. - \beta_6 \left( \frac{1}{35} - \frac{24}{105}I_3 - \frac{4}{35}I_2 + \frac{16}{15}I_2I_3 + \frac{8}{35}I_2^2 \right) \right], \quad (41)$$

$$\beta_2 = \frac{6}{7} \left[ 1 - \frac{1}{5}\beta_3(1+4I_2) + \frac{7}{5}\beta_4\left(\frac{1}{6} - I_2\right) - \beta_6\left(-\frac{1}{5} + \frac{2}{3}I_3 + \frac{4}{5}I_2 - \frac{8}{5}I_2^2\right) \right], \quad (42)$$

$$\beta_5 = -\frac{4}{5}\beta_3 - \frac{7}{5}\beta_4 - \frac{6}{5}\beta_6\left(1 - \frac{4}{3}I_2\right). \quad (43)$$

Due to the reduced number of operations compared to eigen-fitted closures, the IBOF closure is considered as computationally more efficient, while maintaining a similar level of accuracy [41].
## 4. Modeling fiber orientation within the OpenFOAM framework

Tensor-based fiber orientation modeling is a pure advective problem, hence the formulations presented above are better suited for a Lagrangian approach [25]. To calculate the fiber orientation evolution equation in an Eulerian framework, the following definition of material derivative can be employed:

$$ \dot{\underline{A}} = \frac{\underline{D\underline{A}}}{Dt} = \frac{\partial \underline{A}}{\partial t} + \underline{v} \cdot \nabla \underline{A} = \frac{\partial \underline{A}}{\partial t} + \nabla \cdot (\underline{v} \underline{A}) - (\nabla \cdot \underline{v}) \underline{A}. \quad (44) $$

Due to the incompressibility assumption, $(\nabla \cdot \underline{v}) \underline{A}$ is null. The equation that governs the evolution of $\underline{A}$ is a component-wise passive scalar equation, coupled by a source term given by the right-hand side of the fiber orientation model. In order to introduce fiber orientation modeling into OpenFOAM®, some considerations have to be made. In Section 2, the velocity gradient was defined as $L_{ij} = \frac{\partial v_i}{\partial x_j}$, however in OpenFOAM®, $\nabla v_{ij} = L_{ji} = \frac{\partial v_j}{\partial x_i}$, which requires the adaptation of the vorticity tensor, $\underline{W}$, whose definition in the OpenFOAM® framework is the anti-symmetric of the one defined in Section 2 (Eqn. (3)).

Currently, OpenFOAM® does not provide a fourth-order tensor class, and only the *extend* branch offers a simplified fourth-order tensor class, with nine independent components and a limited set of operations. However, to be able to calculate the evolution of fiber orientation, it is necessary to know the product of $D_{kl}A_{ijkl}$. Therefore, $A_{ijkl}$ must be computed based on a closure relation, as described in Section 3. One potential solution for addressing this issue is to implement a comprehensive class to handle the operations. Alternatively, recognizing that $D_{kl}A_{ijkl}$ is a second-order tensor, it can be introduced as a source term in the evolution equation, and symbolic computation can be used to infer the required operations and generate the necessary code. In this work the latter approach was employed.

The framework for fiber orientation modeling was developed considering two abstract classes, **fiberOrientationModel** and **closureModel**. These were developed following the structure of the viscosity models available in `$FOAM_SRC/transportModels/incompressible/viscosityModels`, with run-time selection. For each class the required minimum attributes and member functions were added to model fiber orientation. Accordingly, the **fiberOrientationModel** class contains a dictionary for collecting input data from the user, fields that are common to all models ($\underline{A}$, $\underline{L}$, $\underline{D}$, $\underline{W}$, $\dot{\gamma}$, and $\underline{A} : \underline{D}$), a pure virtual **solve** method for assembling and solving the system of linear equations and some other minor methods to assure its functionality. The **closureModel** class includes one pure virtual function named **computeClosure** to compute the product $D_{kl}A_{ijkl}$ and one virtual function named **computeRSCClosure** to compute $[A_{ijkl} + (1 - \kappa)(L_{ijkl} - M_{ijmn}A_{mnkl})]D_{kl}$, required by the RSC model (Eqn. (14)), together with additional minor functionalities.

Based on the above structure, each fiber orientation model that derives from the **fiberOrientationModel** class should: (i) specify a sub-dictionary to collect model-specific information and any required additional fields (e.g. $\underline{C}$), and (ii) implement the corresponding **solve** method. Additionally, each closure model (hybrid, ORE, IBOF, etc.), which derives from the **closureModel** class, should implement the **computeClosure** method. If one decides to make the closure also available for the RSC approach, the method **computeRSCClosure** should also be defined.

The framework implemented in OpenFOAM® is depicted in Fig. 3. Currently, these classes are wrapped in a **functionObject (FO)** to work as a plug-in for any incompressible flow solver. At the end of each specified time step, upon calling the **execute** method from the **FO**, the **solve** function from a specific **fiberOrientationModel** will be executed, assembling and solving the corresponding system of linear equations. On the user's side, inside the **functions** section in the **controlDict** dictionary, in addition to the usual FO parameters, it is required to define a **dictionary** named **fiberOrientationProperties** that indicates the fiber model (parameter **model**), fiber parameters (subdictionary `<fiberModel>Properties`), and closure model (parameter **closureModel**).

Following the OpenFOAM® modeling approach, the final equation syntax is similar to the mathematical notation. As an example, the final form of the FT model, Eqn. (10), which is included in **solve()** method, is given in Listing 1. The variable `D_doubleDot_A4_` contains the calculation of $D_{kl}A_{ijkl}$, whose calculation is described in Section 5.

The resultant framework for the implemented fiber orientation modeling should be easily extendable to other fiber orientation and closure models. Currently, due to the validity of the implemented formulations, the available closure models are only valid for 3D cases.
<figure><img src="image_4.png" /><figcaption>Figure 3. Concept for the functionObject for modeling fiber orientation.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Listing 1. Folgar-Tucker fiber orientation model.</figcaption></figure>

5. Closure model with symbolic computation

The open-source library Sympy v1.11 [45] was used to perform the operations that OpenFOAM® does not have, namely the contractions between fourth- and second-order tensors. This rich-featured library for symbolic computation is written in Python language and provides functionalities to convert the symbolic code into C-style code, which can be used in OpenFOAM® [46].

In the developed Python scripts provided in the repository folder fiberClosureSymbolic, a vector is treated as a [3 × 1] matrix, a second-order tensor as a [3 × 3] matrix, and a fourth-order tensor as a dense multidimensional array with shape [3 × 3 × 3 × 3].

The Python script required to generate the C-code that calculates the closure method (computeClosure), for the hybrid closure model in the hybrid class, is illustrated in Listing 2. In Lines 1–3 of Listing 2, the Sympy module is imported with the alias sym and the functions symm, fourthOrderIndexPermutation and generateCCode from the module fiberOrientationModelling_Tools. These functions are responsible for generating a symmetric second-order tensor, performing index permutations within the fourth-order tensor and generating C-code, respectively.

Listing 2 comprises two functions, computeHybridClosure and main. The former computes the hybrid closure symbolically, while the latter, creates the symbolic variables, performs the tensor contraction and generates the associated C-code.

In the computeHybridClosure function, a matrix symbol is created representing a second-order tensor (Line 7), symmetry conditions are applied (Line 9), and a second-order identity tensor is defined (Line 11). Finally, in Lines 13–36, the hybrid closure is computed following Eqn. (28), to obtain $A_{ijkl}^{hyb}$.

In the main function, in Lines 41–43, a matrix representing the second-order tensor $\underline{D}$ is created. In Line 45 the computeHybridClosure function is executed, returning the computation of the hybrid
```python
import sympy as sym

from fiberOrientationModelling_Tools import symm, fourthOrderIndexPermutation,
    generateCCode

def computeHybridClosure():
    a = symalletsSymbol('A', 3, 3)
    A2 = symm(symetrixMatrix(a))
    I = sym眼(3)

    A4 linear = (
        symetrixRational(-1, 35) *
        fourthOrderIndexPermutation('ijkl', symetrisproduct(I, I))
        + fourthOrderIndexPermutation('ikjl', symetrisproduct(I, I))
        + fourthOrderIndexPermutation('iljk', symetrisproduct(I, I))
    )

    + symetrixRational(1, 7) *
    (  
        fourthOrderIndexPermutation('ijkl', symetrisproduct(A2, I))
        + fourthOrderIndexPermutation('ikjl', symetrisproduct(A2, I))
        + fourthOrderIndexPermutation('iljk', symetrisproduct(A2, I))
        + fourthOrderIndexPermutation('klij', symetrisproduct(A2, I))
        + fourthOrderIndexPermutation('jlik', symetrisproduct(A2, I))
        + fourthOrderIndexPermutation('jkil', symetrisproduct(A2, I))
    )

    A4 quadratic = symetrisproduct(A2, A2)
    f = 1.0 - 27*A2 det()
    A4 hybrid = f * A4 quadratic + (1.0 - f) * A4 linear

    return A4 hybrid

def main():
    d = symetrixMatrix('D', 3, 3)
    D = symm(symetrixMatrix(d))
    A4 = computeHybridClosure()
    D doubleDot A4 = symetriscontraction(symetrisproduct(D, A4), (0, 4), (1, 5))
    generateCCode(symetrixMatrix(D doubleDot A4))

if __name__ == "__main__":
    main()
```

<caption>**Listing 2.** Symbolic computation of the hybrid closure model.</caption>

closure, and, in Line 47, the contraction $D_{kl}A_{ijkl}$ is performed. To conclude this function, in Line 49, C-code is generated from the symbolic computations.

The derivation of the C-code to implement the `computeRSCClosure` method in the `Hybrid` class is illustrated in Listings 3 and 4. To compute the RSC closure, a script with a general framework was developed through the `RSC` function. This method takes a closure function as an input argument and performs the required symbolic computations. In Lines 2–8 of Listing 3, the slow-down symbolic variable,
$κ$, is created as well as variables that will hold the eigenvalues ($λ_i$) and eigenvectors ($e_i$) of $\underline{A}$. Although not required by the RSC approach, some orientation models and associated closures require the eigenvalues to be sorted in descending order and to have the eigenvectors as column vectors. Within the OpenFOAM® framework, the eigenvalues of a symmetric second-order tensor are computed and sorted in ascending order and the corresponding eigenvectors stored in a row-wise manner. The order of computations within OpenFOAM® was taken into account in the symbolic computations (see Lines 10–21).

In Line 23 the closure relation (the argument of the RSC function) is executed to obtain $A_{ijkl}$. Place holders for tensors $M$ (Line 25) and $L$ (Line 26) are created and the respective expressions are formulated (Lines 28–35) to get the expression between square brackets of the RSC model, Eqn. (14).

```c
1 def RSC(closureFunc):
2     k = sym Symbol('k')
3     e = symallet.SymMatrix('eigenVector', 3, 3)
5     eigVec = symngentMatrix(e)
6     w = sym.SymMatrix('eigenValue', 3, 1)
8     eigVal = sym.SymMatrix(w)
9
10 # Invert order to be consistent with OpenFOAM
11 tmpEigVal = eigVal.copy()
12
13 # sort eigenValues
14 eigVal[0] = tmpEigVal[2]
15 eigVal[2] = tmpEigVal[0]
16
17 # sort eigenVectors
18 eigVec = eigVec.T
19 tmpEigVec = eigVec.copy()
20 eigVec[:, 0] = tmpEigVec[:, 2]
21 eigVec[:, 2] = tmpEigVec[:, 0]
22
23 A4 = closureFunc()
24
25 M = sym.MutableDenseNDimArray(zeros(3,3,3,3))
26 L = sym.MutableDenseNDimArray(zeros(3,3,3,3))
27
28 for i in range(3):
29     evector = symming.SymMatrix(symritisflatteneigVec[i]))
30     M += symetensorproduct(evector, evector, evector, evector)
31     L += eigVal[i]*symetensorproduct(evector, evector, evector, evector)
32
33 tmp1 = symetensorcontraction(symetensorproduct(M, A4), (2,4), (3,5))
34
tensorCombination = A4 + (1-k)*(L- tmp1)
35
return tensorCombination
```

<caption>**Listing 3.** Symbolic computation of the RSC closure model.</caption>

The hybrid closure calculation used in the `computeRSCClosure` method is shown in Listing 4, which comprises only the `main` function. Similarly to the one shown in Listing 2, upon creating $A_{ij}$ (Lines 7–9), and $[A + (1 - \kappa)(L - M : A)]$ (Line 11), following the RSC framework, the contraction with the rate-of-deformation tensor is performed (Line 13). Finally, in Line 15, the associated C-code is generated.

## 6. Case studies

This section presents the analysis performed to verify the implementation of the fiber orientation plug-in, implemented for modeling the evolution of fiber orientation, which was undertaken with two case studies: simple-shear flow and the flow in a center-gated disk. For these, an adaptation of the solver `pimpleFoam` was used. In both cases, because an analytic solution for the velocity field is available, the `pimpleFoam` solver was modified to remove all computations of primitive variables. For assessment
```python
import sympy as sym
from fiberOrientationModelling symbolicComputationTools import symm, generateCCode
from hybricClosure import computeHybridClosure
from rsc import RSC

def main():
    d = sym.MatrixSymbol('D', 3, 3)

    D = symm( sym.Matrix(d) )

    A4 = RSC(computeHybridClosure)

    D_doubleDot_A4 = sym.tensorcontraction( sym.tensorproduct( D, A4 ), (0,4), (1,5) )

    generateCCode( sym.Matrix( D_doubleDot_A4 ) )

if __name__ == "__main__":
    main()
```

<caption>**Listing 4.** Symbolic computation of the hybrid closure model with RSC.</caption>

purposes, the fiber orientation evolution results obtained in OpenFOAM® were compared with reference data obtained from the Python scripts provided in the repository folder `fiberOrientationPython`, which integrate the same fiber orientation evolution equations assuming homogeneous time evolving flows. For that purpose, a Runge-Kutta5(4) time integration scheme [47] was used, with an absolute and relative tolerances of $10^{-12}$ for both quantities.

6.1. Simple-shear case study. The velocity field for a simple-shear flow, with velocity along $x$, and a null velocity at $y = 0$ plane, is given by:

$$ \underline{v} = \begin{bmatrix} \dot{\gamma}y \\ 0 \\ 0 \end{bmatrix}, \qquad (45) $$

which, in accordance with the definition of gradient in OpenFOAM®, will result in the velocity gradient:

$$ \nabla \underline{v} = \begin{bmatrix} 0 & 0 & 0 \\ \dot{\gamma} & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix} \qquad (46) $$

To verify the accuracy of the numerical implementation, a mesh refinement study using 7 progressively refined grids was performed on the three-dimensional geometry shown in Fig. 4. The coarser mesh has 32 cells along the $x$- and $y$-directions and one cell in the $z$- direction. The number of cells in the $x$-direction were doubled for each successively finer grid.

The initial velocity field and boundary conditions were defined with the utilities `setExprFields` and `setExprBoundaryFields` by applying the exact solution in Eqn. (45). The shear-rate and fiber shape factor were assumed to be unitary, $\dot{\gamma} = 1 s^{-1}$ and $\xi = 1$. The initial and Dirichlet boundary conditions for $A_{ij}$ on the `left` patch is given by:

$$ A_{ij} = \begin{bmatrix} 1/3 & 0 & 0 \\ 0 & 1/3 & 0 \\ 0 & 0 & 1/3 \end{bmatrix}. \qquad (47) $$

For the `"(right | top | bottom)"` patches a `zeroGradient` condition was assigned and for the `"(front | back)"` patches a `symmetryPlane` boundary condition was used.

For this case study, three fiber orientation models were tested, namely FT, iARD-RPR and MRD. The IBOF closure relation was employed with all models. For the FT model, two additional closure relations were tested, Hybrid and ORE. The relevant parameters for these models and closure relations are reported in Tab. 1 and were obtained from Favaloro and Tucker [27].

The discretization schemes and solution requirements employed are reported in Tabs. 2 and 3, respectively. The orientation evolution equation was integrated from 0 s until 65 s, with the time step size being controlled by the CFL condition for a maximum Courant number of 0.5.
<table><caption>**Table 1.** Fiber orientation model parameters for the simple-shear case study.</caption><thead><tr><th rowspan="2">Parameter</th><th colspan="3">Model</th></tr><tr><th>FT</th><th>iARD-RPR</th><th>MRD</th></tr></thead><tbody><tr><td>C<sub>I</sub></td><td>0.0311</td><td>0.0562</td><td>0.0198</td></tr><tr><td>C<sub>M</sub></td><td>—</td><td>0.9977</td><td>—</td></tr><tr><td>α</td><td>—</td><td>0</td><td>—</td></tr><tr><td>D<sub>i</sub></td><td>—</td><td>—</td><td>D<sub>2</sub> = 0.4796<br>D<sub>3</sub> = 0.0120</td></tr><tr><td>Closure</td><td>Hybrid<br>IBOF<br>ORE</td><td>IBOF</td><td>IBOF</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Geometry, mesh and patches for the simple-shear case study.</figcaption></figure>

<table><caption>**Table 2.** Temporal and spatial discretization schemes for the simple-shear case study.</caption><thead><tr><th>Term</th><th>Scheme</th></tr></thead><tbody><tr><td>ddtSchemes</td><td>Euler</td></tr><tr><td>gradSchemes</td><td>grad(U) Gauss linear</td></tr><tr><td>divSchemes</td><td>div(φ<sub>i</sub>,A<sub>2</sub>) Gauss upwind</td></tr><tr><td>laplacianSchemes</td><td>none</td></tr><tr><td>interpolationSchemes</td><td>linear</td></tr><tr><td>snGradSchemes</td><td>none</td></tr></tbody></table>

<table><caption>**Table 3.** Solution methods for the simple-shear case study.</caption><thead><tr><th>Parameter</th><th>p</th></tr></thead><tbody><tr><td><i>solver</i></td><td>PBiCGStab</td></tr><tr><td><i>smoother</i></td><td>DILU</td></tr><tr><td><i>tolerance</i></td><td>10<sup>-10</sup></td></tr><tr><td><i>relTol</i></td><td>0</td></tr></tbody></table>
6.2. Center-gated disk case study. The center-gated disk is a representative benchmark case of a typical geometry in the injection-molding (IM) process. Due to the usual small thickness of parts pro- duced by IM, the flow is typically dominated by shear, but extensional flow can also occur [48]. In this case study, assuming constant viscosity (Newtonian behavior) is not realistic for IM of thermoplastic materials. However, it was used for convenience in verifying the implemented fiber orientation add-on, as an analytical solution is available.

The geometry and patches considered for the present case study are illustrated in Fig. 5. Due to the axisymmetric flow conditions just a small angular region of the disk was employed for verification purposes. Accordingly, the geometry in this case study is a 5° slice from a typical center-gated disk with a feeding system, with a radius $r \in [0.01, 0.12]$ m and a thickness of 0.003 m.

<figure><img src="image_7.png" /><figcaption>Figure 5. Part and mesh representation with the corresponding patches for the center-gated disk case study (dimensions in m).</figcaption></figure>

The velocity field, in cylindrical coordinates ($r, \theta$ and $z$) for a Newtonian fluid in a center-gated disk is given by [41]:

$$ \underline{v} = \begin{bmatrix} v_r \\ v_\theta \\ v_z \end{bmatrix} = \begin{bmatrix} \frac{3Q}{8\pi rb} \left(1 - \frac{z^2}{b^2}\right) \\ 0 \\ 0 \end{bmatrix}, \qquad (48) $$

where:

* $Q$ is the volumetric flow-rate,

* $r$ is the radial position,

* $b$ is the half-thickness of the disk, and

* $z$ is the position along the thickness, where $z = 0$ corresponds to the mid-plane.

Taking into account the gradient definition in OpenFOAM®, the velocity gradient tensor is given by [41]:

$$ \nabla \underline{v} = \frac{3Q}{8\pi rb} \begin{bmatrix} -\frac{1}{r} \left(1 - \frac{z^2}{b^2}\right) & 0 & 0 \\ 0 & \frac{1}{r} \left(1 - \frac{z^2}{b^2}\right) & 0 \\ -\frac{2}{b} \left(\frac{z}{b}\right) & 0 & 0 \end{bmatrix}. \qquad (49) $$

Given the cylindrical symmetry of the problem and the Cartesian gradient implemented in OpenFOAM®, to ensure that Eqn. (49) was taken into account, the problem geometry was defined so that the $x$-direction aligns with the radial direction, $y = 0$, and $z$ is perpendicular to the flow direction. For verification purposes a mesh refinement study was carried out, using five progressively finer meshes. The coarser mesh has 128 cells along the radial and thickness directions and one along the angular direction ($\theta$ in the cylindrical coordinates). The cells in the radial direction were doubled for each successively finer grid to increase the accuracy of the fiber orientation evolution.

For this case study, the flux and velocity fields were computed assuming a total flow rate of $Q = 134.774 \text{ cm}^3\text{s}^{-1}$. Given that only a $\alpha = 5^\circ$ slice is being considered, the case study flow rate was adjusted to $Q\alpha/360$.

For this case, the flux and the gradient at the cell center were prescribed through a coded `functionObject` following Eqns. (48) and (49), respectively. For the second-order orientation tensor, the initial condition and Dirichlet boundary condition for the inlet were defined as isotropic state, Eqn. (47). The
**zeroGradient** boundary condition was applied to the walls and outlet, and **symmetry** to the symmWalls patch. The discretization and solution settings were the same as those considered for the simple-shear case study, Tabs. 2 and 3. Moreover, the fiber orientation parameters are described in Tab. 4.

A coded functionObject was used to define a steady-state criteria. The simulation was stopped when the initial maximum residual of $A_{ij}$ was below $10^{-6}$, and during 15 consecutive time steps the absolute difference between the current and previously computed $A_{ij}$ was below $10^{-6}$, in all computational cells.

<table><caption>**Table 4.** Fiber orientation model parameters for the center-gated disk case study.</caption><thead><tr><th></th><th>Model</th></tr></thead><tbody><tr><td>Parameter</td><td>FT</td></tr><tr><td>ξ</td><td>1</td></tr><tr><td>C<sub>I</sub></td><td>0.001</td></tr><tr><td>closure</td><td>IBOF</td></tr></tbody></table>

## 7. Results and discussion

### 7.1. Simple-shear case study.
In the simple-shear case study, taking into account the velocity field given by Eqn. (45), the $A_{xz}$ and $A_{yz}$ components of the fiber orientation tensor are null, whereas the $A_{zz}$ component will have a similar dynamic to the $A_{xx}$ component [39]. Thus, for verification purposes, only the $A_{xx}$ and $A_{xy}$ components will be analyzed. The $A_{xx}$ component indicates the alignment in the x-direction and the off-diagonal component the direction of orientation in the x-y plane [39]. To compare the results from the Python script with the ones computed in OpenFOAM®, and given that all streamlines are parallel to the x-axis, a pseudo-time ($t'$) was calculated based on the cell velocity and the domain length by:

$$t' = \frac{L_x}{u}, \qquad (50)$$

where $L_x$ is the length of the domain along the x-direction, with $L_x = 0$ at the inlet, and $u$ is the x-direction velocity component.

Given that only the velocity component in the x-direction is non-null, each line of cells in the x-direction is expected to match a section in the fiber orientation curve. Figures 6, 7 and 8 depict the comparison between the $A_{xx}$ and $A_{xy}$ components of the fiber orientation tensor with progressive mesh refinement for the FT model using the different closure models tested. The scatter data displayed corresponds to the line of cells at the y-position with the most significant numerical difference compared with the Python reference data. Figure 9 shows the evolution of $A_{xx}$ when using the iARD and MRD models with the IBOF closure.

The parameters defined for the models were set to obtain roughly the same value at a steady-state ($A_{xx} \approx 0.65$) [27]. Comparing the different closures used with the FT model it is clear that the hybrid model over-predicts the value of orientation for this simple-flow. This is a known issue in literature [39,41], however, due to its computational simplicity it is employed in some proprietary software. The IBOF and ORE present very similar transient results. This has also been addressed in literature [41], where the major practical difference is the computational effort. The ORE requires the eigendecomposition of the tensor $A_{ij}$ whilst the IBOF works with its invariants. Comparing the FT model with the iARD and MRD, although the solution at steady-state is the same, the transient behavior differs substantially. The FT model flattens out after $\approx 10$ s while the iARD and MRD model present a pronounced peak. This shows the ability of different phenomenological orientation models to capture different orientation behaviors, whose choice is not trivial.

The results also show that with progressive mesh refinement, the simulated data computed in OpenFoam® approximates the reference data generated with the Python script.

### 7.2. Center-gated disk case study.
The parabolic flow profile given by Eqn. (48) represents a non-homogeneous radial diverging flow field where the velocity magnitude decreases along the radial direction, and, accordingly, the shear ($\dot{\gamma}$) and stretch ($\dot{\epsilon}$) rates ratio are functions of the thickness and radial position, as illustrated in Fig. 10.

Figure 10 depicts the ratio between the shear and extensional components on half-thickness of the disk. For these conditions, null shear and elongation will occur at the mid-plane and wall, respectively.
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Mesh refinement study with the Folgar-Tucker model and the Hybrid closure for the simple-shear case study.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Mesh refinement study with the Folgar-Tucker model and the IBOF closure for the simple-shear case study.</figcaption></figure>

At the mid-plane, where $\frac{\partial u}{\partial x} < 0$ with null shear-rate, the fibers are compelled to orient along the cross-flow direction. As the magnitude of $z$ increases, the shear-rate becomes more pronounced, reaching its maximum at the wall. Consequently, the fibers tend to orient increasingly along the flow direction.

Similarly to the simple-shear case study, the $A_{xy}$ and $A_{yz}$ components of the fiber orientation tensor are always zero, whereas the $A_{zz}$ component will have a similar dynamic to the $A_{xx}$ (since the trace must be unitary). Thus, for verification purposes, only the $A_{xx}$ and $A_{xz}$ components will be presented.

Referring to the results shown in Figs. 10 and 11, near mid-plane, at $z/b = 0.102$, the $A_{xx}$ component initially tends towards zero due to the presence of a negative elongational rate and an insufficient shear-rate to keep the alignment of the fibers in the flow direction. However, as the fluid advances along the
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Mesh refinement study with the Folgar-Tucker model and the ORE closure for the simple-shear case study.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Mesh refinement study with the iARD and MRD models for the $A_{xx}$ component for the simple-shear case study.</figcaption></figure>

radial direction, the small but non-null shear-rate gradually, albeit slowly, aligns the fibers along the flow direction. This effect intensifies with $z/b$, i.e., while approaching the wall.

To compare the results obtained in OpenFOAM®, the FT model was integrated using 3000 linearly spaced points along the radius, and 128 points along the thickness corresponding the z-coordinate of the mesh. The orientation tensor was sampled at the cell centers with the **sampling function** Object utility of OpenFOAM®. The mesh refinement results for the closest ten evenly spaced $z/b$ ratios are displayed in Fig. 11. The mesh refinement study shows that with progressive refinement, the results obtained with OpenFOAM® approach the reference data generated from the *Python* script.
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Map of shear/stretch ratio for the center gated disk.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 11.** Mesh refinement study with the Folgar-Tucker model and the IBOF closure for the center-gated disk case study.</figcaption></figure>

## 8. Conclusions

In this work, a functionObject that works as a plug-in tool for any OpenFOAM® incompressible flow solver was programmed to calculate the evolution of fiber orientation in fiber reinforced thermoplastic materials, following a decoupled approach. Several fiber orientation models and closure relations available in the literature were implemented in the tool, which was verified by comparing the solver predictions with independent results obtained by numerically integrating the associated governing equations.

Due to OpenFOAM®'s limitations in handling higher-order tensors, the developed methodology for modeling fiber orientation, which requires fourth-order tensors, utilizes open-source symbolic computation to perform the necessary calculations. This symbolic approach retains only the non-trivial operations,
thus enhancing efficiency. However, this comes with the cost of code readability, as the functions resulting from the symbolic computations tend to be almost impossible to interpret.

State-of-the-art models and closure relationships currently employed in commercial solutions have been described and implemented in OpenFOAM®. Benchmark test cases, including a simple-shear flow and the flow in a center-gated disk, have been used to verify the numerical implementation, through the comparison with results obtained with a Python routine that integrates the equations for a single material point. A mesh refinement study was carried out for each type of flow. The results show that with progressive mesh refinement the OpenFOAM® results approach the Python counterpart.

The developed code, scripts, and case studies are open to the community, and the authors welcome suggestions and improvements.

## Acknowledgments

This work is funded by FEDER funds through the COMPETE 2020 Programme and National Funds through FCT (Portuguese Foundation for Science and Technology) under the PhD Grant 2020.07424.BD and projects UID-B/05256/2020, UID-P/05256/2020. The authors acknowledge M. Oliveira for the fruitful discussions related to the numerical implementation of the orientation models.

## References

[1] C.-T. Huang and H.-C. Tseng, "Simulation prediction of the fiber breakage history in regular and barrier structure screws in injection molding," *Polymer Engineering & Science*, vol. 58, no. 4, pp. 452-459, 2018. [Online]. Available: https://4spepublications.onlinelibrary.wiley.com/doi/abs/10.1002/pen.24660

[2] F. W. J. Van Hattum and C. A. Bernardo, "A model to predict the strength of short fiber composites," *Polymer Composites*, vol. 20, no. 4, pp. 524-533, 1999. [Online]. Available: https://4spepublications.onlinelibrary.wiley.com/doi/abs/10.1002/pc.10376

[3] M. Sepehr, G. Ausias, and P. Carreau, "Rheological properties of short fiber filled polypropylene in transient shear flow," *Journal of Non-Newtonian Fluid Mechanics*, vol. 123, no. 1, pp. 19-32, 2004. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0377025704001582

[4] S. Mortazavian and A. Fatemi, "Effects of fiber orientation and anisotropy on tensile strength and elastic modulus of short fiber reinforced polymer composites," *Composites Part B: Engineering*, vol. 72, pp. 116-129, 2015. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S1359836814005642

[5] D. Masato, J. Rathore, M. Sorgato, S. Carmignato, and G. Lucchetta, "Analysis of the shrinkage of injection-molded fiber-reinforced thin-wall parts," *Materials & Design*, vol. 132, pp. 496-504, 2017. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0264127517307001

[6] Y. Mu, G. Zhao, A. Chen, Y. Liu, Y. Song, and C. Sun, "Numerical investigation of three-dimensional fiber suspension flow by using finite volume method," *Polymer Bulletin*, vol. 74, no. 11, pp. 4393-4414, 2017.

[7] S. G. Advani and I. Tucker, *Charles L.*, "The Use of Tensors to Describe and Predict Fiber Orientation in Short Fiber Composites," *Journal of Rheology*, vol. 31, no. 8, pp. 751-784, 11 1987. [Online]. Available: https://doi.org/10.1122/1.549945

[8] J. Wang, J. F. O’Gara, and I. Tucker, *Charles L.*, "An objective model for slow orientation kinetics in concentrated fiber suspensions: Theory and rheological evidence," *Journal of Rheology*, vol. 52, no. 5, pp. 1179-1200, 09 2008. [Online]. Available: https://doi.org/10.1122/1.2946437

[9] J. H. Phelps and C. L. Tucker, "An anisotropic rotary diffusion model for fiber orientation in short- and long-fiber thermoplastics," *Journal of Non-Newtonian Fluid Mechanics*, vol. 156, no. 3, pp. 165-176, 2009. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0377025708001602

[10] H.-C. Tseng, R.-Y. Chang, and C.-H. Hsu, "An objective tensor to predict anisotropic fiber orientation in concentrated suspensions," *Journal of Rheology*, vol. 60, no. 2, pp. 215-224, 03 2016. [Online]. Available: https://doi.org/10.1122/1.4939098

[11] H.-C. Tseng, R.-Y. Chang, and C.-H. Hsu, "The use of principal spatial tensor to predict anisotropic fiber orientation in concentrated fiber suspensions," *Journal of Rheology*, vol. 62, no. 1, pp. 313-320, 01 2018. [Online]. Available: https://doi.org/10.1122/1.4998520

[12] A. Bakharev, H. Yu, S. Ray, R. Speight, and J. Wang, "Using new anisotropic rotational diffusion model to improve prediction of short fibers in thermoplastic injection molding," in *SPE ANTEC Conference*, 2018.

[13] S. K. Kugler, A. P. Dey, S. Saad, C. Cruz, A. Kech, and T. Osswald, "A flow-dependent fiber orientation model," *Journal of Composites Science*, vol. 4, no. 3, 2020. [Online]. Available: https://www.mdpi.com/2504-477X/4/3/96

[14] G. B. Jeffery, "The motion of ellipsoidal particles immersed in a viscous fluid," *Proc. R. Soc. Lond. A*, vol. 102, pp. 161-179, 1922.

[15] M. Gupta and K. K. Wang, "Fiber orientation and mechanical properties of short-fiber-reinforced injection-molded composites: Simulated and experimental results," *Polymer Composites*, vol. 14, no. 5, pp. 367-382, 1993. [Online]. Available: https://4spepublications.onlinelibrary.wiley.com/doi/abs/10.1002/pc.750140503

[16] W. Ogierman and G. Kokot, "A study on fiber orientation influence on the mechanical response of a short fiber composite structure," *Acta Mechanica*, vol. 227, no. 1, pp. 173-183, 2016.

[17] H.-C. Tseng, R.-Y. Chang, and C.-H. Hsu, "Numerical prediction of fiber orientation and mechanical performance for short/long glass and carbon fiber-reinforced composites," *Composites Science and Technology*, vol. 144, pp. 51-56, 2017. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0266353816318802

[18] K. Heinen, "Mikrostrukturelle orientierungszustände strömender polymerlösungen und fasersuspensionen," Ph.D. dissertation, Universität Dortmund, 2007.
[19] F. Ospald, "Contributions to the Simulation and Optimization of the Manufacturing Process and the Mechanical Properties of Short Fiber-Reinforced Plastic Parts," Ph.D. dissertation, Technische Universität Chemnitz, 2017.

[20] F. Wittemann, R. Maertens, A. Bernath, M. Hohberg, L. Kärger, and F. Henning, "Simulation of reinforced reactive injection molding with the finite volume method," *Journal of Composites Science*, vol. 2, no. 1, 2018. [Online].
Available: https://www.mdpi.com/2504-477X/2/1/5

[21] F. Wittemann, R. Maertens, L. Kärger, and F. Henning, "Injection molding simulation of short fiber reinforced thermosets with anisotropic and non-Newtonian flow behavior," *Composites Part A: Applied Science and Manufacturing*, vol. 124, no. May, p. 105476, 2019. [Online]. Available: https://doi.org/10.1016/j.compositesa.2019.105476

[22] K. Chiba and F. Chinesta, "Numerical simulation of flow kinematics and fiber orientation for multi-disperse suspension," *Rheologica Acta*, vol. 45, no. 1, pp. 1-13, 2005.

[23] G. Lipscomb, M. Denn, D. Hur, and D. Boger, "The flow of fiber suspensions in complex geometries," *Journal of Non-Newtonian Fluid Mechanics*, vol. 26, no. 3, pp. 297-325, 1988. [Online]. Available: https://www.sciencedirect.com/science/article/pii/0377025788800235

[24] J. Wang, "Improved fiber orientation predictions for injection molded composites," Ph.D. dissertation, University of Illinois at Urbana-Champaign, 2007.

[25] C. Tucker, *Fundamentals of Fiber Orientation: Description, Measurement and Prediction*. Hanser Publishers, 2022.

[26] F. Folgar and I. Charles L. Tucker, "Orientation behavior of fibers in concentrated suspensions," *Journal of Reinforced Plastics and Composites*, vol. 3, no. 2, pp. 98-119, 1984.

[27] A. J. Favaboro and C. L. Tucker, "Analysis of anisotropic rotary diffusion models for fiber orientation," *Composites Part A: Applied Science and Manufacturing*, vol. 126, p. 105605, 2019. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S1359835X19303549

[28] K. Breuer, M. Stommel, and W. Korte, "Analysis and evaluation of fiber orientation reconstruction methods," *Journal of Composites Science*, vol. 3, no. 3, 2019. [Online]. Available: https://www.mdpi.com/2504-477X/3/3/67

[29] S. K. Kugler, A. Kech, C. Cruz, and T. Osswald, "Fiber orientation predictions—a review of existing models," *Journal of Composites Science*, vol. 4, no. 2, 2020. [Online]. Available: https://www.mdpi.com/2504-477X/4/2/69

[30] A. J. Pontes, N. M. Neves, and A. S. Pouzada, "The role of the interaction coefficient in the prediction of the fiber orientation in planar injection moldings," *Polymer Composites*, vol. 24, no. 3, pp. 358-366, 2003. [Online]. Available: https://4spepublications.onlinelibrary.wiley.com/doi/abs/10.1002/pc.10035

[31] "Moldflow ard-rsc," https://help.autodesk.com/view/MFIA/2018/ENU/?guid=GUID-98E7B581-D2E6-4E9A-98C5-B79CAF919035, accessed: 2024-02-27.

[32] H.-C. Tseng, R.-Y. Chang, and C.-H. Hsu, "Phenomenological improvements to predictive models of fiber orientation in concentrated suspensions," *Journal of Rheology*, vol. 57, no. 6, pp. 1597-1631, 11 2013. [Online]. Available: https://doi.org/10.1122/1.4821038

[33] "Fiber function overview," https://support.moldex3d.com/2023/sync/for-nx/en/2-2-5-1_fiberfunctionoverview.html, accessed: 2024-02-27.

[34] S. K. Kugler, G. M. Lambert, C. Cruz, A. Kech, T. A. Osswald, and D. G. Baird, "Efficient parameter identification for macroscopic fiber orientation models with experimental data and a mechanistic fiber simulation," *AIP Conference Proceedings*, vol. 2205, no. 1, p. 020050, 01 2020. [Online]. Available: https://doi.org/10.1063/1.5142965

[35] A. Kech, S. Kugler, and T. Osswald, "Significance of model parameter variations in the pard-rsc model," *Journal of Composites Science*, vol. 4, no. 3, 2020. [Online]. Available: https://www.mdpi.com/2504-477X/4/3/109

[36] A. Scheuer, A. Ammar, E. Abisset-Chavanne, E. E. E. Cueto, A. Scheuer, A. Ammar, E. Abisset-Chavanne, E. E. E. Cueto, F. Chinesta, R. Keunings, and S. G. Advani, "Data-driven upscaling of orientation kinematics in suspensions of rigid fibres," *Computer Modeling in Engineering and Sciences*, vol. 117, no. 3, pp. 367-386, 2018.

[37] M. Yun, C. A. Martin, P. Giormini, F. Chinesta, and S. Advani, "Learning the macroscopic flow model of short fiber suspensions from fine-scale simulated data," *Entropy*, vol. 22, no. 1, p. 30, 2020.

[38] "Fiber orientation prediction theory references," https://knowledge.autodesk.com/support/moldflow-insight/learn-explore/caas/CloudHelp/cloudhelp/2014/ENU/MoldflowInsight/files/GUID-BD0CFAC8-D8A9-42E7-8838-7DD2E67829B1-htm.html, accessed: 2024-02-27.

[39] S. G. Advani and C. L. Tucker, "Closure approximations for three-dimensional structure tensors," *Journal of Rheology*, vol. 34, no. 3, pp. 367-386, 1990.

[40] J. Cintra, Joaquim S. and I. Tucker, Charles L., "Orthotropic closure approximations for flow-induced fiber orientation," *Journal of Rheology*, vol. 39, no. 6, pp. 1095-1122, 11 1995. [Online]. Available: https://doi.org/10.1122/1.550630

[41] D. H. Chung and T. H. Kwon, "Invariant-based optimal fitting closure approximation for the numerical prediction of flow-induced fiber orientation," *Journal of Rheology*, vol. 46, no. 1, pp. 169-194, 2002.

[42] D. Chung and T. Kwon, "Fiber orientation in the processing of polymer composites," *Korea-Australia Rheology Journal*, vol. 14, no. 4, pp. 175-188, 2002. [Online]. Available: http://ftp.cheric.org/research/tech/periodicals/view.php?seq=395367

[43] B. E. VerWeyst, "Numerical predictions of flow-induced fiber orientation in three-dimensional geometries," Ph.D. dissertation, University of Illinois at Urbana-Champaign, 1998.

[44] "Fiber orientation tools," https://github.com/charlestucker3/Fiber-Orientation-Tools, accessed: 2024-02-27.

[45] "Sympy," https://www.sympy.org/en/index.html, accessed: 2024-02-27.

[46] B. Ramoa, R. Costa, F. Chinesta, and J. M. Nóbrega, "A semi-automatic approach based on the method of manufactured solutions to assess the convergence order in openfoam," *OpenFOAM® Journal*, vol. 2, p. 148-165, Nov. 2022. [Online]. Available: https://journal.openfoam.com/index.php/ofj/article/view/75

[47] "Scipy user manual," https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html, accessed: 2024-02-27.

[48] R. Zheng, R. I. Tanner, and X.-J. Fan, *Injection Molding: Integration of Theory and Modeling Methods*. Springer-Verlag Berlin Heidelberg, 2011.