Henry Fool | Betriebsgebäude | Puerto Rico

3D modeling of squeeze flow of multiaxial laminates

3D modeling of squeeze flow of multiaxial laminates

Accepted Manuscript 3D Modeling of Squeeze Flow of Multiaxial Laminates Ch. Ghnatios, E. Abisset-Chavanne, Ch. Binetruy, F. Chinesta, S. Advani PII: ...

4MB Sizes 0 Downloads 0 Views

Recommend Documents

Squeeze flow of a dry foam
Squeeze flow of a dry compressible foam between two circular disks is studied by means of a lubrication analysis. Two si

3D modeling technique for the analysis of delaminated composite laminates
A shell/3D modeling technique was developed for which a local three-dimensional solid finite element model is used only

Optimization of chocolate 3D printing by correlating thermal and flow properties with 3D structure modeling
•Correlate parameter setup based on thermal and flow properties of chocolate.•Additional of MgST aid in lubricating choc

Modeling fatigue behavior of quasi-isotropic laminates
Fatigue behavior of quasi-isotropic carbon/epoxy laminates were studied, using six stacking sequences. Experimental resu

Reverse squeeze film flow in a continous flow system
The continuous-flow squeeze film apparatus has been adapted to permit flow to take place in either direction. This simul

Modeling the impact behavior of AD85 ceramic under multiaxial loading
This paper presents an advanced constitutive model to describe the complex behavior of ceramic materials under impact lo

Model-free inversion of squeeze-flow rheometer data
A method for extracting rheological data from squeeze-flow tests is proposed. The analysis, based upon the lubrication a

Squeeze flow testing of glass mat thermoplastic material
The anisotropic flow behaviour of glass fibre/polypropylene glass mat thermoplastic (GMT) has been investigated using ax

Accepted Manuscript

3D Modeling of Squeeze Flow of Multiaxial Laminates Ch. Ghnatios, E. Abisset-Chavanne, Ch. Binetruy, F. Chinesta, S. Advani PII: DOI: Reference:

S0377-0257(16)30080-5 10.1016/j.jnnfm.2016.06.004 JNNFM 3798

To appear in:

Journal of Non-Newtonian Fluid Mechanics

Received date: Revised date: Accepted date:

30 October 2015 30 May 2016 2 June 2016

Please cite this article as: Ch. Ghnatios, E. Abisset-Chavanne, Ch. Binetruy, F. Chinesta, S. Advani, 3D Modeling of Squeeze Flow of Multiaxial Laminates, Journal of Non-Newtonian Fluid Mechanics (2016), doi: 10.1016/j.jnnfm.2016.06.004

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT Highligths Proposal of a penalized transversally isotropic fluid In-plane-out-of-plane separated representation Solution of a 3D problem with a 2D complexity Consideration of extremely fine thickness descriptions Discussion of many different laminates

AC

CE

PT

ED

M

AN US

CR IP T

• • • • •

Preprint submitted to Elsevier

June 3, 2016

ACCEPTED MANUSCRIPT

3D Modeling of Squeeze Flow of Multiaxial Laminates

1

Notre Dame University-Louaize, P.O. Box: 72, Zouk Mikael, Zouk Mosbeh, Lebanon [email protected]

ESI group chair, GeM UMR CNRS-Centrale Nantes, 1 rue de la Noe, F-44300 Nantes, France {Emmanuelle.Abisset-Chavanne;Christophe.Binetruy;Francisco.Chinesta}@ecnantes.fr University of Delaware, Newark, Delaware 19716-3119, United States [email protected]

ED

M

3

AN US

2

CR IP T

Ch. Ghnatios1,2 , E. Abisset-Chavanne2 , Ch. Binetruy2 F. Chinesta2,∗ , S. Advani3

Abstract

AC

CE

PT

Thermoplastic composites are attractive because they can be recycled and exhibit superior mechanical properties. As part of a more general research effort, in this paper attention is paid to squeeze flow of continuous fiber laminates. In the case of unidirectional prepregs, the ply constitutive equation, when elastic effects are neglected, can be modeled as a transversally isotropic fluid, that must satisfy the fiber inextensibility as well as the fluid incompressibility. When laminate is squeezed, depending on the plies orientation and the lubrication conditions at their interfaces, the flow kinematics exhibits a complex dependency along the laminate thickness requiring one to describe the detailed velocity field through the thickness. Simulations at this scale are computationally intensive or almost impossible with standard FE techniques. In this work the squeeze flow of laminates is revisited from the perspective of the Proper Generalized Decomposition that makes possible extremely detailed 3D calculations with a computational cost characteristic of 2D (sometimes 1D) simulations. Keywords: Squeeze flow, Composite laminates, Sheet Forming, Proper Generalized Decomposition, Ericksen fluid

2

ACCEPTED MANUSCRIPT 1

Introduction

AN US

CR IP T

Thermoplastic composites are preferred structural materials due to their excellent damage tolerance properties, shorter manufacturing cycles and ease of weldability. One of the precursor material to fabricate thermoplastic composite parts is an unidirectional (UD) prepreg which consists of aligned continuous fibers pre-impregnated with thermoplastic resin. In their melt state, UD prepreg can be viewed as inextensible fibers surrounded by an incompressible viscous matrix, and hence can be modeled as a transversally isotropic fluid [21]. These UD laminates are usually stacked in desired orientations to create a composite laminate. The composite laminates or sheets are consolidated to manufacture the net shape part using the forming process [17,16]. During the consolidation and forming process, the thermoplastic laminate is subjected to complex deformation which may include intraply and/or interply shear, ply reorientation and squeeze flow which could result in unacceptable defects such as wrinkles, sliding of prepregs or inefficient fiber distribution due to ply slippage. The squeeze flow behaviour of both unidirectional and multiaxial laminates has been studied in [15,5,14,19,20,13,6].

AC

CE

PT

ED

M

Understanding the influence of ply orientation and stacking on the squeeze flow behaviour of laminates is important as it will affect the material properties of the final part. In particular, experimental investigations have shown that cross-ply laminates produce a different squeeze flow response than unidirectional laminates [20]. An unstable flow response of the cross-ply laminates has been measured and attributed to the lack of material flow in the fiber direction. A stacking sequence effect was also noticed since laminates made by alternating orientation with each successive ply or by alternating lamina composed of many plies oriented in the same direction yielded different responses. The difference between behaviours of unidirectional and cross-ply laminates suggests a strong coupling between the ply surface area and the total volume of the ply. This insight calls for an in depth investigation to obtain a clear understanding of the flow kinematics constrained by the presence of continuous inextensible fibers. Since size effects are suspected from these experimental findings, only detailed enough computational approaches can clarify this. However, the solution of numerical models defined in laminates remains, even nowadays, a tricky issue because well experienced state of the art discretization techniques, like finite elements, fail to capture all the richness of solutions through the thin laminate thickness, and more particularly at the level of the plies composing the laminate. This difficulty is not new, it was encountered many years ago when addressing the solution of problems defined in plate and shell geometries. Plate and shell theories were developed in the structural mechanics framework for addressing quite simple behaviors. They allow, through the introduction 3

ACCEPTED MANUSCRIPT

CR IP T

of reasonable hypotheses, to reduce the 3D complexity to a 2D one related to the problem now formulated by restricting to in-plane coordinates. In the case of fluid flows this dimensionality reduction is known as lubrication theory and it allows efficient solutions of fluid flows in plate or shell geometries for many type of fluids, i.e. linear (Newtonian) and non-linear. In a former work [12] it was proven that lubrication theory cannot be applied in the case of unidirectional and multiaxial laminates involving different behavior across the thickness direction. Thus, 3D simulations are necessary for such cases to predict detailed kinematics to address formation of possible defects. This will increase the number of degrees of freedom dramatically thus requiring efficient discretization techniques to get tractable solutions.

