Nils Meyer

# **Mesoscale simulation of the mold filling process of Sheet Molding Compound**

Nils Meyer

### **Mesoscale simulation of the mold filling process of Sheet Molding Compound**

### **Karlsruher Schriftenreihe Fahrzeugsystemtechnik Band 98**

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.

## **Mesoscale simulation of the mold filling process of Sheet Molding Compound**

by Nils Meyer

Karlsruher Institut für Technologie Institut für Fahrzeugsystemtechnik

Mesoscale simulation of the mold filling process of Sheet Molding Compound

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

von Nils Meyer, M.Sc.

Tag der mündlichen Prüfung: 14. September 2021 Hauptreferent: Dr.-Ing. Luise Kärger Korreferent: Prof. Andrew N. Hrymak Korreferent: Prof. Dr.-Ing. Frank Henning

**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-1173-1 DOI 10.5445/KSP/1000143703

## **Kurzfassung**

Sheet Molding Compounds (SMC) sind diskontinuierlich faserverstärkte Verbundwerkstoffe, die aufgrund ihrer Fähigkeit, Verbundbauteile mit langen Fasern zu geringen Kosten zu realisieren, weit verbreitet sind. Sie ermöglichen Funktionsintegration, wie etwa den Einsatz von Rippen oder metallischen Einsätzen, und können mit kontinuierlichen Kohlenstofffasern gemeinsam verarbeitet werden, um die Formbarkeit von SMC mit den überlegenen mechanischen Eigenschaften von kontinuierlichen Fasern zu kombinieren. Das Streben nach hochintegrierten und komplexeren SMC-Bauteilen erfordert jedoch ein tiefes Verständnis der Verarbeitungsmechanismen und deren Einfluss auf die Leistungsfähigkeit eines Bauteils. Prozesssimulationen adressieren diesen Punkt, indem sie mögliche Fertigungsfehler und Prozessparameter vorhersagen. Diese Ergebnisse können nicht nur zur Prozessauslegung und zur Reduzierung von Trial-and-Error-Phasen genutzt werden, sondern auch für die anschließende Struktursimulation durch eine virtuelle Prozesskette.

In dieser Arbeit wird die Prozesssimulation von SMC zunächst mit einem makroskopischen Referenzmodell auf Basis von Faserorientierungstensoren adressiert. Dies entspricht dem Stand der Forschung, aber die zugrundeliegenden Annahmen von geraden Fasern, die viel kürzer als jedes geometrische Merkmal sind, werden in anspruchsvollen SMC-Anwendungen oft verletzt. Dies führt zu der Hypothese, dass eine direkte Simulation einzelner Faserbündel erforderlich ist, um den SMC-Formfüllprozess komplexer Geometrien genau zu beschreiben. Basierend auf dieser Hypothese wird eine neuartige direkte Bündelsimulationsmethode (DBS) vorgeschlagen, die eine direkte Simulation auf Komponentenebene ermöglicht und dabei die

Beobachtung nutzt, dass Faserbündel während des SMC Fließpressens oft in einer gebündelten Konfiguration verbleiben. Das entwickelte DBS Modell kann mit Patches kombiniert werden, um den Co-Molding-Prozess von SMC mit kontinuierlichen Faserverstärkungen zu simulieren. Daher wird ein Modell zur Beschreibung des Materialverhaltens von unidirektionalen Kohlenstofffaser-Patches einschließlich eines einfachen Schädigungsmodells zur Vorhersage von Defekten entwickelt.

Die Parameter des makroskopischen Referenzmodells, des DBS-Modells und des Patch-Modells werden experimentell bestimmt. Dazu gehören die thermischen Eigenschaften des SMCs, die temperaturabhängige und ratenabhängige Viskosität der SMC-Paste, die Reibung an der Werkzeugwand sowie die Kompressibilität des SMCs. Ebenso werden die temperaturabhängigen und ratenabhängigen mechanischen Eigenschaften der Patches bestimmt, die jedoch große Streuungen zwischen den Proben und Chargen aufweisen.

Schließlich werden die Modelle auf mehrere Validierungsfälle angewandt, um die Anwendbarkeit auf Komponentenebene zu bewerten. Die Beispiele zeigen eine gegenüber dem makroskopische Referenzmodell verbesserte Vorhersage der Faserarchitektur, insbesondere der Faserorientierung in der Nähe von Werkzeugwänden sowie der Vorhersage von Bindenähten und Fließmarken. Zusätzlich bietet das DBS Modell die Option, Krümmungen der Bündel vorherzusagen und den Faservolumenanteil zu berechnen, welche durch Mikro-Computertomographie, thermisch gravimetrische Analysen und Durchleuchtungsbilder validiert werden.

## **Abstract**

Sheet Molding Compounds (SMC) are discontinuous fiber reinforced composites that are widely applied due to their ability to realize composite parts with long fibers at low cost. They enable function integration, e.g. ribs or metallic inserts and can be co-molded with continuous carbon fibers to combine the formability of SMC with the superior mechanical properties of continuous fibers. However, the pursuit for highly integrated and more complex SMC components requires a profound understanding of the processing mechanisms and their influence on the performance of a component. Process simulations address this point by predicting possible manufacturing defects and process parameters. These results can not only be used to configure the process and reduce trial and error phases, but also for subsequent structural simulation by virtue of a virtual process chain.

In this work, the process simulation of SMC is addressed initially with a macroscopic reference model based on fiber orientation tensors. This is in line with the state of research, but the underlying assumptions of straight fibers that are much shorter than any geometric feature are often violated in advanced SMC applications. This leads to the hypothesis that a direct simulation of individual fiber bundles is required to accurately describe the SMC mold filling process of complex geometries. Based on this hypothesis, a novel Direct Bundle Simulation (DBS) method is proposed to enable a direct simulation at component scale utilizing the observation that fiber bundles often remain in a bundled configuration during SMC compression molding. The developed DBS model can be combined with patches to simulate the co-molding process of SMC with continuous fiber reinforcements. Hence, a model is developed to describe the material behavior of unidirectional carbon fiber patches including a damage model to predict defects.

The parameters of the macroscopic reference model, DBS model and patch model are determined experimentally. This includes thermal properties of the SMC, temperature dependent and rate dependent viscosity of the paste, friction at the mold surface as well as the compressibility of the SMC. Likewise, the temperature dependent and rate dependent mechanical properties of the patches are determined, but they show large scatter between samples and batches.

Finally, the models are applied to several validation cases for evaluating the applicability at component scale. The examples show improved prediction of fiber architecture compared to the macroscopic reference model, especially the fiber orientation near mold walls and the prediction of knit lines or flow marks. In addition, the DBS model provides the option to predict bundle curvatures and calculate fiber volume fractions, which are validated by microcomputed tomography, thermal gravimetric analysis and fluoroscopy.

## **Acknowledgments**

The research documented in this thesis was carried out during my time at the Institute of Vehicle System Technology - Lightweight Technology (FAST-LT) at Karlsruhe Institute of Technology (KIT). I would like to thank Dr.-Ing. Luise Kärger for the opportunity to work in her Young Investigator Group at FAST and for the supervision of this thesis. The conditions for research in her group are excellent thanks to her dedication and the atmosphere she creates in the group. I would like to thank her for the helpful discussions and detailed, thoughtful comments on my research. Additionally, I would like to thank Professor Andrew N. Hrymak for the opportunity to conduct parts of the documented research under his supervision at Western University in Canada. His extensive experience in composite processing was immensely helpful to shape the focus of this thesis. Further, I would like to thank Prof. Dr.-Ing. Frank Henning for forming a trustful multi-disciplinary cooperation between FAST-LT and its partner institutions, without which this work would be impossible.

I would also like to thank my colleagues at FAST-LT, Fraunhofer Institute of Chemical Technology and in the International Research Training Group GRK 2078 for the friendly work atmosphere and their immense support. The genuine helpfulness is unparalleled and I would like to highlight the fruitful discussions in our office with Florian Wittemann, Martin Hohberg, Alexander Jackstadt and Louis Schreyer, the support in experimental setups by Ludwig Schöttl, Sergej Ilinzeer, Lucas Bretz and Miriam Bartkowiak as well as the collaborations with Johannes Görthofer, Sven Revfi and Julian Bauer. Additionally, I would like to thank all students, who supported my research and explored additional ideas.

The research documented in this thesis has been funded by the German Research Foundation (DFG) within the International Research Training Group "Integrated engineering of continuous-discontinuous long fiber-reinforced polymer structures" (GRK 2078). The support by the German Research Foundation (DFG) is gratefully acknowledged.

Finally, I would like to thank my family and friends for the support on this journey. This is especially true for Lilian, who gave me emotional support in stressful times, shared joy in successful times and provided very practical advice in linguistic aspects of my research documentation.

Karlsruhe, September 2021 *Nils Meyer*

## **Contents**





## **Acronyms and symbols**

### **Acronyms**



### **Scalars**


xii



xiv



### **Vectors**



### **Tensors (2nd order)**



### **Tensors (4th order)**


### **Sets**


### **Operators and math symbols**



## **1 Introduction and motivation**

Lightweight design is a holistic development approach, which ideally triggers a positive self-reinforcing effect: If mass is saved, inertial loads on a structure and loads from powering a machine are reduced. This scales down requirements on support structures, which then can be designed lighter themselves to reduce mass even further. Ultimately, this reduction in mass reduces primary energy consumption for all mass related loads in a machine and decreases impact energy in a catastrophic failure.

One option to reduce mass is the application of composite materials. Such materials combine properties of multiple materials to obtain a tailored property profile suited to the specific application. Fiber reinforced plastics (FRP) are an important group of composites comprising strong load-carrying fibers (typically made from glass or carbon) and a matrix material (thermoset or thermoplastic polymer), that bonds fibers together, transfers loads to fibers and protects fibers from environmental influences [1]. If such reinforcement fibers span the entire length of a part, the reinforcement type is classified as *continuous FRP* (CoFRP). Continuous fiber reinforcements offer high potential for mechanical performance, especially regarding strength, due to large fiber volume contents and high alignment of fibers. Consequently, they are a popular choice in the aerospace industry, performance cars and sports products requiring high stiffness. However, the necessary raw materials of CoFRPs are rather expensive, the degree of process automation is limited, and careful consideration is needed to make economically viable material choices [2]. Additionally, the geometric shapes are restricted by the underlying architecture of fabrics - sometimes it is simply not possible to place fibers such that they span an entire part. Chopped fibers, on the other hand, offer

much more design freedom as they can be molded together with the matrix material in a flow process such as injection molding or compression molding. Such a material is classified as *discontinuous FRP* (DiCoFRP). DiCoFRPs can be molded to complex geometries at high rates and thus are very cost efficient. However, their mechanical properties are inferior to CoFRPs. Discontinuous fiber reinforced plastics with local continuous fiber reinforcements (CoDi-CoFRP) are a new class of composites aiming to combine the merits of both constituents [3]. They use continuous high performance fibers in critical loading areas and a discontinuous compound to realize function integration with ribs, beads or other complex geometric features.

The mechanical properties of FRPs depend significantly on the manufacturing process, which affects key properties such as fiber orientation and fiber volume fraction. Prediction of these properties by numerical simulations can help to account for these effects in an early product development phase and to save expensive trial and error studies. Process simulations aim to predict manufacturing defects and material properties after processing an FRP. Ideally, these results are fed to subsequent structural simulations in the framework of a continuous CAE chain for more accurate results of structural simulations [4–6].

Therefore, this work addresses the simulation of a Sheet Molding Compound (SMC) CoDiCoFRP manufacturing process. Most current models for SMC compression molding simulations assume a macroscopic behavior without accounting for the underlying fiber bundle architecture of SMC at mesoscale, which can lead to inaccurate predictions of fiber orientations after processing in complex geometries. The objective of this work is an accurate prediction of the manufactured fiber bundle architecture and the possibility to incorporate local unidirectional reinforcement patches in the simulation. The developed methods should result in a better understanding of defects during the early design phase of SMC components and thus pave the way to more structural applications of CoDiCoFRP. These applications have the potential to reduce primary energy consumption in a cost-effective manner.

## **2 Fundamentals and state of research**

## **2.1 Modeling scales**

Fiber reinforced composites can be analyzed and modeled at different scales ranging from the molecular scale at which chemical processes occur (curing, fiber sizing) up to the component scale at which a part is loaded. This work terms the smallest scale *microscale*, which refers to the scale of individual fiber filaments with a diameter of approximately 15 µm. The term *mesoscale* refers to fiber bundles or fiber tows that consist of hundreds of individual fiber filaments. The dimensions of fiber bundles cross sections are in the order of 0.1 mm, while the length is 25 mm to 50 mm. Modeling approaches at mesoscale treat the entire bundle as a single instance with no resolution of the underlying micro-structure. The *macroscale* refers to a homogeneous macroscopic scale, at which the micro- and meso-structure of the composite is completely blurred and described by effective macroscopic properties. The scales in this work are summarized in Figure 2.1.

**Figure 2.1:** Schematic illustration of the modeling scales in fiber reinforced composites. Fibers are at the micro scale, bundles at the mesoscale and the homogeneous material with blurred micro- and meso-structure is a representation at the macroscale.

## **2.2 Fiber suspensions**

This section starts with an introduction to general properties of fiber suspensions and subsequently states deformation measures and balance equations. After introducing stochastic descriptions of the fiber orientation state at macroscale, the section summarizes approaches to model rheology of fiber suspensions at macroscale. Finally, direct numerical models are reviewed and classified.

### **2.2.1 General properties and classification**

A fiber suspension is a material consisting of elongated particles suspended in a viscous matrix material. Such particles are characterized by their aspect ratio

$$r\_{\rm p} = \frac{L}{d} \tag{2.1}$$

with the (equivalent) diameter *d* and length *L*. The aspect ratio describes the dimensionless slenderness of a particle. Another dimensionless shape factor is defined as

$$\xi = \frac{r\_{\rm p}^2 - 1}{r\_{\rm p}^2 + 1} \tag{2.2}$$

and simplifies the notation of fiber orientation evolution equations later in this section.

A suspension of reference volume *V*<sup>r</sup> that contains *N*<sup>p</sup> particles can be described either by the number density

$$m\_{\rm P} = \frac{N\_{\rm P}}{V\_{\rm r}}\tag{2.3}$$

or the fiber volume fraction

$$f = V\_{\rm p} n\_{\rm p} = \frac{\pi d^2 L}{4} n\_{\rm p},\tag{2.4}$$

where *V*<sup>p</sup> described the volume occupied by one cylindrical fiber.

Fiber suspensions are commonly divided in different concentration regimes [7] and different theories apply in each of these regimes:

**Dilute regime** Particles are spaced at a large distance compared to their dimensions. Thus, the particles can be treated individually, regardless of the existence of other particles. The effect on the surrounding fluid flow is typically neglected. The dilute concentration regime is formally bound by

$$m\_{\rm p} \ll \frac{1}{L^3}.\tag{2.5}$$

**Semi-dilute regime** Particles interact with each other through flow field perturbations. They can be distributed isotropic and are well below the nematic transition [7]. The semi-dilute regime is formally bound by

$$\frac{1}{L^3} \ll n\_{\rm p} \ll \frac{1}{dL^2}.\tag{2.6}$$

**Concentrated regime** The average distance between particles is of the order of their diameter. They interact via hydrodynamic effects and contacts [8]. The concentrated regime is formally bound by

$$\frac{1}{dL^2} \ll n\_{\rm p.} \tag{2.7}$$

Typical industrial composites belong to the concentrated regime, as illustrated in Figure 2.2.

**Nematic regime (liquid crystal)** Particles are suspended freely but remain in a constrained alignment, as this is the only possibility to achieve a corresponding high fiber volume fraction.

**Figure 2.2:** Concentration regimes: The black lines separate the dilute regime (white) from the semi-dilute regime (light gray) and the concentrated regime (gray).

The bending stiffness of a flexible cylindrical fiber in a suspension may be described by the dimensionless ratio

$$\mathbf{S}^\* = \frac{E\pi}{4\eta\dot{\eta}^\* r\_\mathrm{p}^4} \tag{2.8}$$

with the Young's modulus of a fiber *E*, scalar shear rate *∞*˙ and matrix shear viscosity *¥* [9]. Fibers bend with an increased aspect ratio, higher shear rates and reduced bending stiffness [10], which is summarized by a small value of *S*§.

### **2.2.2 Rate of deformation measures**

The motion of a fiber suspension is characterized by its differentiable velocity field *<sup>v</sup>* : *<sup>D</sup>* ! <sup>R</sup><sup>3</sup> on a domain *<sup>D</sup>* <sup>Ω</sup> <sup>R</sup>3. The position within that domain is denoted by *x* 2 *D* and differential operators work with respect to *x*, if not mentioned otherwise. Relevant objective measures for the deformation are the symmetric strain rate tensor

$$D = \frac{1}{2} \left( \text{grad}\, v + \text{grad}\, v^{\top} \right) \tag{2.9}$$

and the vorticity tensor

$$\mathbf{W} = \frac{1}{2} \left( \mathbf{grad} \, v - \mathbf{grad} \, v^{\top} \right), \tag{2.10}$$

where grad denotes the gradient operator. The spherical part of the strain rate tensor

$$D^\diamond = \frac{1}{3}\operatorname{tr}\{D\}I = \mathbb{P}\_1 : D \tag{2.11}$$

describes the volumetric rate of change and vanishes in case of incompressible fluids. The tensor <sup>P</sup><sup>1</sup> *<sup>=</sup>* <sup>1</sup> <sup>3</sup> *I* ≠*I* denotes the spherical projection tensor of fourth order. The deviatoric part of the strain rate tensor is denoted as

$$D' = D - D^\circ = \mathbb{P}\_2 : D \tag{2.12}$$

and describes the rate at which the fluids shape changes. The tensor P<sup>2</sup> *=* I*<sup>S</sup>* °P<sup>1</sup> denotes the deviatoric projection tensor of fourth order and I*<sup>S</sup>* is the symmetric part of the fourth order identity tensor. A scalar measure for this deformation rate is the shear rate *<sup>∞</sup>*˙ *<sup>=</sup>* <sup>p</sup> 2*D*<sup>0</sup> : *D*<sup>0</sup> .

### **2.2.3 Balance equations**

The solution fields of the deformation of fiber suspensions are described by several balance equations. These balance equations describe conservation properties observed in nature, e.g. conservation of mass, momentum and energy. A general form of such a balance equation for an arbitrary differentiable scalar field (*•*) : *D* ! R in an Eulerian frame is

$$\frac{\partial(\bullet)}{\partial t} + \text{div}((\bullet)\upsilon) = \text{s} \tag{2.13}$$

with the spatial time derivative *@*(*•*) *@<sup>t</sup>* at a fixed point in space, the divergence operator div(*•*) and source term *s*. The source term *s* summarizes the production of the property, non-convective flux of the property and externally applied sources. The spatial time derivative is related to the material derivative via

$$\dot{\theta}(\dot{\bullet}) = \frac{\partial(\bullet)}{\partial t} + \upsilon \cdot \text{grad}(\bullet) \tag{2.14}$$

and describes the change of a property at a deformed material point.

The first important field of interest is the mass density *Ω* : *D* ! R. As mass can be neither produced (at least not in Newtonian mechanics) nor transported in other ways than by convection, the source term vanishes and leads to

$$\frac{\partial \rho}{\partial t} + \text{div}(\rho \upsilon) = 0 \tag{2.15}$$

for the balance of mass in an Eulerian frame. Application of Equation (2.14) gives the Lagrangian form

$$
\not{\rho} + \rho \operatorname{div}(\upsilon) = 0 \tag{2.16}
$$

with the material derivative of the mass density field *Ω*˙.

Another balanced property is the momentum density field *Ωv* and using the tensor form of (2.13) results in the balance equation for momentum

$$\frac{\partial(\rho v)}{\partial t} + \text{div}\left((\rho v) \otimes v\right) = \text{div}(\sigma) + \rho b \tag{2.17}$$

with the non-convective momentum flux density also known as Cauchy stress tensor and a body force density field *Ωb*. Using the mass balance (2.15) and (2.14), this equation may be recast to the Lagrangian form

$$
\rho \,\dot{\upsilon} = \text{div}(\sigma) + \rho b. \tag{2.18}
$$

Finally, the temperature field *T* : *D* ! R is of interest and described by the balance of internal energy *Ωc*p*T* , where *c*<sup>p</sup> is the specific heat capacity. Without heat sources, this is

$$\frac{\partial(\rho c\_{\rm p}T)}{\partial t} + \text{div}\Big((\rho c\_{\rm p}T)\cdot v\Big) = -\text{grad}(\mathbf{d}) + \boldsymbol{\sigma} : \mathbf{D} \tag{2.19}$$

with the non-convective heat flux *d*, and the dissipated strain energy : *D*. For a constant specific heat capacity, the Lagrangian form of Equation (2.19) is

$$
\rho c\_{\rm p} \dot{T} = -\operatorname{div}(\mathbf{d}) + \boldsymbol{\sigma} : \boldsymbol{D} . \tag{2.20}
$$

### **2.2.4 Stochastic description of fiber orientation states**

The orientation of a single straight fiber, fiber segment or bundle segment is denoted as a unit vector *<sup>p</sup>* <sup>2</sup> *<sup>S</sup>*, where *<sup>S</sup> <sup>=</sup>* {*<sup>p</sup>* <sup>2</sup> <sup>R</sup><sup>3</sup> *<sup>|</sup>* <sup>1</sup> *<sup>=</sup>* ∞ ∞*p* ∞ ∞}. The change of fiber orientation with time *p*˙ is of upmost interest for process simulations of fiber suspensions. The first model to describe the orientation evolution of a single rigid spheroid suspended in an infinitely sized Newtonian fluid domain without buoyancy and inertia was described by Jeffery in 1922 as

$$
\dot{p}\_{\rm l} = W \cdot p + \xi \left( D \cdot p - \left( D : p \otimes p \right) \cdot p \right) \tag{2.21}
$$

with the symmetric strain rate tensor *D* and the skew-symmetric vorticity tensor *W* of the suspending fluid [11]. For pure shear in a 2D domain, this equation reduces to a scalar evolution equation for the fiber orientation angle *µ*

$$
\dot{\theta} = \frac{\dot{\mathcal{Y}}}{2} \left( 1 + \xi \cos 2\theta \right). \tag{2.22}
$$

Strictly speaking, Jeffery's model is only valid for a single ellipsoid in an infinite domain. However, it can be shown, that the equation is also valid for cylindrical fibers, if an equivalent aspect ratio *r*<sup>e</sup> is used [12]. Cox [13] derived such an equivalent aspect ratio based on slender-body theory. Zhang et al. [14] utilized Finite Element Analysis to propose a fitted cubic relation between

*r*<sup>e</sup> and *r*p, which is more accurate for small aspect ratios. Jeffery's model yields reasonable results in the dilute regime, where fibers are dispersed so sparsely that they do not affect each others motion.

Modeling the orientation evolution in semi-dilute fiber suspensions requires a formal treatment of an orientation state consisting of several fiber directions. This can be achieved either by a stochastic description, as described in the following, or modeling of individual fiber instances, as described in the later Section 2.2.6.

The fiber orientation density distribution function ™ : *S* ! *P* gives the probability to find a fiber in direction *p* with *P =* {™ 2 R *|* 0 ∑ ™ ∑ 1}. It fulfills the normalization condition

$$\int\_{\mathcal{S}} \Psi(p) \mathbf{d} \, p = 1 \tag{2.23}$$

and the continuity condition

$$\Psi = -\operatorname{grad}\_{\mathcal{S}}(\Psi(p)\dot{p}).\tag{2.24}$$

Here, d*p* is the surface element on the unit sphere *S* that ensures an invariant integration and grad*<sup>S</sup>* is a gradient on the unit sphere *<sup>S</sup>*. Further, ™(*p*) is a symmetric function, i.e. ™(*p*) *=* ™(°*p*), as the end point of a fiber cannot be distinguished from its starting point.

Folgar and Tucker modified Equation (2.21) to

$$\Psi = -\operatorname{grad}\_{\mathcal{S}}\left(\Psi(\mathcal{p})\dot{\mathcal{p}}\_{\mathbb{I}}\right) + C\_{\mathbb{I}}\dot{\gamma}\operatorname{grad}\_{\mathcal{S}}^2(\Psi(\mathcal{p})) \tag{2.25}$$

by employing the fiber orientation density distribution function and adding a term *C*i*∞*˙ for diffusion [15, 16]. The parameter *C*<sup>i</sup> is termed the interaction coefficient and for *C*<sup>i</sup> *=* 0, Equation (2.25) is identical to Equation (2.21). The diffusion term is inspired by Brownian particle motion and models the tendency of fibers to distribute to an isotropic state due to interactions between individual fibers. The empirical parameter *C*<sup>i</sup> describes the rate at which

fibers move towards a more isotropic state and depends on fiber aspect ratio, fiber volume fraction and other properties of the fiber suspension. [15].

The fiber orientation density distribution function ™(*p*) may have arbitrary complexity. Thus, it is inconvenient to store in numerical simulations and Advani and Tucker [17] suggested to store moments instead. These moments are called fiber orientation tensors and commonly used moments are the second order fiber orientation tensor

$$A = \int\_{S} \Psi(p) p \otimes p \mathrm{d}p \tag{2.26}$$

and fourth order fiber orientation tensor

$$\mathbb{A} = \int\_{S} \Psi(p) p \otimes p \otimes p \otimes p \mathrm{d}p. \tag{2.27}$$

Fiber orientation tensors are fully symmetric and the contractions

$$A \cdot I = 1 \quad \mathbb{A} : I = A \tag{2.28}$$

apply. The set of second order orientation tensors fulfilling such properties is termed *A*. Exemplary fiber orientation states, corresponding second order fiber orientation tensors and plots of the fiber orientation density distribution function are given in Figure 2.3.

Jeffery's model (2.21) and the Folgar-Tucker model (2.25) can be expressed in terms of fiber orientation tensors as

$$\dot{A} = W \cdot A - A \cdot W + \zeta \left( D \cdot A + A \cdot D - 2A : D \right) \tag{2.29}$$

and

$$\dot{A} = W \cdot A - A \cdot W + \xi \left\{ D \cdot A + A \cdot D - 2A : D \right\} + 2C\_l \dot{\gamma} \left\{ I - 3A \right\}, \quad (2.30)$$

respectively. However, the fourth order tensor A is required to solve these equations. One could obtain it from analogous equations for A˙ , but that

**Figure 2.3:** Orientation examples. The distributions (a),(b) and (c) are generated by the procedure described in Section 5.1.2.

would require a sixth order tensor and so forth. Hence, A is typically obtained from a closure approximation. The simplest analytical closure is a quadratic closure

$$A^{Q} = A \otimes A,\tag{2.31}$$

which is correct only for unidirectional fiber orientation. The linear closure

$$\begin{aligned} \mathbb{A}^{\mathsf{L}} = \frac{1}{7} \Big( A \otimes I + I \otimes A + A \square I + I \square A + (A \square I)^{\top\_{\mathsf{R}}} + (I \square A)^{\top\_{\mathsf{R}}} \Big) \\ &- \frac{1}{35} \Big( I \otimes I + I \square I + (I \square I)^{\top\_{\mathsf{R}}} \Big) \end{aligned} \tag{2.32}$$

is exact for a random isotropic fiber orientation state. Hybrid closures try to interpolate these states [17, 18]. More advanced closures approximate the fourth order tensor by fitting eigenvalues such as orthotropic fitted closures [19, 20] and the invariant-based optimal fitting closure (IBOF) [21]. This work proceeds to use the IBOF closure in macroscopic reference solutions, as it is a good compromise between accuracy and computational cost.

The Folgar-Tucker model was extended by several authors and a list (noncomprehensive) is given in Table 2.1. In non-dilute suspensions, Equations

(2.29) and (2.30) overestimate the reorientation rate of fibers in comparison to experimental evidence. Thus, variants with reduced strain rate (RSC) or a retarded principal rate (RPR) have been developed and add up to two additional parameters. Additionally, the isotropic diffusion term has been replaced by an anisotropic term [22], which introduced four additional parameters. A Maier-Saupe potential was suggested to counteract the diffusive process considering nematic conditions in suspensions with high fiber volume fractions [23]. There are also some suggestions to account for confinement of fibers in regions, where the distance between mold walls is smaller than the length of suspended fibers [24, 25]. Recent developments tend to reduce the number of parameters for easier parameter identification [26, 27].

Most of these models are developed primarily for injection molding simulations with a three-dimensional short fiber architecture in mind and the results can be tuned by the choice of parameters. To obtain an unbiased reference solution, the macroscopic models later in this work use Jeffery's basic equation without retarded principal rates and without diffusion parameters.


**Table 2.1:** Overview on stochastic macroscopic fiber orientation models and number of required parameters for the diffusion model and the strain rate reduction. Implementations and examples may be found at

Stochastic fiber orientation descriptions are valid for fibers that are much shorter than any dimension of the flow domain, as they do not account for confinement by mold walls and other non-local effects. A sufficiently high number of fibers is needed at each point, at which an evolution equation is evaluated, to represent a proper sample of the fiber orientation density distribution function. Additionally, these models typically assume a constant fiber volume fraction and constant fiber aspect ratio.

## **2.2.5 Rheology of fiber suspensions**

The total macroscopic stress in an incompressible fiber suspension may be decomposed into

$$
\sigma = -pI + \sigma\_{\rm M}^{\prime} + \sigma\_{\rm FM}^{\prime} + \sigma\_{\rm FF^\*}^{\prime} \tag{2.33}
$$

with a contribution from the volumetric pressure *p*, from the neat matrix fluid <sup>0</sup> <sup>M</sup>, a deviatoric extra stress due to fiber-matrix interactions <sup>0</sup> FM (longrange interactions) and a deviatoric extra stress due to fiber-fiber interactions 0 FF (short-range interactions, such as mechanical contacts and lubrication) [30, 31]. The illustration in Figure 2.4 visualizes the sources of long-range hydrodynamic interactions due to disturbances of the flow field and shortrange interactions at contact points of fiber bundles.

**Figure 2.4:** Illustration of long-range hydrodynamic interactions of bundles with the matrix velocity field and short-range fiber-fiber interactions at contact points between bundles. The hydrodynamic interaction is two-way coupled, i.e. the motion of a bundle depends on the fluid velocity in its environment and the fluid velocity is affected by bundle motion. Adopted from [32].

#### **2.2.5.1 Matrix contribution**

The constitutive equation for the neat matrix fluid

$$\sigma\_{\rm M} = \zeta \text{tr}(\mathbf{D}) \mathbf{I} + 2\eta \mathbf{D}' = \underbrace{\left\{ 3\zeta \mathbb{P}\_1 + 2\eta \mathbb{P}\_2 \right\}}\_{\mathbb{V}\_{\rm kso}} : \mathbf{D} \tag{2.34}$$

computes <sup>M</sup> from a scalar shear viscosity *¥* and a scalar bulk viscosity *≥*, which may be combined to a fourth order isotropic viscosity tensor Viso. The term tr(*D*)*I =* P1*D* vanishes in incompressible fluids and most compressible engineering cases assume *≥ =* 0, as *≥* is only relevant if the reciprocal of the bulk strain rate was at the order of magnitude of molecular movement [33]. Therefore, Equation (2.34) is reduced to

$$
\sigma'\_{\mathcal{M}} = 2\eta \mathcal{D}'\tag{2.35}
$$

for the following sections. Subjecting the isotropic matrix to uniaxial, incompressible elongation (*D*0 *y y = D*<sup>0</sup> *zz =* °*D*<sup>0</sup> *xx* /2) results in the first and second normal stress differences

$$N\_{\rm M,1} = \sigma\_{\rm M,xx}^{\prime} - \sigma\_{\rm M,yy}^{\prime} = 3\eta D\_{xx}^{\prime} \tag{2.36}$$

$$N\_{\rm M,2} = \sigma\_{\rm M,xx}^{\prime} - \sigma\_{\rm M,xz}^{\prime} = 3\eta D\_{\rm x,x}^{\prime} \tag{2.37}$$

