# Compression moulding of SMC: in situ experiments, modelling and simulation

Pierre J.J. Dumont, Laurent Orgéas, Denis Favier, Patrick Pizette, Cécile Venet

► To cite this version:

Pierre J.J. Dumont, Laurent Orgéas, Denis Favier, Patrick Pizette, Cécile Venet. Compression moulding of SMC: in situ experiments, modelling and simulation. Composites Part A: Applied Science and Manufacturing, 2007, 38A (2), pp.353-368. 10.1016/j.compositesa.2006.03.010. hal-00292081

HAL Id: hal-00292081

https://hal.science/hal-00292081v1

Submitted on 26 Feb 2021

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 CC BY 4.0 - Attribution - International License
# Compression moulding of SMC: In situ experiments, modelling and simulation

P. Dumont<sup>a</sup>, L. Orgéas<sup>b,*</sup>, D. Favier<sup>b</sup>, P. Pizette<sup>c</sup>, C. Venet<sup>c</sup>

## Abstract

Compression mouldings of commercial SMC were performed with an instrumented industrial press under various process conditions. Results underline the influence of process parameters such as the initial SMC temperature, the axial punch velocity and the geometry of the mould on local normal stress levels. They also show negligible fibre-bundle segregation in the principal plane of the moulded parts. Thereby, a one-phase plug flow shell model is proposed as a direct extension of the plug flow model proposed by M.R. Barone and D.A. Caulk [J Appl Mech 53(191):1986;361-70]. In the present approach, the SMC is considered as a power-law viscous medium exhibiting transverse isotropy. The shell model is implemented into a finite element code especially developed for the simulation of compression moulding of composite materials. Simulation and experimental results are compared, emphasizing the role of the SMC rheology on the overall recorded stress levels. Despite the simplicity of the model, rather good comparisons are obtained.

**Keywords:** A. Polymer-matrix composites (PMCs); E. Compression moulding; C. Analytical modelling; C. Computational modelling

## 1. Introduction

Compression moulded composites such as sheet moulding compounds (SMC) or alternatively glass mat thermoplastics (GMT) are widely and increasingly used in different industry branches to produce thin semi-structural parts. SMC are made of continuous sheets (thickness ≈2–3 mm) containing ≈15–35 wt% chopped glass bundle mats (length ≈10–50 mm, averaged diameter ≈0.05–1 mm) embedded in a pasty matrix essentially made of CaCO$_3$ (≈60 wt%) and polyester resin (≈20 wt%). A SMC compression moulding process encompasses the different following stages. In a first stage, a charge (≈25 °C) of stacked plies cut from the SMC sheet is placed on the lower platen of a hot mould (≈150 °C) for about 5–20 s before

the mould is closed. Then, in a second stage as short as around 5 s, the charge is squeezed at a closure velocity of ≈0.1–10 mm s⁻¹. This stage involves a considerable deformation of the SMC inside the mould cavity, the usual height reduction being around 30–70%. Once the mould is filled, it is maintained closed during a third curing stage for a period of about 60–120 s. Even if the flow stage is very short, the different flow-induced phenomena, such as bundles orientation evolution and segregation, knit lines, etc., strongly affect the structural and physical properties, the surface aspect and the geometrical stability of produced parts. For this reason, there is a great interest of the industrial SMC manufacturers in simulating this flow stage in order to reduce the long and empirical phases of design and development of the moulds.

The flow of SMC during compression moulding is highly complex because it strongly depends on the rheology of the SMC, its evolving microstructure, the curing of the polyester resin, thermal and mechanical boundary

<footnote><sup>a</sup> Laboratoire de Technologie des Composites et Polymères (LTC), École Polytechnique Fédérale de Lausanne (EPFL), Station 12, CH-1015 Lausanne, Switzerland</footnote>

<footnote><sup>b</sup> Laboratoire Sols-Solides-Structures (3S), CNRS – Universités de Grenoble (INPG-UJF), BP 53, 38 041 Grenoble cedex, France</footnote>

<footnote><sup>c</sup> Schneider Electric Industries SAS, Science and Technology Division – Materials Research, 37 Quai P.L. Merlin, 38 000 Grenoble, France</footnote>

<footnote>* Corresponding author. Fax: +33 4 76 82 70 43.
E-mail address: Laurent.Orgeas@hmg.inpg.fr (L. Orgéas).</footnote>
conditions between the SMC and the mould. Building flow models allowing a complete treatment of the couplings between all these different phenomena is heavy and is actually replaced by some simplification assumptions. For instance, since the flow stage is very short, it is reasonable to neglect the resin cure during this stage [3]. Likewise, paste segregation is known to occur in complex parts of the mould such as ribs [44,11], but is rarely identified in other zones: as a consequence, most of the flowing models outside the ribs are one-phase models.

Other approximations of the SMC flow (outside the ribs) have been stated, taking advantage of the small thickness of the mould cavity with respect to the other in-plane dimensions. They have given rise to simplified descriptions of the flow kinematics of compression moulded materials. Among them, a first group is based on the lubrication approximation, mainly characterized by (i) the shearing of the composite through the thickness of the mould cavity and (ii) the perfect sticking of the composite with the upper and lower parts of the mould. Such a deformation mechanism has been used to establish one-phase generalized Hele-Shaw shell models [41,16,25,2]. A second group of descriptions adopt the “plug flow” approximation, characterized by (i) a uniform in-plane shear and elongational deformation (ii) the absence of out of plane shear in the core of the stacked sheets and (iii) possible sheared boundary layers near the surface of the cavity. This second type of deformation mechanism has given rise to one-phase plug flow shell models [6,5,34]. Also notice that both approaches consider the SMC charge as a one-layer sample and does not account for its stratified structure. This point was studied for GMT materials by Leterrier and G'Sell [28] who developed a multilayer plug flow model.

In practice, the plug flow as well as the through-thickness shear mechanism exist, the predominance of the first or second mechanism depends on processing conditions, thickness of the sheets, length of fibres, etc. In the case of SMC or GMT, which bundles are much longer than the thickness of the mould cavity, experimental evidences show that within a wide range of processing conditions, the plug flow regime leads to an acceptable approximation of the real flow patterns. Indeed, the experiments performed by Barone and Caulk [4] have revealed that “plug flow” is the dominating flow pattern of SMC materials moulded in industrial conditions (mould temperature at 150 °C, thin charges of thickness $h \le 15$ mm, high compression velocities $\approx 1.75–10$ mm s⁻¹). Performing isothermal squeeze flow experiments on SMC at compression velocities ranging from 0.1 to 3 mm s⁻¹, Kotsikos and Gibson [21] have established the same conclusions by analyzing the evolution of the measured normal stress along the radius of samples. Similar trends have been obtained by Servais et al. [39] for the compression moulding of GMT. More recently, Odenberger et al. [33] have analysed more closely the flow of SMC during compression using a camera. The authors have identified three deformation mechanisms. A first and rapid “squish” phase is observed, during which the

flow is very complex, especially at the flow front where the upper and lower plies are often ejected from the stack. This short phase is followed by a steady “plug flow”. At the last step of compression, a “boiling” phase was observed, ascribed to the boiling of styrene. Finally, it is worth noting the plug flow assumption is also supported by the theoretical and numerical work of Lee and Tucker [26]. These authors have shown that the plug flow is a consequence of the interaction of heat transfers localised near the contact with the mould cavity, flow and temperature-dependent viscosity and that it is likely to happen for SMC in realistic processing conditions.

The equations governing the momentum balance in the plug flow model involve a pressure term, related to the incompressibility of the SMC, a term related to the bulk rheology of the SMC and a friction term, introduced to model the interaction between the SMC and the mould cavity:

* Three types of friction terms have been tested by Barone and Caulk [6]: pure and Coulombic dry frictions, and linear hydrodynamic frictions. A better agreement with their experiments was found with the last solution, which physical motivation is the existence of a thin lubricating paste-rich boundary layer adjacent to the upper and lower parts of the mould. More recently, accounting for the non-Newtonian rheology of the paste, the hydrodynamic friction term was modified using a power-law approximation [1].

* The term that accounts for the bulk rheology of the SMC was initially proposed to be transversely isotropic, purely viscous and linear [6]. It requires one additional constitutive parameter, i.e. a constant in-plane extensional viscosity. Different authors have tried to determine this parameter [9,10,1] and showed that such a contribution was not always negligible with respect to the friction term.

The last point calls up the two following comments. Firstly, the rheological model proposed by Barone and Caulk [6] is a first attempt to model the very complex rheology of SMC, strongly affected by the bundles content, bundles orientation, imposed temperatures and strain rates [40,27,32,18,19,29,30,24,13]. It could be interesting to introduce in the plug flow model rheological models more adapted to the rheology of the SMC. Secondly, the rheological term, even in its linear form, is rarely taken into account in commercial simulation softwares, to the best of our knowledge. Instead, a simplified form of the plug flow equations is used, i.e. the very thin charges approximation model, in which only the pressure and friction terms are modelled. The resulting momentum balance equations are formally identical to those developed following lubrication approximations [6]. If this very strong assumption permits to describe very quickly the filling of the mould with a minimum of unknown field (i.e. the pressure only), it may not give correct mould closure force [35].
Hence, the main objectives of the present contribution are (i) to propose a plug flow model following the pioneering work of Barone and Caulk [6] and using a rheological model able to capture the effect of strain rate, temperature, bundle content and orientation on the mechanical response of the SMC, (ii) to test its capability to reproduce compression experiments performed in industrial conditions, and (iii) to implement it in a specially designed simulation code. For that purpose, compressions in a rectangular channel were performed using an instrumented industrial press and with a commercial SMC formulation (Section 2). Results first exhibit little segregation phenomena in the plane of the compressed part and underline the roles of the initial temperature of the SMC, the initial position of the plies in the mould, the mould geometry and the axial punch velocity (Section 3). The plug flow shell model is then formulated: the rheological model used to describe the bulk behavior of SMC materials is a nonlinear and anisotropic viscous model developed in a previous study [13] and the friction model is the power-law friction law proposed by Abrams and Castro [1]. Analytical expressions for compressions in a channel are established and compared with the experiments (Section 4). At last, a finite element (FE) formulation of the developed shell model is proposed and implemented in a FE application, allowing a comparison between experimental and numerical results (Section 5).

