# Squeeze flow in heterogeneous unidirectional discontinuous viscous prepreg laminates: Experimental measurement and 3D modeling

Grégoire Sorba, Christophe Binetruy, Adrien Leygue, Sebastien

Comas-Cardona

► To cite this version:

Grégoire Sorba, Christophe Binetruy, Adrien Leygue, Sebastien Comas-Cardona. Squeeze flow in heterogeneous unidirectional discontinuous viscous prepreg laminates: Experimental measurement and 3D modeling. Composites Part A: Applied Science and Manufacturing, 2017, 103, pp.196-207. 10.1016/j.compositesa.2017.10.007. hal-04325824

HAL Id: hal-04325824

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

Submitted on 17 Oct 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Squeeze flow in heterogeneous unidirectional discontinuous viscous prepreg laminates: Experimental measurement and 3D modeling

G. Sorba *, C. Binetruy, A. Leygue, S. Comas-Cardona

GeM - Research Institute of Civil Engineering and Mechanics, UMR 6183, CNRS, Ecole Centrale de Nantes, Université de Nantes, 1 rue de la Noé, 44321 Nantes, France

The freedom of composite design can be improved by combining continuous and discontinuous prepregs. The forming of a pre heated blank made of optimally oriented and distributed discontinuous prepreg plies may lead to unacceptable defects such as in plane and out of plane wrinkles, sliding of plies over long distance, rotation of adjacent plies, bending of fiber induced by transverse squeeze flow and finally to inappropriate and inefficient fiber distribution. This arises because the individual discontinuous plies are free to move and deform in the mold during the forming step. First, this paper presents some exper iments conducted to identify the behavior of a stack of unidirectional discontinuous viscous prepregs subjected to through thickness compression. Then a model based on a heterogeneous anisotropic fluid approach is gradually developed in accordance with the experimental findings. It is shown that the var ious observed phenomena are retrieved by the numerical model and that the predicted values are in good agreement with measurements, but also that it requires to be solved in 3D with a relatively fine mesh in the thickness to provide good results.

## 1. Introduction

One of the key advantages in designing fiber reinforced composites is the freedom of design. Polymer composites can be shaped in very complex forms to meet design requirements start ing from a prepreg material that is a flat sheet which combines fibers and thermoplastic or uncured thermoset matrix. To get a structural 3D part at a production rate suitable for the automotive industry, one of the available processes is to form the initially pla nar prepreg into a final 3D part through thermo stamping. In order to keep the design freedom, the main concern is to control fiber orientation and reduce the remaining defects below an acceptable level, while keeping the target cycle time. The design freedom can be further developed in thermo forming of thermoplastic or ther mosetting materials combining continuous and discontinuous pre prepregs. Discontinuous prepregs are patches used to locally strengthen the part where necessary. Structural analysis helps to identify locations where a part needs to be particularly reinforced with continuous fiber reinforced composites and where the part is less mechanically loaded, allowing cheaper composites to be used.

A multi thickness/multi material blank made of thermoplastic or thermosetting patches can be manufactured according to this design procedure as shown in Fig. 1. The blank is preheated and formed under press to get the tailored composite [1]. Fig. 1b shows the flat blank designed to be formed to obtain a curved hat shaped framing component. A discontinuous woven prepreg placed on the right hand side of the blank is initially aligned with the component axis. After forming to obtain the hat shape, this patch experienced a rotation as shown in Fig. 1c. No in plane shearing is observed because the forming in this region of the com ponent consists in squeezing and folding the discontinuous pre preg. On the left hand side of the tailored blank, a discontinuous unidirectional prepreg underwent large transformations as seen in Fig. 1c. Fibers are no longer straight and the ply dimensions are not controlled due to the large squeeze flow. Pure sliding of dis continuous prepreg are also observed in vertical zones as a conse quence of the mold closing. The main advantages of this technology are the use of noble materials only where necessary, allowing for fewer scraps while achieving final weight reduction. Beyond these many advantages, there are still some difficulties to overcome as shown in the example in Fig. 1. The pre heated blank is not held inside the tool, which lets the individual discontinuous plies free to move and deform in the mold during the forming step. This can lead to unacceptable defects such as in plane and out of plane wrinkles, sliding of patches over long distance, rotation of

<footnote>* Corresponding author.
E-mail addresses: gregoire.sorba@ec-nantes.fr (G. Sorba), christophe.binetruy@ec-nantes.fr (C. Binetruy), adrien.leygue@ec-nantes.fr (A. Leygue), sebastien.comas@ec-nantes.fr (S. Comas-Cardona).</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1. (a) UD prepreg deformation and sliding, (b) flat multi-material blank before forming, (c) woven prepreg rotation [1].</figcaption></figure>

adjacent plies, bending of fibers induced by transverse squeeze flow and finally to inappropriate and inefficient fiber distribution [1]. All these mechanisms are the result of stresses (compression, tension and shear) which build up during the forming process.

The interlaminar shear effects play a significant role in the forming of multi layered composite parts, especially when discon tinuous patches made of continuous fibers are introduced in the lay up. A resin rich layer exists in between prepreg plies, where the two resin rich layers on either side can represent about 5% of the total thickness [2], with the thickness depending on the tech nology used. It has an impact on the behavior of the plies [3 5]. Due to its important lubrication function, the thermoplastic or thermosetting matrix decreases interly friction and therefore facilitates large amount of slippage, deformation and rotation of discontinuous plies. Interply slip helps to alleviate the in plane stresses due to compression, and fiber waviness. Fiber buckling can then be avoided.

The production of quality parts without extensive experimental investigations is needed. It is therefore important to predict where issues may occur during forming and how they can be alleviated. Designing the process through a trial and error procedure may lead to an acceptable product, but it always brings additional labor and tooling costs, machine time and scraps. Due to the complexity of the problem, the numerical simulation of the forming of complex blanks is necessary. Numerical tools can simulate the production processes and provide to the subsequent Finite Element Analysis (FEA) a material scale description of the as manufactured part. The in plane loading and slippage of a UD prepreg ply has been modeled in [6]. This situation can arise when a prepreg patch is dragged by the tooling in the downward direction during the deep drawing of a blank.

