Florian Wittemann

**Band 101**

F. Wittemann

Methods for fiber-dependent injection molding simulation

**Fiber-dependent injection molding simulation of discontinuous reinforced polymers**

Florian Wittemann

### **Fiber-dependent injection molding simulation of discontinuous reinforced polymers**

### **Karlsruher Schriftenreihe Fahrzeugsystemtechnik Band 101**

Herausgeber

### **FAST Institut für Fahrzeugsystemtechnik** Prof. Dr. rer. nat. Frank Gauterin Prof. Dr.-Ing. Marcus Geimer Prof. Dr.-Ing. Peter Gratzfeld Prof. Dr.-Ing. Frank Henning

Das Institut für Fahrzeugsystemtechnik besteht aus den Institutsteilen Bahnsystemtechnik, Fahrzeugtechnik, Leichtbautechnologie und Mobile Arbeitsmaschinen.

Eine Übersicht aller bisher in dieser Schriftenreihe erschienenen Bände finden Sie am Ende des Buchs.

# **Fiber-dependent injection molding simulation of discontinuous reinforced polymers**

by Florian Wittemann

Karlsruher Institut für Technologie Institut für Fahrzeugsystemtechnik

Fiber-dependent injection molding simulation of discontinuous reinforced polymers

Zur Erlangung des akademischen Grades eines Doktors der Ingenieurwissenschaften von der KIT-Fakultät für Maschinenbau des Karlsruher Instituts für Technologie (KIT) genehmigte Dissertation

von Florian Wittemann, M.Sc.

Tag der mündlichen Prüfung: 3. Dezember 2021 Referent: Prof. Dr.-Ing. Frank Henning Korreferent: Prof. Prof. hon. Dr. Tim A. Osswald Dr.-Ing Luise Kärger

**Impressum**

Karlsruher Institut für Technologie (KIT) KIT Scientific Publishing Straße am Forum 2 D-76131 Karlsruhe

KIT Scientific Publishing is a registered trademark of Karlsruhe Institute of Technology. Reprint using the book cover is not allowed.

www.ksp.kit.edu

*This document – excluding parts marked otherwise, the cover, pictures and graphs – is licensed under a Creative Commons Attribution-Share Alike 4.0 International License (CC BY-SA 4.0): https://creativecommons.org/licenses/by-sa/4.0/deed.en*

*The cover page is licensed under a Creative Commons Attribution-No Derivatives 4.0 International License (CC BY-ND 4.0): https://creativecommons.org/licenses/by-nd/4.0/deed.en*

Print on Demand 2022 – Gedruckt auf FSC-zertifiziertem Papier

ISSN 1869-6058 ISBN 978-3-7315-1217-2 DOI 10.5445/KSP/1000148499

# **Kurzfassung**

Diese Arbeit präsentiert neuartige Simulationstechniken für Spritzgusssimulationen mit faserverstärkten Polymeren (FRPs).

Spritzguss ist einer der meistverbreiteten Prozesse zur Massenproduktion von diskontinuierlich faserverstärkten Polymerbauteilen. Die Prozessparameter (Füllrate, Temperatur, etc.) beeinflussen die Bauteileigenschaften signifikant. Für eine adäquate Vorhersage der finalen Bauteileigenschaften muss eine Simulation alle Prozessschritte (Formfüllung, Nachdruck, Abkühl- /Aushärtungsphase, Abkühlung außerhalb des Werkzeuges) beinhalten.

Während der Formfüllung hat die Strömungsmodellierung oberste Priorität. Das komplexe Matrixverhalten muss unter Beachtung von Scherrate, Temperatur und, falls vorhanden, chemischer Reaktion modelliert werden. Die sich ausprägende Faserorientierung, die von Strömungsfeld, Faserlänge und Volumengehalt abhängt, sollte aus zwei Gründen berechnet werden. Einer ist das Ausprägen von anisotropen Material- und somit auch Bauteileigenschaften aufgrund der Fasern. Zudem rufen die Fasern auch während der Formfüllung anisotropes Verhalten im flüssigen Material hervor. Auch die Faserlänge beeinflusst das mechanische und Fließverhalten des Materials und wird im Umkehrschluss durch das Strömungsfeld während der Formfüllung beeinflusst. Die Faserlänge hat großen Einfluss auf die Schlagzähigkeit des Bauteils, aber auch auf die effektive Viskosität in Faserrichtung im flüssigen Material. Umgekehrt erzeugt das Strömungsfeld aber auch Kräfte auf die Fasern, die diese zum Brechen bringen können. Stand der Technik Simulationen beachten den Einfluss der Faserorientierung und -länge auf das Strömungsfeld nicht. Diese Arbeit präsentiert einen neuartigen Ansatz, in welchem Viskosität, Faserorientierung, Faserlänge und Geschwindigkeit gekoppelt sind.

Zur Berücksichtigung der Fasereigenschaften in der Viskositätsmodellierung und somit auch in der Geschwindigkeit wird die Viskosität als Tensor vierter Stufe, der als Funktion von Matrixviskosität, Faserorientierung, -länge und volumengehalt definiert ist, modelliert. Der Viskositätstensor wird für eine homogenisierte Matrix-Faser-Suspension auf Basis von mikromechanischen Modellen berechnet. Für die Modellierung des Faserbruchs werden die hydrodynamischen Schlepp- und Auftriebskräfte beachtet. Zusätzlich werden makroskopische Ansätze zur Berechnung der Faser-Faser Interaktionskräfte (Schmier- und Reibkraft) gezeigt und verifiziert.

Neben der Formfüllung beeinflussen die weiteren Prozessschritte Nachdruck, Abkühl-/Aushärtungsphase und Abkühlung außerhalb des Werkzeuges ebenfalls die Bauteileigenschaften. Durch das anisotrop visko-elastische Verhalten können Verzug und Eigenspannungen aufkommen. Stand der Technik Software simuliert diese Phänomene in der Regel anisotrop mit linear elastischen Modellen. Diese Arbeit präsentiert einen Ansatz zur Berechnung von Verzug und Eigenspannungen für FRPs mit duromerer Matrix und thermo-visko-elastischen Modellen. Relevante Prozessdaten wie Faserorientierung, Temperatur und Aushärtungsgrad werden übertragen um diese in der Verzugssimulation mit zu betrachten. Faser- und Matrixeigenschaften werden zur Homogenisierung verwendet und unter Beachtung der Faserorientierung wird ein orthotropes Material definiert. Das Matrixverhalten wird als Funktion von Aushärtungsgrad und der Temperatur modelliert. Zusätzlich werden thermische und chemische Schwindung beachtet.

Die vorgestellten Methoden sind für Formfüllsimulationen in der opensource, finite Volumen basierten Software OpenFOAM und für die Verzugsanalyse in die kommerziellen finiten Elemente basierten Software Simulia Abaqus implementiert. Numerische Studien verifizieren die Implementierung und Methoden. Die Formfüllsimulationen zeigen eine gute Übereinstimmung mit experimentellen Ergebnissen, was die neu entwickelten Ansätze validiert.

# **Abstract**

This work presents novel simulation techniques for injection molding of fiber reinforced polymers (FRPs). Injection molding is one of the most applied processes for mass production of discontinuous FRP parts. The process conditions such as filling rate, temperature, etc. have a significant influence on the final part properties. For an adequate prediction of these properties, a process simulation has to depict different aspects, including all process steps, being mold filling, holding, in-mold solidification and out-of-mold cooling.

During the mold filling phase, the flow modeling is of major significance. The complex matrix behavior must be modeled under consideration of shearing, temperature and, if present, chemical reactions. The important aspect of fiber orientation, depending on flow field, fiber length and volume fraction, should be modeled for two reasons. The first one being, that the fiber orientation influences the anisotropic mechanical properties of the material, and therefore, the final part's behavior. Furthermore, the fibers also produce an anisotropic flow behavior in the liquid material during mold filling. Similar to the orientation, the fiber length influences the flow and mechanical behavior of the material and is vice versa influenced by the mold filling process. The fiber length is crucial for the mechanical impact strength and resistance and the effective viscosity in fiber direction. On the opposite the flow field evokes forces on the fibers, leading them to break. In state-ofthe-art simulation techniques, the influence of fiber orientation and fiber length on the flow field is not considered. Therefore, this work presents a novel simulation approach where viscosity, fiber orientation, fiber length and velocity field are coupled.

For consideration of the influence of fiber properties on viscosity and hence velocity, the viscosity is modeled with a fourth order anisotropic viscosity tensor, depending on matrix viscosity, fiber orientation, length and volume fraction. The viscosity tensor is calculated for a homogenized matrix fiber suspension, based on micro mechanical models. For modeling of fiber breakage, hydrodynamic drag and lift are computed. Furthermore, novel macroscopic approaches for fiber-fiber interaction forces (friction and lubrication) and contact points are presented and verified.

Besides the mold filling, the following process steps holding, solidification and out-of-mold cooling also have impact on the final part's properties and geometry. Due to the anisotropic and viscoelastic mechanical behavior, warpage may occur, or residual stresses build up. State-of-the-art software simulates these phenomena anisotropic with linear elastic mechanical models. This work presents an approach to calculate warpage and residual stresses for FRPs with thermoset matrix using thermo-chemo-elastic material models. Temperature and curing fields are mapped to be considered in the warpage simulation. Fiber and matrix properties as well as fiber orientations are used for homogenization to create an orthotropic material model. The matrix behavior depends on the degree of curing and temperature. Thermal and chemical shrinkage are also considered.

The presented methods are implemented in the open-source, finite volume based software OpenFOAM for mold filling simulations and in the commercial finite element based software Simulia Abaqus for warpage simulations. Numerical studies verify the implementations and methodology. The mold filling simulations are validated by comparison to experimental results, showing good agreement.

# **Acknowledgements**

At first, I want to thank Prof. Dr.-Ing. Frank Henning and Dr.-Ing. Luise Kärger for supervising my work and improving the quality of my work and my scientific skills with their support within the last four years, which enabled this thesis.

Secondly, I want to thank all my colleagues for their support on scientific and social level. Hereby, special thanks goes to Dr.-Ing. Nils Meyer, Julian Seuffert, Dr.-Ing. Martin Hohberg, Dr.-Ing. Alexander Bernath and Alexander Jackstadt for great discussions and input in case of flow and fiber orientation modeling as well as modeling of anisotropic visco-elastic material behavior. Additionally, I want to thank Robert Maertens for great support and expertise from the experimental side and provision of experimental data.

At last I want to thank the 'Vector Stiftung Stuttgart', the 'Ministerium für Wissenschaft, Forschung und Kunst Baden-Württemberg' and the 'Ministerium für Wirtschaft, Arbeit und Wohnungsbau Baden-Württemberg' for the financial support of research projects, making this thesis possible.



# **List of Symbols**

#### **Mathematical operations**









# **1 Introduction**

# **1.1 Motivation**

Reduction of costs, weight and CO2-emissions are major priorities in almost every engineering field. Especially in the mobility sector, systems need to be cost-effective and light. Due to this development, new challenges such as functional integration and low-density materials are rising.

The aspect of low-density drives discontinuous fiber reinforces polymers (FRPs) in focus of engineering tasks. These material systems offer good specific mechanical behavior due to fiber reinforcement in combination with low-density of the polymer matrix. Due to combination of fibers, matrix and several available additives, the material systems can be specified in regard of mechanical and thermal properties, chemical resistance, flame protection, etc. Furthermore, FRPs have a high potential for functional integration. The low temperature of the liquid material (compared to metals) creates the possibility of encapsulation of electric components, metal inserts and other functional groups during part manufacturing.

Manufacturing FRPs in combination with functional integration and mass production at low-cost leads to injection molding. Injection molding is one of the most frequently applied processes for manufacturing discontinuous FRPs with thermoplastic or thermoset matrix. Usually, the process is nonisothermal, transient and non-evacuated. Due to these aspects and since the process conditions have crucial impact on the final part properties, especially within the presence of fibers, a high-quality process simulation supports the prediction of the parts properties and geometry. Furthermore, a process simulation shows the general manufacturability of a given geometry under defined process conditions, considering pressure, temperature, air traps etc. Based on process simulations, process conditions can be optimized with respect to production time and energy effort. These aspects help to lower the development time and costs at an early stage of the engineering task.

Although several approaches are already existing to simulate the injection molding process and the flow behavior of FRPs has been in focus of several research projects, a state-of-the-art simulation still includes many simplifications and neglections. Therefore, this work presents a novel approach for all stages of the injection molding process, which takes the influence of fibers and in-mold air on the filling process into account and simulates the following process steps anisotropic and viscoelastic.

### **1.2 Objectives of the thesis**

The major objective of this thesis is increasing the level of detail for an injection molding simulation with FRPs for all process steps. Therefore, two points are addressed in mold filling simulations. One is the simulation of the in-mold air in addition to the FRP phase. The presence of air can lead to air entrapments for complex geometries and has influence on the in-mold pressure and thermal equilibrium. The second point is the influence of fibers on the flow behavior of the FRP and vise versa.

The co-existing FRP and air needs to be modeled during the mold filling simulation. Therefore, the simulation is extended from a single- to a multiphase approach. FRP and air are simulated as immiscible fluids with separate flow and thermal behavior.

By further considering the influence of fibers on the FRP flow behavior and enabling an anisotropic flow, the viscosity is extended from a scalar quantity to a tensor. The tensor is a function of fiber orientation and length. Since the fibers break during mold filling, the fiber length distribution is also considered during mold filling and fiber breakage is calculated, depending on flow induced forces.

For a higher degree of detail in simulations of the following process steps (after mold filling), an anisotropic, non-isothermal structural analysis is performed to predict shrinkage, warpage and residual stresses. Relevant mold filling results are transferred and the calculated fiber orientation is used to create homogenized and orthotropic materials.

### **1.3 Structure of the thesis**

#### **1.3.1 Structure of content**

This work is structured along the course of an injection molding process. In Section 2 the state of the art is given. Section 2.1 gives a short outline of the real process and its different phases and materials. Section 2.2 and Section 2.3 represent the state of the art for the corresponding simulation approaches. Section 1 contains the description of all novel approaches and methods, developed during the thesis. Within Section 1 the approaches are verified with numerical experiments (Section 4.1) and experimental process data (Section 4.2). Section 4.3 gives a short outlook on warpage analysis with the novel approaches. Section 5 includes the conclusion and outlook.

#### **1.3.2 Information about notation**

Within this work, index notation and the Einstein sum convention are used. Counting indices (for example coordinate directions) are visualized with Italian letters (,, , …). Specifications are represented by indices with normal letters (air, FRP, …). For scalars the specification indices are on the bottom (i.e. specification), since they may be potentiated. Specifications of vectors and tensors are represented by indices on the top (i.e. specification) to separate them from the counting indices. All tensors of fourth order in this work show right hand, left hand and main symmetry, so they can be visualized in Voigt notation. For Voigt notation the following order is chosen: [<sup>11</sup> <sup>22</sup> <sup>33</sup> <sup>23</sup> <sup>13</sup> 12]. Therefore, a fourth order tensor in Voigt notation is given by

$$T\_{ijkl} = \begin{bmatrix} T\_{1111} & T\_{1122} & T\_{1133} & T\_{1123} & T\_{1113} & T\_{1112} \\ & T\_{2222} & T\_{2233} & T\_{2223} & T\_{2213} & T\_{2212} \\ & & T\_{3333} & T\_{3323} & T\_{3313} & T\_{3312} \\ & & & T\_{2323} & T\_{2313} & T\_{2312} \\ & \text{sym.} & & & T\_{1313} & T\_{1312} \\ & & & & & T\_{1212} \end{bmatrix}$$

.

# **2 State of the art**

### **2.1 Injection molding process**

### **2.1.1 Machine and process overview**

Injection molding is a discontinuous process, well suited for mass production of complex shaped parts with high tolerances [1]. The major amount of research and industrial applications focuses on thermoplastic injection molding (TIM), while this work focuses on thermoset injection molding, named reactive injection molding (RIM), since thermosets solidify within an irreversible curing reaction (see Section 2.1.2). One manufacturing cycle, also named shot, is divided in plasticization, mold filling, holding, solidification (cooling or curing) and ejection. The cycle can be completely automated.

A general build-up of the injection molding machine and a tool is schematically illustrated in Figure 2.1. The machine is departed in a plastification and clamping unit. The material is feed in granular form through the hopper. Plasticization, transport and mixing is realized with a screw, the plasticization is supported by a heating system. During plastification, the screw retracts and material is accumulated under pressure in a chamber right before the nozzle (screw chamber). In order to inject the material into the mold, the screw moves forward to displace the material. The movement is controlled with hydraulics and usually a constant velocity or a velocity profile is defined. The material is injected via sprue (or runner system) into the mold which are part of the tempered tool. After mold filling the material solidifies under pressure and the part is ejected. During the end of solidification, the plastification unit can already melt material for the next shot to keep cycle times low.

The two major differences between of TIM and RIM, besides the material, are the process temperatures and material plastification. A thermoplastic is processed and injected at high temperatures (160 °C to 400 °C) into a cooled mold (10 °C to 160 °C), while a thermoset is plasticized and injected at cold temperatures (100 °C to 140 °C) into a heated mold (130 °C to 200 °C) [1]. Furthermore, thermoplastics have a more complex screw, containing special compression, mixing and shearing units, as well as a non-return valve. Thermosets are plasticized with a comparatively simple and compressionless screw, having no valve and a constant pitch between the flights [1,2]. The simple screw may lead to material or pressure losses during processing, but enables unscrewing if the material cures within the plastification unit.

Figure 2.1: Schematic illustration of an injection molding machine with tool

#### **2.1.2 Injection molding materials**

This section will only provide a small overview of injection molding materials. More detailed information about polymer materials, the molecular and atomic structure and the resulting properties as well as further applications and relevance in industry can be found in [1–5].

Injection molding can be performed with polymers of all three main categories (thermoplastics, thermosets, elastomers). The RIM process includes all reactive materials, which can be thermosets, elastomers and reactive thermoplastics. Since this work is focused on injection molding with thermosets, the usage of this material group should be considered, although all matrix material independent methods presented within this work are also applicable for thermoplastic materials. Even though the focus is on thermosets, some information about thermoplastics for injection molding should be provided, since TIM is the most important process variant for various industry sectors.

#### **2.1.2.1 Thermoplastics**

Besides the categories amorphous and semi-crystalline, thermoplastics are clustered in commodity (or standard), engineering (or technical) and high temperature (or high performance) thermoplastics, which also illustrates their typical field of usage. Commodity thermoplastics are mostly known from packaging sector due to their low price (< 2 €/kg [5]) and weak mechanical performance, which is a positive aspect in this case. Due to these properties the usage sector is mainly packaging, for example food and medical products. The most prominent thermoplastics in this group are PP, PET, PS and PVC [4].

As indicated by the name, engineering thermoplastics are mostly used in engineering applications. Compared to commodity thermoplastics, they show better mechanical properties, especially at higher temperatures (up to 150 °C) [5]. Due to the application in engineering tasks, the properties are often modified to fulfill special mechanical, thermal, electrical but also optical requirements. Therefore, additives and fillers like for example glass (or carbon) fibers are added to improve stiffness and strength. Typical applications are in the automotive and electronic sector, where these materials are used for support structures, housings, etc. The majority of parts is made of PA (6 and 66), PC and ABS [4].

High performance thermoplastics are characterized by good mechanical properties, especially at elevated temperatures. In order to improve the mechanical behavior, they may be fiber reinforced. Since costs are often a minor aspect for applications with high performance thermoplastics, the reinforcement is with carbon fibers, to unleash more lightweight potential in the most cases. Typical applications are in the aerospace sector and the most used materials are PEEK, PEK and PPS. One of the major material systems, illustrating the lightweight potential of this group, is PEEK reinforced with 30 %-wt. carbon fibers, available from various manufactures. Such a material has a tensile strength comparable to steel and aluminum alloys within a density of only about 1.4 g/cm³ and can be handled by less than 400 °C [6].

#### **2.1.2.2 Thermoset naterials**

Opposite to thermoplastics, thermosets solidify irreversibly and cannot be remolten. During curing (solidification) the polymer chains cross-link covalently, hence the atomic structure is represented by a network and not by individual chains, compared to thermoplastics [3]. Due to the covalent linkage, thermosets offer good mechanical properties and less tendency for creep compared to thermoplastics, only having van-der-Waals forces for chain coupling and physical cohesion.

Thermosets are mainly clustered in four groups: phenol formaldehyde or phenolic (PF), unsaturated polyester (UPE), epoxy (EP) and cross-linked polyurethanes (PU). Thermosets are typically highly filled, with fillers such as calciumcarbonat, wood (flour and fibers) or glass (spheres and fibers). Due to the high stiffness and strength PFs, UPEs and EPs are mostly used in engineering applications such as automotive parts. In the last few years, PFs have also been used for engine parts, caused by the good mechanical properties at high temperatures in combination with low cost production [7,8]. Epoxies show a low viscosity during manufacturing as well as good adhesion and dielectrical properties. Therefore, they are often used for electronic encapsulation and as base material for industrial glue. Typical PU applications are housings of computers, televisions and other commercial electronic products, due to the low cost and high impact strength of this thermoset group [1].

#### **2.1.2.3 Flow behavior of thermoset materials during processing**

The most important material aspect for mold filling behavior and mold filling simulations is the viscosity, defining the materials resistance to flowing. The viscosity of a polymer is influenced by temperature, shear rate and degree of curing/crystallization. Other aspects such as polymer chain architecture (length, side groups, etc.), fibers, pollutions and humidity also influence the material behavior, but are not explicitly mentioned for description of the viscosity.

A general course of polymer's viscosity for different temperatures and shear rates is shown in Figure 2.2. The viscosity is lower for higher temperatures, which is a typical behavior of liquids. Polymers often show a Newtonian behavior for low shear rates, meaning the viscosity is constant (see Section 2.2.2.3). For higher shear rates the materials act shear thinning and the viscosity decreases. The general viscosity level, presence and width of the Newtonian area as well as the change rate in the shear-thinning area depend on material and fillers [3].

Figure 2.2: Qualitive illustration of polymer viscosity over shear rate for two different temperatures

For processing of thermoset (or elastomer) polymers, the irreversible curing reaction needs to be considered for viscosity. Figure 2.3 illustrates the viscosity of a thermoset during processing. Since the material is heating up during mold filling (see 2.1.1) the viscosity would descent continuously during processing, but the cross-linking of the polymer chains has crucial impact on chain movement and hence rises the viscosity. The superimposition of these two effects leads to a parabolic related course of viscosity during processing, with a specific time window, where the viscosity is low and the material can be processed easily and energy efficient [5].

Figure 2.3: Qualitive illustration of thermoset viscosity during processing for temperature only (red), curing only (blue) and combined (green)

The complex viscosity behavior of thermoset materials leads to a complex flow behavior and mold filling. The flow behavior is shown schematically in Figure 2.4.

Figure 2.4: Fiber reinforced thermoset's flow behavior during mold filling. Hot surface layer (red) and cold core region (blue) in tool (grey) with corresponding fiber orientation, temperature and velocity profile

Due to the temperature difference between tool and material, a hot surface layer and cold core region build up. Since the temperature in the surface layer is high, viscosity of the resin is low. For the surface layer a shear dominated laminar flow is observed, while the flow in the core region is more plug-like, dominated by elongation stresses [9,10]. The overlay of these flow phenomena lead to a running ahead of the core region, resulting in a region with partial wall contact and an unstructured flow front [9–12]. Newer studies also show wall-slip effects, which influence the velocity profile, being not zero on the wall and hence also influencing viscous stress, fiber orientation, etc. [11]. This wall-slip effect is often neglected within simulations (see Section 2.2.1.5). Fibers align along the resulting forces of the flow field. Therefore, the alignment is in flow direction in the surface layer and perpendicular in the elongation dominated region [9]. Nevertheless, the fibers influence the viscosity and flow behavior (see Section 3.3). Fiber alignment in flow direction increases the viscosity in this direction and is therefore contrary to the effects mentioned above.

Some studies on the flow behavior of thermosets also show a qualitive distribution of the viscous stress over the different layers [9]. In these cases, the viscous stress is zero in the center of the core region (plane of symmetry) and increases linear towards the tool walls. The linear increase is only true for a Newtonian fluid in a perfect parabolic velocity field (Poiseuille flow), which are both not the case for RIM processes and should therefore not be assumed.

### **2.1.3 Mold filling**

Even though mold filling takes up only a small amount of time in a complete injection molding cycle, it is an important part defining various aspects of the final material state and part. The mold filling has crucial impact on fiber orientation, degree of curing, temperature distribution, air traps and energy requirement. Furthermore, insufficient mold filling may degrade the final part or damage the machine.

The material enters the cavity at the inlet and takes up the void volume until filled, as shown in Figure 2.5. By default, the cavity is not evacuated and the FRP displaces the in-cavity air, which leaves via venting slots, being too small for the liquid FRP to enter. In most cases, the mold filling is volume controlled (constant volume flow or profile), when the cavity is nearly filled, there is a switchover to a pressure-controlled filling. The switchover can be defined differently. Most common methods are position of the screw, hydraulic or screw chamber pressure. Another likely used criterion is reaching a defined mold pressure, dedicated by one or more sensors, ideally placed near the end of the flow path. Of course, these criteria can be combined and some aspect such as maximum screw position and maximum hydraulic pressure are generally present to secure the machine.

Figure 2.5: Schematic illustration of the mold filling process in injection molding. Material (green) enters at the inlet and fills the cavity (white)

The mold may contain functional parts such as metal inserts, threads or electric chips, sensors and connectors. These components get fixed or encapsulated during mold filling to functionalize the part. The filled cavity corresponds to the final geometry of the part, besides eventual deviations due to shrinkage and warpage. The following process steps after mold filling are described in Section 2.1.4.

### **2.1.4 Process steps following mold filling**

#### **2.1.4.1 Holding stage**

After the mold is filled, the next process step is the holding stage, when a high pressure (up to 200 MPa) acts from the plastification unit via the inlet on the material. The primary function of the holding pressure is to press out retaining air and counteract the shrinkage. In general, the shrinkage decreases with increase of holding pressure [13]. The amount of holding pressure depends on the material, part geometry and maximum clamping force. Similar to the filling rate, the holding pressure may be constant or assume a defined profile. The holding time depends on material, geometry and temperature. The holding stage is performed until the inlet is solidified, after this point, the pressure is not present in the cavity anymore. The mold should be designed, so sufficient holding stage is performable, depending especially on wall thickness and position of the inlet [1]. Amount and time can be estimated with simplified equations or determined and optimized with process simulations.