CE

PT

ED

M

AN US

A new discretization technique based on the use of separated representations was proposed some years ago to address multidimensional models suffering the so-called curse of dimensionality, where standard mesh-based techniques fail [3]. The curse of dimensionality was circumvented thanks to those separated representations that transformed the solution of a multidimensional problem into a sequence of lower dimensional problems. The interested reader can refer to the recent reviews [8] and the references therein. A direct consequence of that approach was to separate the physical space. Thus in plate domains an inplane-out-of-plane decomposition was proposed for solving porous media flow models in laminates [8], then for solving elasticity problems [7] and coupled multiphysics problems [9]. In those cases the 3D solution was obtained from the solution of a sequence of 2D problems (the ones involving the in-plane coordinates) and 1D problems (the ones involving the coordinate related to the plate thickness). In [12] authors proposed a numerical model combining the solution of the 3D Stokes problem in cross-ply laminates with different viscosities. Then, a sequence of layers in which the flow was described by the Stokes equations was intercalated within another sequence described using Darcy’s equation. Both models were combined into a unified Brinkman equation. In order to ensure a detailed enough 3D representation despite the domain degeneracy, an in-plane-out-of-plane separated representation within the Proper Generalized Decomposition – PGD – framework was successfully proposed.

AC

As discussed above, thermoplastic composites are usually made of several prepreg plies with a particular orientation sequence to meet design requirements. In the case of unidirectional prepregs, the ply constitutive equation, when elastic effects are neglected, can be assimilated to a transversally isotropic fluid, that must satisfy the fiber inextensibility as well as the fluid incompressibility. When laminate is squeezed, depending on the plies orientation and the lubrication conditions at their interfaces, the flow kinematics exhibits a rich dependency along the laminate thickness requiring detailed 3D descriptions through the full thickness. Simulations at this scale are computationally intensive or almost impossible with standard FE techniques. In this work the squeeze flow of laminates is revisited from the perspective of the Proper Gener4

ACCEPTED MANUSCRIPT alized Decomposition that makes possible extremely detailed 3D calculations with a computational cost characteristic of 2D (sometimes 1D) simulations.

CR IP T

The next section addresses the solution of the Stokes flow problem in a narrow gap making use of the in-plane-out-of-plane separated representation within the Proper Generalized Decomposition framework. This approach will be then extended to a laminate consisting of different plies modeled with a viscous linear fluid of different viscosities. Section 3 extends the approach to the Ericksen fluid flow problem in the case of a monolayer and a laminate. Finally section 4 describes some numerical experiments, addressing: (i) the case of a sequence of Ericksen plies where the continuity of the velocity field is enforced at the ply interfaces; (ii) the case where a friction law at the Ericksen ply interfaces makes possible to introduce sliding of each ply; and (iii) the case where these plies are lubricated from a viscous inter-ply layer that accommodates different kinematic constraints applied to adjacent plies.

3D modeling of Stokes flow in narrow gaps

AC

2

CE

PT

ED

M

AN US

This work is not intended to be exhaustive and cover all scenarios of practical interest. It only aims to introduce a simulation framework able to solve mechanical models defined in laminates exhibiting a very rich behavior in both the in-plane and the through-thickness directions. Most of its applicative consequences will be explored in future works. The main aim of the present work is to prove that high fidelity solutions of complex fluid and flows can be efficiently addressed within the context of separated representations, making possible an extremely fine representation of the flow kinematics occurring during the squeeze flow of laminates, reaching a level of detail rarely attainable by using more experienced state of the art discretization techniques. At present our analysis only concerns the flow kinematics (velocity field). The analysis of other model outputs (pressure field or fibers tension) are under progress, mainly due to the fact that the Ericksen model considered in the present work has difficulties to accommodate flow shear and fibers tension at the plies interfaces. These developments will be deeply analyzed in ongoing publications.

2.1 In-plane-out-of-plane separated representation The in-plane-out-of-plane separated representation allows for solution of full 3D models defined in plate geometries with a computational complexity characteristic of 2D simulations. This separated representation allows for independent representations of the in-plane and the thickness fields dependencies. The main idea lies in the separated representation of the velocity field by using 5

ACCEPTED MANUSCRIPT functions depending on the in-plane coordinates x = (x, y), Pji (x, y), and others depending on the thickness direction z, Tji (z), j = 1, 2, 3 and i = 1, · · · N according to:       

v(x, z) =



u(x, z)      

v(x, z)

w(x, z)





N P

Pi1 (x)



Ti1 (z) 

·    i=1  P   N 2 2 ,  P (x) · T (z) i   i=1 i   N  P 3 3 i=1

(1)

Pi (x) · Ti (z)

Eq. (1) can be rewritten in the compact form N X

Pi (x) ◦ Ti (z),

AN US

v(x, z) ≈

CR IP T

which leads to a separated representation of the strain rate, when introduced into the flow problem weak form allows the calculation of functions Pi (x, y) by solving the corresponding 2D equations and functions Ti (z) by solving the associated 1D equations, as described later.

i=1

(2)

where ”◦” denotes the entry-wise or Hadamard’s product.

M

Remark 1. In the previous separated representation the in-plane coordinates were separated from the out-of-plane one. When considering plate-like domains with rectangular in-plane shape a more efficient separation could be performed by assuming the fully separated representation N X

ED v(x, y, z) ≈

i=1

Xi (x) ◦ Yi (y) ◦ Zi (z),

(3)

AC

CE

PT

that allows computing the 3D velocity field from a sequence of one-dimensional problems. Even if the numerical simulations discussed in the this paper consider this kind of separable domains, for the sake of generality, we do not considered such a fully separated representation, because problems of industrial interest are in general defined in laminates with complex in-plane geometries, making difficult the separation of of coordinates x and y, i.e. for example a circle expressed from a sum of functions of x and functions of y requires an infinite number of terms.

Remark 2. If a and b are vectors of the same dimension, vector c, defined from c = a ◦ b, has as components ci = ai bi , having the same dimension as a and b. If a and b are second order tensors with the same size, tensor c, defined from c = a ◦ b, has components cij = aij bij (no sum with respect to the repeated 6

ACCEPTED MANUSCRIPT indexes). In this case, a : b = c, with the scalar c given by c = aij bij now considering sum with respect to the repeated indexes (Einstein’s summation convention). Using notation in (1), the velocity gradient ∇v(x, z) can be written as: ∂u ∂u ∂u ∂x ∂y ∂z



     ∂v ∂v ∂v     ∂x ∂y ∂z    ∂w ∂w ∂w ∂x ∂y ∂z





 N  X    i=1 

N X i=1

∂Pi1 ∂Pi1 ∂x ∂y ∂Pi2 ∂Pi2 ∂x ∂y ∂Pi3 ∂Pi3 ∂x ∂y

Pi1





1  Ti

     2 ◦ Pi2    Ti  

Pi3

Pi (x) ◦ Ti (z).

Ti3

∂T 1 Ti1 ∂zi ∂T 2 Ti2 ∂zi ∂T 3 Ti3 ∂zi

      

=

CR IP T

∇v =



(4)

2.2

AN US

The solution of full 3D Stokes problem within the in-plane-out-of-plane separated representation is revisited in the next sections. Flow model

M