This paper investigates the through thickness compression of laminated fiber aligned prepreg patches and its numerical model ing. The flow of the viscous matrix between and along the fibers is not considered. The continuous fibers in each ply prevent any extensional flow in the fiber direction [7], although a too low vis cosity might lead to bleeding [8] [5]. Owing to the wide range of prepregs, to the various and complex mechanisms reported in [1], a detailed analysis is required to capture each mechanism indi vidually. To do so, Section 2 presents dedicated experiments with instrumented specimens to measure large deformation including changes in fiber orientation and inter ply rotation. Experimental findings will be used to build a preliminary multi layer heteroge neous 3D model presented in Section 3. Section 4 addresses its numerical solution. Section 5 presents some numerical tests linked to experiments run in Section 2. The numerical predictions are compared to experimental findings to refine the proposed model and make it more realistic.

## 2. Experimental observations

Compaction experiments are run on uncured unidirectional prepregs.

The material used in the experiments is a unidirectional carbon fiber reinforced epoxy prepreg (Hexply M21/35.5 M%/268/T700GC from Hexcel) of 56.9% nominal fiber volume fraction. As it is a material of industrial grade, the given properties are subjected to statistical variations. The experiments are conducted on uncured thermoset unidirectional prepregs because the principal desired characteristic being the material to be a linear viscous fluid. As a monomer is very likely to stay in its linear domain regardless of the rate of shear, this assumption is considered reasonable. Further more, considering that the fillers added in the material are of small dimension compared to the ply, the resin rich layer can be consid ered an homogenized fluid which viscosity is the one provided by the manufacturer. Individual plies have a nominal cured thickness of 0.262 mm. 150 × 150 mm² samples with 3 stacking sequences (UD, cross ply and angle ply) were considered. The ply edges are left unconstrained which allows them to expand and rotate freely. Thin copper threads of 0.1 mm in diameter were added at some interfaces between plies to probe the interface flow. An array of parallel threads aligned with adjacent UD plies is inserted at the interface. To form a grid, a second array of parallel threads can be placed at another interface between adjacent UD plies aligned in the second stacking direction. The two arrays of tracers are not inserted at the same interface as they would significantly increase the interface thickness and possibly interact with each other. A visual description of the configuration is presented in Fig. 2.

The dimension of the tracers is higher than the average thick ness of the resin interface, however it does not have an impact on the flow as their stiffness is very low compared to the viscosity of the matrix. Furthermore they will be convected by the resin in the same fashion as the fibers are, considering that the flow is induced by the through thickness compression. While their exact position in the thickness cannot be controlled, they will penetrate the sheet of fibers with which they are aligned. They therefore remain an interesting solution to inspect the configuration of the fibers.

As these tracers are dragged by the viscous polymer at the inter face between two plies they provide useful experimental data to be compared to model predictions. After placing these tracers in the prepreg assembly, mold plates are pre heated at 260 °C then the stack is debulked in the press under 100 kN during five minutes. Both the initial and final 2D geometry of the grid are obtained from a high resolution X ray computed tomography (CT) scanner (X Radia Carl Zeiss). The high contrast between the copper threads and the carbon fiber polymer composite allows for a fast scanning with a very high accuracy.

The application of a compaction force produces a pressure gra dient within the material that induces transverse flow. This trans verse flow allows fibers to spread sideways under normal compressive load. Because the fibers are inextensible, the sample ply cannot expand along their axis but resin can be squeezed out in that direction if the resin viscosity is low enough. The viscosity evolves with the temperature and the degree of cure, but as the thermo stamping is a fast process there is no evolution of these
<figure><img src="image_3.png" /><figcaption>Fig. 2. Configuration of the copper threads. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

properties. In order to perform experiments at a constant viscosity, preliminary tests are carried out to determine the range of temperatures and time scales where the degree of cure does not evolve and the viscosity is high enough to generate pure squeeze flow without resin bleeding in the fiber direction. The results of these tests are a maximum temperature of 30 min processing time at a temperature of 80 °C maximum to avoid curing of the resin and resin bleeding in the fiber direction. Higher temperatures are close to the transition from squeezing to bleeding flow as observed in [5]. During the actual tests the time limit is never reached as thermo stamping is a very fast process, lasting only a few seconds. Specimens were squeezed in the thickness direction with a slow monotonic loading applied by a press in the range of 250–400 kN.

Table 1 summarizes the testing process parameters.

## 2.1. Compression of a $[0]_6$ stack

Fig. 3 shows the initial and final geometry of the grid for the $[0]_6$ stack. Initial fiber orientation is represented by the red arrow. The array of tracers orthogonal to the fiber direction was in the mid plane, the parallel one on both sides of the fourth ply, counted from the bottom side of the stack. 60 × 60 mm² scans are taken in the center of the sample. The high transverse squeeze flow perpendicular to the fibers induced the breakage of copper threads aligned in that direction as can be seen in the blue circle, while the parallel ones remain intact (Fig. 3b). This provides the evidence that a reaction tension stress develops in fibers to prevent the flow along their direction. The vertical tracers experienced in plane bending, as can be seen in the green rectangle, that is more pronounced close to the lateral edges resulting from the squeeze flow kinematics.

## 2.2. Compression of a $[0/90]_6$ stack

Fig. 4 presents the initial and final geometry of the grid placed in the mid plane of a $[0/90]_6$ cross ply stack. Initial fiber orienta

tions are represented by the red arrows. The horizontal array of tracers was inserted between the fourth and the fifth plies, the second array oriented in the vertical direction was inserted between the eight and the ninth plies. As revealed by these images, the orthogonality of the grid is kept as shown by the blue wedges, but unlike the UD stack, all tracers underwent some bending as can be seen in the green rectangle. The bending and spreading of the composite layers are clearly visible on the composite layers as indicated in Fig. 5. The grid dilated of 15% on average as high lighted by the yellow double headed arrows in Fig. 4, which proves there is actually a viscous flow between the plies. While orthogonally oriented plies impede flow in each direction within composite plies, they do not prevent the flow at their interface. This important observation is consistent with the existence of pure viscous resin layers in prepregs, as mentioned in Section 1. Even for prepregs with high fiber volume fraction, there is flow at ply interfaces.

## 2.3. Compression of a $[+30/ 30]_3$ stack

Fig. 6 shows results obtained for a stack of plies with $[+30/ 30]_3$ lay up. The first array of tracers was in the mid plane, the second one on both sides of the fourth ply, counted from the bottom side of the stack. The initial grid has a shape of a rhombus with a 60° angle corresponding to the total inclination of adja cent plies. A marker crossing one copper thread, highlighted in the blue square, is added to the grid for a more accurate comparison between the undeformed and deformed grids. Measurements made on three repeats indicate that the thin copper threads rotated 8 10° towards [+45/ 45] and kept their straightness in the center of the sample. This experimental finding is confirmed on real composites observations reported in [9] for a composite substitute made of carbon fiber reinforced syrup. The grid dilated of 13 16% of its original spacing for the three repeats. The $[0/45]_3$ stack of 1.7 mm initial thickness, squeezed under 350 kN up to 1.2 mm in thickness exhibited the highest rotation, typically 14°, for the same squeeze rate than the one applied to the $[+30/ 30]_3$ lay up.

