Numerical prediction of flow induced orientation in anisotropic suspensions. Application to injection molding of fibers reinforced thermoplastics

Journal of Materials Processing Technology, 32 (1992)429-438

429

Elsevier

Numerical prediction of flow induced orientation in anisotropic suspensions. Application to injection molding of fibers reinforced thermoplastics. A. POITOU a, F. CHINESTAb, R. TORRES b a Laboratoire de Mtcanique et Technologic, E.N.S. Cachan / C.N.R.S. / Universit6 Pads VI 61, avenue du prtsident Wilson, 94 235 CACHAN Cedex (France) b Escuela Tecnica Superior Ingenieros Industriales. Universidad Politecnica de Valencia 46071 Valencia. Spain Abstract

An approximate method is given to compute the flow of anisotropic suspensions in plane mold with constant thickness. It is shown that the three-dimensional equations reduce into a bidimensional problem which unknown is the average of the velocity field through the thickness. The approximate equations are solved with a pseudo-spectral collocation method using a Chebychev polynomial interpolation. A stream function formulation permits to account for the incompressibility constraint. Numerical results show the limitation of the often used Reynolds approximation to deal with injection molding of fiber reinforced thermoplastics.

1. INTRODUCTION Modeling flow induced orientation of short fibers reinforced thermoplastics is of major interest in composite processing. It permits both to avoid orientation defects in injected molded pieces and to define the mechanical behavior of the composite, which is known to be strongly coupled with the fiber orientation [1-2]. This modeling is commonly achieved in the framework of dilute or semi-dilute suspensions of non spherical particles in a newtonian fluid [3-4]. Resulting equations involve the coupling of an elliptic boundary value problem with a convection type equation. The elliptic problem is associated with the equilibrium equations whereas the convection equation describes the time evolution of the anisotropic viscosity tensor. Computing the 3D solution for arbitrary geometries seems today out of reach, but it is still possible to evidence simplified equations which are valid for thin mold geometries and uniform thickness. After Tucker's paper [4] who showed the existence of different flow regimes according to the relative order of magnitude of the different terms in the equations, we proposed [13] a systematic procedure which enables us to make explicit the simplified equations that govern the different flow regimes. This paper gives some numerical results which are obtained with a multidomain pseudo-spectral method. It is shown that even in thin mold geometries, the velocity field can be strongly modified by the presence of fibers in the suspension. The problem to be solved is shown on figure 1. Of~l denotes the free surface on which a given force is assumed to be prescribed. 3f22 denotes the lateral edges on which the velocity field is 0924-0136/92/$05.00 © 1992ElsevierSciencePublishers B.V. All rights reserved.

430 prescribed (i.e. the mold gate and the lateral walls of the mold). 0f23 and 0f2 4 denote respectively the upper (z = h) and lower surface (z = - h) of the mold. On Of13 and OD.4, a zero velocity field is prescribed. The midplane of the mold is 5"..The trace of 0f21 and Of~2 on E is 0El and OE2 We assume that Of~lUOf~2UOf23UOD.4 = 0f2. The prescribed velocity field on 0f22UO~3UO~4 is vg. The prescribed force on Of21 is Fg. The aim of this paper is to find approximate velocity and stress tensor fields of this problem on ft.

~g i~f24 Figure 1 Mold geometry In section 2, the general equations of anisotropic suspensions are summarized and the approximate model is derived for the flow in thin mold geometries. The numerical procedure is described in section 3. Numerical results are discussed in section 4. Notations: Vectors are written in bold characters Matrices are written in shadow characters Tr denotes the trace operator T denotes the transpose operator

: e.g. : e.g. : e.g. : e.g. : e.g.

u (velocity field) ~ (extra stress tensor) Tr(A) = Aii A T = Aji u T denotes the line vector whereas u denotes the column vector

Scalar product is written : u T v = ui vi Diadic product is written : uvT=uivj Strain rate tensor : d (v) = {grad v + (grad v)T}/2 Every other particular notations are introduced throughout the text. 2. G E N E R A L E Q U A T I O N S

