ELSEVIER

PII:SO955-7997(97)OOO60-X

Engineering Analysis with Bounda~ Elements

20 (1997) 81-89

© 1997 Elsevier Science Ltd

Printed in Great Britain. All rights reserved

0955-7997/97/$17.00

An improved boundary element method for

analysis of profile polymer extrusion

T. Nguyen-Thien a, T. Tran-Cong a'* & N. Phan-Thien b

aFaculty of Engineering and Surveying, University of Southern Queensland, Toowoomba QLD 4350, Australia

bDepartment of Mechanical and Mechatronic Engineering, The University of Sydney, Sydney NSW 2006, Australia

(Received 19 February 1997; accepted 29 May 1997)

An improved boundary element method for non-linear viscoelastic flow analysis is

reported. In this method, the domain integral representing the non-linear effects is

calculated using a more efficient approximating technique. This is achieved by first

transforming the domain integral into a form that can be approximated by particular

solutions to the original problem. These particular solutions are in fact expressed as a

linear combination of radial basis functions, whose coefficients are found by data

fitting technique. As a result, the numerical computation of the volume integral is

eliminated and a significant reduction of CPU time is achieved. The routines are tested

with simple flows and then applied to solve complex three-dimensional direct and

inverse extrusion problems of polymeric fluids, such as thermoplastic melts. Inverse

extrusion process, where an extrusion die profile needs to be computed for a given

profile extrudate, is a very important practical engineering application which is

successfully analysed by the present method. Relative to a previous BEM

implementation where the volume integral is computed directly using numerical

quadratures, a CPU time reduction ranging from 40 to 70% is achieved. © 1997

Elsevier Science Ltd.

Keywords:

Boundary element method, domain integral transformation, particular

solution, profile polymer extrusion.

1 INTRODUCTION

Boundary element methods (BEM) have become popular

techniques for solving boundary value problems in solid

and fluid mechanics. Many linear problems involving partial