The following important conclusions are drawn from this experimental campaign:

* There is always a pure resin layer between the plies in prepreg laminates.

* A reaction tension stress develops in fibers to prevent the flow along their direction.

* There is a complex flow at the plies interface.

<table><caption>**Table 1**

Testing process parameters summary.</caption><tr><td>Material</td><td>Hexcel Hexply M21/35.5M%/268/ T700GC</td></tr><tr><td>Fiber volume fraction</td><td>56.9 %</td></tr><tr><td>Dimensions</td><td>150 × 150 mm²</td></tr><tr><td>Debulking</td><td>260 °C under 100 kN during 5 min</td></tr><tr><td>Maximum testing time to avoid curing</td><td>80 °C during 30 min</td></tr><tr><td>Force range</td><td>250–400 kN - monotonic</td></tr></table>
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** CT-scan of the [0]$_{6}$ stack, 400 kN compressive load. Red arrows represent the initial fiber orientation. (a) Before squeeze flow. Initial thickness: 1.6 mm. Initial grid space: 10 mm. (b) After squeeze flow. Final thickness: 1.3 mm. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** CT-scan of the [0/90]$_{6}$ stack, 400 kN compressive loading. Average expansion of the grid of 15%. (a) Before squeeze flow. Initial thickness: 3.15 mm. Initial grid space: 10 mm. (b) After squeeze flow. Final thickness: 2.38 mm. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

* The kinematics of the copper tracers are not the same as the kinematics of the plies, although the tracers are a good indicator of the movement of the plies.

* The observed phenomena that the model should be able to retrieve are:
  Squeeze flow kinematics perpendicular to the fiber direction, inextensibility in the parallel direction, no bleeding of resin. Rotation and in plane bending of the plies.

## 2.4. Copper threads kinematics

Fig. 3 presents the CT scans of the [0]$_{6}$ stack before and after compression. The vertical threads barrel, as can be expected from a stack of UD under compression, whereas the horizontal threads break, confirming that the copper threads are mostly at the inter face as no carbon fibers were broken during the compression test.

Fig. 4 presents the CT scans of the [0/90]$_{6}$ stack before and after compression, and Fig. 5 presents the photograph of the stack after compression. The Figures show that all threads barrel and spread in the same fashion as the fibers.

It can therefore be concluded that while the copper threads can be thicker than the interface layer they are mostly located in the

inter layer zone, and their kinematics are representative of the kinematics of the fibers.

## 3. Numerical modeling

Different levels of analysis are possible, including the fiber and tow levels, ply by ply analysis in which the prepreg plies are modeled as distinct bodies or composite models that use effective properties. Here the ply by ply analysis is chosen because of the occurrence of large relative displacement of plies allowed by the resin interply. In addition, the primary focus is on fiber aligned prepreg where tows are undistinguishable.

There are two main approaches in composite forming simulations: the geometrical/kinematic approach [10] and the mechanical/Finite Element (FE) approach [11].

The first one is very fast and adapted to simple geometries where only fiber reorientation is considered for preliminary design purposes. However this approach cannot be used here since the effects of inter ply interactions cannot be included. In addition, fiber reorientation can be generated by the transverse squeeze flow in multi axial discontinuous laminates [9]. It is also observed that fibers can deform so that the initial rectangular sample becomes
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Photograph of the $[0/90]_6$ sample after compression. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

barrel shaped [12]. These mechanisms cannot be taken into account in a pure kinematic approach.

Mechanical/Finite Element (FE) simulations are based on solving the equilibrium of the full structure. They can include complex material models and boundary conditions such as tool part friction. Constitutive relations connecting stresses to the kinematics are required in the mechanical approach. They can be generally classified under the assumption that the prepreg material is made of a fluid or a solid. Since our purpose is to model forming processes where discontinuous fiber aligned prepregs can extensively slide, rotate and deform due to transverse squeeze flow, the following focuses on fluid models which seem the most appropriate framework.

The Transversally Isotropic Fluid (TIF) model has been developed on the basis of the elastic ideal fiber reinforced model proposed by Pipkin [13] and Spencer [14]. This work was adapted by Rogers [11] to the case of viscous composites in order to simulate the forming process of a stack of composite plies. The melted resin is considered a viscous incompressible fluid, and the presence of the fibers is taken into account by adding an inextensibility constraint in the fiber direction, leading to an homogenized viscous fluid, transversely isotropic, incompressible and inextensible. It is a model based on physical properties such as fiber volume fraction and resin viscosity. It is attractive because it reduces the need for extensive characterisation campaign. As the rate of deformation during the process is slow enough, non linear dynamical effects are not taken into account, meaning that the Reynolds number $Re \ll 0$. Based on that hypothesis, the homogenized fluid is considered in a succession of quasi static states. The simulation consists in computing the instantaneous velocity, then updating the geometry of the stack of plies and the fiber orientation before going on to the next time step.

Lamers et al. later extended this model with elastic components that allow the model to incorporate fiber stresses and an elastic fabric shear response in drape predictions [15], by assigning a finite stiffness to the fibers and using compressible matrix proper ties, in a way that it is not necessary to use Lagrange multipliers in the FE calculations. In this modified model the interlayer was not modeled as a separate continuum, but was treated with a simple viscous slip law expressed in terms of the velocity differences between adjacent plies.

The following general assumptions are used in the modified TIF model used in this paper:

* Acceleration effects are neglected, i.e. $Re \ll 0$.

* In the plies:
  No distinction is made between fibers and matrix in any way in the continuum model.
  Fibers introduce an inextensible direction in the composite plies.
  Fibers convect with the melted polymer during the flow.
  As a result, the composite ply is modeled as a homogenized linear viscous fluid reinforced by inextensible fibers.
  The strain rate dependency is not included.

* In the interfaces:
  A linear viscous fluid is inserted to model the interface.

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** CT-scan of the $[+30/ - 30]_3$ stack, 350 kN compressive loading. Tracers in the mid-plane. Average expansion of the grid between 13% and 16% on 3 repeats. (a) Before squeeze flow. Initial thickness: 1.6 mm. Initial grid space: 10 mm. (b) After squeeze flow. Final thickness: 1.25 mm. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
### 3.1. Equilibrium of the fluid