If the glass transition temperature (increasing with the degree of cure) reaches the material temperature, residual stresses may build up. Inhomogeneous temperature distribution and shrinkage as well as anisotropy due to fibers influence the evolution of residual stresses [14]. Residual stresses result in warpage after demolding. Several studies investigate these phenomena and their relations. Experimental works show the influence of holding pressure, material temperature and mold temperature [13,15]. The experimental investigations show, that warpage and shrinkage may be decreased by higher holding pressure and longer cooling or curing times, but this of course increases the energy requirements and process time. Therefore, an optimum must be found, which is also addressed by different investigations. Simulation models offer the possibility to perform a DoE with little effort, as for example shown by Huang and Tai [16] and Choi and Im [17]. Other approaches suggest relations with Kriging surrogate models [18] and artificial neural networks [19].

#### **2.1.4.2 In-mold solidification**

The solidification is no serial process step, since the material starts to solidify directly after thermal activation or mixing, including also the steps mold filling and holding. Nevertheless, a general process cycle contains a final solidification step, having the single aim to generate a form stable geometry and part. This step is necessary since the inlet is usually not the area with the greatest thickness, hence there is still liquid (or high viscous) material left in the cavity after the holding step. This step is named cooling for thermoplastics and curing for thermosets and elastomers. For RIM processes the curing step has a crucial impact on the final degree of curing and therefore on the final part's mechanical properties, shrinkage and residual stresses.

#### **2.1.4.3 Ejection and out-of-mold cooling**

After solidification the part is form stable and is ejected from the tool. The ejection is supported by multiple cylinders, pressing the part out. Of course, the mold must be designed the way, that the part can be ejected undamaged. Position, number and velocity of the ejection cylinders influence the final warpage [14].

After ejection the part cools down to room temperature and may post cure in this period. Since the part is no longer bounded to the tool geometry the residual stresses result in warpage. The amount of warpage depends on geometry, fiber orientation and volume fraction, temperature distribution, degree of cure and residual stresses, resulting from prior process steps. Additionally, the cooling in the atmosphere after ejection influences the warpage by temperature and thermal shrinkage. To reduce the effect of warpage or residual stresses, parts may be clamped and/or tempered after ejection.

### **2.2 Mold filling simulation**

#### **2.2.1 Process modeling**

#### **2.2.1.1 Governing equations**

In general, the injection molding process is a transient problem of fluid dynamic type. These can be solved analytically in some cases or under consideration of simplifications. Nowadays, fluid dynamic problems are solved in a computational fluid dynamics (CFD) simulations performed on a discretized mesh [1]. Solving a CFD problem for the mold filling simulations involves solving of the Navier-Stokes equations postulated in the first half of the 19th century [20], including the mass, momentum and thermal energy equation. The mass balance is shown in the form

$$\frac{d\rho}{dt} + \frac{\partial(\rho U\_l)}{\partial \mathbf{x}\_l} = \mathbf{0},\tag{2.1}$$

where is the density, the time, the coordinate vector and is the velocity vector. The momentum equation is re-presented by

$$\frac{d(\rho U\_l)}{dt} + \frac{\partial(\rho U\_l U\_l)}{\partial \mathbf{x}\_l} = -\frac{\partial p}{\partial \mathbf{x}\_l} + \frac{\partial \mathbf{r}\_{lj}}{\partial \mathbf{x}\_l},\tag{2.2}$$

with being the pressure and the viscous stress tensor. Body forces such as gravity would be formulated with terms like , but are often neglected, as in the present work, since they only show minor significance. The thermal energy equation is described with Fourier's Law by

$$\begin{split} \frac{d\left(\rho c\_{\rm p} T\right)}{dt} &+ \frac{\partial \left(\rho c\_{\rm p} U\_{l} T\right)}{\partial \mathbf{x}\_{l}} \\ = \lambda\_{\rm th} \frac{\partial T}{\partial \mathbf{x}\_{l} \partial \mathbf{x}\_{l}} + \mathbf{r}\_{lj} D\_{lj\prime} \end{split} \tag{2.3}$$

where is the temperature, th and <sup>p</sup> are thermal conductivity and specific

heat capacity. The term calculates the viscous heat/dissipation with the strain-rate tensor , being the symmetrical part of the velocity gradient [1,20]. Other terms such as thermal radiation and source terms due to chemical reactions are neglected within this work.

For an isotropic fluid (scalar viscosity) the viscous stress tensor is given by

$$\pi\_{lj} = 2\eta \left( D\_{lj} - \frac{1}{3} \frac{\partial U\_k}{\partial \varkappa\_k} \delta\_{lj} \right), \tag{2.4}$$

where is the dynamic viscosity and is the unity tensor. For anisotropic fluids, the viscosity is described by a tensor of higher order. Similar to anisotropic mechanical behavior, the viscous stress for an anisotropic case can be written as

$$
\pi\_{lj} = \eta\_{ljkl} D\_{kl}.\tag{2.5}
$$

The symmetry rules of continuum mechanics also apply for fourth order viscosity tensor , therefore it contains up to 21 independent entries. Nevertheless, at state of the art the anisotropy of FRPs is usually neglected and the viscous stress is calculated with an isotropic viscosity as described in Eq. (2.4).

#### **2.2.1.2 Analytical approaches**

Since fluid dynamic problems exist far longer than the possibility of solving numerical methods on computers, analytical approaches and solutions exist for special cases. Nevertheless, the solution of the Navier-Stokes equations in a 3D-case is part of the Millennium Problems, therefore all solutions come along with simplification, assumptions and restrictions.

One of the mentioned special cases is the so called Poiseuille flow, describing the laminar and stationary flow of a Newtonian fluid in an infinitly long circular tube with a pressure difference between the two ends and radius . The solution was developed by Hagen and Poiseuille in the mid of the 19th century [20]. By assuming the velocity to be zero at the walls and the pressure to be constant over the tube's cross section, the flow field is described by a parabolic velocity profile as shown in Figure 2.6. The velocity is given by

$$U\_1(r) = \frac{\Delta p}{4\eta l} (R^2 - \chi\_2^2),\tag{2.6}$$

where ∆ is the pressure drop within the distance .

Figure 2.6: Poiseuille flow velocity profile (green) in a circular tube

The equation can be modified for non-circular cross-sections [20]. Today the Poiseuille flow is still an often-used model for flow description and verification and therefore is used in various scientific publications.

Another milestone in flow field modeling is the description of a fluid in a narrow gap between two parallel plates, presented by Hele-Shaw at the end of the 19th century [1,20]. Considering an isotropic, incompressible fluid with no body forces and no wall slip the flow field is described by

$$\begin{aligned} \frac{\partial}{\partial \boldsymbol{x}\_{1}} \left( \boldsymbol{S}\_{\text{HS}} \frac{\partial p}{\partial \boldsymbol{x}\_{1}} \right) + \frac{\partial}{\partial \boldsymbol{x}\_{2}} \left( \boldsymbol{S}\_{\text{HS}} \frac{\partial p}{\partial \boldsymbol{x}\_{2}} \right) &= \mathbf{0} \\ \boldsymbol{S}\_{\text{HS}} = \int\_{0}^{\boldsymbol{\chi}\_{3,\text{max}}} \frac{\boldsymbol{\chi}\_{3}^{2}}{\eta} d\boldsymbol{x}\_{3} \end{aligned} \tag{2.7}$$

with the fluidity HS. Here, the narrow gap between the two plates is in the

range 0 < <sup>3</sup> < 3,max, the two plates are at <sup>3</sup> = 0, <sup>3</sup> = 3,max and perpendicular to the <sup>3</sup> -direction. <sup>3</sup> is zero in the complete domain.

It should be noted, that the viscosity is within the integral since it is not assumed to be constant and may depend on <sup>3</sup> , temperature or shear-rate. For a Newtonian fluid, the integral can be solved independent of the viscosity and the velocity is given by

$$\begin{aligned} U\_1 &= \frac{1}{2\eta} \frac{\partial p}{\partial \boldsymbol{\chi}\_1} \boldsymbol{\chi}\_3 (\boldsymbol{\chi}\_{3,\text{max}} - \boldsymbol{\chi}\_3) \\ U\_2 &= \frac{1}{2\eta} \frac{\partial p}{\partial \boldsymbol{\chi}\_2} \boldsymbol{\chi}\_3 (\boldsymbol{\chi}\_{3,\text{max}} - \boldsymbol{\chi}\_3) \end{aligned} \tag{2.8}$$

resulting in a parabolic velocity profile.

Due to the thin wall character of many injection molding parts, the solution of Hele-Shaw is well suited for injection molding simulation, even though it was developed long before the first injection molding processes. Therefore, many of today's 2D and 2.5D mold filling simulation rely on Hele-Shaw. Nevertheless, it should be handled carefully, since it is only valid for thin areas far away from the inlet and side walls. [1]

#### **2.2.1.3 Numerical methods**

Today it is standard to solve the governing equations for a CFD-simulation under consideration of numerical methods. The most prominent numerical methods are Finite Element Method (FEM), Finite Volume Method (FVM) and Finite Difference Method (FDM). The principle of these three methods is identical, the geometry is discretized into a finite number of elements as shown in Figure 2.7. The problem is then solved for the single elements (or cells) and combined to a global solution. The simulation mesh can be realized with different levels of complexity as illustrated in Figure 2.7. The simplest mesh is a 2D-midplane mesh, representing the geometry by one or more meshed planes. Of course, also a 1D-approximation is possible, but most injection molding simulation software perform at least a 2D-simulation. The 2.5D-mesh discretizes only the geometry's surfaces and is also known as

Boundary Element Method. Additionally, to the boundary surfaces, a finite number of 'inner planes' can be regarded as indicated by the blue elements in Figure 2.7.

The most complex variant is the 3D-mesh, representing the complete geometry with 3D-elements. While the FEM is solved on nodes at the corners and edges of the elements, the FVM focuses on the cell centers and surface fluxes. It must be distinguished between boundary nodes/faces, where boundary conditions are applied and internal nodes/cells, where algorithms solve the problem.

Figure 2.7: Original Geometry and different possibilities to discretize the geometry in a finite element mesh

Every numerical method has its individual advantages and disadvantages. While FEM is well suited for Langrage approaches, with the nodes being able to move, FVM is flux based and therefore well suited for Euler approaches and hence CFD simulations. However, most commercial injection molding software is based on FEM, due to the numerical efficiency, enabling faster solving [21–23]. Nevertheless, injection molding simulation is a CFD problem, which is most likely solved with FVM in other cases. Therefore, some scientific publications offer FVM-based solutions for injection molding simulations [24–26]. Hence, FVM-based simulation approaches are applied within this work.

Besides the three mentioned discretization-based methods CFD problems are also solved with Smooth Particle Hydrodynamics (SPH). This mesh-less approach splits the simulation domain into single particles considering neighbor particles in a defined area for solving. SPH has been successfully used to describe FRP flow processes on microscopic scale in [27].

#### **2.2.1.4 Influence of mesh size**

The mesh size has a crucial impact on the solution quality and numerical effort. One the one side a finer mesh, meaning smaller elements or more discrete points, predicts a more realistic solution, but on the other side this directly correlates to more numerical effort, since more equations must be solved. Nevertheless, deviations from the analytical solution may occur if the mesh is not fine enough. A typical case is the approximation of a parabolic velocity profile, as it may occur in injection molding (see Section 2.2.1.2). The analytical solution and approximations with increasing number of nodes are shown in Figure 2.8.

The parabolic distribution is discretized by three, five and nine nodes with linear interpolation. While the three node solution creates great deviations and is not able to display the analytical solution, the five node solution already shows a quite good agreement. The nine node approach shows only small deviations and displays the analytical solution with good agreement. Of course, the interpolation between the nodes could also be of higher degree. Therefore, the parabolic distribution can be displayed exact by three nodes with interpolation of second order, but the higher order comes along with greater numerical effort. In summary, the mesh fineness is always a compromise between level of detail and numerical effort.

Since the quality of result depends on the mesh, the development of novel methods should include a mesh study, showing a convergence of results for different mesh refinements. A convergence can be reached if the approximated distribution can be represented correctly (for example Figure 2.8 with quadratic interpolation), or if a finer mesh shows no more improvement for the results. A multiphase simulation with immiscible single phases (FRP and air in this case) always depends on the mesh, since the theoretical interface is reconstructed, depending on mesh dimensions [28]. Nevertheless, the flow front shows same tendencies and characteristics for adequatly fine meshes, although the results are not identical. This aspect is further investigated in [26] and Section 4.1.4.

Figure 2.8: Approximation of parabolic distribution with different number of points and linear interpolation between the points

#### **2.2.1.5 Simulation model**

The mold filling simulation is a transient simulation, where the FRP enters the simulation mesh at a defined node for FEM or a defined boundary surface (inlet) in FVM. While FEM simulations with commercial software are singlephase, considering only the FRP, the FVM-based simulations presented in [24–26] perform a multi-phase simulation, considering FRP and air.

The material in-flow is defined as a boundary condition and can be chosen similar to the real process, as described in Section 2.1.3. This also applies for the switchover point. The regarded geometry usually contains only the cavity, not the tool or plastification unit. Figure 2.9 shows different states of filling for a multiphase mold filling simulation of an electrical engine housing. Although there is a wall-slip in RIM processes, as described in Section 2.1.2.3, this effect is often neglected within a simulation, as it is within this work. A no-slip boundary condition is applied at the walls, justified by the need of less parameters and still creating good results [25,26].

Figure 2.9: Different states of a multiphase mold filling simulation of an electrical engine housing. Cavity filled with air in transparent grey and FRP in red

The FRP enters the model at the inlet on the top and the simulation is performed until the cavity is completely filled. One reason for performing this simulation is to verify if the part is manufacturable by the chosen process parameters with respect to maximum machine pressure, clamping force, air traps, material temperature, etc. Another important aspect is information about the material state needed for further analysis, including temperature and fiber orientation distribution as well as solidification state.

### **2.2.2 Material modeling**

#### **2.2.2.1 Thermal properties**

The relevant thermal properties are specific heat capacity and thermal conductivity (Eq. (2.3)). Additionally, the reaction enthalpy may be a material property, relevant for thermal modeling, but the influence of reaction heat is often neglected, as within this work [14,25,26]. This assumption is justified by the thin wall character of the molded parts, so the influence of the reaction heat is neglectable, compared to the thermal influence of the massive and tempered tool.

The specific heat capacity <sup>p</sup> is a scalar property. Independent of polymer type, <sup>p</sup> increases with temperature and shows a jump for the phase change [29,30]. Kamal and Ryan [29] present measurements of the specific heat depending on temperature and degree of cure showing also a significant increase of <sup>p</sup> during curing. Since only the liquid state is relevant for mold filling and the degree of cure has only marginal changes in this process phase, it is often assumed to be constant or only function of temperature.

The thermal conductivity th also increases with temperature and jumps at the phase change. Also similar to p, an increase of th is shown in [29]. For the same reasons compared to specific heat capacity, the thermal conductivity is often assumed to be constant or only function of temperature. The thermal conductivity may be a multi-dimensional tensor for anisotropic materials. The homogenized thermal conductivity can be calculated with the same assumption about serial and parallel connections as performed for mechanical properties [31,32]. Nevertheless, the anisotropy of th is also often neglected during form filling, since the values are generally low (< 0.5 W/(m∙K)) for most polymer composites and the mold filling is performed within a short period of time.

#### **2.2.2.2 Curings kinetics**

The curing reaction, meaning the cross-linking of molecular chains with covalent bonds, of thermoset materials has crucial impact on the viscosity and hence on the mold filling behavior, even though the fraction of crosslinked material is small during mold filling. Therefore, the curing kinetics should be modeled accurately for a mold filling simulation of a RIM process.

The curing of thermosets is an irreversible and exothermal process, which can be modeled by mechanistic or empirical models. Mechanistic models describe the chemical processes during cross-linking and therefore depend on detailed, material specific information, which are hard to determine [14]. Empirical models focus on simpler description with less parameters. Here, the reaction process is described by a simple differential equation, depending on a few empirical parameters, degree of cure and in some cases temperature. The simplest possibility to describe the reaction is given by

$$\frac{dc}{dt} \sim (1 - c)^n,\tag{2.9}$$

with being the degree of cure and the reaction order in this case. One disadvantage of this approach is, that the maximum change of cure is always at = 0, which is no typical behavior of industrially used thermoset materials [14].

Kamal and Sourour [33] extended this approach by

$$\frac{dc}{dt} \sim c^m (1 - c)^n,\tag{2.10}$$

with an empirical parameter in this case, so the change rate is zero at the beginning of the reaction. Finally, the curing kinetics are calculated with the so called Kamal-Malkin (or also Kamal-Sourour) kinetics model [33], determining the rate of change as

$$\frac{dc}{dt} = (K\_1 + K\_2 c^{m\_{\text{KM}}}) (1 - c)^{n\_{\text{KM}}},\tag{2.11}$$

with the empirical parameters KM and KM. The parameters <sup>1</sup> and <sup>2</sup> are often given by an Arrhenius type approximation in the form

$$K\_{1,2} = A\_{\text{KM1},2} \cdot \exp\left(\frac{-E\_{\text{KM1},2}}{R\_{\text{g}} \cdot T}\right) \tag{2.12}$$

with empirical factors KM1 and KM2, activation energies KM1 and KM2 and universal gas constant <sup>g</sup> [14].

The literature offers more complex models for better descriptions, if the temperature is for example near the glass transition temperature [14,34]. Nevertheless, the Kamal-Malkin model is often used and well established for curing kinetics model for different materials and processes [23,25,26,34–36]. Within this work the Kamal-Malkin kinetics model (Eq. (2.11)) in combination with Eq. (2.12) is used to model curing kinetics during the complete process simulation.

#### **2.2.2.3 Viscosity properties**

The viscosity must be modeled to determine the viscous stress tensor, which is needed to solve the momentum equation as shown in Eq. (2.2). The simplest approach to model the viscosity is to assume the FRP as Newtonian, resulting in a constant viscosity

$$
\eta = \text{const.}\tag{2.13}
$$

Meyer et al. [37] and Sommer et al. [38] present good results for simulating compression molding with Newtonian matrix behavior. Nevertheless, polymers show a clearly non-Newtonian behavior, with a viscosity depending on temperature and shear rate, described in2.1.2.3.

The simplest description of a shear rate depending viscosity in case of simple shear load is given by the power-law approach [39], determining the shear stress by

$$
\pi = K\_{\rm PL} ' \cdot \dot{\mathcal{Y}}^{n\_{\rm PL}},\tag{2.14}
$$

and with = ̇ the dynamic viscosity by

$$
\eta = K\_{\rm PL} \cdot \dot{\nu}^{(n\_{\rm PL}-1)} \text{.} \tag{2.15}
$$

with ̇ being the scalar shear rate, defined as ̇ = √ 2 . PL and PL are constant, positive model parameters in this case. The power-law approach models shear thinning ( < 1), shear thickening ( > 1) and Newtonian behavior ( = 1). One disadvantage, is the prognostication of a zero-shear stress and zero viscosity for zero-shear rate, which may be true for the shear stress, but not for viscosity. To overcome this, the approach can be extended to describe a so called Herschel-Bulkley fluid [40] with the viscosity given by

$$\eta = \begin{cases} \tau\_0 \dot{\mathbb{y}}\_0^{-1} + K \cdot \dot{\mathbb{y}}\_0^{n-1} & \dot{\mathbb{y}} < \dot{\mathbb{y}}\_0 \\ \tau\_0 \dot{\mathbb{y}}^{-1} + K \cdot \dot{\mathbb{y}}^{n-1} & \dot{\mathbb{y}} \ge \dot{\mathbb{y}}\_0 \end{cases} \tag{2.16}$$

where ̇<sup>0</sup> and <sup>0</sup> are material specific parameters. For < <sup>0</sup> the material behaves as a solid. ̇<sup>0</sup> indicates the shear rate for a change from Newtonian to shear thinning behavior (cf. Figure 2.2). Today's most common approaches for viscosity modeling of polymer melts are of Cross type. Cross [41] describes the viscosity by

$$\eta = \frac{\eta\_0}{1 + \left(\frac{\eta\_0 \dot{\nu}}{\tau^\*}\right)^{1 - n\_{\rm CWLF}}},\tag{2.17}$$

with zero-shear viscosity <sup>0</sup> and ∗ indicating the critical shear stress for beginning of shear thinning behavior. CWLF is a material specific parameter. In order to respect the temperature dependence, <sup>0</sup> is often described by the Williams–Landel–Ferry (WLF) equation like

$$\begin{aligned} \eta\_0 &= D\_{\text{CWLF1}} \exp\left(-\frac{A\_{\text{CWLF1}}(T - T\_{\text{g}})}{A\_{\text{CWLF2}} + (T - T\_{\text{g}})}\right) \\ T\_{\text{g}} &= D\_{\text{CWLF2}} + D\_{\text{CWLF3}} p \\ A\_{\text{CWLF2}} &= A\_{\text{CWLF3}} + D\_{\text{CWLF3}} p \end{aligned} \tag{2.18}$$

where <sup>g</sup> is the glass transition temperature and CWLF1 , CWLF3 ,CWLF1 ,CWLF2 and CWLF3 are material specific fitting parameters. For most materials, it is assumed to be CWLF3 = 0, so <sup>g</sup> and CWLF2 are constant. The combination of Eq. (2.17) and Eq. (2.18) is the so called Cross-WLF viscosity model, being the most used viscosity model for injection molding simulations with thermoplastic materials.

Although the Cross-WLF model is well suited for thermoplastic materials, there are better approaches for thermoset materials. Castro and Macosko [42] reported a better fitting by describing the temperature dependence with an Arrhenius equation with an additionally term to take the curing kinetics into account. For thermoset material the today's most common viscosity model is the so-called Castro-Macosko (CM) model, combining a Cross type equation with Arrhenius equation and curing dependence by

$$\eta = \frac{\eta\_0}{1 + \left(\frac{\eta\_0 \dot{\mathbf{y}}}{\tau^\*}\right)^{1 - n\_{\rm CM}}} \left(\frac{c\_{\rm g}}{c\_{\rm g} - c}\right)^{(C\_{\rm CM} + C\_{\rm CM} c)},\tag{2.19}$$

$$\eta\_0 = B\_{\rm CM} \exp\left(\frac{T\_{\rm CM}}{T}\right)$$

with material specific fitting parameters CM, CM1, CM2, CM, CM and <sup>g</sup> being the gelation conversation point.

#### **2.2.2.4 PvT modeling**

The aim of pressure-volume-temperature (PvT) modeling is the relation between pressure, specific volume (or density) and temperature, also often titled as equation of state. The simplest way to describe this relation is to assume incompressibility, so

$$
\rho = \text{const.}\tag{2.20}
$$

or equivalent ̇<sup>11</sup> + ̇<sup>22</sup> + ̇<sup>33</sup> = 0, with ̇ being the strain rate tensor. This assumption is valid for polymer injection molding and often chosen in different studies [24,37,38,43].

Since this study includes the air within the process simulation, this phase would also be considered incompressible, which is an unsatisfying simplification. The PvT behavior of gases is often given by the perfect gas law so

$$
\rho = \frac{p}{TR\_s},
\tag{2.21}
$$

with <sup>s</sup> as specific gas constant. The best related description of liquids in this case is to assume a perfect fluid, determining the density with

$$
\rho = \frac{p}{TR\_{\rm s}} + \rho\_0. \tag{2.22}
$$

Here, <sup>0</sup> is the density at = 0 .

During the history of polymer processing several PvT models with different complexity and experimental effort have been created. An overview of different models in relation to different polymers is given by Rodgers [44] and Júnior et al. [45]. Both studies name the so-called Tait model (or Tait equation) as most used and best fitting model for a variety of polymers and pressure/temperature scales. The Tait equation represents the specific volume by

$$\mathbf{u} = \boldsymbol{\nu}\_0(T) \left[ 1 - \mathcal{C}\_t \ln \left( 1 + \frac{\boldsymbol{p}}{B(T)} \right) \right] + \boldsymbol{\nu}\_t(\boldsymbol{p}, T). \tag{2.23}$$

t is an universal constant with value 0.0894. The Tait model is also called two-domain Tait model, since the formulation changes, depending on the volumetric transition temperature, defined as

$$T\_{\mathbf{t}}(p) = b\_5 + b\_6 p. \tag{2.24}$$

The subfunctions, needed for Eq. (2.23) are given in Table 2.1, wherein = − <sup>5</sup> . All mentioned -parameters within this Section are purely empirical and material-specific fitting parameters.

Function <sup>t</sup> () < <sup>t</sup> () > 0 () 1m + 2m 1s +2s () 3mexp(−4m ) 3sexp(−4s ) (, ) 0 7exp(8 − 9)

Table 2.1: Definition of subfunctions for Tait model (Eq. (2.23))

Wang et al. [46] further improved the model, by describing the -parameters as function of the temperature change rate and not as constants.

Due to a lack of available PvT data, the FRPs within this work are described nearly incompressible. Since multi-phase approaches are performed, a complete incompressible simulation would be unsatisfying with respect to the air phase. Therefore, the air will be assumed as perfect gas (Eq. (2.21)) and the FRP is assumed to be a perfect fluid (Eq. (2.22)), where <sup>s</sup> is chosen very high, so /(s) ≈ 0 and ≈ <sup>0</sup> .

#### **2.2.3 Fiber orientation modeling**

#### **2.2.3.1 Movement of a single fiber**

The modeling of fiber orientation has been in the focus of research for several decades. Today's models are based on Jeffery's work to describe the orientation change of a single ellipsoid in a Newtonian fluid with the so called Jeffery's equation [47]. The orientation of a single fiber is described by the nominated orientation vector , as shown Figure 2.10. The orientation can be formulated in Cartesian coordinates or with the orientation angles and , as shown in Figure 2.10. The Jeffery's equation for is

$$\frac{dq\_l}{dt} = \mathcal{W}\_{lj} q\_j + \frac{r\_{\rm f}^2 - 1}{r\_{\rm f}^2 + 1} \Big( D\_{lj} q\_j - q\_l \{q\_l D\_{lj} q\_j\} \Big),\tag{2.25}$$