## 2. Experimental procedure

### 2.1. Material and equipments

The tested material is a standard Low Profile SMC (LP 606) supplied by Mecelec Composite and Recyclage (Tournon, France). Produced sheets are made of a polyester-based paste matrix (see Table 1 for the composition) and 25 wt% of glass fibre bundles (length $l = 25$ mm, composed of ≈200 fibres with a diameter $d = 15$ µm). Note that the formulation characterized in previous works [24,13] was identical to the present one, except it contained an inhibitor to prevent curing during rheological tests.

Compression moulding experiments were performed with an industrial hydraulic press (maximum axial force of 4000 kN or a maximum axial displacement velocity of 40 mm s⁻¹, Compositec, Chambéry, France). The press was equipped with a standard rectangular mould that could be used to produce thin (≈1 mm) or thick (≈15 mm) plates with a in-plane ($e_1, e_2$) surface of 500 × 500 mm². Two pres-

<figure><img src="image_2.png" /><figcaption>Fig. 1. Top views of the lower part of the mould with the initial placements of SMC charges for the (SL, specimen-left) (a) and the (SC, specimen-center) (b) mouldings, side view of the modified mould (c). Notice that the projections of the pressure sensors located in the upper part of the mould have been added in the views (a) and (b).</figcaption></figure>

sure sensors were located in the upper part of the mould (punch). Their locations are given in Fig. 1: the first one (I) lies at the center of the mould whereas the second one (II) is near the lower left corner of the punch.

### 2.2. Compression experiments

All the tests were carried out using a constant mould temperature of 150 °C and two axial constant punch velocities $\dot{h}$ of 1 mm s⁻¹ and 10 mm s⁻¹. The compression began as soon as (≈10 s) the sample was put into the mould. The initial in-plane ($e_1, e_2$) surface of the sample was $S_0 = l_0L_0 = 170 \times 500$ mm². Incomplete compressions showed that the flow front was fairly rectilinear, so that the flow could be considered as parallel to $e_1$ and independent of $x_2$. At last, different initial temperatures of the SMC samples were used, respectively, 30, 40 and 60 °C. For the last two cases, a pre-heating time over approximately a 10 min period was needed to bring the samples to the prescribed temperature.

A first set of experiments was carried out with the standard mould using samples made of three stacked layers (average initial thickness $h_0$ of 8.15 mm). Two initial sample positions were adopted to investigate more precisely the SMC flow. In the first case (SL, sample-left), the samples were initially located at left side of the mould, as depicted in Fig. 1(a), whereas in the second case (SC, sample-center) samples were placed at the center of the mould (see Fig. 1(b)). Note that when $\dot{h} = 10$ mm s⁻¹, it was not possible to achieve the compression mouldings with a good control of $\dot{h}$, so that no reliable results were obtained for this velocity and this initial sample placement.

<table><caption>Table 1
Composition of the paste (wt%)</caption><tr><td>Polyester resin</td><td>19.60%</td></tr><tr><td>Styrene</td><td>1.56%</td></tr><tr><td>Catalytor + inhibitor</td><td>0.56%</td></tr><tr><td>ZnSt</td><td>1.31%</td></tr><tr><td>M additive</td><td>0.76%</td></tr><tr><td>Low profile additive</td><td>13.29%</td></tr><tr><td>MgO</td><td>0.85%</td></tr><tr><td>CaCO<sub>3</sub></td><td>62.07%</td></tr></table>
A second set of compression tests was achieved in order to study the flow of the SMC in a more complex situation. For that purpose, the lower part of the standard mould was modified by adding an intermediate polished aluminum part that contained two successive smooth convergents (height = 2 mm, length = 80 mm), as depicted in Fig. 1(c). For these experiments, only one punch velocity was investigated (h = 1 mm s⁻¹), and the samples were initially located at the left side of the mould (Fig. 1(b)) and formed of six layers (average initial thickness of 16 mm).

During the experiments the axial punch displacement h as well as the signals of the two pressure sensors were systematically recorded. The last two signals directly give the values of local axial stresses σ₃₃ and σ₃₃.

## 2.3. Fibre-bundle segregation measurements

A burning and weighting technique was used to determine the fibre content after the compression moulding. It is very similar to that proposed by Christensen et al. [11] and consists of (i) burning the organic constituents of the matrix and (ii) degrading some of the mineral ones.

In order to detect the calcination and the degradation temperatures of the organic and mineral constituents, a TGA experiment was first performed on a small amount of matrix (volume ≈ 1 mm³). The resulting evolutions of the mass m of the sample and its first derivative dm/dT with respect to the temperature T are plotted in Fig. 2. The dm/dT curve clearly exhibits three peaks. Two of them occur below 500 °C and are attributed to the calcination of the organic constituents (polyester resin + other organic additives). The last one, located between 500 and 800 °C, is mainly induced by the partial degradation of the CaCO₃ fillers (CaCO₃ → CaO + CO₂).

It was then possible to determine the local volume fraction of fibres fᵢ adopting the following procedure:

(1) determine the initial mass of the sample mᵢ<sup>initial</sup>,

(2) heat the sample at 500 °C to burn the organic constituents of mass mᵢ<sup>organic</sup> and measure the mass mᵢ<sup>initial</sup> of the sample to obtain:

<figure><img src="image_3.png" /><figcaption>Fig. 2. TGA measurements. Evolution of the mass m of the sample (%) and of its first derivative dm/dT (% °C⁻¹) as functions of the temperature.</figcaption></figure>

$$m_{i}^{\text{organic}} = m_{i}^{\text{initial}} - m_{i}^{500}, \qquad (1)$$

(3) heat again the sample at 800 °C to degrade a part of the mineral constituents of initial mass mᵢ<sup>mineral</sup> and measure its new mass mᵢ<sup>800</sup>. Introducing k_d, the degradation coefficient of mineral fillers (0 ≤ k_d ≤ 1), one can write:

$$m_{i}^{\text{mineral}} = \frac{m_{i}^{500} - m_{i}^{800}}{k_{d}}, \qquad (2)$$

so that the mass fraction of fibre-bundles τᵢ reads

$$\tau_{i} = \frac{m_{i}^{\text{fibre}}}{m_{i}^{\text{initial}}} = \frac{m_{i}^{800} - (1 - k_{d})m_{i}^{500}}{k_{d}m_{i}^{\text{initial}}} \qquad (3)$$

and the fibre-bundle volume fraction fᵢ becomes

$$f_{i} = \frac{\tau_{i}\rho_{\text{paste}}}{\tau_{i}\rho_{\text{paste}} - \rho_{\text{fibre}}(\tau_{i} - 1)}, \qquad (4)$$

where ρᵢ is the volumetric mass of constituent i. Notice that Christensen et al. [11] used a theoretical value of 0.44 for the degradation coefficient k_d. Here, k_d was determined experimentally, adding a fourth step to the procedure described above on several undeformed samples. It consisted of performing a chemical attack (HCl) of the remaining mineral fillers to make them totally disappear from the burned sample: the mass of the sample then reduces to mᵢ<sup>fibre</sup>. From the knowledge of mᵢ<sup>fibre</sup>, mᵢ<sup>800</sup> and mᵢ<sup>500</sup>, it was then possible to estimate k_d from the measurements:

$$k_{d} = \frac{m_{i}^{500} - m_{i}^{800}}{m_{i}^{500} - m_{i}^{\text{fibre}}} = 0.41 \pm 0.025. \qquad (5)$$

This coefficient, which is slightly different from the theoretical value of 0.44 [11], was used in (3) for the rest of the analyzed samples.

<figure><img src="image_4.png" /><figcaption>Fig. 3. Influence of the sample size a of the cut samples on fᵢ/f₀ (open circles). The black circles represent the mean value of fᵢ/f₀ calculated at fixed value of a.</figcaption></figure>
## 2.4. Initial variability of the fibre-bundle content

The influence of the size of the cut samples on $f_i$ was investigated on undeformed sheets. For this reason, many samples of in-plane dimensions of $a \times a$, ($a = 25, 50$ and $100$ mm) were cut in sheets of mean fibre-bundle volume fraction $f_0$ ($\approx 17.4\%$). The estimated values $f_i/f_0$ as functions of the lengths of the samples are given in Fig. 3. This figure clearly reveals that the calculated mean value (black circles in the figure) weakly depends on the sample dimensions and is equal to the mean fibre-bundle volume fraction $f_0$, provided the number of samples is high enough. Moreover, a higher scattering of $f_i$ is observed for smaller lengths: the estimated standard deviation is $\approx 10\%$, $6.5\%$ and $6\%$ when $a$ equals $25, 50$ and $100$ mm, respectively. This is essentially attributed to the initial local heterogeneity of the fibre-bundle content. Note that similar trends have been observed

