# Finite element implementation of a two-phase model for compression molding of composites

Pierre J.J. Dumont, Steven Le Corre, Laurent Orgéas, Denis Favier, Cyril Gaborit, Pierre Lory

► To cite this version:

Pierre J.J. Dumont, Steven Le Corre, Laurent Orgéas, Denis Favier, Cyril Gaborit, et al. Finite element implementation of a two-phase model for compression molding of composites. Revue Européenne des Éléments Finis, 2005, 14 (6-7), pp.885-902. 10.3166/reef.14.885-902. hal-01581689

HAL Id: hal-01581689

https://hal.science/hal-01581689

Submitted on 5 Sep 2017

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
# Finite element implementation of a two-phase model for compression molding of composites

Pierre Dumont* — Steven Le Corre** — Laurent Orgéas*** — Denis Favier**** — Cyril Gaborit***** — Pierre Lory******

*Laboratoire de Technologie des Composites et Polymères (LTC), Ecole Polytechnique Fédérale de Lausanne, Station 12, CH-1015 Lausanne, Switzerland

**GeM - Institut de Recherche en Génie Civil et Mécanique, CNRS - Ecole Centrale de Nantes, BP 92101, 44321 Nantes cedex 3, France.

***Laboratoire Sols-Solides-Structures (3S), CNRS - Universités de Grenoble BP 53, 38041 Grenoble, Cedex, France, Laurent.Orgeas@hmg.inpg.fr

****Technocentre Renault, TCR Lab045, 1 Avenue du Golf, 78288 Guyancourt Cedex, France

**ABSTRACT.** A two-phase approach is proposed to model the rheology of polymer glass-fiber com-
pounds such as SMC or GMT during processing. The anisotropic behavior of the composite,
which is related to the microstructure of the fiber network, is reduced to the simple case of
transverse isotropy. The rheology of the two media, e.g. the matrix and the fiber network, as
well as their interaction follow non-linear viscous behaviors. The equations of this model are
simplified to the case of the compression of SMC, giving the formulation of a shell model whose
equations are written into a finite element code. Simple simulation examples thus show the
strong influence of material and process parameters on the phenomenon of phase separation.

RÉSUMÉ. Un formalisme biphasique est proposé pour modéliser la rhéologie des composites fibres de verre - matrice polymère de type SMC ou GMT au cours de leur mise en forme. L'ani-sotropie liée à la morphologie du réseau de fibres est simplifiée au cas de l'orthotropie de révolution. Le comportement visqueux non linéaire régit la rhéologie des deux milieux consi-dérés (matrice et réseau de fibres) ainsi que leur interaction. Les équations de ce modèle sont simplifiées et adaptées au cas du procédé de compression des SMC ; ceci aboutit à la formu-lation d'un modèle coque, dont les équations sont introduites dans un code éléments finis. Des exemples simples de simulations permettent alors de souligner l'importance des paramètres « matériau » et « procédé » sur le phénomène de redistribution de phases.

KEYWORDS: SMC, GMT, theory of mixture, two-phase shell model, segregation.

MOTS-CLÉS : SMC, GMT, théorie des mélanges, modèle biphasique coque, ségrégation.
# 1. Introduction

Glass fiber composites processed by compression molding are increasingly used in the automotive and electrical industries as semi-structural parts because of their lightweight and their cost-efficient processing. These composites are formed by a mat of glass fibers or glass fiber bundles, which constitutes the reinforcement of the matrix. The fibers or bundles content is high so that the composites can be seen as highly concentrated suspensions. Their matrix may be either thermoplastic or a thermoset polymer. Thermoplastic based composites are called GMT for Glass Mat Thermoplastics, whereas thermoset based composites are called SMC for Sheet Molding Compounds. All these composites are produced as planar preforms before being processed by compression molding to produce parts. During this compression stage, the initial distributions of fiber orientation and fiber content through the produced parts is not maintained, strongly affecting final geometry and properties of produced parts [OSS 94, YAG 95, THO 96]. These two kinds of phenomena have to be controlled in order to ensure suitable mechanical and geometric properties of the produced parts. Current rheological models entering in the scope of models for semi-dilute or concentrated suspensions can predict fiber orientation evolution during processing of polymer composites [ADV 94]. However, as they assume that the fiber network and the matrix have the same macroscopic velocity, these one-phase models cannot account for fiber segregation, keeping a homogeneous fiber content during the whole compression molding. Therefore, to describe the evolution of the fiber content in SMC or GMT, a two-phase model based on the mixture theory is proposed in this paper. Firstly, the paper briefly shows how the two-phase model is established. Secondly, the model is simplified to the case of compression molding: a formulation for a two-phase "Barone and Caulk" like shell model is proposed. Thirdly, the finite element method for solving the two-phase boundary value problem is formulated. The paper is concluded by few numerical examples where the influence of processing parameters on segregation phenomenon is analyzed.

# 2. Balance equations for a two-phase mixture

The theory of mixture that was initially developed in the pioneering works of Truesdell and Toupin [TRU 60, BOW 76] is a general framework well suited to predict fiber segregation phenomena arising during compression molding. Within such a theoretical framework, the following basic assumptions were stated:

– the composite is seen as the superposition of two continuous media. Each continuous medium represents immiscible constituent or phase of the material, i.e. the fiber network $f$ and the matrix $m$ (polymer + filler). Hence, each material point $M$ of the mixture is simultaneously occupied by material points $M^{(f)}$ and $M^{(m)}$ of the phases $f$ and $m$, respectively<fnref n="1" />.

<footnote n="1">1. In order to be able to distinguish phenomena arising at a microscopic scale (inside the matrix or in the fiber in our case) from those "apparently" resulting at a macroscopic scale, a specific</footnote>
– Each elementary macroscopic volume element $\delta V$ of the mixture (elementary mass $\delta m$ and specific mass $\rho = \frac{\delta m}{\delta V}$) is simultaneously occupied by the phases $\alpha$ whose elementary macroscopic mass $\delta m^{(\alpha)}$ occupies an elementary macroscopic volume $\delta V_{(\alpha)}$ of $\delta V$. This enables to introduce (i) the macroscopic $\rho^{(\alpha)}$ and the microscopic $\rho_{(\alpha)}$ specific masses:

$$ \delta m^{(\alpha)} = \delta V \rho^{(\alpha)} = \int_{\delta V_{(\alpha)}} \rho_{(\alpha)} dV, \qquad (1) $$

and (ii) the volume fraction $f^{(\alpha)}$ of the phase $\alpha$:

$$ f^{(\alpha)} = \frac{\delta V_{(\alpha)}}{\delta V} = \frac{\rho^{(\alpha)}}{\rho}. \qquad (2) $$

In this work, the constituents of the mixture are assumed incompressible at the microscopic scale so that $\rho_{(\alpha)}$ is constant and $\rho^{(\alpha)} = f^{(\alpha)}\rho_{(\alpha)}$.

– The saturation condition for the mixture is stated:

$$ f^{(m)} + f^{(f)} = 1. \qquad (3) $$

– At last, to simplify the analysis, only isothermal situations will be considered in this paper. Hence, only mass and momentum balances have to be considered in the modelling (see next subsection).

## 2.1. Mass and momentum balance equations

Accounting for the previous assumptions the local mass balance equation written at a given material point $M^{(\alpha)}$ for each phase $\alpha$ is

$$ \frac{\mathrm{D}^{(\alpha)}}{\mathrm{D}t} f^{(\alpha)} + f^{(\alpha)} \mathrm{div} \underline{v}^{(\alpha)} = 0, \qquad (4) $$