The Stokes flow model is defined in Ξ = Ω × I, Ω ⊂ R2 and I ⊂ R, and for an incompressible fluid, in absence of inertia and mass terms reduces to:     ∇·σ   

=0

= −pI + 2ηD ,

ED σ      

(5)

∇·v =0

D=

∇v + (∇v)T . 2

(6)

AC

CE

PT

where σ is the Cauchy’s stress tensor, I the unit tensor, η the fluid viscosity, p the pressure (Lagrange multiplier associated with the incompressibility constraint) and the rate of strain tensor D defined as

A penalty formulation is used to circumvent the issue related to stable mixed formulations (that require fulfilling the so-called LBB conditions) within the separated representation, which to this day is an open issue. The mass balance is modified by introducing a penalty coefficient λ whose value is usually very small. This approach was successfully considered in our previous works [12,1,2]. The penalty formulation can be written as: ∇ · v + λ p = 0, 7

(7)

ACCEPTED MANUSCRIPT or more explicitly Tr(D) ∇·v =− , (8) λ λ where Tr() refers the trace operator. The trace of a tensor, and in particular the trace of the rate of strain, can be written as Tr(D) = D : I. p=−

The weak form of the penalized Stokes problem, for a test velocity v∗ vanishing on the boundary in which the velocity is prescribed, and assuming null tractions in the remaining part of the domain boundary, can be written as

Ω×I



1 Tr(D∗ )Tr(D) + 2ηD∗ : D λ



dx dz = 0,

where as proved in the appendix

N X 4   ηX A∗jk (x) : Bjk (z) + Ajk (x) : B∗jk (z) , 2 j=1 k=1

AN US

2ηD∗ : D ≈

(9)

CR IP T

Z

and

N   1 1X Tr(D∗ ) · Tr(D) ≈ F∗j (x) : Gj (z) + Fj (x) : G∗j (z) , λ λ j=1

(10)

(11)

M

with matrices introduced in Eqs. (10) and (11) are defined in the appendix.

Separated representation constructor

ED

2.3

CE

PT

The construction of the solution separated representation is performed incrementally, a term of the sum at each iteration. Thus, supposing that at iteration n − 1, n ≥ 1, the first n − 1 terms of the velocity separated representation were already computed v

n−1

(x, z) =

n−1 X i=1

Pi (x) ◦ Ti (z),

(12)

AC

the terms involved in the weak form (9) are: ∗

n−1

D :D

4   XX 1 n−1 A∗jk (x) : Bjk (z) + Ajk (x) : B∗jk (z) , = 4 j=1 k=1

and Tr(D∗ ) · Tr(Dn−1 ) =

n−1 X



F∗j (x) : Gj (z) + Fj (x) : G∗j (z) ,

j=1

respectively. 8

(13)

(14)

ACCEPTED MANUSCRIPT When looking for the improved velocity field vn (x, z) at iteration n vn (x, z) =

n X i=1

Pi (x) ◦ Ti (z) = vn−1 (x, z) + Pn (x) ◦ Tn (z),

(15)

we consider the test function v∗ (x, z) v∗ = P∗ ◦ Tn + Pn ◦ T∗ ,

(16)

∇v∗ = P∗ ◦ Tn + Pn ◦ T∗ .

(17)

In that case it results n−1 4 XX

j=1 k=1

4 X



(A∗nk (x) : Bnk (z) + Ank (x) : B∗nk (z)) ,

k=1

and



A∗jk (x) : Bjk (z) + Ajk (x) : B∗jk (z) +

Tr(D∗ ) · Tr(Dn ) = F∗n (x)

n−1 X j=1

AN US

4D∗ : Dn =

CR IP T

that implies

(18)



F∗j (x) : Gj (z) + Fj (x) : G∗j (z) +

: Gn (z) + Fn (x) : G∗n (z)

(19)

M

where all the matrices involved in Eqs. (18) and (19) are given in the appendix, with the test function ∇v∗ given by Eq. (17). Z

1 ∗ (F (x) : Gn (z) + Fn (x) : G∗n (z)) dx dz+ λ n

PT

Ω×I

ED

Thus the problem weak form (9) writes at iteration n:

Z

4 η X A∗ (x) : Bnk (z) + Ank (x) : B∗nk (z) 2 k=1 nk

CE

Ω×I

AC



Z

Ω×I

Z

Ω×I





!

dx dz =



 X 1 n−1 F∗j (x) : Gj (z) + Fj (x) : G∗j (z)  dx dz− λ j=1 

4   XX η n−1 A∗jk (x) : Bjk (z) + Ajk (x) : B∗jk (z)  dx dz. 2 j=1 k=1

(20)

The extended weak form (20) becomes nonlinear because it involves the product of unknown functions Pn and Tn . Thus a linearization strategy becomes necessary, the simplest one being an alternating fixed point algorithm that proceeds as follows: 9

ACCEPTED MANUSCRIPT

n

v (x, z) =

n X i=1

CR IP T

(1) Assuming functions Pn (x) are known (arbitrarily chosen at the first iteration of the nonlinear iteration) matrices A∗jk and F∗j vanish. Being all functions depending on x known, integrals in Ω in (20) can be calculated. Thus, it finally results in a one dimensional linear problem that involves the three scalar functions involved in Tn (z), Tn1 (z), Tn2 (z) and Tn3 (z). (2) Then, with the just computed function Tn (z), and with B∗jk and G∗j vanishing, one can proceed to integrate Eq. (20) in I. It finally results in a two-dimensional linear problem for the unknown function Pn (x) that involves the three scalar functions Pn1 (x), Pn2 (x) and Pn3 (x). (3) The convergence is checked by comparing functions Pn and Tn between two consecutive iterations of the nonlinear solver. If both functions are small enough they are used to update the velocity field Pi (x) ◦ Ti (z).

(21)

AN US

If the convergence is not attained, one returns to step 1 with the calculated functions Pn to re-compute Tn

Flow in a laminate

ED

2.4

M

Because of the one-dimensional large scale variation present in the laminate thickness direction one can employ extremely detailed descriptions along the thickness direction without sacrificing the computational efficiency of the 3D solution procedure.

PT

Consider a laminate composed of P layers in which each layer involves a linear and isotropic viscous fluid of viscosity ηi , thus the extended Stokes flow problem in its weak form involves the dependence of the viscosity along the thickness direction.

AC

CE

If H is the total laminate thickness, and assuming for the sake of simplicity and without loss of generality that all the plies have the same thickness h, . Now, from the characteristic function of each ply χi (z), it results h = H P i = 1, · · · , P:  the viscosity reads

χi (z) =

 1  0

if (i − 1)h ≤ z < ih

,

(22)

elsewehere

η(x, z) =

P X i=1

ηi · χi (z),

(23)

where it is assumed, again without loss of generality, that the viscosity does not evolve in the plane, i.e. ηi (x) = ηi . 10

ACCEPTED MANUSCRIPT This decomposition is fully compatible with the velocity separated representation (1) and with the in-plane-out-of-plane decomposition considered for solving Eq. (20).

3

Ericksen fluid flow model in a laminate

AN US

CR IP T

The case of a prepreg ply reinforced by continuous fibres oriented along direction pT = (px , py , 0), kpk = 1, is analyzed here. It is assumed that the thermoplastic resin exhibits Newtonian behaviour. Thus the velocity v(x, z) of the equivalent anisotropic fluid must satisfy the incompressibility and inextensibility constraints ∇ · v = 0, (24) and ∇v : (p ⊗ p) = 0, (25) respectively. Expression (25) can be rewritten using tensor notation as ∇v : a = 0, where the second order orientation tensor a is defined from a = p ⊗ p.

