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

A two-phase approach is proposed to model the rheology of polymer glass-fiber compounds 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.


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 theses 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 (Osswald et al., 1994;Yaguchi et al., 1995;Thomasson et al., 1996). 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 (Advani, 1994). 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.

Balance equations for a two-phase mixture
The theory of mixture that was initially developed in the pioneering works of Truesdell and Toupin (Truesdell et al., 1960;Bowen, 1976) 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 of them represents immiscible phase of the material, i.e. the fiber network f and the matrix m (polymer + filler). Thus, 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.
-each elementary macroscopic volume element δV of the mixture (elementary mass δm and specific mass ρ = δm δV ) is simultaneously occupied by the phases α whose elementary macroscopic mass δm (α) occupies an elementary macroscopic volume δV (α) of δV 1 . This enables to introduce (i) the macroscopic ρ (α) and the microscopic ρ (α) specific masses and (ii) the volume fraction f (α) of the phase α: [1] In this work, the constituents of the mixture are assumed incompressible at the microscopic scale so that ρ (α) is constant and ρ (α) = f (α) ρ (α) .
-the saturation condition for the mixture is stated: -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).
Accounting for the previous assumptions the local mass balance equation written at a given material point M (α) for each phase α is where the notation D (α) Dt is the material time derivative following materials points M (α) of velocity v (α) . Summing equations [3] written for each phase of the mixture, i.e. m and f , and accounting for the saturation of the mixture [2] leads to an incompressibility equation for the mixture: [4] The mixture theory defines the concept of partial stress vector T (α) and partial stress tensor σ (α) for each phase α, writing the total force exerted from the outside on the surface element δS of δV as : where n is an outward normal vector to δ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. σ (α) = σ T (α) . Moreover, rheometry 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 notation for microscopic and macroscopic quantities is introduced: if φ is a scalar physical quantity, then φ (α) (with the subscript (α) ) is the microscopic (or local) value of φ in the phase α, whereas φ (α) (with the exponent (α) ) is the macroscopic (or apparent) value of φ for the same phase in the mixture. 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 (Servais et al., 2002;Le Corre et al., 2002;Dumont et al., 2003). In this case, it can be shown that the partial stress tensors σ (α) can be split into the sum of two terms (Bowen, 1976): where p i is a pressure related to the incompressibility of the mixture, σ e(m) and σ 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 (Bowen, 1976): where p e(m) is a viscous momentum exchange between the two phases.

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 (Hieber et al., 1980;Lee et al., 1984). 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 onephase Barone-Caulk shell models (Barone et al., 1986;Barone et al., 1987;Osswald et al., 1990). 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 (Barone et al., 1985;Servais et al., 2002;Odenberger et al., 2004). Consequently, in order to simplify the full 3D two-phase model introduced in the previous section, a twophase 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 (Barone et al., 1986). 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 (e 1 , e 2 ) plane: the extension to a 3D shell model would be straightforward. Hence, the velocity v (α) of each phase becomes: where the symbol "˜" has been introduced to distinguish 2D-unknown fields and 2Doperators in (e 1 , e 2 ) from 3D ones, and h is the thickness of the sheets. As the shear components D (α) β3 (β ∈ {1, 2}) of the strain rate tensor are constrained to zero, arbitrary reaction terms must be added to the partial stress tensor σ (α) : Hence, the mass balance equation for the fiber network becomes: Moreover, the incompressibility of the mixture now reads: 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: h 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 (Barone et al., 1986) (cf. next subsection).