where the notation $\frac{\mathrm{D}^{(\alpha)}}{\mathrm{D}t}$ is the material time derivative following materials points $M^{(\alpha)}$ of velocity $\underline{v}^{(\alpha)}$. Summing equations (4) written for each phase of the mixture, i.e. $m$ and $f$, and accounting for the saturation of the mixture (3) leads to an incompressibility equation for the mixture:

$$ \mathrm{div}(f^{(f)} \underline{v}^{(f)} + (1 - f^{(f)}) \underline{v}^{(m)}) = 0. \qquad (5) $$

notation for microscopic and macroscopic quantities is introduced: if $\phi$ is a scalar physical quantity, then $\phi_{(\alpha)}$ (with the subscript $(\alpha)$) is the microscopic (or local) value of $\phi$ in the phase $\alpha$, whereas $\phi^{(\alpha)}$ (with the exponent $(\alpha)$) is the macroscopic (or apparent) value of $\phi$ for the same phase in the mixture.
The mixture theory defines the concept of partial stress vector $\underline{T}^{(\alpha)}$ and partial stress tensor $\underline{\sigma}^{(\alpha)}$ for each phase $\alpha$, writing the total force exerted from the outside on the surface element $\delta S$ of $\delta V$ as :

$$ \underline{\delta F} = \underline{\delta F}^{(m)} + \underline{\delta F}^{(f)} = (\underline{T}^{(m)} + \underline{T}^{(f)})\delta S = (\underline{\sigma}^{(m)} + \underline{\sigma}^{(f)}) \cdot \underline{n}\delta S, \quad (6) $$

where $\underline{n}$ is an outward normal vector to $\delta S$. Assuming no external and internal specific moments, the second momentum balance equations written for each phase lead to the symmetry of the partial stress tensor, i.e. $\underline{\sigma}^{(\alpha)} = \underline{\sigma}^{T(\alpha)}$. Moreover, rheometry experiments performed on industrial SMC or GMT clearly reveal that during compression molding these composites can reasonably be approximated as purely viscous and incompressible anisotropic fluids [SER 02, LEC 02, DUM 03]. In this case, it can be shown that the partial stress tensors $\underline{\sigma}^{(\alpha)}$ can be split into the sum of two terms [BOW 76]:

$$ \underline{\sigma}^{(\alpha)} = -f^{(\alpha)}p_i\underline{\delta} + \underline{\sigma}^{e(\alpha)}, \quad (7) $$

where $p_i$ is a pressure related to the incompressibility of the mixture, $\underline{\sigma}^{e(m)}$ and $\underline{\sigma}^{e(f)}$ are "extra" or viscous partial stresses. In this case, neglecting inertial effects and specific external forces, it can be shown that the local form of the first global momentum balance equations for the matrix and the network of fibers are [BOW 76]:

$$ \underline{\text{div}} \underline{\sigma}^{e(m)} - (1 - f^{(f)}) \underline{\text{grad}} p_i + \hat{\underline{p}}^{e(m)} = \underline{0}, \quad (8) $$

$$ \underline{\text{div}} \underline{\sigma}^{e(f)} - f^{(f)} \underline{\text{grad}} p_i - \hat{\underline{p}}^{e(f)} = \underline{0}. \quad (9) $$

where $\hat{\underline{p}}^{e(m)}$ is a viscous momentum exchange between the two phases.

# 3. Towards a two-phase shell model for compression

## 3.1. Experimental observations and consequences

The deformation of composites during industrial compression molding is usually described using two "elementary" mechanisms. The first deformation mode is a "squeeze flow", characterized by the shearing of the composite in the thickness of the molded sheets. It is thus assumed that there is a sticking contact between the composite and the upper and lower parts of the mold: such a mechanism has been used to establish one-phase generalized Hele-Shaw shell models [HIE 80, LEE 84, OSS 90]. The second mechanism is a "plug flow" characterized by an uniform deformation in the thickness of the sheets, with a slipping contact at the interface between the composite and the mold surfaces: this second type of flow has given rise to one-phase Barone-Caulk shell models [BAR 86, BAR 87]. In practice, the predominance of the first or second mechanism depends on processing conditions, thickness of the sheets, length of fibers... In the case of SMC or GMT, where fibers are much longer than the thickness of the sheets, experimental evidences show that within a wide range of processing conditions, the plug flow regime leads to a rather good approximation of the
real flow patterns [BAR 85, SER 02, ODE 04]. Consequently, in order to simplify the full 3D two-phase model introduced in the previous section, a two-phase shell model is proposed in this section, adopting a plug flow kinematics for the composite: it is a direct extension of the one-phase shell model of Barone and Caulk [BAR 86]. Note that for a sake of simplicity, this approach is here restricted to the case of a plate geometry whose midplane is contained in the $(\underline{e}_1, \underline{e}_2)$ plane: the extension to a 3D shell model would be straightforward. Hence, the velocity $\underline{v}^{(\alpha)}$ of each phase becomes:

$$ \underline{v}^{(\alpha)} = v_1^{(\alpha)}(x_1, x_2)\underline{e}_1 + v_2^{(\alpha)}(x_1, x_2)\underline{e}_2 + \frac{\dot{h}}{h}x_3\underline{e}_3 = \tilde{\underline{v}}^{(\alpha)} + \frac{\dot{h}}{h}x_3\underline{e}_3 \quad (10) $$

where the symbol "~" has been introduced to distinguish 2D-unknown fields and 2D-operators in $(\underline{e}_1, \underline{e}_2)$ from 3D ones, and $h$ is the thickness of the sheets. As the shear components $D_{\beta 3}^{(\alpha)}$ ($\beta \in \{1, 2\}$) of the strain rate tensor are constrained to zero, arbitrary reaction terms must be added to the partial stress tensor $\underline{\underline{\sigma}}^{(\alpha)}$:

$$ \underline{\underline{\sigma}}^{(\alpha)} = -f^{(\alpha)}p_i\underline{\underline{\delta}} + T_{\beta 3}^{(\alpha)}(\underline{e}_{\beta} \otimes \underline{e}_3 + \underline{e}_3 \otimes \underline{e}_{\beta}) + \underline{\underline{\sigma}}^{e(\alpha)}, \quad \beta \in \{1, 2\}. \quad (11) $$

Hence, the mass balance equation for the fiber network becomes:

$$ \frac{D^{(f)} f^{(f)}}{Dt} + f^{(f)} \left( \mathrm{div} \tilde{\underline{v}}^{(f)} + \frac{\dot{h}}{h} \right) = 0. \quad (12) $$

Moreover, the incompressibility of the mixture now reads:

$$ \mathrm{div}(f^{(f)} \tilde{\underline{v}}^{(f)} + (1 - f^{(f)}) \tilde{\underline{v}}^{(m)}) + \frac{\dot{h}}{h} = 0. \quad (13) $$

Finally, the two momentum balance equations in the plane of the sheets are integrated over the thickness $h$ of the composite. In the case of isothermal situations they become:

$$ -(1 - f^{(f)}) \underline{\underline{\mathbf{grad}}} p_i + \frac{1}{h} \left( \underline{\underline{\tilde{F}}}_h^{(m)} + \underline{\underline{\tilde{F}}}_0^{(m)} \right) + \underline{\underline{\mathbf{div}}} \tilde{\underline{\underline{\sigma}}}^{e(p)} + \tilde{\underline{\underline{\hat{p}}}}^{e(m)} = \underline{\underline{\mathbf{0}}} \quad (14) $$

and