The orientation tensor a has only planar components (the out-of-plane fiber orientation can be neglected in the case of laminates), it is symmetric and of unit trace, i.e.   axy 0    ayy 0   

M

a=

 axx    ayx  

0

0

0





A 0 

=

0T 0

,

(26)

ED

where A represents the plane component of the orientation tensor a, axy = ayx (i.e. A = AT ) and ayy = 1 − axx .

PT

The simplest expression of the Ericksen’s constitutive equation [11] can be written in the compact form as follows

CE

σ = −pI + T a + 2ηT D + 2(ηL − ηT )(D · a + a · D),

(27)

AC

and introduced in the linear momentum balance ∇ · σ = 0,

(28)

In Eq. (27) p and T represents respectively the Lagrange multipliers related to the incompressibility and inextensibility, and ηL and ηT the longitudinal and transverse shear viscosities respectively. By introducing the incompressibility and inextensibility constraints from a penalty formulation it results: ∇ · v + λp = 0, 11

(29)

ACCEPTED MANUSCRIPT and ∇v : a + ξT = 0, with λ and ξ small enough, leading to: p=−

(30)

1 ∇·v = − Tr(D), λ λ

(31)

and

∇v : a 1 = − D : a, ξ ξ where the fact that a is symmetric ∇v : a = D : a is used. T =−

CR IP T

(32)

AN US

Equations (31) and (32) allow removing pressure and tension fields involved in the flow problem as described below, however both penalty coefficients λ and ξ must be chosen small enough to ensure that the computed velocity remains accurate enough. However, the inverse procedure for obtaining from both expressions the pressure and the tension as soon as the velocity field will be available is not a valuable option because the computed fields remain in general too inaccurate (the results depend too much on the choice of the penalty coefficients). Thus, for calculating accurately pressure and tension an adequate mixed formulation seems the best option.

Z

M

The weak form for a test velocity v∗ (x, z) vanishing at the boundary in which velocity is prescribed and assuming null tractions in the remaining part of the domain boundary can be expressed as D∗ : σ dx dz = 0,

(33)

ED

Ω×I

Z

D∗ : σ dx dz =

Ω×I

D:a Tr(D) I− a + ηT D + η˜ (D · a + a · D) λ ξ

CE

Z

PT

By introducing the Ericksen constitutive equation (27) as well as both penalty expressions (29) and (30), it can be written as



D :

Ω×I

!

dx dz = 0,

(34)

AC

with η˜ = ηL − ηT . This integral form can be rewritten as Z

Ω×I

(

Tr(D∗ ) · Tr(D) (D∗ : a) · (D : a) − + λ ξ

ηT D∗ : D + η˜D∗ : (D · a + a · D)} dx dz = 0.

(35)

At this stage the in-plane-out-of-plane separated representation constructor of v(x, z) proceeds as described in the previous section. 12

ACCEPTED MANUSCRIPT

Remark 3. If a = 0 this formulation reduced to one related to the Stokes flow problem.

a(x, z) =

P X

ai (x)χi (z).

i=1

CR IP T

Remark 4. Laminates can be addressed by associating with each ply the planar fiber orientation pi (x), with its out-of-plane component vanishing, from which the associated orientation tensor ai (x) results in ai (x) = pi (x) ⊗ pi (x). Using again the characteristic function of the i-ply, χi (z), i = 1, · · · , P, the orientation tensor in the laminate, a(x, z), can be expressed as (36)

a(z) =

AN US

Remark 5. If the fiber orientation is constant in each plane, then the laminate orientation tensor can be expressed as P X

ai χi (z).

(37)

i=1

Numerical results

M

4

PT

ED

The numerical results discussed hereafter consider laminates composed of different number of plies, different thicknesses and different behavior laws. The domain occupied by the laminate has a length L, width W and a thickness H, i.e. Ξ = Ω × I, with Ω = [0, L] × [0, W ] and I = [0, H], with x ∈ Ω and z ∈ I.

AC

CE

Prescribed velocities are enforced at the top and bottom boundaries whereas null tractions apply on the lateral boundaries, i.e.:    v(x, y, z = 0) = 0       T   v(x, y, z = H) = (0, 0, −V )       σ(x = 0, y, z) · nx = 0    σ(x = L, y, z) · nx = 0        σ(x, y = 0, z) · ny = 0       σ(x, y = W, z) · n = 0 y

where nx = (1, 0, 0)T and ny = (0, 1, 0)T . 13

,

(38)

CR IP T

ACCEPTED MANUSCRIPT

AN US

Figure 1. Squeezing with prescribed velocity two laminates composed of 5 plies: (left) a [90/0/90/0/90] laminate and (right) [0/90/0] with two intercalated lubricating Stokes layers.

M

It is important to note that only the instantaneous velocity field is calculated. The real process simulation should imply solving the flow problem, then updating the domain and the fiber orientation before recomputing the velocity field and so on. As the magnitude of the applied compression remains very small the main features related to the flow kinematics can be extracted from the first problem solution as discussed later.

AC

CE

PT

ED

It is assumed that during compression molding, the upper mold moves down with a prescribed velocity V during the consolidation. Thus, mass conservation leads to significant velocity variations within Ω, e.g. the central point has a null in-plane velocity because of the symmetry condition whereas the in-plane velocity is maximal at the laminate lateral boundaries ∂Ω×I. Moreover, when the rheology varies significantly along the laminate thickness (e.g. as viscosity depends on temperature, the last evolving throughout the thickness) or to take into account the through-thickness complex kinematics when moving from one layer to the next with different constitutive equation and/or reinforcement orientations, a sufficiently detailed solution is also required along the thickness direction to capture all the characteristics of the solution. Figure 1 depicts the laminate geometry as well as the squeezing conditions.

4.1 Laminate composed of a sequence of Ericksen plies 4.1.1

Cross-ply laminate

The first test case consists of a laminate composed of 5 plies described by the Ericksen constitutive equation, with the unidirectional continuous fiber rein14

ACCEPTED MANUSCRIPT forcement defining the orientation sequence 0/90/0/90/0. The reinforcement orientation is defined by pTi = (1, 0, 0), for i = 1, 3, 5 and pTj = (0, 1, 0), for j = 2, 4, implying the orientation tensors 







1 0 0 



 0 0 0, ai =      000

aj = for j = 2, 4.

0   0  

0 0  

, 1 0   000

CR IP T

for i = 1, 3, 5; and

(39)

(40)

AN US

The following values are used : L = W = 0.5 m, H = 1.04 mm, (0.208 mm each ply), ηL = 200 Pa · s, ηT = 20 Pa · s and both penalty coefficients are kept very small, λ = ξ = 10−8 . The discrete model was defined by using 3600 nodes uniformly distributed in Ω and 500 along I that corresponds with 100 nodes in each ply. The squeezing rate was V = 0.1 mm · s−1 .

AC

CE

PT

ED

M

The separated representation constructor computes progressively the terms of the finite sum. After each calculation (after calculating a term of the finite sum, also called ”mode”) the solution accuracy is checked. The two most usual procedures for checking the solution accuracy are based, the first on the norm of the just calculated mode and the second on the norm of the equation residual. In the first, the norm of the current mode is compared with the norm of the first mode that is usually the one of highest norm. When the relative norm with respect to the first mode norm becomes small enough, it indicates that the current mode contribute very little to the solution improvement and then the iteration process can be stopped. In fact, such a criterion does not guarantee the convergence towards the exact solution. Thus, a second procedure to better quantify the solution accuracy consists of introducing the current solution into the partial differential equation and then computing the L2 norm of its residual. Both strategies are described in detail in [10] and other alternatives based on the calculation of quantities of interest were addressed in [4]. In general terms, the number of required modes depends on the solution separability and its ”a priori” estimation is not an easy task. With respect to the choice of the in-plane and out of plane meshes the procedure is quite similar to the one considered in finite element implementations. After solving a problem, the same problem is solved again by using a refined 15

