# Multiphysics for simulation of forming processes

Luisa Silva, Patrice Laure, Thierry Coupez, Hugues Digonnet

## ► To cite this version:

Luisa Silva, Patrice Laure, Thierry Coupez, Hugues Digonnet. Multiphysics for simulation of forming processes. Nicolas BOYARD. Heat Transfer in Polymer Composite Materials: Forming Processes, Wiley-ISTE, 2016, 978-1-84821-761-4. hal-01284223

HAL Id: hal-01284223

https://hal.science/hal-01284223

Submitted on 7 Mar 2016

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

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

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Multiphysics for simulation of forming processes

Luisa Silva¹, Patrice Laure², Thierry Coupez¹, Hugues Digonnet¹

¹ Institut de Calcul Intensif, Ecole Centrale de Nantes, France

² Laboratoire J.A. Dieudonné, Université de Nice Sophia-Antipolis, France

e-mail: luisa.rocha-da-silva@ec-nantes.fr

**Abstract**

This chapter is dedicated to the simulation of composite forming processes. It concerns numerical methods for molding, at the different material scales, involving multiphysics analysis. From the different possible numerical techniques, we focus on a stabilized finite element method for the resolution of the continuum equations (mass, momentum and energy conservation), in an Eulerian framework, with a level-set approach for interface displacement treatment. Internal variables, such as the reticulation rate or fiber fraction and orientation may be computed using models defined through ordinary differential equations. Examples at the process scale, but also at the composite inner structure scales illustrate the capabilities of such numerical techniques.

**Keywords**

Composite forming simulation, Multiphysics modeling, Finite element method.
# 1 Introduction

In recent years, the use of composite materials has been widely diffused and continues to evolve in the aeronautics, the aerospace and the automotive industries. Most of the composites used are composed of an organic matrix and fiber reinforcements, ensuring a certain level of physical properties, designed taking into account the material's anisotropy, with a large weight reduction when compared to traditionally used metallic alloys. Thermoplastic or thermoset may be used as resins, whereas one may have glass, carbon, or even natural-based fibers.

To produce composite parts, processes may be more generally grouped into three categories [1, 2]. The first one concerns short fiber reinforced matrices, involving fiber transport with resin flow, as a non-Newtonian suspension being injected into a closed mold or through a die. The second one involves highly viscous thermoplastics or thermosets resins containing continuous or long discontinuous fibers, as well as high fiber content rate. In this case, the reinforcement and the resin are heated and deform together under an applied (compression) stress to form the required composite shape like, for example, in compression molding of Sheet Molding Compounds. Finally, the third category concerns porous media processes, when the reinforcement is composed of continuous fiber networks or textiles placed in a (closed or open) mold into which a low viscosity resin is injected or placed, heated and compressed, in such a way that it impregnates the fibrous architecture.

Complex multiphysics phenomena, strongly coupled, arise in these different processes: short fibers orient, interact and may aggregate; preforms are compressed and may locally be subject to large deformations; the fluid flows in a medium that has very high anisotropic features, undergoing simultaneously thermal and chemical solicitations. Control of the process, in all the cases, remains difficult, and inner properties like local fiber fraction and orientation, or residual stresses, strongly vary as a function of process conditions. Numerical techniques may bring answers when comprehension of what happens during the process is difficult, may provide quantitative values for the influence of the different process parameters, allowing its improvement, enhancing its productivity and reducing the development time.

# 2 Multiscale, multiphysics and multidomain modeling

Modeling of composite forming (and in particular when the resin is at the liquid state) is thus very complex and may include simulating the preforming and the flow processing stages. Preforming includes draping, braiding or fiber placement and has also been object of the development of specific numerical applications, like in [3, 4], not illustrated in this chapter. In the last decades, there has been also a huge progress on simulation of the flow step, with different codes simulating the process, at different scales.

To develop numerical models, one needs to first decide what is the scale of description of the physical mechanisms to study: microscopic, mesoscopic or macroscopic
(Fig.1). At the microscopic scale, one distinguishes each fiber and resin flow is computed. It may induce fiber motions, important, for example, in the comprehension of the rheological behavior of suspensions. For continuous fibers assembled in tows, the mesoscopic scale includes not only the study of the resin flow between the tows, but also impregnating them. Most often, the tow is considered as an homogeneous porous media. This dual scale structure is important, since one needs to ensure that the resin fully saturates the inter-tow spaces, but also intra-tow ones, to prevent the appearance of macro or micro voids. At the macroscopic scale, an homogeneous equivalent media being injected or compressed is considered, and simulation occurs for the forming of the whole part. Even if this last is the one of interest for industrial and process development purposes, material parameters and models may be very complex, with data difficult to obtain experimentally. Then, it is useful to couple simulation at all scales: at lower levels, one gets closer to what experiences the material and less inner or homogenized properties are necessary.

<figure><img src="image_2.png" /><figcaption>Figure 1: *Spatial scales for composites: (a) fiber suspensions, with a macroscopic or part scale and a microscopic one, composed of the resin and chopped fibers (b) continuously reinforced composites, with a macroscopic scale, a mesoscopic (tow) scale, and a microscopic (fiber) scale.*</figcaption></figure>

At all scales, the same general continuum mechanics approach is used and detailed in the following. Mass, momentum and conservation equations are solved and constitutive laws and internal variables evolution depend on the scale and on the observed phenomena. Flow equations and classical thermo-rheological-kinetical coupling are illustrated in the next section.
## 2.1 Flow equations

Let us consider the flow of a fluid around a bundle of rigid fibers (chopped or continuous), thus placing ourselves at the microscopic scale. The injected resin may be considered as a Newtonian or non-Newtonian incompressible fluid and the following Navier-Stokes equations, coupled to an appropriate set of boundary conditions, describe the fluid flow:

$$ 
\left\{ 
\begin{aligned}
\rho_f \left( \frac{\partial \vec{v}}{\partial t} + \vec{v} \cdot \nabla \vec{v} \right) - \nabla \cdot \underline{\sigma} &= \rho \vec{g} \\
\nabla \cdot \vec{v} &= 0 \\
\underline{\sigma} &= 2 \eta_f \underline{\varepsilon}(\vec{v}) - p \underline{I} \\
\underline{\varepsilon}(\vec{v}) &= \frac{1}{2} [\nabla \vec{v} + t \nabla \vec{v}]
\end{aligned}
\right. 
\quad (1) 
$$

with $\vec{v}$ the resin velocity, $p$ the pressure, $\rho$ the resin density, $\eta_f$ the dynamic viscosity of the resin, $\underline{\sigma}$ the stress tensor, $\underline{\varepsilon}(\vec{v})$ the strain rate tensor and $\underline{I}$ the identity tensor. In most cases, injection pressures and velocities imply that inertia and gravity terms may be neglected, leading to the Stokes equations. This set of equations needs then to be coupled to a solid behavior model to take into account fiber motion or displacement. Fig.2 shows an example of flow computation of a Newtonian resin between a network of rigid static fibers, to illustrate the type of attended results, at the fiber scale.

<figure><img src="image_3.png" /><figcaption>Figure 2: Illustration of resin flow computation across an array of unidirectional fibers: (a) flow path and (b) velocity norm distribution [5].</figcaption></figure>

To establish conservation equations at the mesoscopic or macroscopic scales, averages or semi-analytic homogenization may be used, leading to two general cases.

In the first one, fibers are chopped, remain rather short and its concentration when impregnated by the fluid is dilute or semi-dilute. Then, the overall media being injected in a part (at the macroscopic scale) can be considered has an homogeneous fluid, with
a stress tensor that may be decomposed in a fluid *f* and a fiber *s* (anisotropic) contributions $\underline{\sigma} = \underline{\sigma}_f + \underline{\sigma}_s$, where $\underline{\sigma}_f$ is defined as in Eq.(1) This last means that one needs to follow, during part forming, the local evolution of fiber orientation and fraction. More details will be given in throughout an illustration in Sec.3.4.1.

In the second one, fibers are gathered in yarns composing a textile reinforcement or chopped but assembled in a dense mat, and the pore scale is very small when compared with the flow length scale. In this case, Darcy's law [6] is the most used to model flow through this porous media, and is established by performing a volume average of the Navier-Stokes equations describing the flow at the microscopic scale, supposing that the solid skeleton composing the reinforcement is static and non-deformable, and that the porous media is saturated. It is given by

$$ \eta_f \underline{K}^{-1} \vec{q} + \nabla p = 0 \qquad (2) $$

Here, $\vec{q}$ is the discharge rate, such that $\vec{q} = \phi \vec{v}$, where $\phi$ is the fluid volume fraction, or porosity. The permeability, $\underline{K}$, is a tensor for an anisotropic porous media. In the case of unidirectional fibers packing, if the *z* axis of the coordinate system is taken in the same direction as the fiber axis, the permeability tensor can be then written as follows:

$$ \underline{K} = \begin{pmatrix} K_{\perp 1} & 0 & 0 \\ 0 & K_{\perp 2} & 0 \\ 0 & 0 & K_{\parallel} \end{pmatrix} \qquad (3) $$

with $K_{\perp 1}$ and $K_{\perp 2}$ the transverse permeability of the packing and $K_{\parallel}$ the permeability along the fiber axis.

Several articles use the Brinkman's equations [7] instead of Darcy's law to describe this type of flow, by generalizing Navier-Stokes equations supposing the flow of a Newtonian fluid through a swarm of fixed particles and writing:

$$ \eta_f \underline{K}^{-1} \vec{q} - \nabla \cdot [2\eta_e \underline{\varepsilon}(\vec{v})] + \nabla p = 0 \qquad (4) $$

where $\eta_e$ is the effective viscosity, which may be different from $\eta_f$. One advantage of Eq. 4 is that it is an extension of the Stokes equation. Thus, numerical tools developed for this latter may be used in the Brinkman's case with a small amount of implementations. However, its physical domain of validity is very restricted (small fiber concentrations, when $\phi \to 1$, corresponding to very large permeabilities), and theoretical analysis shows that it yields to $\eta_e = \eta_f$ [8]. Nevertheless, the Brinkman form will be preserved in the following, since it allows to treat numerically both Stokes and Darcy cases, by considering extreme values of $\underline{K}$ and $\eta_e$. Sec.3.4.2 shows an example of flow on one of these cases.

## 2.2 Thermal-rheological-kinetical coupling

Liquid forming processes involve more complex physical phenomena than "just" solving a fluid (and solid) mechanics problem. In fact, modeling these processes involves
studying the coupling between this problem, the thermal one and, for some resins, like thermosets, the polymerization or degree of curing reaction. In fact, resin cure affects its viscosity, which modifies the flow pattern. Higher viscosities imply also more pressure and higher deformation and motion of the fiber reinforcement, that changes the permeability, influencing the resin flow. Temperature affects directly the rheology, changing the viscosity and modifying the reinforcement geometry by thermal dilatation. In addition, polymerization impacts thermal transfer, since curing is exothermal, producing heat, which also accelerates the reaction.

The temperature, governed by the energy conservation equation, depends of the heat transfer by conduction in the resin, fibers and with the mold, the viscous dissipation (which may often be neglected) and on the heat released by the chemical reaction (function of the reaction kinetics). Let us suppose that one wishes to study directly the homogenized composite media. This will allow us to establish a general form of the energy equation). Then, the heat equation in the absence of viscous dissipation, of the convection-diffusion-reaction type, is written as follows:

$$ \rho c_p (\frac{\partial T}{\partial t} + \vec{v} \cdot \nabla T) - \nabla \cdot (\lambda \nabla T) = H_{\alpha} f(\alpha) \quad (5) $$

Here, $\rho$ and $c_p$ are, respectively, the composite density and specific heat, whereas $\lambda$ is the thermal conductivity, which is, for an anisotropic media, also a tensor. All are dependent on the fiber fraction and the latter most often on other inner geometrical parameters, such as orientation. The source term corresponds to the polymerization or cure reaction and is a function of the cure enthalpy, $H_\alpha$, and of the cure rate, $f(\alpha)$. The evolution of the degree of cure $\alpha$ is obtained by solving a simple advection equation with a source term, the cure rate:

$$ \frac{\partial \alpha}{\partial t} + \vec{v} \cdot \nabla \alpha = f(\alpha) \quad (6) $$

Here, $\alpha$ is the degree of cure and $f(\alpha)$ is function of the cure kinetics model chosen. One of the most general and widely used is the Kamal and Sourour model [9], that has the form

$$ f(\alpha) = \left[ k_1 \exp\left(\frac{E_1}{RT}\right) + k_2 \exp\left(\frac{E_2}{RT}\right) \alpha^m \right] (1 - \alpha)^{1-n} \quad (7) $$

and where $k_i$, $m$, $n$ are material parameters, whereas $E_i$ are activation energies. The degree of cure and the temperature affect the viscosity, in opposite senses. More generally, $\eta = \eta_0 g(\vec{v}) \cdot g(T) \cdot g(\alpha)$. Dependence on $\vec{v}$ through the strain rate tensor is taken into account through purely viscous models for classical non-Newtonian resins, and the influence of temperature may be governed by Arrhenius law [10]:

$$ g(T) = \exp\left[\frac{E_a}{R}\left(\frac{1}{T} - \frac{1}{T_{\text{ref}}}\right)\right] \quad (8) $$
where $E_a$ is the activation energy, $R$ the universal gas constant and $T_{ref}$ the reference temperature. Finally, the degree of cure implies a modification of the viscosity according to, for example, a percolation model, which is widely used for epoxy resins [11],

$$g(\alpha) = \left(1 - \frac{\alpha}{\alpha_{gel}}\right)^p \qquad (9)$$

where $\alpha_{gel}$ is the extent of the reaction at the percolation threshold (the gel point) and $p$ is the percolation exponent, both material properties to be identified experimentally. The reference temperature $T_{ref}$ is usually the glass transition temperature $T_g$, influenced by the cure kinetics (when vitrification occurs) and that may be modeled using the DiBenedetto model [12]

$$T_g = T_g^0 + \frac{(T_g^1 - T_g^0)\kappa\alpha}{1 - (1 - \kappa)\alpha} \qquad (10)$$

where $T_g^0$ and $T_g^1$ are glass transition temperature values of the uncured and fully cured resin, and $\kappa$ is an adjustable structure dependent parameter. To illustrate what type of computations use these features, we consider the overmolding of an electrical part with an epoxy resin, Fig.3 [13].

<figure><img src="image_4.png" /><figcaption>Figure 3: Numerical results obtained in overmolding of an electrical part with an epoxy resin: (a) comparison of experimental end of fill positions and numerical predictions for two different filling configurations; (b) temperature (left) and degree of curing (right) at the end of filling (above) and at the end of the post-filling (below) stages [13].</figcaption></figure>
## 2.3 Orientation and structure development during processing