of the matrix. The factor 3*¥* may be termed elongational viscosity of the isotropic matrix *¥M*,*xx* and the ratio between this elongational viscosity and the shear viscosity *¥M*,*xx* /*¥* is known as Trouton ratio, which is equal to 3 for isotropic incompressible Newtonian media [33, 34].

### **2.2.5.2 Long-range contribution due to hydrodynamic interaction of fibers with the matrix**

Typical fibers in composite applications have a high Young's modulus and thus are considered inextensible compared to the suspending matrix material. Hence, fibers represent a constraint on the fluid field, disturbing the velocity field at fiber scale and leading to additional shear deformation compared to the macroscopically applied deformation rate. This extra shearing of the viscous matrix translates to extra stress at macroscale that depends on the fiber volume fraction, fiber aspect ratio and fiber orientation. This extra stress results purely from the long-range velocity disturbance even in absence of any mechanical contact of fibers.

The general objective form of anisotropic stress in an incompressible fluid derived by Ericksen [35] and can be expressed as

$$\sigma\_{\rm FM} = -P\_0 I + 2\eta\_0 D + \eta\_1 A + \eta\_2 A : D + 2\eta\_3 \left( A \cdot D + D \cdot A \right) \tag{2.38}$$

using orientation tensors [36]. The first two parameters *P*<sup>0</sup> and *¥*<sup>0</sup> describe the isotropic extra stress by fibers. The three parameters *¥*1,*¥*2,*¥*<sup>3</sup> are the internal fiber stress, elongational extra viscosity in fiber direction and difference between shear viscosity along fibers and transverse to the fibers, respectively [37]. The term *¥*<sup>1</sup> is typically neglected under the assumption that there is no internal strain-rate-independent fiber stress at scales, where Brownian motion is insignificant.

The extra stress from fiber-matrix interactions should not depend on the trace of *D* and should be trace-free itself, such that the term *p* describes the complete hydrostatic pressure [38]. Thus, Equation (2.38) is reformulated to the deviatoric stress

$$\sigma'\_{\rm FM} = \left[2\eta\_0 \mathbb{P}\_2 + \eta\_2 \left(\mathbb{A} - \frac{1}{3}I \otimes A\right) + 2\eta\_3 \left(A \Box I + I \Box A - \frac{2}{3}I \otimes A\right)\right] : \mathcal{D}'. \tag{2.39}$$

This equation is obtained by replacing *D* with the deviatoric strain rate *D*<sup>0</sup> , multiplying Equation (2.38) with P<sup>2</sup> from the left and using the identity *I* : A *= A*. For *¥*<sup>0</sup> *=* 0 and *¥*<sup>3</sup> *=* 0, i.e. slender fibers with negligible thickness in a dilute suspension, Equation (2.39) becomes equivalent to the equation derived by

Batchelor [39]. For dilute suspended elongated spheroids, Batchelor obtained the expression

$$\eta\_{2,\text{Ba}} = \frac{2f r\_{\text{p}}^2}{3\left[\ln(2r\_{\text{p}}) - 1.5\right]} \eta \tag{2.40}$$

using slender body theory [40–42]. Dinh and Amstrong suggest the variant

$$
\eta\_{2, \text{DA}} = \frac{2f r\_{\text{p}}^2}{3 \ln(\sqrt{\pi/f})} \eta \tag{2.41}
$$

for semi-dilute suspensions assuming aligned fibers [42, 43]. Shaqfeh and Fredrickson [44] rigorously derived the expression

$$\eta\_{2, \text{SF}} = \frac{4f r\_{\text{p}}^2}{3 \left[ \ln(1/f) + \ln \ln(1/f) + C \right]} \eta \tag{2.42}$$

for slender fibers using multiple scattering, accounting for velocity disturbances at multiple length scales. The factor *C* ranges from 0.1585 for aligned fibers to -0.6634 for isotropic fiber orientation [30, 44]. These authors also imply *¥*<sup>0</sup> *=* 0, since this contribution is small compared to *¥*2.

Tucker [41] recasts Equation (2.38) together with the isotropic matrix component to

$$
\sigma\_{\rm M} + \sigma\_{\rm FM} = 2\bar{\eta} \left[ D + N\_{\rm p} \mathbb{A} : D + N\_{\rm s} \{ A \Box D + D \Box A \} \right] \tag{2.43}
$$

with *¥*¯ containing the entire isotropic contribution from the matrix and suspended particles. The dimensionless particle number *N*<sup>p</sup> describes the importance of extra stress due to particle stretching and is usually dominant over the shear number *N*s, which describes increased shear resistance due to the thickness of fibers [41]. Strictly enforcing deviatoric stresses and strain rate tensors leads to

$$
\sigma'\_{\rm M} + \sigma'\_{\rm FM} = \mathbb{V}\_{\rm a} : \mathbf{D}' \tag{2.44}
$$

with the fourth order anisotropic viscosity tensor

$$\mathbb{V}\_{\mathbf{a}} = 2\bar{\eta} \left[ \mathbb{P}\_2 + N\_{\mathbf{P}} \left( \mathbb{A} - \frac{1}{3} \mathbf{I} \otimes \mathbf{A} \right) + N\_{\mathbf{s}} \left( \mathbf{A} \Box \mathbf{I} + \mathbf{I} \Box \mathbf{A} - \frac{2}{3} \mathbf{I} \otimes \mathbf{A} \right) \right] \tag{2.45}$$

following Bertóti [38]. This version is equivalent to a combination of Equations (2.39) and (2.35) for a proper choice of the parameters, but has a different interpretation. Such parameters can be obtained e.g. from fitting micro-mechanical results [38,45] or by expressing them in terms of Equations (2.40) to (2.42). Equation (2.43) is also commonly used directly to model anisotropic viscosity with *N*<sup>s</sup> *=* 0 for slender fibers [23, 46].

Other authors do not start from the general form given in Equation (2.38), but with the transversely isotropic case of an unidirectional fiber reinforcement in analogy to linear elasticity [47, 48]. For an incompressible material, this is

$$\begin{bmatrix} D\_{\text{xx}}' \\ D\_{\text{yy}}' \\ D\_{\text{zz}}' \\ D\_{\text{zz}}' \\ D\_{\text{yz}}' \\ D\_{\text{zz}}' \\ D\_{\text{yy}}' \end{bmatrix} = \begin{bmatrix} \frac{1}{\eta\_{\text{xx}}} & -\frac{1}{2\eta\_{\text{xx}}} & -\frac{1}{2\eta\_{\text{xx}}} & 0 & 0 & 0 \\ \left(\frac{1}{4\eta\_{\text{xx}}} + \frac{1}{4\eta\_{\text{yz}}}\right) \left(\frac{1}{4\eta\_{\text{xx}}} - \frac{1}{4\eta\_{\text{yz}}}\right) & 0 & 0 & 0 \\\\ \left(\frac{1}{4\eta\_{\text{xx}}} + \frac{1}{4\eta\_{\text{yz}}}\right) & 0 & 0 & 0 \\\\ & \frac{1}{\eta\_{\text{yx}}} & 0 & 0 \\\\ D\_{\text{yy}}' & & \frac{1}{\eta\_{\text{xy}}} & 0 \\\\ \text{symm.} & & & \frac{1}{\eta\_{\text{xy}}} \end{bmatrix} \begin{bmatrix} \sigma\_{\text{FM,xx}} \\ \sigma\_{\text{FM,xx}} \\ \sigma\_{\text{FM,yy}} \\ \sigma\_{\text{FM,zz}} \\ \sigma\_{\text{FM,zz}} \\ \sigma\_{\text{FM,xx}} \\ \sigma\_{\text{FM,xy}} \end{bmatrix} \tag{2.46}$$

in Voigt notation and for fiber alignment in *x*-direction. Pipes et al. [48, 49] determined the corresponding parameters from micro-mechanical analysis to

$$\eta\_{\chi\chi} = \frac{\eta f}{2} \left( \frac{a\sqrt{f}}{1 - a\sqrt{f}} \right) r\_{\text{p}}^2 \tag{2.47}$$

for the elongational viscosity in fiber direction,

$$
\eta\_{\chi\chi} = \frac{\eta}{2} \left( \frac{2 - a\sqrt{f}}{1 - a\sqrt{f}} \right) \tag{2.48}
$$

for the axial shear viscosity and

$$
\eta\_{\mathcal{Y}^Z} = \eta \left( 1 - a^2 f \right)^{-2} \tag{2.49}
$$

for the transverse shear viscosity. A geometrical factor *a* describes the idealized arrangement of fibers (*a*<sup>2</sup> *<sup>=</sup>* <sup>p</sup> 12/*<sup>º</sup>* for hexagonal packing and *<sup>a</sup>*<sup>2</sup> *<sup>=</sup>* 4/*<sup>º</sup>* for square packing). Similar to the volume averaging suggested by Advani and Tucker [17] for elastic stiffness, the result is then used in a volume-averaged equation that depends on fiber orientation tensors [50].

A summary of the parameterizations for anisotropic viscosity is given in Table 2.2. It is noteworthy that all variants rely on three parameters compared to the five parameters of transversely isotropic materials in elasticity due to the assumption of incompressibility. In principal, these formulations are equivalent to each other. The general form with *¥*,*¥*2,*¥*<sup>3</sup> is the preferred parameterization in subsequent chapters, because it is based directly on the matrix viscosity, which can be measured in a simple rheometer.

**Table 2.2:** Parameters for macroscopic viscosity including matrix and long-range contributions.


Figure 2.5 compares predicted first normal stress differences for a suspension of aligned fibers. Batchelor's model and the Dinh-Amstrong model, which are designed for dilute and semi-dilute suspensions, underestimate the increase of elongational viscosity for higher volume fractions. Pipes' model predicts a singularity as soon as the fiber volume fraction approaches the packing limit (square packing is assumed here). The Shaqfeh-Fredrickson model resembles the properties of Batchelor's model for small volume fractions and those of Pipes' model for larger volume fractions, including a singularity at the packing

**Figure 2.5:** Comparison of the first normal stress difference for a suspension of unidirectional aligned, square packed fibers with aspect ratio *r*<sup>p</sup> *=* 25.

limit. It seems applicable to a relatively large range of fiber volume fractions and is therefore used in the reference model in this work. All models scale with *r* <sup>2</sup> p, which results in a significant increase of the elongational viscosity for high aspect ratios of the fibers.

#### **2.2.5.3 Short-range contribution due to fiber-fiber interactions**

Fibers in concentrated suspensions have multiple contact points at which fibers interact through mechanical friction or lubricated contact with a small matrix layer undergoing high shear. Such direct contacts between individual fibers are summarized as short-range interactions or fiber-fiber interactions. The macroscopic stress contribution can be formally expressed as

$$
\sigma\_{\rm FF} = -\frac{1}{V\_{\rm r}} \sum\_{\langle \alpha, \beta \rangle}^{N\_{\rm c}} h\_{a\beta} \otimes f\_{a\beta} \tag{2.50}
$$

with a gap vector *hÆØ* and an interaction force *fÆØ* for all *N*<sup>c</sup> interaction pairs (*Æ*,*Ø*) of a reference volume *V*<sup>r</sup> [31, 52]. The gap vector *hÆØ* describes the shortest path between two straight fiber segment surfaces and is directed from partner *Æ* to partner *Ø*. The scalar distance between two surfaces is denoted *g* and may become negative, if fibers overlap. The interaction force *fÆØ* comprises an elastic contact force in normal direction *f* n,e *ÆØ* and a tangential friction force *f*t,e *ÆØ* as well as lubrication forces in normal direction *<sup>f</sup>* n,l *ÆØ* and tangential direction *f*t,l *ÆØ*. The focus of this section is the effect of averaged short-range interactions at a macroscopic scale.

The average number of contact points per fiber is given as

$$m\_{\mathbb{C}} = 4f \left( \frac{2}{\pi} r\_{\mathbb{P}} \Phi\_1 + \Phi\_2 + 1 \right) \tag{2.51}$$

with orientation functions

$$\Phi\_1 = \int\_S \int\_S \left| \sin \left( \angle \left( p, p' \right) \right) \right| \Psi(p) \Psi(p') \mathrm{d}p' \mathrm{d}p \tag{2.52}$$

$$\Phi\_2 = \int\_S \int\_S \left| \cos \left( \angle \left( p, p' \right) \right) \right| \Psi(p) \Psi(p') \mathrm{d}p' \mathrm{d}p \tag{2.53}$$

as defined by Toll and Månson [53–55]. Exemplary values of the orientation functions are given in Table 2.3. For a planar network of slender elastic fibers (*r*<sup>p</sup> ¿ 1), the number of contact points per unit volume is consequently

$$\frac{N\_{\rm c}}{V\_{\rm r}} \approx \frac{16}{\pi^2} \frac{f^2}{d^3} \Phi\_{\rm 1} \,\tag{2.54}$$

The average normal force at contact points is then given as

$$\left\| \overline{f^{\rm n,e}\_{\alpha\beta}} \right\| = \frac{32}{5\pi^2} E d^2 \Phi\_1^3 f^3,\tag{2.55}$$

where *E* is the Young's modulus of fibers and the averaging is indicated by a bar [53, 54].

Servais et al. [56, 57] build on these results for their investigation of shortrange fiber-fiber interactions for dispersed fibers and bundled configurations.


**Table 2.3:** Exemplary values of the orientation functions [53].

They propose a power-law relation for the tangential hydrodynamic lubrication

$$\mathbf{f}\_{\alpha\beta}^{\mathrm{t,l}} = k\_{\mathrm{S}} \left\| \Delta \boldsymbol{v}\_{\alpha\beta} \right\|^{n-1} \Delta \boldsymbol{v}\_{\alpha\beta} \tag{2.56}$$

with power-law exponent *n*, hydrodynamic friction coefficient *kS* and the relative tangential sliding velocity ¢*vÆØ*. For the tangential friction, they suggest a Coulomb friction model

$$\mathcal{J}^{\rm t,e}\_{\alpha\beta} = -\mu \overline{\left\| \mathcal{J}^{\rm n,e}\_{\alpha\beta} \right\|} \left\| \Delta v\_{a\beta} \right\|\tag{2.57}$$

with a Coulomb friction coefficient *µ*. The friction acts opposed to the direction of the sliding velocity and the operator Ç*•*É *=* (*•*)/k*•*k is used as a convenient notation to compute the normal direction.

Servais et al. [56] showed that Coulombic friction results in a yield stress of the macroscopic suspension with dispersed fibers and thus the result renders a Herschel-Bulkley fluid. For fiber bundles, they proposed a similar model using a Carreau relationship for the shear rate dependency. However, hydrodynamic lubrication dominates in case of fiber bundles, as the sheared surface area is large compared to the effects of the friction at the contact point [57].

A more general formulation of the short range stress contribution was obtained by Djalili-Moghaddam and Toll [30] for linear Newtonian lubrication forces

$$\mathbf{f}\_{\alpha\beta}^{\mathrm{t,l}} = k\_{\mathrm{D}} \frac{4h\_{\mathrm{D}}^3}{d^2} \eta \Delta v\_{\alpha\beta} \tag{2.58}$$

with a constant interaction range 2*h*D, in which short-range interactions occur, and a parameter *k*<sup>D</sup> that accounts for all geometric details of the lubrication flow. Their final result

$$
\sigma\_{\rm FF} = \frac{4k\_{\rm D}}{3\pi^2} r\_{\rm p}^2 \eta f^2 \mathbb{B} : \mathcal{D} \tag{2.59}
$$

is independent of the interaction range [30]. The key term is the fourth order interaction tensor

$$\mathbb{E} = \int\_{S} \int\_{S} \left\| p \times p' \right\| \, \mathbb{V}(p) \mathbb{V}(p') p \otimes p \otimes p \otimes p \mathrm{d}p \mathrm{d}p' \tag{2.60}$$

that evaluates the interaction between fibers and vanishes for perfectly aligned fibers. Djalili-Moghaddam and Toll [30] compute the interaction tensor from discrete fiber directions (compare Section 2.2.6), but Férec et al. [58] suggest to approximately relate the interaction tensor to the known fiber orientation tensors. They denote the second order interaction tensor as

$$B = \int\_{\mathcal{S}} \int\_{\mathcal{S}} \left\| p \times p' \right\| \,\Psi(p)\Psi(p')p \otimes p \mathrm{d}p \mathrm{d}p' \tag{2.61}$$

and obtain the result

$$B \approx \frac{3\pi}{2} \left( A - \mathbb{A} : A \right). \tag{2.62}$$

by replacing the Onsager potential∞ ∞*p*£*p*<sup>0</sup> ∞ ∞ with the Maier-Saupe potential ∞ ∞*p*£*p*<sup>0</sup> ∞ ∞ <sup>2</sup> and an appropriate prefactor [58]. They initially propose a quadratic closure

$$\mathbb{E}^{\mathbb{Q}} = \frac{1}{\Phi\_1} \mathbf{B} \otimes \mathbf{B} \tag{2.63}$$

and a linear closure [58]. The closures were later improved [59], but the effect of short-range interactions has been investigated in more detail with direct numerical models, as precise knowledge of the fiber architecture is necessary to properly account for this effect.

This work will use equations (2.59), (2.62) and (2.63) later in Chapter 4 for the macroscopic reference model. This set of equations has just a single adjustable parameter *k*<sup>D</sup> and provides a closed formulation of the macroscopic stress contribution from fiber-fiber lubrication, which is rarely addressed in literature.

### **2.2.6 Direct numerical models**

Direct numerical models solve the motion of individual fibers or representative sets of fibers in a suspension instead of a stochastic moment. An equation of motion is solved for fiber segments, bead chains or other discretization forms of a fiber instead of an evolution equation for a fiber orientation tensor.

#### **2.2.6.1 Kinematic models**

The first class of direct simulation models transports fiber nodes according to the bulk velocity of the suspension. These models do not include the momentum balance or any forces and are therefore termed *kinematic models* here. One such model is commercialized as 3D TIMON Direct Fiber Simulation (DFS) feature by Toray Engineering Co., Ltd., Japan. A fiber is represented by a chain of rods connected with nodes as shown in Figure 2.6. These nodes are transported with a pre-computed velocity field *v*1(*x*,*t*) of a macroscopic mold filling simulation in a forward Euler scheme as

$$
\hat{x}\_k^{n+1} = x\_k^n + v\_\infty \left( x\_k^n, t^n \right) \Delta t,\tag{2.64}
$$

where *k* 2 N denotes a node index and *n* 2 N denotes a time step index [60]. A regularization is necessary to constrain fiber elongation and is achieved with a correction scheme (*x*ˆ *<sup>n</sup>+*<sup>1</sup> *<sup>k</sup>* ! *<sup>x</sup><sup>n</sup>+*<sup>1</sup> *<sup>k</sup>* ) based purely on the geometric configuration of the fiber [60]. This kinematic model has been used to analyze fiber bending in Long Fiber Thermoplastic (LFT) compression molding [61] and

shows good agreement of fiber orientation in compression molding of carbon prepreg platelets compared to Micro Computer Tomography (µCT) [62].

**Figure 2.6:** The kinematic direct fiber simulation model computes a new temporary node position *x*ˆ*n+*<sup>1</sup> *<sup>k</sup>* from a prescribed velocity field and the current node position *<sup>x</sup><sup>n</sup> <sup>k</sup>* by forward Euler integration. Afterwards, a regularization is applied to ensure constant fiber length.

This model does not account for short-range interactions between fibers nor for long-range interactions by disturbance of the velocity field. This makes the model very efficient, but has the disadvantage that multiple fibers can occupy the same space and therefore leads to unrealistic fiber volume fractions. The effect of fibers on the macroscopic viscosity has to be factored into the macroscopic process simulation, as fibers are advected based on this pre-computed velocity field. Further, the regularization can lead to 'stuck' fibers in regions with high shear rates [63].

#### **2.2.6.2 Stokesian dynamics**

Most direct simulation approaches compute the motion of suspended particles at low Reynolds numbers (*Re* ø 1, Stokes flow) using an equation of the type

$$M\dot{v}\_{\rm p} = f^{\rm H} + f^{\rm I} + f^{\rm B} \tag{2.65}$$

with particle mass *M*, particle velocity *v*p, hydrodynamic interaction force *f*H, non-hydrodynamic interaction forces *f*<sup>I</sup> and stochastic forces causing Brownian motion *f*<sup>B</sup> [64]. The Brownian motion is neglected for fiber suspensions under investigation here and the hydrodynamic interaction is given as

$$f^{\rm H} = -R\_{\rm V} \left( v\_{\rm p} - v\_{\infty} \right) + R\_{\rm s} D' \tag{2.66}$$

with resistance tensors *R*v, *R*<sup>s</sup> and the surrounding fluid velocity *v*<sup>1</sup> [64]. Equivalent equations apply for the evolution of the rotational velocity of a suspended particle. In case of a suspended sphere with radius *R*, the resistance tensors take the form

$$R\_{\rm V} = 6\pi\eta RI, \quad R\_{\rm s} = 0. \tag{2.67}$$

Stokesian dynamics models do not account for two separate phases with an interface, but use the resistance tensors to capture all effects at a length scale smaller than the suspended particle. In most cases, the motion of fibers is driven by discrete evaluation of a given flow field *v*1, which is undisturbed by the presence of fibers (i.e. no long-range interactions).

Yamamoto and Matsuoka [65] developed a beadchain model, in which a flexible fiber is made up of several connected spheres (see Figure 2.7a). The fiber can bend, stretch and twist through linear elastic forces between spheres. Each sphere experiences a hydrodynamic drag according to Equations (2.66) and (2.67) as well as an equivalent torsion moment from a given fluid field [65]. They could demonstrate agreement of rigid fiber motion with Jeffery's Equation (2.21) and correct deformation of flexible fibers compared to Forgacs and Mason's results [10].

The step from single fibers to fiber suspensions requires a model for shortrange fiber-fiber interactions. A lubrication model for the normal direction between two rigid rods (see Figure 2.7b) was introduced by Yamane et al. [66]

in their simulation of a semi-dilute rigid rod suspension. They proposed the analytical solution based on a lubrication approximation

$$f^{\rm n,l}\_{\alpha\beta} = -12\pi\eta \underbrace{\frac{d^2}{\left\|\mathbf{p}\_a \times \mathbf{p}\_\beta\right\|}}\_{\text{contact area}} \underbrace{\left\|h\_{\alpha\beta}\right\|}\_{\text{normal area}} \underbrace{\left\|h\_{\alpha\beta}\right\|}\_{\text{normal}} \tag{2.68}$$

where ∞ ∞ <sup>∞</sup>*h*˙ *ÆØ* ∞ ∞ <sup>∞</sup> describes the rate at which the cylinders approach each other at the contact area. They report a rather weak effect on fiber orientation (small *C*i) and viscosity in the semi-dilute regime. Similarly, Yamamoto and Matsuoka extended their beadchain model to account for multiple fibers in a periodic cell with lubrication forces [67]. Lubrication forces only do not prevent fibers from penetrating each other at high aspect ratios. Therefore, mechanical contact forces have been added by Sundararajakumar and Koch [68]. The interaction of rigid rods and boundaries was investigated by Thomasset et al. [69].

For flexible fibers, the bead chain model was improved with the use of spheres [70,71] (see Figure 2.7e), spheroids [72] (see Figure 2.7d), rods [73, 74] 2.7f) or combined beads [75, 76] (see Figure 2.7c) connected by sockets and joints to reduce the computational effort for fibers of long aspect ratios. Such models were applied to model flocculation [77], to obtain effective suspension properties [78–80], analyze the effect of elastic fiber bending on the suspension elasticity [81], to model fiber jamming [82], to investigate fiber fracture [75] and fiber buckling [83].

Lindström and Uesaka [85] compute the velocity field on a grid instead of using a prescribed velocity field. The velocity field is solved including a body force field that opposes the hydrodynamic interaction forces (momentum conservation) and therefore includes long-range hydrodynamic interaction. They applied their two-way coupled approach to determine effective properties of semi-dilute fiber suspensions [86, 87]. However, they utilize the drag of prolate spheroids, which leads to a total drag force on a fiber that

(spheroids) [72]

**Figure 2.7:** Schematic overview on different approximations of fibers for Stokesian dynamics.

[70, 77]

depends on discretization [88]. These are the only models in this section, that compute a disturbance of the fluid field and therefore account for long-range interactions.

More recent developments are targeting applications beyond representative volume elements further towards the component scale. For example, Kuhn et al. [89, 90] investigated rib filling during LFT compression molding. Hayashi et al. [91] suggested the use of constrained beams at this scale. Sasayama et al. [92] simplified the original beadchain model by dropping all rotational motion components and apply the model on the injection molding process of a center-gated disk with 20 wt% glass fibers [93].

#### **2.2.6.3 Slender body theory**

Models using slender body theory are closely related to Stokesian dynamics. In this classification, such models solve an integral boundary equation for fibers without physical thickness and thus without a moment induced by the fluid flow. Further, slender body theory enables the computation of the

disturbance of the fluid flow by Green's functions of the Stokes equation called Stokeslets.

For small Reynolds numbers *Re* ø 1, an incompressible Newtonian fluid, vanishing velocity at infinite distance and a singular body force field, Equation (2.17) becomes

$$0 = -\operatorname{grad}(p) + \eta \operatorname{div}(D) + \delta(x - x\_0) f\_0,\tag{2.69}$$

with a singular force *f*<sup>0</sup> in position *x*<sup>0</sup> and the Dirac distribution *±*(*x*°*x*0). The fundamental solution to this equation is given as

$$\Gamma(x - x\_0) = \frac{1}{8\pi\eta} \left( \frac{I}{\left\| x - x\_0 \right\|} + \frac{(x - x\_0) \otimes (x - x\_0)}{\left\| x - x\_0 \right\|^3} \right) \tag{2.70}$$

and called the Oseen tensor. It is used to compute the Stokeslet

$$
v = f\_0 \Gamma(x - x\_0) \tag{2.71}$$

as a solution for velocity field [94,95]. An extension of this method for singular forces is obtained by distributing several Stokeslets along a line *x*0(*s*) with arc length coordinate *s*. The solution is then given by the integral

$$w = \int\_0^L \hat{f}\_0(\mathbf{s}) \Gamma(x - x\_0(\mathbf{s})) d\mathbf{s} \tag{2.72}$$

with a force per unit length *f*ˆ 0(*s*) and is called the single-layer formulation for Stokes flow [40,95]. The linear superposition of this formulation is commonly used to represent slender, one-dimensional bodies suspended in a highly viscous flow.

Hinch developed an early model for the motion of a single perfectly flexible, but inextensible thread utilizing slender-body theory [96]. His model

provided early insight in fiber dynamics, but was not able to meet the experimental results by Forgacs and Mason [10] due to neglected width of the thread by slender-body theory.

The effect of hydrodynamic long-range interactions between rigid rods on the fiber orientation evolution was studied by Rahnama et al. [97] and showed good agreement with the results by Stover et al. [98]. Slender body theory faithfully computes the effective velocity increase up to a certain limit of fiber volume fraction, where it underestimates the effect [42] - likely due to neglected short-range lubrication interaction. Fan et al. [99] combined the slender body framework by Rahnama et al. for long-range interactions with the model for short-range lubrication by Yamane et al. They were able to determine effective rheological properties and Folgar-Tucker interaction parameters on periodic domains with rigid suspended fibers. The same method was applied to obtain an empirical fit for the Folgar-Tucker parameter depending on aspect ratio and fiber volume fraction [100]. A recent review on more advanced modeling of flexible fibers suspensions with a focus on slender-body theory is given by du Roure et al. [101].

### **2.2.6.4 Micro-macro approaches for highly concentrated suspensions**

This class of models is applied for the homogenization of highly concentrated fiber suspensions, in which long-rang interactions and viscous matrix stress may be negligible compared to short-range interactions [8, 102]. Le Corre et al. solve the momentum balance in a complex rigid fiber bundle network with viscous tangential lubrication forces and lubrication moments. The lubrication forces are formulated as

$$\mathbf{f}\_{\alpha\beta}^{\mathrm{t,l}} = \eta \frac{d\_{\mathrm{a}}^2}{\left\| \mathbf{p}\_{\alpha} \times \mathbf{p}\_{\beta} \right\|} \left( \frac{\left\| \Delta \mathbf{v}\_{a\beta} \right\|}{G \dot{\gamma}\_0} \right)^{n-1} \frac{\Delta \mathbf{v}\_{a\beta}}{G},\tag{2.73}$$

where the parameter *G* is interpreted as effective sheared gap and *d*<sup>a</sup> is the major diameter of the elliptical bundle cross section [8]. The effective sheared gap *G* is different to the physical gap *g* in direction *hÆØ*, since it describes a thickness that is equivalent to the complex shear zone between bundles and does not describe the physical distance between the bundle surfaces. The authors employ an up-scaling scheme to compute the macroscopic properties of the suspension from the mesoscale model [8, 103]. Their analysis also showed that the lubrication moment is of minor importance for most cases and confirmed the observation by Servais et al. [57] that friction is of minor importance in *bundle* suspensions. Latter was also verified experimentally [104].

### **2.2.6.5 Resolved methods**

Resolved methods are direct numerical simulations that explicitly resolve the particle-fluid interface and solve the equations of motion for these two separate phases. These methods may be further subdivided in immersed boundary methods and particle methods.

Immersed boundary methods use an Eulerian background mesh and a Lagrangian body, enabling the method to properly capture the two-way interaction between those two phases. Such methods have been applied in the field of bio-mechanics to model the whirling motion of flangella to study bacterial locomotion [105], to investigate the flow around phytoplankton employing local mesh refinements [106] and multiple wood pulp fibers [107]. However, immersed boundary methods do not scale very well considering the number of fibers and the requirement of a volume discretization that must be smaller than the fiber diameter [101].

Particle methods are mesh-less Lagrangian methods to solve partial differential equations and are interesting for direct simulations of fiber suspensions, because the phases can be naturally modeled by different types of particles. Kromkamp et al. [108] investigated the dissipation of suspended cylinders

with the Lattice Boltzmann (LB) method. The LB method was also applied to single flexible fibers and good agreement was achieved for bending modes and orientation orbits of stiff fibers [109]. Duong et al. applied dissipative particle dynamics (DPD) to Boger fluids, but several parameters have to be calibrated for correct hydrodynamic forces on fibers [110]. A moving particle semi-implicit method was applied by Yashiro et al. [111, 112] to simulate the injection molding process of dilute short-fiber-reinforced composites with connected particles. They can qualitatively reproduce molding effects, but do not verify the method and show only results with small fiber volume fractions. Smoothed particle hydrodynamics (SPH) was used by Yang et al. to model the motion of a single flapping fiber [113]. They represent the fiber by socalled element bending groups, which ensure an elastic connection between individual SPH particles. SPH was also used to model injection molding with multiple fibers [114, 115], but like Yashiro et al., the reported fiber number density of this work was unrealistically small. A more detailed analysis of SPH for fiber suspensions showed that it can be used to determine Folgar-Tucker constants on periodic domains, if a novel correction term for fiber thickness is added [116].

### **2.2.7 Concluding remarks on suspension models**

Macroscopic models utilizing fiber orientation tensors are numerically efficient models to describe dilute and semi-dilute suspensions with a fine distribution of short fibers. Such a suspension is macroscopically anisotropic due to long-range fiber-matrix interactions and short-range fiber-fiber interactions such as lubrication and contacts. There are different parameterizations for this anisotropic material response and several models are proposed to describe anisotropy based on effective material properties as well as the second and fourth order fiber orientation tensor. However, the application of macroscopic models is more suitable for injection molding process simulations, in which fibers are short in comparison to any geometrical feature.

Direct simulation models renounce the use of fiber orientation tensors and model individual fiber instances instead. The simplest approach within this class of models are kinematic models, where fibers are advected directly with a prescribed fluid velocity in combination with a regularization scheme to maintain a constant fiber length. These models are computationally fast, but do not affect the fluid flow and do not account for interactions between fibers. Resolved computational models with fully modeled individual matrix phase and fiber phase offer the potential to completely describe a suspension by immersed boundary methods or particle methods. However, these models are computationally expensive and limited to small representative volume elements. Slender body theory offers a computationally efficient approach to solve fluid flow and fiber motion simultaneously and is commonly used in the bio-mechanical community. While the method is very suitable for representative volumes with prescribed far-field boundary conditions, it can not be applied easily to component scale compression molding simulations with a deforming cavity and contacts at the mold surface. Similarly, the presented micro-macro up-scaling approaches are restricted to representative volume elements. The most common approach for fiber suspension modeling is Stokesian dynamics which solves the force balance for interlinked suspended bodies due to hydrodynamic forces and non-hydrodynamic interaction forces. Within this class of models, the work by Lindström and Uesaka [85–87] is particularly noteworthy, as it affects the fluid flow to be anisotropic through two-way coupling.

## **2.3 Sheet Molding Compounds**

Sheet Molding Compound (SMC) is a discontinuous fiber reinforced polymer made from uncured thermoset sheets and can be molded to parts under pressure and temperature by compression molding technology. Parts from this material offer a superior strength to density ratio (compared to polymers without reinforcement) at low cost. Furthermore, corrosion resistance and high achievable surface quality make it attractive for outer vehicle parts, such as truck body panels, automotive hoods and trunks. While continuously reinforced polymers have limited forming capabilities due to the inextensibility of fibers in the fabric, SMCs can flow in complex geometries during compression molding. The incorporation of complex geometrical features, carbon fibers, metallic inserts and continuous fiber patches enables the application in structural applications such as subfloor structures in the automotive industry.

## **2.3.1 Production**

The production of SMC parts is subdivided in two major steps: the production of SMC sheets from raw materials and the subsequent compression molding of sheets to a part. The sheet production (see Figure 2.8a )starts with mixing of all ingredients (thermoset resin, filler, thickener, additives, catalyst, mold release agent, inhibitor) except for fibers to prepare the paste that later forms the matrix material of the composite. The paste is then filled to doctor boxes of an SMC sheet production line, which apply the paste on two different plastic foils that run continuously through the machine. The boxes have adjustable plates to ensure a specific paste film thickness on the foil.

At the same time, multi-end fiber rovings are pulled from several bobbins to a cutting unit that continuously cuts the strands in 25 mm or 50 mm long pieces. These bundles fall down on the paste layer of one of the plastic foils in a randomly in-plane orientated manner. The second resin-coated foil is placed on top of the first foil to enclose the resin and fiber bundles between both foils. The compound is then pulled through several (heated) calendering rolls to impregnate the fiber bundles with the thermoset paste and compact the compound.

Finally, the sheet is rolled to a take-up coil and stored for several days and the resin starts a maturing process. During this maturing phase, the viscosity increases several orders of magnitude (compare Figure 2.8b).

**Figure 2.8:** SMC sheet production and maturing phase. The viscosity plot is recorded at a shear rate *<sup>∞</sup>*˙ *<sup>=</sup>* 1.0s°<sup>1</sup> [117].

The sheets are cut, stacked and placed into a hot mold for compression molding. The stack typically covers a fraction of 25% to 100% of the mold and the remaining cavity is filled by flow of SMC, as the press closes the upper mold stamp in the corresponding lower mold half. The compression follows a defined closing velocity profile until the compression force hits a threshold at which the press controller switches to force control. This compression force level is held for a few minutes until the resin is sufficiently cured. Finally, the press opens the mold and the cured part is released.

Possible defects according to Mallick [118] are summarized in Figure 2.9. If air is introduced in the paste during mixing or compounding and cannot be removed by compression between calendering rolls or if air is entrapped between sheets during stacking, pores and pinholes at the surface may form. Entrapped air or gaseous products of the curing reaction, which are closed during the compression may also form blisters that pop up as soon as the mold is opened. Sink marks may occur in resin rich areas due to shrinkage of the resin during curing. Typically, SMC is highly filled with chalk to prevent such surface defects.

Fiber bundles reorient during the mold filling process and this may lead to unwanted fiber miss-orientation. For example, fibers close to mold edges

Shear viscosity at 1 s−<sup>1</sup> shear rate in Pa s

**Figure 2.9:** Possible defects in parts molded from SMC according to Mallick [118].

have a tendency to be oriented parallel to the edge. Additionally, knit lines may form where SMC has to travel around a feature. The line at which two flow fronts of SMC meet is typically weak, as only a limited amount of fibers crosses this knit line. If the paste travels faster than the fiber bundle network during molding, fiber-matrix separation (FMS) occurs. This phenomenon happens often in small ribs, at the flow front and other small features, where fiber bundles may get stuck. Fiber-matrix separation results in a certain area of a part consisting of matrix material only, which has significantly worse properties compared to the composite.

### **2.3.2 Deformation mechanisms**

The micro-structure of SMC sheets implies certain restrictions on the deformation kinematics, which are significantly different to the deformation of an isotropic fluid or solid. The initial stack is characterized by a random in-plane fiber orientation and the thickness is typically small compared to the lateral dimensions of the stack and the fiber length. Marker and Ford [119] present one of the earliest investigations on the deformation mechanisms and heating of SMC sheets. They investigate large stacks of ten sheets thickness in which the interlocking behavior of individual sheets resists transverse

shear deformation and outer sheets close to the mold wall experience larger deformation. Barone and Caulk [120] use this observation to simplify the kinematics of thin stacks (one or two sheets) by proposing that no shear velocity gradient exists along the thickness direction of the stack. This assumption was declared unrealistic by Tucker and Folgar [121], who suggested that a no-slip condition at the surface must apply. Subsequently, Barone and Caulk [122] found experimental evidence for their proposed stack deformation mechanism using a disk shaped initial stack with black and white sheets of SMC. They performed compression experiments at two mold closing speeds (1.75 mm s−<sup>1</sup> and 10.0 mm s<sup>−</sup>1) with multiple stack sizes. Exemplary results for a 5-layer stack are shown in Figure 2.10a and 2.10b. Barone and Caulk concluded that the deformation mechanism is primarily a uniform extension that can be accompanied by slip between layers and mold walls for large stacks and slow deformation speeds. Contrary to the *fountain-flow* observed in thermoplastic compression molding [123], this deformation mode is termed *plug-flow*. Later, Barone and Caulk modeled the slip at the mold surface with a constant friction, Coulomb friction and hydrodynamic friction and concluded that the hydrodynamic friction yields most realistic results for the SMC, they investigated [124].

Later investigations of the stack deformation differentiate a *squish*, *flow* and *boiling* phase during molding [125]. The squish phase refers to the initial complex flow front deformation due to thermal influences and release of air entrapped in the stack. Exemplary flow front shapes are illustrated in Figure 2.10c for 10 mm s−<sup>1</sup> and in Figure 2.10d for 2 mm s−<sup>1</sup> mold closing speed. Especially in the case of faster compression, it becomes apparent that the bottom layer exhibits larger deformation upon start of the compression, as it is hotter due to a longer contact with the mold surface. The subsequent flow phase is characterized by a stable plug-flow and the final boiling phase refers to gas release after complete closing of the mold. A minimization of the squish phase yields to a more homogeneous flow and reduced the volume fraction of pores in the SMC part [126]. Fiber bundles in the SMC core stay

typically intact during molding and only fiber bundles in the lubrication layer near the mold surface disentangle during the flow [117, 127–129].

## **2.3.3 One-phase modeling approaches**

### **2.3.3.1 Generalized Hele-Shaw models**

Hele-Shaw models simplify the incompressible Stokes flow between two plates to a two-dimensional problem, if the lateral dimensions are large compared to the gap between theses plates. They can be applied to model the compression molding process between two flat plates with gap *h*, that close with velocity *h*˙ . The pressure field in such a 2D-domain is given by

$$\text{div}\{\text{S grad}(p)\} = \dot{h} \tag{2.74}$$

where the divergence and gradient operators refer to the 2D-domain and

$$\mathbf{S} = \int\_{-h/2}^{h/2} \frac{(z - z\_0)^2}{\eta} \mathbf{d}z \tag{2.75}$$

is a measure of mobility [130]. Here, the position *z*<sup>0</sup> refers to the position of the flow, where d*vx* /d*z =* d*vy* /d*z =* 0. The corresponding 2D-velocity field is obtained as

$$
\sigma = -\frac{S}{h}\,\mathrm{grad}(p). \tag{2.76}
$$

For the special case of isotropic Newtonian material behavior and no-slip conditions at the mold walls, *<sup>S</sup> <sup>=</sup> <sup>h</sup>*3/12*¥* applies and *<sup>z</sup>*<sup>0</sup> *<sup>=</sup>* 0, as depicted in Figure 2.11a. These assumptions were used in the earliest models of SMC compression molding [131]. However, the thermally induced viscosity change may lead to failure of the parallel squeezing assumption and more complex velocity profiles develop, as the viscosity near the mold surfaces decreases with increasing temperature, as shown in Figure 2.11b [132]. A dimensionless criterion to determine whether a parallel squeezing assumption is valid was developed by Lee et al. [133].

**(a)** Hele-Shaw flow **(b)** Thermal effect in SMC **(c)** Plug-flow assumption

**Figure 2.11:** Velocity distributions through thickness for different mold boundary conditions and thermal conditions.

The Hele-Shaw model for compression molding was solved for arbitrary domains beyond simple analytical cases using the finite element method [121]. The authors indicated the incorporation of the energy balance and Non-Newtonian viscosity in their simulation framework.

Generalized Hele-Shaw models utilize the simplification for the flow between arbitrary curvilinear surfaces with a thin gap between them. They are not restricted to flat plates and the application to compression molding was shown by Lee et al. [130]. The proposed model included heat transfer and isotropic Non-Newtonian viscosity and was solved using finite elements. The comparison to experimental results showed a good agreement for thin charges, but differences for thick charges due to a neglection of extensional stresses. The model was later combined with tensorial orientation models to predict the orientation after compression molding [134]. For a Newtonian material with constant part thickness, the Generalized Hele-Shaw flow may be solved using the boundary element method, which is very efficient and does not require a complicated meshing procedure [135].

### **2.3.3.2 Plug-flow models**

The Generalized Hele-Shaw models employ a no-slip boundary condition at mold walls, which does not comply with the observed SMC deformation behavior of Barone and Caulk [122]. Thus, the boundary element methods were extended to account for the modified flow behavior depicted in Figure 2.11c with slip and hydrodynamic friction at the walls that maintains a velocity profile without bulk shear deformation [136, 137]. Later, non-isothermal modeling (similar to the illustration in Figure 2.11b) confirmed the validity of a flat velocity profile with a thin hot lubrication layer close to the walls [138]. This settled the dispute between Lee and Tucker (advocates of the no-slip model) and Barone and Caulk (advocates of the plug-flow model) in favor of the plug-flow assumption.

As the simplest approach, the hydrodynamic friction can be described by a hydrodynamic friction coefficient *∏*, i.e. the proportionality constant between relative velocity and friction force. For sufficiently thin parts, that is

*<sup>h</sup>¥*/(*B*2*∏*) <sup>ø</sup> 1, friction dominates the momentum equation. In this case, the mobility term in Equation (2.74) can be replaced with