previously concerning the heterogeneity of the local strain field of SMC samples deformed homogeneously at the macroscopic scale [13]. Accounting for the last results and for the studied “1D” type of flow (in $(e_1, e_2)$), the fibre-bundle content will be next analysed using mean values, each of them being calculated with 10 samples of in-plane dimensions of $25 \times 25$ mm$^2$ of identical $x_1$-coordinate.

# 3. Experimental results

## 3.1. Standard mould

Fig. 4 shows a collection of experimental results obtained during the compression tests performed with the standard mould. The graphs plotted in this figure show the evolutions of local axial stresses as functions of the relative height

<figure><img src="image_5.png" /><figcaption>Fig. 4. Compression with the standard mould – evolution with $h - h_f$ of (a) $\sigma_{33}^I$ and $\sigma_{33}^{II}$ during the compression (SL) with $\dot{h} = 1$ mm s$^{-1}$ and for an initial SMC temperature of $T = 60$ °C, (b) $\sigma_{33}^{II}$ recorded during a (SL) compression ($\dot{h} = 1$ mm s$^{-1}$, $T = 40$ °C) and $\sigma_{33}^I$ recorded during a (SC) compression ($\dot{h} = 1$ mm s$^{-1}$, $T = 40$ °C), (c) $\sigma_{33}^{II}$ ($\dot{h} = 1$ mm s$^{-1}$, (SL)) at different initial SMC temperatures and (d) $\sigma_{33}^I$ ((SC), $T = 40$ °C) at different imposed velocities $\dot{h}$.</figcaption></figure>
$h - h_f$, $h_f$ being the final height of the considered sample. Fig. 4(a) gives typical evolutions of $\sigma_{33}^{\text{II}}$ and $\sigma_{33}^{\text{II}}$ during the “standard” (SL) compression experiments performed with an initial SMC temperature of 60 °C and with $\dot{h} = 1$ mm s$^{-1}$. The axial stress $\sigma_{33}^{\text{II}}$ first exhibits a sharp increase until reaching a peak or a threshold stress. During this first stage, air as well as styrene entrapped inside the SMC are expelled, and the wavy surface of the layers is flattened. From the peak, the SMC effectively starts to flow into the mould: stress levels exhibit a slight decrease, which is ascribed to the so-called “squish” deformation mode [33]. Then, axial stresses display a smooth increase that progressively stiffens during the compression. Fig. 4(b) compares the axial stress $\sigma_{33}^{\text{II}}$ recorded during a (SL) compression and the axial stress $\sigma_{33}^{\text{II}}$ recorded during a (SC) compression. The two compressions have been achieved with the same testing conditions, except the initial positions of the sample in the mould. This figure reveals that axial stresses are higher when the sample is initially placed on the left of the mould than at the center.

The role of the initial SMC temperature is underlined in Fig. 4(c), which shows the evolution of $\sigma_{33}^{\text{II}}$ during the standard (SL) compression of three SMC samples, whose initial temperatures were 30, 40 and 60 °C, respectively. As evident from this figure, the higher the initial temperature, the lower the threshold stress. Lastly, the influence of the compression velocity $\dot{h}$ is given in Fig. 4(d), in which the axial stresses $\sigma_{33}^{\text{II}}$ recorded during two compression experiments (SC) performed with two different axial velocities are plotted. The viscous behavior of the SMC samples is emphasized, since the stress levels strongly increase with the imposed axial velocity $\dot{h}$, and in consequence with the axial strain rate $\dot{h}/h$.

As in the case of the previous homogeneous rheometry experiments [13], it is important to note that a large scattering of ±20% on stresses was systematically recorded. This high value is commonly encountered in the study of SMC rheology [40,20] and is probably induced by the rather irregular initial thickness $h_0$ of the samples ($\Delta h_0/h_0 \approx \pm 10\%$), the inaccuracy of $L_0$ and $l_0$ ($\Delta L_0/L_0 = \Delta l_0/$

<figure><img src="image_6.png" /><figcaption>Fig. 5. Fibre-bundle volume fraction ratio $f/f_0$ as a function of $x_1$-coordinate after the compression at $\dot{h} = 1$ mm s$^{-1}$, using the standard mould and an initial SMC temperature of 60 °C (a), using the modified mould and initial SMC temperatures of 30 °C (b) and 60 °C (c), respectively. The dashed-dotted line represents the initial mean value ($f_0 \approx 17.4\%$). The dashed line represents the $\approx \pm 10\%$ deviation from the initial mean value (cf. Section 2.4).</figcaption></figure>
$l_0 \approx \pm 2\%$) and the variability of the initial fibre content (cf. previous section).

### 3.2. Fibre-bundle segregation

Once the samples had been compressed, their fibre-bundle content was analyzed following the methodology given in Section 2.3. Fig. 5 shows typical examples of the evolution along the $x_1$-coordinate of the fibre-bundle volume fraction ratio $f/f_0$, for three moulded parts: a compression with the standard mould (SL) using an initial SMC temperature of 60 °C, and two compressions using the modified mould, at an initial temperature of 30 and 60 °C. For these three examples, the axial punch velocity $\dot{h}$ was set to 1 mm s$^{-1}$. Accounting for the strong initial variability of the fibre content, it is concluded that the fibre-bundle content is approximately kept constant (i.e. close to its initial value $f_0$ of 17.4%), even if a rather weak and wavy increase of the fibre-bundle content, enhanced with the modified mould, is observed just behind the free surface of the sample.

## 4. A plug-flow shell model for the compression of SMC

### 4.1. Basic assumptions and balance equations

From the results obtained in previous studies [24,13] and those obtained in this work, the following assumptions are stated in order to build a plug flow shell model to approximate the mould filling stage:

* The SMC is seen as a one-phase continuum medium: this implies that the fibre network and the matrix have the same velocity $\underline{v}$. As a direct consequence, fibre bundle segregation phenomena will not be accounted for in the present model. Since weak segregation phenomena have been observed during the “simple” compressions investigated in this work (cf. previous subsection), this is a first rough approximation of the SMC rheology. For compressions involving more tortuous flows, such as flow into severe in-plane convergents or ribs, non negligible bundle segregation may occur [44]: for those situations, much more complex and heavy two-phase models have to be formulated [8,20,14].

* The filling phase duration is very short ($<$10 s) so that it is reasonable as a first approximation to consider the flow of the SMC without any curing of the polyester resin. Moreover, the temperature in the stack of plies may be considered as practically isothermal [3], except near the upper and lower parts of the mould.

* The SMC is assumed to be an incompressible and purely viscous body. Its stress tensor can therefore be written as

$$ \underline{\underline{\sigma}} = -p\underline{\underline{\delta}} + \underline{\underline{\sigma}}^{\text{V}}(\underline{\underline{D}}), \quad (6) $$

where $p$ is the incompressibility pressure, $\underline{\underline{\delta}}$ the identity tensor and $\underline{\underline{\sigma}}^{\text{V}}$ the viscous stress tensor, function of the strain rate tensor $\underline{\underline{D}}$:

$$ \underline{\underline{D}} = \frac{1}{2}(\underline{\underline{grad}}\underline{\underline{v}} + \underline{\underline{T}}\underline{\underline{grad}}\underline{\underline{v}}). \quad (7) $$

* External and acceleration forces will be considered as negligible quantities with respect to viscous forces generated by the deformation of the SMC.

Under such assumptions, the unique balance equations that rule the flow of the SMC, i.e. the mass and momentum balance equations, read:

$$ \left\{ \begin{array}{l} \mathrm{div} \, \underline{\underline{v}} = 0, \\ \mathrm{div} \, \underline{\underline{\sigma}} = \underline{\underline{0}}. \end{array} \right. \quad (8) $$

The plug flow kinematics constraints are now introduced, restricting the current presentation to the case where the mould cavity is a plate geometry whose midplane is contained in the $(\underline{\underline{e}}_1, \underline{\underline{e}}_2)$ plane: the extension of this 2D plate model to a 3D shell model would be straightforward. The velocity $\underline{\underline{v}}$ of a material point of SMC composite is supposed to be written:

$$ \underline{\underline{v}} = v_1(x_1, x_2)\underline{\underline{e}}_1 + v_2(x_1, x_2)\underline{\underline{e}}_2 + v_3(x_1, x_2, x_3)\underline{\underline{e}}_3, \quad (9) $$

$x_3$ being the height of the considered point from the mid-plane of the sheet. Introducing the derivative symbol $v_{i,j} = \partial v_i / \partial x_j$, the mass balance Eq. (8a) can be written as

$$ v_{1,1}(x_1, x_2) + v_{2,2}(x_1, x_2) + v_{3,3}(x_1, x_2, x_3) = 0, \quad (10) $$

from which one obtains:

$$ v_{3,33}(x_1, x_2, x_3) = 0. \quad (11) $$

Noting $h(x_1, x_2)$ the current thickness of the sheet, this last equation can be integrated accounting for the following boundary conditions $v_3(x_1, x_2, x_3 = h(x_1, x_2)/2) = \dot{h}(x_1, x_2)/2$ and $v_3(x_1, x_2, x_3 = -h/2) = -\dot{h}(x_1, x_2)/2$:

$$ v_3(x_1, x_2, x_3) = \frac{\dot{h}}{h}(x_1, x_2)x_3 = \frac{\dot{h}}{h}x_3. \quad (12) $$

Hence, the strain rate $\underline{\underline{D}}$ of a material point of SMC composite is supposed to be of the form:

$$ \underline{\underline{D}} = \begin{pmatrix} v_{1,1} & \frac{1}{2}(v_{1,2} + v_{2,1}) & 0 \\ \frac{1}{2}(v_{1,2} + v_{2,1}) & v_{2,2} & 0 \\ 0 & 0 & \frac{\dot{h}}{h} \end{pmatrix} (\underline{\underline{e}}_1, \underline{\underline{e}}_2, \underline{\underline{e}}_3) \\ = \underline{\underline{\tilde{D}}} + \frac{\dot{h}}{h} \underline{\underline{e}}_3 \otimes \underline{\underline{e}}_3, \quad (13) $$

where the symbol “--” has been introduced to distinguish 2D-unknown fields and 2D-operators in $(\underline{\underline{e}}_1, \underline{\underline{e}}_2)$ from 3D ones. As the shear components $D_{\beta 3}$ ($\beta \in \{1, 2\}$) of the strain rate tensor are constrained to zero, arbitrary reaction terms (Lagrange multiplier) must consequently be added to the stress tensor $\underline{\underline{\sigma}}$:

$$ \underline{\underline{\sigma}} = -p\underline{\underline{\delta}} + T_{\beta 3}(\underline{\underline{e}}_{\beta} \otimes \underline{\underline{e}}_{3} + \underline{\underline{e}}_{3} \otimes \underline{\underline{e}}_{\beta}) + \underline{\underline{\sigma}}^{\text{V}}, \quad \beta \in \{1, 2\}. \quad (14) $$

These assumptions are now used to simplify the balance equations of the problem. Hence, the incompressibility of the composite (8a) now reads:
$$ \tilde{\operatorname{div}} \tilde{\mathbf{v}} + \frac{\dot{h}}{h} = 0. \qquad (15) $$

Likewise, the momentum balance Eq. (8b) in $(\mathbf{e}_1, \mathbf{e}_2)$ are integrated over the thickness $h$ of the composite. Accounting for (14), they are now written as

$$ -\overline{\operatorname{grad}} P + \frac{1}{h}(\tilde{\mathbf{T}}^{+} + \tilde{\mathbf{T}}^{-}) + \overline{\operatorname{div}} \tilde{\mathbf{\Sigma}}^{\nu} = \tilde{\mathbf{0}}, \qquad (16) $$

noting

$$ P = \frac{1}{h} \int_{-h/2}^{h/2} p \, dx_3, \quad \tilde{\mathbf{\Sigma}}^{\nu} = \frac{1}{h} \int_{-h/2}^{h/2} \tilde{\mathbf{\sigma}}^{\nu} \, dx_3 \qquad (17) $$

and

$$ (\tilde{\mathbf{T}}^{+} + \tilde{\mathbf{T}}^{-}) \cdot \mathbf{e}_{\beta} = [T_{\beta 3}]_{-h/2}^{h/2} \quad \text{with } \beta \in \{1, 2\}. \qquad (18) $$

$P$ and $\tilde{\mathbf{\Sigma}}^{\nu}$ represent, respectively, the pressure and the in-plane viscous stresses averaged in the thickness of the stack, whereas $\tilde{\mathbf{T}}$ and $\tilde{\mathbf{T}}^{+}$ are the tangential stress vectors exerted on the SMC sample by the upper ($x_3 = h/2$) and lower ($x_3 = -h/2$) parts of the mould.

## 4.2. Constitutive equations

### 4.2.1. Rheology of the composite

Based on previous homogeneous rheometry experiments performed on the same type of formulation (LP606 with inhibitor [24,13]), a continuum one-phase rheological model has been proposed in [13] and will be used in the following. In order to account for both the “planar” microstructure of the fibre bundle network and the nonlinear viscous character of its deformation, the composite is seen as an incompressible power-law fluid displaying a transverse isotropy, which axis is $\mathbf{e}$, normal to the plane of the sheet. The intrinsic form of the viscous stress tensor $\tilde{\mathbf{\sigma}}^{\nu}$ is obtained using the theory of anisotropic tensor functions [7,17]. For the two particular geometries studied in this work ($\mathbf{e} = \mathbf{e}_3$) and the adopted kinematical restrictions, one obtains [13]:

$$ \tilde{\mathbf{\sigma}}^{\nu} = \eta_{\text{eq}} \frac{2}{1+2H} \left(\frac{D_{\text{eq}}}{D_0}\right)^{n-1} \tilde{\mathbf{D}} \qquad (19) $$

and

$$ \sigma_{33}^{\nu} = \eta_{\text{eq}} \frac{2H}{1+2H} \left(\frac{D_{\text{eq}}}{D_0}\right)^{n-1} \frac{\dot{h}}{h}, \qquad (20) $$

with

$$ D_{\text{eq}}^2 = \frac{2}{1+2H} \left( \tilde{\mathbf{D}} : \tilde{\mathbf{D}} + H \left( \frac{\dot{h}}{h} \right)^2 \right). \qquad (21) $$

In these equations, $H$ is a rheological function that accounts for the bundle network anisotropy, $n$ is the power law exponent of the SMC and $\eta_{\text{eq}}$ is its elongational viscosity during a simple compression deformation parallel to $\mathbf{e}$ at an equivalent strain rate $D_{\text{eq}}$ of $D_0$ (arbitrarily fixed to $1 \text{ s}^{-1}$). The power-law exponent $n$ was found to be equal

to 0.44 for the studied SMC (from $f=0.034$ to 0.234) and 0.6 for the matrix without fibre. Likewise, possible expressions for the rheological functions $H$ and $\eta_{\text{eq}}$ were given [13]:

$$ H = \frac{1 + 98f + 980f^2}{0.5 + 67f + 670f^2} - 1, \qquad (22) $$

$$ \begin{aligned} \eta_{\text{eq}} &= \eta_{\text{ps}} \left( 2 \frac{1+H}{1+2H} \right)^{-\frac{n+1}{2}} \\ &= \eta_{\text{ps}}^0 (1 + 98f + 980f^2) \left( 2 \frac{1+H}{1+2H} \right)^{-\frac{n+1}{2}}, \end{aligned} \qquad (23) $$

where the viscosities $\eta_{\text{ps}}$ and $\eta_{\text{ps}}^0$ correspond to the axial viscosities recorded during homogeneous plane strain compressions at an axial strain rate of $1 \text{ s}^{-1}$ performed on the SMC (volume fraction of fibre $f$) and the paste ($f=0$), respectively. The paste viscosity $\eta_{\text{ps}}^0$ is temperature dependent, it is supposed to follow this Arrhenius-like relation:

$$ \eta_{\text{ps}}^0 = \eta_{\text{ps}}^{00} e^{b \left( \frac{1}{T} - \frac{1}{T_0} \right)}, \qquad (24) $$

where $b = 4500 \text{ K}^{-1}$ and $\eta_{\text{ps}}^{00} = 0.18 \text{ MPa s}$ is the plane strain axial viscosity of the paste at $1 \text{ s}^{-1}$ and $T = T_0 = 296 \text{ K}$ (arbitrary choice). Even if the proposed rheological model does not describe the evolution of the in-plane orientation of fibre bundles, it is rather simple, requires few constitutive parameters and gives rather good prediction of the stress levels at the beginning of the SMC flow over a wide range of strain rates and volume fractions of fibres [13].

### 4.2.2. Friction stresses

Friction stresses are assumed to be related to the shearing of thin paste-rich boundary layers during the tangential relative motion of the bulk composite and the mould upper and lower surfaces. The geometry, the microstructure (composition, fibre-bundle content and orientation) as well as the deformation mechanisms of such boundary layers may be very complex. Following the proposition of Abrams and Castro [1], we will simply approximate friction stresses as power-law functions of the tangential SMC velocity:

$$ \tilde{\mathbf{T}}^{+} = \tilde{\mathbf{T}}^{-} = -\lambda \left( \frac{\|\tilde{\mathbf{v}}\|}{v_0} \right)^{m-1} \tilde{\mathbf{v}}, \qquad (25) $$

where $v_0$ is a reference velocity (fixed here to $1 \text{ mm s}^{-1}$), $\lambda$ a hydrodynamic friction coefficient to determine, and where $m$ is a power-law exponent related to the rheology of the boundary layers:

* $m$ should probably range between the power-law exponent of the matrix and that of the SMC at temperatures close to the mould temperature. In the absence of rheo-
logical data in such a temperature range for the tested material, it will be next arbitrarily fixed to 0.6, assuming (i) a low fibre-bundle content in the boundary layers and (ii) a weak variation of the power-law exponent with the temperature.

* Characteristic friction stresses $T_c$ arising during the compression can be roughly estimated considering (i) a characteristic in-plane velocity $v_c$ of the bulk SMC and (ii) characteristic thickness $\delta_c$ and viscosity $\mu_c$ of the boundary layers:

$$ \mu_c \approx \mu_{c0}(T, f, \text{fibre orientation}) \left(\frac{v_c}{\delta_c \dot{\gamma}_0}\right)^{m-1}. \quad (26) $$

This yields:

$$ T_c \approx \lambda \left(\frac{v_c}{v_0}\right)^{m-1} v_c \approx \mu_{c0} \left(\frac{v_c}{\delta_c \dot{\gamma}_0}\right)^{m-1} v_c, \quad (27) $$

so that

$$ \lambda \approx \frac{\mu_{c0}}{\delta_c^m} \left(\frac{v_0}{\dot{\gamma}_0}\right)^{m-1}. \quad (28) $$

In the following and for a sake of simplicity, $\lambda$ will be kept constant, close to this characteristic value.

### 4.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 $(e_1, e_2)$ of the mould 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$ is given, and the initial velocity of the composite 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 mould cavity, accounting for the following set of boundary conditions:

* $\tilde{\mathbf{v}} \cdot \tilde{\mathbf{n}} = 0$ on $\partial_v \Omega = \partial\Omega \cap \partial\Omega_M$, where $\tilde{\mathbf{n}}$ is the unit outward normal to $\Omega$ and $\partial_v \Omega$ is the part of $\partial\Omega$ which is in contact with the mould boundary $\partial\Omega_M$;

* $\tilde{\Sigma} \cdot \tilde{\mathbf{n}} = \tilde{\mathbf{0}}$ on $\partial_\sigma \Omega = \partial\Omega - \partial_v \Omega$, where $\partial_\sigma \Omega$ stands for the free surface of the composite.

### 4.4. Analytical solution for the standard mould – determination of $\lambda$

From the as-described shell model, it is possible to obtain an analytical solution for the compression with the standard mould (cf. Appendix A). Indeed, the local and mean axial stresses at the upper or lower parts of the mould $x_3 = \pm h/2$, i.e. $\Sigma_{33}(x_1)$ and $\langle\sigma_{33}\rangle = \langle\Sigma_{33}\rangle$, can be expressed as the sums of two contributions, i.e. one related to the bulk rheology of the SMC (rheo) and another related to frictions stresses (fric):

$$ \Sigma_{33}(x_1) = \eta_{\text{ps}} \underbrace{\left(\frac{|\dot{h}/h|}{D_0}\right)^{n-1} \frac{\dot{h}}{h}}_{\sigma_{33}^{\text{rheo}}} - \underbrace{\frac{2\lambda}{m+1} \left(\frac{|\dot{h}|}{v_0}\right)^{m-1} \frac{\dot{h}}{h^{m+1}} x_1^{m+1} - \left(\frac{l_0 h_0}{\alpha h}\right)^{m+1}}_{\Sigma_{33}^{\text{fric}}}, \quad (29) $$

$$ \langle\sigma_{33}\rangle = \eta_{\text{ps}} \underbrace{\left(\frac{|\dot{h}/h|}{D_0}\right)^{n-1} \frac{\dot{h}}{h}}_{\sigma_{33}^{\text{rheo}}} + \underbrace{\frac{2\lambda}{m+2} \left(\frac{|\dot{h}|}{v_0}\right)^{m-1} \frac{\dot{h}}{h^{m+1}} \left(\frac{l_0 h_0}{\alpha h}\right)^{m+1}}_{\langle\sigma_{33}^{\text{fric}}\rangle} = \eta_{\text{ps}} \left(\frac{|\dot{h}/h|}{D_0}\right)^{n-1} \frac{\dot{h}}{h}(1+\mathcal{Q}), \quad (30) $$

where $\alpha$, respectively, equals 1 and 2 for the (SL) and (SC) compression mouldings, and where the origin of the $x_1$ axis is, respectively, located at the left extremity or at the center of the mould for the (SL) and (SC) compression mouldings (see Fig. 1).

Since the problem is symmetric about $x_3 = 0$ and $h$ is small, it is expected that the gradient of $\sigma_{33}$ through the thickness is small [6]. Therefore, it is assumed that $\Sigma_{33} \approx \sigma_{33}$.

The dimensionless number $\mathcal{Q}$

$$ \mathcal{Q} = \frac{2}{(m+2)} \frac{\lambda}{\eta_{\text{ps}}} \frac{D_0^{m-1}}{v_0^{m-1}} \left|\frac{\dot{h}}{h}\right|^{m-n-1} \frac{\dot{h}}{h^2} \left(\frac{l_0 h_0}{\alpha h}\right)^{m+1}, \quad (31) $$

gauges the importance of friction stresses “fric” vs. bulk stresses “rheo”. Typically, accounting for (22), the higher the fibre content, or the lower the initial SMC temperature, or the thicker the initial charges, or the smaller their in-plane dimensions, then the lower the friction stresses. Conversely, $\mathcal{Q}$ (or friction stresses) increases during the compression. At last the role of the axial strain rate on $\mathcal{Q}$ also depends on $m-n$.

From (29) and the (SL) experiments, it is now possible to determine $\lambda$, using the difference $\Sigma_{33}(0) - \Sigma_{33}(L_0/2) = \sigma_{33}^{\text{II}} - \sigma_{33}^{\text{I}}$:

$$ \lambda = \frac{2^m (m+1) h^{m+1}}{\left(\left|\frac{\dot{h}}{h}\right|/v_0\right)^{m-1} \dot{h} L_0^{m+1}} \left[\sigma_{33}^{\text{II}} - \sigma_{33}^{\text{I}}\right]. \quad (32) $$

Please note that this form is independent of the bulk rheology of the SMC, that could be different from the viscous model used in this work.

The averaged values of $\lambda$, calculated from the (SL) mouldings, are plotted as functions of the temperature in Fig. 6. Despite the very large scattering of the results, an increase of $\lambda$ with $T$ is observed. In their recent work, Abrams and Castro [1] have assumed that the boundary layers were only made of matrix (i.e. $\mu_{c0} = \mu_{\text{matrix}0}$ in (28)). Following this hypothesis, the recorded increase of
<figure><img src="image_7.png" /><figcaption>Fig. 6. Evolution of $\lambda$ as a function of the initial temperature $T$ of the samples.</figcaption></figure>

$\lambda$ with $T$ should be correlated to a diminishing of the characteristic thickness $\delta_c$ of the lubricating layers. Notice that the other trend could also be a possible explanation. Indeed, an increase of both the thickness $\delta_c$ and the fibre-bundle content $f$ inside the layers with the initial temperature of the SMC ($\mu_{c0}(T, f)$) could also yield an increase of $\lambda$ with the initial SMC temperature, as evident from (28). These two opposite and possible scenarios only prove that the knowledge of the geometry, the microstructure and the rheology of the lubricating layers is still an open problem: further microstructural observations such as X-ray microtomography measurements (for instance) would be very helpful to obtain some answers.

### 4.5. Comparison with the mouldings

Fig. 7 compares the stresses $\sigma_{33}^I$ and $\sigma_{33}^{II}$ recorded during the (SL) or (SC) mouldings with the predictions of (29). As shown from this figure, a rather good correlation is found

<figure><img src="image_8.png" /><figcaption>Fig. 7. (SL) and (SC) compression with the standard mould – comparison between the predictions of the model and the corresponding experiments.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Fig. 8. (SL) and (SC) compression with the standard mould – evolution of the dimensionless number $\mathcal{D}$ with $h - h_f$.</figcaption></figure>

between the model and the experiments for various processing conditions, despite the simplicity of the established plug flow model and accounting for the scattering of the experimental results.

Thereby, we have plotted in Fig. 8 the evolution of the dimensionless number $\mathcal{D}$ during the compression mouldings, for all the investigated processing conditions with the standard mould. As evident from this figure, $\mathcal{D}$ ranges between 0.01 and 10, it never exceeds 10, and its mean value is $\approx 0.8$: the contribution of the SMC rheology in the overall averaged normal stresses (or on the mould closure forces) can therefore not be neglected, at least in the investigated range of processing conditions.

# 5. FE application

It has been shown that the developed shell model was able to give acceptable predictions of simple compression mouldings. This section briefly exposes its implementation in a specially designed Finite Element code, towards the simulation of more complex mouldings.

## 5.1. Resolution scheme

The flow of the composite is described by an Eulerian approach with respect to a reference surface containing $(\mathbf{e}_1, \mathbf{e}_2)$. In the present application, this reference is the mid-surface of the mould $\Omega_M$ that is meshed using triangular elements, each of them having a specific height $h$. During compression moulding, the composite fills the mould 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 mould filling applications [12,37], this was achieved with

an additional scalar variable $\chi$ describing the local volume fraction of the composite in the mould (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 transport equation:

$$ \frac{\mathrm{D}\chi}{\mathrm{D}t} = 0, \qquad (33) $$

where $\frac{\mathrm{D}}{\mathrm{D}t}$ stands for the material time derivative. Thus, the scalar field $\chi$ is added to the unknown fields $\tilde{\mathbf{v}}$ and $P$. All are dependent on time $t$ and space variables $x_1$ and $x_2$. An usual strategy is employed to solve this problem. It consists of splitting the time and the spatial discretizations as in many usual FE treatments. Thereby, the time interval $[0, T]$ of moulding is subdivided into a finite number of time steps $]t^n, t^{n+1}]$ that can have different lengths.<fnref n="1" /> Given the solution at time $t^n$, the step $n+1$ consists of (i) finding the new domain $\Omega^{n+1}$ occupied by the composite, and (ii) determining the unknown fields $\tilde{\mathbf{v}}^{n+1}, P^{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, \tilde{\mathbf{v}}^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 $\tilde{h}^{n+1}$ and $\tilde{h}^{n+1}$ for each element of the mesh, and (d) the determination of the elements that are gained by the composite between $t^n$ and $t^{n+1}$.

(ii) The problem is then subdivided into two sub-problems.

* Sub-problem (SP1), called the pressure–velocity problem. It aims at finding $\tilde{\mathbf{v}}^{n+1}$ and $P^{n+1}$, and couples the incompressibility Eq. (15) and the nonlinear momentum balance Eq. (16). This first sub-problem would be equivalent to the well-known Stokes problem in the particular case where $n = 1$.

* Sub-problem (SP2), called the free surface problem. It consists of solving the transport Eq. (33) on $\Omega^{n+1}$ to compute $\chi^{n+1}$.

## 5.2. FE formulation of the pressure–velocity problem

### 5.2.1. Weak formulation

The pressure–velocities problem SP1 is re-casted into a weak form by multiplying the Eqs. (16) and (15) by a set of test functions $\tilde{\mathbf{v}}^*$ and $P^*$ (belonging to appropriate spaces $V$ and $Q$) and by integrating over $\Omega$ applying the divergence theorem. The problem consists now of finding the set of unknowns functions $(\tilde{\mathbf{v}}, P) \in (V \times Q)$ such as $\forall (\tilde{\mathbf{v}}^*, P^*) \in (V \times Q)$:

<footnote n="1">¹ In the following, we will find the following notation for every function $f: f(\mathbf{x}, t^n) = f^n(\mathbf{x})$.</footnote>
$$
\left\{
\begin{aligned}
\int_{\Omega} \tilde{\underline{\Sigma}}^{\mathbf{V}} : \tilde{\underline{D}}^{*} dV - \int_{\Omega} \frac{\tilde{\underline{T}}^{+} + \tilde{\underline{T}}^{-}}{\tilde{h}} \cdot \tilde{\underline{V}}^{*} dV \\
- \int_{\Omega} P \operatorname{div} \tilde{\underline{V}}^{*} dV = \int_{\partial_{\sigma}\Omega} \tilde{\underline{V}}^{*} \cdot \tilde{\underline{\Sigma}} + \tilde{\underline{n}} dS, \\
\int_{\Omega} P^{*} \operatorname{div} \tilde{\underline{V}} dV = - \int_{\Omega} \frac{\tilde{h}}{\tilde{h}} P^{*} dV,
\end{aligned}
\right.
\quad (34)
$$

where the boundary conditions $\tilde{\underline{\Sigma}} \cdot \tilde{\underline{n}} = \tilde{\underline{0}}$ on $\partial_{\sigma}\Omega$ are taken into account, and where $\tilde{\underline{D}}^*$ is the strain rate tensor associated with test functions $\tilde{\underline{v}}^*$. It is important to underline that this system is highly nonlinear due to the power-laws used to determine $\tilde{\underline{\Sigma}}^V$ and $\tilde{\underline{T}}^+ + \tilde{\underline{T}}^-$ (see Section 4.2).

### 5.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 through-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-velocity problem can be readily obtained via the standard Galerkin discretization process:

$$ \begin{pmatrix} \mathbf{A}(\tilde{V}) & {}^t\mathbf{B} \\ \mathbf{B} & \mathbf{0} \end{pmatrix} \begin{pmatrix} \tilde{V} \\ \mathcal{P} \end{pmatrix} = \begin{pmatrix} F(\tilde{V}) \\ G \end{pmatrix}, \quad (35) $$

where $\tilde{V}$ is the vector of nodal velocities unknowns and $\mathcal{P}$ stands for the vector of integrated pressure nodal unknowns, $F$ is a vector that accounts for the boundary conditions, $G$ contains terms arising from the incompressibility constraint associated to $\tilde{h}/h$, and where the blocks $\mathbf{A}$ and $\mathbf{B}$ accounts for constitutive equations and incompressibility condition, respectively.

In the case of a linear system ($n = m = 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 (35) defined as

$$ \mathcal{L}_r(\tilde{U}, Q) = \frac{1}{2} \tilde{U} \mathbf{A} \tilde{U} + \frac{1}{2} r^t (\mathbf{B} \tilde{U} - G)(\mathbf{B} \tilde{U} - G) - {}^t \tilde{U} F + {}^t (\mathbf{B} \tilde{U} - G) Q, \quad (36) $$

where $r$ is a constant [36,31]. Notice that the problem (36) can be generally written into a linear system of equations:

$$ \begin{pmatrix} \mathbf{A} + r^t \mathbf{B}\mathbf{B} & {}^t \mathbf{B} \\ \mathbf{B} & \mathbf{0} \end{pmatrix} \begin{pmatrix} \tilde{V} \\ \mathcal{P} \end{pmatrix} = \begin{pmatrix} F + r^t \mathbf{B}\mathbf{G} \\ G \end{pmatrix}. \quad (37) $$

If the problem is nonlinear ($n \neq 1$ and/or $m \neq 1$), a Newton-Raphson method is used. Starting with an initial solution $(\mathcal{P}_{nr}, V_{nr})$ at the iteration $nr$, the Newton-Raphson algorithm consists of solving at the iteration $nr+1$ the residual function $R$:

$$ R(\tilde{V}_{nr+1}, \mathcal{P}_{nr+1}) = 0, \quad (38) $$

whose components write

$$ \begin{pmatrix} R_1(\tilde{V}, \mathcal{P}) \\ R_2(\tilde{V}, \mathcal{P}) \end{pmatrix} = \left( \begin{array}{c} (\mathbf{A}(\tilde{V}) + r^t \mathbf{B}\mathbf{B})\tilde{V} + {}^t \mathbf{B}\mathcal{P} - F(\tilde{V}) - r^t \mathbf{B}\mathbf{G} \\ \mathbf{B}\tilde{V} - G \end{array} \right). \quad (39) $$

For that purpose, it is thus required that each Newton-Raphson iterate satisfies the relations:

$$ \left\{
\begin{aligned}
& \left( \begin{array}{cc} \frac{\partial R_1(\tilde{V}_{nr}, \mathcal{P}_{nr})}{\partial \tilde{V}} & \frac{\partial R_1(\tilde{V}_{nr}, \mathcal{P}_{nr})}{\partial \mathcal{P}} \\ \frac{\partial R_2(\tilde{V}_{nr}, \mathcal{P}_{nr})}{\partial \tilde{V}} & \frac{\partial R_2(\tilde{V}_{nr}, \mathcal{P}_{nr})}{\partial \mathcal{P}} \end{array} \right) \begin{pmatrix} \Delta \tilde{V}_{nr} \\ \Delta \mathcal{P}_{nr} \end{pmatrix} \\
& = - \begin{pmatrix} R_1(\tilde{V}_{nr}, \mathcal{P}_{nr}) \\ R_2(\tilde{V}_{nr}, \mathcal{P}_{nr}) \end{pmatrix} \\
& \tilde{V}_{nr+1} = \tilde{V}_{nr} + \Delta \tilde{V}_{nr}, \\
& \mathcal{P}_{nr+1} = \mathcal{P}_{nr} + \Delta \mathcal{P}_{nr}, \\
& nr \leftarrow nr + 1.
\end{aligned}
\right.
\quad (40) $$

The above system can be written as

$$ \left( \begin{array}{c} \frac{\partial \mathbf{A}(\tilde{V}_{nr})}{\partial \tilde{V}} + \mathbf{A}(\tilde{V}_{nr}) + r^t \mathbf{B}\mathbf{B} - \frac{\partial F(\tilde{V}_{nr})}{\partial \tilde{V}} \quad {}^t \mathbf{B} \\ \mathbf{B} \end{array} \right) \left( \begin{array}{c} \Delta \tilde{V}_{nr} \\ \Delta \mathcal{P}_{nr} \end{array} \right) = \dots \\
= -(\mathbf{A}(\tilde{V}_{nr}) + r^t \mathbf{B}\mathbf{B})\tilde{V}_{nr} - {}^t \mathbf{B}\mathcal{P}_{nr} + F(\tilde{V}_{nr}) + r^t \mathbf{B}\mathbf{G} \\
- \mathbf{B}\tilde{V}_{nr} + G \quad (41) $$

To find the solution of this linearized system, we apply the modified Uzawa algorithm used to solve the linear problems. Thus the general algorithm for solving the nonlinear pressure-velocity problem combines the Newton-Raphson and Uzawa algorithms. In practice, the components of the matrices and vectors of the previous system are computed using a standard Gaussian quadrature integration method [45]. The solver of the matrix system given by combining the Newton-Raphson and Uzawa algorithms is based on a biconjugate gradient algorithm.

### 5.3. Free surface transport

The solution of the free surface evolution Eq. (33) is obtained through a finite element discretization in space and a characteristic based method discretization in time. The variational formulation of the Eq. (33) reads as

$$ \int_{\Omega} \frac{\mathrm{D}\chi^{n+1}}{\mathrm{D}t} \chi^{*} \mathrm{d}V = 0, \quad \forall \chi^{*} \in \Phi, \quad (42) $$

where $\chi^{n+1}$ is the solution of the variable $\chi$ field at time $t^{n+1}$ and $\chi^*$ a test function defined in an appropriate space $\Phi$. From (42), the characteristic method is applied to approximate the material derivative $\frac{\mathrm{D}}{\mathrm{D}t} \chi^{n+1}$. It consists of writing the following finite difference scheme:

$$ \frac{\mathrm{D}\chi^{n+1}}{\mathrm{D}t} = \frac{\chi^{n+1}(\tilde{\mathbf{x}}) - \chi^n(\tilde{\mathbf{x}}^n)}{\Delta t}, \quad \text{with } \tilde{\mathbf{x}}^n = \tilde{\mathbf{x}}^{n+1} - \Delta t \tilde{\mathbf{v}}^{n+1}, \quad (43) $$
where it is important to notice that $\tilde{\mathbf{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, i.e. the “foot” of the characteristic line, is the main difficulty of the method. In the current application, a first order intersection method is employed to find this position. The variational formulation finally becomes:

$$ \int_{\Omega} \chi^{n+1} \chi^* dV = \int_{\Omega} \chi^n(\tilde{\mathbf{X}}^n) \chi^* dV, \quad \forall \chi^* \in \Phi. \quad (44) $$

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

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

with $\chi^{n+1}$ the vector of nodal unknowns, $\mathbf{M}^{(\chi)}$ a mass matrix built on the left-hand side of integral (44) and $\Psi^{(\chi)}$ a vector built on the right-hand side of (44). Simple linear P1 elements are used in this discretization scheme.

### 5.4. Simulation of the mouldings

To test the capability of the developed FE code, compression mouldings identical to those achieved experimentally were simulated. For that purpose, the initial SMC charge of dimensions $L_0 \times l_0 = 170 \times 500 \text{ mm}^3$ in $(\mathbf{e}_1, \mathbf{e}_2)$ were initially located on the left side (SL) or at the center (SC) of the rectangular mould, whose dimensions in $(\mathbf{e}_1, \mathbf{e}_2)$ are $L_0 \times L_0$ ($L_0 = 500 \text{ mm}$). The shell mesh and the boundary conditions related to the simulated boundary values problems are given in Fig. 9 in the case of the (SL) mouldings.

A first set of simulations consisted of reproducing all the compressions performed with the standard mould ($h_0 \approx 8.15 \text{ mm}$). The validation of the developed code was achieved with the help of the analytical expressions (29) and (30). It was shown that a very good reproducing (<5%) of the analytical expressions of local stresses $\sigma_{33}^I$ and $\sigma_{33}^{II}$ was obtained whatever the testing conditions and

<figure><img src="image_10.png" /><figcaption>Fig. 9. Reference geometry, boundary conditions and P2+/P1 mesh used to run the simulation.</figcaption></figure>

considering the rather coarse mesh used to run the simulation (cf. Fig. 9, average size of triangle elements 20 mm).

A second type of simulation was then performed in order to model compression mouldings with the mould displaying thickness variations (cf. Fig. 1(c)). In this case, the initial uniform height of the SMC charge was set to $h_0 \approx 16 \text{ mm}$, and a variable height $h(x_1, x_2, t)$ was attributed for each finite element during the simulation. The experimental and predicted time evolutions of the local axial stresses recorded during the compression using this mould ($\dot{h} = 1 \text{ mm s}^{-1}$, initial SMC temperature = 40 and 60 °C) have been plotted in Fig. 10. Even if the simulated stresses at the center of the mould slightly underestimate the experimental data, numerical results give acceptable reproducing

<figure><img src="image_11.png" /><figcaption>Fig. 10. (SL) compression with the modified mould – comparison between the experimental and simulated time evolutions of the axial stresses $\sigma_{33}^I$ and $\sigma_{33}^{II}$: (a) $\dot{h} = 1 \text{ mm s}^{-1}$, initial SMC temperature $T = 40 \text{ °C}$, (b) $\dot{h} = 1 \text{ mm s}^{-1}$, initial SMC temperature $T = 60 \text{ °C}$.</figcaption></figure>
of experimental data, considering the very simple rheological and frictional models used to run the simulation.

## 6. Conclusion

Compression moulding experiments have been performed on an instrumented industrial press under various processing conditions and using a standard SMC formulation. Experimental results have underlined the influence of the position of the stack charge in the mould, its initial temperature as well as the axial punch velocity on recorded local normal stress levels. A destructive technique was then used in order to measure the fibre-bundle content (averaged in the thickness) inside the moulded parts. Small fibre-bundle segregations were observed just behind the flow front of the composite, i.e. ≈10–15% of the initial fibre-bundle content. They seemed to be enhanced using the mould with thickness convergents, but were further assumed to be negligible accounting for the rather strong variability of the fibre-bundle content of the initial SMC. From the experimental results obtained previously in the literature and those gained in this work, a one-phase plug flow shell model was then formulated, extending the plug flow model initially proposed in [6]. A rather simple model was used for the bulk rheology of SMC, that was assumed to behave like an incompressible and purely viscous (in power-law) body displaying transverse isotropy. The constitutive parameters of the model were entirely determined in a previous study over a wide range of strain rates, temperatures and fibre contents, using a SMC formulation very close to that used in this work. Friction stresses between the SMC and the mould were also modelled in a very simple manner, reducing a very complex situation to a phenomenological hydrodynamic friction law with constant friction parameters. From this plug flow model, analytical solution were derived in the case of the standard mould, underlining the role of both the SMC rheology and frictions on the overall stress levels. The unique unknown constitutive parameter, i.e. the friction coefficient, was determined from the analytical solution and the (SL) mouldings. It was found to be temperature dependent. The local physical reasons of such a trend still remain unclear and demand further analysis of the microstructure of the lubricating layers. Nonetheless, the predictions of the as-identified simple plug flow model were satisfactory for the compression mouldings with the standard mould and for a quite wide range of processing conditions. There from, it has been clearly shown that the bulk rheology of the composite should not be neglected with respect to friction stresses.

For that purpose, a specific finite element code able to account for both bulk and friction stresses in the plug flow shell model was then developed within an Eulerian framework and using (i) a nonlinear mixed pressure–velocity formulation of the flow problem and (ii) a volume fraction of SMC advected with the SMC velocity to model the free surface of the composite. Compared with the experiments,

predictions of the simulated compression mouldings using the mould with thickness convergents were found to be acceptable.

Efforts are now focusing in implementing in the plug flow model more sophisticated rheological model based on upscaling approaches and able to account for the strong coupling between rheology and bundle orientation [43,42,15,38,22,23].

## Acknowledgements

The authors are grateful for the support of the Région Rhône–Alpes and the industrial partners Schneider Electric, Renault, Groupe Mecelec-Inoplast and Renault Trucks. The authors also wish to thank Compositec for the help provided during the moulding experiments.

## Appendix A. Analytical solutions for the mean and local axial stresses for the standard mouldings (SL) and (SC)

Accounting for the unidirectional properties of the flow problem through the standard rectangular mould, the velocity fields $\mathbf{v}$ and the strain rate tensor $\underline{\mathbf{D}}$ can be written, respectively:

$$ 
\begin{aligned}
\mathbf{v} &= v_1(x_1)\mathbf{e}_1 + \frac{\dot{h}}{h}x_3\mathbf{e}_3 \quad \text{and} \\
\underline{\mathbf{D}} &= v_{1,1}\mathbf{e}_1 \otimes \mathbf{e}_1 + \frac{\dot{h}}{h}\mathbf{e}_3 \otimes \mathbf{e}_3.
\end{aligned}
\quad (46) $$

The mass balance equation, see Eq. (15), writes:

$$ v_{1,1} + \frac{\dot{h}}{h} = 0. \quad (47) $$

Noting that $\dot{h}/h = c^{\text{ste}}$ and accounting for the boundary conditions $v_1(x_1 = 0) = 0$, the solution of this equation is simply:

$$ v_1(x_1) = -\frac{\dot{h}}{h}x_1. \quad (48) $$

This velocity field when introduced into the rheological and friction laws, see Eqs. (19) and (25), yields the following system of momentum balance equations, see Eq. (16):

$$ 
\begin{cases}
-P_{,1} + \frac{1}{h}(\tilde{T}^{+} + \tilde{T}^{-})_1 = 0, \\
-P_{,2} = 0,
\end{cases}
\quad (49) $$

whose solution $P(x_1)$ is given integrating (49), accounting for the boundary condition at the flow front ($x_1 = L = (l_0h_0)/h$) for the stress state ($\tilde{\Sigma}_{11} = -P(x_1 = L) + \tilde{\Sigma}_{11}^v = 0$):

$$ P(x_1) = -\eta_{\text{ps}} \left( \frac{1}{1+H} \right) \frac{1}{D_0^{n-1}} \left| \frac{\dot{h}}{h} \right|^{n-1} \frac{\dot{h}}{h} + \dots \\ + \left( \frac{2\lambda}{m+1} \right) \frac{1}{v_0^{m-1}} \left( \frac{1}{h} \right) \left| \frac{\dot{h}}{h} \right| \left| \frac{\dot{h}}{h} \right| x_1^{m+1} - \left( \frac{l_0 h_0}{\alpha h} \right)^{m+1}, \quad (50) $$

where $\alpha = 1$ or 2 for the (SL) and (SC) mouldings, respectively.
It is also possible to calculate the mean axial stress $\langle \sigma_{33} \rangle$:

$$ 
\begin{aligned}
\langle \sigma_{33} \rangle &= \frac{F_3}{S} = \frac{1}{L} \int_{0}^{L} \sigma_{33} \, dx_1 \\
&= -\frac{1}{L} \int_{0}^{L} p(x_1, x_3) \, dx_1 + \sigma_{33}^v, \quad \forall x_3,
\end{aligned}
\quad (51) $$

where $F_3$ is the axial compression force and $S = L_0L$ the current surface of the compounds inside the mould and where $\sigma_{33}^v$ has a constant value, see Eq. (20). The mean axial stress writes in particular for $x_3 = h/2$:

$$ \langle \sigma_{33} \rangle = -\frac{1}{L} \int_{0}^{L} p(x_1, h/2) \, dx_1 + \sigma_{33}^{v}. \quad (52) $$

The global form of momentum balance equation along $e_3$ allows to write:

$$ \int_{0}^{L} p(x_1, x_3) \, dx_1 = \int_{0}^{L} p(x_1, h/2) \, dx_1, $$

so that

$$ 
\begin{aligned}
& \frac{1}{h} \int_{-h/2}^{h/2} \left[ \int_{0}^{L} p(x_1, x_3) \, dx_1 \right] \, dx_3 \\
&= \frac{1}{h} \int_{-h/2}^{h/2} \left[ \int_{0}^{L} p(x_1, h/2) \, dx_1 \right] \, dx_3 \\
\Leftrightarrow & \int_{0}^{L} P(x_1) \, dx_1 = \int_{0}^{L} p(x_1, x_3 = h/2) \, dx_1,
\end{aligned}
\quad (53) $$

which leads to the following expression of $\langle \sigma_{33} \rangle$ accounting for (50):

$$ 
\begin{aligned}
\langle \sigma_{33} \rangle = & \frac{2\lambda}{m+2} \left(\frac{|\dot{h}|}{v_0}\right)^{m-1} \frac{\dot{h}}{h^{m+1}} \left(\frac{l_0 h_0}{\alpha h}\right)^{m+1} \\
& + \eta_{\text{ps}} \left(\frac{|\dot{h}/h|}{D_0}\right)^{n-1} \frac{\dot{h}}{h}.
\end{aligned}
\quad (54) $$

In the same way, it is also possible to determine $\Sigma_{33}(x_1)$, which is given by

$$ 
\begin{aligned}
\Sigma_{33}(x_1) &= -\frac{1}{h} \int_{-h/2}^{h/2} p(x_1, x_3) \, dx_3 + \frac{1}{h} \int_{-h/2}^{h/2} \sigma_{33}^v \, dx_3 \\
\Leftrightarrow \Sigma_{33}(x_1) &= -P(x_1) + \sigma_{33}^v \\
\Leftrightarrow \Sigma_{33}(x_1) &= -\frac{2\lambda}{m+1} \left(\frac{|\dot{h}|}{v_0}\right)^{m-1} \frac{\dot{h}}{h^{m+1}} \\
& \quad \times x_1^{m+1} - \left(\frac{l_0 h_0}{\alpha h}\right)^{m+1} + \eta_{\text{ps}} \left(\frac{|\dot{h}/h|}{D_0}\right)^{n-1} \frac{\dot{h}}{h}.
\end{aligned}
\quad (55) $$

## References

[1] Abrams L, Castro J. Predicting molding forces during sheet molding compound (SMC) compression molding. I: Model development. Polym Compos 2003;24:291-303.

[2] Advani S, Tucker C. A numerical simulation of short fiber orientation in compression molding. Polym Compos 1990;11(3):164-73.

[3] Barone M, Caulk D. The effect of deformation and thermoset cure on heat conduction in a chopped-fiber reinforced polyester during compression molding. Int J Heat Mass Transfer 1979;22(7):1021-31.

[4] Barone M, Caulk D. Kinematics of flow in sheet molding compound. Polym Compos 1985;6(2):105-9.

[5] 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 1987;26:185-206.

[6] Barone MR, Caulk DA. A model for the flow of a chopped fiber reinforced polymer compound in compression molding. J Appl Mech 1986;53(191):361-70.

[7] Boelher JP. Applications of tensor functions in solid mechanics. New York: Springer; 1987.

[8] Bowen RM. Continuum physics III: Mixtures and EM field theories, vol. III. London: Academic Press; 1976.

[9] Castro J, Griffith R. Sheet molding compound compression-molding flow. Polym Eng Sci 1989;29(10):632-8.

[10] Castro J, Tomlinson G. Predicting molding forces in SMC compression molding. Polym Eng Sci 1990;30(24):1568-73.

[11] Christensen S, Hutchinson B, Sun E, Osswald T, Davis B. Fiber-matrix separation in ribbed SMC and BMC parts. In: ANTEC, editor. Proc. ann. tech. conf. - ANTEC'. Brookfield (CT): Soc. of Plastics Engineers; 1997. p. 782-7.

[12] Dhatt G, Gao DN, Médal M, Song R. Remplissage et solidification par éléments finis. Euro Finite Element Rev 1992;1:309-32.

[13] Dumont P, Orgéas L, Le Corre S, Favier D. Anisotropic viscous behaviour of sheet molding compounds (SMC) during compression molding. Int J Plast 2003;19(4):625-46.

[14] Dumont P, Le Corre S, Orgéas L, Favier D, Gaborit C, Lory P. Finite element implementation of a two-phase model for compression molding of composites. Euro Finite Element Rev 2005;14:883-900.

[15] Gibson A, Toll S. Mechanics of the squeeze flow of planar fibre suspensions. J Non-Newtonian Fluid Mech 1999;82:1-24.

[16] Hieber C, Shen S. A finite-element/finite-difference simulation of the injection-molding filling process. J Non-Newtonian Fluid Mech 1980;7:1-32.

[17] Hull B, Rogers T, Spencer A. Theoretical analysis of forming flows of continuous-fibre-resin systems. In: Flow and rheology in polymer composites manufacturing. Composite materials series, vol. 10. Amsterdam: Elsevier; 1994. p. 203-56.

[18] Kim J, Shiau Y, Lee L, Im Y. Compression molding simulation of chopped fiber reinforced polymeric composites in plate-rib type geometry. Polym Compos 1992;13(2):97-107.

[19] Kim K-T, Jeong J-H, Im Y-T. Effect of molding parameters on compression molded sheet molding compounds parts. J Mater Process Technol 1997;67:105-11.

[20] Kim S, Im Y. Three-dimensional thermo-viscoplastic analysis of compression molding of sheet molding compounds with fiber volume fraction prediction. J Mater Process Technol 1997;63:631-6.

[21] Kotsikos G, Gibson AG. Investigation of the squeeze flow behaviour of sheet moulding compounds (SMC). Composites Part A 1998; 29(12):1569-77.

[22] Le Corre S, Caillerie D, Orgéas L, Favier D. Behavior of a net of fibers linked by viscous interactions: theory and mechanical properties. J Mech Phys Solids 2004;52:395-421.

[23] Le Corre S, Dumont P, Orgéas L, Favier D. Rheology of highly concentrated planar fiber suspensions. J Rheol 2005;49:1029-58.

[24] Le Corre S, Orgéas L, Favier D, Tourabi A, Maazouz A, Venet C. Shear and compression behaviour of sheet molding compounds. Compos Sci Technol 2002;62(4):571-7.

[25] Lee CC, Folgar F, Tucker CL. Simulation of compression molding for fiber-reinforced thermosetting polymers. J Eng Ind 1984;106: 114-25.

[26] Lee C, Tucker C. Flow and heat transfer in compression mold filling. J Non-Newtonian Fluid Mech 1987;24(3):245-64.

[27] Lee LJ, Marker LF, Griffith RM. The rheology and mold flow of polyester sheet molding compound. Polym Compos 1981;2(4): 209-18.
[28] Leterrier Y, G'Sell C. Multilayer plug flow modeling of the fast stamping process for a polypropylene/glass fiber composite. Polym Compos 1996;17(2):231-40.

[29] Lin C-M, Weng C-I, Ho C-T. Anisotropy in sheet molding compounds during compression molding. Polym Compos 1997; 18(5):613-22.

[30] Lin C-M, Weng C-I. Simulation of compression molding for sheet molding compound considering the anisotropic effect. Polym Compos 1999;20(1):98-113.

[31] Liu W, Xu S. A new improved Uzawa method for finite element solution of Stokes problem. Comput Mech 2001;27:305-10.

[32] Michaeli W, Mahlke M, Osswald T, Nolke M. Simulation of the flow in SMC. Kunststoffe German Plastics 1990;80(6):31-3.

[33] Odenberger P, Andersson H, Lundström T. Experimental flow-front visualisation in compression moulding of SMC. Composites Part A 2004;35:1125-34.

[34] Osswald TA, Tucker CL. Compression mold filling simulation for non-planar parts. Int Polym Process 1990;2:79-87.

[35] Osswald T, Tucker C. Compression mold filling simulation for non-planar parts. Int J Polym Process 1990;5(2):79-87.

[36] Robichaud MP, Tanguy PA, Fortin M. An iterative implementation of the Uzawa algorithm for 3D-fluid flow problems. Int J Numer Meth Fluids 1990;10:429-42.

[37] Scardovelli R, Zaleski S. Direct numerical simulation of free-surface and interfacial flow. Ann Rev Fluid Mech 1999;31:567-603.

[38] Servais C, Luciani A, Månson J-AE. Fiber-fiber interaction in concentrated suspensions: dispersed fiber bundles. J Rheol 1999; 43(4):1005-18.

[39] Servais C, Luciani A, Månson J-AE. Squeeze flow of concentrated long fibre suspensions: experiments and model. J Non-Newtonian Fluid Mech 2002;104:165-84.

[40] Silva-Nieto R, Fisher B. Rheological characterization of unsatured polyester resin sheet molding compound. Polym Eng Sci 1981;21(8): 499-506.

[41] Silva-Nieto R, Fisher B, Birley A. Predicting mold flow for unsatured polyester resin sheet molding compound. Polym Compos 1980;1(1): 14-23.

[42] Souloumiac B, Vincent M. Steady shear viscosity of short fibre suspensions in thermoplastics. Rheol Acta 1998;37:289-98.

[43] Toll S, Månson J-AE. Dynamics of a planar concentrated suspension with non-hydrodynamic interaction. J Rheol 1994;38(4):985-97.

[44] Xu J, Kim J, Ho T, Lee L. Compression molding of sheet molding compounds in plate-rib type geometry. Polym Compos 1993;14(1): 51-8.

[45] Zienkiewicz OC, Taylor RL. The finite element method. The basis, vol. 1, 5th ed. Oxford, UK: Butterworth-Heinemann; 2000.