As in any mechanical problem the first equation to be solved is equilibrium:

$$ \mathrm{div}(\underline{\sigma}) + f_b \underline{0}, \qquad (1) $$

with $\underline{\sigma}$ the stress tensor and $f_b$ the body forces. As the problem is quasi static this equation will be solved at each time step to compute the instantaneous velocity, pressure and fiber tension.

### 3.2. Constitutive model of viscous uniaxial composite

Rogers [11] proposed the following behavior law:

$$ \underline{\sigma} \underline{D} : \dot{\underline{\epsilon}} \underline{P} \underline{1} + T \underline{a} \otimes \underline{a}, \qquad (2) $$

with $\dot{\underline{\epsilon}}$ the rate of strain tensor:

$$ \dot{\underline{\epsilon}} = \frac{1}{2} \frac{\partial \underline{v}^T}{\partial \underline{x}} + \frac{\partial \underline{v}}{\partial \underline{x}}, \qquad (3) $$

$\underline{v}$ the velocity, $P$ the pressure, $T$ the tension in the fibers, $\underline{a}$ the unitary fiber orientation:

$$ \underline{a} \begin{bmatrix} \cos(\theta) \\ \sin(\theta) \\ 0 \end{bmatrix}, \qquad (4) $$

with $\theta$ the angle between the fibers and the horizontal axis and $\underline{D}$ the transversely isotropic viscosity tensor:

$$ \underline{D} = 2\eta_L \underline{1} + 2(\eta_L \eta_T) \underline{a} \otimes \underline{a}, \qquad (5) $$

with $\eta_L$ the longitudinal viscosity and $\eta_T$ the transverse viscosity.

Eq. (5) introduces a fluid of viscosity $\eta_L$ along the fiber direction and $\eta_T$ in the other principal directions, corresponding to the desired behavior for the TIF model.

### 3.3. Kinematic constraints

The modeled fiber aligned prepreg material has a viscous behavior coupled with an incompressibility and an inextensibility in the fiber direction, as discussed by Rogers [11]. The incompressibility and the inextensibility in the fiber direction are respectively modeled as:

$$ \dot{\underline{\epsilon}} : \underline{1} \underline{0}, \qquad (6) $$

with $\underline{1}$ the identity matrix, and

$$ \dot{\underline{\epsilon}} : (\underline{a} \otimes \underline{a}) \underline{0}. \qquad (7) $$

### 3.4. Explicit update

The resolution algorithm consists of considering the fluid in a succession of quasi static states, and to update the geometry and the material properties before proceeding to the next timestep. In this work however, the update will not be performed and only the solution at the first timestep will be computed, as the focus is on the tendency predicted by the model and not the actual full computation.

#### 3.4.1. Mesh position

Once the velocity is known, the position of the mesh can be updated using an explicit scheme:

$$ \underline{x}_{i+1} = \underline{x}_i + \Delta t \underline{v}_i, \qquad (8) $$

with $\underline{x}$ the nodes' position and $\Delta t$ the timestep.

The inextensibility constraint being enforced on the instantaneous velocity, it is not satisfied as soon as the mesh is modified. Therefore, a special attention is paid on the quality of the computation due to the use of an explicit scheme instead of an implicit one, and the stability of the method should be verified. It also constrains the number of time steps in order to keep numerical stability.

#### 3.4.2. Fiber reorientation

As the fibers move with the fluid, their orientation is convected by the flow in the ply. Considering the inextensibility of the fibers, the evolution equation of $\underline{a}$ is:

$$ \frac{D\underline{a}}{Dt} (\nabla \underline{v})\underline{a}, \qquad (9) $$

with $\nabla \underline{v}$ the velocity gradient, as discussed by Rogers [11].

Injecting (4) into (9) and expressing $\dot{\underline{a}}$ and $\underline{a}$ in terms of $\dot{\theta}$ and $\theta$ leads to