with the vorticity tensor , being the unsymmetrical part of the velocity gradient and fiber aspect ratio <sup>f</sup> , defined as

$$r\_{\rm f} = L\_{\rm f} / d\_{\rm f} \tag{2.26}$$

where <sup>f</sup> and <sup>f</sup> are fiber length and diameter.

Figure 2.10: Orientation of a single fiber (blue) represented by vector or angles and

Fiber matrix suspensions are clustered in three regimes: dilute, semi-dilute and concentrated. The regimes are defined by fiber volume fraction <sup>f</sup> or <sup>f</sup> , with being the number density of fibers in the suspension. In dilute suspension it is <sup>f</sup> ≪ 1 and <sup>f</sup> ≪ 1 [48] or <sup>f</sup> ≪ 1/<sup>f</sup> 2 [49]. The space between the fibers is large and fibers rotate without influence by their neighbors, Eq. (2.25) is valid. The stationary solution in this case is a periodic rotation of the fiber. In semi-dilute suspensions (<sup>f</sup> ≪ 1 and <sup>f</sup> ≫ 1 [48] or 1/<sup>f</sup> <sup>2</sup> ≪ <sup>f</sup> ≪ 1/<sup>f</sup> for isotropic orientation [49]) fiber-fiber interactions occur. If also <sup>f</sup> is no longer small, a concentrated regime is reached, fiber interactions significantly influence the movement and fibers are no longer

free to move independently. At state of the art, the movement of fibers in semi-dilute and concentrated regimes cannot be described analytically.

Folgar and Tucker [50] developed a semi-empirical model to take fiber interactions into account. The Folgar-Tucker-model is an extension of the Jeffery's equation with an empirical interaction coefficient <sup>I</sup> . The stationary solution is a final orientation, depending on <sup>I</sup> . The approach still describes the evolution of a single fiber over time. Applying this model, meaning calculate the movement of every single fiber within a mold filling simulation, creates a huge numerical effort, which is unacceptable at state of the art.

#### **2.2.3.2 Macroscopic orientation modeling**

Advani and Tucker [51] present a homogenization scheme, describing the fiber orientation evolution with orientation tensors to reduce the numerical effort. The exact formulations of the second and fourth order tensors and are

$$A\_{lj} = \int \psi(\mathbf{q}) q\_l q\_j \mathbf{dq} \tag{2.27}$$

and

$$A\_{ljkl} = \int \psi(\mathbf{q}) q\_l q\_j q\_k q\_l \mathbf{d}\mathbf{q},\tag{2.28}$$

where ( ) (or (,)) is the probability density function for a specific orientation. Inclusion of into the Folgar-Tucker model leads to

$$\begin{split} \frac{dA\_{lj}}{dt} + \frac{\partial U\_l A\_{lk}}{\partial \boldsymbol{x}\_k} &= -\left(\mathcal{W}\_{lk} A\_{kj} - A\_{lk} \mathcal{W}\_{kj}\right) \\ + \frac{r\_l^2 - 1}{r\_l^2 + 1} \Big\{ \mathcal{D}\_{lk} A\_{kj} + A\_{lk} \mathcal{D}\_{kj} - 2 \mathcal{D}\_{kl} A\_{kllj} \Big\} \\ &+ 2 \mathcal{C}\_l \dot{\boldsymbol{\gamma}} \{ \mathcal{S}\_{lj} - 3 \mathcal{A}\_{lj} \}. \end{split} \tag{2.29}$$

This approach builds the base for several actual orientation models with different focuses [52–56].

Today's most commonly used models are the reduced strain closure (RSC) model by Wang et. al. [52] for slow orientation due to high fiber volume fractions and the RSC-ARD-model (anisotropic rotary diffusion) by Phelps and Tucker [53] for additionally long fibers. The models become more complex to consider more effects by need of more empirical parameters. To keep the empirical parameters low and since small aspect ratios (<sup>f</sup> ≤ 100) being in the focus of this work, the RSC model is used. Based on the Folgar-Tucker model and by considering the eigenvectors and eigenvalues of the orientation tensor, the approach describes the evolution of the second order orientation tensor of by

$$\begin{split} \frac{dA\_{lj}}{dt} + \frac{\partial U\_l A\_{jk}}{\partial \boldsymbol{x}\_k} &= -\{\mathcal{W}\_{lk} A\_{kj} - A\_{lk} \mathcal{W}\_{kj}\} \\ + \frac{r\_{\rm f}^2 - 1}{r\_{\rm f}^2 + 1} \{D\_{lk} A\_{kj} + A\_{lk} D\_{kj} \\ -2[A\_{ljkl} + (1 - \kappa) \\ \{L\_{ijkl} - M\_{ljmn} A\_{mnkl}\}] \boldsymbol{D}\_{kl}\} \\ + 2\kappa \boldsymbol{C}\_{\rm f} \dot{\boldsymbol{y}} \{\mathcal{S}\_{lj} - 3 \boldsymbol{A}\_{lj}\}. \end{split} \tag{2.30}$$

The fourth order tensors and are calculated with the eigenvectors and eigenvalues of as described in [52]. The so-called strain reduction factor is empirical and material specific, with ∈ [0,1]. For = 1 the RSC-model reduces to the Folgar-Tucker-model (Eq. (2.29)).

#### **2.2.3.3 Characteristics of**

Since the orientation tensor is a complex construct, information to interpret the results is required. The side entries reflect the deviation to the coordinate axis, hence it is ∈ [−0.5,0.5] for ≠ . A fiber has no front or back end, therefore, = − is valid and ∈ . The major entries represent the amount of fibers orientated in the respective direction. Due to normalization it is tr() = 1, which leads (in combination with ∈ sym.) to five independent entries [48].

One simple and often used interpretation is to see the eigenvectors as individual fibers with an orientation probability of its corresponding eigenvalue. Therefore, transversely isotropic materials can be created by only regarding the direction with the highest eigenvalue, or the tensor is used for orientation averaging. Further information is given in Section 2.3.2.2.

Figure 2.11 shows specific fiber orientations with corresponding orientation tensor. The three cases quasi-isotropic, aligned and 2D random (Figure 2.11a-c) are often used orientations to verify novel approaches.

Figure 2.11: Different fiber orientations with corresponding orientation tensor. Quasiisotropic/3D random (a), unidirectional/aligned (b), planar isotropic/2D random (c) and fabric like/2D aligned (d)

Comparing Figure 2.11c and Figure 2.11d illuminates one disadvantage of the orientation tensor, since it is identical in both cases, but the orientations are clearly different and the mechanical behavior of the final parts would show significant deviations. Similar configurations can also occur for 3D orientations.

### **2.2.3.4 Determination of**

The interaction coefficient represents fiber-fiber interactions, thus it depends on <sup>f</sup> and fiber orientation [57]. In past works, <sup>I</sup> is determined with experiments [58,59] or mathematically as function of <sup>f</sup> and <sup>f</sup> [48,60]. Heinen compared the experimental and numerical approaches in her PhD thesis [61]. The results show, that a combination of [48,60] such as

$$\mathcal{L}\_{\rm I} = \begin{cases} 0.03(1 - \exp(-0.224\phi\_{\rm f}r\_{\rm f})) & \Phi\_{\rm f}r\_{\rm f} \le 1.3\\ 0.0184 \exp(-0.7148\phi\_{\rm f}r\_{\rm f}) & \Phi\_{\rm f}r\_{\rm f} > 1.3 \end{cases} \tag{2.31}$$

fits best to the experimental work in [58] and [59]. Within this work <sup>I</sup> is calculated according to this approach, which is similar to [61].

#### **2.2.3.5 Closure approximations for**

The change rate of requires the unknown fourth order orientation tensor . Similar to the second order tensor, a change rate equation can be determined for the fourth order tensor, which would create numerical effort and require an unknown orientation tensor of sixth order [51]. Therefore, is determined with a closure approximation.

The first solutions for such closure approximations are the linear approach, which is exact for an isotropic distribution, the quadratic approach, which is exact for aligned fibers, and the hybrid approach, which is a combination of linear and quadratic [62]. Today there are several approaches for closure approximations to determine based on invariants or eigenvalues of [63]. Within this work the IBOF5 closure approximation developed by Chung and Kwon [64] is used to determine , since was found to be a good approximation with quite little numerical effort.

### **2.2.4 Fiber breakage modeling**

#### **2.2.4.1 Modeling of decreasing fiber length**

During processing, or during compounding, fibers may get damaged by interaction with the wall or other fibers as well as by hydrodynamic load of the surrounding matrix. Damage in this case is similar to fiber breakage, hence the average fiber length within a compound or part is at least as long as in the initial material or shorter. Since the fiber length has crucial impact on the mechanical behavior of an FRP, an adequate prediction of the final fiber length distribution is crucial for subsequent material modelling.

One of the first approaches to model fiber length distributions during compounding is presented by Shon et al. [65], describing the length evolution by

$$\frac{dL\_{\rm f}}{dt} = -k\_{\rm Shannon}(L\_{\rm f} - L\_{\rm fo}),\tag{2.32}$$

with kinetic constant Shon, determined for different process devices and the equilibrium fiber length f∞ representing the shortest possible fiber length. Inceoglu et al. [66] extended this approach by considering the specific mechanical energy Πme so

$$\frac{dL\_{\rm f}}{dt} = -K^{\prime\prime} \cdot \Pi\_{\rm me} \cdot \{L\_{\rm f} - L\_{\rm f\alpha}\}.\tag{2.33}$$

Here ´´ represents a change rate constant, determined by data of a batch mixer.

#### **2.2.4.2 Fiber breaking mechanisms**

The approaches of Shon and Inceoglu may describe the shortening of fibers, but in an empirical and linear way. The model parameters may be well suited for compounding processes, but the variety of process conditions and geometries in injection molding is too high to describe the complex process of fiber breaking with such simple models.

In a first step it is crucial to identify the mechanisms, which lead the fibers to break. In principal these are identical to other damage mechanisms, meaning a specific strain or stress is reached. However, the material behavior of fibers and matrix during processing is complex and it is ambiguous which load case or cases are present, so it is unclear if the fiber is breaking due to bending, stretching or buckling. Several studies identify buckling as the major, but not singular, effect for fiber breakage [47,67–70]. Of course, effects such as interaction with walls or other fibers as well as imperfection in material and fiber geometry may favor or complicate this process.

Assuming buckling as breaking criterium, it is necessary to verify that the stress field of the fiber, induced by the flow field, puts the fiber in compression. Therefore, the orientation between the fiber and the flow field must be known. The orientation state, where the fiber is under compression, is known as Jeffery orbit [47] and defined by

$$D\_{lj}q\_iq\_j < 0.\tag{2.34}$$

Figure 2.12 shows the values of for all possible orientations, represented by a unity sphere in case of simple and normed shear flow. According to the color bar the Jeffery orbit is given by the blue area.

After the load case is known, the acting forces and stress distribution must be determined, to evaluate if the fiber is buckling. A first approach is given by Burgers [71], assuming the forces to act from the fiber end to the center along the fiber axis. In a later work Hernandez et al. [68] extended this approach, determining analytical and pseudo-analytical solutions of the force integral along the axis, showing a non-linear distribution. This aspect is quite important for determining the buckling point, meaning the critical force above the fiber buckles. Durin et al. [69] performed simulations for fiber loads with different orientations, flow fields and aspect ratios, also based on the model of Burgers. Within the study it is shown, that for most cases the buckling criterion is reached before the breakage criterium (critical stress). The conclusion is, that for brittle fibers, such as glass or carbon, buckling then will

lead to breakage and is acceptable and meaningful as singular breakage criterium.

Figure 2.12: Visualization of Jeffery orbit on a unity sphere in simple shear flow. Colors represent value of ⟦ ⟧ , blue regions (negative values) represent Jeffery orbit

Defining buckling as the mechanism for fiber breakage leads to the need of defining a point at which the fiber buckles. Therefore, the Euler buckling modes, derived by Euler in the middle of the 18th century are a common model. The critical load of this purely elastic deformation mode can by defined by stress, so

$$
\sigma\_{\rm bu} = E\_{\rm f} \left( \frac{\pi}{4r\_{\rm f}} \right)^2,\tag{2.35}
$$

with <sup>f</sup> being the elastic modulus of the fiber, used by Durin et al. [69]. The criterion can also be defined for forces as

$$F\_{\rm bu} = \frac{\pi^3 E\_{\rm f} d\_{\rm f}^4}{64 L\_{\rm f}^2},\tag{2.36}$$

used by Phelps et al. [70]. Both works identify the first buckling mode as critical, leading to the assumption, that fibers break only at one point, so only two new fibers are created. The work of Meyer et al. [27] shows, that at least in dilute suspensions also other eigenmodes can be reached and fibers may break in three parts.

#### **2.2.4.3 Fiber length distribution models**

To predict an adequate fiber length distribution within the final part, a physically based breakage model considering material properties and process condition must be defined. It can be seen in literature, that a crucial part of fiber breakage takes place within the plastification unit in case of injection molding. Nevertheless, this part of the process is often quite similar and the conditions change only slightly, so a known state of fiber length distribution at the end of plastification can be used as input for the mold filling simulation. Two actual models for fiber length distribution modeling with breakage mechanisms are given by Durin et al. [69] and Phelps et al. [70]. While the first one is built to model breakage in a twin screw extruder and uses single fiber orientations, the latter one is focused on mold filling in macroscopic injection molding simulations. Therefore, only the approach of Phelps et al. will be further discussed.

In a first step, some assumptions must be defined. The most important are, that there is a minimum fiber length, where no more fiber breakage is possible and all possible fiber lengths are a multiple of this length. Furthermore, the flow field is assumed to be pure shear [70]. Additionally, the force acting on the fiber must be determined, to compute if the fiber buckles. Therefore, Phelps et al. use the slender-body analysis presented by Dinh and Amstrong [72], representing the acting force by

$$F\_{\rm f} = \frac{\zeta \eta\_{\rm M} L\_{\rm f}^2}{8} \left( -D\_{lj} q\_l q\_j \right), \tag{2.37}$$

where is a dimensionless drag coefficient and <sup>M</sup> is the matrix viscosity. Buckling is assumed to happen if the critical force, defined in Eq. (2.36) is reached, so

$$\begin{array}{c} F\_n^{\rm f} = \frac{4\zeta \eta\_{\mathcal{M}} \dot{\mathcal{Y}} \Big(L\_n^{\rm f} \Big)^4}{\pi^3 E\_\mathcal{f} d\_\mathcal{f}^4} \Big( -2 \| D\_{l/} \| q\_l q\_\mathcal{l} \Big) \\\ = B\_n^{\rm bu} \Big( -2 \| D\_{l/} \| q\_l q\_\mathcal{l} \Big) \\\ \ge 1, \end{array} \tag{2.38}$$

with bu being the dimensionless buckling index and ⟦⟧ is the normed strain rate defined as ⟦⟧ = /̇. The index corresponds to the possible fiber length f . According to Figure 2.12 it is −1 ≤ −2⟦⟧ ≤ 1 for all possible orientations in case of pure shear. Hence there is no state satisfying Eq. (2.38) if bu < 1. Furthermore, Phelps et al. assume −2⟦⟧ = 1 in regions where the fibers are under compression. Regarding Figure 2.12, this means that the areas with x<sup>3</sup> > 0 and x<sup>1</sup> < 0 or x<sup>3</sup> < 0 and x<sup>1</sup> > 0 are fully and homogenous blue. Within this assumption, the remaining buckling criterium is bu ≥ 1 and no fiber orientations must be determined.

Based on numerical experiments Phelps et al. [70] define a breakage probability br for every possible fiber length as

$$P\_n^{\rm br} = \begin{cases} 0 & \text{for } B\_n^{\rm bu} < 1\\ \mathcal{C}\_{\rm br} \dot{\mathcal{Y}} \left(1 - \exp\{1 - B\_n^{\rm bu}\} \right) & \text{for } B\_n^{\rm bu} \ge 1' \end{cases} \tag{2.39}$$

with the breakage coefficient br. Within [70] it is clearly mentioned, that (1 − exp(1 − bu)) is a fitting approach, and better functions may be found, although the numerical experiments show a low sensitivity towards the exact formulation of this function.

In a constitutive breakage model, fibers must create new, shorter fibers when breaking. Therefore, Phelps et al. [70] introduce a child generation rate defined as

$$R\_{nm}^{\rm br} = \varrho\_m^{\rm br} \text{PDF}\left(L\_{n\nu}^{\rm f} \frac{L\_m^{\rm f}}{Z}, \mathcal{S}\_{\rm br}\right) \tag{2.40}$$

being a Gaussian probability density function (PDF) with standard deviation

br and mean value f /2, so fibers break most likely in the middle. The scaling factor br is defined to fulfill

$$\sum\_{n} R\_{nm}^{\text{br}} = 2P\_m^{\text{br}},\tag{2.41}$$

since every break creates two children fibers. Furthermore, it is

$$R\_{nm}^{\text{br}} = R\_{(m-n)m\prime}^{\text{br}} \tag{2.42}$$

because the length of one child and the parent fiber defines the length of the other child fiber. For example, if one child of a 1.0 mm fiber is 0.7 mm, the other one must be 0.3 mm and hence have the same probability. The last restriction to the child generation rate is

$$R\_{nm}^{\text{br}} = 0 \quad \text{for } n \ge m,\tag{2.43}$$

otherwise recombining fibers to create longer ones would be possible, which is unphysical.

Finally, the evolution equation of fiber lengths is given by

$$\frac{\partial N\_n^{\text{f}}}{\partial t} + U\_l \frac{\partial N\_n^{\text{f}}}{\partial \boldsymbol{\omega}\_l} = -P\_n^{\text{br}} N\_n^{\text{f}} + \sum\_m R\_{nm}^{\text{br}} N\_{m\nu}^{\text{f}} \tag{2.44}$$

where f is the number of fibers with length <sup>f</sup> within the control volume.

In summary, the approach of Phelps et al. [70] is able to capture the evolution of an arbitrary number of different fiber lengths within an injection molding simulation, considering process conditions. Furthermore, only a small amount of material information is needed. Physical parameters such as, matrix viscosity, elastic modulus and diameter of the fiber are usually known. In addition, the three empirical parameters , br and br must be determined. Although this approach is well suited to determine the fiber length distribution with respect to process parameters, there are disadvantages. One is the simplification −2⟦⟧ = 1 and another one is the need of an empirical parameter for force calculation. These two aspects will be addressed with novel approaches in Section 3.4 within this work.

### **2.3 Simulation of process steps following mold filling**

The simulation of following process steps is indispensable to predict shrinkage, residual stress and hence the final part geometry. The range of complexity and numerical effort of such simulations form a wide spectrum. The material model can be isotropic, linear elastic up to anisotropic, timedependent and chemo-thermal viscoelastic. The simulation model itself may be only the filled cavity or the complete tool with integrated cooling/heating system. Several research works investigate in modeling with linear elastic [73], viscoelastic [74] and thermo-viscoelastic approaches [17,75–78].

In general, the material flow during the holding pressure phase is neglected and the simulations are pure solid mechanics. Besides the stress-strain equilibrium, the thermal equation (Eq. (2.3)) still needs to be solved, although the viscous heating can be neglected due to the low strain rates, compared to the mold filling phase. Furthermore, evolution equations for chemical reactions and phase change are still considered, since the majority of such processes take part within these process steps and they have significant influence on the material's mechanical and thermal behavior.

### **2.3.1 Matrix material modeling**

The mechanical behavior of the matrix is an overlay of material aspects such as the exothermal curing reaction, glass transission temperature, material history, but also molecular aspects such as polymer-chain geometry. However, the latter are not considered in macroscopic simulations. Besides material behavior, also the temperature field is significant, since it co-defines the material temperature distribution. Therefore, a wide range of models with different complexity and attendant experimental effort for parameter determination exist. A good overview of polymer viscoelasticity is given by Bergström [79] and in the special case of epoxies with fiber reinforcement by Bernath [14].

#### **2.3.1.1 Viscoelastic modeling**

A viscoelastic material model is able to describe the material behavior within a wide range of frequencies. In addition, polymer phenomena such as relaxing or creep can be modeled. In a general form the stress tensor is given by

$$
\sigma\_{lj} = \int\_0^t G\_{ljkl}(t-\xi) \frac{d\varepsilon\_{kl}(\xi)}{d\xi} d\xi,\tag{2.45}
$$

with strain tensor , time integration variable and relaxation modulus tensor , which is also often represented by a scalar, since the matrix behaves isotropic. Descriptions for in case of injection molding warpage modeling can be found in [75,77] or by Choi and Im [17]. Unfortunately, all these studies focus on thermoplastic materials. According to the author's knowledge during writing this thesis, no work addressing warpage simulation for reactive injection molding with visco-elastic phenolics has been published so far. Although, Mirabedini et al. [80] investigate in phenolic/rubber compounds, the viscoelastic behavior is dominated by the rubber part, and the phenolic is already completely cured within the studies. Nevertheless, viscoelastic modeling of continuous fiber reinforced epoxies is investigated in several studies and the general modeling aspects are transferable, since epoxies show a more comparable behavior to phenolics, than thermoplastics to phenolics. For thermoset materials, the relaxation modulus is often approximated with a Prony series [81]

$$G\_{ljkl}(\mathbf{t}, T, \mathbf{c}) = G\_{ljkl}^{\mathbf{r}}(\mathbf{c}) + \sum\_{p=1}^{N} G\_{ljklp}(\mathbf{c}) \exp\left(\frac{\mathbf{t}}{\Theta\_p(T, \mathbf{c})}\right) \tag{2.46}$$

with r being the elastic modulus in the rubbery state and Θ as actual relaxation time, depending on temperature and state of cure. The suitability of Prony series to model thermoset viscoelasticity is shown in [82–84]. Since the relaxation time strongly depends on the material state [14,83], they are often approximated with a shift function such as

$$
\Theta\_p(T, \mathfrak{c}) = a\_\text{s}(T, \mathfrak{c}) + \Theta\_p^{\text{ref}}(T\_{\text{ref}}, \mathfrak{c}\_{\text{ref}}), \tag{2.47}
$$

where Θ ref are relaxation times referring to a reference state based on a master curve determined with experimental measurements. The shift-factor <sup>s</sup> can be modeled WLF based or with an Arrhenius-type function [14].

#### **2.3.1.2 CHILE model**

Although viscoelastic approaches represent the material adequately, a high experimental effort is necessary to characterize the material and feed the model. Therefore, simplified approaches with accompanying restrictions exist. One example is the class of the cure hardening instantaneously linear elastic (CHILE) models. In the line of simplifications between fully viscoelastic and CHILE, also path-dependent models exist, but since they are not relevant within this work they are not further discussed. More information about path-dependent models is given by Svenberg and Holmberg [85] as well as Bernath [14].

Compared to viscoelastic approaches, the material history is neglected by CHILE models, so

$$
\sigma\_{lj} = \int\_0^t E\_{ljkl}^{\mathbf{c}}(T, \mathbf{c}) \frac{d\varepsilon\_{kl}(\xi)}{d\xi} d\xi,\tag{2.48}
$$

with elastic stiffness tensor defined as

$$E\_{ljkl}^{\mathbf{c}}(T,\mathbf{c}) = \begin{cases} \begin{array}{c} E\_{ljkl}^{\mathbf{c}\mathbf{r}} \\ E\_{ljkl}^{\mathbf{c}\mathbf{r}} \end{array} & \Delta T\_1 \\\end{cases}$$

$$E\_{ljkl}^{\mathbf{c}\mathbf{r}} + \frac{T\_{\mathbf{g}}(\mathbf{c}) - T - T\_{\mathbf{c1}}}{T\_{\mathbf{c2}} - T\_{\mathbf{c1}}} \left( E\_{ljkl}^{\mathbf{c}\mathbf{g}} - E\_{ljkl}^{\mathbf{c}\mathbf{r}} \right), & \Delta T\_2. \\\ E\_{ljkl}^{\mathbf{c}\mathbf{g}} & \Delta T\_3 \\\end{cases} \tag{2.49}$$

Here it is ∆<sup>1</sup> ∈ <sup>g</sup> () − < c1, ∆<sup>3</sup> ∈ <sup>g</sup> ()− > c2 and ∆<sup>1</sup> ≤ ∆<sup>2</sup> ≤ ∆<sup>3</sup> . c1 and c2 are material specific parameters. cr , cg and c are the instantaneous elastic moduli in pure rubbery, pure glassy and actual material state.

One disadvantage of the CHILE approach is the non-release of frozen-in stress, after reheating above <sup>g</sup> + c1, although there is a change of material stiffness. For normal process conditions, this aspect is uncritical, since the material is heated up to mold temperature once and cools down to room temperature afterwards. Both process steps include monotonous temperature changes. During heating, the degree of cure has not reached the point of gelation in most cases. Hence, no residual stress can be frozen in. In the cooling step, the change from rubbery to glassy state should only happen once, so the frozen stresses cannot release.

The approach can be further improved by assuming cr and cg not to be constant. Bogetti and Gillepie [86] describe cr by an interpolation between the values of fully cured and fully uncured, as function of the degree of cure. Adolf and Martin [87] describe the evolution based on the power law. But, due to a lack of experimental data, cr will be assumed to be constant within this study. According to Bernath [14] cg is also assumed to be constant, since the state of cure has neglectable effect on small-strain properties, being relevant for modeling residual stresses and warpage. Furthermore, Bernath verifies the CHILE model to be well suited to model the residual stresses und warpage in case of transverse isotropic material behavior with epoxy matrix in the RTM process.

### **2.3.2 Homogenization and mapping**

#### **2.3.2.1 General procedure**

As mentioned, the parameters have crucial impact on the material during processing, and therefore on the final properties of the part. Hence it is important to describe the material adequately and transform corresponding attributes from the process to the structural simulation, which are often not performed within the same software.

In case of structural simulation of FRPs, an additional material homogenization, based on fiber volume fraction and fiber orientation distribution, is often performed, because, similar to the process simulation, the numerical effort to simulate individual fibers is too high. Therefore, the homogenized FRP is described by an isotropic, transverse isotropic or orthotropic material model. The principle of homogenization is shown in Figure 2.13a.

Since the meshes at start and end of the data transformation (also known as source and target mesh) are not identical, values must be interpolated. This procedure is called mapping, applied for example on temperature and curing field, but also fiber orientation tensors, relevant for the mechanical behavior.