differential equation (PDE) such as potential flow (the

Laplace equation), linear elasticity (Navier's equation)

and viscous creeping flow (the Stokes equation) have been

solved successfully using BEM. Obviously, for the prob-

lems involving homogeneous PDE this technique has cer-

tain advantages over the finite element and the finite

difference methods because it requires only discretization

of the boundary of the domain, thus reducing the dimension-

ality of the problems by one. Unfortunately, this advantage

is greatly offset in non-linear problems such as viscoelastic

flows. In these problems, the non-linearities can be formu-

lated into a body force term and the problem can be solved

*To whom correspondence should be addressed.

81

iteratively. At each stage of the iterative process, a linear

problem is solved, with the non-linear terms estimated from

the results of the previous iteration. 1-3 The resulting integral

equations include domain integrals. These integrals are

usually computed by numerical quadrature techniques

which require domain discretization. 4"5 However, a draw-

back of this procedure is the high demand of CPU time in

the computation of the domain integral, which is about 60-

65% of the total CPU time required for a three-dimensional

extrusion problem using a coarse mesh such as MSHI.

When a fine mesh (such as MSH2) is used, this percentage

is even higher. When the number of volume nodes is large

(for a fine mesh) this procedure is apparently inefficient.

Moreover, it is inconvenient from the numerical point of

view, especially near corner points where the stresses are

infinitely large. 6

Alternative methods of calculating the domain integral

have been used recently in order to make the BEM more

effective. 6-1° Most of the proposed techniques that have

82

T. Nguyen-Thien

et al.

been used so far belong to either methods of particular

solutions or methods of transforming domain integrals to

boundary ones. The concept in both methods is similar in

principle; the difference is in the way the total solution for

the problem is obtained. In the former method, particular solu-

tions satisfying the inhomogeneous PDE are first found, and

the remainder of the solutions satisfying the homogeneous

PDE with appropriate boundary conditions, which must be

adjusted to ensure the correct boundary conditions for the

total solution, is obtained. Then the total solution is found by

adding the particular solutions to homogeneous ones. 6'1°' I I In

the latter method, the divergence theorem or the reciprocity

theorem is applied instead to convert domain integrals (also

called pseudo-body forces) into the boundary ones. 7-932

The functions used to represent the non-homogeneous

part of the problem could be chosen among a great number

of approximants. A good choice, however, is important for

numerical efficiency. Radial basis functions in recent years

have become popular for multidimensional interpolation. 13.

Zheng

et

al. 6"14

and Zheng and Phan-Thien 12 showed some

advantages and successful results of the radial basis func-

tions in the particular solutions methods, as applied to some

inhomogeneous potential problems. We also prefer the use

of radial basis functions in this study.

The main purpose of the present work is to extend the

study of three-dimensional extrusion of viscoelastic

fluids 535 using the methods of particular solutions. 1° In

the next section we briefly recall the basic equations govern-

ing viscoelastic flow problems and their integral equation

representation, followed by a derivation of the approximate

pseudo-body force using radial basis functions. Next the

details of the numerical treatment and solution scheme are

described, which are validated with some test results involving

simple flows. The method is then applied to study complex

three-dimensional direct extrusion and inverse extrusion flows

of viscoelastic fluids. The efficiency of the method is demon-

strated by CPU time-saving ranging from 40 to 70%.

is the upper-convected derivative of the extra stress tensor,

with L the velocity gradient tensor and L x its transpose.

With the introduction of the relaxation time, a new dimen-

sionless group, the Wiessenberg number, arises. It is

defined as

Wi= k'~,

where ~, is a typical shear rate.

In BEM applications, eqn (1) in conjunction with eqn (2)

is written as

o = -

pl +

2~pD + e (4)

where 2~lpD represents the total (arbitrary) linear part of the

stress tensor and e the non-linear part. From this equation

and balance of mass and momentum, we obtain the follow-

ing integral equation:

Cij(x)uj(x) = ~nUij* (x,

y)tj(y) dI'(y)

- f~Dt~(x,y)uj(y)

dr(y)

_ JfDeJk(Y )

Ou~(x,y)

dfl(y) (5)

cgx k

where u is the velocity field, t the traction field, u*(x,y) and

t*(x,y) are known kernels (Stokeslet and its associated trac-

tion field, see, for example, Tran-Cong and Phan-Thien 5'15

and Tran-Cong16). The last integral on the RHS of eqn (5),

considered as a pseudo-body force, does not introduce any

unknown. The use of gaussian quadrature formulas to cal-

culate this integral, however, as mentioned above, could be

time consuming, and sometimes it is practically impossible

for studying complex problems such as extrusion process at

high Wiessenberg numbers, especially at refined mesh.

Alternatively this domain integral can be calculated

approximately using particular soulutions, as described in

the next section.

2 BOUNDARY INTEGRAL FORMULATION

We consider a steady, isothermal flow of viscoelastic fluid

with a single relaxation time, where the stress tensor can be

arbitrarily decomposed as: 5

a= -pl +

2~TsD + T (1)

Here p is the hydrostatic pressure which arises due to the

incompressibility constraint, 1 is the unit tensor, Os is 'sol-

vent' viscosity, D is rate of strain tensor, r is the extra-

stress tensor which is governed by a differential constitu-

tive equation of the type

X AT

~-+R=0 (2)

in which 3, is the relaxation time, R is model dependent, and

AT Jr

-- q- u.Vr - LT - TL + (3)

At at

3 APPROXIMATION OF THE VOLUME

INTEGRAL

The volume integral in eqn (5) can be converted onto the

boundary, using the divergence theorem after integrating by

parts as follows:

f ou~(x,

b i = JneJk(y)

aXk Y) dfl(y)

= [aoejk(y)uij(x,y)nk

dF(y)

fDu~.(x,

Y) 0ej~(x, y)

- Oxk

df/(y)

= J;~iejk(y)uij(x,

y)n~ dF(y) + u p (6)

Analysis of profile polymer extrusion

83

A

i

'- - n '~, x

I Ty=0. ! Vx= I •

I Tz-"O.,I "~Vy=0.

~ '- - ~, %v~.

vy=o. i . ~Tx=O. ~

~,._-o. i IIS,1

, ,trio.,

i

I ' ""~" I ~

k.Tz---O.l h

I, I / \

-.:.2.2.].:.:.:.2.2.:.:.:.:.[.I.:.2.:.2.:[:.:.:::::::~:::.2::::."

". ~./

:i!iiii!i!iiiiiii !!i!i!iii! !ii!ivx o !:!

x -_ "':::::~:2:2:2:2:::I::2::I:::2:2:I:2Vy=0.:22222~22:2:2::2:1 -.

Fig. 1. Boundary conditions for Couette flow. V and T are velocity and traction vectors.

in which u p is defined as

u p = - ~ u*'x

0ej,(x,

Y)

D ij( ,y )

" ~

dfl(y)

= - J V~(x,y)~j(y) dfl(y) (7)

where

1 0ej,(x, y) (8)

~PJ(Y)- 2

Oxk

and

I'ij(x,

y)

=

2uij(x ,

Y)

Parton and Perlin 17 have proved that the solution for u p in

eqn (7) can be found by solving the Navier equations

(

1 p ~2uP )

k, 1 --S-~ vv'u + = - 2~,(x) (9)

#

Here ~ is the Poisson's ratio. We will consider the limiting

case where 1, ---* 1/2 in the final results, since this limit is

regular.

To transform eqn (9) into the biharmonic equation we use

the 'Galerkin' form of solution given by

u p=VV.G --1 V2 G (10)

2(1 - x,)

where G is the Galerkin vector (see, for example, Phan-

Thien and Kim18). Substitution of eqn (10) into eqn (9)

yields the following equation:

V2V2G = f (11 )

where

f= -2 -~

02)

#

If a function G satisfying eqn (11) can be found then the

corresponding u p is determined by eqn (10). Note that

in this final form of the particular solution, the limit of

---, 112 is not singular. One way to find G is to approximate f

X

" ::::::: :::::::::::r::::

============================================

~-iiiiiiiV~-i:i:i!!i!:i:i:!:!:!i~-!:~

~:.:.:.:.:.:.:.:.:.:.:.:.:.:.:.:T-~.:.-:.--:-...-y ~

ITx=O. I )TxffiO. Ii ~ L