$$ \left\{ \begin{array}{l} \dot{\theta} \sin \theta \quad \frac{\partial v_1}{\partial x} \cos \theta + \frac{\partial v_1}{\partial y} \sin \theta, \qquad (a) \\ \dot{\theta} \cos \theta \quad \frac{\partial v_2}{\partial x} \cos \theta + \frac{\partial v_2}{\partial y} \sin \theta, \qquad (b) \end{array} \right. \qquad (10) $$

and combining Eqs. (10a) and (10b) leads to the evolution equation of $\theta$:

$$ \dot{\theta} \cos^2 \theta \left[ \left( \frac{\partial v_1}{\partial x} \right)^2 + \left( \frac{\partial v_2}{\partial x} \right)^2 \right] + \sin^2 \theta \left[ \left( \frac{\partial v_1}{\partial y} \right)^2 + \left( \frac{\partial v_2}{\partial y} \right)^2 \right] + 2 \cos \theta \sin \theta \frac{\partial v_1}{\partial x} \frac{\partial v_1}{\partial y} \frac{\partial v_2}{\partial x} \frac{\partial v_2}{\partial y}. \qquad (11) $$

Eq. (9) can be injected into an explicit update scheme in order to obtain the evolution equation of the fiber orientation:

$$ \theta_{n+1} = \theta_n + \Delta t \cdot \left[ \left( \frac{\partial v}{\partial y} \frac{\partial u}{\partial x} \right) \cos \theta \sin \theta + \frac{\partial u}{\partial y} \cos^2 \theta + \frac{\partial v}{\partial x} \sin^2 \theta \right] \qquad (12) $$

## 4. Numerical simulation

The Finite Element Method (FEM) is used to compute an approximate solution to the equilibrium under constraints presented in Section 3. Bradaigh [16] used a mixed penalty method to compute the velocity and then reconstruct the tension in the fibers when needed. In this paper, a mixed formulation is used to obtain the velocity and the tension in a single step, to avoid the dependency of the solution on the value of a penalty parameter and to provide more accurate results.

### 4.1. Weak form

The calculation of the weak formulation of the governing equations is similar to the one developed by Bradaigh in [16]. The dissipated power is expressed as

$$ E_k \frac{1}{2} \int_{\Omega} \dot{\underline{\epsilon}} : \underline{D} : \dot{\underline{\epsilon}} d\Omega \qquad (13) $$

and the power of the external forces is expressed as

$$ W \int_{\Omega} \underline{v} f_b d\Omega + \int_{\Gamma} \underline{v} f_s d\Gamma, \qquad (14) $$

with $f_b$ the body forces and $f_s$ the surface forces.

To enforce the constraints Lagrange multipliers are used. A Lagrange multiplier $P$ is associated with the incompressibility:

$$ P(\underline{\underline{\epsilon}} : \underline{1}) \underline{0}, \qquad (15) $$
with P the pressure and a Lagrange multiplier T is associated with the inextensibility:

$$T[\dot{\underline{\epsilon}} : (a \otimes a)] = 0, \quad (16)$$

with T the fiber tension.

The principle of stationarity gives:

$$\delta E_k = \delta \left\{ \frac{1}{2} \int_{\Omega} \dot{\underline{\epsilon}} : \underline{D} \equiv \underline{\dot{\epsilon}} d\Omega + \int_{\Omega} P(\underline{\dot{\epsilon}} : \underline{1}) d\Omega + \int_{\Omega} T[\underline{\dot{\epsilon}} : (a \otimes a)] d\Omega \right. \\ \left. \int_{\Omega} \underline{\nu} f_b d\Omega + \int_{\Gamma} \underline{\nu} f_s d\Gamma \right\} = 0. \quad (17)$$

Taking all possible variations of Eq. (17) yields:

$$\left\{ \begin{array}{ll} \delta \underline{\nu}^T \rightarrow \int_{\Omega} \delta \underline{\nu}^T \left[ \underline{\nabla}^{sT} \left( \underline{D} : \underline{\dot{\epsilon}} \right) + (a \otimes a) \underline{\nabla} T + \underline{\nabla} P f_b \right] d\Omega & \\ \int_{\Gamma} \delta \underline{\nu}^T f_s d\Gamma = 0 & (a) \\ \\ \delta P \rightarrow \int_{\Omega} \delta P(\underline{\dot{\epsilon}} : \underline{1}) d\Omega = 0 & (b) \\ \\ \delta T \rightarrow \int_{\Omega} \delta T[\underline{\dot{\epsilon}} : (a \otimes a)] d\Omega = 0 & (c) \end{array} \right. \quad (18)$$

## 4.2. Discretisation

The velocity, the fiber tension and the pressure are discretised using different interpolation functions. This allows to use a higher order of interpolation functions for the velocity than for the fiber tension and the pressure, avoiding numerical issues during the solving of the system. This matter is addressed by Hughes in [17] and will be discussed in Section 4.3. This leads to the forms of the unknown presented in (19).

$$\left\{ \begin{array}{ll} \underline{\nu}(\underline{x}) \approx \underline{\tilde{\nu}}(\underline{x}) & \underline{N}^{\nu}(\underline{x}) \underline{\hat{\nu}} \\ P(\underline{x}) \approx \underline{\tilde{P}}(\underline{x}) & \underline{N}^{p}(\underline{x}) \cdot \underline{\hat{P}} \\ T(\underline{x}) \approx \underline{\tilde{T}}(\underline{x}) & \underline{N}^{T}(\underline{x}) \cdot \underline{\hat{T}} \end{array} \right. \delta \underline{\nu}(\underline{x}) \approx \delta \underline{\tilde{\nu}}(\underline{x}) \underline{N}^{\nu}(\underline{x}) \delta \underline{\hat{\nu}} \quad (19)$$

where $\underline{N}^\nu$, $\underline{N}^p$ and $\underline{N}^T$ are the interpolation functions for the velocity, the pressure and the fiber tension respectively, $\underline{\nu}$, $\underline{\tilde{P}}$ and $\underline{\tilde{T}}$ the numerically approximated fields for the velocity, the pressure and the fiber tension respectively and $\underline{\hat{\nu}}$, $\underline{\hat{P}}$ and $\underline{\hat{T}}$ the nodal values of the velocity, the pressure and the fiber tension respectively.

Introducing these interpolations into Eqs. (18a), (18b) and (18c) and noting that these equations must be verified for all $\delta \underline{\hat{\nu}}$, $\delta \underline{\hat{P}}$ and $\delta \underline{\hat{T}}$, the complete equations system is obtained as in (20). It is important to notice that the system (20) is typical of a saddle point problem, a particular category much discussed in the literature [18,19].

$$\left[ \begin{matrix} \underline{K}^{\nu} \\ (\underline{K}^{p})^{T} \\ (\underline{K}^{t})^{T} \end{matrix} \right] \left[ \begin{matrix} \underline{\underline{\nu}} \\ \underline{0} \\ \underline{0} \end{matrix} \right] \left[ \begin{matrix} \underline{f} \\ \underline{0} \\ \underline{0} \end{matrix} \right] \quad (20)$$

$\underline{K}^\nu$, $\underline{K}^p$, $\underline{K}^t$ and $\underline{f}$ are defined as:

$$\left\{ \begin{array}{ll} K_{ij}^{\nu} & \int_{\Omega} (\underline{B}_{i}^{\nu})^{T} \underline{D} \underline{B}_{j}^{\nu} d\Omega \\ K_{ij}^{p} & \int_{\Omega} (\underline{B}_{i}^{\nu})^{T} \underline{N}_{j}^{p} d\Omega \\ K_{ij}^{t} & \int_{\Omega} (\underline{B}_{i}^{\nu})^{T} \underline{a} \underline{N}_{j}^{t} d\Omega \\ f_{i} & \int_{\Omega} (\underline{N}_{i}^{\nu})^{T} \underline{f}_{\nu} d\Omega \quad \int_{\Gamma} (\underline{N}_{i}^{\nu})^{T} \underline{f}_{s} d\Gamma \end{array} \right. \quad (21)$$

## 4.3. Choice of shape functions

It is discussed in [17] that the shape functions for velocity should be of a higher order than the shape functions for the

Lagrange multiplier associated with the incompressibility constraint in order to avoid locking, and this assumption was made for the inextensibility constraint as well. This consideration leads to the choice of the H27 8 8 elements, containing 27 velocity nodes and 8 pressure/tension nodes with continuous interpolation for all the unknown throughout all the domain. Tri quadratic shape functions are therefore used for interpolating the velocity whereas tri linear shape functions are used for interpolating the pressure and the fiber tension.

The stability of elements for problems of incompressible fluid has been discussed, and a mathematical stability criterion known as the Ladyzhenskaya Babuška Brezzi (LBB) condition was presented in [20 22]. Although it was suggested that the tests used for the stability of elements in an incompressible problem may not suffice for the inextensible problem in [23], it is assumed that an element that verifies the LBB condition for an incompressible problem also verifies it for an incompressible and inextensible problem, as it is shown in Eqs. (15) and (16) that their form is very similar.

## 4.4. Computational cost

As a mixed formulation is used to compute the pressure and the fiber tension the resulting system is rather large. To give the reader an indication of the computational cost, a cube of 24 by 24 by 49 elements, resulting in 775 597 dof, is solved in 5 h using one core of an AMD Opteron 6204 @ 3.3 GHz and required around 350 Gb of memory using Matlab.

## 5. Results

Previous works reported only 2D results. This study will show that a 2D model cannot capture all the distinctive features reviewed in the introduction. Instead, the 3D results obtained show that they are closer to the experimental observations. More over, the results are computed using a mixed formulation, allowing the computation of the Lagrange multipliers and the velocity field with more accuracy. Unless said otherwise only the results for the first timestep are presented. When velocity results are presented, it refers to the instantaneous velocity at the beginning of the timestep.

All computations were run with the resin viscosity $\eta$ 100 Pa.s, longitudinal viscosity of the composite $\eta_L$ 379 Pa.s and transverse viscosity of the composite $\eta_T$ 602 Pa.s. Those values were obtained using a semi empirical model developed by Christensen in [24] with a fiber volume fraction of 56.9%.

The transverse viscosity is higher than the longitudinal one as given by Christensen's model. This means that making the fibers slide along them is easier than making them slide cross wise, which is true when considering the sliding of one ply with respect to the other due to the fibers of one layer falling into the gaps of the other layer.

The first simulations run are not quantitatively representative of the real stack of plies, the first objective being the verification that the model retrieves qualitatively the observed phenomena. Therefore the geometry is not meant to be quantitatively compared against experimental samples.

### 5.1. Cross ply stack without resin layer interface

A stack of two composite plies modeled as two layers of TIF oriented at 0°(bottom ply) and 90°(top ply) is considered, with no sliding at the interface. The height of the cube is 5 mm and the thickness of each TIF layer is 2.5 mm. A slip condition was applied on the top side, a no slip condition on the bottom side and conti
<figure><img src="image_8.png" /></figure>

(a) Velocity magnitude expressed in mm.s⁻¹. The bottom ply is oriented at 0° and the top one at 90°.

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

(b) Deformed interface corresponding to the case in Figure 7a. A multiplicative factor of 3 was applied to the deformation to improve visualization. The fiber tension is expressed in Pa.

<caption>**Fig. 7.** Results for the cross-ply stack without resin layer interface.</caption>

nuity of all the fields is imposed throughout the whole domain. The number of elements in the plane is 14 × 14 and 29 in the thickness, leading to an aspect ratio of 14/29. As seen on Fig. 7a each layer cannot be elongated in the direction of their fibers, resulting in a blocking interface that contains fibers in the two directions and a numerical peak in the fiber tension as it is not allowing any movement as can be seen in Fig. 7b.

The only remaining movement allowed for the interface is the out of plane bending, as can be seen in Fig. 7, with the maximum displacement induced by the bending being equal to 75% of the displacement induced by the compression. This behavior is not consistent with the hypothesis of plane interface made in previous works [25]. There is no indication from the experiments that such a mechanism occurs.

Moreover, it is observed experimentally that the composite layers can slide with respect to each other, this model is therefore not able to reproduce the phenomenon observed. As the 2D model is unable to describe the through thickness behavior, it will behave the same and is therefore not suitable.

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

(a) Velocity magnitude expressed in mm.s⁻¹. The bottom ply is oriented at 0° and the top one at 90°.

## 5.2. Cross ply stack with explicit resin layer interface

To improve the 3D TIF model, layers of fluid were added to the previous cross ply model, two layers at the tool/composite interface and one layer at the composite/composite interface. The mid dle fluid layer will allow the plies to slide, while the side fluid layers allow the sliding between the plies and the tools. Modeling the interply layer as a 3D fluid layer will allow to better capture the interaction between the two plies than the lubrication theory or a friction parameter, although the computational cost will greatly increase.

Those fluid layers are modeled as Stokes fluid with viscosity η equal to the viscosity of the resin.

The cube is 5 mm thick, the thickness of each TIF layer is 1 mm and the thickness of each fluid layer is 1 mm (3 mm in total). No slip conditions were applied on the top and the bottom sides. The number of elements in the plane is 14 × 14 and 29 in the thickness, leading to an aspect ratio of 14/29. Fig. 8a shows the velocity component aligned with the fiber direction of the bottom unidirec

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

(b) Deformed interface. The bottom interface correponds to the tool/composite interface. The top interface if the composite/composite interface. A multiplicative factor of 3 was applied to the deformation to improve visualization. The fiber tension is expressed in Pa.

<caption>**Fig. 8.** Results for the cross-ply stack with explicit resin layer interface. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</caption>
<figure><img src="image_12.png" /><figcaption>**Fig. 9.** Comparison between the simulation with and without the explicit Stokes layer interface. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

tional layer. Due to the viscous fluid layer each ply is now able to slide with respect to each other and spread transversely to the direction of their fibers. As a consequence the numerical peak in the fiber tension is greatly reduced as shown in the top interface in Fig. 8. The bottom interface in Fig. 8 represents the tool/composite interface, whereas the top one is the composite/composite interface. As the interfaces between the plies and the fluid are now less constrained they can elongate and they no longer bend during the transverse squeeze flow.

As can be seen in Fig. 8a the behavior of the fluid layer is entirely 3D and cannot be reduced to 2D without a great loss of information, proving once more the necessity of the 3D model.

Fig. 9 presents a comparison of the fiber tension and the velocity between the simulations with and without the explicit resin layer. Fig. 9a shows that the presence of the Stokes layer dampens greatly the peak in the fiber tension. Fig. 9b shows that the inex tensibility is verified in each ply, and that the Stokes layer accomodates the velocity of each adjacent ply, effectively relieving the coupling between them.

### 5.3. [30, 30] stack with resin layer interface

A stack of two composite plies modeled as two layers of TIF oriented at 30° and 30°, with layers of fluid is considered. The height

of the cube is 5 mm, the thickness of each TIF layer is 1 mm and the thickness of each fluid layer is 1 mm. The number of elements in the plane is 14 × 14 and 29 in the thickness, leading to an aspect ratio of 14/29.

Experimental observations revealed that alternate layers inclined at an angle different from 90°tend to rotate when they are subjected to compression. The computed velocity field for the first time step in [25] suggested the same behavior. Some parametric studies achieved here show that the magnitude of rotation depends on the enforced thickness reduction and the rheological properties of both the composite and resin layers. It is then mandatory to compute the full compression through all the timesteps to determine the final deformation of the composite. The ply re orientation is predicted correctly only by the combination of the TIF model and the presence of the fluid layer between the plies.

As mentioned in Section 2 solid rotation is not the only occur ring phenomenon, as can be seen on Fig. 10 the ply is also subject to in plane shear and in plane bending, although this deformation is somehow hindered by the presence of the fibers.

Fig. 10 presents the decomposition of the velocity field, projected in the plane, taken in the middle of the bottom ply. On Fig. 10a the spreading, i.e. the velocity orthogonal to the fiber direction, is shown. The further from the center, the higher the spreading is, which is consistent with a squeeze flow orthogonal

<figure><img src="image_13.png" /><figcaption>**Fig. 10.** Decomposition of the velocity field in the middle of the lower ply. The velocity is expressed in mm s⁻¹. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
to the fiber direction. On Fig. 10b the rest of the velocity field is shown, and it can be seen that it is a combination of solid body rotation and in plane bending.

## 5.4. Single UD ply with thickness variation of the fluid layer

Due to the variability of the prepreg material, it is very common to see thickness variation of the resin layer. In these situations it is observed that a discontinuous ply will be washed away when subjected to compression.

A single unidirectional ply modeled as a layer of TIF oriented at 0° with two layers of fluid is considered. The base dimensions of the cube are 3 mm, the thickness of the TIF and resin layer is 1 mm and the extra height on the right hand side of the computational domain is 1 mm. The inclination is along the direction of the fibers. The number of elements in the plane is 14 × 14 and 29 in the thick ness, leading to an aspect ratio of 14/29.

As presented on Fig. 11 the TIF layer slides along the fiber direc tion when subjected to transverse compression. The initial unde formed mesh is displayed in the figures to better show the displacement of the TIF layer.

Fig. 11b shows the deformed mesh for the tenth time step, proving that the method does not suffer from numerical instabili ties in the velocity field. This is also true for the pressure and the fiber tension fields. It is recalled here that neither stabilisation method nor implicit update is used to obtain these results.

## 5.5. Comparison with experiment

The objective is now to quantitatively compare the simulation results with the experimental observations. Due to the dimensions of the prepreg material, a finer mesh is required in the thickness.

The modeled stack consists of two layers of TIF of 0.786 mm each surrounded by layers of Stokes of 0.05 mm each, correspond ing approximately to 6% of the composite layer thickness. The ini tial orientation of the bottom layer is 30° and the initial orientation of the top layer is 30°. Due to the use of FEM, the entire stack can not be computed with a sufficient accuracy on the results due to the computational cost. The simulations were therefore carried out on an area of interest in the middle of the stack, and using only one timestep covering the whole process. In this case the area of interest is a cube of dimension 1.7 mm. The number of elements in the plane is 24 × 24, and 49 in the thickness to obtain a suffi cient accuracy without degrading the conditioning of the system, leading to a total of 28,224 elements and 775,597 degrees of free dom. The maximal aspect ratio is therefore 49/24.

Fig. 12 presents the deformation of the middle of the TIF layer, were it is shown that there is a solid body rotation of 6.34° between the two plies. Comparison between the numerical and experimen tal values is given in Table 2.

The differences between the experimental and the numerical values can arise from various origins:

<figure><img src="image_14.png" /><figcaption>**Fig. 11.** Velocity component along the fiber direction. The horizontal velocity is expressed in mm s⁻¹. The plain geometry represents the undeformed shape. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 12.** Deformed geometry of the middle of the TIF layers. The angle is expressed in °. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
<table><caption>**Table 2**
Comparison between the experimental and the numerical values for the (30, -30) stack.</caption><thead><tr><th></th><th>Experimental</th><th>Numerical</th></tr></thead><tbody><tr><td>Solid body rotation</td><td>8-10°</td><td>6.34°</td></tr><tr><td>Spreading</td><td>13-16%</td><td>25%</td></tr></tbody></table>

<table><caption>**Table 3**
Results of the sensitivity study.</caption><thead><tr><th></th><th></th><th></th></tr></thead><tbody><tr><td>Variation of the thickness of the Stokes layers</td><td>+10%</td><td>-10%</td></tr><tr><td>Variation of the rotation</td><td>-22%</td><td>+26%</td></tr><tr><td>Variation of the spreading</td><td>-21%</td><td>+14%</td></tr></tbody></table>

* The assumption was made that the pure resin layer represented 6% of the TIF thickness, and that it is of constant thickness. This is not always verified in industrial prepregs, and the thickness might be not constant in the same sample and might change from one sample to another.

* The resin viscosity was taken from the manufacturer datasheet, and the TIF viscosities calculated using Christensen's semi empirical model presented in [24]. The actual homogenized viscosity of the ply might vary from one sample to another, and Christensen's model gives accurate results considering the material is homogeneous at the scale of the part, which it may not be at the scale of the ply.

* The experimental values are calculated from the displacement of the copper wires which position in the thickness of the inter facial layer is not perfectly controlled. It is worth to note they are located in between the composite layers and not inside them, which could change the results.

* The values are retrieved from a simulation run on a volume of interest in the middle of the sample. A sensitivity study is presented in Section 5.6, which combined with the fact that the material used is of industrial grade lead to variations of the solution that might be important.

## 5.6. Sensitivity study

To assess the impact of the parameters of the model on the results, a sensitivity study has been performed. The viscosity was changed by ±20% and the thickness of the Stokes layers by ±10%.

As the fluid is supposed linear viscous, and Dirichlet boundary conditions are applied, the change of viscosity does not impact the velocity field, therefore computed rotation and spreading are not impacted.

Table 3 presents the impact of the thickness of the Stokes layers on the rotation and the spreading, in which it can be seen that the results are very sensitive to the thickness. This sensitivity explains the differences between the experimental and the numerical values, as the material used is of industrial grade and therefore contains inherent structural variation, such as the average and the spatial distribution of the thickness.

The sensitivity also proves the relevance of the inter ply layers, both in terms of suitable modeling but also in terms of accurate initial geometrical data.

## 6. Conclusion

A stack of discontinuous composite plies can show numerous behaviors, among them a variety of phenomenon of interest regarding the forming process:

* Plies sliding with respect to each other.

* Solid body rotation between plies.

* In plane bending and shearing within a ply.

* Solid body translation.

This paper focused on the through thickness compression of unidirectional discontinuous viscous prepreg laminates. A series of experiments and characterisation has been presented to better understand the flow kinematics induced by the compaction tests. Thanks to thin tracers inserted between plies, the experiments revealed the complexity of the flow at the inter ply. It was demonstrated that the available 2D TIF model was insufficient to correctly predict all these phenomena; a purpose built 3D model was there fore developed and presented. Finally it was proved that in order to reproduce the experiments this model needs a fluid Stokes layer between the TIF layers leading to a heterogeneous 3D model. The latter needs to be solved using a mixed formulation to provide sat isfactory results. The proposed model is finally validated against one experiment using the same range of load rate, ply and pure resin layer thickness and stacking sequence. This model can cap ture the geometrical size effects and use very few material proper ties that can be obtained by direct characterisation or by means of a model.

The computational cost is very high (3 h using 16 cores of two Opteron 6328 @ 3.2 GHz, and 350 Gb of memory), however this is not an obstacle as this model is meant to give reference results, understand the physics observed in experiments and interpret the experimental findings. It is a baseline to elaborate a simplified, more efficient model, particularly needed when simulating the actual parts.

## Acknowledgement

The authors would like to thank the CETIM for the useful dis cussions about mechanisms observed in forming of blanks made of discontinuous prepregs.

## References

[1] Schell JSU, Amory L, Guillon D, Chinesta F, Cueto E, Abisset-Chavanne E. Movement of patches during thermoforming: experiment and simulation. In: AIP conference proceedings, vol. 1769. AIP Publishing; 2016. p. 170032.

[2] Cogswell FN. Thermoplastic aromatic polymer composites: a study of the structure, processing and properties of carbon fibre reinforced polyetheretherketone and related materials. Elsevier; 2013.

[3] McGuinness GB, Brádaigh CMÓ. Development of rheological models for forming flows and picture-frame shear testing of fabric reinforced thermoplastic sheets. J Non-Newton Fluid Mech 1997;73(1):1-28.

[4] McEntee SP, Brádaigh CMÓ. Large deformation finite element modelling of single-curvature composite sheet forming with tool contact. Compos Part A: Appl Sci Manuf 1998;29(1):207-13.

[5] Nixon-Pearson OJ, Belnoue JP-H, Ivanov DS, Potter KD, Hallett SR. An experimental investigation of the consolidation behaviour of uncured prepregs under processing conditions. J Compos Mater 2016:1911-24. https://doi.org/10.1177/0021998316665681.

[6] Sorba G, Binetruy C, Chinesta F, Cueto E, Abisset-Chavanne E. In-plane shearing of a ud prepreg modeled as transversely isotropic fluid: Comparison between continuous and discontinuous fiber tension approaches. AIP conference proceedings, vol. 1769. AIP Publishing; 2016. p. 170008.

[7] Cogswell FN. The experience of thermoplastic structural composites during processing. Compos Manuf 1991;2(3-4):208-16. https://doi.org/10.1016/0956-7143(91)90142-4.

[8] Shuler SF, Advani SG. Transverse squeeze flow of concentrated aligned fibers in viscous fluids. J Non-Newton Fluid Mech 1996;65(1):47-74. https://doi.org/10.1016/0377-0257(96)01440-1.

[9] Jones RS, Roberts RW. Ply re-orientation in compression. Compos Manuf 1991;2(3-4):259-66.

[10] Hancock SG, Potter KD. The use of kinematic drape modelling to inform the hand lay-up of complex composite components using woven reinforcements. Compos Part A: Appl Sci Manuf 2006;37(3):413-22. https://doi.org/10.1016/j.compositesa.2005.05.044.

[11] Rogers TG. Rheological characterization of anisotropic materials. Composites 1989;20(1):21-7. https://doi.org/10.1016/0010-4361(89)90677-0.

[12] Wheeler AB, Jones RS. Numerical simulation of fibre reorientation in the consolidation of a continuous fibre composite material. Compos Manuf 1995;6 (3-4):263-8.
[13] Pipkin AC, Rogers TG. Plane deformations of incompressible fiber-reinforced materials. J Appl Mech 1971;38(3):634. https://doi.org/10.1115/1.3408866

[14] Spencer AJM. Deformations of fibre-reinforced materials. Oxford science research papers. Clarendon Press; 1972.

[15] Lamers EAD, Akkerman R, Wijskamp S. Fibre orientation modelling for rubber press forming of thermoplastic laminates. Int J Form Process 2003;6(3-4):443-63.

[16] Brádaigh CMÓ. Sheet forming of composite materials. In: Flow and rheology in polymer composites manufacturing. Elsevier Science b.v.; 1994. Edition.

[17] Hughes TJR. Finite element analysis of incompressible viscous flows by the penalty function formulation. J Comput Phys 1974;30:1-60.

[18] Fortin M, Glowinski R. Augmented Lagrangian methods: applications to the numerical solution of boundary-value problems (studies in mathematics and its applications). Elsevier Science Ltd; 1983.

[19] Arrow K, Hurwics L, Uzawa H. Programming in linear spaces. In: Arrow K, Hurwics L, Uzawa H, editors. Studies in linear and nonlinear programming. Stanford (California): Stanford University Press; 1958.

[20] Ladyzhenskaya OA, Salkow H. The mathematical theory of viscous incompressible flow. 2nd ed. New York: Gordon and Breach; 1969.

[21] Babuška I. The finite element method with lagrangian multipliers. Numer Math 1973;20(3):179-92. https://doi.org/10.1007/BF01436561

[22] Brezzi F, existence On the uniqueness and approximation of saddle-point problems arising from lagrangian multipliers. ESAIM: Math Model Numer Anal 1974;8(R2):129-51.

[23] Šimáček P, Kaliakin VN, Pipes RB. Pathologies associated with the numerical analysis of hyper-anisotropic materials. Int J Numer Meth Eng 1993;36(20):3487-508. https://doi.org/10.1002/nme.1620362006.

[24] Christensen RM. Effective viscous flow properties for fiber suspensions under concentrated conditions. J Rheol 1993;37(1):103-21. https://doi.org/10.1122/1.550459.

[25] Ghnatios C, Abisset-Chavanne E, Binetruy C, Chinesta F, Advani S. 3d modeling of squeeze flow of multiaxial laminates. J Non-Newton Fluid Mech 2016;234:188-200.