$$S = \frac{h^2}{2\lambda} \tag{2.77}$$

to use an equivalent form of the generalized Hele-Shaw model with a more accurate description of the flow kinematics [139–141]. Such a model was successfully applied to the simulation of complex automotive structures as early as 1990 using finite elements [142]. The boundary element formulation requires too many simplifications for the application to such complex problems. This approach was extended for non-isothermal conditions later on [143]. Kotsikos et al. [144] employed a variational approach to allow shear and extensional deformation components and compared the results to experimental data from squeeze flow compression tests. They confirmed that SMC flow is predominantly of an extensional nature.

After the establishment of basic models for the compression molding simulation, the research focus switched to characterization of the material parameters in order to predict molding forces correctly. Rheometry of SMC is difficult as the material has to be considered inhomogeneous at the scale of typical rheometers (such as plate-plate rheometers). Hence, custom press rheometers were developed to compress a sufficiently large sample and identify rheological parameters with underlying modeling assumptions. Le Corre et al. performed lubricated (i.e. a lubrication film was applied to the mold surface to eliminate friction between SMC and the molds) disk compression experiments and shear experiments with temperatures, deformation rates and fiber volume fractions similar to industrial applications [145]. They showed that the SMC behaves as power-law fluid, if the underlying paste is a power-law fluid and that the material is highly anisotropic. They identify the need to model SMC as a transversely isotropic material. These results were corroborated by plane strain tests to measure the transverse force and a set of constitutive equations was suggested to describe non-linear anisotropic SMC viscosity [146]. First parameters for non-lubricated flow, especially the

hydrodynamic friction parameter, were obtained by a comparison of the compression force between different initial positions of the stack in rectilinear compression [147]. The hydrodynamic friction can also be characterized more accurately by employing pressure sensors along the flow path and it was shown that both, friction and bulk rheology, are relevant to model correct compression forces [148]. Spiral rheometers with pressure sensors along the flow path were suggested as compact alternative to rectilinear press rheometers [149]. However, they may be considered as a rather qualitative method to compare SMC processing capabilities than as an exact tool for characterizing friction parameters.

Several authors argue that the SMC mold filling process can be considered isothermal except for the heating of the small lubrication layer at mold surfaces [148, 150]. Shokrieh et al. [151, 152] disagree with this assumption and suggest a time dependent thickness of the lubrication layer that is controlled by the development of the temperature in the stack and requires no additional experimental parameters.

An additional characterization technique is tensile testing. It enables the characterization of damage and breakage in cases where sheets of SMC are drawn to deep cavities instead of being driven by a squeezing motion between mold surfaces [153].

While most previous models assume SMC to be incompressible, 3D in-situ computer tomography showed that the closing of pores that originate from sheet manufacturing cause a macroscopic compressible behavior [154]. This effect is especially true for structural SMCs and Dumont's equations [146,148] have been modified to account for this compressibility [150, 155].

### **2.3.3.3 3D Models**

Most early simulation approaches simplify the SMC mold filling as a 2D process. This is appropriate for large shell-like structures, but has limitations, when it comes to ribs or other complex 3D features in SMC parts, which are

increasingly important for structural SMC applications. An initial attempt to model SMC mold filling with a 3D CFD model was done by Vahlund [156], even though the method was applied to a shell-like part. Kluge et al. [157] utilized a 3D simulation to accompany their rheological experiments and gain deeper understanding of the material response for different viscosity models. Motaghi et al. [129] used the commercial process simulation software Autodesk Moldflow, and employed an RSC model for fiber orientation. However, they used an isotropic viscosity model and did not compare corresponding compression forces. Hohberg et al. [158] pointed out that boundary conditions used in Autodesk Moldflow at that time do not respect the plug-flow assumption and suggested an elongational viscosity model with appropriate boundary conditions in Simulia Abaqus/Explicit. The model showed potential, but suffered from numerical instabilities. As of 2021, the latest version of Autodesk Moldflow allows the application of wall slip models at all mold surfaces1, which makes Autodesk Moldflow a viable software to predict SMC compression molding at macroscopic scale. The ability to account for wall slip is also available in CoreTech System Moldex3D and has been investigated for thermoset injection molding [159].

### **2.3.3.4 Thermoviscoplastic models**

Initial models for the SMC response during disk compression also suggested the use of a visco-elasto-plastic model [160] instead of the predominantly purely viscous models suggested above. However, this initial work did not account for the plug-flow kinematics described later by Barone and Caulk [122]. Kim et al. [161] suggested a modeling approach similar to that of sheet metal forming describing the SMC as solid anisotropic material with Hill's yield criterion. They characterized parameters from lubricated disk compression tests and validated their model by comparing a finite element

<sup>1</sup> The required improvement was specified by Sven Revfi and Nils Meyer. Sven Revfi suggested the improvement to Autodesk and the issue was resolved in the new release.

simulation with experimental data of filling a T-shaped rib. In contrast to 2D shell simulations, this approach modeled a 2D cross section of the rib and was able to predict the filling behavior. Later, the authors enhanced their model by incorporating heat transfer, a curing model and a simple approach to predict the resulting fiber volume fraction [162] and extended the model to 3D finite elements [163]. The numerical results showed good agreement to molding trials of T-shaped parts, that were molded with multiple colored SMC layers. The characterization of the parameters for the anisotropic flow rules is described by Lin et al. [164, 165]. In most SMC applications however, the elastic component is rather small and the plastic deformation has little effect on the material response. Therefore, the majority of SMC models assumes fluid-like behavior, as described in the previous sections.

## **2.3.4 Two-phase models**

Most SMC models consider the material as a one-phase continuum, as described in the previous Section 2.3.3. However, segregation of fibers and matrix is a relevant phenomenon in SMC compression molding, particularly in advanced applications with ribs, knit line formation and long flow paths. It motivates the development of models that can account for different velocities of fibers and matrix material. Besides the direct numerical models described in Section 2.2.6, this has been addressed by macroscopic SMC models which are shortly described in this section.

Dumont et al. suggested a superposition of two immiscible phases (fiber and matrix) occupying the same volume at isothermal and incompressible conditions [166]. They obtain two momentum balance equations from mixture theory with a viscous momentum exchange term and solve the equations on a shell finite element mesh. The momentum exchange term is chosen such that it reduces to the equations of flow through porous media, if fibers are held in place rigidly. A similar approach was developed by Perez et al. resulting in a general 3D Brinkmann model [167]. However, the weighting factor between the two extrema (flow through porous media and suspension

flow with equal velocities between fiber and fluid) had to be prescribed and a lubrication approximation (compare Figure 2.11a) was used. Later, they related the weighting factor to the fiber-fiber contact density to compute the fiber volume fraction evolution [168].

## **2.4 Unidirectional patch reinforcements for co-molding**

Unidirectional fiber reinforced composites offer superior weight specific mechanical performance in fiber direction compared to discontinuous FRPs and most metals. They are either manufactured by thermoplastic tape laying combined with thermoforming, by consolidation of non crimp fabrics that are pre-impregnated with a thermoset matrix (prepregs) or by infiltration of dry fabrics (RTM, WCM, VARI).

However, these materials are rather expensive and the ability to form complex shapes is limited by the inextensibility of fibers. The combination of discontinuous SMC with local unidirectional reinforced patches in a co-molding process offers potential to benefit from advantages of each constituent [3]. Thus, this section focuses on unidirectional pre-impregnated thermoset patches that can be co-molded with SMC, as developed by David Bücheler in the GRK 2078 [169, 170].

## **2.4.1 Production**

The production of UD-patches for co-molding with SMC is similar to the production of SMC sheets and can be performed on the same production line (see Figure 2.12a). Instead of a cutting unit, the fabric is placed in front of the compounding zone and is impregnated with the matrix material between top and bottom foils. For the application in co-molding, the material is exposed to a special heat treatment after the impregnation to rapidly develop a highly viscous resin B-stage state in the material. The B-staged material has experienced polymer chain growth, but no cross linking of the thermoset. This treatment allows direct cutting and forming after the impregnation process and ensures enough stability of patches during molding [170].

**Figure 2.12:** Patch production and co-molding process.

Patches are stacked while still being tacky directly after impregnation. Subsequently, they are cut, preformed and heated. After 60 min at 80 °C, the patches are in a stable shape close to its final geometry in the co-molded part and are released from the preforming mold. The UD-patches are molded together with SMC made from the same matrix resin in a regular compression molding mold (see Figure 2.12b). During compression molding, patches may be subjected to the defects illustrated in Figure 2.13. Displacement occurs if the SMC excites a net force on the patch that is larger than the friction between patch and mold. The SMC acts upon the patch with normal forces and a sticky tangential friction. Local deformation is caused by a combination of insufficient friction to hold the patch in place and a distributed load by the SMC causing permanent deformation of the patch. If this deformation exceeds the strength of the patch in a particular direction, fracture may occur. This occurs typically perpendicular to the fiber direction by matrix tension or shear, as fibers have a high tensile strength.

**Figure 2.13:** Possible defects in unidirectional patches during co-molding (top view).

## **2.4.2 Deformation mechanisms**

The quasi-inextensibility of fibers and comparably low shear stiffness implies kinematic constraints on the deformation mechanisms of continuous fiber reinforced fabrics. These deformation mechanisms are usually classified as intra-ply mechanisms happening in a single patch layer and inter-ply mechanisms occurring between individual layers of patches. The intra-ply mechanisms can be further categorized in bending behavior, membrane behavior and in-plane shear behavior [171]. Unlike homogeneous isotropic materials, bending stiffness and membrane stiffness of engineering textiles may require a decoupled modeling approach, if the polymer matrix stiffness is extremely small and allows slip between individual fibers during bending load [171]. Inter-ply mechanisms describe slip between plies, slip between mold and plies as well as delamination. The friction between mold and patch has been previously investigated by Bücheler [169] for the material modeled in this work.

The predominant deformation mode of biaxial fabrics is called *pure shear* or *Trellis shear* that is motivated by inextensible fibers in two transverse directions with contacts acting as pivot points. In contrast, UD fabrics naturally deform under so-called *simple shear* assuming inextensibility in fiber direction and incompressibility. The simple shear deformation is depicted in Figure 2.14.

The deformation gradient *F* maps the initial configuration to the current configuration, as indicated in Figure 2.14 by black coordinate systems with and without superscript '0', respectively. The coordinate system, which is transformed by *F* is called the *covariant* reference system. However, the

**Figure 2.14:** Simple shear in a unidirectional patch. The covariant coordinate system is depicted in black (*e*x,*e*y) and the co-rotational system in green ( ˆ*e*x, ˆ*e*y).

deformation gradient includes a rigid body rotation and is therefore not objective. A polar decomposition

$$F = RU = VR\tag{2.78}$$

in a pure rotation *R* and a right stretch tensor *U* or left stretch tensor *V* provides a separation in rigid body rotation and an objective strain measure. However, the by *R* rotated co-rotational frame is typically not suited to describe the deformation of anisotropic composites, because the principal material orientations become non-orthogonal during shearing, as shown with the black coordinate system in Figure 2.14 in comparison to the orthogonal frame depicted as green coordinate system. A possible solution to cope with this is the introduction of additional transformations in a hypoelastic constitutive model (see Section 2.4.3.1) or the formulation of a hyperelastic constitutive model (see Section 2.4.3.2) in the initial configuration. In any way, the constitutive model needs a deformation measure that does not depend on rigid body motion and rigid body rotation (objective strain measure). The definition of the right Cauchy-Green tensor

$$\mathbf{C} = \mathbf{F}^{\top} \cdot \mathbf{F} = \mathbf{U}^{2} \tag{2.79}$$

allows the definition of the Green-Lagrange strain

$$E = \frac{1}{2}(C - I),\tag{2.80}$$

which is a suitable objective deformation measure for modeling the patches.

### **2.4.3 Modeling approaches**

The simplest approach to model patch deformation are kinematic models that assume incompressibility and inextensibility of fibers [172]. In biaxial fabrics the contacts between weft and warp fibers act as pivot points without slippage and in unidirectional reinforcements ideal simple shear with constant fiber spacing is assumed. The forming result is purely controlled by geometrical considerations and without solving a momentum-balance (similar to the kinematic fiber models presented in Section 2.2.6.1). The benefit of such models is the fast computation time, which can be leveraged in optimizations [173], but the prediction capability of defects is very limited [174].

A more rigorous approach is the constitutive modeling of fabric deformation that accounts for material behavior. Such models are developed for the microscale to analyze slippage of fibers and resulting yarn cross sections during weaving [175]. Detailed mesoscale models [176, 177] and discrete models build from 1D or 2D elements representing constituents of the composite [178–180] are used to determine effective properties based on the underlying yarn structure or directly model deformation at component scale. Most commonly, forming behavior of engineering textiles is described by homogeneous macroscale models. This is often done either by formulations that relate the rate of the Cauchy stress to the strain rate (hypoelastic) or by formulating the material behavior with respect to the initial configuration (hyperelastic) [181, 182].

#### **2.4.3.1 Hypoelastic models in co-rotational frame**

Hypoelastic models follow the general form

$$
\sigma^{\nabla} = \mathbb{C} : \mathbf{D} \tag{2.81}
$$

with a tangent stiffness tensor C and an objective stress derivative r, such as the Jaumann rate or Green-Naghdi rate [183]. The constitutive equations are formulated in terms of the current Cauchy stress, which is integrated in time by an appropriate numerical scheme [184]. In general, fiber reinforcement fabrics are described in a non-orthogonal frame and a transformation between the orthogonal objective Cauchy rate (see *e*ˆ*<sup>x</sup>* ,*e*ˆ*<sup>y</sup>* in Figure 2.14 for Green-Naghdi frame) and the non-orthogonal fiber frame (see *e<sup>x</sup>* ,*e<sup>y</sup>* in Figure 2.14) has to be applied [185–187].

#### **2.4.3.2 Hyperelastic models in initial frame**

Hyperelastic models allow for a formulation of constitutive material models with respect to the initial configuration and thus omit the need for difficult coordinate transformations. Such models propose the existence of a strain energy potential *W* and compute the second Piola-Kirchhoff stress as

$$S = 2\frac{\partial W(C)}{\partial C},\tag{2.82}$$

where *C* is the right Cauchy-Green tensor [183]. The strain energy potential is typically formulated in terms of invariants of the right Cauchy-Green tensor and includes energy stored by tension and shear [188,189] as well as coupling terms [190, 191].

The Cauchy stress can be computed by a *push forward transformation* with the deformation gradient as

$$
\sigma = \frac{1}{J} \boldsymbol{F} \cdot \boldsymbol{S} \cdot \boldsymbol{F}^{\top}, \tag{2.83}
$$

where the Jacobian determinant *J =* det(*F*) denotes the volume change [183].

Given the non-orthogonal deformation behavior illustrated in Figure 2.14, a hyperelastic formulation w.r.t the initial configuration seems most appropriate to model the deformation of unidirectional patches during co-molding.

## **3 Objective and outline**

## **3.1 Research hypothesis**

Fibers in SMC are typically 25 mm to 50 mm long and are required to fill geometric features of similar dimensions in advanced (semi-)structural applications. This violates several conditions for the application of tensor based fiber orientation descriptors, e.g. fibers are not straight, fibers are not much shorter than geometric features and the effects are not local, because the motion of one fiber end affects the other end up to 50 mm away. The question is in which scenarios these violations may still lead to reasonable results and which alternative model could be used to alleviate these issues while still being applicable at component scale with reasonable computational times.

The main hypothesis is that a mesoscale process model that accounts for the motion of individual fiber bundles improves the prediction quality of fiber architecture in SMC compression molding simulations at component scale. The improved prediction of the architecture can be qualitatively and quantitatively validated by comparison to µCT scans.

## **3.2 Objective**

The objective of this work is the development of a 3D numerical simulation method for the compression molding process of SMC in regions, where fiber orientation tensors may not be suitable descriptors of the fiber orientation state. These regions are geometrical features such as narrow ribs and beads which are common in advanced structural SMC applications as well as the area near mold walls. The goal is the prediction of structural defects during SMC compression molding, such as knit line formation, fiber-matrix separation and fiber misalignment. Additionally, accurate predictions of necessary compression forces and fiber orientation states are desired. Further, the simulation method should be able to incorporate a model for continuous fiber patches to predict processing effects in parts made from CoDiCoFRP.

The process step under consideration is the mold filling process starting from positioning of the initial SMC stack and patches until complete fill of the mold. Preforming of the material placed in the mold is not considered, nor is subsequent curing, cooling and warpage of the part.

The detailed sub-objectives are summarized as follows:


## **3.3 Outline**

The work is laid out as follows: The next Chapter 4 formulates a 3D macroscopic reference model for SMC employing fiber orientation tensors and anisotropic viscosity. It is solved using the Coupled Eulerian Lagrangian framework in Simulia Abaqus/Explicit and a custom VUMAT subroutine. The equations are reduced to the 1D case in a press rheometer to investigate principal effects of individual parameters and provide a tool for rapid computation of reference results.

Chapter 5 includes the main contribution of this work, which is the development of a mesoscale model that is capable of running at component scale to provide more detailed prediction of defects and features during compression molding. The chapter is subdivided in the derivation of long-range bundlematrix interactions and short-range bundle-bundle interactions. Sub-models are used in each case to derive phenomenological expressions for relevant terms at mesoscale, such as the drag coefficients of bundles and the effective sheared gap between bundles. Finally, a verification ensures that the model is able to reproduce analytical reference solutions.

A model for unidirectional reinforcement patches is described in Chapter 6. The model describes each patch ply with two shell layers with shared nodes. One layer represents the unidirectional bundles and the matrix, the other layer models the unidirectional stitching yarn. Both layers are hyperviscoelastic and a Hashin damage criterion is employed to describe failure of patches during molding.

Chapter 7 describes the characterization of necessary parameters for the SMC model and the patch model in the uncured state. This includes the transverse heat conductivity of SMC sheets and patch plies as well as the heat conductance from a steel mold to the SMC stack and patch. Further, the rate dependent and temperature dependent viscosity of the SMC paste and the compressibility of a uncured stack of SMC are determined. For patches, the relevant anisotropic visco-elastic properties are determined.

In Chapter 8, several experimental results and corresponding simulation results are compared to analyze the merits of direct mesoscale simulations and macroscale models. The first case investigates the general ability of the DBS to predict fiber architectures under quasi-incompressible isothermal Newtonian conditions in a double curved dome. The resulting fiber orientation, curvature and knit line formation agree well with computer tomography scans of experimental specimens. The second case investigates fiber-matrix separation at ribs with quasi-incompressible isothermal Newtonian modeling conditions and compares the results to pyrolysis results of specimens. The full non-isothermal non-Newtonian compressible model is applied to a press rheometer simulation in the next example. The resulting compression force and pressure at sensor positions is compared to experimental data to evaluate the model's capability to not only predict fiber architectures, but also to predict realistic compression forces. The full model is also applied to a more complex small demonstrator part featuring ribs, beads and unidirectional reinforcement patches to show the ability to integrate the patch model with DBS. The final validation example is a large part to demonstrate that the suggested approach can be applied at the scale of larger SMC components.

The modeling domains and investigated features are visually summarized in Figure 3.1.

**Figure 3.1:** Graphical overview of the modeling domains and corresponding simulation features. The molds are represented by rigid shells (R3D3/R3D4) at a homogeneous isothermal temperature. The motion is controlled by a virtual press controller that mimics the controller of the physical press with a switch from displacement control to force control. A key feature is the explicit model of fiber bundles by flexible thread-like truss elements (T3D2) at component scale and with realistic fiber volume fractions. The fiber motion is governed by drag forces from the fluid phase and fiber bundles interact with each other as well as with the mold surfaces. The fluid phase is modeled in an Eulerian frame (EC3D8RT) and is generally non-Newtonian, non-isothermal and compressible. The drag of fiber bundles is subjecting the fluid phase to a body force field that causes anisotropic flow behavior. The simulation framework may incorporate continuously reinforced patches modeled as shells (S3RT/S4RT), which are modeled as hyperviscoelastic material with a Hashin criterion for damage initialization.

## **4 Macroscopic reference model**

## **4.1 General three-dimensional model**

Typically, SMC process simulation is conducted at a macroscopic continuum scale, which efficiently describes the flow of the SMC fiber suspension as a homogeneous material with effective material properties. To verify the mesoscale model developed subsequently in Chapter 5, such a macroscopic model is formulated here. It is implemented in the commercial finite element software Simulia Abaqus/Explicit and provides a reference solution to discuss differences in the modeling scales.

### **4.1.1 Assumptions and scope**

The homogeneous description of SMC implies some assumptions, which are introduced to simplify the chemo-thermo-elasto-visco-plastic behavior of SMC:


The first two assumptions are reasonable for simple plate-like structures with in-plane dimensions that are much larger than fiber length. However, these assumptions are highly questionable in confined regions, knit lines and any feature at fiber length scale. The third assumption is founded in the hypothesis that flow occurs on a time scale much faster than curing (10<sup>0</sup> s vs. 10<sup>2</sup> s). Thus, curing is neglected to reduce the number of parameters in the model and characterization efforts. However, even a small onset of curing might affect the viscosity significantly and thus introduce errors in the simulation results. The fourth assumption reduces the thermo-elasto-visco-plastic material behavior to a thermo-viscous behavior. There is certainly some elasticity in the composite (fibers are elastic and the polymer matrix induces entropy-elasticity as well as some early degree of cross-linking). Even so, the elastically stored energy in the SMC is assumed to be small compared to the dissipated viscous energy during the molding process. Not all dissipated energy is necessarily rate dependent. Internal friction, such as between fibers, causes some degree of plastic behavior. Contrary to fiber suspensions with dispersed bundles, which show a Herschel-Bulkley behavior with yield stress, fiber bundle suspensions are dominated by viscous fiber interactions [56, 57]. Hence, plasticity is neglected here for simplicity.

### **4.1.2 Governing equations**

The initial boundary value problem is formulated on a domain *<sup>x</sup>* <sup>2</sup> *<sup>D</sup>* <sup>Ω</sup> <sup>R</sup>3. The solution fields are the velocity *<sup>v</sup>* : *<sup>D</sup>* ! <sup>R</sup>3, the mass density *<sup>Ω</sup>* : *<sup>D</sup>* ! R, the temperature *T* : *D* ! R and the inner variables describing the fiber orientation state *A* : *D* ! *A*. Considering symmetries, this leads to a total of eleven independent variables © *v*1, *v*2, *v*3,*Ω*,*T*, *A*11, *A*22, *A*33, *A*12, *A*23, *A*13™ .

The solution is obtained from the macroscopic mass balance

$$\frac{\partial \rho}{\partial t} + \text{div}(\rho v) = 0,\tag{4.1}$$

the macroscopic momentum balance

$$\frac{\partial(\rho v)}{\partial t} + \text{div}\left((\rho v) \otimes v\right) = \text{div}(\sigma),\tag{4.2}$$

the macroscopic balance of inner energy

$$\frac{\partial(\rho c\_{\rm p}T)}{\partial t} + \text{div}\Big((\rho c\_{\rm p}T)\cdot v\Big) = -\text{div}(d) + \sigma :D,\tag{4.3}$$

and an evolution equation for the second order fiber orientation tensor

$$\frac{\partial \mathbf{A}}{\partial t} + \text{div}\,(\mathbf{A} \otimes v) = \mathbf{W} \cdot \mathbf{A} - A \cdot \mathbf{W} + \lambda \left\{ D \cdot A + A \cdot D - 2 \mathbb{A} : D \right\}.\tag{4.4}$$

The employed orientation model is the simplest possible version (see Section 2.2.4) without diffusive interaction parameters. This could be easily replaced by a more complex orientation model, if the required additional parameters are available. However, it is shown later in this work that the simple base model shows remarkable agreement with the detailed mesoscale model for planar SMC compression molding.

The governing equations result in a total number of eleven equations for the unknowns. The heat flux *d* and momentum flux are related to the solution variables via the constitutive equations outlined in the next section to close the problem.

### **4.1.3 Constitutive equations**

The Cauchy stress is computed according to Equation (2.33) from several individual contributions as

$$
\sigma = -p(\Delta \varepsilon)I + \sigma\_{\text{M}}^{\prime} + \sigma\_{\text{FM}}^{\prime} + \sigma\_{\text{FF}}^{\prime},\tag{4.5}
$$

where the pressure *p* is no individual solution variable, but modeled with an equation of state that depends on the current mass density through the Hencky strain

$$
\Delta \varepsilon = \log \left( \frac{V}{V\_0} \right) = \log \left( \frac{\rho\_0}{\rho} \right). \tag{4.6}
$$

A relation between mass density and pressure is determined in Chapter 7. Alternatively, a sufficiently large bulk modulus may be employed to model quasi-incompressible SMC. The isotropic deviatoric matrix stress

$$
\sigma'\_{\mathcal{M}} = 2\eta(T, \dot{\gamma}) : \mathcal{D}' \tag{4.7}
$$

uses a matrix viscosity that generally depends on shear rate *∞*˙ and temperature *T* . The extra stress from long-range fiber-matrix interactions is

$$\begin{split} \sigma\_{\rm FM}^{\prime} &= \frac{4f r\_{\rm p}^{2} \eta(T, \dot{\gamma})}{3 \left[ \ln(1/f) + \ln \ln(1/f) + C \right]} \bigg( \mathbb{A} - \frac{1}{3} I \otimes A \Big) : \mathcal{D}^{\prime} \\ &+ 2n\_{\rm s} \eta(T, \dot{\gamma}) \left( A \Box I + I \Box A - \frac{2}{3} I \otimes A \right) : \mathcal{D}^{\prime}, \end{split} \tag{4.8}$$

as given in Equation (2.39) in combination with the parameter *¥*2, SF from Equation (2.42). Strictly speaking, the resulting Equation (4.8) is only valid for Newtonian fluids with constant viscosity due to the assumptions made in the underlying derivations. It may be applicable, if the non-linearity introduced by *¥*(*T*,*∞*˙) is not too severe (approximately constant local viscosity at the length scale of fibers), but results have to be reviewed critically when using Non-Newtonian viscosity. The shear factor *n*<sup>s</sup> is included to improve numerical stability in 3D models and might be considered a numerical parameter that is sufficiently small to not affect the result but large enough to suppress spurious shear. The last term in Equation (4.5) was obtained by projecting Equation (2.59) with the deviatoric projector P<sup>2</sup> to ensure that it does

not generate volumetric stresses. Hence, this extra stress from short-range fiber-fiber interactions is formulated as

$$
\sigma\_{\rm FF}^{\prime} = \frac{4k\_{\rm D}}{3\pi^2} r\_{\rm p}^2 \eta(T, \dot{\gamma}) f^2 \left( \mathbb{B} - \frac{1}{3} I \otimes B \right) : \mathcal{D}^{\prime}. \tag{4.9}
$$

The thermal flux for the energy balance in Equation (4.3) is modeled according to Fourier's law

$$d = -\kappa \,\mathrm{grad}(T)\tag{4.10}$$

where represents the heat conductivity tensor. The conductivity is reduced to a scalar value *= ∑I* here, to reduce the number of parameters to be characterized. This scalar model is a simplification, as in-plane heat conductivity is different to the transversal conductivity. However, the conduction process in SMC is controlled by the heating between hot molds transverse to the sheets, while the majority of in-plane heat transport is of convective nature due to the material flow.

Finally, the fourth order orientation tensor A is computed from the IBOF closure [21] and the fourth order interaction tensor B is computed from the quadratic closure provided in Equation (2.63). A summary of the parameters for the macroscopic model is given in Table 4.1.


**Table 4.1:** Parameters for the macroscopic SMC model. Additional parameters are needed to describe the relations for *¥*(*T*,*∞*˙) and *p*(*Ω*).

### **4.1.4 Friction model at mold surface**

As mentioned in Section 2.3.3, Barone and Caulk investigated Coulomb friction, constant friction and hydrodynamic friction as candidates for the thin lubrication layer between SMC core and the mold surface [124]. Hydrodynamic friction became the most common model and is generally described by

$$\tau\_{\mathsf{M}} = -\lambda \left( \frac{\left\| \boldsymbol{v}\_{\mathsf{s}} \right\|}{\nu\_0} \right)^{m-1} \boldsymbol{v}\_{\mathsf{s}},\tag{4.11}$$

where ⌧<sup>M</sup> is the stress at the mold surface, *∏* is a hydrodynamic friction coefficient, *m* is a power-law coefficient, *v*<sup>0</sup> is an arbitrary reference velocity for non-dimensionalization of the power-law term and *v*<sup>s</sup> is the slip velocity at the mold surface [148, 150]. The determination of friction parameters in literature focused mostly on a flow phase, where a stable plug-flow has formed and the initial squish flow has not be investigated in the same depth. During the squish phase, hydrodynamic friction models have the disadvantage to prescribe no friction stress at zero velocity. However, the initial contact may be rather described by a sticking behavior in some SMCs, which would favor a model with constant friction stress *ø*0. Such a model restricts the slip velocity to *v*<sup>s</sup> *=* 0 for wall shear stresses smaller than *ø*0.

In this work, a transition function ° is suggested to smoothly transition between a constant value *ø*0, which is applicable in absence of relative motion and a hydrodynamic friction model for relative motion. The transition function is defined as

$$\Gamma = e^{-\left(\frac{\|\|\mathbf{w\_k}\|\|}{\nu\_\mathbf{t}}\right)^2} \tag{4.12}$$

with a parameter *v*<sup>t</sup> that describes the transition velocity between sticking and slipping states. The maximum shear stress at the mold surface is then formulated as

$$\tau\_{\mathsf{M}} = -\Gamma \tau\_0 \left[ \upsilon\_{\mathsf{s}} \right] - (1 - \Gamma) \lambda \left( \frac{\left\| \upsilon\_{\mathsf{s}} \right\|}{\nu\_0} \right)^{m-1} \upsilon\_{\mathsf{s}}.\tag{4.13}$$

Relative velocity at mold surface k*v*sk

**Figure 4.1:** Comparison of friction models depending on relative velocity. While hydrodynamic models are widely used in literature [124,148,150], the initial squish phase may require a sticking constant friction at low slip velocities. Therefore, a transition model for sticking SMC is proposed in this work.

A graphic summary of constant friction (*v*<sup>t</sup> ! 1), linear hydrodynamic friction (*m =* 1, *v*<sup>t</sup> ! 0), power-law hydrodynamic friction (*m* 6*=* 1, *v*<sup>t</sup> ! 0) and a transition model is presented in Figure 4.1.

### **4.1.5 Virtual press controller**

The physical press follows a press profile given as *M* pairs of mold gap and corresponding closing velocity, as shown in Figure 4.2. Eventually the press controller switches to force-control in order to limit stresses on mold and press. A virtual press-controller is used to mimic this behavior as boundary condition during the compression molding simulation.

As long as the compression force is below the force at switch-over *F*max, the profile is linearly interpolated to obtain the current press velocity. Regions outside the prescribed profile are linearly extrapolated. The linear interpolation/extrapolation w.r.t the gap leads to a non-linear velocity profile w.r.t compression time.

**Figure 4.2:** The press profile is given as *M* pairs of mold gap and press speed and is linearly interpolated w.r.t the gap.

After the switch-over at time increment *ns*, a simple discrete PI-controller is employed to determine the current press speed

$$\dot{h}\_{n+1} = \dot{h}\_n + P\_{\rm p} e\_n + P\_{\rm l} \sum\_{q=n\_s}^{n} \frac{\Delta E\_q}{2} \Delta t \tag{4.14}$$

from the normalized errors

$$
\Delta E\_q = \frac{F\_{\text{max}} - F\_q}{F\_{\text{max}}} \dot{h}\_M,\tag{4.15}
$$

where *h*˙*<sup>M</sup>* describes the final compression speed of the press. The normalization ensures reliable force-control through a wide range of simulation parameters with constant control parameters *P*<sup>p</sup> *= P*<sup>i</sup> *=* 0.5.

### **4.1.6 Solution procedure**

The Coupled-Eulerian-Lagrangian framework of the commercial Finite Element solver Simulia Abaqus/Explicit is employed to solve the equations, as the deformations during compression molding are too large to be described

in a Lagrangian frame. This operator splitting method divides the equations in a Lagrangian part and an Eulerian part [192, 193]. The Lagrangian part

$$\left.\frac{\partial\rho}{\partial t}\right|\_{\mathcal{L}} = 0\tag{4.16}$$

$$\left.\frac{\partial(\rho v)}{\partial t}\right|\_{\mathcal{L}} = \text{div}(\sigma) \tag{4.17}$$

$$\left.\frac{\partial(\rho c\_{\rm p}T)}{\partial t}\right|\_{\rm L} = -\operatorname{div}(\mathbf{d}) + \boldsymbol{\sigma} : \mathbf{D} \tag{4.18}$$

$$\left. \frac{\partial \mathbf{A}}{\partial t} \right|\_{\mathbf{L}} = \mathbf{W} \cdot \mathbf{A} - \mathbf{A} \cdot \mathbf{W} + \lambda \left\{ \mathbf{D} \cdot \mathbf{A} + \mathbf{A} \cdot \mathbf{D} - 2 \mathbb{A} : \mathbf{D} \right\} \tag{4.19}$$

contains only the source terms and neglects the convective terms. These equations are solved in a traditional Lagrangian frame by explicit time integration and enable to account for contact with other Lagrangian bodies. A subsequent step is used to solve the convection problem

$$\left.\frac{\partial\rho}{\partial t}\right|\_{\mathcal{E}} + \text{div}(\rho v) = 0 \tag{4.20}$$

$$\left.\frac{\partial(\rho v)}{\partial t}\right|\_{\mathcal{E}} + \text{div}\left((\rho v)\otimes v\right) = 0\tag{4.21}$$

$$\left.\frac{\partial(\rho c\_{\rm p}T)}{\partial t}\right|\_{\rm E} + \text{div}\Big((\rho c\_{\rm p}T)\cdot\upsilon\Big) = 0\tag{4.22}$$

$$\left.\frac{\partial \mathbf{A}}{\partial t}\right|\_{\mathcal{E}} + \text{div}\,(\mathbf{A}\otimes\boldsymbol{\upsilon}) = \mathbf{0}.\tag{4.23}$$

This step effectively moves back nodes to remain an Eulerian frame and computes the flux through element boundaries by a second order advection scheme [192, 193]. An illustration of the process is given in Figure 4.3, where a circular shape is stretched horizontally. The initial shape is represented by a volume fraction of SMC *e* 2 © R *|* 0 *< e <* 1 ™ in each element and properties of partially filled elements are volume averaged.

The contact between Lagrangian bodies and the Eulerian material surface utilizes a material reconstruction algorithm [194]. The surface is approximated

**Figure 4.3:** Illustration of the CEL steps. The first step is equivalent to a conventional Lagrangian step, if the spatial time derivative would be replaced by a material time derivative. The second step moves the mesh back to its original configuration, computes the flux of material through element faces and adjusts the Lagrangian variables [193].

by a set of planes for each element that does not necessarily form a continuous surface. The reconstruction allows a definition of contact overclosures for a set of seed points on the contacting surfaces and thus allows to compute contact forces from a penalty approach. The resulting local contact force is then distributed on the elements nodes with their corresponding shape functions [193].

The constitutive model is implemented in a VUMAT subroutine, the virtual press controller is implemented as VUAMP subroutine and the friction is modeled with a pre-computed look-up-table depending on the absolute value of slip∞ ∞*v*<sup>s</sup> ∞ ∞.

## **4.2 Reduction to one-dimensional plug-flow in a press rheometer**

In some cases, such as the elongational plug-flow in a one-dimensional press rheometer, the previously described equations may be further simplified to a one-dimensional set of equations. It is desirable to be able to compute fast solutions to such a problem to verify more complicated three-dimensional implementations and to identify parameters with a press rheometer. Therefore, this section states a set of one-dimensional PDEs to approximate the compressible, non-Newtonian plug-flow in a press rheometer. Subsequently, an analytical solution for the further simplified incompressible Newtonian case is derived.