Figure 2.13: Schematic illustration of homogenization (a). Highlighting the influence of different meshes for change of values during mapping (b). Colors indicate different values of attributes

Within the mapping process information may vanish or blur, as schematically shown in Figure 2.13b. The differences of source and target mesh have a significant influence within the mapping, as indicated in the intermediate state in Figure 2.13b, although this influence may become irrelevant, as indicated on the right side. While the mapping of scalars and vectors is explicit, this is not the case for tensors.

#### **2.3.2.2 Homogenization schemes for FRPs**

The aim of homogenization is to calculate effective material parameters with accompanying minimal effort for experimental data or empirical model parameters. Especially in case of FRPs, homogenized material models are often used, since on the one hand, the numerical effort to simulate single fibers and matrix separately is too high and on the other hand the experimental effort to get all necessary material parameters for every material component is also unrealistic.

The literature offers several homogenization schemes for different material groups and load scenarios. In general, homogenization schemes are separated in mean field and full field approaches. The latter are based on representative or statistical micro models to determine the mechanical answer of a macroscopic load, as for example presented by Müller et al. [88,89]. Due to the experimental and numerical effort, these approaches will not be further discussed. Mean field methods can be divided in bounding and estimating methods. The most prominent boundary methods may be the upper bound by Voigt [90] in combination with the lower bound by Reuss [91]. Since these approaches only take fiber volume content into account, but no aspect ratio or fiber orientation distribution, the usage for discontinuous reinforced material is not suitable. Another often used approach for second order bounds is provided by Hashin and Shtrikman [92–94].

Additionally, some often used, well suited and more detailed approaches are the approach by Mori and Tanaka [95] and by Tandon and Weng [96]. The approach of Tandon and Weng considers fiber volume content and aspect ratio and does not need further empirical parameters. Based on the assumption of aligned spheroidal inclusions, a transverse isotropic stiffness tensor is determined for equal strain.

Thermal properties can also be homogenized with different approaches. The specific heat can be simply volume averaged, since it is a scalar and isotropic property. The thermal conductivity can be determined with the approach given by Clayton [32], representing an extension of the Maxwell approach for dilute suspension and valid for higher concentrations. For the coefficient of thermal expansion, Schapery presents an energy based approach [97]. Another often used model is presented by Rosen and Hashin [98], determining the effective properties by superposition of an isothermal problem with surface displacement and an uniform temperature problem with homogenous boundary displacement.

Nevertheless, these approaches only describe transverse isotropic behavior, which is insufficient in most cases for injection molded FRPs. The material properties must be further orientation averaged. Besides the orientation tensors, Advani and Tucker [51] present an orientation averaging scheme for tensor properties, based on a transverse isotropic microstructure. According to [51] a transverse isotropic second order Tensor is fully described by

$$T\_{lj} = B\_1 q\_i q\_j + B\_2 \delta\_{lj\text{\textquotedblleft}}\tag{2.50}$$

with the material specific constants <sup>1</sup> and <sup>2</sup> . Similar to the orientation tensor, described in Section 2.2.3.2 the orientation average of is given by

$$\begin{split} \langle T\_{lj} \rangle &= \int T\_{lj}(q\_l) \psi(q\_l) dq\_l \\ &= B\_1 \langle q\_l q\_j \rangle + B\_2 \langle \mathcal{S}\_{lj} \rangle \\ &= B\_1 A\_{lj} + B\_2 \mathcal{S}\_{lj}. \end{split} \tag{2.51}$$

Here, 〈∙〉 indicates the orientation averaging. Hence, the orientation average of a transverse isotropic tensor can be described by the orientation tensor and the corresponding unity tensor.

A fourth order, transverse isotropic tensor is completely defined by five independent constants 1−5 , as for example a unidirectional continuous FRP, where the mechanical behavior is represented by the five engineering constants. Therefore, the orientation averaging is given by

$$\begin{aligned} \langle T\_{ljkl} \rangle &= B\_1 A\_{ljkl} \\ &+ B\_2 \Big( A\_{lj} \delta\_{kl} + A\_{kl} \delta\_{lj} \Big) \\ &+ B\_3 \Big( A\_{lk} \delta\_{jl} + A\_{ll} \delta\_{jk} + A\_{jl} \delta\_{lk} + A\_{jk} \delta\_{ll} \Big) \\ &+ B\_4 \Big( \delta\_{lj} \delta\_{kl} \Big) + B\_5 \Big( \delta\_{lk} \delta\_{jl} + \delta\_{ll} \delta\_{jk} \Big). \end{aligned} \tag{2.52}$$

47

The combination of the Tandon-Weng homogenization and the Advani-Tucker orientation averaging is successfully applied by Hohberg [99] for SMC material and also used within this work.

#### **2.3.2.3 Mapping schemes**

Several commercially available mapping tools exist to transform data between a process simulation and a structural analysis. In case of fiber reinforced polymers, these are for example MPCCI MapLib [100,101] or Digimat [102]. Further approaches are presented in literature, for example by Reclusado and Nagasawa [103], mapping fiber orientation based on only the first eigenvalue with corresponding eigenvector. MPCCI MapLib has been successfully used for SMC material by Hohberg [99] and Görthofer et al. [104]. The principal of these mapping processes is always the same and shown in Figure 2.14.

Figure 2.14: Schematic illustration of mapping process. Colors indicate different values of attributes. Due to source mesh, every element in the target mesh shows a different color / value

MPCCI MapLib maps independent of the software, which created the source and target mesh. Therefore, the mapping is purely based on geometry and the results and meshes. In case of mapping between different software packages, the meshes must be converted in a neutral data structure [100], being for example VTK within this work. In case of tensors, MapLib interpolates each tensor component independently. This method may lead to a loss of shape of the tensor. Nevertheless, the only tensor mapped within this work is the orientation tensor, being positive definite and having a normed diagonal. Therefore, this mapping strategy is acceptable.

# **3 Development of simulation methods**

### **3.1 Considered interactions in process and material**

A comprehensive process simulation of FRP injection molding includes more than solving the Navier-Stokes equation with standard material models, as described in the Sections 2.2.1.1, 2.2.2 and 2.2.3. Figure 3.1 shows all considered interactions within this work. At first, interactions between the FRP and the tool must be considered. This aspect is solved via boundary conditions, for example the tool temperature as boundary condition for FRP's temperature at the boundary faces. Such interactions are already considered in a state-of-the-art single-phase and isotropic injection molding simulation.

Due to the single-phase approach, the air in the mold is neglected in state-ofthe-art simulations. Within this work, the air is considered as separate phase in the mold, leading to a multi-phase approach with two immiscible fluids for the simulations. This leads to two more interaction pairs, which need to be mentioned. One is, similar to the FRP, the interactions between the air and the tool, which are also solved via the boundary conditions. The other one is the interaction between FRP and air, which is solved within the multi-phase approach as described in Section 3.2 and [25,28].

Furthermore, internal material interactions are regarded within this work. The interactions between fibers and matrix are separated, so the influence from fibers on the matrix and vice versa are modeled individually and with different approaches.

Due to their non-spherical shape, the fibers induce stresses in the matrix in an anisotropic way. Therefore, the FRP's viscosity is modeled by a fourth order viscosity tensor, depending on non-Newtonian matrix viscosity, fiber orientation, length and volume fraction. The viscosity tensor and the underlying micro-mechanical models are described in Section and [26,38,105].

Figure 3.1: Macroscopic interactions between FRP (green), air (light blue) and tool (shaded) with corresponding simulation aspects. Material internal interactions between fibers (blue) and matrix (yellow) with corresponding simulation approaches

On the opposite, the relative velocity between matrix and fibers induces hydrodynamic forces in the fibers, leading to fiber re-orientation and may also lead to fiber breakage. The hydrodynamic forces are calculated to support fiber breakage simulations, depending on matrix viscosity, flow field and fibers. The corresponding methods are described in Section 3.4.1 and [37,106].

Besides the matrix, fibers interact with other fibers. These interactions induce lubrication, friction and normal forces in the fibers. Hence, these forces are approximated as function of matrix viscosity, fiber orientation, geometry and volume fraction, as described in Section 3.4.2 and [106–110].

To keep the numerical effort in an acceptable range to perform a macroscopic process simulation of a whole part, the fiber orientation is described with state-of-the-art models and orientation tensors, described in Section 2.2.3. Hence, for the hydrodynamic forces on fibers, only the portion in fiber direction is considered, leading to buckling and breaking. The non-fiberdirection component of the forces leads to re-orientation and would therefore be an overshoot, because the orientation is calculated separately.

### **3.2 Multiphase approach<sup>1</sup>**

### **3.2.1 Volume-of-fluid method**

A multi-phase approach requires a methodology to differentiate between the phases and separate material properties in the calculations. Many different approaches to separate two, or more, phases in a CFD simulation have been published and can be found in several fluid dynamic books. Within this work, the FRP and air are considered as homogeneous fiber matrix suspension and gas mixture. Therefore, the multi-phase approach distinguishes only between two phases, which are immiscible. A well-suited approach for these conditions is the volume-of-fluid method (VOF) with one VOF-factor ∈ [0,1].

Within this work, = 0 equals pure air and = 1 pure FRP. Cells with ∈ (0,1) are partially filled. Effective material specific properties can be calculated by volume averaging

$$\cdot \mathbf{i}\_{\rm eff} = \alpha |\cdot|\_{\rm FRP} + (1 - \alpha) |\cdot|\_{\rm air} . \tag{3.1}$$

<sup>1</sup> The methods presented in this Section are published in [25].

#### **3.2.2 Phase-dependent boundary condition**

Due to the presence of air within the mold, a phase-dependent boundary condition must be formulated, the way that air can leave the mold but FRP not. The boundary condition is defined as function of to consider the phasedependence and applies for . The air is assumed to leave the mold without a change of state, so

$$\frac{\partial U\_l}{\partial \mathbf{x}\_j} = \begin{pmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix} \text{ for } a = \mathbf{0}, \tag{3.2}$$

being a Neumann boundary condition named *zeroGradient*. Furthermore, the boundary condition should act like a wall for the FRP resulting in a Dirichlet boundary condition with

$$U\_l = \begin{pmatrix} 0 \ 0 \ 0 \end{pmatrix} \text{ for } \alpha = 1,\tag{3.3}$$

named zero. Hence a no-slip boundary condition is applied, which is a simplification for RIM, but creates good results with less need of model parameters (see Section 2.2.1.5).

There are different options to define the condition for ∈ (0,1). One is to define a fixed switchover for = 1, leading to a loss of material caused by outflow, while cells are partially filled. To reduce this loss, the switchover can be defined at any point ≤ 1, which leads to in-mold air, unable to leave. Within this work, a compromise between these two cases is chosen. Similar to material properties and according to [25,26] the two boundary conditions *zeroGradient* and zero are interpolated with . The final velocity boundary condition is described schematically with

$$\mathcal{U}\_{\text{bound}} = \mathcal{a}U\_{\text{zero}} + (1 - \alpha)zeroGradient. \tag{3.4}$$

Numerical studies and a verification of the boundary condition are presented in Section 4.1.1. To minimize the material loss caused by this formulation, the condition is only applied for a specific region of the mold as shown in Figure 3.2. This region, called outlet, is defined at the end of the flow path, where the venting slots are positioned in the real process. The other boundary areas (besides the inlet) are impermeable walls, so neither FRP nor air may leave the mold.

Figure 3.2: Different boundary regions for an arbitrary mold filling model. Inlet in blue, impermeable walls in grey and outlet with phase-dependent boundary condition in red

## **3.3 Anisotropic viscosity modeling<sup>2</sup>**

### **3.3.1 Theory**

#### **3.3.1.1 Transversely isotropic fluids**

The real viscosity of a FRP is a complex combination of temperature, chemical reaction, shearing and, due to presence of fibers, anisotropy. Within this Section the focus is on description of anisotropic viscosity due to fibers.

Since fibers are non-spherical, stress builds up along the surface in an anisotropic way. This effect is similar to FRPs in solid mechanics and hence, the simplest anisotropic behavior is transversely isotropic. This behavior is present if all fibers are aligned parallel in one direction (cf. Figure 2.11b and Figure 3.3). Gibson [43] describes a transversely fluid analogy to solid mechanics by

<sup>2</sup> The methods presented in this Section are published in [26].

$$\mathcal{E}\_{lj} = \mathcal{S}\_{ljkl}^{\text{transverse}} \sigma\_{kl}$$

$$\mathcal{S}\_{ljkl}^{\text{transverse}} = \begin{bmatrix} \mathcal{S}\_{1111} & \mathcal{S}\_{1122} & \mathcal{S}\_{1122} & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ & \mathcal{S}\_{2222} & \mathcal{S}\_{2233} & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ & & \mathcal{S}\_{2222} & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ & & & \mathcal{S}\_{2222} & -\mathcal{S}\_{1111} & \mathbf{0} & \mathbf{0} \\ & & & \text{sym.} & & \mathcal{S}\_{1212} & \mathbf{0} \\ & & & & & \mathcal{S}\_{1212} \end{bmatrix}. \quad (3.5)$$

with as forth order fluidity tensor.

Additionally, Gibson assumes incompressibility, so it is

$$
\varepsilon\_{11} + \varepsilon\_{22} + \varepsilon\_{33} = 0,\tag{3.6}
$$

leading to

$$\begin{aligned} \left(\sigma\_{22} + \sigma\_{33}\right) \left(\mathbb{S}\_{2222} + \mathbb{S}\_{1122} + \mathbb{S}\_{2233}\right) \\ = -\sigma\_{11} \left(\mathbb{S}\_{1111} + 2\mathbb{S}\_{1122}\right) .\end{aligned} \tag{3.7}$$

Since Eq. (3.7) must be valid for any arbitrary stress state, two more conditions are given by

$$\begin{aligned} \{\mathcal{S}\_{2222} + \mathcal{S}\_{1122} + \mathcal{S}\_{2233}\} &= 0, \\ \{\mathcal{S}\_{1111} + 2\mathcal{S}\_{1122}\} &= 0, \end{aligned} \tag{3.8}$$

hence the number of independent material constants reduces to three and the fluidity tensor can be rewritten as

$$
\begin{bmatrix}
\mathbf{S}\_{1111} & \frac{-S\_{1111}}{2} & \frac{-S\_{1111}}{2} & \mathbf{0} & \mathbf{0} & \mathbf{0} \\
\frac{S\_{1111} + S\_{2323}}{4} & \frac{S\_{1111} - S\_{2323}}{4} & \mathbf{0} & \mathbf{0} & \mathbf{0} \\
\frac{S\_{1111} + S\_{2323}}{4} & \mathbf{0} & \mathbf{0} & \mathbf{0} \\
& & \frac{S\_{1111} + S\_{2323}}{4} & \mathbf{0} & \mathbf{0} \\
& & & \mathbf{S}\_{2323} & \mathbf{0} & \mathbf{0} \\
& & & & \mathbf{S}\_{1212} & \mathbf{0} \\
\end{bmatrix},\tag{3.9}
$$

Within [43], these three constants are independent of fiber attributes.

#### **3.3.1.2 Dependence on fiber attributes of transversely isotropic fluids**

A solution to describe the three constants of a transversely fluid as function of fiber attributes is given by Pipes et al. [105]. Besides the full alignment and the incompressibility Pipes et al. define a few more assumptions. The fiber ends are touching and the fibers are arranged in square or hexagonal packaging. Fiber ends in one row are next to fiber centers in neighbor rows, also in hexagonal packaging, which is a simplification. All fibers have the same diameter <sup>f</sup> , length <sup>f</sup> and neighbor distance ℎ<sup>f</sup> . The rows have an offset of f⁄2. The variation of velocity in fiber direction is linear (<sup>1</sup> = ̇11<sup>1</sup> +const.), so ̇<sup>11</sup> = const. Fibers move with the velocity of their center. A schematic visualization of such an assembly with resulting shearing is given in Figure 3.3.

Figure 3.3: Aligned fibers in a fluid. Shearing due to 1,2 > 1,1 .

Based on the assumptions and information in Figure 3.3 it is

$$U\_{1,2} - U\_{1,1} = \varepsilon\_{11} L\_{\text{f}} / 2 \tag{3.10}$$

and

$$
\dot{\varepsilon}\_{12} = \dot{\chi} = \dot{\varepsilon}\_{11} L\_{\rm f} / \left( 2(h\_{\rm f} - d\_{\rm f}) \right). \tag{3.11}
$$

For a fiber formation as shown in Figure 3.3, the fiber volume fraction is given by

$$
\Phi\_{\rm f} = \Phi\_{\rm max} \left( \frac{d\_{\rm f}}{h\_{\rm f}} \right)^2,\tag{3.12}
$$

with Φmax being the maximum possible fiber volume fraction, which is /√12 for hexagonal and /4 for square packing. Since hexagonal packing represents a stable equilibrium it is the more realistic case. Therefore, hexagonal packaging is assumed within this work.

According to [105] the average force ̂ acting in fiber direction at the fiber midpoint is given by

$$
\hat{F} = \frac{2\sigma\_{11}}{\Phi\_{\rm f}} \left(\frac{\pi d\_{\rm f}^2}{4}\right). \tag{3.13}
$$

On the opposite, the force due to shear stress <sup>12</sup> on half of the fiber surface is

$$
\hat{F} = \pi\_{12} \left( \frac{\pi \ L\_{\text{f}} \, d\_{\text{f}}}{2} \right) = \eta\_{\text{M}} \dot{\mathcal{Y}} \left( \frac{\pi \ L\_{\text{f}} \, d\_{\text{f}}}{2} \right). \tag{3.14}
$$

Building up the force equilibrium leads to

$$
\eta\_{\mathsf{M}} \dot{\mathsf{y}} = \frac{\sigma\_{11}}{\Phi\_{\mathsf{f}}} \left( \frac{d\_{\mathsf{f}}}{L\_{\mathsf{f}}} \right) = \frac{\sigma\_{11}}{\Phi\_{\mathsf{f}} \mathsf{r}\_{\mathsf{f}}}.\tag{3.15}
$$

Finally, the combination of Eq. (3.11), (3.12) and (3.15) yields to

$$\eta\_{11} = \frac{\sigma\_{11}}{\varepsilon\_{11}} = \frac{\eta\_{\text{M}} \Phi\_{\text{f}}}{2} \left( \frac{\sqrt{\Phi\_{\text{f}} / \Phi\_{\text{max}}}}{1 - \sqrt{\Phi\_{\text{f}} / \Phi\_{\text{max}}}} \right) r\_{\text{f}}^2. \tag{3.16}$$

The axial elongational viscosity <sup>11</sup> is identical to the entry <sup>1111</sup> in the fourth order viscosity tensor for a complete micro-mechanical approach. Nevertheless, later in this work the formulation of <sup>1111</sup> will be defined

differently due to homogenization based on this micro mechanical approach (see Section 3.3.2). Therefore, the notation <sup>11</sup> is chosen.

Similar to <sup>11</sup> the axial shear viscosity <sup>12</sup> and transverse shear viscosity <sup>23</sup> can be determined, as shown in [105]. The results are

$$\eta\_{12} = \frac{\eta\_{\rm M}}{2} \left( \frac{2 - \sqrt{\Phi\_{\rm f}/\Phi\_{\rm max}}}{1 - \sqrt{\Phi\_{\rm f}/\Phi\_{\rm max}}} \right) \tag{3.17}$$

and

$$
\eta\_{23} = \frac{\eta\_{\rm M}}{(1 - \Phi\_{\rm f}/\Phi\_{\rm max})^2}.\tag{3.18}
$$

With these three parameters (11, <sup>12</sup> and 23) it is possible to fully describe the viscous behavior of a transversely fluid with fiber arrangement as shown in Figure 3.3. The matrix viscosity can be assumed constant, as in [38], or be calculated with an viscosity model such as in [26]. Within this work, the matrix viscosity is calculated with the Castro-Macosko model described in Section 2.2.2.3, if not explicitly mentioned differently.

Pipes et al. published further investigations to describe the material behavior on microscopic scale for a fiber-fiber overlap ≠ f⁄2 [111] or within a power-law matrix [112]. Nevertheless, these further approaches are not suitable for a macroscopic process simulation with orientation tensors, since too much information, for example the overlap length of individual fibers, is not defined.

#### **3.3.2 Viscosity tensor of discontinuous FRPs**

#### **3.3.2.1 Viscosity tensor of transversely isotropic FRPs**

Based on the assumptions and formulations of Section 3.3.1 and by setting <sup>1111</sup> = 1/11, <sup>2323</sup> = 1/<sup>23</sup> and <sup>1212</sup> = 1/<sup>12</sup> the viscosity would be defined as inverse of the fluidity tensor (Eq. (3.9)), but, due to the incompressibility assumption, the fluidity tensor is singular. Sommer et al. [38]

present a solution to formulate a fourth order viscosity tensor for the assumption of incompressibility. In a first step, a pseudoinverse viscosity tensor for a transversely isotropic fluid is built. The methodology is given by Loredo and Klöcker [113], determining the stiffness tensor of an incompressible solid material. The pseudoinverse is defined as

$$\mathcal{S}\_{ljkl}^{-1} = \left(\mathcal{S}\_{ljkl} + \left(\frac{1}{3}\right)\delta\_{lj}\delta\_{kl}\right)^{-1} - \left(\frac{1}{3}\right)\delta\_{lj}\delta\_{kl}.\tag{3.19}$$

Applying this for the fluidity tensor as defined in Eq. (3.9) and further replacing <sup>1111</sup> = 1/11, <sup>2323</sup> = 1/<sup>23</sup> and <sup>1212</sup> = 1/<sup>12</sup> leads to

$$
\eta\_{ijkl}^{\text{transverse}} = \frac{1}{\eta} \begin{bmatrix}
4\eta\_{11} & -2\eta\_{11} & -2\eta\_{11} & 0 & 0 & 0 \\
\eta\_{11} + 9\eta\_{23}\eta\_{11} - 9\eta\_{23} & 0 & 0 & 0 \\
& \eta\_{11} + 9\eta\_{23} & 0 & 0 & 0 \\
& & \eta\_{23} & 0 & 0 \\
& & & \eta\_{12} & 0 \\
& & & & \eta\_{12}
\end{bmatrix} .\tag{3.20}
$$

#### **3.3.2.2 Orientation averaging of the viscosity tensor**

The tensor given by Eq. (3.20) is only valid for transversely isotropic, meaning fully aligned, materials. However, this state is never reached in a real process, although it might be a good approximation near the surface and in thin part regions. To close this gap and apply the tensor also for less orientated regions, the orientation averaging by Advani and Tucker, introduced in Section 2.3.2.2 (Eq. (2.52)) is used, similar to [26,38]. Finally, the orientation averaged fourth order viscosity tensor is defined as

$$\begin{split} \langle \eta\_{ljkl}^{\text{FRP}} \rangle &= (\eta\_{11} - 4\eta\_{12} + \eta\_{23}) A\_{ljkl} \\ &+ \left( -\frac{\eta\_{11}}{3} + \eta\_{23} \right) \{ A\_{lj} \delta\_{kl} + A\_{kl} \delta\_{lj} \} \\ &+ (\eta\_{12} - \eta\_{23}) \begin{pmatrix} A\_{lk} \delta\_{lj} + \\ A\_{lk} \delta\_{lk} + A\_{jl} \delta\_{lk} + A\_{jk} \delta\_{ll} \end{pmatrix} \\ &+ \left( \frac{\eta\_{11}}{9} - \eta\_{23} \right) \{ \delta\_{lj} \delta\_{kl} \} \\ &+ \eta\_{23} \{ \delta\_{lk} \delta\_{jl} + \delta\_{ll} \delta\_{jk} \} .\end{split} \tag{3.21}$$

This formulation fulfills Eq. (3.9) and generates similar results compared to the formulations of Ericksen [114] and Hinch and Leal [115]. The viscosity is calculated anisotropic, depending on fiber orientation, fiber volume fraction, fiber aspect ratio and matrix viscosity, represented by the three parameters 11, <sup>12</sup> and <sup>23</sup> given in Eq. (3.16)-(3.18).

#### **3.3.2.3 Extension to multi-phase modeling**

The work of Sommer et al. [38] focusses on modeling a single-phase compression molding process. The approach must be extended to be applicable for a multi-phase simulation with FRP and air, as described in [26]. The viscosity is a material specific variable, interpolated with the VOF method (cf. Section 3.2.1). Hence the effective viscosity eff is

$$
\eta\_{ljkl}^{\rm eff} = a \cdot \eta\_{ljkl}^{\rm FRP} + (1 - a) \cdot \eta\_{ljkl\nu}^{\rm air} \tag{3.22}
$$

with

$$
\eta\_{ijkl}^{\text{air}} = \eta\_{\text{air}} \begin{bmatrix}
4/3 & -2/3 & -2/3 & 0 & 0 & 0 \\
& 4/3 & -2/3 & 0 & 0 & 0 \\
& & 4/3 & 0 & 0 & 0 \\
& & & 1 & 0 & 0 \\
& & & & 1 & 0 \\
\end{bmatrix},
\tag{3.23}
$$

where air is the scalar viscosity of air, assumed to be constant within this work. The fourth order viscosity tensor of air is isotropic and the formulation fulfills the condition

$$\begin{aligned} \eta\_{ljkl}^{\text{air}} D\_{kl} &= \\ 2\eta\_{\text{air}} \left( D\_{lj} - \frac{1}{3} \frac{\partial U\_k}{\partial \varkappa\_k} \delta\_{lj} \right), \end{aligned} \tag{3.24}$$

so the results for the viscous stress tensor are identical for the scalar and fourth order tensor formulation (see Section 2.2.1.1, Eq. (2.4) and Eq. (2.5)).

#### **3.3.2.4 Implementation in OpenFOAM**

The momentum equation (Eq. (2.2)) is solved implicitly within OpenFOAM 4.1. Since the viscous stress tensor directly depends on the velocity gradient, Eq. (2.4) is used during solving, so the viscous stress can change and the equilibrium is reached. This term cannot be directly replaced by eff , since the implicit solving algorithm is only defined for scalars, vectors and second order tensors in Cartesian coordinate systems. Nevertheless, implicit solving is more stable and robust. To gain more stability during solving, the effective viscosity tensor is split into an isotropic and anisotropic part so that

$$
\eta\_{ijkl}^{\rm eff} = \eta\_{ijkl}^{\rm eff,iso} + \eta\_{ijkl}^{\rm eff,aniso},\tag{3.25}
$$