= I,,.._~ I ivy-~., I I ^'1 F

I'ql

,v~,t

IVy=°'I /

i*~:1 '- - -,i ~Vz---0.1 I~

B .............................~..:.:.......:.:........~_~j [ h

• .....-...:.:.:.:. ==========================================

Fig. 2. Boundary conditions for Poiseuille flow. V and T are velocity and traction vectors.

84

T. Nguyen-Thien

et al.

Table 1. Maximum differences of velocity and stresses between solutions obtained from the two methods and analytical solutions for

Couette flow using a coarse mesh MSI-I1 (note that P1 = NTDBEM and P2 = NTDBEM96)

Wi

% Av .u max

I %Ar~,I ma~,

%At ~:(

m~ I % S

AV

P 1 P2 P 1 P2 P 1 P2

0.1 0.51 0.50 6.65 6.51 2.48 2.62 40

0.2 0.55 0.61 5.46 5.56 2.67 2.86 42

0.3 0.67 0.71 5.26 5.35 3.12 3.27 40

0.4 0.74 0.81 5.37 5.45 3.55 3.85 43

0.5 0.93 0.92 6.05 6.25 4.07 4.26 42

0.6 1.05 1.15 6.35 6.41 5.02 5.2 40

in terms of radial basis functions if:

N

f/(r) = y.

Olin~/(~'n)

(13)

n=]

in which i = 1, 2, 3 for three-dimensional problems, N is

the number of distributed points in the domain and

Ir-r~l

rn = (14)

where r~, ~n are suitably chosen constants and ~b is a func-

tion of a single variable] 2. Now eqn (1 l) reduces to, if its

right-hand side is replaced by a radial basis function (and

therefore a particular solution is also a radial basis func-

tion), the following:

1 O(~2O( 1 0(~20~)))

~2~\ ~\g~ =~(~) (15)

When the radial basis function ~b is chosen to be exp( - ~2)

for three-dimensional problems, a particular solution of eqn

(15) is given by

1(( ' ) v/~erf(~) + exp( _ ~2) _ 2)

(16)

The Galerkin vector in (10) corresponding to this particular

solution will be given by

N

Gi(r) -= E 4ffr.)ai (17)

n=l

(see Coleman

et al. 1°

for more details).

Finally, with the substitution of the derivatives of

G i

from

eqns (15), 0 and (17) into eqn (10), the displacement u p can

be written as

1

u p ~- aiq~l - 1 - p(o~i4~2 +

~,i~,k~k4~3)

where

4,~ - q~"+ 2_4¢

r

(18)

(19)

1

~b 2 ~- -~b' (20)

F

1

~3~ ''- 7~ (21)

r

in which ? i =

O~'lOx,.

4 NUMERICAL IMPLEMENTATION

A decoupled technique, similar to that reported by Bush 3 is

implemented here. The procedure for finding the extra-

stress (the non-linear part) is similar to that of Tran-Cong

and Phan-Thien 5 and the details will not be repeated here.

The new feature of the present method is the way in which

the domain integral is computed. The data fitting techniques

will be used instead of numerical quadrature techniques.

Ifain, fin

and rn in eqns (13), (14), 0 and (17) are known

from eqns (19)-(21) and (18), the pseudo-body force b in

eqns (6) and (7) is also determined. In principle, the three

parameters can be found by using available non-linear data-

fitting techniques. However the calculation is performed

only by iteration and therefore such a scheme could be

Table 2. Maximum differences of velocity and stresses between solutions obtained from the two methods and analytical solutions for

Couette flow using a fine mesh MSH2 (note that P1 = NTDBEM and P2 = NTDBEM96)

Wi

%AY ~(maxl

%Ar~(,,a~l %Arx:,

,,ax~ %SAV

PI P2 PI P2 PI P2

0.1 0.13 0.14 2.25 2.41 1.48 1.51 55

0.2 0.19 0.21 1.55 1.70 1.66 1.81 55

0.3 0.28 0.30 1.86 2.15 2.08 2.12 56

0.4 0.37 0.39 2.92 3.02 2.25 2.32 58

0.5 0.41 0.42 3.42 3.61 2.55 2.67 58

0.6 0.52 0.58 3.55 3.82 2.80 2.93 58

Analysis of profile polymer extrusion

85

Table 3. Maximum differences of velocity and stresses between solutions obtained from the two methods and analytical solutions for

Poiseuille flow using a coarse mesh MSH1 (note that PI = NTDBEM and P2 = NTDBEM96)

Wi

%Avx(max) %ATxx(max) %AT"xz(max ) %SAV

P1 P2 P1 P2 P1 P2

0.1 1.43 1.54 1.10 1.05 1.81 2.01 41

0.2 1.42 1.50 1.05 1.06 1.85 2.2 42

0.3 1.30 1.43 1.02 1.06 1.80 2.3 42

0.4 1.10 1.05 1.04 1.05 1.90 2.3 43

0.5 1.00 0.92 1.02 1.04 1.85 2.1 43

0.6 0.85 0.71 1.00 1.02 1.85 2.1 43

time consuming. Alternatively we can use linear data-fitting

techniques by choosing/3, and rn

apriori

and leave o~in to be

determined. Zheng

et

al. 6"14 and Coleman

et al.~°

obtained

good results following this approach. Recommendations for

/3~ and r, are also given. In this work rn are chosen to be

suitably distributed points in the domain (here we use the

points resulted from a typical FE-type discretization). The

/3n for each point are chosen to be the average distance from

its immediate neighbours multiplied by a weighting factor,

and the weighting factor can vary in a range from 0.5 to 2.0

by numerical experiments. ~° In order to obtain accurate

solutions, i3~ (or the weighting factor) must be carefully

chosen. This matter will be discussed in the next section.

To find oti, one can use either a collocation method or a

least-square fitting method. Both of those methods result in

a set of linear algebraic equations in which c~i, are unknown.

The use of the gaussian distribution as the basis functions

makes the system matrix effectively sparse, for which either

the Gauss elimination (GE) or an iterative solution method

can be used, such as the conjugate gradient method

(CGM) 19 which has been known to be more efficient for

large systems. The results shown in the next section were

obtained with both GE and CGM methods. The differences

in the results by the two methods are in the range of 0.01-

0.05%, which must be due to the iterative nature of the

CGM method.

5 NUMERICAL RESULTS

5.1 Test problems

In this section we show some test results obtained from the

procedure described above. The code was tested in some

simple flows (Couette and Poiseuille flows) for which

analytical solutions are known. The pseudo-body force

was calculated using gaussian quadrature formulas; 5 the

same procedure is adopted here. For each problem a

coarse mesh, MSH1 (with 73 boundary nodes and 79

domain nodes) and a fine mesh, MSH2 (with 709 boundary

nodes and 1549 domain nodes) were used.

Couette and Poiseulle flows of an upper-convected Max-

well fluid are used to test the program. The boundary con-

ditions for the two flows are illustrated in Figs 1 and 2. The

velocity and stress fields are obtained analytically by sol-

ving the field equations with appropriate boundary con-

ditions. 21 Field variables are non-dimensionalized

according to

t~ t X'---- X_, u'----U, 7"' 7"

-----~' a -------U

7--

a

where t is time, x is the position vector, u is the velocity

vector, r is the extra stress tensor and ~ is the viscosity. X is

the fluid relaxation time, a is a typical length and U is a

typical flow speed. Then, for Couette flow, the velocity

field v and stress field 7 are given as

V~

[!

t/p

, 7 ~

L2wi%~

10

For Poiseuille flow, we have

V=

~h 2

~-~(1 - z 2)

0

0

[

2~pWiz 2 0 -z%

0 0

L - Z~p 0 0

In these formulas 5' denotes shear rate,

Wi

= ~,'i' is the

Table 4. Maximum differences of velocity and stresses between solutions obtained from the two methods and analytical solutions for

Poiseuille flow using a fine mesh MSH2 (note that P1 = NTDBEM and P2 = NTDBEM96)

Wi

%Avx(max) °~AT"xx(max ) °~mTxz(max ) %SAV

P1 P2 PI P2 PI P2

0.1 0.83 0.92 0.23 0.23 0.65 0.59 55

0.2 0.80 0.83 0.26 0.25 0.69 0.63 55

0.3 0.78 0.75 0.25 0.26 0.71 0.65 56

0.4 0.74 0.67 0.27 0.28 0.71 0.66 56

0.5 0.60 0.58 0.30 0.31 0.73 0.66 56

0.6 0.58 0.41 0.31 0.30 0.81 0.67 56