ACCEPTED MANUSCRIPT mesh and the convergence is checked. In our simulations, we considered from the beginning an extremely fine mesh in the thickness directions because as it involves the solution of one-dimensional problems the associated computational complexity is negligible with respect to the one related to the in-plane problems solution. More sophisticated choices for checking the convergence and adapting the meshes within the PGD framework can be found in [4,18].

AN US

CR IP T

The separated representation of the velocity field involves 6 modes (N = 6) that requires the solution of about 20 two-dimensional problems of size 3×3600 (the number of nodes multiplied by the number of components of the velocity field – 3 in 3D –) and the same number of one-dimensional problems of size 3×500. The number of problem solutions (about 20) is higher than the number of modes (N = 6) because to compute each term of the finite sum one must solve a nonlinear problem that requires a few number of iterations. Thus the solution of about 20 problems of size 3×3600 allows computing the 3D problem solution that would require if using a 3D mesh the solution of a linear problem of size 1.8 106 , and such a discretization would involved very distorted elements (with a ratio of the in-plane to the thickness characteristic lengths of about 3000).

CE

PT

ED

M

In order to check quantitatively the solution obtained using the PGD separated representation, a 3D finite element solution was performed in an equivalent discretization making use of Q8 finite elements (the standard 8-node velocity based brick element with trilinear shape functions) consisting of 60 × 60 × 500 nodes. When considering the finite element solution as reference, the difference with respect to the PGD solution was computed and the error using a L2 norm obtained (the L2 norm is preferred because we are interested in attaining an acceptable accuracy along the domain thickness and the maximum norm is mostly controlled by the field at the plies interfaces). The L2 norm was lower than 0.1% proving the efficiency of the proposed approach (the maximum norm was lower than 1%, decreasing with the number of modes considered in the velocity separated representation). Obviously, when considering richer configurations the comparison was not performed because it was impossible solving the resulting the 3D finite element problem.

AC

Figure 2 depicts the velocity magnitude. It can be noticed that the velocity components along the reinforcement directions vanish, and only velocities transverse to the fiber direction exist. Moreover, when assuming the continuity of the velocity across the plies, the in-plane velocity becomes zero because any non-zero velocity (different to a rigid motion difficult to imagine by symmetry considerations) is incompatible with the inextensibility of fibers oriented along the x and the y-directions in the plies sharing each interface. As a consequence the shear stress is maximum at the plies interfaces, which has important consequences as discussed later. 16

CR IP T

ACCEPTED MANUSCRIPT

AN US

Figure 2. Velocity magnitude (m · s−1 ) of a 0/90/0/90/0 laminate with plies modeled by the Ericksen constitutive model. As the in-plane components of the velocity vanish at both laminate surfaces, and the through-thickness component of the velocity at the top surface is two orders of magnitude smaller than the maximum in-plane velocity, the magnitude of the velocity at the top surface seems being zero (according to the color) but it is in fact the one related to the compression velocity, i.e. V = 10−4 m · s−1 . The same consideration applies for most figures in the present paper.

AC

CE

PT

ED

M

In what follows we illustrate the progressive construction of the separated representation solution. As described below, at the first iteration we compute the first functional products given an approximation of the solution, but as a single product is not enough for approximating the exact flow problem solution, other terms are then calculated to enrich the solution ensuring progressively the convergence towards the problem solution. Figures 3-8 depict the most significant modes involved in the in-plane-out-of-plane separated representation of the velocity field. Using the notation introduced in Eq. (1): P11 (x) and T11 (z) in Fig. 3; P12 (x) and T12 (z) in Fig. 4; P13 (x) and T13 (z) in Fig. 5; P21 (x) and T21 (z) in Fig. 6; P22 (x) and T22 (z) in Fig. 7; and finally P23 (x) and T23 (z) in Fig. 8. The first functional product P11 (x) · T11 (z) depicted in Fig. 3 is, as expected, higher at plies oriented along 90◦ , being minimum at the three plies oriented along 0◦ direction, due to the fact that the x-component of the velocity vanishes in plies oriented at 0◦ . However, a single product does not suffice to ensure a null x-component of the velocity at these plies. Other modes are then required for that purpose. We can notice in Fig. 6 that the correction associated with the second functional product P21 (x) · T21 (z) has opposite sign in plies oriented along 0◦ and it tends to enforce a null velocity along the x-direction in plies oriented along 0◦ direction. A similar reasoning applies for the y-component of the velocity. Moreover, Fig. 5 proves that the z-component of the velocity related to the first mode is almost linear, being maximum at the top and null at the bottom surface, in agreement with the expected behavior. 17

ACCEPTED MANUSCRIPT

0.04

0.03

1.2

0.02

1

×10-3

0.01

0.8

0

0

-0.01

-0.02

-0.02

-0.04 0.6

-0.03

0.4 0.2 Y(m)

0

0

0.2

z (m)

P 11( x)

0.02

0.4 0.2

0.6

0.4

0.6

0 -2

-1.5

-1 T1(z)

X(m)

1

-0.5

0 ×10-3

CR IP T

Figure 3. Modes involved in the separated representation: P11 (x) (left) and T11 (z) (right).

When considering all the modes involved in the separated representation and particularizing for xM = (0.48, 0.48), Fig. 9 depicts the solution u(xM , z) and v(xM , z) whose expressions read 

 u(xM , z)   

v(xM , z)



  N P 1 1 P (x ) · T (z) M i   i=1 i .    P N 2 2

AN US



i=1

(41)