as presented in [26]. The isotropic part can be represented by a scalar value and hence be solved implicitly without loss of information. According to Bertóti and Böhlke [116] the effective scalar shear viscosity of an anisotropic fluid is

$$
\eta\_{\rm eff,lso} = \frac{1}{10} \eta\_{ljkl}^{\rm eff} P\_{ljkl}^2 \tag{3.26}
$$

where 2 is the second projector tensor of fourth order. The corresponding fourth order tensor is defined similar to the viscosity of air as

$$
\eta\_{ijkl}^{\text{eff,iso}} = \eta\_{\text{eff,iso}} \begin{bmatrix}
4/3 & -2/3 & -2/3 & 0 & 0 & 0 \\
& 4/3 & -2/3 & 0 & 0 & 0 \\
& & 4/3 & 0 & 0 & 0 \\
& & & 1 & 0 & 0 \\
& & & & 1 & 0 \\
& & & & & 1
\end{bmatrix} . \tag{3.27}
$$

The anisotropic part can be determined by eff,aniso = eff − eff,iso . Finally, the viscous stress is calculated with

$$\begin{split} \pi\_{lj} &= \eta\_{ljkl}^{\mathrm{eff}} D\_{kl} = \left( \eta\_{ljkl}^{\mathrm{eff,iso}} + \eta\_{ljkl}^{\mathrm{eff,aniso}} \right) D\_{kl} \\ &= \underbrace{2 \eta\_{\mathrm{eff,iso}} \left( D\_{lj} - \frac{1}{3} \frac{\partial U\_k}{\partial \boldsymbol{x}\_k} \delta\_{lj} \right)}\_{\text{solvent implicitly}} + \underbrace{\eta\_{\boldsymbol{l}jkl}^{\mathrm{eff,aniso}} D\_{kl}}\_{\text{solvent explicitly}}. \end{split} \tag{3.28}$$

### **3.4 Fiber forces and breakage modeling<sup>3</sup>**

#### **3.4.1 Hydrodynamic forces from fluid on fibers**

The anisotropic viscosity tensor, presented in Section 3.3 represents the fiber to fluid part in the stress equilibrium, but not the fluid on fiber part. The matrix introduces forces in the fibers, due to relative velocity, leading to reorientation or in some cases to fiber damage. This fluid dynamic phenomenon is known as hydrodynamic forces. In principal the effect is similar to aerodynamics and therefore, the same approaches can be adapted.

Since fibers are non-spherical bodies, the hydrodynamic force hyd is a combination of drag force d and lift force li. The first appears whenever a fluid has a relative velocity and contact to a rigid body, while the latter occurs due the inhomogeneous stress distribution on the surface, caused by the nonspherical geometry [20].

The literature offers several studies on the modeling of hydrodynamic forces on fibers, such as for example work of Phan-Thien et al. [117], Lindström and Uesaka [118], Dinh and Amstrong [72] and Meyer et al. [37]. Within this work, the approach of Meyer et al. is adapted, since it has proven to be well suited for modeling forces on glass fibers, or fiber bundles, in a macroscopic process simulation, although Meyer et al. focus on compression molding simulations. Furthermore, the approach provides the possibility to distinguish

<sup>3</sup> The methods for fiber force modeling presented in this Section are published in [106].

between drag and lift force, which is not directly possible in other approaches.

#### **3.4.1.1 Modeling of hydrodynamic drag and lift force**

The drag force of a sphere is given by Stokes law, so

$$F\_l^d = 3\pi \eta\_{\mathsf{M}} d\_{\mathsf{sp}} \Delta U\_{l\prime} \tag{3.29}$$

with sp being the diameter of the sphere and Δ the relative velocity. According to Batchelor [20], the absolute hydrodynamic resistance is proportional to Msp‖Δ‖ and independent of the body's shape, if the flow is Newtonian, incompressible and the internal forces are neglectable compared to the viscous forces so ≪ 1. Hence, Meyer et al. define an equivalent diameter

$$d\_{\rm sp} = k\_{\rm d}(r\_{\rm f}, \phi) d\_{\rm f} \tag{3.30}$$

with the dimensionless correction coefficient d, depending on fiber aspect ratio and relative angle between fiber and relative velocity .

Besides the drag force, the lift force should also be considered. Meyer et al. [37] present a similar approach using a dimensionless correction coefficient li. Of course, the lift force acts not in direction of the relative velocity, but perpendicular to the relative velocity and the fiber axis. Hence, the lift force direction is defined as

$$p\_l = (q\_l \times \mathbb{I} \| \Delta U\_l \|) \times \| \Delta U\_l \|. \tag{3.31}$$

Similar to the drag force, but acting in a different direction, the lift force is then defined by

$$F\_l^{\rm ll} = 3\pi \eta\_{\rm M} k\_{\rm ll}(r\_{\rm f}, \phi) d\_{\rm f} \|\Delta U\_l\| \|p\_l\|. \tag{3.32}$$

Finally, the complete hydrodynamic force on one fiber is given by

$$\begin{aligned} \boldsymbol{F}\_{l}^{\text{hyd}} &= \boldsymbol{F}\_{l}^{\text{d}} + \boldsymbol{F}\_{l}^{\text{li}} \\ = 3\pi\eta\_{\text{M}}\boldsymbol{d}\_{\text{f}}(\boldsymbol{k}\_{\text{d}}\Delta\boldsymbol{U}\_{l} + \boldsymbol{k}\_{\text{li}} \|\Delta\boldsymbol{U}\_{l}\| \|\boldsymbol{p}\_{l}\|). \end{aligned} \tag{3.33}$$

The ⟦∙⟧ brackets represent a normed vector, defined as ⟦ ⟧ = /‖‖.

#### **3.4.1.2 Definition of correction coefficients and**

Within [37], <sup>d</sup> and li are approximated, based on numerical experiments computing the stress on single fibers in a simple shear flow with different orientations and aspect ratios. Since the hydrodynamic force is the result of an interaction between fluid and fiber on the fiber surface, it is also given by

$$F\_l^{\text{hyd}} = \int \sigma\_{lj} n\_j \text{d}A\_{\text{f}} \,\tag{3.34}$$

with <sup>f</sup> describing the fiber surface and being the corresponding normal vector. Here, the fiber ends are neglected. Meyer et al. [37] performed numerical experiments with a fixed fiber within a flow field with known viscosity and velocity. Combining Eq. (3.33) and (3.34) leads to

$$k\_{\rm d} = \frac{1}{3\pi\eta\_{\rm M}d\_{\rm f} \|\Delta U\_{\rm l}\|} \int \sigma\_{\rm l} n\_{\rm l} \mathbf{d}A\_{\rm f} \tag{3.35}$$

and

$$k\_{\rm li} = \frac{1}{3\pi\eta\_{\rm M}d\_{\rm f} || \Delta U\_{\rm l}||} \int \sigma\_{\rm f} n\_{\rm f} \mathrm{d}A\_{\rm f} \tag{3.36}$$

where the index and indicate the components vertical and horizontal on the surface here. Meyer et al. solve Eq. (3.35) and Eq. (3.36) for Δ = (1 0 0) mm/s and ∈ [0°, 90°] as well as different aspect ratios. The angle interval represents all possible angles between fiber and relative velocity, since a fiber has no front or back end. Within [37] numerical fits are performed to approximate <sup>d</sup> and li as function of aspect ratio and angle. The results are given by

$$k\_{\rm d} = 1 - \alpha\_{\rm hyd}(r\_{\rm f} - 1)\cos(2\phi) + \beta\_{\rm hyd}(r\_{\rm f} - 1) \tag{3.37}$$

and

$$k\_{\rm li} = a\_{\rm hyd} (r\_{\rm f} - 1) \sin(2\phi),\tag{3.38}$$

with hyd = 0.09 and hyd = 0.3125. By these approximations <sup>d</sup> shows a steady increase for <sup>f</sup> > 1 and for an increase of , while it is constant for <sup>f</sup> = 1 , applying [20] in case of spherical inclusions. Whereas li is zero for <sup>f</sup> = 1 or = {0° 90°} and has a maximum for = 45°, since fibers have a symmetric surface.

#### **3.4.1.3 Application of hydrodynamic force modeling in homogenized material**

The work of Meyer et al. [37] includes a material model approach with individual fibers (or fiber bundles) enabling the calculation of the hydrodynamic force as given by Eq. (3.33). Within a homogenized material three important aspects are unknown. The first one is that the relative velocity between fibers and fluid is unknown. Secondly, the angle between the fiber and Δ is unknown. Consequently and thirdly, the direction of the lift force is also unknown. In a first step, the relative velocity is approximated. Similar to [37] and [106], Δ can be approximated by

$$
\Delta U\_l = \sum\_n \frac{w\_{nm}}{W\_m} (U\_{lm} - U\_{lm}) \tag{3.39}
$$

with = exp (−(9 2 )/(2<sup>f</sup> 2 )) and = ∑ . The indices and indicate different cells and is the distance between the centers of these cells. Contrary to [37] is not the velocity of the bundle part within this work, but the velocity in cell at the cell center and is assumed that all fibers within this cell have on average the same velocity as the cell center. This simplification is justified by homogenization and the fact that no fiber matrix separation is assumed [106]. Hence, Δ depends on the velocities of the considered neighbor cells. The value of is chosen so 1/ ∑ ≥ 1.5 <sup>f</sup> , but at least one generation of cell neighbors is always taken into account.

The next step is to approximate the angle between Δ and the fibers. Without knowledge of single fiber orientation, the next meaningful approximation is to determine some kind of average angle. This angle can be defined in relation to a reference fiber. In structural analysis, the material is often described transverse isotropic with fibers along the first eigenvector, as described in Section 2.3.2 and [103]. Hence, the eigenvector can be interpreted as such a representative fiber. Since there are always three eigenvectors, there are three reference fibers with exact known orientation and the orientation probability of the corresponding eigenvalue can be determined for every known second order orientation tensor. This assumption fits to the definition of the orientation tensor, since it is

$$A\_{lj} = \int \psi(q\_l) q\_l q\_j \mathrm{d}q\_l = \sum\_{n=1}^{3} \lambda\_n \nu\_{ln} \nu\_{jn} \tag{3.40}$$

in case of three fibers. Here, is the -th eigenvector of and the corresponding eigenvalue.

By knowing the orientation of the eigenvectors, the angle to the relative velocity is given by

$$\phi\_n = \arccos\left(\frac{\Delta U\_l \nu\_{ln}}{\|\Delta U\_l\|}\right) \tag{3.41}$$

and the direction of the lift force by

$$p\_{ln} = (\nu\_{ln} \times \mathbb{[\Delta U\_l]}) \times \mathbb{[\Delta U\_l]}.\tag{3.42}$$

Based on this information, Eq. (3.33) can be solved and hyd can be approximated in a discretized mesh with fiber orientation tensors [106]. Due to the dependence of and on the three eigenvectors, hyd is calculated three times in every cell. Furthermore, <sup>d</sup> and li depend on the fiber aspect ratio, so in case of fiber length distributions, hyd is determined three times for every possible fiber length.

#### **3.4.2 Fiber-fiber interactions**

#### **3.4.2.1 Approximation of fiber-fiber contact points**

During processing, fibers are in contact, which gets more significant for more highly filled suspensions. Whenever there is a contact between two fibers, contact forces act at these contacts. Since the forces appear on every contact point, in a first step, the number of these contact points needs to be determined.

Toll [119] describes that the average number of fiber centerlines, intersecting an average test fiber is exactly given by

$$N\_{\rm fl} = n\_{\rm f} L\_{\rm f}^2 d\_{\rm f} f + \,^1/\_{\rm f} \pi n\_{\rm f} L\_{\rm f} d\_{\rm f}^2 (g+1),\tag{3.43}$$

with <sup>f</sup> being the number of fibers per volume and and are the scalar invariants of the fiber orientation distribution, given by

$$f = \int \int |q\_l^a \times q\_j^\beta| \psi(q\_l^a) \psi(q\_j^\beta) \mathrm{d}q\_l^a \mathrm{d}q\_j^\beta \tag{3.44}$$

and

$$g = \int \int \left| q\_l^\alpha q\_l^\beta \right| \psi(q\_l^\alpha) \psi(q\_l^\beta) \mathrm{d}q\_l^\alpha \mathrm{d}q\_l^\beta \,. \tag{3.45}$$

Here, the indices α and β indicate different individual fibers. In a later study, Toll shows, that the average number of fiber contact points fc per fiber is determined by replacing <sup>f</sup> within Eq. (3.43) with 2<sup>f</sup> [107], so

$$N\_{\rm fc} = 2n\_{\rm f}L\_{\rm f}^2d\_{\rm f}f + \pi n\_{\rm f}L\_{\rm f}d\_{\rm f}^2(g+1). \tag{3.46}$$

The fiber volume fraction <sup>f</sup> is given by (1/4)ff<sup>f</sup> 2 , so the contact points may also be calculated as function of the volume fraction by

$$N\_{\rm fc} = \,^{\mathcal{B}}\!/\_{\mathcal{B}}\Phi\_{\rm f}r\_{\rm f}f + 4\phi\_{\rm f}(g+1). \tag{3.47}$$

Finally, the number of contact points is given, depending on fiber volume fraction, aspect ratio and orientation. Although this approach is a great advance in FRP process modeling, one disadvantage is the dependence of and on single fiber orientations, being an information, which is not available in a macroscopic process simulation.

Therefore, Férec et al. [120] present an approach to approximate the number of contacts points, with information provided by the orientation tensors. Within the study, new rheological models for fiber orientation are developed, including a so-called interaction tensor I , defined as

$$
\delta b^1\_{lj} = \int \int q^a\_l q^a\_j |q^a\_k \times q^\beta\_l| \psi(q^a\_l) \psi(q^\beta\_j) \mathrm{d}q^a\_l \mathrm{d}q^\beta\_j. \tag{3.48}
$$

Due to the normed orientation vectors, it is

$$f = b^{1}\_{ll}.\tag{3.49}$$

Furthermore, | <sup>α</sup> × β | = 1− α β α β , so

$$\begin{split} & \quad \mathbf{b}\_{lj}^{\mathbb{I}} \\ &= \int \int \left( q\_l^\alpha q\_j^\alpha - q\_k^\alpha q\_k^\beta q\_l^\alpha q\_l^\beta q\_l^\alpha q\_j^\alpha \right) \\ & \quad \quad \quad \quad \psi(q\_l^\alpha)\psi\left(q\_j^\beta\right) \mathrm{d}q\_l^\alpha \mathrm{d}q\_j^\beta \\ &= \int \int \left( q\_l^\alpha q\_j^\alpha - q\_l^\alpha q\_j^\alpha q\_k^\alpha q\_l^\alpha q\_l^\beta q\_j^\beta \right) \\ & \quad \quad \quad \quad \psi(q\_l^\alpha)\psi\left(q\_j^\beta\right) \mathrm{d}q\_l^\alpha \mathrm{d}q\_j^\beta \\ & \quad \quad \quad = A\_{lj} - A\_{ljkl}A\_{kl}. \end{split} \tag{3.50}$$

Additionally, Férec et al. [120] introduce a fitting factor, so values for calculated with Eq. (3.44) and Eq. (3.49) are identical in case of quasiisotropic fiber orientation and the final interaction tensor is finally given by

$$\mathbf{b}\_{lj}^{\mathbf{l}} = {}^{\mathbf{3}\pi} \slash\_{\mathbf{3}} \{ A\_{lj} - A\_{ljkl} A\_{kl} \kern-1.74(-1.74){\times}}.\tag{3.51}$$

According to Toll [107] the term 4<sup>f</sup> ( + 1) in Eq. (3.47) can be neglected in case of high aspect ratios <sup>f</sup> , only appearing in the first term, which becomes dominant in this case. Férec et al. neglect this term too and no approximation for the scalar invariant is given in [120]. Nevertheless, for high orientations decreases, while increases and becomes more important as it is also the case in short fiber materials. Within this work, the two novel approaches to determine and presented in [106] are used, also only depending on information, given by . In this way, the average contact points of a fiber within a fiber network can also be determined adequately in case of high orientations and short fiber materials in a macroscopic simulation.

In the first approach the eigenvectors and eigenvalues of are used to determine and . Similar to Eq. (3.40) and can be reformulated as

$$f = \sum\_{n,m=1}^{3} |\nu\_{ln} \times \nu\_{jm}| \lambda\_n \lambda\_m \tag{3.52}$$

and

$$g = \sum\_{n,m=1}^{3} |\nu\_{ln}\nu\_{lm}| \lambda\_n \lambda\_m. \tag{3.53}$$

Since ∈ sym. it is ⊥ , so || = 0 for ≠ . By further regarding | | = 1, and simplify to

$$f = \,^3\pi \left/ \_\otimes \left( 2\lambda\_1 \lambda\_2 + 2\lambda\_1 \lambda\_3 + 2\lambda\_2 \lambda\_3 \right) \right. \tag{3.54}$$

and

$$g = \lambda\_1 \lambda\_1 + \lambda\_2 \lambda\_2 + \lambda\_3 \lambda\_3. \tag{3.55}$$

As shown in [106] and 4.1.3.3 the calculation of with the interaction tensor and Eq. (3.54) creates identical results. Therefore, the correction factor 3/8 is also used for this approach, to fit the value in case of quasi-isotropic orientation. In the case of full alignment, it is = 1 for Eq. (3.45) and Eq. (3.55) and since this invariant is more important for more highly orientated states, no fitting factor is used. Eq. (3.54) and Eq. (3.55) show, that and can be approximated as function of the eigenvalues with corresponding factors in a polynomial way [106]. Hence, in a second approach, the invariants can also be approximated by

$$f, g = \sum\_{n,m=1}^{3} M\_{nm} \tilde{\lambda}\_n \tilde{\lambda}\_m \tag{3.56}$$

with containing the polynomial coefficients and ̃ = (<sup>1</sup> <sup>2</sup> 1), since the eigenvalues are normed and all information is provided by two eigenvalues. The entries of are determined in [106] and given in Table 3.1.


Table 3.1: Coefficients of polynomial fit of scalar invariants and [106]

#### **3.4.2.2 Approximation of fiber contact forces**

The three forces, present at every fiber contact point, are the normal force n , the friction force fr and the lubrication force lu. The direction of the normal force depends on the orientations of the two fibers in contact. Since the reference fibers are represented by the eigenvectors within this work, the normal force cannot be approximated adequately due to ∈ sym. and ⊥

 for ≠ . Nevertheless, the amount of the average normal force, as for example used in [121,108,109] and relevant for friction, can be given by

$$\|\|F\_{\rm l}^{\rm n}\|\| = \frac{32}{5\pi^2} E\_{\rm f} d\_{\rm f}^2 \Phi\_{\rm f}^3 f^3. \tag{3.57}$$

One benefit presented within [106] is the possibility to model the invariant and not only . Therefore, when modeling the normal force, and hence friction force, can also be considered, as shown in [106]. According to Toll and Månson [108] the average normal force within a volume is given by

$$\|\|F\_{l}^{\rm n}\|\| = \frac{\left\|F\_{l}^{\rm node}\right\|}{n\_{\rm node}d\_{\rm f}},\tag{3.58}$$

with the node force node and the number of nodes per volume node defined as