**Figure 4.4:** Dimensions of a press rheometer for one-dimensional elongational flow in the initial configuration (*t =* 0).

The dimensions of a generic press-rheometer are sketched in Figure 4.4. The solution domain is then described by a non-dimensional scalar position *x*§ *= x*/*X*, where the position *x* is normalized with the flow front position *X* and *x*§ 2 © R *|* 0 *< x*§ *<* 1 ™ . The fiber orientation is assumed to be planar, such that only the *A*xx, *A*xy, *A*yx and *A*yy components of *A* remain. Further, all entries of A except for *A*xxxx, *A*yyxx, *A*xyxx and corresponding symmetries vanish. The vertical strain rate is prescribed by the deformation speed *<sup>D</sup>*zz *<sup>=</sup> <sup>h</sup>*˙/*<sup>h</sup>* in terms of the mold gap *h*.

Isothermal bulk flow is a common assumption in literature (see Section 2.3). However, even small changes in temperature can have significant effect on the matrix viscosity. Hence, an analytical solution is employed to compute an

approximate average mold temperature during compression. The temperature between two closing plates with constant closing speeds and ideal heat transfer at the mold surfaces (Dirichlet boundaries) can be expressed as

$$T(t, \mathbf{z}) = T\_\Gamma + (T\_\Gamma - T\_0) \left[ \frac{2}{\pi} \sum\_{q=1}^\infty \frac{\cos(q\pi - 1)}{q} \frac{\sin(q\pi z)}{h} \exp\left(\frac{-q^2 \pi^2 \kappa t}{h\_0 h \rho c\_\mathbf{p}}\right) \right] \tag{4.24}$$

with a constant mold surface temperature *T*T, initial homogeneous stack temperature *T*<sup>0</sup> and *q* 2 N [132, 140]. Averaging the temperature distribution over the thickness

$$\bar{T}(t) = \frac{1}{h} \int\_0^h T(t, z) dz\tag{4.25}$$

yields

$$\bar{T}(t) = T\_{\Gamma} + (T\_{\Gamma} - T\_0) \left[ \frac{4}{\pi^2} \sum\_{q=1}^{\infty} \frac{\cos(q\pi - 1)}{q^2} \exp\left(\frac{-q^2 \pi^2 \kappa t}{h\_0 h \rho c\_{\rm p}}\right) \right] \tag{4.26}$$

as an approximation of the average temperature in the SMC stack.

### **4.2.1 Scalar set of PDEs**

For a model without short-range interactions, the constitutive relation may be formulated as

$$
\sigma = -p(\rho)I + \mathbb{V} : D' \tag{4.27}
$$

with an anisotropic viscosity tensor

$$\mathbb{V} = 2\eta(\bar{T}, \dot{\chi})\mathbb{P}\_2 + \eta\_{2, \text{SF}}(\bar{T}, \dot{\chi})\left(\mathbb{A} - \frac{1}{3}\mathbf{I} \otimes \mathbf{A}\right) \tag{4.28}$$

from equations (2.35), (2.39) and (2.42). With the assumptions stated above, this reduces to

$$V\_{\rm xxxx} = +\frac{4}{3}\eta(\bar{T}, \dot{\gamma}) + \eta\_{2, \rm SF}(\bar{T}, \dot{\gamma}) \left(A\_{\rm xxxx} - \frac{1}{3}A\_{\rm xx} \right) \tag{4.29}$$

$$V\_{\rm xzzz} = -\frac{2}{3}\eta(\bar{T}, \dot{\gamma}) + \eta\_{2, \rm SF}(\bar{T}, \dot{\gamma}) \left(\mathcal{A}\_{\rm X\overline{xzz}} - \frac{1}{3}\mathcal{A}\_{\rm Z}\right) \tag{4.30}$$

$$V\_{\rm Zxxx} = -\frac{2}{3}\eta(\bar{T}, \dot{\gamma}) + \eta\_{2, \rm SF}(\bar{T}, \dot{\gamma}) \left(A\_{\rm Z\widetilde{x}\widetilde{\rm X}} - \frac{1}{3}A\_{\rm xx}\right) \tag{4.31}$$

$$V\_{\rm ZZZZ} = +\frac{4}{3}\eta(\boldsymbol{T}, \dot{\boldsymbol{\gamma}}) + \eta\_{2, \rm SF}(\boldsymbol{T}, \dot{\boldsymbol{\gamma}}) \left(\mathcal{A}\_{\rm ZZZZ} - \frac{1}{3}\mathcal{A}\_{\rm EZ} \right) \tag{4.32}$$

for the relevant terms. The crossed out components vanish due to the assumption that all fibers are in the *x y*-plane.

Mass balance, momentum balance and the orientation equations are then formulated on the normalized domain as

$$
\dot{\rho} = -\rho D\_{\text{xx}} - \rho D\_{\text{ZZ}} \tag{4.33}
$$

$$
\rho \,\dot{\nu} = \frac{\partial}{\partial x^\*} \left( -p(\rho) + V\_{\text{XXX}} D\_{\text{XX}} + V\_{\text{XXZZ}} D\_{\text{ZZ}} \right) - 2 \frac{\pi}{h} \tag{4.34}
$$

$$\dot{A}\_{\text{xx}} = 2(A\_{\text{xx}} - A\_{\text{XXxx}})D\_{\text{xx}} \tag{4.35}$$

$$\dot{A}\_{\rm yy} = -2A\_{\rm yy}D\_{\rm xx} \tag{4.36}$$

$$
\dot{A}\_{\rm xy} = (A\_{\rm xy} - 2A\_{\rm xyxx})D\_{\rm xx},
\tag{4.37}
$$

where domain remains undeformed and progression of the flow front is account for by *@*(*•*) *@x*§ *<sup>=</sup>* <sup>1</sup> *X @*(*•*) *@<sup>x</sup>* , e.g. *<sup>D</sup>*xx *<sup>=</sup> @<sup>v</sup> @x*§ . The additional source term in the momentum equation describes a friction stress *ø* that acts on both mold surfaces (hence the factor 2) against the flow direction. The unknown fields are the mass density *Ω*, velocity *v* and the orientation components *A*xx, *A*xy and *A*yy.

### **4.2.2 Initial conditions**

The stack is considered to be initially at rest and at rest density *Ω*<sup>0</sup>

$$\left. \nu \right|\_{t=0} = 0, \quad \left. \rho \right|\_{t=0} = \rho\_0,\tag{4.38}$$

because it lays unconstrained on the lower mold half. The orientation is considered planar isotropic

$$\left.A\_{\rm{xx}}\right|\_{t=0} = 0.5, \quad \left.A\_{\rm{YY}}\right|\_{t=0} = 0.5, \quad \left.A\_{\rm{yx}}\right|\_{t=0} = 0.0,\tag{4.39}$$

if the SMC sheet production process (see Section 2.3) does not induce any preferential fiber orientation.

### **4.2.3 Boundary conditions**

The Dirichlet boundary condition at *x*§ *=* 0 is

$$\left.\nu\right|\_{\chi^\*=0} = 0\tag{4.40}$$

at all times. The boundary condition at the flow front *x*§ *=* 1 is given by the flux boundary

$$\left. -p(\rho) + V\_{\text{XXX}} D\_{\text{XX}} + V\_{\text{XXZZ}} D\_{\text{ZZ}} \right|\_{X^\* = 1} = 0 \tag{4.41}$$

during flow and by the Dirichlet boundary

$$\left.\nu\right|\_{\mathcal{X}^\*=1} = 0\tag{4.42}$$

as soon as the mold is completely filled. All other boundaries are zero-flux boundaries.

### **4.2.4 Solution procedure**

The set of coupled PDEs, initial conditions and boundary conditions is completed with algebraic expressions for utility variables. The first expression is the evolution of the flow front position

$$
\dot{X} = \nu|\_{X^\*} = 1 \tag{4.43}
$$

for *X < X*max. The second expression is the compression force

$$F\_{\mathbb{C}} = B \int\_{0}^{1} -p(\rho) + V\_{\text{ZZxx}} D\_{\text{XX}} + V\_{\text{ZZZZ}} D\_{\text{ZZ}} \text{d.x}^{\*},\tag{4.44}$$

which is integrated from the normal stress at the mold surface.

The set of coupled PDEs is solved using two *pdepe* solvers in MathWorks Matlab. The first solver integrates the equations until the mold is completely filled and a second solver restarts at that point with modified boundary conditions at *x*§ *=* 1. The algebraic expression for the flow front is solved simultaneously with a forward Euler time discretization and the compression force is evaluated using trapezoidal integration. The computation of the compression force is needed at run time, as it is used by a virtual press controller to switch from displacement control to force control, as soon as the maximum compression force is reached.

### **4.2.5 Analytical solution**

Further simplification of the one-dimensional flow allows to derive a simple analytical solution. Assuming incompressibility (*D*zz *=* °*D*xx *= D*) leads to a direct solution for the flow front position

$$X = \frac{h\_0}{h} X\_0 \tag{4.45}$$

and the velocity field

$$
\nu = -D\mathbf{x}.\tag{4.46}
$$

If the fiber orientation state is approximately constant, i.e. *A*xxxx *=* 0.375 and *A*xx *= A*yy *=* 0.5, the set of equations simplifies to just the momentum balance

$$\frac{\partial p}{\partial \mathbf{x}} = -2\rho D^2 \mathbf{x} + 2\frac{\pi}{h}. \tag{4.47}$$

Integrating this equation once for Newtonian viscosity and a constant friction *ø*<sup>0</sup> yields the pressure field

$$p = -\rho D^2(\mathbf{x}^2 - X^2) + 2\frac{\tau\_0}{h}(\mathbf{x} - X) - \underbrace{2\eta D - \eta\_{2,\text{SF}}(A\_{\text{xxx}} - \frac{1}{3}A\_{\text{xx}})D}\_{\text{from boundary condition}} \tag{4.48}$$

where the last term originates from the boundary condition (4.41). Finally, the compression force can be computed from Equation (4.44) as

$$F\_{\rm C} = B \left( \underbrace{\frac{2}{3} \rho D^2 X^3}\_{\text{inertia}} - \underbrace{\frac{\pi\_0}{h} X^2}\_{\text{friction}} - \underbrace{4 \eta D X - \eta\_{2, \text{SF}} D X A\_{\text{xxx}}}\_{\text{matrix}} \right) \tag{4.49}$$

with a contribution from inertia, friction, matrix viscosity and extra viscosity due to fiber-matrix long-range interactions. The deformation rate *D* is typically negative (compression) and thus each component adds a positive contribution to the overall compression force.

In case of a linear hydrodynamic friction *ø = ∏v*, the solution is

$$F\_{\rm c} = B \left( \frac{2}{3} \rho D^2 X^3 - \frac{\lambda}{3h} DX^3 - 4 \eta DX - \eta\_{2, \rm SF} DX A\_{\rm xxxx} \right), \tag{4.50}$$

which reduces to the result obtained by Castro [139, 140], if the inertial term and the fiber-matrix term are dropped.

## **4.3 Verification of one-dimensional plug-flow**

An exemplary solution of the 1D PDE model is given in Figure 4.5. In this generic example, the mold is closed at a rate of 1 mm s−<sup>1</sup> from an initial gap of 9 mm to a final gap of 2 mm. Due to compressibility of the material, the flow front advances in a non-linear fashion until the mold is completely filled after 8.5 s. The compression force reaches its maximum force at 6.1 s and the moldprofile switches to force control. The solution to the fields *Ω*, *A*xx, *A*yy, *A*xy is plotted as mean values over the mold length. Additionally, the pressure at several positions along the flow path is computed in a post-processing step to enable comparison with sensors in a press rheometry experiment. The exact sensor positions depend on the employed tool and are given for example in Section 7.1.4, where the one-dimensional model is used to determine friction parameters.

A simple example without a switch-over to force control is used to illustrate the general effect of friction, viscosity and compressibility. The resulting pressures are plotted in Figure 4.6 and parameters for the frictionless base model example are given in Table A.1. For these parameters, the pressure sensors are triggered one after another and jump to identical stress levels until the mold is filled after triggering the last sensor. The introduction of hydrodynamic friction (*∏*=50 N s m<sup>−</sup>3) leads to a pressure difference between sensors, as friction introduces a pressure gradient that increases the pressure level from the flow front towards the inward direction. An increase of the viscosity (*¥*=2000 Pa s) simply shifts the general pressure levels up. The times, at which sensors are triggered, remains unchanged in these three configurations, as the bulk modulus is chosen high enough to represent quasi-incompressible behavior. Lowering the bulk modulus (*K*0=1900 Pa) results in a compression of the stack and thus a delayed flow front advancement, which can be observed by the delayed time at which sensors are triggered.

Thus, the pressure sensors in the press rheometer deliver an unambiguous understanding of the SMC flow. The simple 1D model allows for a quick computation of the underlying effects and can be used to fit parameters to experimentally recorded results from the press rheometer.

**Figure 4.5:** Exemplary solution of the one-dimensional compressible flow in a press rheometer. The maximum compression force is reached after 6.1 s in this example before the filling completes at 8.5 s. The different colors in the pressure plot represent the pressures at eleven different sensor positions along the flow path.

**Figure 4.6:** Parameter variation in press rheometer model.

## **5 Direct Bundle Simulation**

## **5.1 Method development**

The fundamental idea of the proposed Direct Bundle Simulation (DBS) approach utilizes the observation that fiber bundles often remain in their bundled configuration during SMC compression molding. Hence, fiber bundles can be represented as individual instances in a direct Stokesian dynamics simulation with a reasonable number of fiber bundle instances. The number of bundles is roughly two or three orders of magnitude smaller than the number of fibers and thus computationally accessible.

### **5.1.1 Assumptions and scope**

Fiber bundles in SMC can be considered very flexible with a dimensionless bending stiffness *S*§ ø 1 due to the high aspect ratio and high viscosity of the suspending matrix (compare Equation (2.8)). The assembled structure of a bundle decreases its bending stiffness even further compared to a homogeneous bundle with the same cross section, because individual fibers may slip relative to each other during bending. Therefore, fiber bundles are modeled with truss elements, which are one-dimensional elements that only have a stretch and compression elasticity, but offer no resistance to bending load. Upon bending, individual truss elements remain straight and curvature is approximated by angles at the connecting nodes. As discussed previously in Chapter 4, the matrix is considered purely viscous and curing is neglected

during the flow process. The assumptions for the DBS model are summarized as


## **5.1.2 Generation of an initial bundle configuration**

First of all, the DBS requires an initial bundle configuration that resembles the configuration in the SMC sheets at the beginning of the compression molding process. The physical distribution results from bundles falling down on the lower foil in a random process (see Section 2.3), thus generating a random transversely isotropic in-plane orientation.

**Figure 5.1:** In a first step, a bundle direction *p*ˆ<sup>0</sup> *<sup>j</sup>* is drawn from a uniform isotropic distribution. It is then projected with a second order fiber orientation tensor and normalized to *p*<sup>0</sup> *j* . Finally, the bundle is translated randomly within the initial stack and protruding bundle elements are removed.

For the simulation, such a distribution is generated by first drawing random directions

$$\hat{\mathbf{p}}\_{\cdot j}^{0} = \sin(\hat{\theta})\cos(\hat{\phi})\mathbf{e}\_{\mathbf{x}} + \sin(\hat{\theta})\sin(\hat{\phi})\mathbf{e}\_{\mathbf{y}} + \cos(\hat{\theta})\mathbf{e}\_{\mathbf{z}},\tag{5.1}$$

where *¡*<sup>ˆ</sup> <sup>2</sup> *<sup>U</sup>*(0,2*º*) and *<sup>µ</sup>*<sup>ˆ</sup> <sup>2</sup> *<sup>U</sup>*(°1,1) with *<sup>U</sup>*(*a*,*b*) being a uniform probability distribution between *a* and *b*. This direction is then mapped to a prescribed initial fiber orientation state *A*0, which is often a planar isotropic state in SMC, by

$$\mathbf{p}\_{\rangle}^{0} = \left[\mathbf{A}^{0}\hat{\mathbf{p}}\_{\rangle}^{0}\right].\tag{5.2}$$

A fiber bundle comprised of several 1D-elements is generated with this direction and randomly translated within the rectangular domain of the initial stack by a shift vector

$$
\Delta \hat{x} = \Delta \hat{x} e\_{\mathbf{x}} + \Delta \hat{y} e\_{\mathbf{y}} + \Delta \hat{z} e\_{\mathbf{z}}, \tag{5.3}
$$

where ¢*x*ˆ 2 *U*(¢*x*min,¢*x*max), ¢*y*ˆ 2 *U*(¢*y*min,¢*y*max) and ¢*z*ˆ 2 *U*(¢*z*min,¢*z*max). The variables ¢*x*min and ¢*x*max are the minimum and maximum allowed shifts such that at least one node of the bundle resides within the initial stack domain. All elements outside the initial stack are deleted. This procedure is illustrated in Figure 5.1 and is repeated until the total fiber volume fraction equals the prescribed SMC fiber volume fraction.

The placement does not consider contacts between bundles. These initial contacts are resolved during an initial analysis step in the solver Simulia Abaqus/Explicit. Overclosures that cannot be resolved are stored as offsets and these offsets are resolved during the simulation.

### **5.1.3 Matrix model**

The matrix is described by similar equations as the macroscopic model

$$\frac{\partial \rho}{\partial t} + \text{div}(\rho v) = 0 \tag{5.4}$$

$$\frac{\partial(\rho v)}{\partial t} + \text{div}\left((\rho v) \otimes v\right) = \text{div}(\sigma) + \hat{f}^{\text{H}} \tag{5.5}$$

$$\frac{\partial(\rho c\_{\rm p}T)}{\partial t} + \text{div}\Big((\rho c\_{\rm p}T)\cdot v\Big) = -\text{div}(d) + \sigma : D.\tag{5.6}$$

However, the stress tensor is

$$\sigma = -p(\rho)I + 2\eta(T, \dot{\gamma})D'\tag{5.7}$$

with an equation of state for the entire composite material and the matrix shear viscosity *¥*(*T*,*∞*˙). The volumetric stress describes the entire SMC compression behavior including fiber bundles and voids, the deviatoric part describes only the matrix response without fiber bundles. A hydrodynamic body force field *f*ˆH has been added to model the influence of fibers on the matrix motion and is described in the subsequent section. The body force field on the matrix is denoted with a hat symbol (*•*ˆ) to distinguish it from the hydrodynamic force on bundles *f*H. The evolution equation for the fiber orientation state is dropped, because orientation tensors can be computed in post-processing from the simulated fiber bundle architecture.

### **5.1.4 Bundle-Matrix interaction<sup>1</sup>**

Bundles are represented as a chain of one-dimensional truss elements (=bundle segments) that do no transfer bending loads and the stiffness in elongational direction is determined by the fiber material. These truss elements do not share an interface with the Eulerian domain of the matrix, but float in the

<sup>1</sup> Parts of this section are based on [195].

same domain. The fiber-matrix interaction follows the concept of Stokesian dynamics (see Section 2.2.6.2) and computes the hydrodynamic force on each bundle segment from its velocity relative to the Eulerian environment. The resulting hydrodynamic force of each bundle segment *f*<sup>H</sup> is then distributed similarly to the approach by Lindström and Uesaka [85] as body force field *f*ˆH on the bulk matrix.

### **5.1.4.1 Parameterization of hydrodynamic resistance from micro model**

The hydrodynamic resistance of idealized cylindrical bundle segments is computed from a numerical study of micro-simulations. The fluid flow in a cube is simulated, while a cylinder is placed in its center. The fluid is subjected to a unidirectional flow with inlet boundary condition *v*<sup>1</sup> normal to the inlet face and a zero-pressure outlet condition, as shown in Figure 5.2. All other walls of the cube are allowed to slip, but the normal velocities are constrained to zero. A no-slip condition is applied at the cylinder surface.

The cylinder aspect ratio *r*<sup>p</sup> is varied in the range {1,2,3,5,8,13,25} and its orientation angle relative to the unidirectional flow *¡* is varied in the range {0,15±,30±,45±,60±,75±,90±}. The incompressible steady Stokes flow problem

$$0 = \operatorname{div} \left( -pI + \eta D' \right) \tag{5.8}$$

$$0 = \text{div}(\upsilon). \tag{5.9}$$

is solved for each configuration utilizing the commercial FEM solver Comsol Multiphysics.

To evaluate the resulting forces on the cylinder segment, its surface is parameterized as *r = r*r*e*<sup>r</sup> *+r*t*e*<sup>t</sup> *+r*p*p*, where {*e*r,*e*t,*p*} describes the local cylinder

**Figure 5.2:** Two-dimensional section through the cube, in which a Stokes flow is solved to compute the drag force on bundle segments. A constant velocity is prescribed at the inlet and a zero-pressure condition at the outlet. Other walls of the cube slip, while the cylinder interface is subjected to a no-slip condition. The cylinder is parameterized by its orientation angle *¡* and aspect ratio *r*p.

coordinate system. With this parameterization, the lateral cylinder surface is defined as

$$\mathcal{C} := \left\{ (r\_{\mathbf{r}}, r\_{\mathbf{t}}, r\_{\mathbf{p}}) \in \mathbb{R}^3 \mid r\_{\mathbf{r}} = R, 0 < r\_{\mathbf{t}} < 2\pi, 0 < r\_{\mathbf{p}} < l \right\},\tag{5.10}$$

where *R* is the cylinder radius and *l* is the length of a bundle segment. The surfaces at the cylinder ends are excluded, as these are typically connected to the next segment of the bundle and thus do not contribute to the overall resistance. The total hydrodynamic force exerted on the cylinder can be determined using an integral over the lateral surface *C* as

$$\mathbf{f}^{\mathrm{H}} = \int\_{\mathcal{C}} \boldsymbol{\sigma} \cdot \boldsymbol{n} \, \mathrm{d}A \tag{5.11}$$

with surface normal *n*.

A dimensionless measure for the hydrodynamic force in flow direction and the transverse direction is obtained with the non-dimensional forms

$$k\_{\rm d} = \frac{1}{6\pi\eta R\nu\_{\infty}} \int\_{\mathcal{C}} \boldsymbol{\sigma}\_{\rm x} \cdot \boldsymbol{n} \, \mathrm{d}A \quad \text{and} \quad k\_{\rm l} = \frac{1}{6\pi\eta R\nu\_{\infty}} \int\_{\mathcal{C}} \boldsymbol{\sigma}\_{\rm y} \cdot \boldsymbol{n} \, \mathrm{d}A \qquad (5.12)$$

from the vertical and horizontal surface stress components <sup>x</sup> and y. The non-dimensionalization is such that the drag factor *k*<sup>d</sup> becomes one for a suspended sphere and the lift factor *k*<sup>l</sup> becomes zero for a suspended sphere. Hence, the factors describe the resulting hydrodynamic force components in flow direction and perpendicular to the flow direction relative to a suspended sphere.

**Figure 5.3:** Dimensionless drag coefficient *k*<sup>d</sup> and lift coefficient *k*<sup>l</sup> from computation (dots) and fit according to Equation (5.13) and Equation (5.14) [195].

The resulting factors from the computation are plotted in Figure 5.3 as dots. For small aspect ratios, the factors approach those for a sphere independent of the orientation. The drag *k*<sup>d</sup> increases with larger aspect ratios and with an orientation closer to 90° due to an increased cross section normal to the flow direction. The lift component *k*<sup>d</sup> peaks at *¡ =* 45<sup>±</sup> and vanishes for a fiber aligned with the flow *¡ =* 0<sup>±</sup> or perpendicular to it with *¡ =* 90±. As for the drag, the severity of the peak increases with the aspect ratio. The computational results are approximated by a fit for the drag factor

$$k\_{\rm d}(r\_p, \phi) = 1 - a(r\_p - 1)\cos(2\phi) + \beta(r\_p - 1) \tag{5.13}$$

and the lift factor

$$k\mathbf{k}(r\_\mathbf{p}, \phi) = \alpha (r\_p - 1) \sin(2\phi) \tag{5.14}$$

with *Æ =* 0.09 and *Ø =* 0.3125. These fits allow efficient evaluation of the factors for a given orientation and aspect ratio. They are plotted as solid lines in Figure 5.3.

#### **5.1.4.2 Computation of the relative velocity**

The application of a Stokesian hydrodynamic force requires knowledge of the relative velocity between bundle segments and their environment. Figure 5.4 illustrates how a fiber segment of length *l* with index *j* 2 *B* at position *x<sup>j</sup>* and orientation *p<sup>j</sup>* floats through the Eulerian domain. The neighborhood of such a segment is defined as

$$\mathcal{N}\_{l} := \{ i \in \mathcal{E} \mid 0 < \left\| \Delta x\_{lj} \right\| < l \}, \tag{5.15}$$

where *E* describes the set of all indices of Eulerian elements that have an element volume fraction *e >* 0.5. The condition describes a sphere with a radius equal to the segment length around the segment center. The neighborhood relation is dynamic, as fiber bundles move independently from the

Eulerian mesh. It has to be updated every time step and the performance of the neighborhood search is critical to the overall performance of the DBS.

**Figure 5.4:** A fiber bundle with index *j* and orientation *pj* is positioned at *xj* . The neighborhood is called *N<sup>j</sup>* and one exemplary neighborhood element with index *i* is highlighted at position *xi* . The velocity at this element is *vi* and its relative position to the bundle segment is termed ¢*xi j* . [195]

A naive neighborhood search would include a check of condition (5.15) for each bundle and for each Eulerian element. This simple approach is displayed in Algorithm 1 and the time complexity of single run is *O*(*NBN<sup>E</sup>* ) with the number of bundle segments *<sup>N</sup><sup>B</sup>* and the number of filled Eulerian elements *N<sup>E</sup>* .

**Algorithm 1** Pseudo code for naive neighborhood search


For the number of elements in consideration here, this algorithm is too slow and some additional time to build more sophisticated data structures is negligible compared to potential gains during search. Such a more efficient data structure is a kd-tree, which is a space-partitioning binary search tree. For a 1d-tree, i.e. a tree that structures one-dimensional data, the data structure is built by sorting the data and taking the median as a node with a left and right node containing data that is arranged in the same way. A recursive binary bisection allows to find data points very effectively, as each choice for a next node rules out approximately half of the nodes at that stage [196]. A kd-tree is a generalization, in which a k-dimensional tree is build and searched, as indicated in Algorithm 2. The function *search* is evaluated recursively starting from a root node and decides which of its two child nodes is closer to the target position. If the node farther away is still within potential proximity of the target value, it is searched, too. Since the function *search* has *O*(log*N<sup>E</sup>* ) complexity due to the effective data structure, the overall performance is *O*(*N<sup>B</sup>* log(*N<sup>E</sup>* )), which is a massive improvement over the naive search method.

The tree is then used to determine *N<sup>j</sup>* for each segment during a time step and compute the relative velocity. The relative velocity for each bundle segment is determined as

$$
\Delta \boldsymbol{\sigma}\_{j} = \sum\_{i \in \mathcal{N}\_{j}} \frac{\boldsymbol{\mu}\_{ij}}{\boldsymbol{W}\_{j}} \left( \boldsymbol{v}\_{i} - \boldsymbol{v}\_{j} \right) \tag{5.16}
$$

**Algorithm 2** Pseudo code for kd-tree search

1: Build tree with nodes of *E* 2: **for** *j* 2 *B* **do** 3: *N<sup>j</sup>* = SEARCH(root, *x<sup>j</sup>* ) 4: **for** *i* 2 *N<sup>j</sup>* **do** 5: Compute relative velocity *v<sup>i</sup>* °*v<sup>j</sup>* for Equation (5.16). 6: **end for** 7: **end for** 8: 9: **function** SEARCH(node, *x*) 10: **if** node has no children **then** 11: Process terminal node 12: **return** results 13: **else** 14: Determine closer node and farther node 15: SEARCH(closer node, *x*) 16: **if** farther node is closer than *l* **then** 17: SEARCH(farther node, *x*) 18: **end if** 19: **end if** 20: **end function**

with

$$\|w\_{lj} = e^{-\frac{3}{2}\frac{\left\|\Delta x\_{lj}\right\|^2}{l^2}} \quad \text{and} \quad W\_{lj} = \sum\_{l \in \mathcal{N}\_f} w\_{lj}. \tag{5.17}$$

The weighting factors *wi j* and normalization with *Wj* represent a Gaussian weighting, where the variance is chosen to be *l*/3, such that the sum of weights outside *N<sup>j</sup>* represents less than 1%. The Gaussian weighting is a more accurate scheme than simply taking the nearest Eulerian element for the relative velocity and ensures that ¢*v<sup>j</sup> = v<sup>j</sup>* °*v* Ø Ø *xj* for *l* ! 0.

#### **5.1.4.3 Body force field evaluation**

The computed factors *k*d(*r*p,*¡*), *k*l(*r*p,*¡*) and the relative velocity ¢*v<sup>j</sup>* are used to compute the force on each bundle segment *f*<sup>H</sup> and the corresponding body force field *f*ˆH. The hydrodynamic force on each element is

$$\mathbf{f}\_{\dot{f}}^{\rm H} = 6\pi R\eta(T,\dot{\gamma}) \left[ k\_{\rm d}(r\_{\rm p},\phi\_{\dot{f}})\Delta v\_{\dot{f}} + k\_{\rm l}(r\_{\rm p},\phi\_{\dot{f}}) \left\| \Delta v\_{\dot{f}} \right\| q\_{\dot{f}} \right] \tag{5.18}$$

with the angle between the relative velocity and the fiber orientation

$$\phi\_f = \arccos\left(\frac{\Delta v\_f \cdot p\_f}{\left\|\Delta v\_f\right\|}\right) \tag{5.19}$$

and with the direction perpendicular to *v<sup>j</sup>* computed by

$$\mathbf{q}\_{\slash} = -\text{sgn}(\mathbf{p}\_{\slash} \cdot \Delta \mathbf{v}\_{\slash}) \left[ \mathbf{p}\_{\slash} - \left( \mathbf{p}\_{\slash} \cdot \left[ \Delta \mathbf{v}\_{\slash} \right] \right) \left[ \Delta \mathbf{v}\_{\slash} \right] \right] \tag{5.20}$$

in the plane that is formed by the two unit vectors *<sup>p</sup><sup>j</sup>* and <sup>á</sup> ¢*v<sup>j</sup>* à (details may be found in Appendix A.1).

Finally, the same Gaussian weights from the determination of the relative velocity are used to distribute the body force contribution of each bundle. Hence, the contribution of bundle *j* to the body force field in Eulerian element *i* is

$$\hat{\mathbf{f}}\_{ij}^{\text{H}} = -\frac{1}{V\_l} \frac{w\_{lj}}{W\_j} \mathbf{f}\_j^{\text{H}} \tag{5.21}$$

with the volume of the Eulerian element *Vi* . These contributions are summed up for each Eulerian element.

#### **5.1.4.4 Implementation**

The model is implemented with a combination of several subroutines in the commercial FEM program Simulia Abaqus/Explicit.

A VEXTERNALDB subroutine is utilized to control the entire analysis workflow. The subroutine parses the input deck before the start of the analysis to extract parameters such as index ranges of bundle elements and Eulerian elements. During the analysis it builds kd-trees at each time increment for all elements that have material volume fraction *e >* 0.5. The implementation for building and searching the kd-trees itself is taken from Kennel [196], who built one of the most efficient kd-trees in Fortran. A tree is build at the beginning of each time step using new memory such that the previous tree is still in memory for delayed parallel threads working on this tree. The memory is freed such that always two alternating trees, an old one and a new one, exist at all times (see Figure 5.5). Additionally, the VEXTERNALDB subroutine computes an approximate stable time increment for the bundles

$$
\Delta t\_{\rm B} = 0.25 \sqrt{\frac{R}{2a\_{\rm max}}} \tag{5.22}
$$

based on the idea that a bundle should travel no more than a fraction of a bundle radius during a step with the maximum acceleration of all bundles *a*max.