Pi (xM · Ti (z)

It can be noticed that the reconstructed solution when adding all the modes verifies the kinematic constraints previously discussed.

CE

PT

ED

M

To further emphasize the capabilities of the proposed computational approach, the consolidation of a cross-ply laminate composed of 70 plies ([0/90]35 ) is analysed. We consider the geometry that was addressed in [20] consisting of L = W = 1 m, H = 7.28 mm and use the same rheology that in the previous test case. The discrete model consists of 3600 nodes uniformly distributed in Ω and 1400 nodes uniformly distributed in I, that corresponds to 20 nodes in each ply. The in-plane-out-of-plane separated representation in this case involves 9 modes. Figures 10 and 11 depict respectively the velocity magnitude and a detail of it to better understand the velocity evolution throughout the laminate thickness.

Multiaxial laminate

AC

4.1.2

To further emphasize the capabilities of the proposed computational approach, in what follows we analyse the consolidation of a laminate consisting of 37 nonorthogonal plies, where each ply is inclined at an angle of 10◦ to each other, to give the [0/10/20/ · · · /350/360] stacking sequence (from top to the bottom). The geometrical dimensions of the laminate are L = W = 1 m, H = 3.64 mm and the same rheology as in the previous test case is used. The discrete model consists of 3600 nodes uniformly distributed in Ω and 1400 nodes uniformly 18

ACCEPTED MANUSCRIPT ×10 -3

×10-3 4

2

0.8 z (m)

0

0

-1

-2

×10-3

1

1

2 P 21( x)

1.2

0.6 0.4

-2

-4 0.6 0.4 0.2 0

Y(m)

0.2

0

0.2

0.6

0.4

0 -0.03

-0.025

-0.02

X(m)

-0.015 T2(z)

-0.01

-0.005

0

1

×10 -6 2

×10-6 3

1.5

0.8

0.5

1

z(m)

0

0

-0.5

0.4

-1

-2 0.6

-1.5

0.4 0.2 0

Y(m)

0.4

0.2

0

0.6

AN US

-1

×10-3

1

1

2 P 31( x)

1.2

CR IP T

Figure 4. Modes involved in the separated representation: P12 (x) (left) and T12 (z) (right).

0.2

0.6

0 -0.5

X(m)

-0.4

-0.3

3

-0.2

-0.1

0

T1(z)

M

Figure 5. Modes involved in the separated representation: P13 (x) (left) and T13 (z) (right). 0.08 0.06 0.04

P 12( x)

0.05 0

-0.1 0.6 0.4

PT

-0.05

0.2

0

0

-0.02 -0.04

0.4

1

0.6 0.4

-0.06 -0.08

0.2

×10-3

0.8

0

0.6

0.2 0 -0.02

X(m)

CE

Y(m)

0.02

z (m)

ED

0.1

1.2

-0.015

-0.01

-0.005 T12(z)

0

0.005

0.01

AC

Figure 6. Modes involved in the separated representation: P21 (x) (left) and T21 (z) (right).

distributed in I. The in-plane-out-of-plane separated representation involves 9 modes. Figures 12, 13 and 14 depict respectively the velocity magnitude, a close-up view of its evolution through the thickness direction and the same close-up view but magnifying the velocity magnitude to emphasize its 3D fine resolution. That goal was successfully achieved thanks to the use of the in-planeout-of-plane separated representation. 19

ACCEPTED MANUSCRIPT 0.15

1.2

×10-3

0.1

1

0.2

0.05

0.8

0

0

z (m)

P 22( x)

0.1

-0.05

-0.1

0.4

-0.1

-0.2 0.6

0.2

-0.15

0.4 0.2 0

Y(m)

0

0.2

0.6

0.4

0.6

0 -6

-4

-2

X(m)

T2 (z)

0

2

2

4 ×10-3

×10 -6

×10-3

2.5

-6

×10 3

1

2

1

1 0.5

0

0

-1

-0.5

-2 0.6

-1

0.4 0.2 0

0

0.2

0.4

X(m)

0.6 0.4

0.6

0.2

0 -3

-2

-1

3

0

1

2

T2 (z)

M

Y(m)

0.8

z (m)

2

1.5

AN US

P 32( x)

1.2

CR IP T

Figure 7. Modes involved in the separated representation: P22 (x) (left) and T22 (z) (right).

ED

Figure 8. Modes involved in the separated representation: P23 (x) (left) and T23 (z) (right).

AC

CE

PT

Fig. 15 depicts the through-thickness component of the velocity on the midplane of the half of plies located at the bottom of the laminate, with fibers oriented along directions 360◦ , 330◦ , 300◦ , 270◦ , 240◦ and 210◦ . It can be observed the complex kinematics successfully captured by the in-plane-outof-plane decomposition. As expected the lowest velocity occurs in the bottom plane (360◦ ). We also verified that at the mid-plane of the top ply, with fibers along 0◦ , the z-component of velocity is very close to the compression rate of 10−4 m · s−1 . In order to check that the velocity field ensures the fibers inextensibility in each ply, the projection of the velocity in the fiber direction, i.e. v · p, is presented in Fig. 16. As expected the scalar product remains almost constant along the fiber direction in each ply, in agreement with the fiber inextensibility. The directional gradient was found to be of the order 10−6 for velocities of the order of 10−2 m · s−1 . 20

ACCEPTED MANUSCRIPT

1.2

×10 -3

1

0.6

0.4

0.2

0

0

0.005

0.01

1.2

×10 -3

1

M

0.8

0.6

ED

z (m)

0.015

AN US

u (m · s-1 )

CR IP T

z (m)

0.8

PT

0.4

0.2

0

CE

0

0.002

0.004

0.006

0.008

0.01

0.012

0.014

0.016

0.018

v (m · s-1 )

AC

Figure 9. Velocity (m · s−1 ) evolution through the thickness: u(xM , z) (top) and v(xM , z) (bottom); xTM = (0.48, 0.48) for the 0/90/0/90/0 laminate

. 4.2

Laminate composed of a sequence of Ericksen plies with friction at the plies interfaces

In this case the laminate is composed by two 0/90 Ericksen plies and viscous friction (the shear stress at the commun inter-plies interface is proportional to 21

CR IP T

ACCEPTED MANUSCRIPT

ED

M

AN US

Figure 10. Velocity magnitude (m · s−1 ) of a [0/90]35 laminate with plies modeled by the Ericksen constitutive model. The norm of the velocity at the top surface coincides with the compression velocity V = 10−4 m · s−1 .

PT

Figure 11. Detail of the solution depicted in Fig. 10 in the proximity of the laminate corner.

AC

CE

the relative velocity between them) is allowed at the interface between them, where the friction coefficient is supposed low, then moderate and finally high. It was noticed that increasing the friction coefficient makes the in-plane kinematics more and more constrained with the expected decrease of the maximum in-plane velocities. 4.3 Laminate composed of Ericksen plies with intercalated lubricating Stokes layers In this case the objective is to investigate the effect of intercaling between two consecutive plies modeled by the Ericksen constitutive model, a polymer layer accurately described by the Stokes flow model. 22

CR IP T

ACCEPTED MANUSCRIPT

PT

ED

M

AN US

Figure 12. Velocity magnitude (m · s−1 ) of a [0/10/20/ · · · /350/360] laminate with plies modeled by the Ericksen constitutive model). Again, the velocity magnitude at the top surface coincides with the compression velocity.

CE

Figure 13. Close-up of the solution depicted in Fig. 12.

AC

The case of a laminate composed of 0/90/0 Ericksen plies with two intercalated Stokes lubrication layers of viscosity η = 2 Pa · s is considered. When comparing the lubricated laminate with the one in which the lubricating layers are absent it was noticed that the in-plane velocities in the Ericksen ply located in the middle of the laminate are significantly higher than the ones in absence of lubrication. This behavior can be explained by the fact that the in-plane kinematics in the Stokes layer is not constrained and then it transmits the shear to the layers in contact, which in turn enhanced the in-plane velocity transverse to the fiber orientation. 23

CR IP T

ACCEPTED MANUSCRIPT

5

AN US

Figure 14. Close-up of the solution depicted in Fig. 12 with the velocity magnitude magnified in order to better appreciate the fast evolution of the velocity magnitude, justifying the necessity of high resolution approximations as the one provided by the separated representation here employed.

Conclusions

PT

ED

M

A new numerical procedure is proposed to simulate the squeeze flow of multiaxial laminates, able to present resolution levels never envisaged until now. It is based on the use of an original in-plane-out-of-plane separated representation of the velocity field, that allows calculating extremely detailed 3D solutions while keeping the computational complexity characteristic of 2D problems. Thus, extremely high resolutions can be attained in the thickness direction and able to capture localized behaviors, as the ones occurring at the ply interfaces when adjacent plies have different fiber orientations, when viscous lubrication occurs or when the interface sliding with viscous friction occurs.

AC

CE

Numerical simulations allowed us to identify very localized behaviors, and to understand the role of lubrication and friction, thanks to an extremely detailed resolution in the thickness direction. In this paper we only discussed kinematics because flow incompressibility and fiber inextensibility were considered within a penalty formulation. In our recent researches both constraints were considered by using appropriate Lagrange multipliers and the associated mixed weak form. Thus, we have access to the fluid pressure as well as to the fibers tension, leading to more valuable physical outputs (e.g. rheological behavior of laminates or fibers behavior that as expected exhibits different behavior depending on the solicitation mode extension or compression -). These results constitute an ongoing publication. 24