$$n\_{\rm node} = \frac{4\phi\_{\rm f}}{\pi d\_{\rm f}^2 \overline{a\_{\rm f}}'} \tag{3.59}$$

where ̅<sup>f</sup> is the average node space. One node is defined by four fibers with three contact points, as schematically shown in Figure 3.4. By assuming the contact points to be homogenously distributed along the fibers and according to Figure 3.4 the average node space is given by

$$
\overline{a\_{\rm f}} = \frac{2L\_{\rm f}}{N\_{\rm fc}} = \frac{2L\_{\rm f}}{\phi\_{\rm f} N\_{\rm fc} \phi},
\tag{3.60}
$$

with the specific number of contacts fc (cf. Eq. (3.47)) defined as

$$N\_{\rm fc\phi} = \,^\mathsf{B}\!/\_\mathsf{T} r\_\mathsf{f} f + 4 \,(g+1) = \frac{N\_{\rm fc}}{\phi\_\mathsf{f}}.\tag{3.61}$$

According to [108], node is given by

$$\left\| \left| F\_{l}^{\rm node} \right| \right\| = \frac{4}{\pi} \int\_{0}^{\Phi\_{\rm f}} \frac{1}{\overline{a\_{\rm f}} \overline{\mathbf{S}\_{\rm n}}} \mathrm{d}\phi\_{\rm f} \, ^{\prime} , \tag{3.62}$$

where <sup>n</sup> ̅̅̅ is the average nodal compliance, defined by

$$
\overline{S\_{\rm n}}^{\cdots} = \frac{2\overline{a}\_{\rm f}^{\overline{3}}}{\pi E\_{\rm f} d\_{\rm f}^4}.\tag{3.63}
$$

Figure 3.4: Schematic illustration of one node as defined in [108], consisting of four fibers and with node space . The contact points are assumed to be homogeneously distributed along the fiber

Solving Eq. (3.62) by considering Eq. (3.60) and Eq. (3.63) leads to

$$\begin{split} \left\| \left| F\_{\mathrm{f}}^{\mathrm{node}} \right| \right\| &= \int \frac{E\_{\mathrm{f}} d\_{\mathrm{f}}^{4} \phi\_{\mathrm{f}} \prime^{4} N\_{\mathrm{f} \mathrm{c} \phi}}{8 L\_{\mathrm{f}}^{4}} \mathrm{d} \phi\_{\mathrm{f}} \prime \\ &= \frac{E\_{\mathrm{f}} d\_{\mathrm{f}}^{4} \phi\_{\mathrm{f}} \prime^{5} N\_{\mathrm{f} \mathrm{c} \phi}\_{\mathrm{f} \mathrm{c} \phi}}{40 L\_{\mathrm{f}}^{4}} . \end{split} \tag{3.64}$$

Further combining Eqs. (3.58), (3.59) and (3.64) the average normal force is

$$\|\|F\_{\ell}^{\mathrm{n}}\|\| = \frac{\pi E\_{\mathrm{f}} d\_{\mathrm{f}}^{5} \phi\_{\mathrm{f}}^{3} N\_{\mathrm{f}\mathrm{c}\phi}^{3}}{80 L\_{\mathrm{f}}^{3}}.\tag{3.65}$$

For fc = 8f/, meaning neglecting the -term of the contact points,

which is the assumption within [108], Eq. (3.65) and Eq. (3.57) are identical. Within this work and similar to [106] Eq. (3.65) is used for normal and hence friction force.

Regarding the friction force, Coulomb friction is assumed in most cases and the force is given by

$$F\_l^{\rm fr} = k\_{\rm fr} \| F\_l^{\rm n} \| \| \[ \Delta U\_l^{\rm ff} \] \| \tag{3.66}$$

with fr being the friction coefficient and Δ ff is the relative velocity between the two fibers. Of course, the real relative velocity of the fibers is unknown due to the homogenization and the relative velocity Δ between matrix and fibers defined in Eq. (3.39), used for calculating the hydrodynamic force, is not valid in this case. The relative velocity of two fibers within one cell is given by

$$
\Delta U\_l^{\text{ff}} = D\_{lj} h\_l^{\text{ff}}, \tag{3.67}
$$

where ℎ ffis the distance vector between the fiber centers. Assuming the distance to be equal in every direction [106], Eq. (3.67) can be reformulated to

$$
\Delta U\_l^{\rm ff} = \sqrt{\frac{h\_j^{\rm ff} h\_j^{\rm ff}}{3}} \sum\_j D\_{lj} = \sqrt{\frac{h\_j^{\rm ff} h\_j^{\rm ff}}{3}} \Delta \mathcal{U}\_l^{\rm ff}. \tag{3.68}
$$

Since the friction force only depends on the normed relative velocity, ℎ ff does not need to be determined and the direction can be approximated by summing up the corresponding entries of the strain rate tensor.

Regarding the lubrication force, the literature contains several approaches with different level of detail and complexity. Yamane et al. [122], Bounoua et al. [123] or Meyer et al. [110] present detailed approaches depending on single fiber orientations and therefore not suitable for a macroscopic approach. Approaches suitable for a homogenized material are usually assumed to be linear proportional to a dimensionless coefficient lu. The linear dependence is further on the relative velocity Δ ff and matrix viscosity, as for example presented by Servais et al. [109]. Djalili-Moghaddam and Toll [124] further consider the surface area, but the area is assumed to be constant and the dependence is represented by <sup>f</sup> . Within this work and identical to [106], the lubrication force is assumed to be proportional to the relative velocity, matrix viscosity, overlapping area and reciprocal proportional to the fiber distance, since it vanishes with rising distances. Finally, the lubrication force is given by

$$\begin{split} F\_l^{\rm lu} &= \frac{k\_{\rm lu}}{\sqrt{h\_f^{\rm ff} h\_f^{\rm ff} / 3}} \eta\_{\rm M} A\_{\rm ff} \Delta U\_l^{\rm ff} \\ &= k\_{\rm lu} \eta\_{\rm M} A\_{\rm ff} \Delta \tilde{U}\_l^{\rm ff}, \end{split} \tag{3.69}$$

with ff being the average overlap area of two fibers within one cell.

#### **3.4.2.3 Approximation of fiber-fiber overlap area**

To approximate the fiber-fiber overlap area, it is assumed, that the fibers do not overlap near the ends [106]. Therefore, the overlap area is completely defined by the fiber diameter and the overlap angle . Two different cases, depending on orientation state must be separated. The first one is ω ≥ ωcrit, where the overlap area is a parallelogram, as shown in Figure 3.5a and Figure 3.5b. The second one is ω < ωcrit, where the area is a hexagon, as shown in Figure 3.5c.

If it is ω ≥ ωcrit the projected overlap area is given by

$$A\_{\rm ff} = \frac{d\_{\rm f}^2}{\sin(\omega)}.\tag{3.70}$$

For realistic aspect ratios (<sup>f</sup> > 10) and material orientation states, it is ω ≥ ωcrit in most processes. Nevertheless, the case ω < ωcrit can be reached, and must be considered. In this case, the assumption of small angles is an adequate approximation and the simplification sin(ω) ≈ ω is valid, hence the critical angle is given by

Figure 3.5: Two overlapping fibers with highlighted overlap area (green). Arbitrary angle > (a), critical angle = (b) and over-critical angle < (c). Red area is subtracted for calculation of overlap [106]

$$
\omega\_{\rm crit} = \frac{d\_{\rm f}}{L\_{\rm f}} = \frac{2}{r\_{\rm f}}.\tag{3.71}
$$

The overlap area is given by subtracting the not overlapping areas from the complete projected surface. Regarding Figure 3.5c, the red areas are subtracted from the complete area and result is the green area. So, for ω < ωcrit the overlap area is given by

$$A\_{\rm ff} = d\_{\rm f} L\_{\rm f} - \frac{L\_{\rm f}^2}{4} \omega. \tag{3.72}$$

Combining Eq. (3.70) and (3.72) the projected overlap area is given by

$$A\_{\rm ff} = \begin{cases} \frac{d\_{\rm f}^2}{\sin\{\omega\}} & \text{for} \quad \omega \ge \omega\_{crit} \\\frac{L\_{\rm f}^2}{\omega\_{\rm f} - \frac{L\_{\rm f}^2}{4}} \omega & \text{for} \quad \omega < \omega\_{crit} \end{cases} \tag{3.73}$$

The angle of the individual fibers is simply given by

$$
\omega = \cos^{-1} \{ q\_l^a q\_l^\beta \}. \tag{3.74}
$$

Unfortunately, the eigenvectors are always perpendicular to each other, so computing the individual angle would not be meaningful or representative. Furthermore, the average ω is unknown in the homogenized material, but as shown in [106] the averaged angle within one cell can be approximated with and is given by

$$
\omega \approx \frac{4f}{\pi}.\tag{3.75}
$$

Hence, the averaged, projected overlap area of the fibers within one cell can be determined depending on fiber geometry and orientation state, with information provided by the second order orientation tensor, fiber length and diameter.

#### **3.4.3 Constitutive fiber breakage model**

#### **3.4.3.1 Breaking criterium based on hydrodynamic forces**

As explained in Section 2.2.4.2, buckling is one of the major mechanisms for fiber breakage. In a first step, it is determined, whether the fibers within one cell are under pressure and buckling in general is possible. This state is reached if < 0 (Eq. (2.34)). Similar to calculation of hydrodynamic forces, explained in Section 3.4.1, the eigenvectors of the orientation tensor are used as reference fibers, and the Jeffery orbit is reformulated to

$$D\_{lj} \nu\_{ln} \nu\_{jn} < 0,\tag{3.76}$$

where, again the index indicates the number of the eigenvector and ∈ {1,2,3}. A summation of the corresponding eigenvalues approximates the amount of fibers, that are able to buckle, within one cell.

In a next step, it is checked, whether the force, acting on the fiber is high enough to cause buckling. This concerns only the eigenvectors fulfilling Eq. (3.76). Therefore, the hydrodynamic forces are used, as described in Section

3.4.1. For buckling, only the amount of force in fiber direction is relevant, which is given by the dot product of the force and the corresponding fiber direction. So, the relevant force is given by

$$
\hat{F}\_n^{\text{hydd}} = \nu\_{\text{ln}} F\_l^{\text{hydd}} \tag{3.77}
$$

and the buckling criterium is reached if

$$B\_{nm}^{\text{hydd}} = \frac{\mathcal{F}\_n^{\text{hydd}}}{F\_m^{\text{bu}}} \ge \mathbf{1},\tag{3.78}$$

with bu being the critical force for fibers with length f and defined in Eq. (2.36). Similar to [69] and [70] transverse forces, which favor buckling, are not considered, which is a huge simplification due to a lack of information in an macroscopic and homogenized approach.

#### **3.4.3.2 Fiber length evolution**

The evolution of fiber lengths is strongly based on the model of Phelps et al. [70] explained in Section 2.2.4.3. The breaking probability is defined as

$$P\_{nm}^{\text{br}} = \begin{cases} 0, & \text{caseA} \\ \mathcal{C}\_{\text{br}} \dot{\mathcal{Y}} \left( 1 - \exp\{1 - B\_{nm}^{\text{hyd}} \} \right), & \text{caseB'} \end{cases} \tag{3.79}$$

for every eigenvector and fiber length f . The cases are

$$\text{caseA} = B\_{nm}^{\text{hyd}} < 1 \quad \text{or} \quad D\_{lj} \nu\_{ln} \nu\_{ln} > 0 \tag{3.80}$$

and

$$\text{caseB} = \, B\_{nm}^{\text{hyd}} \ge 1 \quad \text{and} \quad D\_{lf} \nu\_{ln} \nu\_{ln} \le 0. \tag{3.81}$$

The final breakage probability for a specific length is weighted with the eigenvalues and given by

$$P\_m^{\text{br}} = \lambda\_n P\_{nm}^{\text{br}}.\tag{3.82}$$

The subsequent fiber breakage simulation is similar to Phelps et al. [70]. The child generation rate is defined similar to Eq. (2.40) with the same restrictions as described in Section 2.2.4.3 and [70]. The time evolution is given by Eq. (2.44), simulating the fiber length distribution by number of different possible fiber lengths in every cell.

### **3.4.3.3 Model restrictions and separation to state-of-the-art modeling**

Within Section 3.4.3.1 only the hydrodynamic forces from Section 3.4.1 are mentioned for fiber breakage, but not the friction and lubrication force explained in Section 3.4.2.2. Although these forces can be approximated, their influence on fiber breakage is unclear [68,70]. The fiber breakage is based on buckling, and therefore on pressure. The contact force may be set in relation to the direction of the eigenvectors and the contact points may be assumed equidistant along the fiber, but still it is unclear if the forces favor or hinder buckling and breakage. Additionally, the Jeffery orbit is only based on hydrodynamic forces and this aspect would not be valid any longer. Future studies have to focus on the role of contact forces in fiber breakage to enable a meaningful use in such simulations. Furthermore, the information may be used for a dynamic description of the breakage coefficient br, which is assumed to be constant to this point of time.

The main difference of the novel approach presented within this work and the state-of-the-art, is the calculation of the forces and the derived buckling criterium. In Phelps et al. [70] the forces are calculated by the approach of Dinh and Armstrong [72] (see Eq. (2.37)). In this thesis, the forces are determined by considering hydrodynamic drag and lift force, as described in Section 3.4.1. These forces act on reference fibers, represented by the eigenvectors of the orientation tensor.

Within [70], −2⟦⟧ = 1 is assumed, so the buckling criterium is independent of orientation and determined by matrix viscosity, shear rate, fiber geometry and a fitting parameter. Furthermore, the assumption −2⟦⟧ = 1 leads to an overestimation of br , since it is −2⟦⟧ < 1 for almost all possible orientations (cf. Eq. (2.38) and Figure 2.12). Additionally, this approach leads to a model where all fibers within one cell offer potential to break. The novel approach distinguishes between three reference fibers in every cell, where the Jeffery orbit and forces are calculated independent for every reference fiber. The final breaking probability is weighted with the corresponding eigenvalues. Consequently, this novel approach offers the following two advantages. One is, that the amount of fibers, able to break is determined more accurate. The other one is, that the forces are determined without material specific fitting parameters, so the breakage modeling needs one empirical parameter less compared to the state-of-the-art.

### **3.5 Modeling of warpage and residual stresses**

#### **3.5.1 Schematic procedure**

The simulation of warpage and residual stresses, needs information of the mold filling simulation and material state. Therefore, pre-processing algorithms, providing and transforming the relevant data are necessary. The general procedure of pre-processing is schematically shown in Figure 3.6.

The results of the mold filling simulation are mapped, including the fiber orientation tensor, for orientation averaging as well as the temperature and curing distribution, for the simulation itself. The input of material data, including information of the matrix in rubbery and glassy state, as well as mechanical and geometrical information of the fibers and the fiber volume content are used for the homogenization. The results of the homogenization are the effective mechanical and thermal properties as well as the coefficient vectors of thermal expansion and chemical shrinkage th and ch. The exact methods and models, which are used for the pre-processing and the simulation, are mentioned and explained in the following sections and a summary is given in Table 3.2 at the end of Section 3.5.2.3.

Figure 3.6: Schematic illustration of data transformation as pre-processing of structural analysis

### **3.5.2 Material model**

#### **3.5.2.1 Mechanical model**

To perform an adequate simulation of residual stresses and warpage, in a fiber reinforced part, the material model of the solid-mechanical structural simulation should be anisotropic. The material model for the fiber material is elastic, isotropic and independent of temperature, since glass fibers are simulated. For modeling the mechanical behavior of the matrix, the CHILE model, introduced in Section 2.3.1.2 is used. The anisotropy arises due to the homogenization of fibers and matrix to a homogeneous FRP. Such an approach was already successfully applied by Bernath [14] in case of continuous fiber reinforced epoxy parts, with a transversely isotropic behavior. Within this work, the model of Bernath is extended to represent orthotropic material behavior, being more meaningful for discontinuous reinforced injection molded parts, which are not fully aligned and hence not transversely isotropic.

The mechanical behavior of the matrix is given by Eq. (2.49). The elastic modulus in the glassy state is assumed to be constant. In the rubbery state, the elastic modulus is given by

$$E\_{\rm cr} = E\_{\rm cr,max} \left( \frac{c^2 - c\_{\rm g}^2}{1 - c\_{\rm g}^2} \right)^{8/3},\tag{3.83}$$

with cr,max representing the elastic modulus in the fully cured rubbery state (c = 1) and <sup>g</sup> being the point of gelation. This power law based approach is presented by Adolf and Martin [87] and also used by Bernath [14].

Besides the mechanical load, mainly based on holding pressure, the part undergoes chemical and thermal change, evoking thermal and chemical shrinkage. Hence, the complete change of strain is given by

$$
\Delta \varepsilon\_{lj} = \Delta \varepsilon\_{lj}^{\text{mech}} + \Delta \varepsilon\_{lj}^{\text{th}} + \Delta \varepsilon\_{lj}^{\text{ch}}, \tag{3.84}
$$

with change of thermal strain Δ th being defined by

$$
\Delta \varepsilon\_{ll}^{\text{th}} = \mathcal{O}\_l^{\text{th}} \Delta T \tag{3.85}
$$

and change of chemical strain Δ ch being defined by

$$
\Delta \varepsilon\_{ll}^{\text{ch}} = \mathcal{o}\_l^{\text{ch}} \Delta c,\tag{3.86}
$$

where Δ ch = Δ th = 0 for ≠ .

For the homogenization of the mechanical properties the method presented by Tandon and Weng [96] is applied to create a transversely isotropic material, which is orientation averaged afterwards with the method presented by Advani and Tucker [51]. This procedure is similar to Hohberg [99]. Description of th and ch are given in Section 3.5.2.2 for thermal expansion and 3.5.2.3 for chemical shrinkage.

#### **3.5.2.2 Thermal model**

The thermal properties, which need to be considered for the warpage analysis are the specific heat capacity, thermal conductivity, coefficient of thermal expansion and glass transition temperature. The latter depends on the matrix material and, since thermoset materials are regarded, on the degree of cure. The cross-linking due to curing complicates the sliding of the polymer molecules and hence <sup>g</sup> rises with the degree of cure. A modeling approach for the relationship between <sup>g</sup> and the degree of cure is given by DiBenedetto [125], by

$$T\_{\rm g} = T\_{\rm g.0} + \frac{(T\_{\rm g.o} - T\_{\rm g.0})\kappa\_{\rm Tg.}c}{1 - (1 - \kappa\_{\rm Tg})c} \tag{3.87}$$

with g,0 and g,∞ being <sup>g</sup> for = 0 and = 1 and Tg as material specific dimensionless modeling parameter. This approach is well established and applied in several studies [14,126–128].

The effective thermal conductivity is determined according to the work of Clayton [32]. Contrary to the mold filling phase, the holding pressure and curing phase take up a longer time period, with a significantly higher amount of heat flow. Therefore, it is reasonable to model the thermal conductivity as anisotropic vector th .

Within a transversely isotropic FRP, the conductivity is only distinguished in fiber direction <sup>1</sup> th = <sup>∥</sup> th and perpendicular <sup>2</sup> th = <sup>3</sup> th = <sup>⊥</sup> th. In fiber direction the property is simply given by volume averaging so

$$
\lambda\_1^{\text{th}} = \Phi\_{\text{f}} \lambda\_{\text{th,f}} + (1 - \Phi\_{\text{f}}) \lambda\_{\text{th,M}} \tag{3.88}
$$

with th,f and th,M representing the thermal conductivity of fibers and matrix, which are assumed to be isotropic.

For the perpendicular direction, Clayton [32] extends the general form for dilute dispersion given by

$$\frac{\lambda\_2^{\rm th} - \lambda\_{\rm th,M}}{\lambda\_2^{\rm th} + \mu \lambda\_{\rm th,M}} = \phi\_{\rm f} \left( \frac{\lambda\_{\rm th,f} - \lambda\_{\rm th,M}}{\lambda\_{\rm th,f} + \mu \lambda\_{\rm th,M}} \right) \tag{3.89}$$

and the corresponding differentiation

$$\frac{\left(\mu + 1\right)\lambda\_{\text{th},\text{M}}}{\left(\lambda\_2^{\text{th}} + \mu\lambda\_{\text{th},\text{M}}\right)^2} \text{d}\lambda\_2^{\text{th}} = \left(\frac{\lambda\_{\text{th},\text{f}} - \lambda\_{\text{th},\text{M}}}{\lambda\_{\text{th},\text{f}} + \mu\lambda\_{\text{th},\text{M}}}\right) \text{d}\Phi\_{\text{f}}.\tag{3.90}$$

with the dimensionless shape factor µ. It is

$$\frac{d\lambda\_2^{\rm th}}{(1+\mu)\lambda\_2^{\rm th}} = \left(\frac{\lambda\_{\rm th,f} - \lambda\_2^{\rm th}}{\lambda\_{\rm th,f} + \mu\lambda\_{\rm th}}\right) \frac{d\phi\_{\rm f}}{1 - \phi\_{\rm f}}\tag{3.91}$$

assuming <sup>2</sup> th = th,M for small additions of <sup>f</sup> and 1/(1 − <sup>f</sup> ) as correction factor for this assumption. Further the integration from <sup>2</sup> th = th,M at <sup>f</sup> = 0 to <sup>2</sup> th = th,f at <sup>f</sup> = 1 leads to

$$1 - \Phi\_{\rm f} = \left(\frac{\lambda\_{\rm th,f} - \lambda\_2^{\rm th}}{\lambda\_{\rm th,f} + \lambda\_{\rm th,M}}\right) \left(\frac{\lambda\_{\rm th,M}}{\lambda\_2^{\rm th}}\right)^{\frac{1}{1+\mu}}.\tag{3.92}$$

The shape factor is given by µ = 1 for cylindrical inclusions, such as fibers [32], so Eq. (3.92) leads to

$$\frac{\lambda\_2^{\rm th}}{\lambda\_{\rm th,M}} = \left(\frac{\lambda\_{\rm th,f}}{\lambda\_{\rm th,M}}\right) (1 - \Phi\_{\rm f}) \left(\frac{\lambda\_{\rm th,f}}{\lambda\_{\rm th,M}} - 1\right) \left(\frac{\lambda\_2^{\rm th}}{\lambda\_{\rm th,M}}\right)^{\frac{1}{2}}.\tag{3.93}$$

Finally, the thermal conductivity perpendicular to the fiber direction is given by

$$\begin{split} \lambda\_{2}^{\text{th}} &= \\ \frac{\lambda\_{\text{th,M}}}{4} \Biggl( \sqrt{(1-\Phi\_{l})^{2} \left(\frac{\lambda\_{\text{th,f}}}{\lambda\_{\text{th,M}}}-1\right)^{2} + 4\frac{\lambda\_{\text{th,f}}}{\lambda\_{\text{th,M}}}} \\ &- (1-\Phi\_{l}) \left(\frac{\lambda\_{\text{th,f}}}{\lambda\_{\text{th,M}}}-1\right) \Biggr)^{2} .\end{split} \tag{3.94}$$

The last property needed for warpage analysis is the thermal expansion coefficient. Here, the energy based approach of Schapery [97] is used. Based on the energy potential, bounds are defined for a linear elastic surface traction problem with space-wise constant temperature change. Within [97] it is shown, that the upper bound shows the minimal error for volume averaging, so

$$\mathfrak{d}\_{1}^{\text{th}} = \frac{\Phi\_{\text{f}} E\_{\text{f}} \theta\_{\text{th,f}} + (1 - \Phi\_{\text{f}}) E\_{\text{M}} \vartheta\_{\text{th,M}}}{\Phi\_{\text{f}} E\_{\text{f}} + (1 - \Phi\_{\text{f}}) E\_{\text{M}}},\tag{3.95}$$

where <sup>M</sup> is the elastic modulus of the matrix. The same procedure for the perpendicular direction results in

$$\begin{aligned} \boldsymbol{\vartheta}\_{2}^{\rm th} = \boldsymbol{\vartheta}\_{3}^{\rm th} = (1 + \nu\_{\rm f}) \boldsymbol{\Phi}\_{\rm f} \boldsymbol{\vartheta}\_{\rm th, f} &\quad + (1 + \nu\_{\rm M}) \boldsymbol{\Phi}\_{\rm M} \boldsymbol{\vartheta}\_{\rm th, M} \\ &\quad - \overline{\boldsymbol{\vartheta}\_{\rm th}^{\rm th} \vec{\nu}} \end{aligned} \tag{3.96}$$

with <sup>f</sup> and <sup>M</sup> being the Poisson ratio of fibers and matrix. Similar to the mechanical properties, the thermal conductivity and expansion coefficient are only determined for a transversely isotropic material. Therefore, the properties are further orientation averaged to create an orthotropic material with the orientation averaging scheme presented by Advani and Tucker [51] and explained in Section 2.3.2.2.

#### **3.5.2.3 Model for chemical shrinkage**

The homogenization of the chemical shrinkage is not implemented within homogenization the algorithm. Therefore, it is determined by simple volume averaging as

$$
\boldsymbol{\theta}\_{l}^{\text{ch}} = \boldsymbol{\Phi}\_{\text{f}} \boldsymbol{\theta}\_{l}^{\text{ch}, \text{f}} + (1 - \boldsymbol{\Phi}\_{\text{f}}) \boldsymbol{\theta}\_{l}^{\text{ch}, \text{M}},\tag{3.97}
$$

with ch,,f and ch,,M being the coefficients of chemical shrinkage for pure fibers and pure matrix. Since there is no chemical shrinkage in the fibers ( ch,,f = 0), Eq. (3.97) simplifies to

$$
\partial\_l^{\text{ch}} = (1 - \Phi\_{\text{f}}) \partial\_l^{\text{ch}, \mathcal{M}}.\tag{3.98}
$$

The chemical shrinkage is assumed to be isotropic, so

$$\mathfrak{v}\_{l}^{\text{ch}} = \mathbb{1}\_{\mathfrak{J}} \left( \left\| \mathfrak{v}\_{l}^{\text{ch}} \right\| \quad \left\| \mathfrak{v}\_{l}^{\text{ch}} \right\| \quad \left\| \mathfrak{v}\_{l}^{\text{ch}} \right\| \right). \tag{3.99}$$

Table 3.2 gives a summary and overview of the homogenized mechanical and thermal material properties with corresponding homogenization method and whether it is orientation averaged or not.

Table 3.2: Summary of homogenized material properties with corresponding homogenization methods


# **4 Verification and validation**

### **4.1 Numerical verification**

### **4.1.1 Phase-dependent boundary condition**

To verify the boundary condition described in Section 3.2.2 a stair-like model, as shown in Figure 4.1, is simulated, so different amounts of material reach the outlets at different points of time. At the beginning, the cavity is filled with air and the FRP enters the model at the inlet (blue) with a constant velocity of = (8.0 ∙ 10−3 0 0) m/s. Hence the cavity should be completely filled after 5 s, if no material leaves the model. The model is meshed with cubic hexahedrons, having an edge length of 1 mm.

Figure 4.1: Stair model for verification of phase-dependent boundary condition. Inlet in blue and three outlets with different distances to the inlet in red. Walls in transparent grey

To keep the simulation as simple as possible, the material is modeled isotropic and Newtonian with = 10 Pa∙s and a no slip boundary condition is chosen for the walls. Four different formulations of the phase-dependent boundary condition at the outlet are compared, as shown in Figure 4.2.

Figure 4.2: Simulation results after 5 s for different formulations of the phase-dependent boundary condition. Cut through the <sup>2</sup> -plane of symmetry (see Figure 4.1). Colors visualize the value of the VOF factor . Switchover means change from zeroGradient to (Eq. (3.3))

The case 'ZeroGradient' applies Eq. (3.2) all the time and independent of VOF factor , so the FRP can leave the cavity without any resistance. 'Switchover for > 0' means a spontaneous transition from *zeroGradient* to zero (Eq. (3.3)) as soon as > 0 without any intermediate state. Similar 'Switchover for = 1' means a spontaneous transition from *zeroGradient* to zero as soon as = 1. At last the interpolated formulation, as given by Eq. (3.4), is regarded. Figure 4.2 shows the comparison of the simulations at the 2 -plane of symmetry after 5 s, so the cavities should be completely filled with FRP. The blue areas for 'ZeroGradient' and 'Switchover for = 1' indicate a relative high material loss compared to the other two formulations, where the cavity is almost completely filled. Since all finite volumes of the

simulation mesh have identical dimensions, the results can be quantified by summing up all values of and divide the result by the total number of finite volumes of the domain. The results would be one for a completely filled cavity. The results are shown in Table 4.1.

Table 4.1: Normalized filled volume for different formulations of the phase-dependent boundary condition at the outlet


As expected, 'ZeroGradient' generates the highest material loss. The results of 'Switchover for > 0' and Eq. (3.4) are quite similar, although the interpolated boundary condition is slightly better. A reason for the similarity is the sharp flow front with only one or two partially filled elements along the flow path. A more complex material behavior resulting in a torn and chaotic flow front would lead to more partially filled elements at the flow front and hence to more remaining air in the cavity for 'Switchover for > 0'. Nevertheless, Eq. (3.4) generates the best results and is chosen for further mold filling simulations.

### **4.1.2 Analytical verification of isotropic flow modeling with anisotropic viscosity tensor**

The isotropic solution of is given by an isotropic distribution of fibers, represented by = 1/3 and = 0 for ≠ (see Figure 2.11), resulting in

$$
\eta\_{ljkl}^{\text{FRJ,iso}} = \begin{bmatrix}
\eta\_{\parallel\parallel}^{\text{iso}} & \eta\_{\parallel\perp\perp}^{\text{iso}} & \eta\_{\parallel\perp\perp}^{\text{iso}} & 0 & 0 & 0 \\
& \eta\_{\parallel\parallel}^{\text{iso}} & \eta\_{\parallel\perp\perp}^{\text{iso}} & 0 & 0 & 0 \\
& & \eta\_{\parallel\parallel}^{\text{iso}} & 0 & 0 & 0 \\
& & & \eta\_{\perp\perp}^{\text{iso}} & 0 & 0 \\
& & & & \eta\_{\perp}^{\text{iso}} & 0 \\
& & & & & \eta\_{\perp}^{\text{iso}}
\end{bmatrix}. \tag{4.1}
$$

The values for || iso , ||⊥ iso and <sup>⊥</sup> iso are given in Table 4.2.

Table 4.2: Values for the entries of , (Eq. (4.1)), calculated with Eq. (3.21)


Here, is determined with the linear closure approximation, representing the exact solution for an isotropic fiber distribution [62]. For an isotropic and Newtonian fluid, it would be <sup>12</sup> = <sup>23</sup> = and <sup>11</sup> = 3, known as the Trouton ratio as relation between shear and elongation viscosity, first described by Trouton in 1906 [129]. Assuming this, ,iso can be rewritten as

$$
\eta\_{ijkl}^{\text{FRP,iso}} = \eta \begin{bmatrix}
4/3 & -2/3 & -2/3 & 0 & 0 & 0 \\
& 4/3 & -2/3 & 0 & 0 & 0 \\
& & 4/3 & 0 & 0 & 0 \\
& & & 1 & 0 & 0 \\
& & & & 1 & 0 \\
& & & & & 1
\end{bmatrix},
\tag{4.2}
$$

being the correct fourth order viscosity tensor of an isotropic and Newtonian

fluid. Nevertheless, with respect to Eqs. (3.16)-(3.18) the only two solutions for <sup>12</sup> = <sup>23</sup> are Φ<sup>f</sup> = 0 and Φ<sup>f</sup> ≈ 3.383Φmax. Since the latter is unrealistic, Φ<sup>f</sup> = 0 is chosen to compare the anisotropic tensor to an isotropic solution, therefore, it is <sup>11</sup> = 0. Of course, this is also unrealistic, since the assumption and micro-mechanical models in Section 3.3.1.2 are not valid for Φ<sup>f</sup> = 0, but it is the only way to compare the approach to one of the analytical approaches mentioned in Section 2.2.1.2. By assuming Φ<sup>f</sup> = 0 it is <sup>12</sup> = <sup>23</sup> = <sup>M</sup> and <sup>11</sup> = 0. Hence the scalar effective viscosity according to Eq. (3.26) is given by eff,iso <sup>=</sup> <sup>4</sup> 5 ⁄ <sup>M</sup> and the complete fourth order tensor is given by

$$
\eta\_{ijkl}^{\text{FRP,iso}} = {}^4\!/\_5\eta\_{\text{M}} \begin{bmatrix}
4/3 & -2/3 & -2/3 & 0 & 0 & 0 \\
& 4/3 & -2/3 & 0 & 0 & 0 \\
& & 4/3 & 0 & 0 & 0 \\
& & & 1 & 0 & 0 \\
& & & & 1 & 0 \\
& & & & & 1
\end{bmatrix} .\tag{4.3}
$$

Hence, if implemented correctly, the simulation creates similar results to an analytical approach, calculated with a viscosity being 80 % of the viscosity in the simulations. As reference case, a Poseuille flow as described in Section 2.2.1.2 is chosen. The parameters are given in Table 4.3.

Table 4.3: Parameters for comparison of the viscosity tensor to a Poiseuille flow


The simulation model is built up in 3D. The boundary conditions are defined to fit the assumptions of a Poiseuille flow. The fiber aspect ratio is set to one to suppress fiber re-orientation. The mesh is created with 21 hexahedral elements along the diameter and 200 along the flow path. Since the simulation method is created for injection molding simulation, it is a transient solver. The simulation time is set to 10 s, according to the mean volume flow, the tube should be 50 times completely flown through at this point of time and the steady state is reached.

The results are shown in Figure 4.3, represented by the velocity profile along the diameter. The simulation results are slightly lower than the analytical ones, but fit well. The maximum velocity in the center is 1.5625 m/s for the analytical and 1.4991 m/s for the simulation approach, giving a maximum deviation of about 4 %. Regarding the complete profile the mean square error (MSE) is 0.0026, which represents a good agreement. One reason for the difference is the numerical solving and discretization in the simulation. Additionally, the analytical approach solves only a 1D problem, while the simulation is still solved in full 3D. 

Figure 4.3: Comparison of velocity profile between analytical approach (black) and simulation (red) in case of a Poiseuille flow

The analytical solution should not be seen as correct, since it also underlies assumptions and simplification, which are not all identical to the ones of the simulation. One example is the assumption of no radial flow in the analytical solution, which is not suppressed in the 3D-simulation model, although the velocity components are very low, they are not exactly zero. In summary, the simulation creates meaningful results, close to the analytical approach. This verifies the implementation in case of pure flow modeling, independent of anisotropic effects.

#### **4.1.3 Numerical verification of modeled forces and interaction points 4**

#### **4.1.3.1 Creation of individual fibers**

The approaches for calculating fiber forces, contact points and angles presented in Sections 3.4.1 and 3.4.2 are all based on orientations of individual fibers and modified to apply for a homogenized material with no information of individual fibers. Therefore, it is meaningful to verify these approaches by comparing them to results, calculated with individual fibers. The second and fourth order orientation tensors, needed for verification, are directly computed with individual fibers with Eq. (2.27) and Eq. (2.28).

For verification, 22 different orientation states (OS), consisting of 500 individual fibers each, are considered. The OS of a single fiber is defined by the two angles and as illustrated in Figure 2.10. In OS 1, the possible values are ∈ [0, ] and cos() ∈ [0,1] representing quasi-isotropic orientation. For OS with higher numbers, the window with possible values is reduced by /20 for and 1/20 for cos(), so the fibers become more highly orientated towards the <sup>1</sup> -direction. Due to this rate of change, full orientation is given for OS 21, represented by = 0 and cos() = 0. OS 22 represents planar isotropic orientation, so ∈ [0, ] and cos() = 0. Although planar

<sup>4</sup> Similar numerical studies and results are published in [106].

isotropic orientation is not typical in RIM processes, this OS is often used for verifications.

The possible values of and for the different orientation states are shown in Figure 4.4.

Figure 4.4: Possible values of and for different orientation states illustrated by colored areas (cf. Figure 2.10)

Due the definition of changing the value for cos() in a linear way, the change of is non-linear. By this definition the amount of fibers in <sup>3</sup> direction decreases faster, compared to the other directions <sup>1</sup> and <sup>2</sup> . This is also the case in many real fiber orientations in injection molding, due to the thin walled part design. Figure 4.5 shows three examples of the OS 1, 10 and 17 with corresponding orientation tensors. The single fiber orientations are visualized by the red dots on a unity sphere. The different fibers are created by using the 'rand()' function of Matlab R2019, creating a random value

between zero and one, which is multiplied with the maximum value of and cos() in the corresponding OS.

Figure 4.5: Visualization of different orientation states (red dots) on a unity sphere with corresponding orientation state number and second order orientation tensor

#### **4.1.3.2 Comparison of hydrodynamic forces for single fibers and homogenized material**

The hydrodynamic drag and lift force are calculated with Eq. (3.33), as explained in Section 3.4.1.1. For calculations with individual fibers ( hyd\_f ) the relation between orientation and relative velocity is given by . In case of a homogenized material ( hyd\_EV), the eigenvectors of are used as described in Section 3.4.1.3. Since 500 individual fibers exist for every OS, hyd\_f is calculated 500 times for each OS, while hyd\_EV is calculated three times, due to the three eigenvectors. To verify the novel approach, the average force of each OS is compared, being defined as

$$F\_{l}^{\text{hydd\\_av\\_f}} = \frac{1}{500} \sum\_{n=1}^{500} F\_{ln}^{\text{hydd\\_av\\_f}} \tag{4.4}$$

for the individual fibers. The averaged force of the homogenized approach is weighted with the eigenvalues so

$$F\_l^{\text{hydd}\_\text{av}\text{EV}\_\text{-EV}} = \sum\_{n=1}^3 \lambda\_n F\_{ln}^{\text{hydd}\_\text{-av\\_EV}}.\tag{4.5}$$

The comparison of the individual and homogenized average force is shown in Figure 4.6, split into the drag and lift part. Two relative velocities are regarded, the first one is chosen to be U<sup>i</sup> = 1/√3(1 1 1) m/s (Figure 4.6a) and the second one is U<sup>i</sup> = (1 0 0) m/s (Figure 4.6b). The material parameters are <sup>f</sup> = 0.015 mm, <sup>f</sup> = 100 and <sup>M</sup> = 20 Pa ∙ s.

The drag force components in Figure 4.6a are all identical, since the direction weighting of the drag force is given by , whose entries are identical. The lift force is zero for a quasi-isotropic orientation (OS 1), because the average fiber direction and are the same. Of course, there are lift forces, acting on the single and reference fibers, but they cancel each other during the summation for the average force. This is also the fact for planar quasi-isotropic orientation (OS 22), where <sup>1</sup> lift\_av ≈ <sup>2</sup> lift\_av ≈ 0 and only <sup>3</sup> lift\_av exists.

Figure 4.6: Comparison of average drag and lift force for individual fibers (f) and homogenized material (based on eigenvalues, EV) with different relative velocities = 1/√3 (1 1 1) / (a) and = (1 0 0) / (b)

In detail, the force components are not exact zero, since the fibers are created randomly, so the orientations are not exact (planar) quasi-isotropic.

In Figure 4.6b it is <sup>2</sup> = <sup>3</sup> = 0, hence <sup>2</sup> drag\_av <sup>=</sup> <sup>3</sup> drag\_av <sup>=</sup> <sup>0</sup>. <sup>1</sup> drag\_av decreases with higher degree of orientation since the projected area of the fibers in the <sup>1</sup> -plane decreases, being relevant for this and is represented by the change of d, as explained in Section 3.4.1.2. The lift force is zero for (planar) quasi-isotropic and full alignment, since these orientations are symmetric to . Hence <sup>3</sup> lift\_av = 0 for all orientations, since the orientation in <sup>3</sup> -direction is always symmetric to in this case. Furthermore, it is 1 lift\_av = 0, because the lift force acts perpendicular to .

In general, the homogenized results are in exact agreement with the individual solution as shown within [106], verifying the novel approach to determine the average force on an arbitrary orientation state. Nevertheless, an inference on single fiber forces is not possible, being a disadvantage since forces on individual fibers may exist, but cancel each other while calculating the average. Such phenomena require an approach on microscopic scale and modeling of individual fibers. However, the average force within one cell is well approximated.

#### **4.1.3.3 Comparison of contact points for single fibers and homogenized material**

Three approaches are compared to the individual fiber solution in case of fiber-fiber contact points. One is the interaction tensor as explained in Eqs. (3.48) - (3.51) and [120], representing the state of the art. The other two are the eigenvector based (Eq. (3.54) and Eq. (3.55)) and the polynomial fit (Eq. (3.56)) with the values given in Table 3.1 [106]. For all approaches the contact points are determined with Eq. (3.47), the difference is in the calculation of the invariants and of the fiber orientation distribution. For individual fibers Eq. (3.44) and Eq. (3.45) are used to determine and . The creation of the individual fibers and different orientation states is described in Section 4.1.3.1 and the second and fourth order orientation tensors are directly computed with the individual fibers.

low,×) and polynomial approximation (blue,+) to solution with individual fibers (green,□) and interactions tensor based 

Figure 4.7: Values of the fiber orientation invariants

method presented in [120]

 (red,○)

(a) and 

(b) for the different approaches. Comparison of eigenvector approach (yel-

In a first step, the results of and are compared for the 22 different orientation states, shown in Figure 4.7. Due to the correction factor (see Section 3.4.2.1) the interaction tensor and eigenvector approach fit exactly, which is also shown analytically within [106] for specific cases. Furthermore, both approaches predict a too low for higher OS, while the polynomial approach fits well besides OS 21. The maximum deviation of the eigenvector and interaction tensor approach is about 90 % at OS 20. Besides full orientation (OS 21), which is an unrealistic OS and will never be reached in the real process, the polynomial approach has a maximum deviation of about 8 % for OS 10. The MSE of the eigenvector and interactions tensor approach is 0.01 and 0.001 for the polynomial approach, which is 10 times smaller.

Regarding the approximation of (Figure 4.7b) the interaction tensor is no longer considered, since only an approximation for is given within [120]. Contrary to the eigenvector method shows a maximum deviation for quasiisotropic orientation of about 19 % and fits better for higher orientation. Again, the eigenvector approach fits perfectly for full alignment. The polynomial fit approach shows again good agreement for all OS and has maximum deviation of about 2.4 % for planar isotropic orientation (OS 22). The MSE of the polynomial approach is about 0.000035 and therefore about 484 times lower than the MSE of the eigenvector approach, being 0.015. In general, the polynomial fit shows the best results with the lowest MSEs and should also create the best results in case of calculating contact points.

For comparison of contact points per fiber, according to Eq. (3.47), a theoretical material with <sup>f</sup> = 0.35 and the aspect ratios <sup>f</sup> = 10 (Figure 4.8a) and <sup>f</sup> = 100 (Figure 4.8b) is regarded [106]. As expected, the polynomial approach creates the best fitting results, with a maximum deviation of 29.4 % for <sup>f</sup> = 10 and 292 % for <sup>f</sup> = 100 in the unrealistic case of full alignment.

The results of the eigenvector approach are slightly better than the interaction tensor for all OS, due to the consideration of , creating an offset of exactly 4f( +1) (see Eq. (3.47)).

(b)

The relative importance of this offset becomes smaller for higher aspect ratios, since the other term to determine the contact points is linear in and f , which is the reason why 4f( + 1) is often neglected. The maximum deviation of the eigenvector approach is about 29.8 % for <sup>f</sup> = 10 and 65 % for <sup>f</sup> = 100 at OS 17 and OS 19.

The MSEs of the polynomial and eigenvector approach are 0.0917 and 1.013 for <sup>f</sup> = 10 and 8.83 and 82.82 for <sup>f</sup> = 100, so the MSE of the polynomial approach is about 10 times smaller. In general, and as expected, the polynomial approach predicts the contact points most accurately, with respect to the other approaches.

#### **4.1.3.4 Comparison of average fiber-fiber angle and overlap area for single fibers and homogenized material**

In the next step, the approach to approximate the average fiber-fiber angle presented in [106] and Eq. (3.75) is verified. Again the 22 different orientation states described in Section 4.1.3.1 are used. Firstly, the individual angles between the 500 single fibers are determined with Eq. (3.74), resulting in 124500 angles, since one angle always includes two fibers and a fiber has no angles with itself (500<sup>2</sup> /2− 500 = 124500). Afterwards, the average of these 124500 angles is calculated and compared to the -based approximation given in Eq. (3.75). The value of is computed directly with the individual fibers (Eq. (3.44)) since the aim is to verify Eq. (3.75) and therefore should be as exact as possible and no approximation should blur the results.

The comparison is given in Figure 4.9, showing good agreement. Due to the fitting factor 4/, there is a good agreement for (nearly) quasi-isotropic orientation states. The maximum deviation is about 25.8 % for OS 20 and the MSE is 8.45. The angle deviates about 2.25° on average. Due to the good agreement, the results verify Eq. (3.75) to be well suited to approximate the average angle within an arbitrary fiber distribution. Hence, the next step is to verify the modeling of the average fiber-fiber overlap area, which is based on this angle approximation. The comparison of the average fiber-fiber overlap area is given in Figure 4.10 for the aspect ratios of <sup>f</sup> = 10 (Figure 4.10a) and <sup>f</sup> = 100 (Figure 4.10b). For the reference case, 124500 overlap areas,

according to the 124500 angles are determined with Eq. (3.73) and averaged for every orientation state. 

Figure 4.9: Computed average fiber-fiber angle by averaging single angles of individual fibers (black, □) and approximation based on f given by Eq. (3.75) (green, )

Since ∈ sym., the eigenvectors are always perpendicular to each other and computing the individual overlap is not meaningful for the approximations. Therefore Eq. (3.75) is used to determine the average angle and afterwards the average overlap with Eq. (3.73). Only the eigenvector and polynomial approach are compared, since the average angle only depends on , the results of the eigenvector and interaction tensor based approaches would be identical.

Both approximations predict a slightly too low overlap for both aspect ratios in case of lower orientation states (OS < 15). In higher orientated states, the eigenvector approach overshoots the results, while the polynomial approach fits well, except for the unrealistic state of full alignment.

Here, the polynomial approach shows its maximum deviation with about 27.5 % for <sup>f</sup> = 10 and 870 % for <sup>f</sup> = 100, while the eigenvector approach fits perfectly. Again, this case will not be reached in a real process and therefore will not be included in the further discussion. The eigenvector approach has a maximum deviation of 64.5 % for <sup>f</sup> = 10 (OS 17) and 162 % for <sup>f</sup> = 100 (OS 20). The MSEs are 1.48 ∙ 10−20 for <sup>f</sup> = 10 and 3.36 ∙ 10−19 for <sup>f</sup> = 100 in case of the polynomial approach and 6.1 ∙ 10−20 for <sup>f</sup> = 10 and 9.1 ∙ 10−18 for <sup>f</sup> = 100 in case of the eigenvector approach. The MSE of the polynomial approach is about 4.1 times smaller for <sup>f</sup> = 10 and about 27 times smaller for <sup>f</sup> = 100. Again, the polynomial fit approximation creates the best results with the smallest MSE and should be used for further simulations, since the approximation of the lubrication force will be the most adequate, compared to the other mentioned approximations.

### **4.1.4 Mesh study**

The interpolation of significant different material properties of FRP and air with the VOF in the flow front region will always create different results for different meshes. In general, simulations with VOF will never be independent of the mesh [28]. Nevertheless, the flow front shows identical tendencies, when reaching a corresponding mesh refinement level, such as in a mesh study presented in [26]. Here, the flow front of a standard isotropic Newtonian and the anisotropic non-Newtonian approach are compared for different mesh refinements. The simulated geometry is a 3 mm thick square plate with 100 mm edge length and a circular inlet with a 15.5 mm diameter, as shown in Figure 4.11. The model is meshed with hexahedral elements with 1 mm edge length in <sup>1</sup> - and <sup>2</sup> -direction and 5, 10 and 15 elements in <sup>3</sup> -direction (plate thickness).

The cavity is filled with a constant volume flow of 50 cm³/s, always perpendicular to the inlet surface. Similar to [26], two different flow modeling cases are compared. One is isotropic and Newtonian with = 100 Pa∙s, the other one is with anisotropic viscosity tensor (Eq. (3.21)) and the Castro-Macosko model (Eq. (2.19)) for matrix viscosity and Kamal-Malkin model (Eq. (2.11)) for curing kinetics . The model parameters are given in Table 4.4 and Table 4.5.

Figure 4.11: CAD-model and geometry information for the mesh study. Square plate with circular inlet in blue and outlet in red

The fiber orientation at the inlet is <sup>11</sup> = 0.8 at the walls and <sup>11</sup> = 0. 3̅ at the inlet center with parabolic interpolation. The further entries are <sup>22</sup> = <sup>33</sup> = (1 −11)/2 and <sup>12</sup> = <sup>13</sup> = <sup>23</sup> = 0, if not explicitly mentioned differently. The fiber aspect ratio is 75. The fiber orientation is determined using the RSC-model (Eq. (2.30)) with =0.03.


Table 4.4: Model parameters for Castro-Macosko matrix viscosity model [26]


Table 4.5: Parameters for the Kamal-Malkin kinetic model [25]

The simulation results of the six different cases are shown in Figure 4.12. A significant difference of the flow fronts between 5 and 10 elements can be detected in the isotropic and anisotropic case.

Figure 4.12: Simulated flow fronts for different simulation methods and different mesh sizes with increasing number of elements over plate thickness. Detail area of the complete geometry shown in Figure 4.11 [26]

Further refinement from 10 to 15 elements has almost no effect on the isotropic Newtonian results. In the anisotropic non-Newtonian case, the flow fronts are slightly different for 10 and 15 elements, but show the same tendencies, deviating from the sharp result for 5 elements.

Both flow fronts are torn and chaotic as expected in RIM process, even though the mesh is structured. The torn and chaotic area is in a similar scale for both cases. Therefore, a mesh with 10 elements over plate thickness seems to reach some kind of convergence state. Although, the exact formation of the flow front will always scatter for different meshes, due to minimal mesh inaccuracies and numerical discretization. According to [26] a minimum number of 10 elements is applied between two walls in the simulations.

#### **4.1.5 Verification of flow and force modeling 5**

In the following Sections 4.1.5.1 to 4.1.5.3 different configurations are simulated to highlight the benefits of anisotropic flow coupling [26] and show the newly gained information about the calculated forces [106]. Therefore, fiber aspect ratios and initial orientations are varied. The simulation model is identical in all cases, and similar to the model used in Section 4.1.4 (Figure 4.11) with 10 elements in <sup>3</sup> -direction. The material is injected with a constant volume flow of 50 cm<sup>3</sup> /s, perpendicular to the inlet's surface and a no slip boundary condition is applied for the walls. The material has temperature of 110 °C at the inlet and the tool temperature is 170 °C. The Castro-Macosko viscosity model (Eq. (2.19)) is used for the matrix viscosity and the Kamal-Malkin model for curing kinetics (Eq. (2.11)). The parameters are given in Table 4.4 and Table 4.5. The fiber aspect ratio is constant <sup>f</sup> = 20 or <sup>f</sup> = 100, depending on the simulation. Fiber orientation model and initial conditions are identical to Section 4.1.4.

<sup>5</sup> Similar numerical studies and results published in [26] and [106].

#### **4.1.5.1 Influence of fiber length**

Figure 4.13 compares two different fiber lengths with resulting aspect ratios of r = 20 and r = 100. Due to the anisotropic viscosity tensor, depending on fiber length, the flow front builds up in a different way. The higher aspect ratio leads to higher viscosities and hence to higher viscous stress, resulting in a higher gradient of viscosity along the <sup>3</sup> -direction. Consequently, the flow front in the simulation with longer fibers is more torn and has a larger area with partial wall contact. Such a flow front is a typical phenomenon of RIM processes as explained in Section 2.1.2.3 and would not be predicable within an isotropic, single-phase simulation.

Figure 4.13: Comparison of flow fronts at the wall for two different fiber lengths: = 20 (left) and = 100 (right) at = 0.2 (a) and = 0.3 (b)

#### **4.1.5.2 Influence of fiber orientation**

Within this Section, two different initial conditions of the fiber orientation state are compared. The first one is the one explained in Section 4.1.4. The second one is the reverse case, so it is <sup>11</sup> = 0. 3̅ at the walls and <sup>11</sup> = 0.8 at the inlet center. Both simulations are performed with r = 100.

Figure 4.14 shows the results for state of filling (Figure 4.14a) and cavity pressure (Figure 4.14b) after 0.5 s of filling. The different orientations have no major influence on the flow front evolution, but there is a significant influence on the in-mold pressure. Due to the high degree of orientation near the walls, the viscosity rises, leading to a higher pressure. In an isotropic standard simulation, the viscosity does not depend on the orientation state and therefore, no difference between the two initial states would be visible. This highlights a benefit of the anisotropic flow coupling, with focus on better pressure prediction especially near gates and other regions with significant change of fiber orientation.

Figure 4.14: Comparison of state of filling (a) and corresponding cavity pressure (b) at the wall for different initial fiber orientations. 80% fibers aligned at the wall (left) compared to 80% fibers aligned in the symmetry plane (right). The aspect ratio is = 100 in both cases

#### **4.1.5.3 Force distributions**

Figure 4.15 shows the <sup>11</sup> component results of the injection molding simulation with <sup>f</sup> = 100. The resulting average hydrodynamic force (Eq. (3.33)), friction force (Eq. (3.66)) and lubrication force (Eq. (3.69)) are shown in Figure 4.16, Figure 4.17 and Figure 4.18. The average hydrodynamic force (Figure 4.16) is quite independent of the fiber orientation. Due to the high relative velocity, the force is higher near the inlet, compared to the remainder part. Regarding a path along the <sup>3</sup> -direction, the force is low at the wall elements and in the core region, which is also traceable to the low relative velocity in these regions.

Figure 4.15: Injection molding simulation result of fiber orientation tensor component 11. Cut through the <sup>2</sup> -plane of symmetry (cf. Figure 4.11)

The comparably high forces near the wall, but not directly at the wall correspond to the regions where fiber breakage is observed in the real process. Of course, fiber damaging is also observed directly at the walls, but other phenomena such as fiber-wall interactions, which are not considered within this simulation are more important here.

Figure 4.16: Injection molding simulation result of average hydrodynamic force ‖ ℎ\_‖ on a single fiber according to fiber orientation shown in Figure 4.15. Cut through the <sup>2</sup> -plane of symmetry (cf. Figure 4.11)

Figure 4.17: Injection molding simulation result of average friction force on contact point ‖ ‖ according to fiber orientation shown in Figure 4.15. Cut through the <sup>2</sup> -plane of symmetry (cf. Figure 4.11)

Contrary to hydrodynamic and lubrication forces, the friction force (Figure 4.17) is not depending on an absolute value of any relative velocity. Hence it is not higher near the inlet. The friction force is depending on the specific number of contacts fc, which is higher in less orientated regions, as shown in Section 4.1.3.3. Therefore, the friction is higher in the core region, where the fibers are less orientated and lower at the walls, especially in regions with high degree of orientation (cf. Figure 4.15).

The lubrication force distribution, shown in Figure 4.18, is mainly depending on fiber-fiber overlap and velocity gradient. Hence, the high regions are regions with high degree of orientation (cf. Figure 4.15), due to the larger overlap area.

Figure 4.18: Injection molding simulation result of average lubrication force on contact point ‖ ‖ according to fiber orientation shown in Figure 4.15. Cut through the <sup>2</sup> plane of symmetry (cf. Figure 4.11)

Most high values are reached by a combination of high degree of orientation and high velocity gradients as detectable at walls near the inlet. Another region with slightly higher lubrication forces is near the walls and near the flow front, where the fiber orientation is almost isotropic, but the velocity gradient assumes high values.

In general, all three computed forces show meaningful results along the flow path with respect to fiber orientation and flow field. The friction and lubrication force show both a significant dependency on orientation. Therefore, an orientation dependent modeling of these forces is essential when using them in further modeling approaches. Knowledge about the mentioned forces may help to improve state of the art models of fiber orientation, breakage, fibermatrix separation or viscosity and flow modeling itself. The hydrodynamic forces, used to calculate fiber breakage, are experimentally validated in Section 4.2.3.

### **4.2 Experimental validation of filling simulation**

#### **4.2.1 Anisotropic flow modeling 6**

The novel approach for anisotropic flow modeling should be experimentally validated. The validation should include fiber materials with fiber dependent quantifiable flow effects. Therefore the experimental work of Chiba and Nakamura [130] is referred. Here, the steady state flow of a Newtonian fiber suspension in a 1:4 backward-facing step channel is determined for different Reynolds numbers (). Caused by the backward step, eddies built up for higher , influencing the flow field and fiber orientation. The suspension is a Newtonian water-syrup-solvent with <sup>M</sup> = 0.0463 Pa∙s, filled with Φ<sup>f</sup> = 2 ∙ 10−5 vinylon fibers, having an aspect ratio of r<sup>f</sup> = 283. Due to the low fiber volume content, the Folgar-Tucker model (Eq. (2.29)) is used for fiber orientation modeling. The simulation model with geometrical information is shown in Figure 4.19. Due to the large model, the mesh is created with cubic elements having an edge length of 2 mm. Hence there are areas with less than 10 elements between two walls, which is acceptable in this case, since the complete cavity is filled with suspension and no flow front is regarded.

<sup>6</sup> Results also published in [26].

Figure 4.20 represents the simulation results in case of streaming lines and fiber orientation ellipsoids, calculated with the eigenvalues and eigenvectors of the second order orientation tensor.

Figure 4.19: Simulation model for the experimental work of [130], 1:4 backward-facing step channel for creating a recirculating flow. Inlet in blue, outlet in red and green frame, representing the area of Figure 4.20

The results are obtained within the green frame in Figure 4.19 ±2 mm in x<sup>3</sup> direction, so some streamlines disappear, since they leave this control volume. The streamlines visualize no eddy for =1.9, but for =11.1 and =14.9. The vortex increases with . This effect and the dimensions of the eddies correspond to the numerical and experimental results of Chiba and Nakamura. Also similar to [130], the fibers align along the streaming lines with higher orientation near the walls and more randomly orientation in the vortexes. The good agreement validates the simulation approach in case of flow modeling and fiber orientation for low fiber volume content. [26]

#### **4.2.2 Fiber orientation 7**

#### **4.2.2.1 Material, models and process conditions**

For validation of fiber orientation, the work of Englich [9] is considered, containing experimental data of a 35 %-vol short glass fiber filled phenolic compound. The fiber orientation is determined with microscopic images, presented in [9]. The images are created of RIM experiments with a square plate, having an edge length of 150 mm, a thickness of 2 mm and a 135 mm long and 1.5 mm high line inlet along one edge. The complete cavity with sprue is given in Figure 4.21. The sprue is meshed with 27 hexahedral elements along the diameter and 115 in <sup>3</sup> -direction. The plate is meshed with 300 hexahedral elements in <sup>1</sup> - and <sup>2</sup> -direction and 10 in <sup>3</sup> -direction. The simulation is performed anisotropic, non-Newtonian and non-isothermal. The FRP is injected with 16 cm³/s and 110 °C into the mold with a constant surface temperature of 175 °C. Similar to [26] and Section 4.1.4, the Castro-Macosko viscosity model (Eq. (2.19)) is used for the matrix viscosity and the Kamal-Malkin model for curing kinetics (Eq. (2.11)), the parameters are given in Table 4.4 and Table 4.5.

The fiber orientation is determined using the RSC-model (Eq. (2.30)) with strain reduction factor =0.03 and =0.08 for a reference case with isotropic viscosity modeling. Both strain reduction factors are optimized to fit the orientation distribution [26]. The fiber orientation at the inlet is chosen to be <sup>33</sup> = 0.8 at the walls and <sup>33</sup> = 0. 3̅ at the inlet center with parabolic interpolation. The further entries are defined as <sup>11</sup> = <sup>22</sup> = (1 − 33)/ 2 and <sup>12</sup> = <sup>13</sup> = <sup>23</sup> = 0, being also similar to Section 4.1.4, but the fibers are orientated dominantly in <sup>3</sup> -direction, not <sup>1</sup> , due to the direction of the sprue. The fiber aspect ratio is constant <sup>f</sup> = 15 [26]. The isotropic reference case is also simulated with the Castro-Macosko viscosity model, the model parameters are given in [25]. The curing modeling of the isotropic case is identical to the anisotropic one, so the parameters are given in Table 4.5.

<sup>7</sup> Results also published in [26].

Figure 4.21: Square plate mold model for fiber orientation simulations, derived from [9]. Circular inlet in blue and outlet in red [26]

#### **4.2.2.2 Results for fiber orientation distribution**

Figure 4.22 shows the results, represented by distribution of <sup>11</sup> over plate thickness in the plate's center. 

Figure 4.22: Second-order orientation tensor component <sup>11</sup> over plate thickness for RSC-model with isotropic (red) and anisotropic flow modeling (green). Experimental result (black) derived from microscopic images in [9] [26]

The experimental results are given by the numerical mean with corresponding standard deviation. Both simulation approaches create good results validating the anisotropic approach in case of fiber orientation modeling. The anisotropic results fit better near walls, while the isotropic result is better in the core region. Near the walls (first two and last two points) the MSE is 0.0042 for the isotropic and 0.0016 for the anisotropic simulation. In the core region (three middle points) the mean square error is 0.0069 for the isotropic and 0.019 for the anisotropic simulation.

Regarding the complete region the mean square error is 0.0054 for the isotropic and 0.0091 for the anisotropic simulation. Although the results are quite similar, the strain reduction factor is more than 2.5 times higher in the isotropic simulation, highlighting the influence of fiber-flow coupling for fiber re-orientation.

### **4.2.3 Fiber length**

#### **4.2.3.1 Material, models and process conditions**

For the experimental investigations focusing on fiber length and cavity pressure within this study, an experimental long fiber phenolic molding compound of the novolac type is considered. The compound Porophen GF9201L12a by Sumitomo Bakelite Co., Ltd. was supplied as unidirectionally glass fiber reinforced flakes with an initial fiber length of f0 = 12 mm [131]. According to the measurements, the fibers are assumed with a constant diameter of <sup>f</sup> = 0.0017 mm. This compound is typically used for compression molding applications, not for injection molding, due to the processing difficulties. However, it is used for injection molding in a study to develop a long fiber RIM process in a DFG-project running since 2018, with project number 400343062. These experiments have been performed at the Fraunhofer ICT in 76327 Pfinztal, Germany [132].

Processing the experimental long fiber phenolic compound has proven to be difficult in the injection molding process, due to the strong adhesion of the resin to the screw and the barrel of the injection molding machine, resulting in step-like and abrupt screw movement during plasticizing. Previous research studies by other authors reported similar problems [133]. A higher back pressure and higher barrel temperatures had to be used to counteract the adhesion, compared to conventional short fiber molding compounds. The higher temperatures reduce the adhesion to the barrel, but of course also support the curing of the resin. In order to avoid premature curing, plasticizing was only performed with the plasticizing unit retracted from the hot mold [132]. The process parameters found in this way are given in Table 4.6. A stable and repeatable process was possible.

The injection molding experiments were performed on a KraussMaffei 550/2000 GX injection molding machine equipped with a standard 60 mm thermoset screw and no non-return valve. The temperature control of the plasticizing unit is realized by four individually controlled, oil-tempered zones. The clamping unit has a maximum clamping force of 5500 kN. For inmold pressure measurement, two capacitive sensors of the type 6163 manufactured by Kistler Instrumente GmbH (Sindelfingen, Germany) are placed within the mold. Plates with a constant wall thickness of four millimeters and a size of 480 mm × 190 mm are manufactured. The mold has a central sprue with a cone shape and 185 mm height, further defined by a start diameter of 9 mm, an end diameter of 15.5 mm. The simulation model of the plate with sprue and part of the screw chamber is shown in Figure 4.23.

In the simulation, the front part of the plasticizing unit barrel and the nozzle are considered (Figure 4.23b). By doing so, a more realistic modeling of fiber orientation and length distribution at the beginning of the sprue is possible, since the initial fiber length distribution is determined in the frontal part of the plasticizing unit barrel (see Section 4.2.3.2). The model is meshed with hexahedral cells. The cells in the plate have a length of 2.3 mm in <sup>1</sup> -, 1.75 mm in <sup>2</sup> - and 0.4 mm in <sup>3</sup> -direction. The sprue and nozzle are built up with 24 cells along the diameter, 32 along the circumference and 1.5 mm edge length in <sup>3</sup> -direction.


Table 4.6: Process parameters for injection molding experiments of the long fiber phenolic molding compound [132]

The screw chamber is meshed with 72 elements along the diameter, 64 along the circumference and 1.57 mm edge length in <sup>3</sup> -direction.

At the beginning of the simulation, the screw chamber is filled with FRP, while the rest of the cavity contains air. In the following, FRP enters the model at the top surface of the screw chamber, which is highlighted in blue in Figure 4.23. Both short edges of the plate are defined as outlet, but only one is visible in Figure 4.23a. The matrix viscosity, curing and fiber orientation model are identical to Section 4.2.2.1. Since the resin system is quite similar to the one considered in Section 4.2.2.1 the same parameters are used for the viscosity model, but the viscosity is lowered to 50 %, since the resin system of the long fiber material shows a lower viscosity with same tendencies on temperature and shearing.

Figure 4.23: Cavity for validation (rectangle plate with central sprue) with positions of pressure sensors <sup>1</sup> and <sup>2</sup> , area for fiber length evaluation and inlet area in blue and outlet in red. Complete cavity (a) and cut through <sup>2</sup> -plane of symmetry for detailed visualization of screw chamber and nozzle (b)

The thermal boundary conditions are of Dirichlet type and fit to the process boundary conditions (see Table 4.6). For the velocity, a no-slip boundary condition at the walls is assumed, which is a simplification in RIM simulation, but creates good results as mentioned in Section 2.2.1.5. At the inlet the velocity is inlet = (0 0 −0.425) m/s. For the velocity at the outlet, the phase-dependent boundary condition given by Eq. (3.4) is applied. The material has an initial cure of 0.01 % and the fiber orientation distribution at the inlet and in the initially stored material is quasi-isotropic. The fiber length distribution at the inlet and in the initially stored material is identical to the measurements in the screw chamber (see Section 4.2.3.2). All other boundary fields (fiber orientation at the walls, curing rate, etc.) are set to *zeroGradient*  at the walls*.* Only the mold filling with a time of 3.3 s is simulated, so the mold is filled. Since no relative velocity appears in the completely filled mold, no fiber breakage will occur in the simulation. Therefore, ongoing process steps like holding and curing have no effect on the fiber length distribution within the simulations.

#### **4.2.3.2 Results for fiber length distribution**

For the fiber length analysis, specimens are extracted from the frontal part of screw chamber and the molded part (position <sup>f</sup> Figure 4.23). The measurement is carried out by and according to the method described and validated by Maertens et al. [132]. The fiber length measurement principle is destructive and based on the commercially available FASEP system (IDM Systems, Darmstadt, Germany). The measurement process works in four steps according to the work of Goris et al. [134]. For the first step (matrix removal), a sample of approximately 2.5 g is extracted from the plate and the matrix is removed by for 18 h in pyrolysis at 650 °C and air atmosphere using a TGA701 device by LECO Corporation, (St. Joseph, Michigan, USA). Afterwards, the ash residue is transferred into an aqueous solution, which is diluted using an apparatus consisting of a stirrer and beaker. Therefore, a repeatable and controlled sample-taking out of the dilution is enabled. Contrary to current state-of-the-art fiber dispersion and sample-taking methods, the influence of the operator on the measurement results is reduced [132]. The third and fourth step are the image acquisition and the image processing (fiber detection) by using the algorithms provided by the FASEP software.

For the simulation results, the fiber breakage is simulated with the method described in Section 3.4.3. For the child generation rate and the length evolution, the state-of-the-art approaches given by Eq. (2.40) and Eq. (2.44) are used. The force is determined with eigenvectors and hydrodynamic forces as described in Section 3.4.1 and based on these forces, the breakage probability is determined with Eq. (3.79). 24 different fiber lengths between 0.25 mm and 11.75 mm (in 0.5 mm steps) are possible, corresponding to the clustering of the fiber length measurement, being also 0.5 mm. The breakage coefficient is set to be br = 0.025 and the standard deviation for breaking point is br = 1, being identical to the values used by Phelps et al. [70]. The difference is the force modeling, which is based on hydrodynamic forces and set in relation to the eigenvectors of the orientation tensor in this work (Eq. (3.33)).

The experimental results of the screw chamber (black) and plate (green), as well as the simulation results in the plate (blue) are shown in Figure 4.24 for fiber lengths from 5.25 mm to 11.75 mm and Figure 4.25 for 0.25 mm to 4.75 mm. The experimental results are the average of five measurements, illustrated with corresponding standard deviation and results are weighted by fiber length.

Only slightly more than 5 % of the initially 12 mm fibers remain at the end of plastification, meaning that most of the fiber breakage happens within the plastification and not during mold filling. There are some other longer fibers left with an amount of 2 % around 6 mm and 8 mm length, but most of the fibers (about 70 %) are already shorter than 1 mm before entering the mold. Regarding the plate measurement, nearly no fibers longer than 3.25 mm are left, since the amount is about 0.6 % in sum. Therefore, the amount of fibers shorter than 1 mm has risen to about 78.5 %.

Regarding the complete measurement of screw chamber and plate, the average fiber length in the experiments decreases by about 14 % from 0.434 mm to 0.38 mm. The simulation corresponds to the measurements, predicting about 79.6 % to be shorter than 1 mm, being a deviation of only 1.1 %. As shown in Figure 4.24 and Figure 4.25 the simulated fiber length distribution is in good agreement with the experimental data, showing a slight rise of proportions from 11.25 mm to 1.75 mm and a following faster rise of proportions up to 0.25 mm. 

Figure 4.24: Length weighted amount of fibers with lengths of 5.25 mm to 11.75 mm. Experimental results of screw chamber (black) and in plate (green) with corresponding standard deviation. Simulation results in the plate in blue ×

One exception is 11.75 mm, where the predicted amount of fibers is about 5.6 times higher, compared to the measurements. Nevertheless, compared to the initial state, nearly 87 % of the 11.75 mm fibers which break due to the measurement are also broken within the simulation.

Figure 4.26 shows the simulated average fiber length for the beginning and end of the sprue, being critical transit areas for fiber breakage.

Figure 4.25: Length weighted amount of fibers with lengths of 0.25 mm to 4.75 mm. Experimental results of screw chamber (black) and in plate (green) with corresponding standard deviation. Simulation results in blue ×

By transition from the screw chamber to the sprue (Figure 4.26a), the average fiber length decreases already about 5 %. A state-of-the-art simulation often starts at the beginning of the sprue and the screw chamber and nozzle are neglected. As can be seen, a proportion of fibers is already broken at this point and the assumption of a fiber length distribution similar to the stored material is not valid at this point.

Figure 4.26b further highlights the transit from sprue to mold as an important area of fiber breakage, where the average fiber length is reduced by 10 %, due to the sharp edge, creating high shearing rates and relative velocities.

Figure 4.26: Simulated average fiber length at the end of filling, weighted by number. Image detail of transition from screw chamber to sprue (a) and sprue to mold plate (b)

The simulation results are in good agreement with the experimental data, validating the novel approach of simulating fiber breakage based on eigenvectors and hydrodynamic forces on macroscopic scale. The fiber breakage is predicted at meaningful areas within the part. Of course, shortening of the fibers affects the anisotropic viscosity modeling, depending on the fiber aspect ratio. Hence, the simulated cavity pressure is also influenced by the fiber breakage, which will be discussed in the next Section.

### **4.2.4 Cavity pressure**

Figure 4.27 shows the experimental and simulation results of cavity pressure for two positions within the mold, as illustrated in Figure 4.23. The experimental pressure data is again the average of five measurements with corresponding standard deviation. The data fits to the fiber length measurements presented in Section 4.2.3.2 so the process parameters are given by Table 4.6. Two simulations are compared. Both are based on the anisotropic viscosity tensor presented in Section 3.3. One simulation is calculated with fiber breakage, where the fiber length distribution is given in Section 4.2.3.2. 

Figure 4.27: Cavity pressure during mold filling at sensor <sup>1</sup> (solid lines) and <sup>2</sup> (dotted lines) as shown in Figure 4.23. Comparison of experiments (black) and anisotropic modeling without fiber break (red) and with fiber break (green). Experimental lines are the average of five measurements with corresponding standard deviation

The other simulation assumes a constant fiber length of 0.434 mm, based on the measurements in the screw chamber. Since the simulation model is identical to the fiber length simulation, all parameters and boundary conditions are given in Section 4.2.3.1. Both simulations predict the cavity pressure well at position <sup>1</sup> and slightly too high at position <sup>2</sup> being still within the standard deviation of the experiments. The results match the previous work [26], where the anisotropic flow coupling with constant fiber length is validated to be suitable for in-mold pressure simulation during mold filling. The focus point of this comparison is the deviation at <sup>1</sup> of the two simulations after 4.5 s, where the simulation without fiber breakage predicts a higher pressure. The maximum deviation is about 0.43 MPa at 5.45 s, being about 10 %. The deviation is within an area, where simulated pressure is too high, compared to the experimental results, so in summary, the results with fiber breakage are slightly better than without.

As expected, the predicted pressure, or pressure drop between <sup>1</sup> and <sup>2</sup> is higher without fiber breakage, due to the longer fibers, causing a higher viscosity. Although the difference is quite marginal with about 10 %, it corresponds to only a slight change of fiber length with about 14 % (see Section 4.2.3.2) and highlights the benefit of a fiber length dependent viscosity modeling with parallel fiber breakage calculation. However, the predicted pressure at <sup>2</sup> is too high for both simulations. One reason may be the no-slip boundary condition, which should be further investigated. Sliding along the wall, would lower the pressure by lowering the kinetic energy, needed for material transport. Of course, this would affect both regarded positions <sup>1</sup> and <sup>2</sup> , but <sup>2</sup> will be affected stronger, due to the longer flow path.

### **4.3 Outlook on prediction of warpage simulation**

This Section will give a short outlook on the potential of predicting warpage and residual stresses, presented in Section 3.5. Due to a lack of information about material data and experimental data, the approach cannot be validated within this work. But at least the general feasibility of such a simulation will be shown.

### **4.3.1 Model and procedure**

#### **4.3.1.1 Initial state of warpage analysis**

The warpage analysis begins after mold filling, including the process steps holding, curing, ejection and out-of-mold cooling. Four different cases are considered. On the one side, the time period of the curing step is performed for 10 s and 60 s, on the other side both configurations are simulated with an ejection step and without. The warpage analysis is performed with the FEM software Abaqus 2018 (Dassault Systémes, Johnston, RI, USA).

The simulation model is a rectangular plate, similar to the one shown in Figure 4.23, except the thickness of the plate is only 3 mm not 4 mm so the model is similar to [26], where the anisotropic flow modeling without fiber breakage is validated. Here, fiber breakage is ignored, since the algorithm for homogenization and material creation can only handle constant fiber aspect ratios. Therefore, the short fiber material and mold used in [26] are regarded. As input for the warpage analysis, fiber orientation, temperature and degree of curing distribution at the end of the mold filling simulation are mapped, using MPCCI MapLib [100,101]. Therefore, all relevant parameters for material modeling are identical to the mold filling simulation in the initial state.

#### **4.3.1.2 Material modeling and boundary conditions**

During the holding and curing step, displacement is blocked on the whole part's surface, since it is still in mold in the real process. The surface temperature is set constant to 170 °C. The holding stage applies for 40 s, the curing step for 10 s or 60 s. Afterwards, an ejection step of 3 s is performed, followed by an out-of-mold cooling step of 6000 s. During the latter two steps, only the top end of the sprue is blocked for displacement to fix the model for numerical stability. Over the course of the ejection step, the corner nodes of the plate are displaced 20 mm in <sup>3</sup> -direction during the complete ejection step. In the cooling step, the corner nodes have no specific displacement boundary condition, similar to the rest of the model. During ejection and cooling, a convection boundary condition for temperature is applied on the

whole part. The corresponding film coefficient and sink temperature are 15 W/(m²K) and 20 °C. The effective material properties are determined with the schemes listed in Table 3.2 and explained in Section 2.3.2 and Section 3.5.2. For the mechanical model, the CHILE model (Section 2.3.1.2) is used in an orthotropic formulation. The model parameters are given in Table 4.7 and Table 4.8. The part is clustered in 30 different material sections, depending on orientation state, so similar orientation states are summarized. The sections must not be contiguous, but have similar fiber orientations, which are normally nearby. According to the 30 material sections, 30 orthotropic materials are created, describing the material behavior also for the corresponding fiber orientation state. This procedure is performed to lower the numerical effort and increase numerical stability.


Table 4.7: Material parameters of fiber and matrix for material modeling of warpage analysis


Due to a lack of material data, the parameters are not correct for the rubbery state and glass transition of the matrix for the corresponding material, but orientated at the work of Bernath [14]. Nevertheless, the aim is a proof of concept, not generating quantifiable results. The values for chemical shrinkage and specific heat capacity are known for the compound and therefore chosen identical for matrix and fiber for simplification, to fit to the value of the compound. For the homogenization, a fiber volume content of 0.35 and an aspect ratio of 15 is applied.


Table 4.8: Glass transition parameters of the matrix for material modeling of warpage analysis

### **4.3.2 Results of the warpage analysis**

The results of the warpage analysis are shown in Figure 4.28 for 60 s curing time and Figure 4.29 for 10 s curing time. In both figures, the deformation on the part is scaled up five times, so the general shape of the final part is detectable. After 60 s of curing, the degree of cure is above 0.75 at every point of the part and above 0.9 within the plate. Therefore, it is <sup>g</sup> > 190 °C within the plate and the behavior is linear elastic at the point of ejection.

Figure 4.28: Magnitude of displacement after 40 s holding, 60 s curing, ejection and 6000 s cooling. Simulation without ejection deformation (a) and with ejection deformation (b)

Figure 4.29: Magnitude of displacement after 40 s holding, 10 s curing, ejection and 6000 s cooling. Simulation without ejection deformation (a) and with ejection deformation (b)

Due to the linear elastic behavior for the complete part, the simulation results with ejection deformation (Figure 4.28b) are identical to the results without ejection deformation (Figure 4.28a), because the ejection deformation simply goes back to its initial state during cooling. Therefore, the ejection does not influence the warpage in this case and the deformations must have a different source. The displacement shown within Figure 4.28 is caused by frozen-in residual stresses during curing and thermal shrinkage in combination with anisotropic behavior. The maximum displacement is about 7.5 mm at the edge with maximum <sup>1</sup> -value and a convex shell is built.

In case of only 10 s curing time, the results with ejection deformation (Figure 4.29b) and without (Figure 4.29a) are different, since the material does not behave linear elastic during ejection. When ejecting the part after only 10 s curing it is <sup>g</sup> < 170 ° in most regions of the part, therefore the behavior is entropy elastic (rubbery state). Hence, a part of the ejection deformation is captured due to the change of material behavior during ejection and the final deformation is higher in Figure 4.29b, compared to Figure 4.29a.

Comparison of Figure 4.28a and Figure 4.29a, shows lower deformations for lower curing times, resulting in a maximum displacement of only about 5.5 mm. This is caused by less frozen-in stresses at the point of ejection, since the material is in rubbery state and no residual stresses can build up. The general distribution of deformation and final part shape is quite similar in both cases, since is mainly driven by thermal shrinkage and anisotropy, where the latter is identical in all simulations. In reality it would be unlikely, that a part with such small curing time, and hence degree of cure will deform less. Effects such as movement of the part and gravity will deform the still visco-elastic and not completely form-stable part, but are neglected within the simulation.

In summary, the orthotropic CHILE model in combination with mapping of process fields is able to capture warpage of fiber reinforced injection molding parts, depending on process parameters, fiber orientation distribution and curing reaction. The displacements are in a realistic scale. Nevertheless, the material parameters in the area of the rubbery state and glass transition are only estimated and not based on experimental data. Further investigations on material parameters and warpage experiments with injection molding parts must be performed to validate the method.

# **5 Conclusion and outlook**

## **5.1 Conclusion**

Novel simulation approaches for fiber reinforced injection molding have been presented. The state-of-the-art is extended by a multiphase approach, considering air and FRP within the simulation. A phase-dependent boundary condition enables that the air may leave the mold, while the FRP is blocked. The boundary condition is formulated to minimize the material loss during simulation.

To consider fiber length and orientation within the flow simulation, the viscosity is modeled non-Newtonian with a fourth order viscosity tensor, taking fiber orientation, length and volume fraction into account. Numerical studies show the influence of fiber length and orientation on flow front evolution and on in-mold pressure, highlighting the benefits of the tensorial approach.

New approaches approximate hydrodynamic forces, fiber-fiber contact points and fiber-fiber forces within the homogenized material with information provided by the second order orientation tensor. The contact point and force modeling is verified with numerical experiments, including randomly generated individual fibers. The novel approaches show good results for different orientations states and fiber lengths.

The hydrodynamic forces are used to model fiber breakage with orientation and flow dependency during mold filling. The breakage modeling is validated with fiber length distributions and pressure data of injection molding experiments. The numerical results are in good agreement with the experimental data. In combination with the viscosity tensor, the fiber breakage directly influences the flow modeling, which is also detectable for the simulated in-mold pressure.

After mold filling, relevant process data is mapped to perform a nonisothermal structural analysis with respect to fiber orientation, length and curing distribution. Therefore, the part is clustered into a defined number of different materials, which are all orthotropic and the values are determined by homogenization and orientation averaging. The simulation results are meaningful and show a prediction of warpage with respect to non-isotropic material behavior, curing and glass transition in a realistic scale.

In summary, the novel approaches improve the quality and level of detail for injection molding simulations with fiber reinforced polymers. A more detailed prediction of air traps, flow front evolution, cavity pressure and material internal forces, depending on fiber orientation and length distribution and fiber volume fraction is enabled. In combination with the subsequent warpage prediction, these approaches help to improve the final part's quality. Furthermore, the higher degree of knowledge during the form filling supports the process and product design in an early stage of development, which decreases the loops of part design optimization and tool revision and therefore has a positive effect on costs, economy and environment.

### **5.2 Outlook**

Although the presented novel approaches improve the quality of an injection molding simulation, there is still unused potential for more improvement which needs to be regarded and validated. This concerns all process steps mentioned within this work.

The calculated fiber-fiber interaction points and forces have no effect on the simulation at this point of time, which is of course not the case in the real process. The mentioned aspects may be used to improve fiber breakage calculations by more detailed force modeling or dynamic description of model parameters, used within the breakage model. Furthermore, they offer potential for a more detailed fiber orientation modeling and even flow modeling itself. This applies for viscosity modeling and for the momentum balance equation itself, since the interaction forces represent material internal forces. Nevertheless, the qualitive and quantitative influence of fiber-fiber contact is not well studied at this point of time and the resulting consequences on fiber breakage and material behavior are not completely known. Therefore, further investigations from experimental and numerical side are needed to improve the state of knowledge on these phenomena. Of course, these effects are on microscopic scale and cannot be separated from each other, creating a high experimental effort.

Additionally, some boundary conditions can be improved, by a more detailed approach. This concern the wall-slip effect, which should be further regarded experimental and simulative. Influence of tool, temperature, pressure, viscosity and of course velocity and material on this friction contact should be quantified and taken into account in the simulations. Another aspect is the temperature boundary condition, which is often chosen homogeneous and constant at state of the art, or the tool with cooling/heating channels is also part of the simulation model, creating higher numerical effort. A more detailed description of heat flow as boundary condition improves the temperature modeling without significant rise of numerical effort, since no tool must be regarded. Here also, the heat transfer coefficient should be modeled with respect temperature, pressure, curing state, etc. Both extensions of the boundary conditions may improve the simulation results, especially in case of inmold pressure and temperature modeling.

Furthermore, more investigation on the early state of warpage analysis shown within this work should be done. In a first step the shown approaches must be validated (or modified) based on experimental data, which also includes experiments to determine relevant material properties. Furthermore, the method can be improved by more detailed material modeling. The fiber aspect ratio should be included into the mapping and homogenization process, so the warpage analysis can consider the length distribution.

# **6 Literature**


domain Tait equation of state for injection molding. Materials & Design 2019;183(1):108149.


### **Karlsruher Schriftenreihe Fahrzeugsystemtechnik FAST Institut für Fahrzeugsystemtechnik (ISSN 1869-6058)**

Eine vollständige Übersicht der Bände finden Sie im Verlagsshop



Die Bände sind unter www.ksp.kit.edu als PDF frei verfügbar oder als Druckausgabe bestellbar.


### Karlsruher Schriftenreihe Fahrzeugsystemtechnik

This work presents novel simulation techniques for injection molding of fiber reinforced polymers. Injection molding is one of the most applied processes for discontinuous reinforced polymers parts. The process conditions such as filling rate, temperature, etc. have a significant influence on the final part properties. For an adequate prediction of these properties, a process simulation has to depict different aspects, including all process steps, being mold filling, holding, in-mold solidification and out-of-mold cooling. To enable a fiber-dependent mold filling simulation, an anisotropic forth order viscosity tensor, considering fiber orientation, length and volume fraction is used within this work. Novel approaches to approximate the hydrodynamic and fiber-contact forces within the homogenized material are presented. The hydrodynamic forces are used in a novel fiber breakage modeling approach to predict the transient fiber length distribution during mold filling. Due to the anisotropic fiber flow coupling, the fiber breakage directly influences the flow modeling and hence the modeled cavity pressure. The fiber breakage and cavity pressure are validated with experimental data, showing good agreement. Furthermore, novel approaches for anisotropic and cure-dependent material modeling during holding and cooling stage are presented, to predict residual stresses and warpage.

F. Wittemann

Methods for fiber-dependent injection molding simulation

**Band 101**