Short or long fiber resin reinforced polymers are mostly applied in injection or compression molding. In particular, injection molding induces an oriented layered structure for this type of composites, where orientation depends on the flow pattern. For example, in the vicinity of the injection gate, the material’s flow front adopts a radial flow extension which is divergent and one gets a transverse orientation compared to the flow direction. Far from the gate, the shear rate is weaker and fibers preserve a state of orientation close to the initial one. Solidification occurring near the wall contributes to the formation of a frozen-in oriented layer. Near this layer, the shear rate is very important and fibers orient in the flow direction. Beyond the ones referred, other parameters influence the orientation state such as the injection speed or the holding pressure [14, 15, 16].

At the macroscopic scale, the most general descriptor of an oriented state is the probability distribution function of the orientation, $\psi(\vec{p}, t)$, that represents the probability to find, at time $t$, a fiber with axis oriented in the direction of the unit vector $\vec{p}$. Computation of $\psi$ can be performed using the Fokker and Plank equation [17], but this quantity is difficult to handle numerically, and thus not commonly used. Other authors [18] introduced the second order $a_2$ (and more) orientation tensors, expressed as a function of $\psi$ and $\vec{p}$, as a quantitative measure of the orientation state, such that:

$$ \underline{a}_2 = \int_{\vec{p}} \psi(\vec{p}, t) \vec{p} \otimes \vec{p} \qquad (11) $$

$$ \underline{a}_4 = \int_{\vec{p}} \psi(\vec{p}, t) \vec{p} \otimes \vec{p} \otimes \vec{p} \otimes \vec{p} \qquad (12) $$

$a_2$ is symmetrical and positive definite, and its trace is equal to 1. For a random orientation, its diagonal terms are equal to 1/3 and, for a random planar orientation, to 1/2. For a perfect orientation in the direction $i$, $a_{ii} = 1$, and the other values are zero. Lipscomb [19] combined Jeffery [20] and Fokker and Plank equations to represent the evolution of the second order orientation tensor by a convection-reaction equation, containing also a fourth order tensor, which can be expressed as a function of the second order one using a closure approximation. Folgar and Tucker [21] extended the previous model to take into account the interaction between fibers. They obtained

$$ \frac{\partial \underline{a}_2}{\partial t} + \vec{v} \cdot \nabla \underline{a}_2 = \omega \underline{a}_2 - \underline{a}_2 \omega + \beta (\underline{\epsilon} \underline{a}_2 + \underline{a}_2 \underline{\epsilon} - 2 \underline{\epsilon} : \underline{a}_4) + 2C_i |\underline{\epsilon}| (\underline{I} - 3\underline{a}_2) \qquad (13) $$

where $\omega$ is the vorticity tensor, $\beta$ is a function of the fiber’s aspect ratio, $|\underline{\epsilon}|$ is the equivalent strain rate and $C_i$ an empirical constant, named usually the interaction coefficient. Experimental data has shown that, for concentrated suspensions, kinetics of orientation is much slower than the one predicted by Folgar and Tucker’s model. More recently, other authors introduced modified versions such as the reduced strain closure
model [22], the anisotropic rotary diffusion model [23] or a model taking into account fiber-fiber interactions [24]. Most theories have derived the following expression for the stress tensor [25]:

$$ \underline{\sigma} = -pI + 2\eta_f [\underline{\varepsilon} + N_p\underline{\varepsilon} : \underline{a}_4 + N_s(\underline{\varepsilon}\underline{a}_2 + \underline{a}_2\underline{\varepsilon})] \quad (14) $$

where $N_p$ and $N_s$ are parameters depending on the fluid's viscosity, fiber aspect ratio, orientation and concentration. A review is given in [14].

Fig.4(a), (b) and (c) illustrate, for a complex geometry and at the end of the filling stage, the principal components of $\underline{a}_2$ [15, 16]. It allows a precise knowledge of the fiber orientation distribution everywhere in the part and, starting from this orientation distribution, we may deduce the thermo-elastic properties of the part.

As a first constitutive law, one may use a linear anisotropic thermo-elastic relation (supposing that each point of the part cools from $T$ to $T_0$):

$$ \underline{\sigma}(T_0) = \underline{C}(\underline{a}_2, T_0) [\underline{\varepsilon}(T_0) - \underline{\alpha}(\underline{a}_2, T_0)\Delta T] \quad (15) $$

Components of both the stiffness tensor $\underline{C}$ and the dilatation tensor $\underline{\alpha}$ can be expressed, according to [26], as a function of the components of $\underline{a}_2$ and $\underline{a}_4$:

$$ C_{ijkl} = b_1 a_{ijkl} + b_2(a_{ij}\delta_{kl} + a_{kl}\delta_{ij}) + b_3(a_{ik}\delta_{jl} + a_{il}\delta_{jk} + a_{jk}\delta_{il} + a_{jl}\delta_{ik}) + b_4\delta_{ik}\delta_{kl} + b_5(\delta_{ik}\delta_{jl} + \delta_{il}\delta_{jk}) \quad (16) $$

$$ \alpha_{ij} = P_1 a_{ij} + P_2 \delta_{ij} \quad (17) $$

The different constants $b_i$ are identified for an unidirectional composite using, for example, Mori-Tanaka homogenization theory [27], being the stiffness tensor of the composite $\underline{C}$ expressed as a function of the stiffness tensors of the matrix, $\underline{C}_m$, and of the fiber $\underline{C}_f$. In the same way, parameters $P_1$ and $P_2$ can be expressed as a function of the longitudinal and transverse thermal dilatation coefficients, $\alpha_1$ and $\alpha_2$.

Fig.4 shows, for the previous example, the distribution of the elastic modulus in the different directions, computed from the stiffness tensor, and using the described semi-analytic methodology [15, 16].

## 3 Advanced numerical techniques and macroscale simulations