A VUFIELD subroutine is used to extract positions and velocities at nodes at each time step. The values are saved to field variables to make them accessible to a VUSDFLD subroutine, which interpolates position and velocity at integration points of Eulerian elements as well as truss elements. Each Eulerian element (EC3D8R or EC3D8RT) and each truss element (T3D2) has only one unique integration point. The variables *e*,*T*,*V*,*¥*,*x*,*v*,*f*H/*f*ˆH at the integration points are saved to global arrays with a length equal to the sum of truss and Eulerian elements for further processing in other subroutines. The global arrays are accessed with the index of elements and due to the

**Figure 5.5:** The first tree is built using allocated memory for the maximum tree size at the beginning of the first time increment. This tree should not be overwritten at the beginning of the next time step, as some threads may still operate on this data structure. Hence, the second tree is built in different pre-allocated memory. Then, the third tree overwrites the memory of the first tree to save time on allocation/deallocation and the second tree remains in the second allocated memory space. This procedure with alternating trees omits any waiting for threads, which is potentially disastrous for threaded parallel performance.

parallel domain decomposition, the access happens in parallel while still being thread safe.

Finally, a VDLOAD subroutine evaluates equations (5.16), (5.18) and (5.21). It performs parallel tree searches in the current kd-tree, which do not alter the tree, but only obtain the neighborhood relation.

## **5.1.5 Bundle-Bundle interaction<sup>2</sup>**

The objective of this section is a model for short-range interaction of bundlebundle contact points. The model needs to determine the mechanical contact in normal direction and tangential friction as well as lubrication in normal an tangential direction.

A simple mechanical contact can be modeled quite easily as repulsive force

$$f^{\rm n,e}\_{\alpha\beta} = \begin{cases} 0 & \text{g} > 0 \\ K\_{\rm c} h\_{\alpha\beta} & \text{g} \le 0 \end{cases} \tag{5.23}$$

<sup>2</sup> Parts of this section are based on [32].

with a contact stiffness *K*<sup>c</sup> and the gap between bundle surfaces *g* . The force is absent, if fibers are separated by a gap, and penalizes overclosure linearly. The corresponding tangential friction force is

$$\mathcal{f}^{\mathrm{t,e}}\_{\alpha\beta} = \begin{cases} 0 & \mathbf{g} > 0 \\ -\mu \left\| \mathcal{f}^{\mathrm{n,e}}\_{\alpha\beta} \right\| \left\| \Delta v\_{\alpha\beta} \right\| & \mathbf{g} \le \mathbf{0} \end{cases} \tag{5.24}$$

with a Coulomb friction coefficient *µ*. The friction acts opposed to the direction of the relative tangential velocity ¢*vÆØ*.

The Newtonian form of the tangential lubrication given in (2.73) is

$$\mathbf{f}\_{\alpha\beta}^{\text{t,l}} = \begin{cases} \eta \frac{d\_{\text{a}}^2}{\left\| \underbrace{\mathbf{p}\_{\alpha} \times \mathbf{p}\_{\beta}}\_{\text{contact area}} \right\|} & \underbrace{\Delta \mathbf{v}\_{\alpha\beta}}\_{\text{eff}} & \mathbf{g} > \mathbf{0} \\\\ \mathbf{0} & \mathbf{g} \le \mathbf{0} \end{cases} \tag{5.25}$$

and other models [30, 56, 57] may be expressed with it using a proper choice of the effective sheared gap *G*. This work uses the effective sheared gap as a single adjustable parameter, because it can be interpreted as the mapping from a physical distance between bundles *g* to an equivalent sheared effective gap *G* accounting for all geometrical details of the contact area. This relation is obtained in the following two sections using analytical considerations and a parametric simulation study.

#### **5.1.5.1 Analytical considerations for the effective sheared gap**

The effective sheared gap *G* is equal to the separation distance *g* in case of two flat fiber bundles with parallel faces. In this case, the interaction force becomes infinite, if the separation distance between two such bundles approaches zero. This section aims at finding an expression of *G* for bundles

with elliptical cross sections and investigates the behavior when the separation distance becomes close to zero.

**Figure 5.6:** Schematic illustration of the contact area between two fiber bundles. The first bundle *Æ* is aligned with the *x*-direction and forms the angle *√* between this bundle and a second bundle *Ø*. [32].

Two identical fiber bundles with constant elliptical cross sections are considered. The bundle cross section is described by a minor axis diameter *d*<sup>b</sup> and a major axis diameter *d*a. Both bundles are positioned with their minor axes aligned in *z*-direction, as depicted in Figure 5.6. Bundle *Æ* is aligned with the *x*-direction and a second bundle *Ø* crosses it such that an angle *√* is formed between them in the *x y*-plane. Bundle *Ø* is placed such that a separation distance *g* remains at the closest contact point between both bundles. The coordinates (*x*, *y*, *z*), the separation distance *g* and the effective sheared gap *G* are non-dimensionalized with the bundle dimensions as

$$\mathbf{x}^\* = \frac{\mathbf{x} \| \sin \psi \|}{d\_\mathbf{a}}, \quad \mathbf{y}^\* = \frac{\mathbf{y}}{d\_\mathbf{a}}, \quad \mathbf{z}^\* = \frac{\mathbf{z}}{d\_\mathbf{b}}, \quad \mathbf{g}^\* = \frac{\mathbf{g}}{d\_\mathbf{b}}, \quad G^\* = \frac{G}{d\_\mathbf{b}} \tag{5.26}$$

to obtain generalized results. The dimensionless *z*§-component of bundle surface *Æ* may be expressed as

$$z\_{\alpha}^{\*} = -\frac{1}{2} + \sqrt{{\+y^{\*}(1 - y^{\*})}}\tag{5.27}$$

with *y*§ 2 [0,1]. The dimensionless *z*§-component of bundle surface *Ø* is consequently

$$z\_{\beta}^{\*} = \frac{1}{2} - \sqrt{\frac{1}{4} - \left(\mathbf{x}^{\*} - \mathbf{y}^{\*}\cos\psi\right)^{2}} + \mathbf{g}^{\*} \tag{5.28}$$

with *x*§ 2 h °1 <sup>2</sup> *<sup>+</sup> <sup>y</sup>*§ cos(*√*), <sup>1</sup> <sup>2</sup> *+ y*§ cos(*√*) i . The upper and lower bound depend on the *y*-position as this bundle is placed with an angle *√* relative to bundle *Æ*. The separation distance is added to the parameterization of this bundle surface to place the point closest to the surface of bundle *Æ* at distance *g* relative to bundle *Æ*.

A strong simplification is introduced by assuming that the velocity profile between both surfaces is approximately linear. Hence, the velocity gradient at each infinitesimal surface element of the interaction zone is given as *<sup>∞</sup>*˙ *<sup>=</sup>* ¢*vÆØ <sup>z</sup>Ø*°*z<sup>Æ</sup>* with the velocity difference of the bundles ¢*vÆØ*. Assuming Newtonian viscosity *¥*, an integration over the sheared domain yields

$$\mathcal{J}\_{\alpha\beta}^{\mathrm{l,t}} = \eta \Delta v\_{\alpha\beta} \frac{d\_{\mathrm{b}} d\_{\mathrm{a}}^2}{|\sin(\psi)|} \underbrace{\int\_0^1 \int\_{-\frac{1}{2} + y^\* \cos(\psi)}^{\frac{1}{2} + y^\* \cos(\psi)} \frac{1}{z\_{\beta}^\* - z\_{\alpha}^\*} \mathrm{d}x^\* \mathrm{d}y^\*}\_{1/G^\*}. \tag{5.29}$$

Inserting equations (5.27) and (5.28) and using the substitution *s*§ *= x*§ ° *<sup>y</sup>*§ cos*√<sup>+</sup>* <sup>1</sup> <sup>2</sup> , this can be simplified to

$$\frac{1}{G^\*} = \int\_0^1 \int\_0^1 \frac{1}{1 - \sqrt{\mathbf{s}^\*(1 - \mathbf{s}^\*)} - \sqrt{\mathbf{y}^\*(1 - \mathbf{y}^\*)} + \mathbf{g}^\*} \mathbf{ds}^\* \mathbf{d} \mathbf{y}^\*,\tag{5.30}$$

where it becomes apparent that the result is independent of *√*.

Solving the inner integral leads to

$$\frac{1}{G^\*} = \int\_0^1 -\pi + \frac{4k\_1(\mathbf{g}^\*, \mathbf{y}^\*)}{\sqrt{4k\_1(\mathbf{g}^\*, \mathbf{y}^\*)^2 - 1}} \left(\frac{\pi}{2} + \text{arccot}\sqrt{4k\_1(\mathbf{g}^\*, \mathbf{y}^\*)^2 - 1}\right) \text{d}\mathbf{y}^\* \tag{5.31}$$

with *k*<sup>1</sup> *=* 1*+g* §° <sup>p</sup>*y*§(1° *<sup>y</sup>*§). After a change of integration variables, with another substitution *k*<sup>2</sup> *=* q <sup>4</sup>*<sup>g</sup>* §<sup>2</sup> °8*<sup>g</sup>* §*u*§ *<sup>+</sup>*4*u*§<sup>2</sup> *<sup>+</sup>*8*<sup>g</sup>* § °8*u*§ *<sup>+</sup>*<sup>3</sup> and the firstorder approximation arctan*<sup>x</sup>* <sup>º</sup> *<sup>º</sup>* <sup>4</sup> *x* for small values of *g* §, this is further simplified to

$$\frac{1}{G^\*} = -4 \int\_0^{\frac{1}{2}} \frac{\mu \left[ (2 + \mathbf{g}^\* - \boldsymbol{\mu}^\*) \pi k\_2(\mathbf{g}^\*, \boldsymbol{\mu}^\*) - (4 \mathbf{g}^\* - 4 \boldsymbol{\mu}^\* + 4) \pi \right]}{k\_2(\mathbf{g}^\*, \boldsymbol{\mu}^\*) \sqrt{-4 \boldsymbol{\mu}^{\*2} + 1}} \mathrm{d}\boldsymbol{\mu}^\*. \tag{5.32}$$

Unfortunately, the result of this integration can only be expressed in terms of incomplete elliptic integrals. However, a series expansion leads to an analytical expression that is exact if *g* § approaches zero. Employing only the first series term, the effective sheared gap for two bundles can be expressed as

$$G^\* = \frac{8}{\pi (-16 \arctan \frac{\sqrt{3}}{3} + 16\sqrt{3} - 48 + 24 \ln 2 - 8 \ln \text{g}^\* + \pi)} + O(\text{g}^\*). \tag{5.33}$$

The dimensionless effective sheared gap depends only on the dimensionless separation distance *g* § and is independent of the contact angle *√* and the elliptical cross section. The effective sheared gap approaches zero, if the separation distance approaches zero. This means that interaction forces would become infinite for infinitesimal close contact, just as in the case of flat bundles. However, the assumption of linear velocity profiles is strong and should be compared to full solutions of Stokes' equations on the sheared domain. Therefore, the next section presents a parametric study of the sheared gap that solves the creeping flow between two bundles numerically.

### **5.1.5.2 Parametric study of the effective sheared gap with a micro model**

The sheared fluid between two fiber bundles is simulated to determine the effective sheared gap *t* between bundles. A parametric model of the sheared domain is built using the following parameters for a parametric sweep:


A small buffer zone around the domain is added to the sheared domain to prevent poor element quality at the sharp boundaries. An illustration of the parametric model for one exemplary configuration is given in Figure 5.7.

Stoke's equations for creeping flow are solved on the domain highlighted green in Figure 5.7. The momentum balance is

$$\mathbf{0} = \operatorname{div}\left(-p\boldsymbol{I} + \eta \left(\operatorname{grad}(\boldsymbol{v}) + \operatorname{grad}(\boldsymbol{v})^{\top}\right)\right) \tag{5.34}$$

**Figure 5.7:** Two fiber bundles with a sheared domain (green) between them. The elliptical cross section is characterized by minor and major diameter *d*<sup>b</sup> and *d*a, respectively. The separation distance between bundles is denoted as *g* and the upper bundle is placed with an angle *√* towards the lower bundle. The load direction *m* (=direction of relative motion) is characterized by an angle *±* ranging from −90° to 90° due to symmetry. [32]

with pressure *p*, Newtonian viscosity *¥* and velocity *v*. The fluid is assumed incompressible, hence the mass balance gives

$$\text{div}(\upsilon) = 0.\tag{5.35}$$

The boundary conditions are illustrated in Figure 5.8. The wall in contact with bundle *Æ* is fixed with zero velocity. The wall in contact with bundle *Ø* is subjected to a tangential velocity with magnitude *v*<sup>0</sup> and direction *±* within the shear plane. The remaining walls are subjected to a zero-pressure boundary condition and the inflow or outflow of fluid is limited to a direction *m* that corresponds to the load direction.

Figure 5.9 visualizes some computed effective sheared gaps *G*§ for different aspect ratios of the cross section in each sub-figure. Each set of lines indicated by the same plot style is composed from simulation results for contact angles *√* 2 [15±,30±,45±,60±,75±,90±]. Conversely to the analytical solution, these results vary with load direction and small peaks occur if the load direction matches the orientation angle.

**Figure 5.8:** Boundary conditions for a sheared volume with *r*<sup>c</sup> *=* 2 and *¡ =* 45±. The blue areas indicate the surfaces at which boundary conditions are applied. The vectors *v* denote the velocity boundary condition at the corresponding surfaces. The vector *m* denotes the constrained direction of velocity at the Neumann boundaries and depends on the load angle *±*. [32]

Figure 5.10 shows the mean value of the dimensionless effective sheared gap *G*§ plotted against the dimensionless contact gap *g* § for all contact angles and load directions. The colored lines represent the simulation results for different aspect ratios of the cross section, which correspond to the colors in Figure 5.9. The error bars indicate the standard deviation within each set due to varied contact angles and load angles. The analytical result of Equation (5.33) is plotted as a solid black line.

The mean values of the simulation results fall on one curve in Figure 5.10 regardless of the cross section aspect ratio *r*c, even though the models vary significantly in size and shape. The simplified analytical solution seems to overestimate the effective sheared gap. This may be attributed to the pressure distribution and resulting non-linear velocity profiles across the gap, while the velocity profile is assumed linear in the analytical solution.

Literature values for the effective gap *G*§ are given for example by Le Corre et al. [8] and Guiraud et al. [104]. Le Corre et al. estimate a value *G =* 2.0 £ <sup>10</sup>°<sup>3</sup> mm from compression experiments for bundles with minor diameter *db =* 0.06 mm, which is equivalent to *G*§ *=* 0.033. Guiraud et al. report *<sup>G</sup> <sup>=</sup>* 1.5 £ <sup>10</sup>°<sup>2</sup> mm from pull-out experiments and *db <sup>=</sup>* 0.2mm, which is equivalent to *G*§ *=* 0.075. Both literature results are added to Figure 5.10

**Figure 5.9:** Dependence of *G*§ on load angle for exemplary simulation results with varying contact angles *√* (line style) and varying aspect ratios (color) at different separation distances *g*§. Not all simulations are included to enhance visibility. Mean values of all load angles and contact angles are used to generate the diagram in Figure 5.10. [32]

at *g* § *=* 0 for reference. The magnitude of these reported values is in agreement with the results of the analytical considerations and parametric study presented here.

The analytical considerations showed that the interaction force between two elliptical bundles becomes infinite for infinitesimal small separation distance. However, this singularity cannot be introduced to numerical simulations. Therefore, an offset for the effective sheared gap *G*§ <sup>0</sup> º 0.04 is proposed and

**Figure 5.10:** Comparison of several approaches to determine the dimensionless effective sheared gap *G*§ as a function of the dimensionless separation distance *g*§. [32]

used in a simple fit of the results shown in Figure 5.10. One possible fit for the simulation results is given as

$$\mathbf{G}^\* = \mathbf{0}.027 \arctan(65.5 \mathbf{g}^\*) + \mathbf{g}^\* + \mathbf{G}\_0^\*. \tag{5.36}$$

This approach fulfills the condition

$$\left.G^\*\right|\_{\mathbb{g}^\*}\right|\_{\mathbb{g}^\*} = G\_0^\* \tag{5.37}$$

to achieve a finite separation for bundle contact and the condition

$$\mathbf{g}^\*|\_{\mathbf{g}^\* \to \infty} = \mathbf{g}^\* \tag{5.38}$$

for large separations where the surface curvature becomes irrelevant compared to the separation.

**Figure 5.11:** The variables passed into the subroutine (here *g*ˆ) assume a cylindrical shape (dashed lines), which overlap in this case. However, bundles of thickness *d*<sup>b</sup> may not overlap and a gap remains.

#### **5.1.5.3 Implementation**

Bundle-bundle interactions are implemented with a VUINTERACTION subroutine in Simulia Abaqus/Explicit. By default, the contact variables of any truss or beam contact imply a cylindrical shape, but one may modify these parameters in the subroutine to account for a different shape. In this case, the gap is defined as

$$\mathbf{g} = d + \hat{\mathbf{g}} - d\_{\mathbf{b}} \tag{5.39}$$

where *<sup>d</sup> <sup>=</sup>* <sup>p</sup>*d*a*d*<sup>b</sup> is the equivalent cylinder diameter, *<sup>g</sup>*<sup>ˆ</sup> is the gap between cylinders as provided by Simulia Abaqus/Explicit' contact solver and *d*<sup>b</sup> is the bundle thickness. This definition is a simple approximation to correct the contact thickness from a cylindrical surface to an elliptical shape, as illustrated in Figure 5.11. Strictly, this correction applies only for bundles that are perfectly aligned with their minor axis and neglects contact of bundle edges along the major axis. However, a contact of bundle edges is rare in SMC compared to the contact of the relatively flat bundle surfaces. Further, the rheology is likely dominated by the large sheared zones between bundles, which share a contact area in the plane of their major diameters.

In this implementation, the normal contact stress is defined as

$$\mathbf{s}\_{\mathrm{n}} = \begin{cases} -K\_{\mathrm{c}} \mathbf{g} + K\_{\mathrm{d}} \Delta \boldsymbol{\nu}\_{\mathrm{n}} & \mathbf{g} \le \mathbf{0} \\ \eta \frac{d\_{\mathrm{d}}^2}{|\sin \boldsymbol{\psi}|} \frac{\Delta \boldsymbol{\nu}\_{\mathrm{n}}}{\overline{\boldsymbol{G}} A\_{\mathrm{c}}} & \mathbf{g} > \mathbf{0}, \end{cases} \tag{5.40}$$

where *K*<sup>c</sup> is the contact stiffness (typically Young's modulus), *K*<sup>d</sup> is a contact damping parameter (typically *¥*), *A*<sup>c</sup> is the contact area and ¢*v*<sup>n</sup> is the relative velocity of the contact partners in the normal direction. The tangential contact stresses are defined as

$$\mathbf{s}\_{\mathbf{t}} = \begin{cases} -\mu \left\| \begin{smallmatrix} \mathbf{f}^{\mathrm{n},\mathrm{e}}\_{\alpha\beta} \end{smallmatrix} \right\| \left\| \Delta \mathbf{v}\_{\mathbf{t}} \right\| & \mathbf{g} \le \mathbf{0} \\ -\eta \frac{d\_{\mathrm{a}}^{2}}{|\sin\psi|} \frac{\Delta v\_{\mathbf{t}}}{GA\_{\mathrm{c}}} & \mathbf{g} > \mathbf{0} \end{smallmatrix}, \tag{5.41}$$

where ¢*v*<sup>t</sup> is the relative velocity of the contact partners in tangential direction. The dependence of *G* on the contact gap *g* is modeled with Equation (5.36).

### **5.1.6 Post-processing of bundles**

The DBS results in a bundle configuration represented by a set of nodes and line elements for each simulation time step. Post-processing is required to derive macroscopic properties of the result and compare it to macroscopic models and experimental results.

#### **5.1.6.1 Bundle curvature**

The bundle curvature is a measure for buckling of fiber bundles and is typically high in areas that experience compression during the process. The curvature is approximated as

$$c\_k = \frac{2}{l} \tan\left(\frac{1}{2}\arccos\left(p\_{\cdot j} \cdot p\_{\cdot \prime}\right)\right) \tag{5.42}$$

103

at each node *k* connecting two neighboring bundle segments of length *l* with unit directions *p<sup>j</sup>* and *p<sup>j</sup>*<sup>0</sup> . This evaluation assigns each nodal position an approximate curvature value, which can than be compared to curvature results from computer tomography.

#### **5.1.6.2 Contact number**

To compute the total number of contacts between bundles, a distance map between bundle element centers is computed. This distance map describes for each bundle element center the distance to each other bundle element center. Potential contacts occur, where the distance between two centers *j* and *j*0 is smaller than the segment length *l*. For these possible contact points, the normal to both directions is

$$m = \left[p\_{j'} \times p\_{j'}\right] \tag{5.43}$$

and the relative position can be determined by solving the linear system of equations

$$\begin{aligned} \mathbf{x}\_{J} - \mathbf{x}\_{J'} &= \begin{bmatrix} p\_{J,\mathbf{x}} & n\_{\mathbf{x}} & -p\_{J',\mathbf{x}} \\ p\_{J,\mathbf{y}} & n\_{\mathbf{y}} & -p\_{J',\mathbf{y}} \\ p\_{J,\mathbf{z}} & n\_{\mathbf{z}} & -p\_{J',\mathbf{z}} \end{bmatrix} \begin{bmatrix} u\_{1} \\ u\_{2} \\ u\_{3} \end{bmatrix} \end{aligned} \tag{5.44}$$

for *<sup>u</sup>*<sup>1</sup> to *<sup>u</sup>*3. A contact occurs, if<sup>Ø</sup> Ø*u*<sup>1</sup> Ø <sup>Ø</sup> *<sup>&</sup>lt; <sup>l</sup>*/2 and<sup>Ø</sup> Ø*u*<sup>3</sup> Ø Ø *< l*/2 (that is, the intersection is within the range of both bundle lengths) as well as Ø Ø*u*<sup>2</sup> Ø Ø *< d* (that is, the intersection is close enough to be considered a contact of overlapping cylinders). The evaluation of contacts during the simulation allows to estimate the influence of bundle-bundle interactions on the resulting compression force.

### **5.1.6.3 Mapping of discrete orientation tensors and fiber volume fraction**

If the fiber architecture has to be compared to a macroscopic model or if the results of the process simulation are required as input for further structural simulation, a mapping of line elements to an arbitrary mesh is necessary. Hence, for each cell *i*, the length of each potential bundle segment *j* in this cell ¢*li j* is computed. If bundles cut a face of the cell, the length within the cell up to the intersection point is computed, as indicated in Figure 5.12.

**(a)** Both nodes within cell **(b)** Both nodes outside cell **(c)** Partial intersection **Figure 5.12:** Case distinction for computation of the length ¢*li j* in a cell.

The local fiber volume fraction is then approximated as

$$f\_l = \pi \frac{d\_\mathbf{e}^2}{4} \sum\_{j \in \mathcal{B}\_l} \frac{\Delta l\_{lj}}{V\_l},\tag{5.45}$$

which assumes that the remaining cell volume is filled completely by the composite. The local discrete second order fiber orientation tensor is

$$\mathcal{A}\_{l} = \frac{1}{\sum\_{j \in \mathcal{B}\_{l}} \Delta I\_{lj}} \sum\_{j \in \mathcal{B}\_{l}} \Delta l\_{lj} \mathbf{p}\_{j} \otimes \mathbf{p}\_{j},\tag{5.46}$$

where *B<sup>i</sup>* Ω *B* is the set of bundles with a length ¢*li j >* 0 in cell *i*. The procedure is implemented as a ParaView filter and allows generic mapping between CAE files in the VTK file format. <sup>3</sup>

## **5.2 Verification**

## **5.2.1 Orientation evolution of a single bundle<sup>4</sup>**

The motion of a single bundle in shear flow must be faithfully reproduced to qualify the method for computation of fiber reorientation. Therefore, a fiber bundle with a length of 25 mm and an aspect ratio of *r*<sup>p</sup> *=* 25 is subjected to a shear rate *<sup>∞</sup>*˙ *<sup>=</sup>* 10s°<sup>1</sup> in this verification. The domain for this simulation is illustrated in Figure 5.13. The bundle is placed at the center, discretized with ten segments and positioned vertically, such that the initial orientation is *µ =* 0.

**Figure 5.13:** The contour plot shows a fiber bundle discretized with ten segments in a shear flow. The color codes indicate the velocity in x-direction. The fluctuations at both fiber bundle ends show how the two-way coupling influences the macroscopic velocity field [195].

<sup>3</sup> The filter is available at

<sup>4</sup> Parts of this section are based on [195].

Figure 5.13 shows bundle position and velocity in x-direction shortly after starting the simulation. The contour plot of the horizontal velocity component depicted in Figure 5.13 indicates the two-way coupled nature of the presented approach, particularly visible at the two ends of the fiber bundle. Although the bundle is flexible, it behaves like a rigid body until alignment with the flow due to the positive normal stress in the direction of the bundle axis.

Figure 5.14 compares the orientation evolution of the DBS with ten truss elements and two truss elements to the solution of Equation (2.22) with aspect ratio an equivalent aspect ratio computed from the formula by Cox [13]. The simulation is in good agreement with the reference solution for both discretizations.

**Figure 5.14:** Comparison of bundle orientation angle computed from DBS and Jeffery's Equation with an bundle aspect ratio *r*<sup>p</sup> *=* 25 [195].

There is a small difference between simulation and analytical solution at the almost horizontal state in Figure 5.14. At this point, torque induced by friction at the lateral surface dominates bundle motion. In SMC, bundles are heavily confined by other bundles and the mold. It is assumed that the torque that spins a free bundle in a dilute situation is small compared to the confinement effects and therefore is neglected here.

Additionally, a bundle with a bundle aspect ratio *r*<sup>p</sup> *=* 25 is placed 90° to the flow under the same conditions as in the parameter identification (see Figure 5.2 in Section 5.1.4.1) and meshed with one and ten segments. The resulting drag force normalized with 6*º¥Rv*<sup>1</sup> is 9.31 and 9.55, respectively. This is close to each other, but slightly smaller than the drag coefficient shown in Figure 5.3, because the averaged velocity around the bundle is smaller than the nominal velocity far away. Nevertheless, the orientation result and the total drag indicate that bundle motion is generally only slightly affected by discretization. The choice of bundle discretization may be considered the choice of which length scale of velocity disturbance is computed by the Eulerian mesh and which length scale is included in the drag coefficients.

### **5.2.2 Anisotropic flow front progression**

The presence of fiber bundles leads to anisotropic flow and the extreme case of unidirectional orientation is considered in this verification. Therefore, an initial rectangular stack of the dimensions 25mm£25mm£15mm and fiber bundle length 10mm is compressed with 1mms°<sup>1</sup> between two plates and simulation results of the macroscopic reference model, the DBS and an isotropic macroscopic model are qualitatively compared. Additional simulation parameters can be found in Table A.2.

Tables 5.1 and 5.2 show the flow front progression of the suspension for 0.25% fiber volume fraction and 2.5% fiber volume fraction, respectively. The left columns show isotropic results, where the stack is evenly stretched in horizontal and vertical directions. There is a slight tendency towards a more circular shape, as the edges of the initial stack are rounded during compression.

The anisotropic models consider unidirectional horizontal fiber bundle alignment. For 0.25% fiber volume fraction, there is a slight preference for flow perpendicular to fiber bundles, as the elongation in fiber bundle direction is prohibited by quasi-inextensible fiber bundles. The DBS shows rounder

**Table 5.1:** Semi-dilute anisotropic flow with unidirectional fiber orientation and 0.25% fiber volume fraction.

corners than the mascroscopic model, which is probably due to the inhomogeneous dispersion of bundles in the domain (there are just 125 bundles in the model) and specifically due to the absence of long bundles in the corners, as they are cut outside the stack. In contrast, the macroscopic model assumes equal fiber aspect ratios at every material point without considering the stacks boundaries. The flow front shape is similar for both models and the anisotropy is described qualitatively correct. For a higher fiber volume

fraction of 2.5% the anisotropy becomes more severe and there is much more elongation perpendicular to fiber bundles than in fiber bundle direction. The round corners of the DBS model likely arise from a lack of long fibers in the corners. The overall comparison shows that DBS can describe the anisotropic motion of the fluid qualitatively correct.

### **5.2.3 Frictionless linear compression**

To verify quantitatively correct compression forces due to anisotropic viscosity, a quasi-incompressible press rheometer configuration (see Figure 4.4) without friction is considered in this section. Therefore, an initial rectangular stack of the dimensions *X*<sup>0</sup> *=* 50mm,*B =* 50mm,*h*<sup>0</sup> *=* 5mm is compressed with a constant compression speed *<sup>h</sup>*˙0 *<sup>=</sup>* 1.66mms°<sup>1</sup> between two rigid plates. The mold length is *X*max *=* 75mm, which results in an initial mold coverage of 66% and a volumetric fill time of 1 s. Snapshots of the mold filling process of such a configuration with 25% fiber volume fraction using DBS are shown in Figure 5.15.

**Figure 5.15:** Snapshots of the compression molding process with 25% fiber bundles. The top mold is not displayed for visualization purposes.