Constitutive equations
-It seems reasonable to consider that the polymer matrix behaves like a powerlaw fluid at the microscopic scale, with a power-law exponent n (m) (Le Corre et al., 2002;Dumont et al., 2003). Following this idea, the extra macroscopic stress tensor of the matrix is described by a simple constitutive law: , µ 0(m) being the viscosity of the polymeric matrix at a reference shear rateγ 0 of 1 s −1 . Accounting for the current kinematical restrictions, -The extra stress tensorσ e(f ) of the fiber network is a direct extension of the onephase mechanical model proposed by Dumont et al. (2003). Hence, the fiber network is seen as a compressible power-law fluid displaying a transverse isotropy whose axis is e 3 :σ where with n (f ) the power-law exponent, η (f ) ps a plane strain compression viscosity at a characteristic strain rate D 0 of 1 s −1 , and H a rheological function that accounts for the fiber network anisotropy. η (f ) ps and H strongly depend on f (f ) . For example, using the experimental results obtained on a standard SMC formulation (Dumont et al., 2003), possible expressions for η (f ) ps and H are where η (f ) is a constant.
-The momentum exchange˜ p is induced when there is a relative motionṽ r = v (m) −ṽ (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 [7] 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 (Auriault et al., 2002;Idris et al., 2004): [20] 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 κ * is a constitutive parameter to determine. According to the experimental data collected for various fibrous media: 0.1 < κ * < 100 (Jackson et al., 1986). It can be shown that et al., 1993;Idris et al., 2004): where f was approximated by the analytical lubrication model (Keller, 1964), that gives rather good predictions in a wide range of fiber volume fraction: [22] -In their one-phase shell-model, Barone et Caulk (Barone et al., 1986) have given a physical meaning for the forcesF appearing in [13][14]: 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: where v 0 is a reference velocity and κ (α) H are friction parameters to determine.

Initial and boundary conditions
At the beginning of the compression, i.e. for t = 0, the composite occupies a surface Ω(x 1 , x 2 , t = 0) = Ω 0 of boundary ∂Ω(x 1 , x 2 , 0) = ∂Ω 0 in the principal plane (e 1 , e 2 ) of the mold whose total surface and boundary are Ω M (Ω 0 ⊂ Ω M ) and ∂Ω M , respectively. The initial height h(x 1 , x 2 , 0) = h 0 and fiber volume fraction 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 Ω(x 1 , x 2 , t) of boundary ∂Ω(x 1 , x 2 , t) fills the mold cavity, accounting for the following set of boundary conditions: where ∂ σ Ω stands for the free surface of the mixture.

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 Ω 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 Ω grows from Ω 0 to Ω M and the evolution of the free surface ∂ σ Ω has to be determined. As in many mold filling applications (Dhatt et al., 1992;Scardovelli et al., 1999;Souli et al., 2001), this was achieved with an additional scalar variable χ describing the local volume fraction of the composite in the mold (i.e. χ = 0/1 if the considered material point is empty/full of composite). The dynamic of the variable χ(x 1 , x 2 , t) is ruled by the following equation: where D (mix) Dt stands for the material time derivative following the mixture of velocitỹ v (mix) . As defined in (Bowen, 1976): Thus, the scalar field χ is added to the unknown fields f (f ) ,ṽ (m) ,ṽ (f ) andp 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 Ω n+1 occupied by the mixture, and (ii) determining the unknown fields in Ω n+1 . These two points are detailed below.
(i) The calculation of the new domain Ω n+1 from the knowledge of Ω n , v (mix) n and χ n requires (a) the computation of the time increment ∆t = ∆t n→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 + ∆t, (c) the calculation of h n+1 andḣ 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 andp n+1 i , and couples the incompressibility equation [12] and the two momentum balance equations [13][14]. This first sub-problem would be equivalent 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 [11] 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ṽ (m) n+1 andṽ (f ) n+1 , pressurep n+1 i and fiber volume fraction fields f (f ) n+1 .
In what follows, the different numerical methods used to solve these problems are presented.