ACCEPTED MANUSCRIPT

×10 -7

1

×10 -5 -3.8

1

-2

0.9

0.9

-3.9

-2.5

0.8

0.8 -4

-3

0.7

0.7 -3.5

-4.1

0.6 -4

0.5

y(m)

y(m)

0.6

0.5

-4.2

-4.5

0.4

0.4 -4.3

-5

0.3

0.3 -5.5

0.2

-4.4

0.2 -6

0.1

0.1

-4.5

-6.5

0

0.1

0.2

0.3

0.4

0.5 X(m)

0.6

0.7

0.8

0.9

0

1

×10 -5

1

0

0.1

0.2

0.3

0.4

0.5 X(m)

0.7

0.8

0.9

1

×10 -5

1

-2.8

0.6

CR IP T

0

-1.9

0.9

0.9 -2.9

-2

0.8

0.8 -3

-2.1

0.7

0.7 -3.1

-2.2

0.6 -3.2

0.5

y(m)

y(m)

0.6

-2.3

0.5

-2.4

-3.3

0.4

0.4

-2.5

0.3

-3.4

0.3

0.2

-3.5

0.2

-2.6

-2.7

0

0.1

0.2

0.3

0.4

0.5 X(m)

0.6

0.7

0.8

0.9

×10 -5 -1

0.9

-1.1

0.8

-1.2

0.7

-1.3

0.6

-1.4

0.5 -1.5

0

0.1

0.2

0.3

0.4

0.5 X(m)

0.6

0.7

-2.8

0.8

0.9

1

×10 -6

1

0.9

-3

0.8

-4

0.7

-5

0.6

-6

0.5

M

y(m)

0

1

1

0.1

AN US

0

-3.6

y(m)

0.1

0.4

0.4

-7

-1.6

0.3

0.3 -8

-1.7

0.2

0.2

-1.8

0.1

-9

0.1

0

0

0.1

0.2

0.3

0.4

0.5 X(m)

ED

-1.9

0.6

0.7

0.8

0.9

0

1

-10

0

0.1

0.2

0.3

0.4

0.5 X(m)

0.6

0.7

0.8

0.9

1

Separated representation of the Stokes weak form

AC

A

CE

PT

Figure 15. Through-thickness component of the velocity (m · s−1 ) in the mid-plane of plies with fibers oriented along directions: 360◦ (top-left); 330◦ (top-right); 300◦ (middle-left); 270◦ (middle-right); 240◦ (bottom-left) and 210◦ (bottom-right).

The efficient computer implementation of the separated representation constructor discussed in Section 2.3 needs a separated form of the flow problem expressed in a weak form (9). For that purpose we first consider the second term D∗ : D, that takes into account expression (6) as follows 4D∗ : D = ∇v∗ : ∇v + ∇v∗ : (∇v)T + (∇v∗ )T : ∇v + (∇v∗ )T : (∇v)T . (A.1) 25

ACCEPTED MANUSCRIPT

×10 -5

1 0.9

×10 -4

1 0.9

6

0.8

1

0.8 4

0.7

0.7 2

0.5

0

0.4

0.6 y(m)

y(m)

0.6

0.5

0

0.4

-2

0.3

0.3 -4

0.2

0.2 -6

0

0

0.1

0.2

0.3

0.4

0.5 X(m)

0.6

0.7

0.8

0.9

0

1

×10 -4

1 0.9

-1

0.1

0

0.1

0.2

0.3

0.4

0.5 X(m)

1.5

0.7 0.6 y(m)

y(m)

0

0.4

0.5

0

0.4

-0.5

-0.5

0.3

0.3 -1

0.2

0.1

0.2

0.3

0.4

0.5 X(m)

0.6

0.7

0.8

0.9

0

1

×10 -4

1

2

0.9 0.8

1

0.6

0.1

0.2

0.3

0.4

0.5 X(m)

0.6

0.7

-1.5

0.8

0.9

1

×10 -4

1

4

0.9

3

0.5

0

0.7

2

0.6

1

0.5

0

M

y(m)

0

0.8

0.7

0.4 0.3

-1

0.2

0.4

-1

0.3

-2

0.2

0.1

0

0.1

0.2

0.3

0.4

0.5 X(m)

0.6

0.7

0.8

0.9

-3

0.1

-2

ED

0

AN US

0

0.1

y(m)

0

-1

0.2 -1.5

0.1

1

0.5

0.5

0.5

0.9

1

1

0.6

0.8

1.5

0.8

0.7

0.7

×10 -4

1 0.9

0.8

0.6

CR IP T

0.1

0

1

-4

0

0.1

0.2

0.3

0.4

0.5 X(m)

0.6

0.7

0.8

0.9

1

PT

Figure 16. Projection of the velocity field in the fiber direction in the mid-plane of plies with fibers oriented along directions: 360◦ (top-left); 330◦ (top-right); 300◦ (middle-left); 270◦ (middle-right); 240◦ (bottom-left) and 210◦ (bottom-right).

AC

CE

The simplest choice of the test function v∗ (x, z) is v∗ (x, z) = P∗ (x) ◦ T(z) + P(x) ◦ T∗ (z),

(A.2)

from which the velocity gradient is: ∇v∗ = P∗ ◦ T + P ◦ T∗ .

(A.3)

The choice of P and T in Eq. (A.3) is discussed in Section 2.3. Developing the first term in Eq. (A.1) (the other terms follow the same ratio26

ACCEPTED MANUSCRIPT nale) taking into account Eq. (4) results 

∇v∗ : ∇v ≈ (P∗ ◦ T + P ◦ T∗ ) : 

N X

j=1



Pj ◦ T j  .

(A.4)

CR IP T

It is easy to note that if matrices M(x) and N(x) depend on the in-plane coordinates x, and matrices U(z) and V(z) depend on the out-of-plane coordinate z, we have (M ◦ U) : (N ◦ V) = (M ◦ N) : (U ◦ V) . (A.5) Using this equality, Eq. (A.4) can be written as N X



∇v : ∇v ≈

j=1

{(P∗ ◦ Pj ) : (T ◦ Tj ) + (P ◦ Pj ) : (T∗ ◦ Tj )} ,

and the other terms involved in Eq. (A.1) as

j=1

(∇v)∗T : ∇v ≈

 

N n X

j=1

: (∇v) ≈

(∇v)

 

 





o

 

,

P∗T ◦ Pj : TT ◦ Tj + PT ◦ Pj : T∗T ◦ Tj

N n X

j=1

 





(A.7) o

,

(A.8)

 

P∗T ◦ PTj : TT ◦ TTj + PT ◦ PTj : T∗T ◦ TTj

ED

T



P∗ ◦ PTj : T ◦ TTj + P ◦ PTj : T∗ ◦ TTj

and ∗T



AN US

N n X

M

∇v∗ : (∇v)T ≈

(A.6)

o

.

(A.9)

PT

Thus, we finally obtain

AC

CE

4D∗ : D ≈ N n X

j=1

N X

j=1

N n X

j=1

 





 

 





 

 





 

P∗ ◦ PTj : T ◦ TTj + P ◦ PTj : T∗ ◦ TTj

N n X

j=1

{(P∗ ◦ Pj ) : (T ◦ Tj ) + (P ◦ Pj ) : (T∗ ◦ Tj )} + o