First, the resulting compression force of the pure matrix without any fibers is computed by the analytical Equation (4.49), the 1D PDE model given in Section 4.2, the 3D VUMAT model 4.1 and a DBS without any bundles. The results are depicted in green colors in Figure 5.16 and show congruent compression forces. The compression force increases non-linearly, as the magnitude of the compression strain rate *<sup>D</sup> <sup>=</sup> <sup>h</sup>*˙0/*<sup>h</sup>* increases non-linearly, since *<sup>h</sup>* is reduced with a constant rate. The results at the very beginning are left out of this visualization for all numerical models, because the material is initially at rest and has to accelerate after establishing a contact with the mold.

To verify the viscosity increase due to added bundles, 300 bundles with a nominal length of *L =* 25mm are added employing the procedure described in Section 5.1.2. After the removal of protruding bundle elements, this equates to approximately 1% fiber volume fraction in the stack. The DBS model is then solved with disabled bundle-bundle interactions and with enabled bundle-bundle interactions assuming a bundle width *d*<sup>b</sup> *=* 0.5mm. As indicated by orange bullets in Figure 5.16, the introduction of 1% fiber bundles approximately triples the reaction force compared to pure matrix. The incorporation of short-range bundle-bundle interactions hardly changes the result, as interaction between bundles is relatively rare. On average, the 1610 bundle segments experience 131 contact pairs resulting in only 0.16 contacts per bundle segment, or in other words, every sixth fiber bundle is in potential contact with another bundle at all.

For reference, the analytical solution, 1D model and 3D VUMAT model are solved for 1% fiber volume fraction as well. These models require an estimation of a fiber aspect ratio that is comparable to the generated architecture with a length distribution due to the cutting process. The average fiber bundle length is computed as

$$\bar{L} = \frac{1}{X\_0 B\_0} \int\_0^{B\_0} \int\_0^{X\_0} L(\mathbf{x}, \mathbf{y}) \mathrm{d}x d\mathbf{y} \tag{5.47}$$

**Figure 5.16:** Comparison of analytical, 1D macro, 3D macro and DBS models to verify the implementations and models. The abbreviation *M*, *BM*, *BB* stand for a matrix component, twoway-coupled bundle-matrix interaction and bundle-bundle interactions, respectively.

and results to *<sup>L</sup>*¯ *<sup>=</sup>* 7/12*<sup>L</sup>* in this case, because the initial stack is small in comparison to the bundle length. The aspect ratio definition is not unambiguous for fiber bundles and candidates are the ratio between some average

fiber length and fiber diameter, equivalent bundle diameter, major bundle diameter. For simplicity, the aspect ratio between average fiber length and equivalent bundle diameter is chosen here, which is *r*<sup>p</sup> *=* 72.9. The analytical solution depicted as solid orange line in Figure 5.16 uses the constant values *A*xx *=* 0.5, *A*xxxx *=* 0.375 for the fiber orientation state, which represents an ideal planar isotropic distribution. The solution of the 1D model is depicted as dotted orange line and includes the computation of fiber orientation evolution. It starts at the same compression force level as the analytical equation, but the force increase exceeds it, as fibers align with the flow. The 3D VU-MAT solution (orange dashed) matches the 1D solution precisely for disabled short-range bundle-bundle interactions. With bundle-bundle interactions (*kD =* 5, solid light orange line), the compression force is increased by a rather small amount, which is in line with the expected behavior at this low fiber volume fraction and with the increase introduced by bundle-bundle interactions in DBS.

At 25% fiber volume fraction, the compression force level of DBS without bundle-bundle interactions is significantly increased compared to the previous results due to the increased in-plane elongational viscosity. Again, the 1D solution starts at the same compression force level as the analytical solution, but deviates with increasing time due to fiber alignment (violet solid line and dotted line). The 3D VUMAT solution shows a slightly higher compression force response, as this simulation requires usage of a shear number *Ns =* 10 to stabilize the numerical solution procedure. Very high ratios of elongational viscosity to shear viscosity tend to cause transverse shear bands that destabilize the solution process.

In the DBS model with 25% fiber volume fraction, 87120 contacts pairs form between the 37577 involved bundle segments on average, which means that each segment is in contact with 4.63 other bundles. Therefore, a significant contribution of short-range bundle-bundle interactions is expected. Hence, bundle-bundle interactions are introduced with 0.25 mm, 0.5 mm and 1.5 mm bundle width in the DBS model and indicated with light violet stars, squares and triangles in Figure 5.16, respectively. The cross section area

of bundles is not modified, such that wider bundles are proportionally flatter. It becomes apparent, that at such a high fiber volume fraction, short-range hydrodynamic bundle-bundle interactions indeed lead to a significant increase of the compression force. The compression force increases with flatter and wider bundles. For a width of 0.5 mm, the hydrodynamic interaction factor *k*<sup>D</sup> *=* 5 yields similar results to the DBS.

**Figure 5.17:** Evolution of discrete orientation tensor components during the linear compression example. The abbreviations BM and BB stand for bundle-matrix interaction and bundle-bundle interaction, respectively.

While the compression force can be altered significantly by the presence of bundle-bundle interactions, the orientation evolution experiences only little change. This is demonstrated in Figure 5.17 for the DBS with 25% fiber volume fraction, where one-way coupling refers to an application of Equation (5.18), but not Equation (5.21), which results in the same compression force as the pure matrix behavior. The globally evaluated discrete second order fiber orientation tensor components are hardly influenced by two-way coupling and bundle-bundle interactions.

This section verified, that the two-way coupled bundle-matrix interaction increases the compression force significantly. Without bundle interactions, the analytical solution, 1D model, 3D VUMAT model and DBS predict similar

orders of compression force magnitudes. The effect of bundle-bundle interaction further increases the compression force. Bundle-bundle interaction introduces locally large interaction forces between bundles, which affect the time step of the explicit DBS negatively. While bundle-bundle interactions have a large effect on compression force, they have a weak effect on fiber re-orientation, as many bundle-bundle interactions likely cancel each other w.r.t the overall motion of a bundle. If the fiber architecture is of primary interest in a simulation, bundle-bundle interactions could be omitted for the benefit of better computational times.

However, the choice of *r*<sup>p</sup> and *k*<sup>D</sup> is ambiguous. The choice of *r*<sup>p</sup> *=* 72.9 is reasonable, but it cannot be considered the definitive result, because the reference model assumes rods with constant length, while it is applied here for bundles with a length distribution. Additionally, the application of Shaqfeh and Fredrickson's model at 25% is certainly beyond the semi-concentrated regime and the aspect ratio is not uniform for all fiber bundles. A validation with experimental results of a press rheometer is necessary to validate the results and evaluate the importance of contributions from two-way bundlematrix interaction, bundle-bundle interaction and friction at the mold surface. This will be addressed in Section 7.1.4 and Section 8.3.

## **6 Patch model**

## **6.1 Method development**

The patches are local continuous carbon fiber reinforcements that are preformed to a desired shape and subsequently co-molded with SMC. They can be placed in critical areas of a part to provide a tailored improvement of mechanical performance while keeping the amount of costly unidirectional carbon fibers at a minimum. This chapter addresses the formulation of a hyperviscoelastic model including a simple damage model.

### **6.1.1 Assumptions and scope**

Typical shell formulations of UD fabric forming assume decoupled membrane and bending behavior of the fabric, to account for the relative slip between fibers during bending [171]. However, the investigated material is much stiffer than fabrics in wet compression molding or thermoforming at process conditions and does not show a clear separation between membrane and bending properties. This assumption is verified in the characterization experiments in Chapter 7. Given the high stiffness and the initial temperature of the material in the process, the preform's ability to deform further during the molding process is limited. Patches rather disintegrate upon large deformations instead of being draped into a complex shape during co-molding. Consequently, the main objective of the co-molding simulation in this context is a prediction of the damage initialization in patches under certain process

conditions. Hence, a simple Hashin criterion is implemented to describe patch damage by matrix tension during co-molding.

The patch material considered here is made from a stitch-bonded unidirectional carbon fiber fabric. It is stitched with glass fiber yarn as sketched in Figure 6.1. As a simplification of the stitching pattern, it is assumed that the stitching can be represented by a single direction that is initially transverse to the carbon fiber.

**(a)** Front **(b)** Back **Figure 6.1:** Sketch of the stitching pattern of the fabric used for fabrication of the unidirectional reinforcement patches (ZOLTEK PX35 UD300).

Considering the observed material behavior, the assumptions are summarized as follows:


### **6.1.2 Hyperviscoelastic model**

The material response (membrane, bending, shear) is modeled based on an *Ideal Fiber Reinforced Material* (IFRM). The IFRM model assumes inextensible fibers and incompressibility, which form a reaction term. The stress can be expressed as

$$\sigma = \underbrace{-pI + S\_{\text{F}} \left(p\_{\text{F}} \otimes p\_{\text{F}}\right)}\_{\text{reaction term}} + \tau \tag{6.1}$$

with the hydrostatic pressure *p*, a fiber stress *S*<sup>F</sup> in fiber direction *p*<sup>F</sup> and an extra-stress ⌧ , which is modeled via a constitutive equation [197, 198].

Following Dörr [199], a hyperviscoelastic model based on the IFRM is formulated in the initial configuration as

$$\mathbf{S} = \mathbf{S}\_{\mathbf{e}} + \mathbf{S}\_{\mathbf{V}} + \mathbf{S}\_{\mathbf{r}} \tag{6.2}$$

with *S* denoting the second Piola-Kirchoff stress tensor. It comprises an elastic (e), viscous (v) and reaction (r) part. The isotropic elastic part is a hyperelastic St. Venant-Kirchoff material with

$$\mathbf{S}\_{\mathbf{e}} = \left(\frac{\mathbf{\mathring{E}\_{\mathbf{M}}} \mathbf{\mathring{G}\_{\mathbf{M}}}}{3\mathbf{\mathring{G}\_{\mathbf{M}} - \mathbf{\mathring{E}\_{\mathbf{M}}}} \mathbb{P}\_1 + 2\mathbf{\mathring{G}\_{\mathbf{M}} \mathbb{P}\_2} \right) \colon \mathbf{E},\tag{6.3}$$

where *E*ˇM and *G*ˇM are the isotropic matrix Young's modulus and shear modulus, respectively. The notation *•*ˇ denotes that these properties depend on the matrix damage variable. The viscous part is

$$\boldsymbol{S}\_{\text{V}} = \underbrace{\boldsymbol{F}^{-1} \cdot \left(2\check{\boldsymbol{\eta}}\_{\text{M}} \underbrace{\boldsymbol{\mathcal{F}}^{-\top} \cdot \dot{\boldsymbol{E}} \cdot \boldsymbol{F}^{-1}}\_{\text{current configuration}}\right) \cdot \boldsymbol{F}^{-\top}}\_{\text{current configuration}},\tag{6.4}$$


where the rate of deformation *E*˙ is transformed to the current frame, then multiplied by twice the patch viscosity *¥*ˇ *<sup>M</sup>* and then transformed back to the initial configuration. The reaction term is expressed as the regularization

$$\mathcal{S}\_{\mathbf{f}} = \check{E}\_{\mathbf{F}} \varepsilon\_{\mathbf{F}} \left( p\_{\mathbf{F}} \otimes p\_{\mathbf{F}} \right) \quad \text{with} \quad \varepsilon\_{\mathbf{F}} = E : \left( p\_{\mathbf{F}} \otimes p\_{\mathbf{F}} \right) \tag{6.5}$$

representing the extension in fiber direction. Finally, the second Piola-Kirchoff stress is transformed to the current configuration by Equation (2.83).

### **6.1.3 A simple damage model**

Large deformations of the unidirectional patches are often accompanied by damage in this application. Therefore, the prediction of damage initialization during the molding process is critical to ensure that patches are designed such that they stay intact.

Hashin developed a model specifically for unidirectional fiber reinforced composites and distinguished the four failure modes matrix tension, matrix compression, fiber tension and fiber compression [200, 201]. Previous observations of UD patches in co-molding show that matrix tension and tension of stitching yarn are the dominating failure mechanisms. Tension failure of the UD carbon fibers is hardly observed under processing conditions and compression results more likely in instabilities than actual failure. Hence, this section only considers damage initialization and progression under tension loading. Patches are considered shell-like structures, allowing to apply plane stress criteria.

#### **6.1.3.1 Damage initialization**

Hashin suggests a smooth piece-wise failure criterion for unidirectional composites [200, 201]. The initialization of matrix damage in plane stress is

$$
\left(\frac{\sigma\_{22}}{Y\_{22}}\right)^2 + \left(\frac{\sigma\_{12}}{Y\_{12}}\right)^2 = 1 \quad , \quad \sigma\_{22} > 0 \tag{6.6}
$$

where *æ*<sup>22</sup> denotes the elastic in-plane transverse Cauchy normal stress and *æ*<sup>12</sup> denotes the elastic in-plane Cauchy shear stress. The tensile strength and shear strength of the matrix are given by *Y*<sup>22</sup> and *Y*12, respectively. Fiber tension damage is initialized at

$$\left(\frac{\sigma\_{11}}{Y\_{11}}\right)^2 = 1 \quad , \quad \sigma\_{11} > 0,\tag{6.7}$$

where *Y*<sup>11</sup> is the tension strength of the fiber and *æ*<sup>11</sup> is the tensile elastic Cauchy stress in fiber direction. Carbon fiber damage is neglected due to their superior strength. However, the stitching yarn represents also a layer of unidirectional fiber reinforcement and its failure is relevant to the characteristic failure for UD patches during co-molding.

#### **6.1.3.2 Damage progression**

Strain-softening damage behavior can lead to mesh-dependent energy release in Finite Elements. To alleviate this drawback, the following damage progression is formulated in terms of displacement instead of strain utilizing the characteristic length *l*<sup>c</sup> of an element [202]. The equivalent displacement and the initial equivalent stress at onset of matrix tension damage are

$$\delta\_{\rm M0} = l\_{\rm c} \sqrt{\langle E\_{22} \rangle^2 + 2E\_{12}^2} \quad \text{and} \quad \sigma\_{\rm M0} = \frac{l\_{\rm c}}{\delta\_{\rm M0}} \left\{ \langle E\_{22} \rangle \langle \sigma\_{22} \rangle + 2E\_{12} \sigma\_{12} \right\}, \tag{6.8}$$

respectively. Here, *l*<sup>c</sup> denotes a characteristic element length and h*•*i *=* 1/2(*•+ |•|*) is the Macaulay bracket, which is used to consider tension stress only. The initial equivalent displacement and the corresponding initial equivalent elastic stress at onset of fiber damage are

$$\mathcal{S}\_{\rm F0} = l\_{\rm c} \sqrt{(E\_{11})^2} \quad \text{and} \quad \sigma\_{\rm F0} = \frac{l\_{\rm c}}{\delta\_{\rm F0}} \left( E\_{11} \sigma\_{11} \right), \tag{6.9}$$

respectively. After damage initialization, the current matrix damage variable *D*<sup>M</sup> is defined as

$$D\_{\rm M} = \frac{\delta\_{\rm Mf}(\delta\_{\rm Mt} - \delta\_{\rm M0})}{\delta\_{\rm Mt}(\delta\_{\rm Mf} - \delta\_{\rm M0})} \tag{6.10}$$

with

$$\delta\_{\rm Mt} = l\_{\rm c} \sqrt{\langle E\_{22} \rangle^2 + 2E\_{12}^2} \quad \text{and} \quad \delta\_{\rm Mf} = 2 \frac{W\_{\rm M}}{\sigma\_{\rm M0}} \tag{6.11}$$

representing the current equivalent displacement *±*Mt and the final equivalent displacement *±*Mf at which the material is completely damaged. The final equivalent displacement is computed from the tensile fracture energy of the matrix *W*<sup>M</sup> (per area due to the displacement formulation) that is dissipated during the fracture. The matrix Young's modulus *E*M, shear modulus *G*<sup>M</sup> and viscosity *¥*<sup>M</sup> are degraded by

$$
\check{E}\_{\rm M} = (1 - D\_{\rm M}) E\_{\rm M} \quad , \quad \check{G}\_{\rm M} = (1 - D\_{\rm M}) G\_{\rm M} \quad \text{and} \quad \check{\eta}\_{\rm M} = (1 - D\_{\rm M}) \eta\_{\rm M}. \tag{6.12}
$$

Analogously, the evolution of fiber tension damage is described by

$$D\_{\rm F} = \frac{\delta\_{\rm Ff}(\delta\_{\rm Mt} - \delta\_{\rm F0})}{\delta\_{\rm Mt}(\delta\_{\rm Ff} - \delta\_{\rm F0})} \tag{6.13}$$

with

$$\delta\_{\rm Ft} = l\_{\rm c} \sqrt{\langle E\_{22} \rangle^2} \quad \text{and} \quad \delta\_{\rm Ff} = 2 \frac{W\_{\rm F}}{\sigma\_{\rm F0}}. \tag{6.14}$$

The fiber stiffness *E*<sup>F</sup> is degraded by

$$
\check{E}\_{\rm F} = (1 - D\_{\rm F}) E\_{\rm F}.\tag{6.15}
$$

This results in the behavior illustrated in Figure 6.2. The equivalent stress increases linearly until it reaches *æ•*0, where *•* is either F for fiber or M for matrix. Further loading results in damage and linear softening, where the final point *±•*<sup>f</sup> is defined by the total fracture energy, i.e. the total area under the solid curve. Unloading at partial damage results in the dashed line, where the stress reduces with a flatter slope. The damage variable can never be reduced, such that a subsequent loading would result in the identical slope.

**Figure 6.2:** Relation between equivalent stress and equivalent displacement with linear softening.

### **6.1.4 Implementation**

The model described above is implemented in a Simulia Abaqus/Explicit subroutine VUMAT. A single ply of patch is represented by two element layers, one representing unidirectional carbon fibers and the matrix, and a second one representing the stitching yarn in transverse direction. This model allows for damage in the matrix due to shear, while the transverse stitching yarn remains intact and can be damaged separately due to transverse tension. Both layers share all their nodes, as illustrated in Figure 6.3.

If the damage variable approaches 1, elements may experience significant deformation due to the stiffness degradation. Hence, elements are deleted, if *D*<sup>M</sup> *=* 1 or *D*<sup>F</sup> *=* 1.

**Figure 6.3:** A single ply of patch is represented by two layers of elements. The fiber orientation *p*F of the UD layer is initially perpendicular to the orientation of stitching yarn.

## **6.2 Verification**

The implementation of the patch model is verified for a single element with generic parameters to ensure that the model behaves as expected. A single layer with a single element 1 mm x 1 mm x 0.1 mm is loaded transverse to fiber direction with different fracture energies. The element is loaded beyond its strength and then unloaded back to its initial configuration.

For a transverse stiffness *E*<sup>M</sup> *=* 10GPa and tensile matrix strength *Y*<sup>22</sup> *=* 500MPa, results are shown in Figure 6.4. The matrix fracture energy *W*<sup>M</sup> *=* 25mJmm°<sup>2</sup> is twice the nominal elastic energy stored after reaching *Y*22. Therefore the slope after damage initiation should be identical apart from the sign. If the fracture energy is reduced to half its value, the damage should occur instantaneous after reaching the strength. Increasing the fracture energy should result in a flatter slope after reaching the strength limit. The depicted behavior in Figure 6.4 agrees with these expectations. Equivalently, correct behavior of tension in fiber direction and shear is verified.

**Figure 6.4:** Single element verification of damage implementation for different matrix fracture energies. In all cases, the stress increases with the slope prescribed by the Young's modulus up to the prescribed tensile matrix stress. For *<sup>W</sup>*<sup>M</sup> *<sup>=</sup>* 12.5mJmm°2, the energy (area under the curve) is just equivalent to the fracture energy at this point and the fracture occurs suddenly. For *<sup>W</sup>*<sup>M</sup> *<sup>=</sup>* 25mJmm°2, further displacement reduces the force with the same slope, as it previously increased. For *<sup>W</sup>*<sup>M</sup> *<sup>=</sup>* 50mJmm°2, the corresponding slope is flatter, since the energy dissipation until complete failure is larger. Unloading results in the expected reduced stiffness for all cases.

## **7 Characterization**

## **7.1 SMC**

The primary material under investigation is a glass fiber reinforced SMC based on an unsaturated polyester-polyurethane hybrid resin (UPPH) with 25 mm long glass fibers. The material was developed to improve co-molding with unidirectional carbon fiber patches [170]. The employed glass fiber multi-end roving has a total Tex of 4800 and consists of 80 strands. Hence, cutting one multi-end roving once yields 80 fiber bundles with 200 fibers with 14 µm diameter each, that fall on the SMC carrier foil. The complete composition of the material is given in Table 7.1. The material is produced at Fraunhofer ICT, Pfinztal, Germany.


**Table 7.1:** Composition of UPPH glass fiber Sheet Molding Compound

It has to be noted that the material under investigation is not subject to a rigorous quality control. Hence, the local composition and thickness varies within a single SMC bobbin and especially between different batches. Additionally, the properties depend on manufacturing conditions, storage time, storage conditions etc. Therefore, the goal of this characterization is to find a reasonable approximate description of the material behavior, but not *the* definitive parameterization of the material.

## **7.1.1 Thermal properties**

The balance equations for the inner energy of the macroscopic model (4.3) and the DBS model (5.6) are closed with the transverse heat conductivity *∑* and specific heat capacity *c*<sup>p</sup> of SMC. Additionally, the solution requires a thermal gap conductance *k*<sup>T</sup> to model the heat flow at the contact between SMC and mold surfaces. These parameters are identified in this section for uncured UPPH-GF SMC material.

**Method.** The thermal properties are determined from temperature measurements in a stack of SMC sheets and were conducted at Fraunhofer ICT, Pfinztal, Germany. Ten sheets of 50 mm x 50 mm x 1.1 mm are stacked and thermocouples are placed centrally between each layer. The first sensor *T*<sup>1</sup> is placed on the surface of the stack, such that it is positioned between mold and stack as soon as the stack comes in contact with the mold. The stack is then surrounded with glass wool insulation and placed on a mold surface heated to 145 °C. A small weight on top of this configuration ensures proper contact. The described setup is illustrated in Figure 7.1.

**Results.** The measured temperatures *T* of three samples are shown in Figure 7.2 as light colored lines for sensors T1 to T8 for the first 60 s after contact initialization. The transient heat transfer is approximated as a 1D process according to

$$
\rho c\_p \frac{\partial T}{\partial t} = \kappa \frac{\partial^2 T}{\partial z^2} \tag{7.1}
$$

**(a)** Sketch of the setup **(b)** Photo of the prepared stack **Figure 7.1:** Setup for the evaluation of the transverse heat conductivity. Ten sheets of SMC are stacked with temperature sensors between each layer and placed on a heated plate, while being insulated on all other boundaries.

with *z* 2 [0,*D*s] and the stack height *D*<sup>s</sup> *=* 11mm. The specific heat capacity *c*<sup>p</sup> is approximately 1530 J kg−<sup>1</sup> K<sup>−</sup>1, which is estimated from data of the resin system [203] and glass by a rule of mixture. The boundary conditions are

$$\kappa \left. \frac{\partial T}{\partial z} \right|\_{t,z=0} = -k\_\Gamma (T\_\mathcal{M} - T) \quad \text{and} \quad \kappa \left. \frac{\partial T}{\partial z} \right|\_{t,z=D\_\delta} = 0 \tag{7.2}$$

with a constant mold temperature *T*<sup>M</sup> *=* 145±C. Initially, the temperature is homogeneous at *T |t=*0,*<sup>z</sup> =* 24±C.

An optimal fit to the measured data is obtained by solving the initial value problem for an initial guess of *∑* and *k*T, computing a scalar squared error and minimizing the error iteratively. The solution of Equation (7.1) with optimal parameters (summarized in Table 7.2) is plotted at the sensor positions in Figure 7.2.

**Discussion.** The results for senors *T*<sup>9</sup> and *T*<sup>10</sup> are not considered, as they do not show any relevant change during the time frame investigated here. The experimental data deviates from the numerical solution of (7.1) due to manufacturing related inhomogenities and uneven sheet thickness. The fitted parameters are considered an averaged idealization of the heat transfer process. The specific heat capacity just affects the heat flux into the SMC, but

**Figure 7.2:** Measured temperatures (light colored lines) and best fit based on one-dimensional heat transfer equation (dark colored lines) for sensors T1 to T8


**Table 7.2:** Thermal properties of UPPH-GF SMC in B-staged state.

not the resulting temperature. Hence, an error in the approximation from the rule of mixture does not affect the resulting temperature profiles in this fit or the non-isothermal simulations presented later, if the same specific heat capacity is used in both cases. The heat flux might be computed inaccurately, if the specific heat capacity is inaccurate, but this is of no interest in this work.

### **7.1.2 Paste viscosity**

The viscosity of the pure paste *¥*(*T*,*∞*˙) with additives, but without fibers, is of major interest to this work. It is used in the macroscopic reference model from Chapter 4 (see Table 4.1) and it is needed in the DBS from Chapter 5. For example, the DBS requires the paste viscosity for the matrix contribution

in Equation (5.7), to compute drag forces in Equation (5.18) and to compute lubrication forces between bundles in Equation (5.25).

**Method.** Specimens for viscosity measurement of the SMC paste are prepared by filling the paste in a mold instead of processing it on the SMC line after mixing. The mold is filled up to approximately 1 mm thickness and sealed with a styrene tight foil. The paste is matured for two weeks at room temperature similar to the SMC for molding. Round coupons of 25 mm diameter are cut from the cured paste and placed in a Anton Paar MCR 501 rheometer at Fraunhofer ICT, Pfinztal, Germany in plate-plate configuration. The viscosity is measured in oscillatory mode at 20 °C, 40 °C and 80 °C for three specimens each.

**Results.** The measured viscosity shows typical power-law behavior, as depicted in Figure 7.3. The displayed viscosity refers to the norm of the measured complex viscosity.

**Figure 7.3:** Measured viscosity (points) and best fit.

As the power-law behavior has a singularity at zero shear rate, a Cross-WLFlike model

$$\eta(T,\dot{\gamma}) = \frac{\eta\_0(T)}{1 + \left(\frac{\dot{\gamma}}{\dot{\gamma}\_0}\right)^{1-n}} \quad \text{with} \quad \eta\_0(T) = D\_1 e^{\frac{-a\_1(T-T^\*)}{a\_2 + (T-T^\*)}}\tag{7.3}$$

is used to model the paste viscosity. A fixed transition shear rate *∞*˙0 *=* 0.1 is employed to limit the viscosity at small shear rates. The power-law coefficient *n*, transition shear rate *∞*<sup>0</sup> and parameters *T* §, *D*1, *Æ*1, *Æ*<sup>2</sup> are fitted to the experimental results. The fit is obtained by minimizing the squared sum of the normalized least-squares error at each measured temperature for all specimens. The results of the best fit are illustrated in Figure 7.3 and the optimal parameters are summarized in Table 7.3.


**Table 7.3:** Viscous properties of UPPH paste in B-staged state.

**Discussion.** All computed values of the fit displayed in Figure 7.3 are within the measured data range, except for low shear rates at 80 °C. The parameters should not be considered strictly Cross-WLF parameters, as the experimental data in Figure 7.3 shows no plateau. The model is mainly chosen to smoothly transition to a finite viscosity as the shear rate approaches zero.

### **7.1.3 Compaction**

In both models, the macroscopic reference model from Chapter 4 and the DBS from Chapter 5, pressure is related to the volumetric strain via an equation of state *p*(¢*"*). This relation is obtained in this section as tabulated data for subsequent simulations.

**Method.** SMC stacks of 800 mm x 450 mm size, which represents 100% mold coverage, are consolidated in a Dieffenbacher DYL630/500 hydraulic press with parallelism control at Fraunhofer ICT, Pfinztal, Germany. The mold temperature is 145 °C and the mold closes with a constant speed of 1 mm s<sup>−</sup>1. Automated cutting of SMC sheets ensures an accurate initial mold coverage and compression data (force, mold displacements) is recorded during the trials. The expected relation between volumetric strain and pressure is illustrated in the schematic in Figure 7.4: Upon first contact, the bulk material offers small resistance to compression, as trapped air is compressed and released. The resistance increases, as an increasing amount of pores is closed and air escapes until the maximum compression force is reached. At this constant compression force, the material first expands due to heating and consequently shrinks due to cross linking of the thermoset. Finally, the cured part expands elastically as the compression force is relaxed and the mold opens. The final recorded gap is the part thickness.

**Results.** The measurement in Figure 7.5 shows the behavior sketched in Figure 7.4 for three samples molded with the identical configuration of three SMC sheets each and is plotted with light colors. The thickness of the demolded part is indicated by dots and is in agreement with the measured displacement at release. However, the entire mold and press are elastic and an empty stroke is used to subtract the pressure-dependent mold displacement. The so corrected measurements are plotted in dark colors in Figure 7.5. Only the region of the rising flank is of interest for the compression molding simulation here. Therefore, these points are extracted and plotted over the Hencky strain ¢*" =* log (*h*/*h*0) in Figure 7.5. An averaged relation between strain and pressure is obtained by averaging the strains at various pressure

**Figure 7.4:** Schematic relation between pressure and mold gap for compression. 1) First contact between mold and SMC stack. 2) The maximum compression force is reached. 3) The maximum thermal extension is reached. 4) The part is fully cured. 5) The part is demolded.