$$ -f^{(f)} \underline{\underline{\mathbf{grad}}} p_i + \frac{1}{h} \left( \underline{\underline{\tilde{F}}}_h^{(f)} + \underline{\underline{\tilde{F}}}_0^{(f)} \right) + \underline{\underline{\mathbf{div}}} \tilde{\underline{\underline{\sigma}}}^{e(f)} - \tilde{\underline{\underline{\hat{p}}}}^{e(m)} = \underline{\underline{\mathbf{0}}}, \quad (15) $$

where

$$ \left( \underline{\underline{\tilde{F}}}_h^{(\alpha)} + \underline{\underline{\tilde{F}}}_0^{(\alpha)} \right) \cdot \underline{\underline{e}}_{\beta} = \left[ T_{\beta 3}^{(\alpha)} \right]_0^h \quad \text{with} \quad \beta \in \{1, 2\}. $$

$\tilde{\underline{F}}_0^{(\alpha)}$ and $\tilde{\underline{F}}_h^{(\alpha)}$ account for possible friction effects that are generated by the slipping on the upper ($x_3 = 0$) and lower ($x_3 = h$) surfaces of the mold parts [BAR 86] (cf. next subsection).

## 3.2. Constitutive equations

– It seems reasonable to consider that the polymer matrix behaves like a power-law fluid at the microscopic scale, with a power-law exponent $n_{(m)}$ [LEC 02, DUM 03].
Following this idea, the extra macroscopic stress tensor of the matrix is described by a simple constitutive law:

$$ \tilde{\underline{\underline{\sigma}}}^{e(m)} = 2\mu^{(m)} \left( \frac{\dot{\gamma}_{eq}^{(m)}}{\dot{\gamma}_0} \right)^{n_{(m)}-1} \tilde{\underline{\underline{D}}}^{(m)}, \quad \sigma_{33}^{e(m)} = 2\mu^{(m)} \left( \frac{\dot{\gamma}_{eq}^{(m)}}{\dot{\gamma}_0} \right)^{n_{(m)}-1} \frac{\dot{h}}{h} \quad (16) $$

with $\mu^{(m)} = (1 - f^{(f)})\mu_{0(m)}$, $\mu_{0(m)}$ being the viscosity of the polymeric matrix at a reference shear rate $\dot{\gamma}_0$ of 1 s⁻¹. Accounting for the current kinematical restrictions,

$$ \dot{\gamma}_{eq}^{(m) 2} = 2 \left( \tilde{\underline{\underline{D}}}^{(m)} : \tilde{\underline{\underline{D}}}^{(m)} + \left( \frac{\dot{h}}{h} \right)^2 \right) \quad (17) $$

– The extra stress tensor $\tilde{\underline{\underline{\sigma}}}^{e(f)}$ of the fiber network is a direct extension of the one-phase mechanical model proposed by Dumont et al. [DUM 03]. Hence, the fiber network is seen as a compressible power-law fluid displaying a transverse isotropy whose axis is $\underline{e}_3$:

$$ \tilde{\underline{\underline{\sigma}}}^{e(f)} = \eta^{(f)} \frac{2}{1+2H} \left( \frac{D_e}{D_0} \right)^{n-1} \tilde{\underline{\underline{D}}}, \quad \sigma_{33}^{e(f)} = \eta^{(f)} \frac{2H}{1+2H} \left( \frac{D_e}{D_0} \right)^{n-1} \frac{\dot{h}}{h}, \quad (18) $$

where

$$ D_e^{(f) 2} = \frac{2}{1+2H} \left( \tilde{\underline{\underline{D}}}^{(f)} : \tilde{\underline{\underline{D}}}^{(f)} + H \left( \frac{\dot{h}}{h} \right)^2 \right), \quad \eta^{(f)} = \eta_{ps}^{(f)} \left( 2 \frac{1+H}{1+2H} \right)^{-\frac{n+1}{2}}, \quad (19) $$

with $n_{(f)}$ the power-law exponent, $\eta_{ps}^{(f)}$ a plane strain compression viscosity at a characteristic strain rate $D_0$ of 1 s⁻¹, and $H$ a rheological function that accounts for the fiber network anisotropy. $\eta_{ps}^{(f)}$ and $H$ strongly depend on $f^{(f)}$. For example, using the experimental results obtained on a standard SMC formulation [DUM 03], possible expressions for $\eta_{ps}^{(f)}$ and $H$ are

$$ \eta_{ps}^{(f)} = \eta_{(f)} (98f^{(f)} + 980f^{(f)2}), \quad H = \frac{1 + 98f^{(f)} + 980f^{(f)2}}{0.5 + 67f^{(f)} + 670f^{(f)2}} - 1, \quad (20) $$

where $\eta_{(f)}$ is a constant.

– The momentum exchange $\tilde{\underline{\underline{p}}}^{e(m)}$ is induced when there is a relative motion $\tilde{\underline{\underline{v}}}^r = \tilde{\underline{\underline{v}}}^{(m)} - \tilde{\underline{\underline{v}}}^{(f)}$ of the matrix with respect to the fiber network. The expression of this term is such that the momentum balance for the matrix phase (8) reduces to the permeation law of a power-law fluid trough a rigid network of fibers or bundles or fibers (in the latter case, no permeation inside the bundles is assumed), when the matrix does not deform at the macroscopic scale. A possible momentum exchange is [AUR 02]:

$$ \tilde{\underline{\underline{p}}}^{e(m)} = -(1 - f^{(f)}) \frac{2^{n+1} \mu_{0(m)}}{\kappa^* k_{sq}^* (f^{(f)}, n_{(m)}) d^2} \left( \frac{\|\tilde{\underline{\underline{v}}}^r\|}{\dot{\gamma}_0 d} \right)^{n_{(m)}-1} \tilde{\underline{\underline{v}}}^r, \quad (21) $$
where d is the average diameter of the fibers (or bundles), $k_{sq}^*(f^{(f)}, n_{(m)})$ is the dimensionless permeability for transverse flow of a power-law fluid through a regular square array of parallel fibers of diameter d and fiber volume fraction $f^{(f)}$, and where $\kappa^*$ is a constitutive parameter to determine. According to the experimental data collected for various fibrous media: $0.1 < \kappa^* < 100$ [JAC 86]. It can be shown that $k_{sq}^*(f^{(f)}, n_{(m)})$ can be directly deduced from $k_{sq}^*(f^{(f)}, 1)$ [BRU 93, IDR 04]:

$$k_{sq}^*(f^{(f)}, n_{(m)}) = \left( \frac{\sqrt{\Phi}}{0.8(1 - \sqrt{\Phi})^2} \right)^{1-n_{(m)}} k_{sq}^*(f^{(f)}, 1), \quad (22)$$

where $f_{max}^{(f)} = \pi/4$, $\Phi = f^{(f)}/f_{max}^{(f)}$. In this work, $k_{sq}^*(f^{(f)}, 1)$ was approximated by the analytical lubrication model [KEL 64], that gives rather good predictions in a wide range of fiber volume fraction:

$$k_{sq}^*(f^{(f)}, 1) = \frac{(1 - \Phi)^2}{3\sqrt{\Phi}^3} \left( 3\sqrt{\Phi} \frac{\arctan\left(\sqrt{\frac{1+\sqrt{\Phi}}{1-\sqrt{\Phi}}}\right)}{\sqrt{1-\Phi}} + \frac{1}{2}\Phi + 1 \right)^{-1}. \quad (23)$$

– In their one-phase shell-model, Barone et Caulk [BAR 86] have given a physical meaning for the forces $\tilde{F}$ appearing in (14-15): they are related to the presence of a thin amount of viscous matrix entrapped between the compressed sheets and the mold, that is sheared during the relative motion of the sheets with respect to the mold. Adopting a similar reasoning for the two-phase model, we propose to write:

$$\tilde{\underline{F}}_h^{(\alpha)} = \tilde{\underline{F}}_0^{(\alpha)} = -f^{(\alpha)} \mu_{0(m)} \kappa_H^{(\alpha)} \left( \frac{\tilde{\underline{v}}^{(\alpha)} \cdot \tilde{\underline{v}}^{(\alpha)}}{v_0^2} \right)^{\frac{n_{(m)}-1}{2}} \tilde{\underline{v}}^{(\alpha)}, \quad (24)$$

where $v_0$ is a reference velocity and $\kappa_H^{(\alpha)}$ are friction parameters to determine.

### 3.3. Initial and boundary conditions

At the beginning of the compression, i.e. for $t = 0$, the composite occupies a surface $\Omega(x_1, x_2, t = 0) = \Omega_0$ of boundary $\partial\Omega(x_1, x_2, 0) = \partial\Omega_0$ in the principal plane ($\underline{e}_1, \underline{e}_2$) of the mold whose total surface and boundary are $\Omega_M$ ($\Omega_0 \subset \Omega_M$) and $\partial\Omega_M$, respectively. The initial height $h(x_1, x_2, 0) = h_0$ and fiber volume fraction $f^{(f)}(x_1, x_2, 0) = f_0^{(f)}$ are given, and the initial velocity of each phase is zero.

During the compression, i.e. for $t > 0$, the local thickness $h(x_1, x_2, t)$ is imposed, so that the composite of surface $\Omega(x_1, x_2, t)$ of boundary $\partial\Omega(x_1, x_2, t)$ fills the mold cavity, accounting for the following set of boundary conditions:

- $\tilde{\underline{v}}^{(\alpha)} \cdot \tilde{\underline{n}} = 0$ on $\partial_v \Omega = \partial\Omega \cap \partial\Omega_M$, where $\tilde{\underline{n}}$ is the unit outward normal to $\Omega$ and $\partial_v \Omega^{(\alpha)}$ is the part of $\partial\Omega$ which is in contact with the mold boundary $\partial\Omega_M$;

- $\tilde{\underline{\sigma}}^{(\alpha)} \cdot \tilde{\underline{n}} = \underline{0}$ on $\partial_\sigma \Omega = \partial\Omega - \partial_v \Omega$, where $\partial_\sigma \Omega$ stands for the free surface of the mixture.
# 4. FE formulation of the two-phase boundary value problem

## 4.1. Resolution scheme

The flow of the two-phase mixture is described by an Eulerian approach with respect to a reference surface containing ($e_1, e_2$). In the present application, this reference is the mid-surface of the mold $\Omega_M$ that is meshed using triangular elements, each of them having a specific height $h$. During compression molding, the composite fills the mold cavity so that $\Omega$ grows from $\Omega_0$ to $\Omega_M$ and the evolution of the free surface $\partial_\sigma\Omega$ has to be determined. As in many mold filling applications [DHA 92, SCA 99, SOU 01], this was achieved with an additional scalar variable $\chi$ describing the local volume fraction of the composite in the mold (i.e. $\chi = 0/1$ if the considered material point is empty/full of composite). The dynamic of the variable $\chi(x_1, x_2, t)$ is ruled by the following equation:

$$ \frac{D^{(\text{mix})}\chi}{Dt} = 0, \qquad (25) $$

where $\frac{D^{(\text{mix})}}{Dt}$ stands for the material time derivative following the mixture of velocity $\tilde{v}^{(\text{mix})}$. As defined in [BOW 76]:

$$ \tilde{v}^{(\text{mix})} = \frac{1}{f^{(f)}\rho_{(f)} + (1 - f^{(f)})\rho_{(m)}} \left( f^{(f)}\tilde{v}^{(f)} + (1 - f^{(f)})\tilde{v}^{(m)} \right). \qquad (26) $$

Thus, the scalar field $\chi$ is added to the unknown fields $f^{(f)}, \tilde{v}^{(m)}, \tilde{v}^{(f)}$ and $p_i$. All are dependent of time $t$ and space variables $x_1$ and $x_2$. An usual strategy is employed to solve this problem. It consists in splitting the time and the spatial discretizations as in many usual FE treatments. Thereby, the time interval $[0, T]$ of molding is subdivided into a finite number of time steps $]t^n, t^{n+1}]$ that can have different lengths $^2$. Given the solution at time $t^n$, the step $n + 1$ consists in (i) finding the new domain $\Omega^{n+1}$ occupied by the mixture, and (ii) determining the unknown fields $v^{(m)^{n+1}}, v^{(f)^{n+1}}, p_i^{n+1}, f^{(f)^{n+1}}$ and $\chi^{n+1}$ in $\Omega^{n+1}$. These two points are detailed below.

(i) The calculation of the new domain $\Omega^{n+1}$ from the knowledge of $\Omega^n, v^{(\text{mix})}^n$ and $\chi^n$ requires (a) the computation of the time increment $\Delta t = \Delta t^{n \to n+1}$ (it is determined such as the flow front progression is of the order of a typical element dimension), (b) the determination of $t^{n+1}$ such as $t^{n+1} = t^n + \Delta t$, (c) the calculation of $h^{n+1}$ and $\dot{h}^{n+1}$ for each element of the mesh, and (d) the determination of the elements that are gained by the mixture between $t^n$ and $t^{n+1}$.

(ii) The two-phase problem is then subdivided into three sub-problems:

- Sub-problem (SP1), called the pressure-velocities problem. It aims at finding $v^{(m)^{n+1}}, v^{(f)^{n+1}}$ and $p_i^{n+1}$, and couples the incompressibility equation (13) and the

<footnote n="2">2. In the following, we will find the following notation for every function $f$: $f(\underline{x}, t^n) = f^n(\underline{x})$.</footnote>
two momentum balance equations (14-15). This first sub-problem would be equiv-
alent to the well-known Stokes problem in the case of a one-phase approach with a
Newtonian incompressible fluid.

- Sub-problem (SP2), called the fiber fraction problem. It consists in solving the mass balance equation (12) to determine $f^{(f)^{n+1}}$. Notice that the resolution of SP1 and SP2 is coupled: indeed, at $t^{n+1}$, the pressure-velocities problem SP1 and the fiber volume fraction evolution problem SP2 are successively solved up to convergence, i.e. up to finding a fixed point for the velocities $\tilde{\underline{v}}^{(m)^{n+1}}$ and $\tilde{\underline{v}}^{(f)^{n+1}}$, pressure $p_i^{n+1}$ and fiber volume fraction fields $f^{(f)^{n+1}}$.

- Sub-problem (SP3), called the free surface problem. It consists in solving the transport equation (25) on $\Omega^{n+1}$ to compute $\chi^{n+1}$.

In what follows, the different numerical methods used to solve these problems are presented.

## 4.2. FE formulation of the two-phase pressure-velocities problem

### 4.2.1. Weak formulation

The pressure-velocities problem SP1 is re-casted into a weak form by multiply-
ing the equations (14-15 and 13) by a set of test functions $\tilde{\underline{v}}^{(m)*}$, $\tilde{\underline{v}}^{(f)*}$ and $p^*$ (be-
longing to appropriate spaces $V_0^{(m)}$, $V_0^{(f)}$ and $Q$) and by integrating over $\Omega$ apply-
ing the divergence theorem. The problem consists now in finding the set of un-
knowns functions $(\tilde{\underline{v}}^{(m)}, \tilde{\underline{v}}^{(f)}, p_i) \in (V^{(m)} \times V^{(f)} \times Q)$ such as $\forall (\tilde{\underline{v}}^{(m)*}, \tilde{\underline{v}}^{(f)*}, p^*) \in (V^{(m)} \times V^{(f)} \times Q)$:

$$
\left\{
\begin{aligned}
\int_{\Omega} \tilde{\underline{\sigma}}^{e(m)} : \underline{\underline{D}}^{(m)*} d\Omega - \int_{\Omega} \frac{1}{h} (\underline{\underline{F}}_h^{(m)} + \underline{\underline{F}}_0^{(m)}).\underline{\underline{v}}^{(m)*} d\Omega - \int_{\Omega} \tilde{\underline{\underline{p}}}^{e(m)}.\underline{\underline{v}}^{(m)*} d\Omega - \dots \\
\qquad \qquad \qquad \qquad \int_{\Omega} \mathrm{div} \left[ (1-f^{(f)})\underline{\underline{v}}^{(m)*} \right] p_i d\Omega = \int_{(\partial_v \Omega)} \underline{\underline{v}}^{(m)*} \cdot \underline{\underline{\sigma}}^{e(m)}. \underline{\underline{n}} d(\partial_v \Omega), \\[2em]
\int_{\Omega} \underline{\underline{\sigma}}^{e(f)} : \underline{\underline{D}}^{(f)*} d\Omega - \int_{\Omega} \frac{1}{h} (\underline{\underline{F}}_h^{(f)} + \underline{\underline{F}}_0^{(f)}).\underline{\underline{v}}^{(f)*} d\Omega + \int_{\Omega} \tilde{\underline{\underline{p}}}^{e(m)}.\underline{\underline{v}}^{(f)*} d\Omega - \dots \\
\qquad \qquad \qquad \qquad \int_{\Omega} \mathrm{div} \left[ f^{(f)} \underline{\underline{v}}^{(f)*} \right] p_i d\Omega = \int_{(\partial_v \Omega)} \underline{\underline{v}}^{(f)*} \cdot \underline{\underline{\sigma}}^{e(f)}. \underline{\underline{n}} d(\partial_v \Omega), \\[2em]
\int_{\Omega} \mathrm{div} \left[ f^{(f)} \underline{\underline{v}}^{(f)} \right] p^* d\Omega + \int_{\Omega} \mathrm{div} \left[ (1-f^{(f)})\underline{\underline{v}}^{(m)} \right] p^* d\Omega = - \int_{\Omega} \frac{\dot{h}}{h} p^* d\Omega,
\end{aligned}
\right.
\tag{27}
$$

It is important to underline that this system is highly non-linear due to the power-laws used to determine $\underline{\underline{\sigma}}^{e(\alpha)}$, $\tilde{\underline{\underline{p}}}^{e(m)}$ and $\underline{\underline{F}}_h^{(\alpha)} + \underline{\underline{F}}_0^{(\alpha)}$ (see subsection 3.2).
### 4.2.2. FE approximations

A mixed finite element method is used to discretize the previous weak formulation. The considered elements are P2+ for the two velocity fields and P1 for the in-thickness integrated pressure. The interpolation of the velocity fields is quadratic (six degrees of freedom plus an internal velocity node) whereas the interpolation of the pressure is piecewise linear. The use of these elements allows to circumvent the Brezzi-Babuška compatibility condition. The following matrix formulation of the pressure-velocities problem can be readily obtained via the standard Galerkin discretization process:

$$ \left( \begin{array}{cc} \mathbf{A}(\tilde{\mathbf{V}}) & {}^{t}\mathbf{B} \\ \mathbf{B} & \mathbf{0} \end{array} \right) \left( \begin{array}{c} \tilde{\mathbf{V}} \\ \mathbf{P} \end{array} \right) = \left( \begin{array}{c} \mathbf{F}(\tilde{\mathbf{V}}) \\ \mathbf{G} \end{array} \right) \qquad (28) $$

where $\tilde{\mathbf{V}}$ is the vector of nodal velocities unknowns (m and f), i.e.

$$ \tilde{\mathbf{V}} = \left( \frac{\tilde{\mathbf{V}}^{(m)}}{\tilde{\mathbf{V}}^{(f)}} \right), \qquad (29) $$

P stands for the vector of integrated pressure nodal unknowns, F is a vector that accounts for the boundary conditions and where G contains terms arising from the incompressibility constraint associated to $\dot{h}/h$. More precisely, one can write:

$$ \left( \begin{array}{cc} \mathbf{A}(\tilde{\mathbf{V}}) & {}^{t}\mathbf{B} \\ \mathbf{B} & \mathbf{0} \end{array} \right) = \dots $$

$$ \left( \begin{array}{c|c|c} \mathbf{A}_1^{(m)} + \mathbf{A}_2^{(m)} + \mathbf{A}_3^{(f/m)} & -\mathbf{A}_3^{(f/m)} & {}^t\mathbf{B}^{(m)} \\ \hline -\mathbf{A}_3^{(f/m)} & \mathbf{A}_1^{(f)} + \mathbf{A}_2^{(f)} + \mathbf{A}_3^{(f/m)} & {}^t\mathbf{B}^{(f)} \\ \hline \mathbf{B}^{(m)} & \mathbf{B}^{(f)} & \mathbf{0} \end{array} \right) \qquad (30) $$

where the matrices $\mathbf{A}_1^{(\alpha)}$, $\mathbf{A}_2^{(\alpha)}$ and $\mathbf{A}_3^{(f/m)}$ are "viscous" matrices arising respectively from the rheology of the phases, the friction terms, and the momentum exchange, and where $\mathbf{B}^{(\alpha)}$ stands for the incompressibility condition.

In the case of a linear system ($n_{(\alpha)} = 1$), the problem can be solved by applying the Uzawa algorithm for decoupling the momentum equations and the incompressibility constraint. To reduce the number of iterations up to convergence, the Uzawa algorithm was modified by solving the saddle-point problem for the augmented Lagrangian function of the system (28) defined as

$$ \mathcal{L}_r(\tilde{\mathbf{U}}, \mathbf{Q}) = \frac{1}{2} {}^t\tilde{\mathbf{U}}\mathbf{A}\tilde{\mathbf{U}} + \frac{1}{2} r^t(\mathbf{B}\tilde{\mathbf{U}} - \mathbf{G})(\mathbf{B}\tilde{\mathbf{U}} - \mathbf{G}) - {}^t\tilde{\mathbf{U}}\mathbf{F} + {}^t(\mathbf{B}\tilde{\mathbf{U}} - \mathbf{G})\mathbf{Q}, \quad (31) $$
where $r$ is a constant [ROB 90, LIU 01]. Notice that the problem (31) can be generally written into a linear system of equations:

$$ \left( \begin{array}{cc} \mathbf{A} + r^t \mathbf{B} \mathbf{B} & \mathbf{t} \mathbf{B} \\ \mathbf{B} & \mathbf{0} \end{array} \right) \left( \begin{array}{c} \tilde{\mathbf{V}} \\ \mathbf{P} \end{array} \right) = \left( \begin{array}{c} \mathbf{F} + r^t \mathbf{B} \mathbf{G} \\ \mathbf{G} \end{array} \right). \quad (32) $$

If the problem is non-linear, a Newton-Raphson method is used. To find the solution of this linearized system, we apply the modified Uzawa algorithm used to solve the linear problems. The general algorithm for the non-linear pressure-velocities problem combines the Newton-Raphson and Uzawa algorithms. Given $r$, $nr_{max}$, $\varepsilon_{NR}$, $\rho$, $k_{max}$, $\varepsilon_{UZ}$, an arbitrary choice of $P_{nr=1}$ for $P$ and $\tilde{V}_{nr=1}$ for $\tilde{V}$:

Begin Newton-Raphson

repeat

Compute $\Delta P_{nr}$, $\Delta V_{nr}$:

$$ \Delta P_{nr}^{k=1} \leftarrow 0 $$

Begin Uzawa

repeat

Compute $\Delta \tilde{V}_{nr}^k$ by solving the following equation:

$$ \begin{aligned} & \left( \frac{\partial \mathbf{A}(\tilde{\mathbf{V}}_{nr})}{\partial \tilde{\mathbf{V}}} + \mathbf{A}(\tilde{\mathbf{V}}_{nr}) + r^t \mathbf{B} \mathbf{B} - \frac{\partial F(\tilde{\mathbf{V}}_{nr})}{\partial \tilde{\mathbf{V}}} \right) \Delta \tilde{\mathbf{V}}_{nr}^k \\ & = -(\mathbf{A}(\tilde{\mathbf{V}}_{nr}) + r^t \mathbf{B} \mathbf{B}) \tilde{\mathbf{V}}_{nr} - \mathbf{t} \mathbf{B} P_{nr} + F(\tilde{\mathbf{V}}_{nr}) + r^t \mathbf{B} \mathbf{G} - \mathbf{t} \mathbf{B} \Delta P_{nr}^{k-1} \end{aligned} $$

Then compute $\Delta P_{nr}^k$:

$$ \begin{aligned} \Delta P_{nr}^k & \leftarrow \Delta P_{nr}^{k-1} + \rho (\mathbf{B}(\tilde{\mathbf{V}}_{nr} + \Delta \tilde{\mathbf{V}}_{nr}^k) - \mathbf{G}) \\ k & \leftarrow k + 1 \end{aligned} $$

until ($\|\mathbf{B}(\tilde{\mathbf{V}}_{nr} + \Delta \tilde{\mathbf{V}}_{nr}^k) - \mathbf{G}\| \le \varepsilon_{UZ}$ or $k = k_{max}$)

End Uzawa

$$ \tilde{\mathbf{V}}_{nr+1} \leftarrow \tilde{\mathbf{V}}_{nr} + \Delta \tilde{\mathbf{V}}_{nr}^k $$

$$ P_{nr+1} \leftarrow P_{nr} + \Delta P_{nr}^k $$

$$ nr \leftarrow nr + 1 $$

until ($\|\mathbf{R}(\tilde{\mathbf{V}}_{nr}, P_{nr})\| \le \varepsilon_{NR}$ or $nr = nr_{max}$)

End Newton-Raphson

$$ \tilde{\mathbf{V}}^{n+1} \leftarrow \tilde{\mathbf{V}}_{nr} $$

$$ P^{n+1} \leftarrow P_{nr} $$

where $nr_{max}$, $\varepsilon_{NR}$ are respectively the maximum number of iterations and the convergence criterion, where $k_{max}$ and $\varepsilon_{UZ}$ are respectively the maximum number of iterations and the convergence criterion of the Uzawa algorithm, and where $0 < \rho < 2r$ is the Uzawa parameter. In practice, the components of the matrices and vectors of the previous system are computed using a standard Gaussian quadrature integration method [DHA 81, ZIE 00]. The solver of the matrix system given by combining the Newton-Raphson and Uzawa algorithms is based on a biconjugate gradient algorithm.
### 4.3. Mass balance equation for the fiber phase and free surface transport

– The solution of the fiber volume fraction evolution equation (12) is obtained trough a finite element discretization in space and a characteristic based method discretization in time. The variational formulation of the equation (12) reads $\forall \phi^* \in X$:

$$ \int_{\Omega} \frac{D^{(f)} f^{(f)^{n+1}}}{Dt} \phi^* d\Omega + \int_{\Omega} \left( \tilde{\mathrm{div}} \tilde{\underline{v}}^{(f)^{n+1}} + \frac{\dot{h}}{h} \right) f^{(f)^{n+1}} \phi^* d\Omega = 0, \quad (33) $$

where $f^{(f)^{n+1}}$ is the solution of the fiber volume fraction field at time $t^{n+1}$ and $\phi^*$ a test function defined in an appropriate space $X$. From (33), the characteristic method is applied to approximate the material derivative $\frac{D^{(f)}}{Dt} f^{(f)^{n+1}}$. It consists in writing the following finite difference scheme:

$$ \left\{ \begin{array}{l} \displaystyle \frac{D^{(f)} f^{(f)^{n+1}}}{Dt} = \frac{f^{(f)^{n+1}}(\tilde{\underline{x}}) - f^{(f)^n}(\tilde{\underline{X}}^n)}{\Delta t} \\ \text{with} \\ \displaystyle \underline{\tilde{X}}^n = \tilde{\underline{x}}^{n+1} - \Delta t \tilde{\underline{v}}^{(f)^{n+1}}, \end{array} \right. \quad (34) $$

where it is important to notice that $\tilde{\underline{X}}^n$ is the first-order approximation of the position of the particle at the time $t^n$, $\Delta t = t^{n+1} - t^n$. Finding this position, which is called the departure point (or the "foot") of the characteristic line, is the main difficulty of the method. In the current application, an intersection method [FOU 02] is employed to find this position. Notice that other approximation schemes based on higher orders characteristic methods could be used [MAG 94, FOU 02, KAA 03]. The variational formulation finally becomes $\forall \phi^* \in X$:

$$ \int_{\Omega} \left[ 1 + \Delta t \left( \tilde{\mathrm{div}} \tilde{\underline{v}}^{(f)^{n+1}} + \frac{\dot{h}}{h} \right) \right] f^{(f)^{n+1}}(\tilde{\underline{x}}) \phi^* d\Omega = \int_{\Omega} f^{(f)^n}(\tilde{\underline{X}}^n) \phi^* d\Omega. \quad (35) $$

A standard Galerkin finite element scheme yields to a matrix global formulation, which writes as

$$ (\mathbf{M}^{(f)})(F^{(f)^{n+1}}) = (\Psi^{(f)}), \quad (36) $$

where $F^{(f)^{n+1}}$ is the vector of nodal unknowns, $\mathbf{M}$ the mass matrix built on the left-hand side integral of (35), and $\Psi^{(f)}$ a vector based on the right-hand side integral of (35). In this characteristic finite element formulation, quadratic P2 elements are used. In practice the elements of these matrix $\mathbf{M}$ and vector $\Psi^{(f)}$ are approximated using Gaussian quadrature. More precisely, to compute $\Psi^{(f)}$, the Gaussian quadrature points of the elements have to be backtracked: the background elements that contain the departure points of these quadrature points have to be located. Then the scalar values $f^{(f)^n}(\tilde{\underline{X}}^n)$ at these departure points can be determined by interpolation from the nodal values of the background elements at time $t^n$. Once the system is built, it is solved using a conjugate gradient method.
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Reference geometries and boundary conditions used for the simulation examples 1 (a) and 2 (b).</figcaption></figure>

A numerical scheme which is identical to the previous one has been used to compute $\chi^{n+1}$. Therefore, introducing $\varphi^*$ a test function defined in an appropriate space $Y$, the variational formulation of equation (25) is:

$$ \int_{\Omega} \chi^{n+1} \varphi^* d\Omega = \int_{\Omega} \chi^n(\tilde{\mathbf{X}}^n)\varphi^* d\Omega, \quad \forall \varphi^* \in Y. \quad (37) $$

The characteristic Galerkin method yields the following matrix form of this problem:

$$ (\mathbf{M}^{(\chi)}) (\chi^{n+1}) = (\Psi^{(\chi)}), \quad (38) $$

with $\chi^{n+1}$ the vector of nodal unknowns, $\mathbf{M}^{(\chi)}$ a mass matrix built on the left-hand side of integral (37) and $\Psi^{(\chi)}$ a vector built on the right-hand side of (37). Simple linear P1 elements are used in this discretization scheme. Here, the mixture velocity $\tilde{\mathbf{v}}^{(\text{mix})}$ is used to determine the departure point $\tilde{\mathbf{X}}^n$ of the characteristic lines.

## 5. Numerical examples

The solver is implemented in Fortran 95 programming language and uses some routines from the SLATEC mathematical library written in Fortran 77. This program runs either on Windows or Unix systems. The meshes are computed following the standard ".unv" norm available in ProEngineer or I-DEAS softwares. The meshes as well as the input data are transmitted to the solver through a F95 interface. Results are post-treated in Matlab interface. To illustrate the influence of material and process parameters on the segregation phenomenon, we have considered two different numerical tests described below.

In the first example, the different phenomena arising at the beginning of a compression molding are analyzed, so that only the pressure-velocities sub-problem (SP1) is considered. The simulation consists in submitting a rectangular sample of composite of dimensions $L \times l \times h = 100 \times 10 \times 10 \text{ mm}^3$ in dimension to a plane strain compression deformation ($D_{22}^{(\alpha)} = 0$) at a constant axial strain rate $D_{33}^{(\alpha)} = D_{33} = \dot{h}/h$. The sample has a uniform fiber volume fraction $f_0^{(f)}$.
A reference test is defined with the following parameters: fiber volume fraction $f_0^{(f)} = 0.2$, fiber diameter $d = 0.6$ mm, viscosities $\mu_{0(m)} = 0.055$ MPa s and $\eta_{ps}^{(f)} = 0.18$ MPa s, power-law exponents $n_{(m)} = 0.58$ and $n_{(f)} = 0.44$ [DUM 03], permeability coefficient $\kappa^* = 10$, friction coefficients $\kappa_H^{(\alpha)} = 0$, and axial compression strain rate $D_{33} = -0.1$ s$^{-1}$. Figure 1a displays the geometry, the boundary conditions and the mesh used to perform the simulations. The various tests that are carried out in the following are all based on this set of reference parameters and will aim at revealing the influence of each parameter one by one. At last, in order to characterize the fiber segregation phenomenon, a local segregation rate $D_{seg}$ is defined such as

$$D_{seg} = \frac{\mathrm{div}\underline{\mathbf{v}}^{(f)}}{D_{33}} = -\left(\frac{1}{D_{33}f^{(f)}}\right) \frac{D^{(f)}f^{(f)}}{Dt}, \quad (39)$$

so that the higher $D_{seg}$, the higher the segregation phenomenon. On the contrary, an one-phase behavior of the composite is observed when $D_{seg} = 0$.

Figure 2 shows the influence of the axial compression strain rate $D_{33}$ on the fiber segregation rate $D_{seg}$ along the dimensionless X-coordinate $x_1/L$ of the sample for various pairs $n_{(m)}$ and $n_{(f)}$. It can be observed that the "contrast" between the rheology of the two phases plays a key role in the segregation phenomenon. In the case where $n_{(m)} \neq n_{(f)}$, figures 2a and 2b depict that the segregation increases as the axial compression strain rate decreases. This means that a large amount of matrix is expelled from the fiber network. Conversely, for high strain rates, the mixture behaves merely like an one-phase medium. The comparison of figures 2a with 2b underlines the strong influence of $\Delta n = n_{(m)} - n_{(f)}$: the segregation is high when the difference $\Delta n$ is large. If $\Delta n = 0$, there is no influence of the axial strain rate (figures 2c). Figure 3a shows the strong influence of the parameter $\kappa^*$ on the segregation rate $D_{seg}$. This latter parameter depends on the geometry of the microstructure of the fiber network. Low values of $\kappa^*$ are related to a low interaction between the two phases of the mixture and to a two-phase behavior of the composite. On the contrary, for high values of this parameter, the composite behaves merely as an one-phase medium: the permeation of the polymer matrix through the fiber network is hindered. Figure 3b shows that fiber volume fraction $f_0^{(f)}$ also affects $D_{seg}$. With the current set of constitutive equations, more pronounced segregation phenomenon is observed when the fiber volume fraction $f^{(f)}$ is high. Such an influence results from the competition between the reinforcement of the fiber network and the increase of the interaction term as $f^{(f)}$ is increased: such result could not be easily intuited without a FEM simulation.

The second example consists in showing the evolution of segregation during the filling of a channel. Thus, the whole two-phase problem (SP 1-2-3) is solved. The sample is $L_0 \times l \times h = 170 \times 50 \times 10$ mm$^3$ in dimension and is located on the side of the rectangular mold, whose dimensions in $(e_1, e_2)$ are $L_M \times l$ ($L_M = 500$ mm). The mesh and the boundary conditions are given in figure 1b. The reference parameters are the same than those from the previous example except $n_{(f)} = n_{(m)} = 1$. Figure 4 shows the evolution of the fiber volume fraction $f^{(f)}$ along the $x_1$ coordinate of
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Influence of the strain rate $D_{33}$ on the fiber segregation rate $D_{seg}$ for various pairs $n_{(p)}$ and $n_{(f)}$, $L = 100$ mm.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Influence of the permeability parameter $\kappa^*$ (a) and the fiber volume fraction $f^{(f)}$ (b) on the segregation rate $D_{seg}$.</figcaption></figure>

the mold for three different time steps corresponding to three different positions of the front of the mixture, respectively 200, 300 et 450 mm. Two compression molding simulations have been carried out using two values of $\mu_{0(m)}$: the reference one, i.e. $\mu_{0(m)} = 55000$ Pa s that corresponds to the viscosity of SMC matrixes (figures 4a-4c), and $\mu_{0(m)} = 550$ Pa s, corresponding to a very small viscosity of the matrix (figures 4d-4f). The figure brings up the following comments:

- At a given time, the evolution of the fiber volume fraction has a wavy form that increases as the compression is pursued: at the core of the sample the fiber volume fraction $f^{(f)}$ is higher than its initial value $f_0^{(f)}$ (fiber segregation), increases as $x_1$ increases up to a maximum value $f_{max}^{(f)}$, and finally decreases to values lower than $f_0^{(f)}$ near the flow front (matrix segregation).

- Fiber and matrix segregation phenomena are more pronounced when the contrast between the matrix viscosity and the fibre network is higher: for $\mu_{0(m)} = 550$ Pa s, an
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Fiber volume fraction profiles during the channel compression tests at various filling times and for two values of the matrix viscosities: $\mu_{0}(m) = 55000$ Pa s (a-c) and $\mu_{0}(m) = 550$ Pa s (d-f).</figcaption></figure>

important segregation phenomenon is observed, whereas it is much lower for $\mu_{0}(m) = 55000$ Pa s, for which the composite behaves almost as an one-phase medium.

**6. Conclusion**

A two-phase model based on the mixture theory has been proposed to simulate the segregation phenomenon arising during the compression molding of composites such as SMC or GMT. These composites are seen as being formed by two interacting viscous porous continuous media: the polymer matrix phase and the fiber network phase. This model has been simplified to the case of a two-phase shell model, assuming a plug flow for the two-phases of the composite. The model is restricted to isothermal situations and has been implemented in a finite element software. A 2D characteristic/FEM scheme with a moving free surface was developed. Numerical examples given in the present work first illustrate the influence of parameters such as the axial compression strain rate $D_{33}$, the difference $\Delta n = n_{(m)} - n_{(f)}$, the initial fiber volume fraction $f^{(f)}$ and the permeation coefficient $\kappa^*$ on the fiber segregation phenomenon at the beginning of a compression stage. They also illustrate the complex fiber volume fraction $f^{(f)}$ profiles that can be expected during the compression molding of a composite in a simple channel. Future efforts will focus on the comparison between
the predictions of the two-phase model and results of experiments performed either on SMC or GMT compounds.

## Acknowledgements

The authors wish to thank the Région Rhône-Alpes, and the industrial partners Renault, Schneider Electrics, Mecelec-Inoplast and Volvo Trucks for their support.

## 7. References

[ADV 94] ADVANI S. G., Flow and Rheology in Polymer Composites Manufacturing, vol. 10 of Composite Materials Series, Elsevier, 1994.

[AUR 02] AURIAULT J.-L., ROYER P., GEINDREAU C., "Filtration Law for Power-Law Fluids in Anisotropic Porous Media", Int. J. Eng. Sci., vol. 40, 2002, p. 1151-1163.

[BAR 85] BARONE M. R., CAULK D. A., "Kinematics of Flow of Sheet Molding Compounds", *Polymer Composites*, vol. 6, num. 2, 1985, p. 105-109.

[BAR 86] BARONE M. R., CAULK D. A., "A Model for the Flow of a Chopped Fiber Reinforced Polymer Compound in Compression Molding", *J. Appl. Mech.*, vol. 53, num. 191, 1986, p. 361-371.

[BAR 87] BARONE M., OSSWALD T., "Boundary integral equations for analyzing the flow of a chopped fiber reinforced polymer compound in compression molding", *J. Non-Newtonian Fluid Mech.*, vol. 26, 1987, p. 185-206.

[BOW 76] BOWEN R. M., *Continuum Physics III: Mixtures and EM Field Theories*, vol. III, Academic Press, 1976.

[BRU 93] BRUSCHKE M. V., ADVANI S. G., "Flow of Generalized Newtonian Fluid Across a Periodic Array of Cylinders", *J. Rheol.*, vol. 37, num. 3, 1993, p. 479-498.

[DHA 81] DHATT G., TOUZOT G., *Une Présentation de la Méthode Des Éléments Finis*, Maloine S.A. Editeur, Paris, et Les Presses de l'Université Laval, Québec, 1981.

[DHA 92] DHATT G., GAO D. N., MÉDAL M., SONG R., "Remplissage et Solidification Par Éléments Finis", *Rev. Eur. Elém. Finis*, vol. 1, 1992, p. 309-332.

[DUM 03] DUMONT P., ORGÉAS L., CORRE S. L., FAVIER D., "Anisotropic Viscous Behavior of Sheet Molding Compounds (SMC) During Compression Molding", *Int. J. Plasticity*, vol. 19, 2003, p. 625-646.

[FOU 02] FOURESTEY G., "Une Méthode Des Caractéristiques D'ordre Deux sur Maillages Mobiles Pour la Résolution Des Équations de Navier-Stokes Incompressible Par Éléments Finis", report num. 4448, 2002, INRIA.

[HIE 80] HIEBER C., SHEN S., "A finite-element/finite-difference simulation of the injection-molding filling process", *J. of Non-Newtonian Fluid Mech.*, vol. 7, 1980, p. 1-32.

[IDR 04] IDRIS Z., ORGÉAS L., GEINDREAU C., BLOCH J. F., AURIAULT J. L., "Microstructural Effects on the Flow Law of Power-Law Fluid Through Fibrous Media", *Modelling Simul. Mater. Sci. Eng.*, vol. 5, 2004, p. 995-1116.

[JAC 86] JACKSON G. W., JAMES D. F., "The Permeability of Fibrous Porous Media", *Can. J. Chem. Eng.*, vol. 64, 1986, p. 364-374.
[KAA 03] KAAZEMPUR-MOFRAD M. R., MINEV P. D., ETHIER C. R., "A Characteristic/Finite Element Algorithm for Time-Dependent 3-D Advection-Dominated Transport Using Unstructured Grids", *Comput. Meth. Appl. Mech. Eng.*, vol. 192, 2003, p. 1281-1298.

[KEL 64] KELLER J. B., "Viscous Flow Through a Grating or Lattice of Cylinders", *J. Fluid Mech.*, vol. 18, num. 1, 1964, p. 94-96.

[LEC 02] LE CORRE S., ORGÉAS L., FAVIER D., TOURABI A., MAAZOUZ A., VENET C., "Shear and Compression Behaviour of Sheet Moulding Compounds", *Compos. Sci. Technol.*, vol. 62, 2002, p. 571-577.

[LEE 84] LEE C. C., FOLGAR F., TUCKER C. L., "Simulation of Compression Molding for Fiber-Reinforced Thermosetting Polymers", *Polymer Composites*, vol. 106, 1984, p. 114-125.

[LIU 01] LIU W., XU S., "A New Improved Uzawa Method for Finite Element Solution of Stokes Problem", *Comput. Mech.*, vol. 27, 2001, p. 305-310.

[MAG 94] MAGNIN B., "Modélisation Du Remplissage Des Moules D'injection Pour Polymères Thermoplastiques Par Une Méthode Eulérienne-Lagrangienne Arbitraire", PhD thesis, École Nationale Supérieure des Mines de Paris, 1994.

[ODE 04] ODENBERGER P., ANDERSON H., LUNDSTRÖM T., "Experimental Flow Front Visualisation in Compression Moulding of SMC", *Composites: Part A*, vol. 35, 2004, p. 1125-1134.

[OSS 90] OSSWALD T. A., TUCKER C. L., "Compression mold filling simulation for non-planar parts", *Intern. Polymer Processing*, vol. 2, 1990, p. 79-87.

[OSS 94] OSSWALD T. A., TSENG S. C., "Flow and Rheology in Polymer Composites Manufacturing", vol. 10 of *Composite Materials Series*, chapter C8 - Compression Molding, p. 361-413, Ed. S. G. Advani - Elsevier, 1994.

[ROB 90] ROBICHAUD M. P., TANGUY P. A., FORTIN M., "An Iterative Implementation of the Uzawa Algorithm for 3-Fluid Flow Problems", *Int. J. Num. Methods in Fluids*, vol. 10, 1990, p. 429-442.

[SCA 99] SCARDOVELLI R., ZALESKI S., "Direct Numerical Simulation of Free-Surface and Interfacial Flow", *Annu. Rev. Fluid Mech.*, vol. 31, 1999, p. 567-603.

[SER 02] SERVAIS C., LUCIANI A., MÁNSON J.-A. E., "Squeeze flow of concentrated long fibre suspensions: experiments and model", *J. Non-Newtonian Fluid Mech.*, vol. 104, 2002, p. 165-184.

[SOU 01] SOULI M., ZOLESIO J. P., "Arbitrary Lagrangian-Eulerian and Free Surface Methods in Fluid Mechanics", *Comput. Methods Appl. Mech. Eng.*, vol. 191, 2001, p. 451-466.

[THO 96] THOMASSON J. L., VLUG M. A., "Influence of Fibre Length and Concentration on the Properties of Glass Fibre-Reinforced Polypropylene: 1. Tensile and Flexural Modulus", *Composites: Part A*, vol. 27A, 1996, p. 477-484.

[TRU 60] TRUESDELL C., TOUPIN R., vol. III1, S Flugge, Berlin, 1960.

[YAG 95] YAGUCHI H., HOJO H., LEE D. G., KIM E. G., "Measurement of planar orientation of fibers for reinforced thermoplastics using image processing", *Intern. Polymer Processing*, vol. 3, 1995, p. 262-269.

[ZIE 00] ZIENKIEWICZ O. C., TAYLOR R. L., *The Finite Element Method, Fifth Edition. Volume 1: The Basis*, Butterworth-Heinemann, Oxford, UK, 2000.