For multiphase and multiphysics flows as those encountered in the referred processes, either with short or textile reinforcements, different numerical techniques have been proposed to solve the conservation equations, coupled to the constitutive ones. The most popular ones are the finite difference method (FDM) [28], the finite element method (FEM) [29, 30, 31, 32, 14], the boundary element method (BEM) [33, 34, 35,
<figure><img src="image_5.png" /><figcaption>Figure 4: Illustration of fiber orientation and elastic modulus calculation on a prototype part (courtesy of Schneider Electric), where the material used was a glass fiber reinforced polycarbonate, with 20% fiber content, showing the results on half of the geometry. Distribution of the different diagonal components of $a_2$ at the end of the filling stage: (a) $a_{11}$, flow direction, with the injection gate represented on the right (b) $a_{22}$, height direction (c) $a_{33}$, width direction (d) $E_{11}$ (e) $E_{22}$ and (f) $E_{33}$ [16].</figcaption></figure>

[36] the smooth particle hydrodynamics (SPH) method [37], as well as reduced order methods [38, 39]. Review of the literature for injection molding (mostly used for short fiber composites) may be found in [32] and for liquid composite molding of continuous reinforcements in [40]. Even if computational cost may be greater, finite elements remain the most popular solution due to its capability to handle complex geometries and are the ones used and described in this section.

In most cases, the resin does not completely fill the cavity at the beginning of the process and it is necessary to simulate an impregnation or injection step, in particular to determine void formation and vent placement. This mold filling simulation is a typical problem of a moving boundary (fluid-air interface) with the computed velocity field. To solve it, one may find in the literature two main approaches: the computation on a moving domain or an Eulerian framework, where there is one computational domain composed of all the phases; the displacement of the phases as well as their deformations, which is a Lagrangian computation of the velocity field, meaning generally that there is one computational domain *per* phase. Numerical simulation is still strongly dependent on the meshing capabilities of complex geometries and one of the drawbacks of Lagrangian approaches is that it often constrains the enforcement of the boundary or interface meshes in the volume mesh. This task becomes more and more difficult since, due to domain deformation, remeshing is required at almost each time
step.

In Eulerian methods, the whole domain contains all the phases (fluid, air, mold, ...) and boundaries are represented in an implicit way. It means that the boundary is not given anymore by a surface mesh or any explicit representation but implicitly by a scalar field which value can be accessed anywhere in the domain. Two representative methods using such a field are the Volume Of Fluid (VOF) [41] and the level-set [42] methods. In the following, a modified version of the latter will be detailed [43].

## 3.1 Implicit boundaries

To solve conservation equations in an Eulerian framework, one needs to define material properties everywhere in the computational domain. To detail this multiphase context and as a matter of simplification, let us consider a domain $\omega$ totally included in the larger one, $\Omega$, being $\Gamma = \delta\omega$ the boundary of $\omega$. $\eta$, $\rho$, $K$ or $\lambda$, among others, are heterogeneous fields, defined by their values inside and outside $\omega$. In fact, these properties are discontinuous from one phase to the other, and may be represented in a general way as a variable field $\xi$ within the mesh and through a function $H$, such that

$$ \xi = \xi_{\omega} H + \xi_{\Omega \setminus \omega} (1 - H) \qquad (18) $$

where $H$ is an Heaviside function, usually equal to 0 outside $\omega$ and 1 inside (Fig.5(b)). Dealing with such discontinuity is not affordable by any standard numerical method and a smooth transition around the interface has been proposed [43].

<figure><img src="image_6.png" /><figcaption>Figure 5: Phase functions used to describe, in an implicit way, the different phases of implicit boundaries computation. In this example, two phases are present [44].</figcaption></figure>

## 3.2 Immersed subdomains and regularisation

From the previous description, one may say that $\omega$ is immersed into $\Omega$. To represent $\omega$, an implicit function has been chosen, being the immediate candidate the signed distance function, Fig.5(a):
$$ \bar{d}(x, \Gamma) = \begin{cases} d(x, \Gamma) & \text{if } x \in \omega \\ -d(x, \Gamma) & \text{if } x \notin \omega \end{cases} \qquad (19) $$

Reciprocally, the interface is completely defined by the zero value of the distance function, and the knowledge of $\omega$ is ensured by the definition of this function everywhere in $\Omega$. From $\bar{d}$, one may calculate the Heaviside function, Fig.5(b), by:

$$ H(\bar{d}) = \begin{cases} 1 & \text{if } \bar{d}(x, \Gamma) \ge 0 \\ 0 & \text{if } \bar{d}(x, \Gamma) < 0 \end{cases} \qquad (20) $$

The Heaviside function enables to separate strictly the domains, but the counterpart is that it reintroduces the discontinuity which has been left by the use of an implicit function. The following implicit function:

$$ u_{\epsilon} = u(d, \epsilon) = \epsilon \tanh\left(\frac{d}{\epsilon}\right) \qquad (21) $$

enables the introduction of a thick interface which thickness is related to the parameter $\epsilon$. Moreover, $u_{\epsilon}$ approximates $d$ and we recover $\Gamma$ from it, for whatever $\epsilon$. Then, a smoothed Heaviside function can be introduced as well by:

$$ H_{\epsilon} = \frac{1}{2}\left(1 + \frac{u_{\epsilon}}{\epsilon}\right) \qquad (22) $$

### 3.3 Multiphase flow and thermokinetic numerical resolution

Let us come back to the resolution of the conservation equations, and firstly the flow ones, under the Brinkman form, Eq.(4). For the case of a resin impregnating a porous media (one fluid, $f$, entering the cavity and air, $a$, venting it), material properties are controlled by the thickness parameter $\epsilon$, such that

$$ \eta = \eta_f H_{\epsilon} + \eta_a (1 - H_{\epsilon}) \\ \underline{K} = \underline{K}_f H_{\epsilon} + \underline{K}_a (1 - H_{\epsilon}) \qquad (23) $$

Brinkman's problem may be solved using a mixed finite element method [45]. The authors proposed a stabilized formulation which approximates $\vec{v}$ and $p$ with finite dimensional spaces of continuous piecewise polynomials, such that $\vec{v}_h|_K \in P^1(K)^d$ and $p_h|_K \in P^1(K)$, $\forall K \in \mathcal{K}$. It enriches the velocity space with a space of bubbles (the finer scale), using a polynomial function of order one, with a value equal to unity at the center of the element and vanishing at its boundary. The enrichment performed stabilizes the formulation for both the limit Stokes and the Darcy cases, with a modification of bubble stabilization for both regions, to be compatible with stability and unicity of the solution conditions. The element contribution to the linear system arising can be written, in the matrix form [45]:
$$ \begin{bmatrix} A_{vv} - A_{vb}(A_{bb}^{-1})^t A_{vb} & A_{vp} - A_{vb}(A_{bb}^{-1})^t A_{bp} \\ A_{vp} - A_{bp}(A_{bb}^{-1})^t A_{vb} & -A_{bp}(A_{bb}^{-1})^t A_{bp} \end{bmatrix} \begin{bmatrix} V \\ P \end{bmatrix} = \begin{bmatrix} F_v \\ P \end{bmatrix} \quad (24) $$

In what concerns the solution of the energy equation, as well as the ODE representing the evolution of the different internal variables (like $\alpha$ or $a_2$), these are generally equations of the convection-diffusion-reaction type. Classical techniques generate numerical instabilities, while diffusion problems may cause oscillations during the treatment of thermal shocks.

To avoid stability difficulties due to convection, some techniques have been developed that behave well on highly convective problems but are limited by their tendency to denature the solution. Techniques to avoid oscillations during the treatment of thermal shocks are well-known: one must adapt the mesh size or the time step in the direction of the gradient. To be optimal on the solution of the two issues, two different stabilization techniques are widely used in the literature: the SUPG (Streamline Upwind Petrov Galerkin Method [46]) and the RFB (Residual Free Bubbles [47]) method.

In all the cases, the proposed formulation leads to the resolution of a sparse linear system, which represents an important fraction of the overall simulation time. To solve it, the PETSc library [48] may be used, with a classical preconditioned (ILU, Block Jacobi,...) iterative method of the Krylov type (Conjugate Residual or GMRES).

## 3.4 Composite forming simulation illustrations

### 3.4.1 Compression molding

Compression molding is a high-volume, high-pressure method which is suitable for molding complex and high performing parts. Thermosets (mainly) or thermoplastics may be compression molded in different types of reinforcements, like unidirectional fabrics, woven textiles, random mats or chopped fibers. The use of composite materials like Sheet Molding Compounds (SMC), which are compressed, has increased in the last years, in particular in the automotive industry.

In the literature, several authors have modeled compression molding, but assuming important hypothesis: a 2 – 1/2D model; a dilute regime for fiber concentration; decoupling between cure, temperature and flow; no anisotropy effects, in particular in the thermal conductivity; perfect contact between the SMC and the mold (and between the different layers of SMC). More recently, 3D simulations involving full coupling have been shown [49], with an illustration given in Fig.6. In this case, a Sheet Molding Compound has been compressed. At the beginning, several layers where cut and disposed in three different packs. Before curing begins to influence flow, one observes the high deformation and the flow of the samples, under heating conditions. On the right, complexity of the orientation obtained is plotted at the end of the compression.
<figure><img src="image_7.png" /><figcaption>Figure 6: Simulation of SMC compression molding using Rem3D software [50]: (a) evolution of the material inside the cavity (b) orientation distribution at the end of filling, represented by the eigenvecteurs of $\underline{a}_2$ multiplied by its respective eigenvalues.</figcaption></figure>
### 3.4.2 Resin transfer molding

For the production of structural composites, one of the most promising routes is resin transfer molding (RTM). In this process, the reinforcement is laid in the mold cavity after which the mold closes. Once the textile has been enough heated, resin is injected from one or several gates. Consolidation occurs before the mold opens.

Several variations of this process exist. One of the most used, in particular for very large parts, is infusion: when the clamping forces required to keep the mold closed are very high, a plastic film is used as an upper mold and vacuum is imposed, so that the resin is injected. One other variation is HP-RTM, high pressure resin transfer molding, used to produce parts with very high fiber volume fraction, i.e. very low permeability, enabling a reasonable injection time [51].

HP-RTM (and also Compression RTM, C-RTM) simulation and experiments have been studied in [51] using the above described finite element method. The first test carried out was the injection of a resin into a small disk and a comparison between analytic and numerical solutions was made, Figs.7(a) and (b), showing a good agreement. In what concerns temperature, a strong gradient appears near the injection gate, and a strong skin flow phenomenon is developed, as seen in Figs.7(c) and (d). Temperature difference between a simulation made with or without thermo-rheo-kinetical coupling is also illustrated, showing its importance since, in this last case, skin flow does not occur.

## 3.5 Parallel mesh adaptation and high performance computing

Solutions obtained, for multiphase flows with very different properties, converge towards a sharp or exact interface solution when $\epsilon \to 0$. This means that the proposed methodology is based on generating a thickness $\epsilon$ small enough so that computations are accurate. Since to represent $u_\epsilon$ one needs also several layers of elements within $\epsilon$, this may only be attained using anisotropic adaptive meshing and parallel implementations.

### 3.5.1 Automatic anisotropic adaptation

Anisotropic mesh adaptation enables to capture discontinuities or gradients of the solution in flow, thermal, kinetical situations [52, 43], with a good accuracy at a very low extra number of elements. Building unstructured anisotropic meshes may be based on local mesh modifications, using a metrics field to redefine the lengths [53] and a recently developed procedure optimal in terms of remeshing is here described.

Let $\mathcal{N}$ be a finite set of node numbers and $\mathcal{K}$ a mesh topology or a set of elements that provides the element-node connectivity. Let us consider only the mesh topologies composed of simplex elements (triangles, tetrahedra,...), $K$ being a $d$-simplex ($d$ is the spatial dimension) composed of $D = d + 1$ vertices. The mesh generator described adapts an initial mesh by iteratively modifying and optimizing it locally [53]. This technique consists in the improvement of the quality of local cavities composed of
<figure><img src="image_8.png" /><figcaption>Figure 7: Simulation of resin transfer molding using Rem3D software [50]: (a) evolution of the material inside the cavity (b) pressure evolution throughout filling (c) velocity norm in a thickness cut (d) temperature in a thickness cut, with or without thermo-rheo-kinetical coupling [51].</figcaption></figure>
clusters of elements by remeshing. Two principles are enforced: the minimal volume, which gets the conformity of the mesh, with no element overlaps; the geometrical quality distribution, ensuring that, if the minimal volume provides several possible cavity re-triangulations, then the one chosen is such that it provides the best worst element. The geometrical quality, $Q(K)$, is evaluated for each element $K$ and varies from 0 (worst) to 1 (best quality). It depends on $M_K$, the metric of the element $K$, which is thus the information needed by the mesher to generate a new mesh.

Construction of $M$ may be done using the interpolation error analysis on a solution field, for example, $u$, which may be a velocity, a phase function or other variable, or a vector containing all. Following [52], let us define the edge vector as $\vec{X}^{ij} = \vec{X}^j - \vec{X}^i$, connecting nodes $i$ and $j$, and the edge length distribution tensor at $i$ as

$$ \underline{X}^i = \frac{1}{|\Gamma(i)|} \sum_{j \in \Gamma(i)} \vec{X}^{ij} \otimes \vec{X}^{ij} \qquad (25) $$

where $\Gamma(i)$ represents the nodes that share an edge with $i$. It may be established that the metrics computed at the node is: $\underline{M}^i = \frac{1}{d}(\underline{X}^i)^{-1}$.

Let us suppose that one wishes to generate a new mesh that equidistributes in each edge a constant target interpolation error, $e_{\text{target}}$. It has been previously shown [54] that this target error can be directly linked to a fixed number of nodes $N$ (which may also be related with the available CPU power), $e_{\text{target}} = e_{\text{target}}(N)$. An a posteriori interpolation error analysis states that the error along an edge $ij$ on the solution field $u$, $e^{ij}$, is $e^{ij} = \nabla u^{ij} \cdot \vec{X}^{ij}$, where $\nabla u^{ij} = \nabla u^j - \nabla u^i$, nodal gradients in $i$ and $j$ obtained with a gradient recovery procedure [52].

The optimal stretch factor, $s^{ij}$ that must be applied to $\underline{X}^{ij}$ to attain $e_{\text{target}}$ may now be used to compute the target metric given to the mesher, which allows the new mesh (sequential) generation:

$$ s^{ij} = \frac{e_{\text{target}}(N)}{e^{ij}} \quad \text{and} \quad \tilde{M}^i = \frac{1}{d} \left[ \frac{1}{|\Gamma(i)|} \sum_{j \in \Gamma(i)} (s^{ij})^2 \vec{X}^{ij} \otimes \vec{X}^{ij} \right]^{-1} \qquad (26) $$

An illustration of the type of result obtained in adaptation around one or several layers of fibers is illustrated in Fig.8, where one may observe the mesh density around the tows in (a) and (c) or the accuracy in the representation of the interface in (b).

### 3.5.2 Parallel computing

For complex 3D geometries and complex physical laws as the ones concerning composites, parallel computation remains essential. It allows us to run complete simulations with a reasonable precision, inaccessible (for memory and time limits) with sequential runs. It also makes possible to speed up the simulation run such that the global computation time will still be acceptable. In the last years, processor's performance has not increased by improving the clock rate but by multiplying the number
<figure><img src="image_9.png" /><figcaption>Figure 8: Adapted mesh to accurately represent, in an implicit way, the interface between tows and fluid: (a) mesh (b) isozero value of $u_{\epsilon}$, distance to the fiber/resin interface and (c) mesh and $u_{\epsilon}$ representation [55].</figcaption></figure>

of cores in a CPU. Actual top supercomputers contain several hundreds of thousands to millions of cores with hundreds of TB to PB of memory. It is thus necessary to develop fully parallel applications that follow, at least, this multi-core CPU evolution. Between the existing parallelization techniques, we have chosen to refer those developed and used in the former and following examples [56]. They are based on the decomposition of the computational domain into several sub-domains, where each is affected to one core. In the simulations illustrated above, most of the computation time is spent on remeshing and on the resolution of several large linear systems.

In what concerns meshing, the chosen parallelization procedure uses the sequential mesher described before in a massive parallel context, by following an iterative two successive steps procedure [56]: independent adaptive meshing of each subdomain with blocked interfaces; constrained repartitioning and interface displacement. Even if several iterations are done, work per iteration decreases quickly. When a good quality mesh is obtained, a finite element load-balancing repartitioning step is performed. All fields defined on the initial mesh may then be mapped on the new one.

For linear system resolution, parallelization has been performed by interfacing the parallel partitionner with the PETSc library [48] and also by applying a multi-grid preconditionner. An example illustrating more in detail the capabilities of parallel computing and the performances of this method is shown in Sec.4.2. It has already allowed solving direct flow simulations on fiber samples that required meshes of around 780 millions of nodes and linear systems of 3.6 billions of unknowns [57].
# 4 Determination of equivalent properties and microscale simulations

To obtain accurate models at the macroscopic scale, one often performs a homogenization step. Semi-analytical homogenization, as the one presented in Sec.2.3, is a well established and practical method to compute the effective properties of fiber reinforced materials, in particular at the solid state.

However, if one wishes to later extend it to the non-linear case (for example, by considering that the matrix is viscoelastic) or to the flow regime, these methods may not give enough accuracy. The use of micro-mechanical and micro-rheology finite element computations has thus increased in the last decade [58, 59, 60, 61]. In these computations, properties are obtained by performing finite element simulations on three dimensional 'Representative Elementary Volumes' (REV) of the microstructure. Generating several REVs allows to perform a sensitivity analysis to the topology of the REV (size, number and spatial distribution of fibers or tows), as well as to the mechanical or rheological solicitations. From all the numerical techniques used to generate REVs, direct methods from 3D images will be the focus of this section.

## 4.1 Generation of representative numerical samples

Recently, imaging techniques, like X-ray tomography, have been used to obtain precise geometrical and topological information on microstructures, which allowed extraordinary contributions in material science. The passage from the image to a finite element mesh allows, on one hand, the construction of a numerical representation but, on the other hand, may also reduce the size of the stored data. The first step (before building a mesh) is usually to perform the image segmentation, that is, to partition the image into its multiple objects, so that it becomes easier to analyze it [62]. The result of this segmentation can be regions of the image where the pixels are similar with respect to a characteristic or property (like, for example, the density), or extracted contours. Secondly, a mesh is built from the segmented image by using different techniques.

Recent work [63] proposes a novel technique to simultaneously segment and construct a finite element mesh, using directly the 3D image data, and applies it to fibrous materials. For that, the Immersed Image Method is built, which is the interpolation of the image's pixel/voxel values (a sort of topographic distance function) on the nodes of an initial mesh. The obtained field can approximately represent the original image, depending on the discretization and on the overall number of nodes. Optimization of the number of necessary nodes can be performed through anisotropic mesh adaptation using the above described procedure, so that the obtained field accurately represents the original image. To further progress towards segmentation, the mesh is adapted not on the interpolated topographic distance function (the pixel/voxel value), but on $u_e$. To build this function, a reinitialization method has been implemented. This procedure, coupled to mesh adaptation, allows segmentation but allows also to build up a mesh
with a smooth representation of a phase function distribution, usable in the aimed type of application, Fig.9.

<figure><img src="image_10.png" /><figcaption>Figure 9: Illustration of the results obtained for mesh adaptation directly from an image: (a) original composite image, with 200 × 200 × 220 voxels (b) computed distance function (c) computed regularized Heaviside function (d) final adapted mesh, with 50 000 nodes.</figcaption></figure>

## 4.2 Permeability of a composite

With computed velocity and pressure on a porous REV, permeability can be determined (at the micro or meso scales). For semi-analytic homogenization, three different approaches are most often used: the capillary model; models arising from the hydrodynamic lubrication theory; the cell model. Hydrodynamic lubrication theory models suppose the fibrous media as a regular cylinder spacing geometry and compute an analytic solution for the Stokes equations inside this REV. The most popular one for the transverse permeability of the array of cylinders is Gebart’s expression [64].

In numerical homogenization, Darcy’s equation, defining permeability is directly considered. Thus, for at least three principal directions, one must compute
$$K = \eta_f \left( \frac{\langle \| \vec{v} \| \rangle}{\langle \| \nabla p \| \rangle} \right) \qquad (27)$$

where $\langle \| \vec{v} \| \rangle$ and $\langle \| \nabla p \| \rangle$ are the averaged velocity and pressure gradient fields for the considered direction, only computed in the fluid domain.

At the microscopic scale, numerical results can be compared with Gebart's ones, by considering a regular array of fibers, as illustrated in [55], for square and triangular packings. A pressure gradient is imposed from left to right and the normal velocity in the upper and lower planes is zero. Inside the fibers, velocity is also zero. Fig.10 illustrates a comparison between numerical and semi-analytic results, for different fiber volume fractions. Overall results show that the computed values are closer to lubrication model results for a low fiber volume fraction and closer to cell models for high fiber volume fraction.

<figure><img src="image_11.png" /><figcaption>Figure 10: Permeability determination for regular arrays of fibers and comparison with semi-analytic solutions for (a) a square packing and (b) a triangular packing, for different volume fraction of fibers [55].</figcaption></figure>

However, fibers may be randomly placed in the fabric and do not generally follow a strict organized structure (even for the case where one wishes to impose a certain alignment). To assess the feasibility and capability of the overall procedure, a 3D test case is presented, Fig.11 [63], which complexity implies that no statistic or equivalent organized REV reconstruction is possible.

An image (900 × 900 × 220 voxels) issuing from 3D X-Ray tomography [65] concerning the mat-reinforced sample is treated and a phase function equal to 1 inside the fibres (solid) and 0 in the resin (fluid) phases is built. Finite element resolution of these equations is done using the above described mixed stabilized finite element method, with a generated mesh of 5 millions of nodes, from the image. In terms of computational time, mesh generation, with the overall adaptation procedure, was performed in 3 hours, whereas the Navier-Stokes resolution took only 10 min, on 96 CPUs. Results detailed in Fig.11 were obtained by imposing a pressure gradient through the sample
on one direction, and imposing no-velocity in the fibers, showing that: the mesh is well adapted at the fluid-solid interface; the definition of this interface is very accurate, as illustrated by the isovalue zero of $u_{\epsilon}$; the computed velocity field enhances well the acceleration in the interstitial spaces.

<figure><img src="image_12.png" /><figcaption>Figure 11: Permeability determination for a 3D complex sample: (a) original image with 178.2 millions of voxels, (b) von Mises equivalent of the computed metrics field, (c) final mesh of about 5 millions of nodes, (d) isovalue zero of the regularized Heaviside, (e) computed velocity and (f) pressure fields [63].</figcaption></figure>
4.3 Stiffness tensor determination

A discussed previously, direct numerical thermoelastic simulations on REVs also al-
low the determination of the stiffness and dilatation tensors of the sample. Like before,
let us consider a single mesh which contains both fibers and matrix, Fig.12(a).
A tied contact between fiber and matrix is supposed, as well as between two inde-
pendent fibers. Knowing the stiffness and the thermal dilatation for fibers and resin,
Eq.18 is used to derive the stiffness and the thermal dilatation everywhere in the REV.
Several elastic thermo-mechanical tests are applied on this REV and numerically ho-
mogenized stiffness and thermal dilatation tensor are obtained. Numerical resolution
of this heterogeneous elastic problem may be performed using a mixed finite element
method [16], with forces imposed on each side of the REV. To avoid the imposition
of periodic boundary conditions, homogenization is performed on a smaller included
REV than the one where computations are done.

This type of computations have been performed by [16], who has shown that one
important issue was the optimization of the number of fibers in the REV in order to
obtain a stabilized result. For unidirectional aligned fibers, Fig.12(b) plots the com-
puted modulus in one of the principal directions, as a function of the number of fibers.
It is compared with Mori-Tanaka semi-analytic result showing that, for the sample
represented, 22 fibers provided a convergent result, using a very fine adapted mesh.

<figure><img src="image_13.png" /><figcaption>Figure 12: Stiffness and dilatation tensor determination (a) typical Representative Elementary Volume: the cube represents the computational domain (meshed) and in red, the fibre-matrix interface (b) Young modulus as a function of the number of fibres in the REV (the Mori-Tanaka modulus is 6.21 GPa), for a given mesh size (very fine, in that case) [16].</figcaption></figure>
# 5 Conclusions

In this chapter, we have detailed conservation and constitutive equations involved in multiphysics, multidomain and multiphase modeling of composites, which may be either short, long or continuously reinforced. Advanced numerical techniques involved in flow, thermo-rheo-kinetical and structure development computations were presented. They include: mixed stabilized finite elements; Eulerian approaches and modified level-set techniques; automatic anisotropic mesh adaptation and parallel computing. Illustrations from different authors were presented, for forming at the macroscopic scale and property's direct numerical homogenization at the macroscopic or mesoscopic scales. Examples included: short fiber reinforced thermoplastics injection molding, Sheet Molding Compound compression molding, resin transfer molding; permeability determination at the microscopic and mesoscopic scales; stiffness tensor determination. As a conclusion, huge steps forwards have been achieved to consider complex physics in molding simulations, but other advanced phenomena need to be accounted which will surely induce new mathematical and numerical developments.

# References

[1] S. G. Advani. Opportunities and challenges of multiscale modeling and simulation in polymer composite processing. *International Journal of Material Forming*, 2(1):39–44, 2009.

[2] S.G. Advani and E. M. Sozer. *Process Modeling in Composites Manufacturing, Second Edition*. CRC Press, 2010.

[3] M.A. Khan, T. Mabrouki, E. Vidal-Sall, and P. Boisse. Numerical and experimental analyses of woven composite reinforcement forming using a hypoelastic behaviour. application to the double dome benchmark. *Journal of Materials Processing Technology*, 210(2):378 – 388, 2010.

[4] S. Allaoui, P. Boisse, S. Chatel, N. Hamila, G. Hivet, D. Soulat, and E. Vidal-Salle. Experimental and numerical analyses of textile reinforcement forming of a tetrahedral shape. *Composites Part A: Applied Science and Manufacturing*, 42(6):612 – 622, 2011.

[5] K. Hitti. *Direct numerical simulation of complex Representative Volume Elements (RVEs): Generation, Resolution and Homogenization*. PhD thesis, MINES ParisTech, 2011.

[6] H. Darcy. *Les fontaines publiques de la ville de Dijon*. Dalmont, Victor, 1856.

[7] H.C. Brinkman. A calculation of the viscous force exerted by a flowing fluid on a dense swarm of particles. *Applied Scientific Research*, 1(1):27–34, 1949.
[8] J.-A. Auriault. On the domain of validity of Brinkman's equation. *Transport in Porous Media*, 79:215-223, 2009.

[9] S. Souror and M.R. Kamal. Differential scanning calorimetry of epoxy cure: isothermal cure kinetics. *Thermochimica Acta*, 14(1-2):41-59, 1976.

[10] J.F. Agassant, P. Avenas, J.P. Sergent, B. Vergnes, and M. Vincent. *La mise en forme des matires plastiques*. Lavoisier, 1996. **matières**

[11] D. Serrano, J. Peyrelasse, C. Boned, D. Harran, and P. Monge. Application of the percolation model to gelation of an epoxy resin. *Journal of Applied Polymer Science*, 39(3):679-693, 1990.

[12] A.T. DiBenedetto. Prediction of the glass transition temperature of polymers: A model based on the principle of corresponding states. *Journal of Polymer Science*, 25(9):1949-1969, 1987.

Venet, C. and Silva, L.

[13] E. Foudrinier and L. Venet, C. and Silva. 3d computation of reactive moulding processes. *International Journal of Material Forming*, 1(1):735-738, 2008.

[14] P. Laure, L. Silva, and M. Vincent. Modelling short fibre polymer reinforcements for composites. In P. Boisse, editor, *Composite reinforcements for optimum performance*, pages 619-650. Woodhead Publishing Limited, 2011.

[15] H. Miled, L. Silva, J.F. Agassant, and T. Coupez. Numerical simulation of fiber orientation and resulting thermo-elastic behavior in reinforced thermo-plastics. In *Mechanical Response of Composites*, volume 10 of *Computational Methods in Applied Sciences*, pages 293-313. Springer Netherlands, 2008.

[16] H. Miled, L. Silva, T. Coupez, and J.F. Agassant. Injection molding of fibre reinforced thermoplastics: Integration of fibre orientation and mechanical properties computations. *International Polymer Processing*, 27(5):547-556, 2012.

[17] A. Ammar and F. Chinesta. A particle strategy for solving the fokker-planck equation modelling the fiber orientation distribution in steady recirculating flows involving short fiber suspensions. In M. Griebel and M.A. Schweitzer, editors, *Meshfree Methods for Partial Differential Equations II*, volume 43 of *Lecture Notes in Computational Science and Engineering*, pages 1-15. Springer Berlin Heidelberg, 2005.

[18] G.L. Hand. A theory of anisotropic fluids. *Journal of Fluid Mechanics*, 13(1):33-46, 1962.

[19] G.G. Lipscomb, R. Keunings, G. Marrucci, and M.M. Denn. Continuum theory for fiber suspensions. In *Advances in rheology: Polymers : IX International Congress on Rheology*, 1984.

[20] G. Jeffery. The motion of ellipsoidal particles immersed in viscous fluid. *Proceedings of the Royal Society A*, 102(715):161-179, 1922.
[21] F. Folgar and C.L. Tucker. Orientation behavior of fibers in concentrated suspensions. *Journal of Reinforced Plastics and Composites*, 3(2):98-119, 1984.
correct "O'Gara"

[22] J. Wang, J.F. OGara, and C.L. Tucker. An objective model for slow orientation kinetics in concentrated fiber suspensions: Theory and rheological evidence. *Journal of Rheology*, 52(5):1179-1200, 2008.

[23] J.H. Phelps and C.L. Tucker III. An anisotropic rotary diffusion model for fiber orientation in short- and long-fiber thermoplastics. *Journal of Non-Newtonian Fluid Mechanics*, 156(3):165–176, 2009.

[24] J. Frec, G. Ausias, M. C. Heuzey, and P. J. Carreau. Modeling fiber interactions in semiconcentrated fiber suspensions. *Journal of Rheology*, 53(1):49–72, 2009.

[25] C.L. Tucker. Flow regimes for fiber suspensions in narrow gaps. *Journal of Non-Newtonian Fluid Mechanics*, 39(3):239–268, 1991.

[26] S.G. Advani and C.L. Tucker. The use of tensors to describe and predict fiber orientation in short fiber composites. *Journal of Rheology*, 31(8):751–784, 1987.

[27] T. Mori and K. Tanaka. Average stress in matrix and average elastic energy of materials with misfitting inclusions. *Acta Metallurgica*, 21(5):571–574, 1973.

[28] F. Trochu and R. Gauvin. Limitations of a boundary fitted finite difference method for the simulation of the resin transfer molding process. *Journal of Reinforced Plastics and Composites*, 11:772–786, 1992.

[29] M. V. Bruschke and S. G. Advani. A finite element/control volume approach to mold filling in anisotropic porous media. *Polymer Composites*, 11(6):398–405, 1990.

[30] F. Trochu, E. Ruiz, V. Achim, and S. Soukane. Advanced numerical simulation of liquid composite molding for process analysis and optimization. *Composites Part A: Applied Science and Manufacturing*, 37(6):890–902, 2006.

[31] P. Simacek and S.G Advani. Modeling resin flow and fiber tow saturation induced by distribution media collapse in vartm. *Composites science and technology*, 67(13):2757–2769, 2007.

[32] L. Silva, J.F. Agassant, and T. Coupez. Three-dimensional injection molding simulation. In M.R. Kamal, A.I. Isayev, and S.J. Liu, editors, *Injection molding, technology and fundamentals*, pages 599–651. Carl Hanser Verlag, 2009.

[33] M.K. Um and W.I. Lee. A study on the mold filling process in resin transfer molding. *Polymer Engineering and Science*, 31(11):765–771, 1991.

[34] Y.E. Yoo and W.I. Lee. Numerical simulation of resin transfer molding process using boundary element method. *Polymer Composites*, 17(3):368–374, 1996.
[35] S. Soukane and F. Trochu. Application of the level set method to the simulation of resin transfer molding. *Composites Science and Technology*, 66(7):1067–1080, 2006.

[36] R. Gantois, A. Cantarel, G. Dusserre, J.N. Felices, and F. Schmidt. Bem-based models to simulate the resin flow at macroscale and microscale in lcm processes. *Polymer Composites*, 34(8):1235–1244, 2013.

[37] S. Comas-Cardona, P. Groenenboom, C. Binetruy, and P. Krawczak. A generic mixed fe-sph method to address hydro-mechanical coupling in liquid composite moulding processes. *Composites Part A: Applied Science and Manufacturing*, 36(7):1004–1010, 2005.

[38] C. Ghnatios, F. Chinesta, E. Cueto, A. Leygue, A. Poitou, P. Breitkopf, and P. Villon. Methodological approach to efficient modeling and optimization of thermal processes taking place in a die: Application to pultrusion. *Composites Part A: Applied Science and Manufacturing*, 42(9):1169–1178, 2011.

[39] E. Abisset-Chavanne and F. Chinesta. Toward an optimisation of the reactive resin transfer molding process: thermo-chemico-mechanical coupled simulations. *International Journal of Material Forming*, 7(2):249–258, 2014.

[40] C.H. Park. 15 - numerical simulation of flow processes in composites manufacturing. In Philippe Boisse, editor, *Advances in Composites Manufacturing and Process Design*, pages 317–378. Woodhead Publishing, 2015.

[41] C.W. Hirt and B.D. Nichols. Volume of fluid (vof) methods for the dynamics of free boundaries. *Journal of Computational Physics*, 39(1):201–225, 1981.

[42] S. Osher and J.A. Sethian. Fronts propagating with curvature-dependent speed: Algorithms based on hamilton-jacobi formulations. *Journal of Computational Physics*, 79(1):12–49, 1988.

[43] T. Coupez, L. Silva, and E. Hachem. Implicit boundary and adaptive anisotropic meshing. In *New Challenges in Grid Generation and Adaptivity for Scientific Computing*. SEMA SIMAI Springer Series, Vol. 5, 2015.

[44] J.X. Zhao, T. Coupez, E. Decencire, D. Jeulin, D. Cardenas, and L. Silva. Multiphase mesh generation from 3d image by anisotropic mesh adaptation and re-distancing equation. *Computer Methods in Applied Mechanics and Engineering*, submitted, 2015.

[45] L. Silva, G. Puaux, M. Vincent, and P. Laure. A monolithic finite element approach to compute permeability at microscopic scales in lcm processes. *International Journal of Material Forming*, 3(1):619–622, 2010.
[46] A.N. Brooks and T.J.R. Hughes. Streamline upwind/petrov-galerkin formulations for convection dominated flows with particular emphasis on the incompressible navier-stokes equations. *Computer Methods in Applied Mechanics and Engineering*, 32(13):199–259, 1982.

[47] F. Brezzi, L.P. Franca, and A. Russo. Further considerations on residual-free bubbles for advective-diffusive equations. *Computer Methods in Applied Mechanics and Engineering*, 166(12):25–33, 1998.

[48] S. Balay, S. Abhyankar, M.F. Adams, J. Brown, P. Brune, K. Buschelman, L. Dalcin, V. Eijkhout, W.D. Gropp, D. Kaushik, M.G. Knepley, L.C. McInnes, K. Rupp, B.F. Smith, S. Zampini, and H. Zhang. PETSc Web page. http://www.mcs.anl.gov/petsc, 2015.

[49] L. Silva, L. Ville, and P. Laure. Numerical simulation of compression moulding using an adaptative monolithic finite element approach. In *17th ESAFORM Conference*, 2014.

[50] L. Ville. REM3D Web page. http://www.transvalor.com, 2015.

[51] A. Bissuel. Process modelling of HP-RTM and Compression RTM using Rem3D. Technical report, Numerical simulation of flows at microscopic scale and mesoscopic scale in the RTM process, 2013.

[52] T. Coupez. Metric construction by length distribution tensor and edge based error for anisotropic adaptive meshing. *Journal of Computational Physics*, 230:2391–2405, 2011.

[53] T. Coupez. A mesh improvement method for 3D automatic remeshing. In N.P. Weatherill, P.R. Eiseman, J. Hause, and J.F. Thompson, editors, *Numerical Grid Generation in Computational Fluid Dynamics and Related Fields*, pages 615–626. Pineridge Press, 1994.

[54] T. Coupez and E. Hachem. Solution of high-reynolds incompressible flow with stabilized finite element and adaptive anisotropic meshing. *Computer Methods in Applied Mechanics and Engineering*, 267:65–85, 2013.

[55] G. Puaux. *Numerical simulation of flows at microscopic scale and mesoscopic scale in the RTM process*. PhD thesis, École Nationale Supérieure des Mines de Paris, 2011.

[56] T. Coupez, H. Digonnet, and R. Ducloux. Parallel meshing and remeshing. *Applied Mathematical Modelling*, 25:153–175, 2000.

[57] L. Silva, T. Coupez, and H. Digonnet. Massively parallel mesh adaptation and linear system resolution for multiphase flows. In *27th International Conference on Parallel Computational Fluid Dynamics*, 2015.
[58] J. LLorca and J. Segurado. Three-dimensional multiparticle cell simulations of deformation and damage in sphere-reinforced composites. *Materials Science and Engineering: A*, 365(12):267–274, 2004.

[59] J. Segurado and J. LLorca. A computational micromechanics study of the effect of interface decohesion on the mechanical behavior of composites. *Acta Materialia*, 53(18):4931–4942, 2005.

[60] C.J. Sun, P. Saffari, K. Sadeghipour, and G. Baran. Effects of particle arrangement on stress concentrations in composites. *Materials Science and Engineering: A*, 405(12):287–295, 2005.

[61] E. Totry, J.M. Molina-Aldaregua, C. Gonzlez, and J. LLorca. Effect of fiber, matrix and interface properties on the in-plane shear deformation of carbon-fiber reinforced composites. *Composites Science and Technology*, 70(6):970–980, 2010.

[62] L. Shapiro and G. Stockman. *Computer Vision*. Prentice Hall, 2001.

[63] L. Silva, J.X. Zhao, H. Digonnet, and T. Coupez. Flow simulations on 3d segmented images using reinitialization and anisotropic mesh adaptation. In *Computational Modeling of Objects Presented in Images. Fundamentals, Methods, and Applications*, pages 339–350. Springer International Publishing, 2014.

[64] B. R. Gebart. Permeability of unidirectional reinforcements for RTM. *Journal of Composite Materials*, 26(8):1100–1133, 1992.

[65] L. Orgéas, P.J.J. Dumont, J.P. Vassal, O. Guiraud, V. Michaud, and D. Favier. In-plane conduction of polymer composite plates reinforced with architecture networks of copper fibres. *Journal of Materials Science*, 47(6):2932–2942, 2012.