levels. This tabulated data (see Table 7.4) is then used to interpolate the equation of state *p*(¢*"*) in subsequent simulations. In rare cases, where the tabulated data range is exceeded, the data is extrapolated linearly.


**Table 7.4:** Relation between Hencky strain and pressure for UPPH-GF SMC in B-staged state.

**(a)** Measurement **(b)** Interpolation **Figure 7.5:** Compressive behavior of the SMC. (a) The recorded data (light color) is corrected with the stiffness of the mold and press, which is determined from an empty stroke. The corrected data is shown in dark colors and the final part thickness agrees with the recorded mold gap measurements. (b) The rising flank is extracted and an interpolation is used as tabulated data to describe the relation between Hencky strain and pressure.

**Discussion.** The resin system under investigation shows compressible behavior, as previously reported in [155]. The measurement agrees with the proposed mechanisms, although the dependence on these specific mechanisms is not explicitly proven. The result may be negatively affected by spillage through the mold gap, but the amount of material lost by spillage is typically rather small compared to the bulk material. Only the rising flank is considered relevant for compression molding simulations, because the loading in this application is monotonic. However, the overall behavior is irreversible, as seen in Figure 7.5a. Finally, the temperature dependence and rate dependence of the paste viscosity may introduce an error to this quasi-static analysis, if the ability to close pores is affected by resin flow.

### **7.1.4 Friction at mold surface**

The friction between mold and SMC is an important factor in any SMC compression molding model and is introduced in Section 4.1.4. It can be determined experimentally from pressure differences in a press rheometer, which has been introduced in Section 4.2 and Figure 4.4. This section reports pressure results from molding trials with a press rheometer and discusses several friction models in a parametric study with the one-dimensional model from Section 4.2.

**Method.** The press rheometer is equipped with up to 13 Kistler 6167ASP1.6 pressure sensors, which can measure pressures up to 200 bar at mold temperatures up to 200 °C. The placement of these sensors in the mold is illustrated in Figure 7.6 and the outer dimensions of the mold are 800 mm x 450 mm.

**Figure 7.6:** The pressure sensors P1 to P10 are placed along the center line at 32 mm, 146 mm, 248 mm, 350 mm, 450 mm, 552 mm, 604 mm, 654 mm, 709 mm and 764 mm. The sensors P11 to P13 are positioned in a perpendicular line at 500 mm.

Trials are performed in collaboration with Sergej Ilinzeer at Fraunhofer ICT, Pfinztal, Germany with stack dimensions 600 mm x 450 mm (3 layers, 75% mold coverage) and 400 mm x 450 mm (4 layers, 50% mold coverage) on a

Dieffenbacher COMPRESS PLUS DCP-G 3600/3200 AS hydraulic press with parallel cylinder control. The mold is heated to 145 °C for all trials and is closed with a constant speed *<sup>h</sup>*˙ *<sup>=</sup>* °1mms°<sup>1</sup> until the compression force reaches 4400 kN. After reaching the maximum compression force, the mold is kept closed for 120 s to cure the resin and plates are subsequently demolded. The pressures of sensors, displacement of hydraulic cylinders and compression force are recorded during the trials. The recorded data of five specimens in each configuration is aligned at the switch-over point to force control. This point is *t =* 0 in the evaluation, hence *t <* 0 means that the press operates in displacement control and *t >* 0 means that the press tries to keep the force to a constant value of 4.4 MN.

The press rheometer is modeled by the one-dimensional model developed in Chapter 4, which describes the SMC as an anisotropic, compressible, non-Newtonian, non-isothermal material with the parameters determined in the previous sections. All parameters of the numerical model except for the friction model are determined independently from the press rheometer. Hence, the one-dimensional model is used to efficiently vary frictional parameters to find a qualitatively and quantitatively fitting model. The objective is to find a single parameterization that can describe the effects observed in the 50% coverage and 75% coverage trials.

**Results.** First, the friction coefficient of a hydrodynamic power-law friction model is varied in the range 1.0 MN s m−<sup>3</sup> to 4.0 MN s m−<sup>3</sup> with a constant power-law coefficient *<sup>m</sup> <sup>=</sup>* 0.6 and reference velocity *<sup>v</sup>*<sup>0</sup> *<sup>=</sup>* 1mms°1. These are typical values reported for SMC in literature measured after the squish phase for stable plug-flows [3, 147, 148, 150].

The simulated press data and experimental press data of the 75% mold coverage configuration are shown in Figure 7.7. Regardless of friction parameters, the computed final gap is larger than the experimentally observed thickness. The initial compression force is similar for all tested friction coefficients, but the final time at which the virtual press controller switches to force control is earlier for higher hydrodynamic friction coefficients. In the force controlled

**Figure 7.7:** Press data and pressures at 75% mold coverage for different hydrodynamic friction coefficients *∏* - experiment (light colors) and 1D model (dark colors).

phase, the numerical compression force oscillates. The numerical compression forces for higher friction are closer to the experimental force and the earlier switching point leads to a reduced slope of the gap, which resembles

the experiments well, except for the different final thickness. However, compression force and mold displacement alone do not allow for an unambiguous calibration of the friction parameters. The compression force is a holistic measure that does not distinguish the contributions from friction and the behavior of the SMC itself. Therefore, the pressures at the sensor positions of the model are compared to the measured pressures, too. At *<sup>∏</sup> <sup>=</sup>* 1.0MNsm°3, the spread between the highest pressure and lowest pressure is much smaller than in the experiments and the last sensor is triggered too early. An increase to *<sup>∏</sup> <sup>=</sup>* 2.0MNsm°<sup>3</sup> increases the spread, but still does not result in a sufficiently large initial slope of the leftmost sensors. For *<sup>∏</sup> <sup>=</sup>* 3.0MNsm°<sup>3</sup> and *<sup>∏</sup> <sup>=</sup>* 4.0MNsm°3, the computed pressures agree well with the experiments in terms of spread, slope, absolute pressure level and compression time.

The results for 50% mold coverage are depicted in Figure 7.8. The resulting gap and part thickness agree and larger friction coefficients result in a displacement profile closer to the experiment. The compression force increases with a similar slope for all tested friction parameters during the initial squish phase. Beyond that, the compression force is closer to the experimental value for larger values of *∏*, but does not describe the initial hump in the force profile. The compression time from first contact to the switch over point and complete fill agrees well with the experiment for *<sup>∏</sup> <sup>=</sup>* 3.0MNsm°<sup>3</sup> and *<sup>∏</sup> <sup>=</sup>* 4.0MNsm°3. The pressure plots corroborate the observations of press data. The hydrodynamic friction model is not able to predict the initial pressure peak of the experiment with 50% initial mold coverage. However, the solution with *<sup>∏</sup> <sup>=</sup>* 3.0MNsm°<sup>3</sup> predicts the slopes correctly in the force controlled phase (*t >* 0). The time at which each pressure sensor is triggered agrees with experiments as well, which indicates a correct compressibility model.

The hydrodynamic friction model does not apply any resistance to a stack motion relative to the molds at absent relative velocity. The initial hump in the experimental data might be caused by such a sticking behavior, though. Therefore, also constant friction is implemented in the one-dimensional

**Figure 7.8:** Press data and pressures at 50% mold coverage for different hydrodynamic friction coefficients *∏* (*ø*<sup>0</sup> *=* 0Pa) - experiment (light colors) and 1D model (dark colors).

model from Section 4.2 (see also Equation (4.49) for the corresponding analytical solution) and varied between 30 kPa and 200 kPa. A constant friction value of 30 kPa has been previously employed successfully to describe the friction of the considered SMC material [6]. The resulting simulated press

**Figure 7.9:** Press data and pressures at 50% mold coverage for different constant friction stresses *<sup>ø</sup>*<sup>0</sup> (*<sup>∏</sup> <sup>=</sup>* 0Nsm°3) - experiment (light colors) and 1D model (dark colors).

data is compared in Figure 7.9. The compression force increases rapidly until the capability of the sticking friction is exceeded and then increases slower with an increase of friction caused only by the increased contact area between SMC and molds. The experimental kink in the compression force curve is between the simulated compression forces with 100 kPa and 200 kPa.

**Discussion.** At 75% initial mold coverage, the computed final gap is larger than experimentally observed gap. As the initial thickness is a fixed value and the compressibility is known from Section 7.1.3, this indicates a loss of material through the mold gap in the experiment, while the simulated mass remains constant throughout the process. The numerical compression force oscillates in the force controlled phase, as the solver uses large time steps in this phase, which makes it more difficult for the controller to maintain a constant force.

The simulation model with 50% initial mold coverage predicts a monotone increase of pressure from the flow front to the other end of the mold due to its plug-flow assumption. This is in disagreement with the experiment, where the inner layers get squeezed out of the stack and hence lower pressures are observed at the first sensor position (green curve) compared to the second sensor (orange curve). This will be discussed in more detail in Section 8.3.

The constant friction allows to modify the compression force at which the flow starts, but at high constant friction values, the mold would need much more time to fill than experimentally observed. This can be seen by a flat displacement slope in the force controlled phase. The pressure sensor results show how a constant friction is able to describe the initial pressure level. However, without reduction of this value after the initiation of flow, the pressure spread is vastly overestimated.

This parameter study leads to the conclusion that both, an initial sticking behavior with *ø*<sup>0</sup> *=* 100 kPa to 200 kPa as well as a hydrodynamic friction, are required to model the full SMC compression molding process including squish flow and a stable plug flow. Hence, the friction model from Equation (4.13) is parameterized according to Table 7.5.


**Table 7.5:** Mold friction parameters

## **7.2 Patches**

The unidirectional carbon fiber reinforced patches are produced from ZOLTEK UD300 stitch bonded non-crimp fabric UD300 with 50K PX35 rovings at Fraunhofer ICT, Pfinztal, Germany [170]. Compared to the SMC formulation, the composition (see Table 7.6) is designed such that the B-stage is reached faster. This enables direct cutting after production of the semi-finished patches without a long maturing phase. Like the SMC production, the manufacturing process of semi-finished patch material is not subjected to quality control. Local dry spots due to insufficient impregnation, wet spots due to insufficient mixing of the components and uneven patch thickness distribution occur frequently in the semi-finished material. This is accompanied by a dependence on individual batch production conditions and the storage time until processing or testing. These conditions lead to large scatter of the semi-finished patches. Nonetheless, an attempt is made to approximately determine the properties of this material.

### **7.2.1 Thermal properties**

Similar to the SMC, transverse heat conductivity *∑*, specific heat capacity *c*<sup>p</sup> and thermal gap conductance *k*<sup>T</sup> are required to solve the balance of inner energy in patches. The parameters are determined analogously to the


**Table 7.6:** Composition of UPPH patches

procedure described in the previous Section 7.1 at Fraunhofer ICT, Pfinztal, Germany.

**Method.** Six patches of 50 mm x 50 mm x 0.45 mm are stacked and thermocouples are placed centrally between each layer. The first sensor *T*<sup>1</sup> is placed on the surface of the stack, such that it is positioned between mold and stack as soon as the stack comes in contact with the mold. The stack is then surrounded with glass wool insulation and placed on a mold surface heated to 145 °C. A small weight on top of this configuration is used to lightly press the stack against the mold surface.

**Results.** The measured temperature profiles are shown in Figure 7.10 with light colored lines. The resulting properties of the best fit and the specific heat capacity from rule of mixtures are listed in Table 7.7. The solution of the one-dimensional heat transfer Equation (7.1) with these optimal parameters is plotted in dark colors in Figure 7.10.

**Discussion.** The heat transfer from mold to the stack is faster than in SMC, but the transverse conduction within the stack is similar to the SMC. A possible source of error is the adhesion quality of the patches, which is not as good as the adhesion between SMC sheets. Therefore, the transverse heat

**Figure 7.10:** Measured temperatures (light colored lines) and best fit based on one-dimensional heat transfer equation (dark colored lines) for sensors T1 to T6.


**Table 7.7:** Thermal properties of UPPH-CF UD patches in B-staged state.

conductivity could be underestimated due to imperfect interface bonding between patch plies.

### **7.2.2 Tensile stiffness and strength**

The patch model in Chapter 6 requires the Young's modulus in carbon fiber direction *E*0 <sup>F</sup> and in stitching yarn direction *E*<sup>F</sup> as well as the matrix Young's modulus *E*M. Further, the corresponding strengths *Y* <sup>0</sup> <sup>11</sup>, *Y*<sup>11</sup> and *Y*<sup>22</sup> are needed to parameterize the model. Hence, these values are approximated via tensile tests at different temperatures in this section.

**Method.** Pre-impregnated stitched specimens with dimensions 20 mm x 160 mm x 0.45 mm are tested under tensile loading with a clamping length of 80 mm at KIT IAM-WK, Karlsruhe, Germany. At room temperature, specimens are tested on a ZwickRoell Z2.5 universal testing machine with 3 mm min−<sup>1</sup> and 300 mm min−<sup>1</sup> deformation velocities. For specimens with fiber orientation in loading direction, testing is started after applying a tension load of 100 N and the stiffness is evaluated at 0.3% strain. For specimens with fiber orientation perpendicular to the loading direction, testing is started after applying a tension load of 2.5 N and the stiffness is evaluated at 0.3% strain. Five specimens were tested in each configuration, but only four experiments can be considered valid for transverse testing, as other specimens experienced damage within the clamping area (see Figure 7.11).

**Figure 7.11:** Tensile test with fibers oriented perpendicular to the loading direction are considered valid if failure occurs outside the clamped region. If progressive damage occurs in the clamped region, the test is considered invalid.

Figure 7.10 shows that the top part of the patch closest to the mold surface (orange sensor) reaches a temperature of 70 °C to 90 °C after 10 s. Only the sensor in direct contact to the mold exceeds this temperature. As co-molding processes are typically shorter than 10 s, a temperature of 80 °C is considered as the upper limit for a homogeneously heated patch for characterization. At this temperature, specimens are tested on an INSTRON ElektroPuls E3000 universal testing machine equipped with a temperature chamber. The same deformation rates and initial forces as for room temperature are applied.

**Results.** The results of five specimens per rate tested at room temperature and in fiber direction are shown in Figure 7.12a. The strength in fiber direction exceeds the capabilities of this machine. The resulting stiffness with perpendicular fiber orientations at room temperature is shown in Figure 7.12b and the strength is shown in Figure 7.12c.

**Figure 7.12:** Tensile stiffness and strength at 20 °C.

The resulting tensile stiffness in fiber direction at 80 °C is shown in Figure 7.13a. Tensile stiffness and strength at 80 °C for fibers oriented perpendicular to the load direction are shown in Figure 7.13b and 7.13c, respectively.

**Discussion.** The tensile stiffness in fiber direction does not depend on deformation rate and temperature. It is considered a constant value of *E*0 <sup>F</sup> *=* 14GPa subsequently. The tensile strength in fiber direction is high enough to consider fracture in the direction of carbon fibers negligible. The tensile stiffness perpendicular to fibers shows no statistically significant dependence on the deformation rate, but it is temperature dependent. The summed stiffness in perpendicular direction is approximately *E*<sup>F</sup> *+E*<sup>M</sup> *=* 1.3GPa at 20 °C and 0.4 GPa at 80 °C. To parameterize the model and separate the yarn stiffness from the matrix stiffness, some rough approximations are made: It is assumed that the transverse stiffness at higher temperatures is dominated by

the stitching yarn and hence the stitching yarn stiffness is approximated by *E*<sup>F</sup> *=* 0.3GPa. Similar to the carbon fibers, this fiber stiffness and strength is assumed to be independent of temperature, and the strength is conservatively estimated as 3 MPa. The matrix is assigned a tensile stiffness of 1.0 GPa and strength of 16 MPa at 20 °C, which are linearly discounted by a factor of 1/10 for an increase to 80 °C. This is certainly just a very rough estimation of the tensile properties, but large scatter has been observed between individual batches of the patch material depending on the impregnation quality and storage conditions. A robust quality control in the production process is necessary to determine exact properties. Nonetheless, the obtained values give a reasonable estimation to demonstrate the co-molding model. The set of parameters in Table 7.8 is applied in later computer simulations for the tensile parameterization of the patches.

### **7.2.3 Bending stiffness**

Bending in a rheometer enables access to the bending properties, if bending stiffness and membrane stiffness are decoupled, as observed in many prepregs [199, 204]. This section investigates whether such a decoupling has


**Table 7.8:** Tensile properties used for simulation of carbon fiber UD patch in B-stage

to be considered for the deformation in the SMC co-molding process, where relatively stiff patches are placed initially at room temperature into the mold.

**Method.** A photo of the bending apparatus in an Anton Paar MCR 501 rheometer, which is applied for this investigation at Fraunhofer ICT, Pfinztal, Germany, is shown in Figure 7.14a. The specimen dimensions are 25 mm x 50 mm x 0.45 mm. They are subjected to nominal rotation rates of 0.1 min<sup>−</sup>1, 1 min−<sup>1</sup> and 10 min−<sup>1</sup> at 25 °C as well as 0.1 min−<sup>1</sup> and 10 min−<sup>1</sup> at 80 °C. The prescribed deflection angle *Æ*<sup>b</sup> ranges from 0° to 60°. The kinematics of the apparatus lead to a nominal curvature of the specimen

$$\kappa\_{\rm b} = r\_0 \cos a\_{\rm b} + r\_0 (1 + \cos a\_{\rm b}) \cot a\_{\rm b} \tag{7.4}$$

with the distance from clamps to the pivot point *r*<sup>0</sup> *=* 9.5mm (details may be found in Appendix A.1). The Young's modulus for bending may then be approximated by beam theory

$$E\_{\rm F;b} = -\frac{M\_{\rm b}}{I\_{\rm b} \varkappa\_{\rm b}}\tag{7.5}$$

with the bending moment *M*<sup>b</sup> and the moment of area *I*<sup>b</sup> for deflection angles *Æ*<sup>b</sup> 6*=* 0.

**Results.** The results are shown in Figure 7.14b. The initial large values are caused by the singularity in Equation (7.5). The subsequent increase of values is the build up of tension closing the clearance between the apparatus and specimen, which was larger during the trials at 80 °C. The maximum and minimum Young's modulus determined from tensile testing are added as gray dashed lines for reference.

**(a)** Photo of the bending setup **(b)** Results of the bending tests **Figure 7.14:** Rheometer bending test of the patch material. The solid lines in (b) indicate results at 25 °C, the dashed lines indicate results at 80 °C.

**Discussion.** The approximated stiffness at 25 °C is even larger than during the tensile tests. This can be attributed to either friction in the apparatus, or the fact that these specimens come from a different batch. The stiffness at 80 °C is about 5 MPa, which is rather an underestimation due to the play in the apparatus during these tests, which leads to an overestimation of the curvature *∑*<sup>b</sup> in Equation (7.5). In any way, the bending stiffness is in the same order of magnitude as the tensile stiffness. The conclusion is that the behavior of these patches seems not decoupled and can be described by conventional shell elements.

### **7.2.4 Shear stiffness and strength**

In addition to the tensile properties, the shear stiffness *G*M, shear strength *Y*12, shear viscosity *¥*<sup>M</sup> and fracture energy *W*<sup>M</sup> are required to complete the patch model from Chapter 6. These parameters are obtained in this section from bias extension tests.

**Method.** The shear stiffness and strength are obtained from bias extension tests at KIT IAM-WK, Karlsruhe, Germany with a setup similar to that used by Schirmaier [205]. Specimens with dimensions 160 mm x 80 mm x 0.45 mm and with 45° fiber orientation are prepared on an automated cutting table by Zünd Systemtechnik AG, Altstätten, Switzerland. They are subjected to 5 mm min<sup>−</sup>1, 50 mm min−<sup>1</sup> and 500 mm min−<sup>1</sup> deformation rates on a Zwick-Roell Z2.5 universal testing machine at room temperature (see Figure 7.15). Additional tests are performed at 80 °C on a ZwickRoell Zmart.Pro universal testing machine equipped with a temperature chamber at 5 mm min−<sup>1</sup> and 400 mm min−<sup>1</sup> deformation rates. The free length is 135 mm and adhesive tape is employed to improve the fixation in the clamped region.

The deformation is measured via Digital Image Correlation (DIC) with an unstructured pattern of white marker points. The marker points are painted on the rovings of the patches to enable an evaluation of the deformation between individual rovings (see Figure 7.16a). Contrary to dry fabrics, which can be evaluated with a regular grid [205], the points are distributed in an unstructured grid, as rovings are deformed from the manufacturing process of the patch and not equally spaced. The deformation is recorded with a Canon EOS 70D DSLR camera at 30 frames per second with a resolution of 4 pixels mm<sup>−</sup>1. The gray-value image is seeded with tracking points based on regional brightness maxima of the image (see Figure 7.16b). A MathWorks Matlab DIC tool [206] is used to compute the deformation of these tracking points with a correlation window size of 21 pixels x 21 pixels. The resulting deformation of the tracking points is then imported to Paraview for further processing. A Delaunay triangulation is employed to generate a mesh between tracked points (see Figure 7.16c). After removing strongly distorted

**(a)** Sketch of the bias extension setup **(b)** Photo of the bias extension setup **Figure 7.15:** Bias extension test. The region without clamped carbon fibers is subjected to shear deformation.

elements, the displacement can be interpolated with shape functions (see Figure 7.16d) and the deformation gradient of each cell is computed. The employed deformation measure is the Green-Lagrange strain (see Equation (2.80)), which is computed from an element-size weighted average of all elements.

**Results.** The resulting shear stress at room temperature is shown in Figure 7.17a for all tested samples. It is plotted over the Green-Lagrange shear strain obtained from DIC.

The scatter between individual specimens in terms of stiffness and especially in terms of strength is large. There is a statistically significant rate dependency with a steeper initial slope at higher deformation rates (see Figure 7.18a), but not distinct rate effect on strength (see Figure 7.18b). The results at 80 °C are shown in Figure 7.17b for fewer samples. There is also large scatter between samples and no distinct rate effect (see Figure 7.19).

**(a)** Marker points **(b)** Initialization **(c)** Mesh **(d)** Evaluation **Figure 7.16:** Steps for DIC analysis. (a) The specimen is prepared with white marker points. (b) The tracking points are identified in the initial frame. (c) The resulting deformed points are interpolated by Delaunay triangulation (d) The displacement field is interpolated with shape functions and deformation measures are evaluated on the evaluation mesh.

**Figure 7.17:** Shear test results. The deformation of the highlighted dark colored specimen is displayed in Figure 7.20.

The failure at 20 °C varies between abrupt fracture in well impregnated patches to a progressive damage in patches with lower impregnation quality. A representative sample for such a progressive failure is highlighted in dark color in Figure 7.17a and the deformation is shown in Figure 7.20. The displayed strain is the maximum principal value of the Green-Lagrange strain q (*E*xx °*E*yy)<sup>2</sup> *<sup>+</sup>E*<sup>2</sup> xy obtained from DIC. The specimen fails under shear load first. This corresponds to the first peak in Figure 7.17a and can be seen as a

**Figure 7.19:** Shear stiffness and strength at 80 °C.

shift between bundles at the localized shear zone in Figure 7.20b. Then the stitching yarn fails progressively with the formation of a single shear band parallel to the rovings.

**Discussion.** The shear stiffness at 5 mm s−<sup>1</sup> test velocity is 0.110 GPa at 20 °C and increases apparently with increased deformation rates. For 500 mm s−<sup>1</sup> test velocity, the rate *<sup>E</sup>*˙12 is nominally 0.06 s−<sup>1</sup> and at *<sup>E</sup>*<sup>12</sup> *<sup>=</sup>* 0.01 the shear stress is º 3.5 MPa. This suggests a shear viscosity of approximately 10 MPa s with the model in Equation (6.4). The value *¥*<sup>M</sup> *=* 10MPas is subsequently used as a rough approximation of the matrix viscosity of the patch without considering shear-thinning effects. The fracture energy of 5 mJ mm−<sup>2</sup> is estimated from the enclosed area under the stress-displacement curve. Similar to the tensile properties, all matrix properties are simply linearly discounted by a factor 1/10 for a temperature increase from 20 °C to 80 °C. The results may

**Figure 7.20:** Experimental deformation at 20 °C for different crosshead travel states. The failure causes local loss of some facets, if these cannot be accurately correlated in a picture after cracks occur.

be negatively affected by the introduction of adhesive tape in the clamping region, even though the evaluation of strains via DIC negates some of this error. The set of parameters in Table 7.9 is used in later simulations for the shear parameterization of patches.


**Table 7.9:** Shear properties used for simulation of carbon fiber UD patch in B-stage

**Simulation.** To verify correct qualitative and quantitative behavior of the parameterized model, the bias extension test is simulated with the patch model described in Chapter 6 and using the material parameters from Table 7.7, Table 7.8 and Table 7.9. The modeling domain is partitioned in clamping regions and free region (see Figure 7.21). The free region is further partitioned with 5 mm wide sections representing the width of rovings in the patch. As this is a characteristic length of the specimen, the entire domain is meshed with

triangular elements of approximately 5 mm edge length. The partitioning ensures that one edge of each element is always oriented in fiber direction.

**Figure 7.21:** Boundary conditions and partitioning of the simulated bias extension test.

The patch clamping is critical and boundary conditions are chosen carefully. The clamping regions are coupled to reference points with kinematic constraints that prohibit any motion in the direction of carbon fiber rovings. However, the transverse motion is not constrained, as the experimental clamping does not reliably prohibit this motion. Instead, the edges marked at the top right and bottom left in Figure 7.21 have constrained horizontal displacement to ensure that the upper clamping does not deviate from the loading axis. The recorded time-displacement curve from the testing machines are then applied to the upper reference point, while the lower reference point remains in a fixed position.

The resulting deformation behavior is illustrated in Figure 7.22. First, the shear band forms and leads to a localized shift between bundles in Figure 7.22b. The position of the shear band is determined by numerical or experimental imperfections. Its formation is observed either at the upper position or at the lower position. Even after complete failure of the matrix material, the stitching yarn still holds both halves of the specimen together (see Figure 7.22c). Finally, the stitching breaks completely and the specimen is separated (see Figure 7.22d). This behavior qualitatively matches the observed experimental behavior shown in Figure 7.20.

**Figure 7.22:** Simulated deformation of a patch under bias extension loading at 20 °C and 500 mm min<sup>−</sup>1. The colormap shows the maximum principal value of the Green-Lagrange strain with pink color representing values that exceed the color range.

The corresponding force-displacement curves are shown in Figure 7.23, where the extracted positions of Figure 7.22 are marked by dots. The simulation agrees with those experimental results showing a progressive damage and the images at the extracted positions show that the same mechanisms apply. The initial slope, maximum force and damage energy are close to the experiments. Evaluations at 80 °C and different rates show similar results and it is concluded that the model described in Chapter 6 and the calibration of parameters from this chapter allow a reasonable prediction of patch deformation and damage initialization.

**Figure 7.23:** Force-displacement curves of the bias extension test at 20 °C and 500 mm min<sup>−</sup>1. Simulation results with the parameterization of this chapter are plotted on top of the experimental results (light colors) in dark colors and the points indicate the position of the snapshots shown in Figure 7.22.

## **8 Validation cases**

In contrast to the previous verification in Chapter 5 and Chapter 6, the validation compares the model predictions to experimental results. Therefore, parts are molded with a wide range of geometries to explore merits and weaknesses of the proposed mesoscale modeling approach. An overview over the validation cases is given in Table 8.1.


**Table 8.1:** Overview over the chosen validation cases.

## **8.1 Double curved dome with isothermal Newtonian matrix1**

The first validation example should show the general ability of DBS presented in Chapter 5 to predict the fiber architecture in a 3D geometry under simplified Newtonian isothermal conditions. Therefore, a small geometry is molded from SMC with different initial stack configurations and the resulting parts are scanned by computer tomography. The real fiber architecture is then compared to DBS and to the results of a simplified version of the macroscopic model presented 4.

## **8.1.1 Experimental setup**

The material under investigation is the glass fiber reinforced UPPH SMC that has been described in Section 7.1. The specimen under investigation is a curved dome with outer dimensions 120 mm x 94 mm and a nominal thickness of 2 mm. The specimen is molded from two different initial stack configurations. The split configuration is molded from two stacks with dimensions 80 mm x 30 mm x 5.3 mm, which are placed at both ends of the mold, as illustrated in Figure 8.1a. This configuration leads to a knit line in the center, where both stacks meet. The bundle architecture in vicinity of this knit line is of particular interest, as it would represent a weak spot in a structural application. The asymmetric configuration is molded from a single stack with dimensions 80 mm x 60 mm x 5.3 mm, which is placed on one side of the mold and enables a longer flow path, as shown in Figure 8.1b. The mold is heated to 145 °C and closed with a Lauffer hydraulic press in collaboration with Lucas Bretz at KIT wbk, Karlsruhe, Germany. The maximum press force was limited to 50 kN.

<sup>1</sup> Parts of this section are based on [195].

**(a)** Split stack configuration **(b)** Asymmetric configuration **Figure 8.1:** The molded part has outer dimensions 120 mm x 94 mm. For the split stack configuration, two SMC stacks are placed at both ends of the mold. For the asymmetric configuration, a single stack is placed on one side of the mold.

The molded samples were analyzed by volumetric imaging using an Yxlon X-ray µCT system with a Perkin Elmer flat panel Y.XRD1620 detector and a reflection tube by Comet. A volumetric gray scale image is reconstructed based on the Feldkamp, Davis and Kress (FDK) algorithm [207] with a resulting voxel size of 69 µm. µCT scans and volumetric image reconstruction were performed by Ludwig Schöttl at KIT IAM-WK, Karlsruhe, Germany.

## **8.1.2 Numerical setup**

The experimental process is modeled with DBS and a macroscopic model according to Chapters 4 and 5. The molding domain is represented by Eulerian elements and only those Eulerian elements occupied by initial stack positions are initially filled with material.

Both mold halves are represented by rigid shell elements. They interact with the SMC paste through normal contacts and hydrodynamic friction (see Equation (4.11)). While the lower mold is constrained at a fixed position, the upper mold is closed with the profiles given in Figure 8.2. These profiles are an idealization of experimentally recorded press profiles. The variation in the experimental profiles can be attributed partly to a non-uniform thickness of SMC sheets and to the reaction time of the press control unit. The simulation stops after a complete fill with the final part height and does not include the subsequent holding and curing process.

**Figure 8.2:** Distance between upper and lower mold during the compression flow phase of SMC. Six parts of the split configuration were produced and are shown with light green lines. Four parts of the asymmetric configuration were produced and are shown with light orange lines. Additionally, the idealized mold profiles for simulations are shown in solid green and solid orange for split stack configuration and asymmetric configuration, respectively. [195]

The macroscopic reference solution was obtained with a simplified (Newtonian, no anisotropic interaction terms) version of the model described in Chapter 4 that resembles the state of technology. The initial fiber orientation is described by a planar isotropic fiber orientation tensor.

For the DBS, bundles are generated with the procedure described in Section 5.1.2, assuming a planar isotropic orientation distribution in the stacks. The total number of bundles per stack is 7600, which is computed from the bundle properties and a nominal bundle volume in the part of 5844 mm3 at the given nominal fiber volume fraction. Each bundle is discretized with ten linear truss elements.

The explicit time integration requires a small time increment due to the high resin viscosity. The mass of the entire model was therefore scaled by a factor *∑*<sup>m</sup> to improve the time increment, while ensuring that kinetic energy remains negligible small compared to the external work. A summary of the simulation parameters is listed in Table A.6.

## **8.1.3 Results and discussion**

### **8.1.3.1 Flow front progression**

Figure 8.3 provides an overview on the compression molding process simulation of the split stack configuration. The initial mold gap at *t =* 0s is 20 mm and the lowest part of the upper mold is just not touching the SMC stacks. Closing the mold initially deforms the stacks into the curved geometry, but does not start material flow. During forming, the two-way coupled approach pulls the stack sideways in the dome-shaped mold. This can be observed by the lateral deformation of the stack tips depicted at *t =* 2s in Figure 8.3. The mold gap is reduced to the initial stack height of 5.3 mm after approximately two seconds. From there on, flow dominates the mold filling process and fiber bundles are carried with the SMC until the final part thickness of 2 mm is reached.

### **8.1.3.2 Orientation and separation effects**

Figure 8.4 shows slices through the midplane of the upper and lower planar regions of the scanned part in split stack configuration. The white strands represent fiber bundles, which remain in their bundled structure even for the applied high degree of deformation. The knit line features a severe fibermatrix separation and only a small amount of fiber bundles bridges the gap in this zone. The inner slice in Figure 8.4 even shows some pores. Regions close to the mold boundaries and the knit line show a bundle alignment parallel to the boundary. Bundles perpendicular to the boundary are likely pulled out of

**Figure 8.3:** Snapshots of the molding process for the split stack configuration. The compression molding process starts with a deformation of the two initial stacks. Subsequently, the SMC is forced to flow until the part reaches its final thickness of 2 mm. The DBS approach lets bundles deform and flow with the matrix material while enforcing two-way coupling. Therefore, the flow is naturally anisotropic and depends on the current bundle configuration. [195]

this region by forces acting over the entire length of the bundle and parallel bundles remain close to the boundaries (see also Figure 2.9). Regions farther away from boundaries show a regular random in-plane orientation.

The DBS result is sliced in the same planes and the result is depicted in Figure 8.5. The simulation results show a slightly larger area of fiber-matrix separation and bundles do not bridge the resin-rich knit line. This might be caused either by the experimental setup, because the part was compressed further than the nominal thickness, or by an underestimation of the drag forces on bundles in the model. Similar to the µCT-scan, boundary regions show a predominant orientation parallel to the boundaries.

**Figure 8.4:** Slices through the upper and lower planar regions of the µCT Scan. Fiber bundles stay intact during molding and fiber-matrix separation can be observed at the knit line. The knit line region includes pores marked with red circles. [195]

For a quantitative comparison of the DBS to the macroscopic simulation and the µCT scans, bundle orientations are evaluated on a uniform 12 x 16 x 4 grid of sub-volumes. The discrete second-order fiber orientation tensor for each of the sub-volumes is computed with the procedure described in Section 5.1.6.3. The slices of the µCT scan shown in Figure 8.4 are analyzed in 2D using the OrientationJ plugin of the image processing software Fiji [208]. This procedure assigns a major direction to each 10 x 10 pixel area. The discrete fiber orientation tensor is evaluated on an equivalent 12 x 16 grid to represent the orientation state as tensor components.

A comparison of the DBS approach, µCT scan and the conventional macro fiber orientation model is depicted in Figure 8.6 for the split stack configuration. The *A*xx-component of the µCT-analysis features three significantly higher oriented vertical stripes at both ends of the mold and the knit line. Conversely, the *A*yy-component of the µCT-analysis indicates a dominant orientation in horizontal direction at the top and bottom mold boundaries with lower values at the vertical mold boundaries to the left and right of the

**Figure 8.5:** Slices through the planar regions of the DBS result. Each gray cylinder represents a bundle segment consisting of 200 individual fibers. The knit line at the center is matrix rich and no bundles gap this region. Bundles close to the boundaries show a reduced fiber volume fraction and more bundles are oriented parallel to the boundary. [195]

figure. The corresponding DBS is able to reproduce these three stripes of higher vertical orientation at the correct positions. Characteristic gradients and the level of orientation is predicted well. The macroscopic model using fiber orientation tensors and Jeffery's equation does not account for the constraints at mold walls and shows a homogeneous orientation distribution. In homogeneous regions, such as the inner slice with some distance to the knit line, Jeffery's equation leads to a reasonable prediction of the orientation state.

The DBS limits any bundle orientation normal to the molds, because bundle segments cannot be physically arranged in normal direction in the constrained mold gap. Thus, the *A*zz-component is small in the planar regions of the part. An investigation of a magnified µCT scan with higher resolution confirms that fiber bundles at the knit line are primarily oriented in-plane. The computation based on fiber orientation tensors shows a dominant normal component of fiber orientation at the knit line.

**Figure 8.6:** Comparison of DBS results with µCT Analysis and fiber orientation tensor based computation utilizing Jeffery's equation for the split stack configuration. The first row shows orientation tensor component *A*xx which indicates vertical fiber orientation in this representation. The second row shows the *A*yy-component representing horizontal fiber orientation. The third row shows the *A*zz-component representing fiber orientation normal to the observation plane. The orientation analysis of the µCT image slices is limited to two dimensions. Thus, the central image in the third row shows a high resolution µCT scan of the region indicated in the illustration above. The magnified view reveals a dominant in-plane orientation of bundles. [195]

Figure 8.7 is analogous to Figure 8.6, but describes the evaluation of the asymmetric stack configuration with a maximum flow path of 60 mm in *y*direction. This configuration confirms observations of the previous case with significantly higher orientations parallel to mold walls that are not described by tensor based theory. Despite a longer flow path, the magnitude of reorientation is similar to the split stack configuration due to a similar stretch in *y*-direction (50% initial mold coverage each).

**Figure 8.7:** Comparison of DBS results with µCT Analysis and fiber orientation tensor based computation utilizing Jeffery's equation for the asymmetric stack configuration. Refer to Figure 8.6 for a detailed explanation of the layout. [195]

### **8.1.3.3 Bundle curvature**

The curvature of bundles is evaluated with the procedure described in Section 5.1.6.1 at each node *k*. A color-coded scatter plot of the curvature for the split stack configuration is plotted in Figure 8.8. It shows that the largest curvatures occur at corners and close to the knit line. The curvature at the knit line originates probably from a flow in *x* direction compressing bundles to a zigzag shape. The curvature in the µCT scan is obtained only for the central region in order to have sufficient resolution for tracking bundle curvature [209].

Curvature of each node *ck* in mm−<sup>1</sup>

**Figure 8.8:** Simulation results of bundle curvature. The highest values occur at the corners of the mold and at the knit line. The part's three-dimensional shape is visible in this plot due to the bending of bundles at curvatures of the geometry. High-resolution µCT curvature data for the central area marked with a black rectangle is obtained by Ludwig Schöttl at KIT IAM-WK, Karlsruhe, Germany. [195]

The projection of curvature values in the µCT scan region on the *y* direction is plotted in Figure 8.9. The maximal values of the µCT scan agree well with the maximal curvatures computed from the DBS. The mean curvature of DBS slightly underestimates the curvature in the outer regions, but agrees well around the knit line. It should be mentioned that simulated curvature might depend on the segment length of bundles.

#### **8.1.3.4 Number of contacts**

An a priori estimate for the number of contacts per bundle segment is given by equation (2.51). This estimate predicts about 6.2 contacts per bundle segment. The total number of contacts (*g <* 0 in Equation (5.23)) in the

**Figure 8.9:** Curvatures projected to the *y* axis in the scanned central region (see rectangular section in Figure 8.8) encompassing the knit line.

DBS is evaluated for each frame of the simulation results and is plotted in Figure 8.10. This averages to approximately 1.8 × 10<sup>5</sup> contact *pairs* for the split stack configuration and 2.2 × 10<sup>5</sup> contact *pairs* for the asymmetric flow, which has a slightly increased fiber volume fraction compared to the nominal value. Considering the total amount of 77,438 and 87,950 bundle segments, this evaluates to 4.6 and 5.0 contacts per bundle segment, respectively. This evaluation fits well with the estimate given in equation (2.51).

### **8.1.4 Concluding remarks**

The mold filling process of a double curved dome structure is described with DBS using an isothermal Newtonian paste and an isotropic Newtonian reference model. The predicted fiber architecture agrees well with measured µCT data w.r.t. fiber miss-alignment at mold walls, knit line formation and curvature. The flow of material is assumed to be isothermal in this validation case. This assumption is quite common for the bulk material of SMC, as the time scale of thermal diffusivity in SMC is large compared to the time it takes the material to flow (less than 5 s). However, even small changes in temperature can cause a relevant change in viscosity and is therefore

**Figure 8.10:** Number of bundle-bundle contact pairs during the molding process. The number of contact pairs decreases during the forming phase of the stack and increases during flow, when the entire stack is compressed. [195]

introduced to later validation examples. Additionally, lubrication between bundles is neglected in this validation example for computational efficiency, since it has little effect on fiber architecture (see Section 5.2.3). The DBS is able to predict fiber-matrix separation effects at the knit line and thus enables a better description of structural weak spots in such areas. At regions close to the mold walls and the knit line, DBS accounts for spatial constraints of the fiber orientation due to mold boundaries and leads to more accurate fiber orientation results than orientation tensor based macroscopic simulation models. Nonetheless, Jeffery's equation leads to reasonable results in planar, homogeneous regions at up to ten times faster computational times.

## **8.2 Honeycomb ribs with isothermal Newtonian matrix**

Fiber-matrix separation is investigated in a plate mold with honeycomb shaped ribs. The experiments were conducted as part of the Master's thesis by Florian Rothenhäusler at Volkswagen AG, Wolfsburg, Germany using a VE SMC with 50 mm fiber length [63].<sup>2</sup> This section reports the application of the DBS model to this geometry and material system. The resulting fiber-matrix separation and compression forces are compared to the experimental results.

## **8.2.1 Experimental setup**

The material in this example differs from the default material described in Section 7.1. Here, the commercial SMC HUP EJ 43529 by Polynt Composites Germany GmbH is used for compression molding. It has a fiber volume fraction of 36% and fibers are 50 mm long with a fiber diameter of 20 µm. The mold is a plate with four insert positions, of which one is filled with a rib structure consisting of *thick* ribs (2.5 mm width at top, 2° draft angle), one is filled with a *thin* rib structure (2 mm width at top, 1° draft angle) and the remaining two insert positions are filled with dummy plugs. The SMC is molded in the two different configurations shown in Figure 8.11 on an RSP 160.42SZ hydraulic press by Röcher GmbH&Co KG Maschinenbau, Germany. In configuration A, the stack is placed directly under the honeycomb ribs and the ribs are filled at the beginning of the compression. In configuration B, the stack is placed opposite to the honeycomb ribs and the material flows along the mold before filling the ribs. Both configurations use a mold closing speed of 10 mm s<sup>−</sup>1, mold temperature of 140 °C and maximum compression force of 1500 kN. Press forces, velocity and displacement are recorded during trials.

The fiber volume fraction in the ribs is analyzed by thermal gravimetric analysis (TGA). To limit the experimental effort, only the ribs in the position marked with labels 1 to 6 in Figure 8.11 are analyzed. This position is chosen, because these ribs show a challenging fiber-matrix separation pattern, but are also reliably filled during the process. Each of the six walls is subdivided in a lower, central and upper partition and TGA is performed following procedure

<sup>2</sup> Parts of the Master's thesis by Florian Rothenhäusler are included here for comparative purposes. They will be published soon under the title *Experimental and Numerical Analysis of SMC Compression Molding in Confined Regions — A Comparison of Simulation Approaches*.

B of DIN EN ISO 1172:1998. The specimen is pyrolyzed for 30 min at 650 °C and subsequently the mineral filler is dissolved by hydrochloric acid [63]. This procedure is repeated three times for each configuration (324 samples in total) and the fiber volume content is calculated from the fiber mass after the extraction procedure.

**Figure 8.11:** Mold configurations of the honeycomb experiments.

In addition to the fully molded samples, short shot moldings in configuration A were produced by Simon Wehler at Volkswagen AG, Wolfsburg, Germany. The short shots were obtained by inserting 1.6 mm and 4.8 mm thick rigid spacers in the mold gap prohibiting complete closure of the mold. The surface of the resulting parts were scanned with a 3D optical measurement system ATOS by GOM GmbH, Braunschweig, Germany.

### **8.2.2 Numerical setup**

The process is modeled with the Direct Bundle Simulation technique explained in Chapter 5. The upper and lower mold surfaces are extracted from the CAD model of the part, placed at an initial distance of 12 mm and are meshed with rigid body elements of 2 mm average edge length. The lower

mold remains fixed, while the upper mold is either controlled by a virtual press controller (see 4.1.5) following the press profile (configuration B), or subjected to a constant force (configuration A). Configuration A is simulated with a constant compression force, because the press comes effectively to a stop after reaching the stack and proceeds in a force controlled mode in the trials with this configuration. This is likely caused by the ribs clinging to the SMC stack, which delays the start of the flow process. The Eulerian domain is expanded upward and downward beyond the mold surfaces and meshed with 1 mm x 1 mm x 1 mm EC3D8R elements in the region of the honeycombs and 1 mm x 2 mm x 2 mm else. The in-plane boundaries are not implemented via rigid body contact, but by constraining the velocity in normal direction to zero. The regions of initial stacks are filled with material in the beginning and 17000 bundles of 50 mm are generated in the stack to represent the desired fiber volume fraction. The bundle cross section area is set to 1.768 × 10−<sup>7</sup> m2 assuming 560 fibers of 20 µm thickness each.

The simulation is isothermal and quasi-incompressible due to a lack of detailed material characterization of the material employed in this example. The paste viscosity is set to a constant Newtonian value of 25 kPa s, which was obtained as the zero-shear viscosity for the paste employed in this SMC and applicable for shear rates up to approximately 10 s−<sup>1</sup> [63]. Bundles interact with each other, but tangential lubrication forces are disabled as the resulting time step would become prohibitive. The friction between mold and SMC is modeled with a conventional hydrodynamic friction model that has been parameterized with a press rheometer at Fraunhofer ICT, Pfinztal, Germany following the procedure described by [150]. A summary of all simulation parameters can be found in Table A.7.

Both configurations are simulated three times each, with different initial random fiber bundle distributions. This allows to investigate the effect of the initial configuration on process results and allows to quantify the scatter of computed fiber volume fraction.

## **8.2.3 Results and discussion**

### **8.2.3.1 Flow front progression**

The snapshots of the flow front progression in Figure 8.12 show the rib filling process for both configurations.

**(a)** Configuration A **(b)** Configuration B **Figure 8.12:** Snapshots of the honeycomb mold filling process.

In configuration A, ribs are filled shortly after the start of the compression process. There is a slight tendency for ribs down the flow path to fill later, but the overall filling process occurs quite homogeneous. After 0.8 s simulation time, the ribs are completely filled and the flow front in the plate region flows until it reaches the end of the mold. Long fibers, which are entangled in the ribs, influence the flow and lead to a significant alignment of fiber bundles in flow direction. This can also be seen in Figure 8.13, where the matrix of an entire plate was burned off and is compared to the simulated fiber bundle architecture. Both, the simulation and pyrolyzed part, feature a distinct fiber bundle alignment left to the ribs.

**(a)** Pyrolysis [63] **(b)** Simulation **Figure 8.13:** Comparison of the fiber architecture in configuration A.

In configuration B, the stack elongates first similar to the linear compression flow in a press rheometer. As soon as the flow front reaches the ribbed section, a portion of the flow diverges in the ribs and the ribs closer to the initial stack are completely filled before the ribs further away even begin to fill. The diverged flow causes the formation of a knit line behind the ribs, which can be also seen in the bundle architecture in Figure 8.14.

To validate the flow front progression, the filled domain is compared to 3D optically measured surfaces of short shots of configuration A, which are shown in Figure 8.15. Figure 8.15a is obtained from a specimen with 4.8 mm thick spacers in the mold gap stopping the compression flow soon after the stack is touched. Overlaying the simulation time step closest to this thickness (see Figure 8.16a) shows that the flow front progression is captured

**(a)** Pyrolysis [63] **(b)** Simulation **Figure 8.14:** Comparison of the fiber architecture of the thick honeycomb in configuration B.

accurately. For a later compression state obtained with a 1.6 mm thick blank, the simulation and scan also agree (see Figure 8.16b). The simulation predicts equal flow front progression in both, thin and thick, ribs. However, the thinner ribs appear to be filled slightly faster in the experiment. This is likely caused by higher temperatures in the thin region that result in lower viscosity and a faster filling progress, which is not accounted for in this isothermal simulation. Further, the filling process may be affected by missing parallelism control of the press, or the compressibility of the stack.

**(a)** Short shot 4.8 mm **(b)** Short short 1.6 mm **Figure 8.15:** Surface scans of short shots obtained by Simon Wehler from Volkswagen AG, Wolfsburg, Germany.

**(a)** Short shot overlay 4.8 mm **(b)** Short short overlay 1.6 mm **Figure 8.16:** Comparison of scanned short shots (gray surface) with paste fill region (green volume) of thick ribs in configuration A.

### **8.2.3.2 Fiber volume fraction**

The focus of this evaluation are the thick honeycomb ribs labeled with numbers in Figure 8.11. These are regularly filled and the observations are qualitatively equivalent to the mirrored ribs in the thin honeycomb. An evaluation mesh representing the experimental partitions is used to compute fiber volume fractions of the DBS (see Section 5.1.6.3).

A comparison between thermal gravimetric analysis and the simulation results for configuration A is shown in Figure 8.17. The numbers at the bottom indicate the rib label according to Figure 8.11 and the *±* values indicate the standard deviation between samples at each position. The experimental results show a homogeneous fiber volume fraction distribution in the ribs with small deviation between the three analyzed samples. The ribs, which are placed directly above the initial stack position (1,2,3) exhibit slightly higher fiber volume fractions than the other ones. Contrary, the corresponding simulation predicts severe fiber-matrix separation with local fiber volume fractions well above the experimentally observed values. An investigation of the simulation progress reveals that the ribs fill homogeneously first, but fiber bundles are pulled out subsequently. Fiber bundles entangle at the bottom of ribs with multiple adjacent ribs (e.g. 1,3,6) and lead to a concentration at the

**(b)** Fiber volume fraction obtained by Direct Bundle Simulation **Figure 8.17:** Fiber volume fraction distribution in thick ribs of configuration A.

bottom of these ribs. The discrepancy between simulation and experiment is likely caused by missing shear thinning properties of the matrix material in the model, which lead to smaller forces on the bundle ends extending to the base plate. Another reason could be the gap correction described in Section 5.1.5.3, which is just a coarse approximation of the actual bundle shape.

The experimental result of configuration B predicts a gradient in volume fraction from the top left to the bottom right in Figure 8.18a. The top sections of ribs 1 to 3 exhibit large experimental uncertainties - up to 13.3% standard deviation at a mean value of 16.6%. The fiber architecture in these regions depends strongly on individual realizations of randomness in the mold filling process. The numerical model predicts overall slightly higher fiber volume fractions and severe separation in Rib 2, that occurs reliably in all three repetitions.

**(a)** Fiber volume fraction measured by TGA [63] 1 2 3 4 5 6 Top Center Base 32.6 *±*1.60 9.8 *±*0.60 50.1 *±*9.40 46.8 *±*3.10 44.4 *±*2.60 53.5 *±*3.60 54.0 *±*4.30 24.9 *±*1.20 50.2 *±*4.60 48.2 *±*1.70 37.8 *±*2.00 60.4 *±*3.00 54.6 *±*2.10 42.5 *±*2.30 47.9 *±*2.70 47.7 *±*0.90 44.5 *±*3.40 52.2 *<sup>±</sup>*5.70 0% 25% 50% 75% 100%

**(b)** Fiber volume fraction obtained by Direct Bundle Simulation **Figure 8.18:** Fiber volume fraction distribution in thick ribs of configuration B.

#### **8.2.3.3 Compression forces**

The simulated compression force, velocity and mold gap are compared to press data in order to evaluate the model's ability to predict not just the fiber architecture, but also resulting press parameters. The press used in the experiments is not designed for precise measurements, but for industrial manufacturing processes, and does not have parallelism support. For example, the closing velocity differs from the prescribed value of 10 mm s−<sup>1</sup> during the displacement control and the compression force is not reliably held at the nominal value of 1500 kN. This can be seen from the gray experimental recordings shown in Figure 8.19 for configuration A and in Figure 8.20 for configuration B. The dots in both figures indicate the snapshots positions for an easier comparison to Figure 8.12. In configuration A, the press effectively stops at contact with the stack and tries to proceed force controlled. Hence, the DBS model also utilizes a constant compression force in this case. This results in a rapid initial compression similar to the experimental velocity. However, the subsequent filling phase is faster than the experimental

records indicate, which is partially caused by the lower than nominal compression force in the experiment and partially caused by skipping the initial pre-compression phase. In configuration B, the virtual press controller follows the averaged press velocity profile and switches to force control, after the maximum force of 1500 kN is reached. The resulting displacement during the force control phase is comparable to the experimentally observed displacements.

**Figure 8.19:** Compression press data in comparison to DBS for configuration A. Dots mark the snapshot positions shown in Figure 8.12.

**Figure 8.20:** Compression press data in comparison to DBS for configuration B. Dots mark the snapshot positions shown in Figure 8.12.

### **8.2.4 Concluding remarks**

The molding of a plate with honeycomb shaped ribs has been simulated with DBS employing isothermal, Newtonian paste viscosity. The flow front progression is reasonable and verified by short shorts. The thin honeycomb is filled slightly faster, which can be likely attributed to the neglected temperature change. The predicted fiber volume fraction in configuration A disagrees with the experimental result, as the effect of fiber bundle pull out is too large compared to the experiments. Conversely, the fiber volume prediction in configuration B is able to correctly identify critical regions at the top of Rib 1 and Rib 2. A comparison to other simulation tools (Moldex, 3DTimon DFS) showed that macroscopic models are not able to predict the fiber-matrix separation in this example and that a neglection of fiber-fiber interactions can lead to unreasonably high fiber volume fractions, because multiple fibers can occupy the same space [63]. The consideration of non-isothermal conditions and a non-Newtonian paste seems required to accurately predict compression forces as well as bundle architectures in thin structures.

## **8.3 Press rheometer**

A press rheometer is a plate mold for rheology experiments at component scale and has been introduced in Section 4.2, Figure 4.4 and the characterization of friction in Section 7.1.4. Typically, such a mold is used to characterize the macroscopic elongational viscosity and hydrodynamic friction properties of the stable plug-flow with an array of pressure sensors that is distributed within the mold. The macroscopic elongational viscosity can be measured only in such a large mold because only this way it can be ensured that the fiber length is sufficiently short compared to the dimensions of the measurement instrument. Conversely, the proposed DBS model only requires the paste rheology and thus the press rheometer is used in terms of validation to investigate whether the observed pressures agree with the bottom-up DBS

model. This is done with precise parallelism controlled presses and the full non-isothermal, non-Newtonian, compressible material model.

### **8.3.1 Experimental setup**

The experimental setup is already described in Section 7.1.4 and the experiments are extended by additional trials with the stack configuration 200 mm x 450 mm (8 layers, 25% mold coverage) on a Dieffenbacher DYL630/500 hydraulic press with parallel cylinder control in collaboration with Sergej Ilinzeer at Fraunhofer ICT, Pfinztal, Germany. The part and the initial stack configurations are sketched in Figure 8.21. The pressures of sensors, displacement of hydraulic cylinders and compression force are recorded during the trials. The recorded data of five specimens in each configuration is aligned at the switch-over point to force control. This point is *t =* 0 in the evaluation, hence *t <* 0 means that the press operates in displacement control and *t >* 0 means that the press tries to keep the force to a constant value of 4.4 MN.

**Figure 8.21:** Mold configurations of the press rheometer experiments.

### **8.3.2 Numerical setup**

The SMC flow is dominated by a horizontal one-dimensional elongation and the variance in transverse direction is expected to be rather small. Hence, for computational efficiency, the DBS domain width is reduced to *B =* 50mm, which is twice the fiber length. The length is *X*max *=* 800mm and the inplane Eulerian mesh size is 2.5 mm x 2.5 mm. The Eulerian mesh size in thickness direction (*z*) ranges from 0.6 mm in the region of the final part geometry to 2 mm at the top. The velocity in normal direction is set to zero at all boundaries of the domain. The domain of the initial stack is filled with 40000 fiber bundles with a length of *L =* 25mm, which are discretized by ten truss elements each. Bundles interact with each other and with the mold surfaces through contact forces, but tangential lubrication is disabled due to its prohibitive restriction on the stable time step during explicit integration. The initial temperature distribution is precomputed assuming 10 s contact with the lower mold between stack placement and the upper mold arriving at the stack. The pre-computed temperature and Eulerian mesh profile are depicted in Figure 8.22.

**Figure 8.22:** Initial stack at the flow front. The colors indicate the initial temperature distribution in all elements that are initially filled with material. Both mold halves are heated to a fixed value of 145 °C.

The SMC paste is a compressible, non-Newtonian and non-isothermal material with properties listed in Table A.8. Both mold halves are rigid body shells (element size 2.5 mm x 2.5 mm) with constraint degrees of freedom and a

prescribed temperature of 145 °C. The upper mold's vertical motion is controlled by a VUAMP subroutine that mimics the press controller, as described in Section 4.1.5.

The Eulerian domain is subdivided at the sensor positions to evaluate the pressure as average of the nodal pressures which are filled with material at any given time. The compression force of the reference point linked to the upper rigid body mold is smoothed with a Butterworth filter, to eliminate spurious noise that occurs during the explicit time integration.

## **8.3.3 Results and discussion**

### **8.3.3.1 75% mold coverage**

The recorded press data for 75% mold coverage is depicted in Figure 8.23 in light colors. The data lines represent the mean of five repetitions and the small standard deviation is depicted as light colored area behind the lines. All data is temporally aligned with the switch to force control, such that *t <* 0 refers to the displacement controlled profile and *t >* 0 refers to the force controlled profile. The mold is closed with a constant velocity of 1 mm s−<sup>1</sup> until it is controlled by the constant force resulting in a flatter slope and finally reaches its final thickness. For reference, the final thickness of the part is indicated with a dashed line and shows that the press displacements recordings are well calibrated. However, the final parts are approximately 0.25 mm thinner than the expected thickness from volumetric considerations accounting for the measured compressibility. This is most likely caused by loss of material through the mold gap due to the high pressure experienced in this configuration. The pressure sensors show the expected behavior for a plug-flow with the left most sensor reaching the highest pressures and a decrease in pressure towards the flow front. However, the pressure of the left most sensor might reach its maximum capability of 200 bar and therefore might underestimate the maximum pressure in the mold. The resulting plates are homogeneous without any flow marks.

**Figure 8.23:** Press data at 75% mold coverage - experiment (light colors) and DBS model (dark colors).

The DBS model is used to simulate the compression molding process in the press rheometer with 75% mold coverage and the result is shown in Figure 8.23 in dark colors. As described previously, the numerical model predicts a thicker part than produced in the experiments, which is likely caused by leakage in the experiments. The displacement slope reduces correctly as soon as the virtual press controller reaches the force controlled phase. The time until the switch agrees with the experiments and the compression force is controlled reliably at the prescribed level. The spread between pressure sensors is similar to the experimental results and the general pressure levels are predicted correctly. However, the exact values are not met because of the earlier completion of the filling process in the simulation due to leakage, the maximum sensor capability of 200 bar and idealizations in the model.

#### **8.3.3.2 50% mold coverage**

The experimental results with 50% initial mold coverage are depicted in Figure 8.24 in light colors. The compression force features a hump at which the rate of compression force increase slows down for a moment before reaching the maximum force. This phenomenon can be likely attributed to a transition from a sticking friction to a hydrodynamic lubricated friction. The final thickness of the parts agrees with the measured mold deformation and is in line with the expected thickness from mass conservation considering compressibility.

**Figure 8.24:** Press data at 50% mold coverage - experiment (light colors) and DBS model (dark colors).

The measured pressures exhibit a larger variance between individual plates, which is indicated by the colored background areas representing the difference between the upper and lower quartile of results. Contrary to the results

of 75% mold coverage, the left most pressure sensor P1 does not record the highest pressures and the molded plates show evenly sized zones with disturbed fiber bundles exhibiting higher curvature at both ends of the molds (see Figure 8.25). The separation lines between these zones are displayed as gray lines in the mold coverage sketch in Figure 8.24.

**Figure 8.25:** Photos of exemplary plates molded with 50% initial mold coverage. Both ends feature a darker zone with highly curved compressed fibers.

The results of the DBS simulation with 50% mold coverage are displayed in Figure 8.24 in dark colors. The time to reach the compression force switch over point and the total compression time to fill the mold agree with the experimental data. The initial compression force increase and plateau are not met by the model even after accounting for an initial sticking behavior in the friction model. The used wall stress of 150 kPa should be sufficiently large to cause a steep compression force increase in the beginning according to the previous parametric study with one-dimensional model. An evaluation of the tangential mold stress reveals that this tangential stress is not properly enforced by Simulia Abaqus/Explicit at the reconstructed surface between molds and SMC, because it is only enforced if a surface overclosure is present. However, the contact initially oscillates between a closed and open state in many regions and thus allows a stack motion, whenever the contact is open. Hence, the initial compression force and the initial pressures are underestimated. The DBS does not explicitly account for individual sheets and is therefore not capable to replicate the shear induced stack deformation that

causes the pressure sensor P1 to report pressures below the values obtained by P2.

#### **8.3.3.3 25% mold coverage**

The experimental results for 25% initial mold coverage are depicted in Figure 8.26 in light colors. This data was recorded on a smaller press, which deforms noticeably at the given load. Hence, the plot distinguishes the gap computed from press displacement and a true gap, which is corrected by the stiffness of mold and press. The true gap agrees with the measured final part thickness. The compression force increases in a relatively flat slope until the mold is almost filled and then raises quickly to the maximum compression force.

**Figure 8.26:** Press data at 25% mold coverage - experiment (light colors) and DBS model (dark colors).

The measured pressure does not adhere to the expected behavior of an ideal plug-flow (compare Figure 4.5), as the pressure does not monotonically decrease from the wall towards the mold front. Instead, the leftmost pressure sensor P1 records small pressures and at times the sensor P3 records the highest pressure. This is corroborated by the formation of large zones with disturbed fiber bundles, which are indicated by gray lines in Figure 8.26. Despite the asymmetrical initial stack placement, these zones are approximately even in size at both ends of the mold.

The formation of the zones is likely caused by a shear deformation of the stack. To analyze this effect, a sheet at the center of the nine-layered stack is colored by black paint to track its deformation. The black sheet is clearly visible after molding through the other sheets (see Figure 8.27) and its shape is similar to the central zone observed without coloring. This suggests that the central sheets are stretched less than predicted by the ideal plug-flow model and the upper and lower sheets fill the room at both ends, which originates from the reduced stretch of the central sheets. This mechanism is in line with the observations by Hohberg [158], but it is a stark contrast to the assumptions of the classical plug-flow models. Temperature gradients must affect more than just a very thin lubrication layer at the mold surface.

**Figure 8.27:** Photos of exemplary plates molded with 25% initial mold coverage. The left plate is molded from nine sheets, were the central sheet was colored black.

The DBS result with 25% initial mold coverage is displayed in Figure 8.26 in dark colors. The experimental results of this case were obtained on a smaller

press, which is compliant under the given loads and results in 0.2 mm apparently measured displacement at 4400 kN, even if there is none. Hence, the experimental gap determined from press displacement is corrected for the compliance by an empty reference press stroke. To model the compression process accurately, the virtual press controller mimics the press compliance. The resulting compression time and compression force agree with the experimental data, but similar to the case with 50% initial mold coverage, the initial increase and pressure drop of the left-most sensors due to the stack deformation are not captured.

The orientation evolution in the plate is shown in Figure 8.28. The bundle orientations are mapped to an evaluation cell encompassing the entire domain using the procedure described in Section 5.1.6.3. The resulting orientation tensor represents the global orientation state in the cavity and its components are plotted as solid lines in Figure 8.28. The orientation has an initial preferred direction in the flow direction, because the removal of truss elements at boundaries during the cutting procedure affects this direction less than the transverse direction with its longer boundary. Fiber bundles orient themselves in flow direction during the process, as expected. The solution of the one-dimensional reference model based on Jeffery's equation and the IBOF closure approximation is plotted as dotted lined for comparison. This solution assumes an ideal transversal isotropic initial fiber orientation state and is therefore shifted in comparison to the DBS result. However, except for the shift, this simple model agrees well with the mesoscale model. This suggests that the orientation evolution in a sufficiently large planar SMC region can be described equivalently with this macroscopic model and no further interaction parameters.

### **8.3.4 Concluding remarks**

Press rheometer trials with distributed in-mold pressure sensors allow for an in-depth analysis of the compression molding process that allows to distinguish contributions from compressibility, friction and the elongation of SMC

**(a)** 75% mold coverage **(b)** 50% mold coverage **Figure 8.28:** Orientation evolution for the entire cavity computed by DBS (solid) and with the macroscopic model using Jeffery's model (dotted).

itself. The experimental results show that the SMC under investigation does not follow the ideal plug flow assumption for thick initial stacks. Contrary, the introduction of a colored sheet in the center of the stack proves that the central sheet is not stretched to the entire mold length and the outer layers fill the remaining region with a disturbed bundle structure. The bundle architecture is not just a question of whether it is the region of the initial stack, but the result of complex flow mechanisms that can transport entire sheets to the center of the mold.

The novel mesoscale simulation method can predict the compression time and occurring compression forces accurately with viscous behavior based only on a paste viscosity characterization at coupon scale. This coupon scale characterization is much simpler than press rheometry, but has the disadvantage that the pure paste might not be available for commercial SMC prepregs. The contributions from friction and anisotropic viscosity of the bundle suspension are proportionate. The non-isothermal simulation model accounts for compressibility of SMC, non-Newtonian viscosity, bundle-matrix interactions, bundle-bundle interactions (excluding lubrication), press compliance and press control. The predicted pressure differences between sensor positions and overall pressure levels are predicted correctly, even if the exact pressures of the experiment are not met exactly. Reasons for the difference are

on the one hand experimental deficiencies, such as the maximum pressure capability of the sensors, leakage through the mold gap or the point-wise nature of the measurement, as well as on the other hand limitation in the simulation model, such as no explicit resolution of individual sheets and imperfect contact enforcement of the Simulia Abaqus/Explicit solver. Even though the initial sticking behavior cannot be accurately enforce in the contact, it is concluded that such a sticking friction is needed to model the entire compression process including squish and flow phase. Modeling the pressure drop of the first sensor and the motion of the inner sheet likely requires different fluid phases that represent individual SMC sheets. However, this would require a much finer Eulerian mesh and a solution for merging these phases with progressing simulation time, as heating and deformation remove any distinguishability of the sheets.

## **8.4 Complex small part and co-molding**

A successor to the double curved dome in Section 8.1 has been designed in cooperation with Lucas Bretz and manufactured at KIT wbk, Karlsruhe, Germany. The part is used to validate the combined application of all DBS model features (non-isothermal, non-Newtonian, compressible, co-molding with patches) in a single part that contains beads and ribs.

### **8.4.1 Experimental setup**

The part and initial stack position are shown in Figure 8.29. Two beads with different flank angles are positioned in the part and a rib is placed centrally in each bead. Both ribs are 2 mm wide at the top and have a draft angle of 3°. The rib height varies from 15 mm at one end to 0 mm at the other end. A lever is integrated into the part for easier recovery from the mold. The SMC stack has the initial dimensions 100 mm x 75 mm x 7 mm and is placed centrally into the mold, which is heated to 145 °C. Optionally, a single ply patch with

dimensions 80 mm x 7 mm x 0.45 mm can be placed in the narrow bead and a patch with dimensions 80 mm x 15 mm x 0.45 mm in the wide bead. The press profile is prescribed as a table with a maximum compression force of 130 kN and parts are molded on a Lauffer hydraulic press in collaboration with Lucas Bretz at KIT wbk, Karlsruhe, Germany.

**Figure 8.29:** Complex small part geometry comprising two beads with different flank angles and centrally placed ribs with decreasing height. The molded part has outer dimensions 120 mm x 94 mm. Patches are optionally placed inside the beads.

The molded samples were analyzed by volumetric imaging using the same setup as described in Section 8.1, but with the resolutions 31 µm and 50 µm for the high resolution and low resolution scans, respectively. The µCT scans and volumetric image reconstruction were performed by Ludwig Schöttl at KIT IAM-WK, Karlsruhe, Germany.

## **8.4.2 Numerical setup**

The process is modeled with the Direct Bundle Simulation technique explained in Chapter 5. The upper and lower mold surfaces are extracted from the CAD model of the part, placed at an initial distance of 22.4 mm and are meshed with rigid body elements of 1 mm average edge length. The lower

mold remains fixed, while the upper mold is controlled by a virtual press controller (see 4.1.5) following the press profile listed in Table A.9. The Eulerian domain is expanded upward and downward beyond the mold surfaces and meshed with 0.8 mm x 0.8 mm x 0.8 mm EC3D8RT. The in-plane boundaries are implemented by constraining the velocity in normal direction to zero. The regions of initial stacks are filled with material in the beginning and 18000 bundles of 25 mm are generated in the stack to represent the desired fiber volume fraction. The paste material is compressible, non-isothermal and non-Newtonian with parameterization from Section 7.1. A hydrodynamic friction model with parameters from Section 8.3 is employed at the mold walls. All simulation parameters are summarized in Table A.9.

The basic configuration does not include unidirectional reinforcement patches. But patches can be added via 3.5 mm x 4 mm and 3.75 mm x 4 mm S4RT shell elements following the model described in Chapter 6 and with the parameterization in Section 7.2. The interaction between the patch surfaces and the SMC paste is modeled via no-slip conditions, as both utilize the same sticky matrix material. The interaction between the molds and patches is modeled with a Coulomb contact with a friction coefficient *µ =* 0.3 according to the measurements by Bücheler [169].

## **8.4.3 Results and discussion**

### **8.4.3.1 Flow front progression**

The deformation of the SMC stack is shown in Figure 8.30 for the simulation and short-shots at corresponding thicknesses. First, the stack is formed to the shape of the mold without a flow process. As soon as the mold gap decreases below the thickness of the formed stack, a flow process starts and fills the ribs. It is noticeable that fiber bundles are pressed against the stamps, because both ends of fibers oriented in *x* direction are pulled in different directions, while the stamps press them down. Both ribs are filled equally from the bottom and the side connected to the bead wall. Hence, the last

filled region is not equivalent to the maximum rib height, but the volume farthest away from any entry to a rib. This behavior is coherent between simulation and experimental short-shots. As the mold is completely filled, a fraction of the simulated paste penetrates the rigid body walls and floats through the domain as droplets.

**Figure 8.30:** Snapshots of the compression process with two beads and ribs. The left column shows the simulation results and the right column shows short shots of the molded parts. The complex geometry leads to leakage of the paste through the rigid mold surfaces due to inaccurate contact handling between the Eulerian surface and Lagrangian surface in Simulia Abaqus/Explicit.

### **8.4.3.2 Fiber-matrix separation**

Fiber-matrix separation can be observed with bare eyes in the molded parts due to the opaque resin. A exemplary part is scanned via µCT and a section across the rib in the wider bead is shown in Figure 8.31. Figure 8.31a shows a volumetric image reconstruction of the fiber architecture in the wider bead. The fiber bundles displayed here stay intact even after entering a rib. The slice of a lower resolution scan shown in Figure 8.31c reveals that a resinrich region forms in the region that is filled last according to the simulated flow front progression. A higher resolution slice of this region is shown in Figure 8.31b. The mesoscale simulation approach predicts a similar pocket of resin in the same region, even though the actual detailed realization of the fiber architecture depends on the randomness of the initial charge in the simulation as well as the experiment.

**(a)** 3D µCT view **(b)** High resolution µCT slice

**(c)** Low resolution µCT slice.

**(d)** Direct Bundle Simulation

**Figure 8.31:** Fiber-matrix separation in the ribs. The µCT images were obtained by Ludwig Schöttl at KIT IAM-WK, Karlsruhe, Germany.

### **8.4.3.3 Patch damage**

The quality of parts molded with single patch layers varies significantly between parts with intact patches and parts with defects such as fracture, displacement and deformation of patches. This is shown in Figure 8.32 and the variance is in line with the large scatter observed in the characterization experiments in Section 7.2.

**Figure 8.32:** Patch defects observed in molded parts.

**(a)** Temperature **(c)** Matrix damage **Figure 8.33:** Simulation results of the section points adjacent to the mold surface.

The result of a simulation with patches and fiber bundles is shown in Figure 8.33. The flow front progress and resulting fiber architecture is unaffected by the presence of patches, except for a slightly faster filling process due to the volume occupied by patches in the cavity. The maximum temperature of patches adjacent to the mold surface is approximately 75 °C and therefore within the temperature range investigated during characterization. The patch model does predict damage and deformation in the patches, but the deformation mode does not strictly match the experimental results. Some paste

material reaches between patch and mold surface, because leakage at convex mold surfaces allows material to get to the other side of the patch.

## **8.4.4 Concluding remarks**

The application example demonstrates the ability of the developed models to describe non isothermal compression molding and co-molding at component scale. The mold filling process and formation of resin-rich areas is predicted qualitatively correct and the µCT images illustrate that fiber bundles stay intact even after entering a rib. The partial penetration of Eulerian material through mold surfaces is a well known problem in compression molding simulations with the Coupled Lagrangian Eulerian framework of Simulia Abaqus/Explicit [158, 210]. It can occur at convex curved mold surfaces due to the employed surface reconstruction method [194] of the internal solver. Large pressures and low viscosity foster the formation of droplets that penetrate the wall. Even though a finer mesh and smaller time step can reduce the effect, it always remains a disadvantage of the employed surface reconstruction technique. The integration of patches in the mesoscale compression molding process simulation is possible and the occurrence of damage in patches is predicted correctly. However, the accurate prediction of patch defects is difficult due to the inhomogeneous patch quality, which results in different experimental outcomes and uncertain model parameters.

## **8.5 Demonstrator compression molding**

The demonstrator component of IRTG GRK 2078 serves as a use case that validates the applicability of the developed mesoscale DBS at larger scale components. The component is designed collaboratively within IRTG GRK 2078 to incorporate three beads that reinforce the part under bending load.

## **8.5.1 Experimental setup**

The mold consists of a base mold with dimensions 800 mm x 250 mm and an insert geometry with dimensions 230 mm x 230 mm. This insert geometry is visualized in Figure 8.34a. The upper insert can be replaced with a ribbed structure that has varying rib features in each bead (see Figure 8.34b). The outer beads have six crossing rib pairs each, with a thickness of 3 mm and 2 mm, respectively. The radius at the rib base varies between 0.5 mm and 1.5 mm in these outer beads. The central bead has a higher rib density with a total of eight crossing rib pairs, but constant radii of 1.5 mm and a constant thickness of 2 mm.

**<sup>(</sup>b)** Insert with ribs and full mold coverage **Figure 8.34:** Molds and initial stack configurations of the IRTG demonstrator component.

The components are manufactured from UPPH-GF by Sergej Ilinzeer at Fraunhofer ICT, Pfinztal, Germany on a Dieffenbacher COMPRESS PLUS DCP-G 3600/3200 AS and a Dieffenbacher DYL630/500 hydraulic press with parallel

cylinder control. Compression force and cylinder displacements are recorded during these manufacturing molding trials. Parts without ribs are molded from a 320 mm x 240 mm x 7 mm initial stack that is placed centrally over the insert in the mold and parts with ribs are manufactured with full mold coverage of two SMC sheets.

The demonstrator component's size prohibits an efficient evaluation of the fiber orientation via µCT for the entire component. However, the opacity of the employed resin system enables fluoroscopy with visible light. Hence, a fluoroscopy workflow was developed in collaboration with Sven Revfi from KIT IPEK, Karlsruhe, Germany and the workflow is depicted in Figure 8.35.

**Figure 8.35:** Image acquisition and image processing workflow. The specimens are placed on an overhead projector and images of the shine through region are taken. The images are merged and the lighting is homogenized in an image processing software to obtain a gray valued image of the entire specimen with 12.5 pixels mm<sup>−</sup>1. The gray valued image is then analyzed with orientation analysis software. (The procedure was developed and applied together with Sven Revfi from KIT IPEK, Karlsruhe, Germany and Ludwig Schöttl from KIT IAM-WK, Karlsruhe, Germany.)

Specimens are placed on an overhead projector and partial images are acquired with a Canon EOS 70D DSLR due to size restrictions and to improve the image resolution. These partial images are digitally processed to correct inhomogeneous lighting, distortions and improve contrast. Subsequently

they are aligned and merged to a single gray value image of the entire component. The procedure is repeated for the top and bottom side of three parts in total. Each image is analyzed with OrientationJ [208] utilizing a smoothing length of 1 pixel and a grid size of 5 pixels at which orientation vectors are determined. In addition, the orientation on image data is analyzed using structure tensors [211]. The fluoroscopy procedure is corroborated by µCT with 42 mm resolution and 56 µm resolution. The µCT scans, volumetric image reconstruction and fiber orientation analysis by structure tensors were performed by Ludwig Schöttl at KIT IAM-WK, Karlsruhe, Germany.

### **8.5.2 Numerical setup**

The production process is modeled with the Direct Bundle Simulation technique explained in Chapter 5. The upper and lower mold surfaces are extracted from the CAD model of the part, placed at an initial distance of 25.2 mm and are meshed with rigid body elements of 2.5 mm average edge length in the planar region and 1.0 mm average edge length in the bead region. The lower mold remains fixed, while the upper mold is controlled by a virtual press controller (see 4.1.5) following the press profile listed in Table A.10. The Eulerian domain encompasses the entire region between the molds and is meshed with 2.5 mm x 2.5 mm x 1.0 mm EC3D8RT elements in the planar regions and 1.0 mm x 1.0 mm x 1.0 mm EC3D8RT elements in the region of the beads. The normal velocity of the boundaries of the Eulerian domain is constrained to zero and the region of the initial stack is filled with material in the beginning. For the unribbed components, a total amount of 150000 bundles with *L =* 25mm length with element size *l =* 2.5mm are generated in the stack to represent the desired fiber volume fraction. As the ribbed part is molded from a full initial mold coverage and the plate region has only little effect on the rib filling, the simulated domain is reduced to the insert area for this case. This also reduces the number of simulated fiber bundles to 30000, which are resolved finer with *l =* 1mm. The paste material is compressible,

non-isothermal and non-Newtonian following the parameterization of Section 7.1. A hydrodynamic friction model with parameters from Section 8.3 is employed at the mold walls. All simulation parameters are summarized in Table A.10.

### **8.5.3 Results and discussion**

### **8.5.3.1 Flow front progression**

Snapshots of the filling process are depicted in Figure 8.36. Initially, the stamps push the stack into the beads and fiber bundles are pulled into the shape leading to a higher orientation perpendicular to the beads. The outer stack contour is affected by this pull-in and the stack experiences out-of-plane deformation. At 4.8 s, the beads are fully filled and the flow front extends uniformly in both directions. However, leakage occurs at some convex mold surface areas similar to the previous example in Section 8.4. Due to the asymmetric placement on top of the insert, the flow front reaches one end of the mold first (7.2 s) and the flow direction reverts in some regions to completely fill the mold.

Flow marks with highly aligned fibers form at the ends of the beads during molding (see Figure 8.37a, where flow marks are highlighted with a light green color). For an insert-centered initial stack, these flow marks are symmetric (equal length at both ends of the beads). The mesoscale DBS is able to predict the formation of these flow marks as a result of bundles with restricted motion due to the beads. As shown in Figure 8.37c, fiber bundles at the end of the beads travel a shorter distance than those surrounding them, which leads to an alignment of fibers. An adjustment of the initial stack to the mold center results in longer flow marks that are restricted to one side of the part (see Figure 8.37b, where flow marks are highlighted with a light green color). The effect of this small shift by 6 cm can be reproduced by DBS, where similar long flow marks are predicted (see Figure 8.37d). The fiber architecture shows artifacts in the region of the beads, which are caused by leakage of paste

**Figure 8.36:** Snapshots of the simulated mold filling process for the IRTG demonstrator component without ribs and an insert-centered stack. The left column shows the top view on the component and the right column shows a section through the center bead. The last snapshot shows some droplets leaving the cavity due to leakage.

droplets through the mold surfaces. These artifacts can occur, if a droplet leaves the domain rapidly and affects the neighborhood of a fiber bundle, with which it interacts.

**Figure 8.37:** Formation of flow marks in the demonstrator part without ribs. The displacement describes the distance that each bundle traveled relative to its initial position. The insertcentered stack results in symmetric flow marks at all ends of beads. The mold-centered stack results in long one-sided flow marks, which are observed in experiments and the simulation.

The ribbed demonstrator component molded with full mold coverage shows severe fiber-matrix separation in the ribs with resin-rich regions at the crossing points of each rib pair. The separation can be observed visually, as depicted in Figure 8.38a, and seems to occur independent of rib width and radii at the rib entrance. The ribs are completely filled with paste in the mesoscale simulation, but fiber bundles are also segregated in the ribs, as shown in Figure 8.38b. The segregation is caused by entangled bundles as well as bundles that have a limited capability to enter multiple ribs due to the inextensibility. Additionally, the drag force on bundle ends that just reach into a rib is not sufficient to move the entire bundle into the rib and increases separation.

The central bead of a exemplary ribbed part is scanned by µCT and depicted in Figure 8.39a. The contour surface separating fiber bundles and matrix agrees with the visually observed separation. The DBS also shows a clear lack of bundles in the crossing points of ribs, as depicted in Figure 8.39b.

**(a)** Photo **(b)** DBS **Figure 8.38:** Overview on fiber-matrix separation in the insert region of the demonstrator component. (a) The dark translucent shade at the crossing points of ribs can be clearly attributed to fiber-matrix separation. (b) The crossing points of ribs are lacking fiber bundles. The matrix is not shown for visualization purposes, but the simulation predicts a complete filling of the ribs with matrix material and hence a fiber-matrix separation similar to the experiments.

**(b)** DBS

**Figure 8.39:** Overview on the fiber architecture of the central bead obtained from µCT and DBS.

A more detailed view on the rib pair highlighted by a circle in Figure 8.39a is displayed in Figure 8.40 and compared to the corresponding simulation result. The flow front progression towards the center of the bead and resulting fiber architecture is predicted correctly. This is remarkable, as a recent study on very similar structures with commercial macroscopic models falsely predicted a preferred filling of the wider central crossing points of ribs [212]. The

authors assume that this stems from a lack of anisotropy in their models, but the results presented here suggest that the mesoscale structure plays an important role in the fill pattern of such ribs.

**(a)** µCT scan **(b)** DBS **Figure 8.40:** Zoom into the fiber architecture in a central rib pair. The section is highlighted by a circle in Figure 8.38

### **8.5.3.2 Fiber orientation**

The *A*<sup>11</sup> component of the fiber orientation tensor is depicted in Figure 8.41 for the analysis with OrientationJ on fluoroscopy images, structure tensors on fluoroscopy images, structure tensors on µCT scanned sections and the DBS result. The image data is evaluated on an exemplary sample and the DBS result is mapped on a 2.5 mm x 2.5 mm x 2 mm evaluation mesh with the method described in Section 5.1.6.3.

Both results based on image data show a predominant orientation in *x*direction with a major perpendicular orientation in the beads and at the ends of the mold. This perpendicular orientation is more prominent in the structure tensor based computation. It is likely that this structure tensor based computation is more accurate due to better image pre-processing. Additionally, the structure tensor based computation leaves out the region with the sample label to avoid misinterpretation of the orientation in this region. Compared to image data, the evaluation of regions analyzed via µCT show a higher orientation in terms of absolute values, but the same trend of

**Figure 8.41:** Orientation tensor component *A*11 for an exemplary part analyzed with all described methods. The computation in (b) leaves out the region with the sample label to avoid misinterpretation of the orientation in this region.

perpendicular orientation in the beads. As the µCT scan contains information through the thickness of the sample, this suggests that the fluoroscopy method underestimates strong alignment due to its thickness averaging property. The DBS features distinct highly oriented regions at the end of each bead that correspond to the flow marks described in the previous section. These are less pronounced in the experimental data, even though strong fiber alignment is visually observed in the molded parts. The mesoscale simulation does predict a strong perpendicular orientation in the beads and even between beads. This is reasonable, as fibers are pulled into the beads such that the initial stack is stretched ª140% in *y*-direction at the central area. Pulling fibers into the beads affects also the region between beads due to the fiber length of 25 mm. However, the effect is less pronounced in the experimental data and only visible in the central section of the µCT scan. More µCT scans of the other beads would be necessary to investigate the fiber architecture of this region in more detail.

#### **8.5.3.3 Compression force**

The compression force of the demonstrator part is plotted in Figure 8.42. The plot shows the mean force of five trials in dark gray as well as the area between the lower quartile and upper quartile in light gray. The compression force features a plateau just after the initial forming phase, in which the stamps deform the stack to the beads. The most likely reason for this plateau is the initial sticking friction at the mold surface. The DBS simulation features a similar plateau, but it is less pronounced. This is likely affected by the poor contact enforcement in Simulia Abaqus/Explicit between the Eulerian phase and complex Lagrangian molds. Accordingly, the final thickness of the part is also underestimated due to the leakage at convex mold features.

**Figure 8.42:** Compression force and mold displacement during the demonstrator compression molding process. Dots mark the snapshot positions depicted in Figure 8.36.

### **8.5.4 Concluding remarks**

The mesoscale Direct Bundle Simulation approach is applied to a large-scale component of dimensions 800 mm x 250 mm with ribs and beads. The model is able to predict molding defects such as the formation of flow marks and fiber-matrix separation in ribs. In contrast to recent reports on macroscopic models in similar structures [212], the rib filling process is described qualitatively correct from the outer borders towards the center and not with a preferred flow in the wider center of crossing rib pairs. The comparison of fiber orientation with experimental results shows a correct prediction of the predominant fiber orientation in *x*-direction of the plate with perpendicular orientation in the beads and at the end of the mold. However, the coarse mesh and previously described penetration of droplets through the mold walls cause artifacts in the bundle architecture and an underestimation of the final part thickness. The mapping procedure of the orientation evaluation can be also used to transfer the process simulation results to a mesh for subsequent structural simulation. Contrary to regular macroscopic models applied in the state of research, the new mesoscale model allows to transfer fiber volume fraction and evaluate higher order fiber orientation tensors, if desired. A comparison to results obtained with the commercial software Autodesk Moldflow may be found in the thesis of Revfi [213].

## **9 Summary**

After a summary of the state of research in SMC compression molding and stating the objective, this work proposes a macroscale reference model for SMC process simulation. The model is supplemented by utility functions such as a virtual press controller and a friction model, which are also used for the mesoscale simulations. The main contribution is the development of a mesoscale Direct Bundle Simulation (DBS) approach that models SMC paste as a combination of paste moving in an Eulerian frame and Lagrangian truss elements representing individual fiber bundles. The model includes bundlebundle interactions as well as bundle-matrix interactions. A collection of pre- and post-processing tools for productive application at component scale completes the mesoscale model. Additionally, an anisotropic patch model with individual damage variables for stitching yarn and infiltrated UD carbon fibers is developed. The required parameters of the models are characterized and applied to several validation examples.

## **9.1 Conclusions**

The prediction of fiber architecture via the novel mesoscale Direct Bundle Simulation is improved in comparison to macroscale models in confined regions and close to the mold walls. The method is able to predict the formation of knit lines and flow marks in the same positions as experimentally validated. All this can be applied for entire components and is to the authors best knowledge the first method resolving individual bundles enabling component simulations.

For large planar SMC regions, Jeffery's model without interaction parameters seems to be a good choice for a macro model, as it agrees with the more detailed DBS in such large planar regions. This is remarkable, because it does not require additional interactions parameters and retarding principal rate parameters as they have been developed for injection molding. The macroscopic model features up to ten times faster computation times than the mesoscale model, but it is not able to predict fiber-matrix separation and bundle curvatures.

The proposed mesoscale method allows the prediction of fiber-matrix separation and bundle curvature, because flexible bundles and paste are able to move at different speeds. The predicted locations of resin-rich regions agree with experimental results.

Lubrication forces between individual fiber bundles are incorporated into the model, but are numerically challenging due to the resulting small time step during explicit integration. The introduction of lubrication affects the compression force during molding of planar SMC parts, while the orientation is hardly affected. In most cases, mold friction is a dominant contributor to compression force in comparison to lubrication and lubrication may be neglected in favor of faster computation times.

The ideal plug-flow assumption has limited validity for the considered UPPH-GF SMC and sheets do slip relative to each other, which has been proven by a colored sheet. Even with mesoscale DBS and accounting for mold friction, press compliance, compressibility, non-Newtonian non-isothermal paste, it is not possible to describe this layer-induced mechanism with a single matrix phase. Likely, it would require separate matrix phases, which arises new challenges w.r.t mesh resolution and the ability to merge such phases.

To model the squish and flow phase of SMC, sticking friction should be added to the hydrodynamic friction model. Therefore, a friction model is suggested that transitions from an initial sticking behavior to a hydrodynamic powerlaw friction model. However, the friction model relies on the implementation of the surface reconstruction from the Eulerian phase. The employed Simulia

Abaqus/Explicit CEL method shows some weaknesses in this area that lead to an imperfect contact and leakage of paste material through the mold surfaces.

Patches from the in-house production suffer from significant scatter due to a missing series production quality control. They can be integrated to the simulation and the occurrence of damage is predicted correctly, but the exact deformation mode cannot be predicted reliably due to the scatter in patch properties.

## **9.2 Outlook and recommendations**

For macroscopic SMC orientation models, the mesoscale simulation results suggest that a focus on further development of refined fiber orientation parameters is not necessary for large planar SMC regions. Instead, it seems more promising to further develop models that account for confinement in narrow regions following for example the work of Perez et al. [25].

The mesoscale process model can enrich the previously developed CAE chain for SMC [6] with higher-order fiber orientation tensors, fiber volume fraction, bundle curvature and other information given by the more detailed prediction of fiber bundle architecture. Given the probabilistic nature of some molding mechanisms, the scatter should be assessed with a continuous probabilistic CAE chain that leverages the ability of the developed model to predict different realizations for different initial randomly generated stacks under identical processing conditions.

The mesoscale model itself could be improved by extensions to model flattening of bundles and disintegration of bundles that can occur in rare cases. For example, a damage variable could describe the disintegration of the bundle character based on the subjected shear stress on bundles. Finally, the use of an improved implementation could address the leakage problem and generally improve the computation performance in Simulia Abaqus/Explicit CEL.

## **A Appendix**

## **A.1 Sketches for geometric considerations**

Equation (5.20) was obtained by a projection, which is illustrated in Figure A.1a. The normal direction is computed from the distance between the tip of *p<sup>j</sup>* and the projection of *p<sup>j</sup>* on the unit vector in direction of the velocity [[¢*vj*]].

The curvature given in Equation (7.4) is the inverse radius of the bent section. The radius *R*<sup>b</sup> comprises a lower section with length *r*<sup>0</sup> sin*Æ*<sup>b</sup> and an upper section with length (*r*<sup>0</sup> *+r*<sup>0</sup> cos(*Æ*b))cot*Æ*.

**(a)** Normal direction **(b)** Bending **Figure A.1:** Sketches for the derivation of Equation (5.20) and Equation (7.4)

## **A.2 Simulation parameters**


**Table A.1:** Parameter variation of the reference model presented in Figure 4.6 in Section 4.3.


**Table A.2:** Simulation parameters for anisotropic flow verification in Section 5.2.2.

**Table A.3:** Simulation parameters for linear compression example with 0% fiber volume fraction.



**Table A.4:** Simulation parameters for linear compression example with 1% fiber volume fraction.


**Table A.5:** Simulation parameters for linear compression example with 25% fiber volume fraction.


**Table A.6:** Simulation parameters for simple double curved dome.


**Table A.7:** Simulation parameters for the honeycomb part.


**Table A.8:** Simulation parameters for the press rheometer.


**Table A.9:** Simulation parameters for the complex small part.


**Table A.10:** Simulation parameters for the IRTG demonstrator part.

## **List of Figures**





A.1 Sketches for the derivation of Equation (5.20) and Equation (7.4) . 215

## **List of Tables**



## **Bibliography**


microstructure of SMC composites," *Composites Part A: Applied Science and Manufacturing*, vol. 39, no. 1, pp. 91–103, 1 2008.


uncured sheet moulding compounds: Rheology and flow-induced microstructures," *Composites Part A: Applied Science and Manufacturing*, vol. 101, no. July, pp. 459–470, 10 2017.


compound (SMC) processes," *International Journal of Material Forming*, pp. 1–11, 8 2019.


tension coupling," *Textile Research Journal*, vol. 89, no. 3, pp. 434–444, 2 2019.


## **List of own Publications**

## **Journal articles**


## **Book contributions**

• Johannes Görthofer, Nils Meyer, Ludwig Schöttl, Anna Trauth, Malte Schemmann, Pascal Pinter, Benedikt Fengler, Sergej Ilinzeer, Martin Hohberg, Tarkes Dora Pallicity, Luise Kärger, Kay A. Weidenmann, Peter Elsner, Frank Henning, Andrew Hrymak, and Thomas Böhlke. Compression Molding of the Demonstrator Structure. In *Continuous–Discontinuous Fiber-Reinforced Polymers*. Carl Hanser Verlag GmbH & Co. KG, Munich, Germany, 8 2019.

## **Conference proceedings**


## **Presentations**


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

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




## Karlsruher Schriftenreihe Fahrzeugsystemtechnik

Sheet Molding Compounds (SMC) are discontinuous fiber reinforced composites that are widely applied due to their ability to realize composite parts with long fibers at low cost. They enable function integration and forming of complex components, but require a profound understanding of the processing mechanisms and their influence on the performance of a component. Process simulations address this point by predicting possible manufacturing defects and process parameters. These results can not only be used to configure the process and reduce trial and error phases, but also for subsequent structural simulation by virtue of a virtual process chain.

The hypothesis in this work is that a direct simulation of individual fiber bundles is required to accurately describe the SMC mold filling process of complex geometries. Based on this hypothesis, a novel Direct Bundle Simulation (DBS) method is proposed to enable a direct simulation at component scale utilizing the observation that fiber bundles often remain in a bundled configuration during SMC compression molding. The developed DBS model can be combined with patches to simulate the co-molding process of SMC with continuous fiber reinforcements.

ISSN 1869-6058 ISBN 978-3-7315-1173-1