P∗T ◦ Pj : TT ◦ Tj + PT ◦ Pj : T∗T ◦ Tj

o

P∗T ◦ PTj : TT ◦ TTj + PT ◦ PTj : T∗T ◦ TTj N X 4 n X

o

A∗jk (x) : Bjk (z) + Ajk (x) : B∗jk (z) ,

j=1 k=1

27

+ +

o

= (A.10)

ACCEPTED MANUSCRIPT

A∗jk

   P∗        P∗

=

◦ Pj ,

◦ PTj , if k = 2

  P∗T       P∗T

◦ Pj , if k = 3

j

Bjk =  

 TT       TT

if k = 1 ,

AN US

if k = 2

◦ Pj , if k = 3

(A.12)

(A.13)

◦ PTj , if k = 4

◦ Tj ,

if k = 1

◦ TTj , if k = 2

   T∗T       T∗T

ED

B∗jk =

,

◦ TTj , if k = 4

j

   T∗        T∗

if k = 2

M

and

(A.11)

if k = 1

◦ Tj , if k = 3

    P ◦ Pj ,       P ◦ PT ,    PT       PT

,

◦ PTj , if k = 4

   T ◦ Tj ,        T ◦ TT ,

Ajk =

if k = 1

CR IP T

where ∀j, j = 1, · · · , N

◦ Tj , if k = 3

,

(A.14)

◦ TTj , if k = 4

On the other hand, the first term in Eq. (9) can be expressed as:

CE

with

PT

Tr(D∗ ) · Tr(D) = Tr(∇v∗ ) · Tr(∇v), Tr(∇v) ≈

N X i=1

!

Pi ◦ Ti : I.

(A.15)

(A.16)

AC

Thus, a generic term in Eq. (A.15) can be written as ((P∗ ◦ T + P ◦ T∗ ) : I) · ((Pj ◦ Tj ) : I) .

(A.17)

Now, by defining V(J) the vector form of the diagonal matrix J, Eq. (A.17) results ((P∗ ◦ T + P ◦ T∗ ) : I) · ((Pj ◦ Tj ) : I) = (V(P∗ ◦ I) ⊗ V(Pj ◦ I)) : (V(T ◦ I) ⊗ V(Tj ◦ I))+ 28

ACCEPTED MANUSCRIPT (V(P ◦ I) ⊗ V(Pj ◦ I)) : (V(T∗ ◦ I) ⊗ V(Tj ◦ I)),

(A.18)

that allows finally casting the first term in the weak form (9) as Tr(D∗ ) · Tr(D) ≈

N X

F∗j (x) : Gj (z) + Fj (x) : G∗j (z)

(A.19)

j=1

with F∗j = V(P∗ ◦ I) ⊗ V(Pj ◦ I),

(A.20)

Fj = V(P ◦ I) ⊗ V(Pj ◦ I),

(A.22)

Gj = V(T ◦ I) ⊗ V(Tj ◦ I),

G∗j = V(T∗ ◦ I) ⊗ V(Tj ◦ I).

(A.23)

AN US

References

CR IP T

and

(A.21)

[1] S. Aghighi, A. Ammar, C. Metivier, M. Normandin, F. Chinesta. Non incremental transient solution of the Rayleigh-B´enard convection model using the PGD. Journal of Non-Newtonian Fluid Mechanics, 200, 65-78, 2013.

M

[2] M.S. Aghighi, A. Ammar, C. Metivier, F. Chinesta. Parametric solution of the Rayleigh-B´enard convection model by using the PGD: Application to nanofluids. International Journal of Numerical Methods for Heat and Fluid Flows. In press.

ED

[3] A. Ammar, B. Mokdad, F. Chinesta and R. Keunings, A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids, J. Non-Newtonian Fluid Mech., 139, 153-176, 2006.

PT

[4] A. Ammar, F. Chinesta, P. Diez, A. Huerta . An Error Estimator for Separated Representations of Highly Multidimensional Models. Computer Methods in Applied Mechanics and Engineering, 199, 1872-1880 2010.

CE

[5] J.A. Barnes, F.N. Cogswell. Transverse flow processes in continuous fibrereinforced thermoplastic composites. Composites, 20/1, 38-42, 1989.

AC

[6] H.E.N. Bersee, A. Beukers. Consolidation of thermoplastic composites. Journal of Thermoplastic Composite Materials, 16/5, 433-455, 2003. [7] B. Bognet, A. Leygue, F. Chinesta, A. Poitou and F. Bordeu, Advanced simulation of models defined in plate geometries: 3D solutions with 2D computational complexity, Computer Methods in Applied Mechanics and Engineering, 201, 1-12, 2012. [8] F. Chinesta, A. Ammar, A. Leygue and R. Keunings, An overview of the Proper Generalized Decomposition with applications in computational rheology, J. NonNewtonian Fluid Mech., 166, 578-592, 2011.

29

ACCEPTED MANUSCRIPT [9] F. Chinesta, A. Leygue, B. Bognet, Ch. Ghnatios, F. Poulhaon, F. Bordeu, A. Barasinski, A. Poitou, S. Chatel and S. Maison-Le-Poec, First steps towards an advanced simulation of composites manufacturing by automated tape placement, International Journal of Material Forming, http://www.springerlink.com/index/10.1007/s12289-012-1112-9 [10] F. Chinesta, R. Keunings, A. Leygue. The Proper Generalized Decomposition for advanced numerical simulations. Springer International Publishing Switzerland, 2014.

CR IP T

[11] J.L. Ericksen. Anisotropic Fluids. Archive for Rational Mechanics and Analysis, 231-237, 1959. [12] Ch. Ghnatios, F. Chinesta, Ch. Binetruy. The Squeeze Flow of Composite Laminates. International Journal of Material Forming, 8, 73-83, 2015.

[13] J.A. Goshawk, V.P. Navez, R.S. Jones. Squeezing flow of continuous fibrereinforced composites. Journal of non-newtonian fluid mechanics, 73/3, 327-342, 1997.

AN US

[14] D.J. Groves, A.M. Bellamy, D.M. Stocks. Anisotropic rheology of continuous fibre thermoplastic composites. Composites, 23/2, 75-80, 1992.

[15] T.G. Gutowski, Z. Cai, S. Bauer, D. Boucher, J. Kingery, S. Wineman. Consolidation experiments for laminate composites. Journal of Composite Materials, 21/7, 650-669, 1987.

M

[16] T.C. Lim, S. Ramakrishna. Modelling of composite sheet forming: a review. Composites Part A: Applied science and manufacturing, 33/4, 515-537, 2002.

ED

[17] S.R. Morris, C.T. Sun. Analysis of forming loads for thermoplastic composite laminates. Composites Part A: Applied Science and Manufacturing, 27/8, 633640, 1996.

PT

[18] E. Nadal, A. Leygue, F. Chinesta, M. Beringhier, J.J. Rodenas, F.J. Fuenmayor. A separated representation of an error indicator for the mesh refinement process under the Proper Generalized Decomposition framework. Computational Mechanics, 55/2, 251-266, 2015.

CE

[19] S.F. Shuler, S.G. Advani. Transverse squeeze flow of concentrated aligned fibers in viscous fluids. J. Non-Newtonian Fluid Mech., 65, 47-74, 1996.

AC

[20] S.F. Shuler, S. Advani. Flow instabilities during the squeeze flow of multiaxial laminates. Journal of Composite Materials, 31/21, 2156-2160,1997. [21] Spencer AJM, Theory of fabric-reinforced viscous fluids, Composites Part A, 2000, 31, 1311-1321.

30