2-1 CONSTITUTIVE LAWS FOR ANISOTROPIC SUSPENSIONS Molten thermoplastics composites are commonly modelled as suspensions of axisymetric particles. A complete derivation of the equations can be found in Giesekus [5], Hinch and Leal [6], Batchelor [7-8], Dinh and Armstrong [9]. Two main steps permit to derive macroscopic constituve equations for anisotropic suspensions: (i) a now standard homogeneisation

431 procedure which leads volume average quantities if the orientation of the particles is perfectly defined and (ii) a statistic description of the orientation. Volume average

Basic equations are obtained with a volume average procedure. Let ~ be a representative volume of our macroscopic scale, which contains many particles located in Clfi. Let 'L u and d(u) denote the microscopic extra stress tensor, velocity vector and strain rate tensor and n be an outwards normal vector to the particles' boundary 3ClPi. For a given orientation of every fiber in the suspension and if inertia terms can be neglected, the corresponding macroscopic variables ']1",v, d(v), are defined by: /.

v = 1---/clpu dv ,/,¢

(2-1)

q,e

.y~, i

Ct'PJ
q'P

-5~

C~ i

i

If the ambient fluid is assumed to be newtonian, the microscopic constitutive law within the fluid reads: 'l = 2 rl d(u)

(2-4)

This leads to the expression of the macroscopic pressure and extra stress tensor: /, P=~-

Io pdv

and

T = 2rl d(v) +:gp

(2-5)

P is the volume averaged pressure and ~p describes the contribution of the particles to the macroscopic stress: t"

Zp = _L_ ~ C~P

/

'l n x T dv

(2-6)

i

Determination of~,p requires to solve a microscopic problem (e.g., Jeffery problem [10] for dilute suspension theory or interacting slender bodies problem for semi-dilute theories [8-9]). Orientation distribution

Particle orientation is described through the probability distribution of fiber orientation ~g(p). If p denotes a unit vector aligned along the symmetrical axis of the particle, if the particles have a quasi-infinite aspect ratio and if brownian motion can be neglected, equations of evolution for p and ~g read: d__pp= g r a d (u) p - p V d ( u ) p p

dt

(2-7)

432 The macroscopic and orientation averaged extra stress tensor ~ is then defined by: /.

/ T V(P)dp

(2-9)

d

L e t , and a4 be the second order and forth order orientation tensor. •

A straightforward calculation yields: ~ t = (8radv) • + ,(grad v)T - 2 , 4 d

(2-11)

One can approximate a4 with a quadratic closure relation [11]: •4 = z ® a

(2-12)

Then a general expression for ~ can be written which is the constitutive relation of the suspension: = 2rlI{d + Nptr(ad)a+Ns{da+ad}}

(2-13)

ddat= ( g r a d v ) , + a (grad v)T- 2 tr (ad)d

(2-14)

This constitutive equation has the same basic structure as that of anisotropic fluids introduced by Hand [12]. Thus, though the microscopic model is generally not valid for thermoplastics composites because the fiber concentration is too large for the composite to be considered as a dilute or even a semi-dilute suspension, equations (2-13), (2-14) remain valid as a first linear approximation, provided that the rheological parameters rib Np and Ns are determined experimentally. It is to be noted that for concentrated suspensions, Ns can usually be neglected before Np [4]. This approximation is made in the following. 2.2 SIMPLIFIED EQUATIONS FOR THIN MOLD GEOMETRIES In this section, starting from the virtual work theorem and assuming that a shear strain kinematic hypothesis is valid, a simplified solution is evidenced. Details of the method can be found in [14] where a comparison with an asymptotic expansion analysis is also given. Because of the non newtonian behavior of the fluid, there is no minimum principle that expresses both the constitutive relation and the momentum equation. However, the momentum equation can be written:

V u*¢ el/'*= {u*, div (u*) = 0, u* = 0 on ~f~2u~f~3u~f24}

(2-15)