Weak formulation
The pressure-velocities problem [SP1] is re-casted into a weak form by multiplying the equations [13-14] and [12] by a set of test functionsṽ (m) * ,ṽ (f ) * and p * (belonging to appropriate spaces V It is important to underline that this system is highly non-linear due to the power-laws

FE approximations
A mixed finite element method is used to discretize the previous weak formulation. The considered elements are P 2+ for the two velocity fields and P 1 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: where V is the vector of nodal velocities unknowns (m and f ), i.e.
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ḣ/h. More precisely, one can write : where the matrices A 1 (α) , A 2 (α) and A 3 (f /m) are "viscous" matrices arising respectively from the rheology of the phases, the friction terms, and the momentum exchange, and where B (α) stands for the incompressibility condition. In the case of a linear system (n (α) = 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 [27] defined as where r is a constant (Robichaud et al., 1990;Liu et al., 2001). Notice that the problem [30] can be generally written into a linear system of equations: 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 pressurevelocities problem combines the Newton-Raphson and Uzawa algorithms. Given r, nr max , ε N R , ρ, k max , ε U Z , an arbitrary choice of P nr=1 for P and V nr=1 for V: Compute ∆ V k nr by solving the following equation: where nr max , ε N R are respectively the maximum number of iterations and the convergence criterion, where k max and ε U Z are respectively the maximum number of iterations and the convergence criterion of the Uzawa algorithm, and where 0 < ρ < 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 (Dhatt et al., 1981;Zienkiewicz et al., 2000). The solver of the matrix system given by combining the Newton-Raphson and Uzawa algorithms is based on a biconjugate gradient algorithm.

Mass balance equation for the fiber phase and free surface transport
-The solution of the fiber volume fraction evolution equation [11] is obtained through a finite element discretization in space and a characteristic based method discretization in time. The variational formulation of the equation [11] reads ∀ φ * ∈ X: where f (f )n+1 is the solution of the fiber volume fraction field at time t n+1 and φ * a test function defined in an appropriate space X. From [32], the characteristic method is applied to approximate the material derivative D (f ) Dt f (f ) n+1 . It consists in writing the following finite difference scheme: where it is important to notice thatX n is the first-order approximation of the position of the particle at the time t n , ∆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 (Fourestey, 2002) is employed to find this position. Notice that other approximation schemes based on higher orders characteristic methods could be used (Magnin, 1994;Fourestey, 2002;Kaazempur-Mofrad et al., 2003). The variational formulation finally becomes ∀ φ * ∈ X: [34] A standard Galerkin finite element scheme yields to a matrix global formulation, which writes as where F (f ) n+1 is the vector of nodal unknowns, M the mass matrix built on the lefthand side integral of [34], and Ψ (f ) a vector based on the right-hand side integral of [34]. In this characteristic finite element formulation, quadratic P2 elements are used. In practice the elements of these matrix M and vector Ψ (f ) are approximated using Gaussian quadrature. More precisely, to compute Ψ (f ) , the Gaussian quadrature points of the elements have to be backtracked: the background elements that contain the departure points of theses quadrature points have to be located. Then the scalar  values f (f ) n (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.
-A numerical scheme which is identical to the previous one has been used to compute χ n+1 . Therefore, introducing ϕ * a test function defined in a appropriate space Y , the variational formulation of equation [24] is : [36] The characteristic Galerkin method yields the following matrix form of this problem: with χ n+1 the vector of nodal unknowns, M (χ) a mass matrix built on the left-hand side of integral [36] and Ψ (χ) a vector built on the right-hand side of [36]. Simple linear P1 elements are used in this discretization scheme. Here, the mixture velocitỹ v (mix) is used to determine the departure pointX n of the characteristic lines.

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  (Dumont et al., 2003), permeability coefficient κ * = 10, friction coefficients κ (α) H = 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 or-der to characterize the fiber segregation phenomenon, a local segregation rate D seg is defined such as 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) = 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 ∆n = n (m) − n (f ) : the segregation is high when the difference ∆n is large. If ∆n = 0, there is no influence of the axial strain rate (figures 2c). Figure 3a shows the strong influence of the parameter κ * on the segregation rate D seg . This latter parameter depends on the geometry of the microstructure of the fiber network. Low values of κ * 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 (f ) 0 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 × l × h = 170 × 50 × 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 ×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 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 µ 0(m) : the reference one, i.e. µ 0(m) = 55000 Pa s that corresponds to the viscosity of SMC matrixes (figures 4a-4c), and µ 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 (f ) 0 (fiber segregation), increases as x 1 increases up to a maximum value f (f ) max , and finally decreases to values lower than f (f ) 0 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 µ 0(m) = 550 Pa s, an important segregation phenomenon is observed, whereas it is much lower for µ 0(m) = 55000 Pa s, for which the composite behaves almost as an one-phase medium.

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 ∆n = n (m) −n (f ) , the initial fiber volume fraction f (f ) and the permeation coefficient κ * 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.