The other equations governing a steady solution are written at a local stage: n ~ elf = {u, div (u)= 0 in f2, u = ug on

~£"~2k.l~23k.)~k'~4}

(2-16)

433 0

= 2nI{d + Nptr(ad)a}

(2-17)

(u T grad) a = (gradu) a + a (grad u)T- 2 tr (ad) d

(2-18)

A simplified solution is looked for in a subset q/'l of q/' in writing the momentum equation for every u* in q/'l*: q/'l ={u~R3, u =f(z)v(x'Y)= 311 z 2[/ / V h2][Vy] x l 2 ' div(v)= 0, v = l f f

2 1 h2]~v} ' div(v*)= 0,

Vg dzon OZ2}

= 0 on OY'2

A straightforward calculation shows that the locally aligned orientation is an exact solution of equation 2-18 for every u in %°1. The stress tensor then reads: o

tT ~

with "1"=2 f(z)Th D + N p T r [ V V T D / v v T I ~vXv / v T v I '

t='qI

v, ~ = 0

(2-19)

The momentum equation is reduced into the following form:

+

~

f(z) FdV*

=0

Vv ° ~

~*

(2-20)

I

Defining the pressure p as a Lagrange multiplier associated with the incompressibility condition and calculating in the above equation the integrals with respect to z leads to the formulation of the simplified equations: Weak form: -[~ThhTr{{ +

¢-

D+NpTrlvvTD/vvT/D*[+6~Ivv',vTv , vTv, , h

2h Fg v* = 0 V v* satisfaing v* = 0 on ~E2 and div (v*) = 0 in ~;

1

(2-21)

434 Local equations:

11I h div {D + NpTr IvvTv vy DII vvYv vy ~_ 6 1]I J - ~ v = 2h g r a d ~ and div (v) = 0 inZ J

3

1]I h [ D + NpTr Iv vTD] v vT[ n - 2h ~ n = 2h F--g on ~Z1

t

~vTv

I VTV J

v = vg on 8~2

(2-22)

Discussion: Let e be the mold aspect ratio (i.e. typical length / thickness) I f N p < < e -2 ,

5~-1]ihdiv{D+NpTrfV vTDIv vTl<< 111 [vTv t v T v l 6~--v

(2-23)

The flow is then governed by shear stresses and it follows the Reynolds equation. if Np___E_2,

~ 1]i h div {D + NpTr Iv vTDI v vTI = 1]I 1VTV ! vTv I 6-~-V

(2-24)

Tension terms and shear stress terms are of the same order of magnitude. This case yields an anisotropic equivalent of the Barone ~indCaulk model [ 14], which does not imply a flat velocity profile through the thickness. If Np >> e -2 ,

52~4-111h div {D + NpTr [v vTDI v vT I >> 6 TII v [ vTv / vT v J h

(2-25)

Only tension terms are to be accounted for. This leads to an anisotropic membrane model. The physical meaning of this situation would be that the fibers are so intimately intricate that the existence of solid upper and lower walls in the mold does not alter the flow significantly. If the suspension is assumed to slide along the upper and lower mold surfaces instead of the sticking contact that is considered in this paper, this leads to the fourth flow regime described by Tucker

[4].

3. N U M E R I C A L

ALGORITHM

3.1. STREAM FUNCTION FORMULATION Since the approximate model is bi-dimensional and because of the incompressibility condition, a stream function formulation has been chosen. Though this method does not permit to deal with arbitrary boundary conditions, it is appropriate to evaluate the model behavior in simple reference geometry. Let V denote the stream function, the velocity field can be written: 8~F 8V u =By' v = - 8y (3-1) Equation (2-22) then reads:

A2V--2--A

2h 2 ~

=

0

(3-2)

435 Where G is a non linear function of partial derivatives of V up to order 3. We only considered the case for which the boundary condition is a prescribed velocity field on the boundary: = V,

dv ~

_

=q

on bY.

(3-3)

On the solid lateral edges, ~t is assumed to be constant to account for the non sliding condition. At the entrance and at the exit V is an arbitrary function provided that: I z a~-~--ds= 0,

where t denotes a unit vector tangent to the boundary

(3-4)

3.2. PSEUDO SPECTRAL METHOD: CHEBYSHEV COLLOCATION The pseudo spectral method has been chosen for two main reasons: (i) this method is easy to implement for high order equations in simple geometries and permits to introduce one degree of freedom per point only, (ii) its order of convergence is exponential for smooth solutions. It is thus possible to check easily the convergence of the method. The equations are solved on a square reference domain [-1,1]x[-1,1]. This requires to modify the equation on E through the Jacobean transformation. A basis of the approximation space is chosen to be the tensorial product of Chebyshev polynomials up to order N in x direction and M in y direction: N

M

1q/(x,y)N,M= ~ ~ Cij Ti(x) Tj{y) i=o j=o

(3-5)

Or in an equivalent nodal form: N

M

lqJ(x,y)N,M= ~ Y. ~xi,Yj) Ni(x) Nj(y) i=0 j=0

(3-6)

Where Ni (or Nj) is a Chebyshev shape function of order N (or M). If xj denotes the-collocation point defined by: xj=-co~-~-}C~TN(Xj)=0,

then

Ni(xj)=Sij

(3-7,

The expression of the derivatives at each collocation point is given by:

Lax

iffi~"r~"YJ4"-?"~ I

i=0j=

FaqNAy)] t - - t

L Oxp ..kx=x,iL ~ I _ky=y,)

(3-8)

Where the derivatives of the shape functions at each collocation point are known and tabulated

number [15]. Equation (3-2) is solved with a collocation method in requiring that this equation is satisfied at each point ( xi, yj ). This method is equivalent with a variational formulation (Chebyshev shape

436

functions and Gauss Lobato Chebyshev quadrature). The non linearities are solved using a standard iterative Newton-Raphson scheme. The convergence is increased by a relaxation and step by step technic. The initialization is achieved in solving the associated problem with N o = 0 (Stokes problem). If the two boundary conditions are exactly prescribed at each boundary node, the resulting problem is ill posed. The inner equations located just near the boundary are thus eliminated. This leads to a well posed problem which has been shown to be convergent for a similar biharmonic equation. Multi-domain analysis To account for more complex geometry a multi-domain analysis has been used. On each subdomain, the stream function is approximated as above. The compatibility condition between to neighboring sub-domains has been written in enforcing the continuity of the function as well as its derivatives up to order 3. 4. RESULTS On a rectangular domain, at the first step of the Newton Raphson scheme, the residual reaches high values near the boundary and particularly at the comer of the domain because of the existence of the singularity which results from the discontinuity of the velocity field at these points. The number of iterations depends essentially on Np. For example, if Np is of order 1/e 2, the convergence is obtained after 5 or 6 iterations. On non rectangular domains, with the multi-domain procedure, the convergence cannot be reached directly for large Np values (Np of order l/e2). However, this convergence is made possible with a step by step method. The values of N 0 are increased progressively and at each step, the initial value of the stream function is taken to be the solution of the problem at the previous step.

A 0.00

0.06

0.12

0.18

0.24

0.30

0.36

0.42

0.54

0.48

i i i i i i I I I I l l l l l l l l l l l

0.10

I

I

G60 t

0.C8 --

0.14

0.0'7

0.06

-

- - - - - - ~ - -

0.66 I

I

I

0.72 I

I

I

0.78 I

I

I

0.84- 0.90 I

I

01

0.06

0.06

-0.14

--

I

I

0.96 I

-

1

0.06

-0.02

-0.06

-0.06

--

- -0.1

-0.1+

I

0.02

002

-0.1

I

0.1 -

-0.02

-0.06 -

I

0 . I 'i - ~

0 14

-

0.02

0.03

I

" O. I

-0.02

0.04

i

-0.1 -

-0.14-

0.01 O.OQ

= l l l l l l l l l l l l l l l l l l l l l i

0.00

0.06

0.12

0.18

0.24

0.30

0.36

0.42

0.48

I I ~ t il 0.54 0 6 0

III 0.66

t , I I I I I I I I I I 0,72 0 . 7 8 0.84- 0 . 9 0 0.96

8 Figure 2 Contour of the stream function

I--J

437 NP 0

,.

//

/), <2

~

J

r,l r,l

//

~

.

o

4

-,.~z-e.e,-o.6a-~.E'2-0.36-e.2o

.o.e4

0.~2

o.l~a

m44

o.~

e.r6

e.9~

A8 Figure 3 Average velocity through the thickness, along the line AB Figure 2 and 3 show the results that were obtained for a uniform inlet and outlet velocity profile on a rectangular (i.e. parallelepipedic) domain. Degree of interpolation is 9 in each direction. If this degree is increased, the solution is not altered significantly but the little oscillations of the stream function are smoothed. Figure 3 shows that if Np is of order 1/e2, the solution is different from the Reynolds solution which is computed with N o = 0. Figure 4 and 5 give the same results with the multidomain analysis. 3 spectral elements tiave been chosen in which the degree of interpolation is 9 in each direction. The convergence is also obtained for Np of order 1/e2. Figure 5 illustrates the significant difference between the profile of the average of the velocity through the thickness in the case Np = 0 and Np = 1/~2 (ie Np = 100 for this geometry). 0.00

0.17

0.2 ~-

0,48

0.72

CI

0.96

1.20

1.4.

1.68

1.92

~ i I I I I I I i i i 018i i~1i4..1I I . . I. . [ . I. . I I I i i i

0.11

~

2.1e

~,40

I I I 11

2.64

2.86

3.12

3.36

,~.~0

I I i i i i i i i i i i i i I

3.84

i i ~0.1711

-

=

0,06

2,;~

;.2 D.O=t

i i i i i 0"000.00

0.2l

0,4~

0:/2

, l.~....L_l.._r J i i i ~J i i~"" 0..06 t20 1.44, 1,88 1.92 2.16

l_-~::3z2AO

2,64

o.oo 2.86

D= Figure 4 - Contour o f the stream function

3.12

3.36

3.ElO

3.84

438 r,p

//

//

/

\\ oD

t'' o.~1

I'''1'''1 il. ~;l ~.o'3

IJ.~i

''l'''l e.~

'l' ~.~

r' I~l. ll~

' I ' ' ' I ' ' ' I ' 'll' ' ' 1 ' ' I' II. I I t~, i:1 o . 1 ~ i1.16 i~,16 0.19

Figure 5 - Average velocity through the thickness, along the line CD 4. CONCLUSION A simplified model has been introduced, which takes into account two significant properties of short fibers composites, namely the existence of non zero normal stress differences in shear strain and the flow induced anisotropy. The kinematical hypothesis which leads to assume a shear strain velocity profile through the thickness enables us to reduce the three-dimensionnal equations into a two-dimensionnal flow problem. This two-dimensionnal problem has been solved numerically with a spectral element method for simple geometries and boundary conditions. However, extensions of this method to more general conditions (e.g. existence of a free surface), are possible. References 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

G.G. Lipscomb, M.M. Denn, J. Non Newtonian Fluid Mech. 26 (1988) 297-325 K. Chiba, K. Nakamura, J. Non Newtonian Fluid Mech 35 (1990), 1-14 M.A. Bibbo, S.M. Dinh, R.C. Armstrong, J. Rheol. 29 (1985) 905-929 C.L. Tucker, J. non newt. fluid mech. 39 (1991) 239-268 H. Giesekus, Rheol. Acta 2 (1962) 50-62 E.J. I-linch and L.G. Leal, J. Fluid Mech. 57 (1973) 753-767 G.K. Batchelor, J. Fluid Mech. 41 (1970) 545-570 G.K. Batchelor, J. Fluid Mech. 46 (1971) 813-829 S.M. Dinh and R.C. Armstrong, J. Rheol. 28 (1984) 207-227 G.B. Jeffery, Proc. Roy. So(:. A-102 (1922) 161-179 E.J. Hinch and L.G. Leal, J. Fluid Mech. 76 (1976) 187-208 G.L. Hand, Arch. Rat. Mech. Anal. 7 (1961) 81-86 A. Poitou, J. non newt. fluid mech., (submitted) M.R. Barone, D.A. Caulk, J. Appl. Mech. 53 (1986) 361-371 D. Gottlieb, Y. Hussaini, S. Orszag, in "spectral methods for partial differential equations". SIAM 1984