• No results found

BEM simulation for the pressing of glass

N/A
N/A
Protected

Academic year: 2021

Share "BEM simulation for the pressing of glass"

Copied!
15
0
0

Bezig met laden.... (Bekijk nu de volledige tekst)

Hele tekst

(1)

BEM simulation for the pressing of glass

Citation for published version (APA):

Wang, K., Mattheij, R. M. M., & Morsche, ter, H. G. (1999). BEM simulation for the pressing of glass. (RANA : reports on applied and numerical analysis; Vol. 9923). Technische Universiteit Eindhoven.

Document status and date: Published: 01/01/1999

Document Version:

Publisher’s PDF, also known as Version of Record (includes final page, issue and volume numbers)

Please check the document version of this publication:

• A submitted manuscript is the version of the article upon submission and before peer-review. There can be important differences between the submitted version and the official published version of record. People interested in the research are advised to contact the author for the final version of the publication, or visit the DOI to the publisher's website.

• The final author version and the galley proof are versions of the publication after peer review.

• The final published version features the final layout of the paper including the volume, issue and page numbers.

Link to publication

General rights

Copyright and moral rights for the publications made accessible in the public portal are retained by the authors and/or other copyright owners and it is a condition of accessing publications that users recognise and abide by the legal requirements associated with these rights. • Users may download and print one copy of any publication from the public portal for the purpose of private study or research. • You may not further distribute the material or use it for any profit-making activity or commercial gain

• You may freely distribute the URL identifying the publication in the public portal.

If the publication is distributed under the terms of Article 25fa of the Dutch Copyright Act, indicated by the “Taverne” license above, please follow below link for the End User Agreement:

www.tue.nl/taverne

Take down policy

If you believe that this document breaches copyright please contact us at: openaccess@tue.nl

(2)

BEM Simulation for the Pressing of Glass

K. Wang, R.M.M. Mattheij & H.G.ter Morsche EMail: wang@win.tue.nl

Abstract

The main objective of this research is to develop reliable numerical methods to simulate the pressing of glass. The glass can be modelled as a Stokes ow with three types of boundaries: free, xed and \moving" boundaries. Because of axi{symmetry we actually have a 2-D problem. This problem is solved by rewriting it as a boundary integral equation. A boundary element method is then employed to solve the integral equation of the Stokes ow and time integration is carried out by a kind of predictor-corrector scheme.

1 Introduction

In Mattheij at el. 12] a model was given to describe the ow of glass in a mould with partially free and partially prescribed boundaries. Here we shall consider a numerical approach to simulate the actual solution by a boundary element method.

If we let v= (vi) be the velocity vector, p the pressure and  = (ij) stress

tensor, then the ow of the glass can be seen to be described by the continuity

equation

rv= 0 (1)

and themomentum equation

v;rp= 0: (2)

Whereas boundary conditions are given by

n=;p 0

n on the free boundary, (3)

(

vn= 0

nt=;mvt

on the mould boundary, (4)

and (

vn=vpn

nt=;p(v;vp)t

on the plunger boundary. (5)

Here n = (ni) is the outward unit normal, t= (ti) is the unit tangential

direc-tion, vp is the velocity of the plunger, m and p are parameters indicating the

roughness of the mould and the plunger respectively.

(3)

free boundary plunger boundary mould boundary

Figure 1: The boundaries

dx dt

=v(x(t)): (6)

Now our objective is to solve problem (1) { (6). In section 2 we solve the Stokes problem in an axi-symmetric domain by a Boundary Element Method and In section 3 we further deal with equation (6), by developing a predictor-corrector scheme. Finally we will assess the method by considering some actual simulations in section 4.

2 BEM for an Axisymmetric Stokes Flow

There are two approaches for deriving the governing integral formulation for axi{ symmetric problems based on hydrodynamic potentials of single{ and double{ layers. Both methods are leading to the same equation. The rst one is to obtain the integral equation by using the axi{symmetric fundamental solution based on ring forces(cf. Becker2]). The second approach is to apply the fun-damental solution derived from a point force to obtain the a Cartesian version of the three{dimensional integral equation. Subsequently cylindrical coordinates are then substituted in this formulation. Here we use the latter method.

2.1 The Boundary Integral Equations

Let us denote the uid region by  and its surface by S. The fundamental

(4)

8 < : uki(xy) = 1 8 h ik jx;y j+ (x i ;y i )(x k ;y k ) jx;y j 3 i  qk(xy) = xk ;y k 4jx;y j 3  (7) i.e. (ukqk) satisfy ( ruk = 0 uk;rqk =;(xy)ek:

where ek is thek{th unit vector of an Cartesian coordinate system.

By using the divergence theorem, we can obtain theGreen's formula for the

Stokes problem, i.e., if v and v

0 are divergence free,

p and p

0 are smooth scalar

then R ( vi; @p @xi) v 0 i;( v 0 i; @p0 @xi) vi]d =R ; ij(v p)njv 0 i;ij(v 0 q 0) njvi]d;

where ij(v p) =;pij+vij+vji and ij(v

0 p 0) = ;p 0 ij+v 0 ij+v 0 ji: By replacing v 0 i and p

0 by the fundamental solutions

uki and qk, identifying

vi and p with the solution to Eqns. (1) and (2), and letting the source point

approach the boundary, we obtain the following boundary integral equation(cf.

Ladyzhenskaya7], Pozrikidis14] and Brebbia et al. 3])

cij(x)vj(x) + Z Sqij(xy)vjdS y= Z Suij(xy)bjdS y  (8)

where the kernels qij is equal to :

qij(xy) = 3(

xi;yi)(xj;yj)(xk;yk)nk

4jx;y j

5

: (9)

The coecientscij depend on the position of the pointx, but in the BEM, we will

show that it is not necessary to know the analytical expressions of cij although

they are available (cf. Brebbia at el. 3]).

To obtain the integral equation for the axi-symmetric case, we reformulate

the representation above by employing cylindrical coordinates (rz), i.e.

y= (y 1 y 2 y 3) T = (rcosrsinz)T:

Because of the rotational symmetry, we only have to determine vr and vz at

the intersection of the surfaceS and (say) the half{space= 0. This intersection

curve is denoted by ;(see Figure 2), therefore dS=rdd;. Letx= (R 0Z)T 2

;. After successive substitution of cylindrical coordinates and integration along

the {direction of Eqn. (8) we obtain

ccij(xc)vcj(xc) + Z ; rqcij(xcyc)vcj(yc)d; = Z ; rucij(xcyc)bcj(yc)d; (10)

(5)

r

Γ

z

S

Figure 2: ;(the solid lines) and S

where the superscriptcstands for cylindrical andi,jare now either 1 or 2 source

point xc= (R Z)T, eld pointyc= (rz)T velocity vc= (vci) = (vrvz)T, stress

vector bc= (bci) = (brbz)T. the kernelucij can be written as

ucij = 1 2 p a+b h E(k) a;b  A 0 ij+abA 1 ij+2a 2 ;b 2 b2 A 2 ij  ; K(k) b  A 1 ij+2 abA 2 ij i  (11) where a=r 2+ R 2+ c 2  b= 2rR  c=Z;z k = s 2b a+b :

K(k) and E(k) are the complete elliptic integrals of the rst and second kind

respectively, i.e. K(k) := Z  2 0 d q 1;k 2sin 2  E(k) := Z  2 0 q 1;k 2sin 2 d : (12)

The coecient matrices An are dened by

A 0 =  ; 1 2 b cR ;cr a+c 2 ! A 1=  2a;c 2 ;cr cR ;b ! A 2 =  ; 3 2 b 0 0 0 ! :

While the kernel qcij can be represented by

qcij = 1 (a 2 ;b 2 ) p a+b h K(k)  ;B 0 ij;abB 1 ij+2a 2 ;3b 2 b2 B 2 ij+a(8a 2 ;9b 2 ) b3 B 3 ij  +E(k) a;b  4aB 0 ij+a2 +3b 2 b B 1 ij+2a(3b 2 ;a 2 ) b2 B 2 ij; 8a 4 ;15a 2b2 +3b 4 b3 B 3 ij i :

Here coecient matrices Bn are dened by

B 0 =  ; 1 2 db R dc ;rdc dc 2 !  B 1=  de; 1 2 bR nr (R 2 nr;rd)c (d;rnr)R c R c 2 nr !  B 2 =  (enr;rd)R ; 1 2 bcnr R 2 cnr 0 !  B 3 =  ; 1 2 bR nr 0 0 0 ! 

(6)

where d=;rnr+cnze=R 2+

r 2 and

n= (nrnz) is the unit outward normal

of boundary ;.

We remark that if the point x is lying at the z{axis, i.e. R = 0, the above

formulae are not applicable because of b= 0. However, direct computation leads

to ucij= 1 8a 3 2 h 2(A 0 ij+A 2 ij) i  and qcij = 3 4a 5 2 h 2(B 0 ij+B 2 ij) i :

Note that if the domain includesz{axis then ; is an open curve.

2.2 The Boundary Element Method

Although physically 3{D, the BEM implementation of Eqn. (10) is similar to a 2{D problem. To obtain a BEM formulation, we rewrite the boundary integral

Eqn. (10) (we remove the superscript cfor simplicity)

c(x)v(x) + Z ; q(xy)v(y)d; = Z ; u(xy)b(y)d;: (13)

The boundary is divided into segments ;k, i.e.

; =X

k ;k

 (14)

and the velocity and stress vector are represented by shape functions

v= X j vj j b= X j bj j: (15)

By substituting Eqns. (14) and (15) into Eqn. (13) we have

c(x)v(x) + P jP kR ; k q(xy) j(y)d;vj = P jP kR ; k u(xy) j(y)d;bj

For a particular pointxi, denotingc(xi),v(xi) byci andvi respectively, we thus

have civi + P jP kR ; k q(xiy) j(y)d;vj = P jP kR ; k u(xiy) j(y)d;bj: (16)

Now dene the matrixHby

Hij := X k Z ; k q(xiy) j(y)d; +ciij (17)

(7)

and the matrix Gby Gij:= X k Z ; k u(xiy) j(y)d;: (18)

This results in the following square full rank system of linear algebraic equations, denoted by

HV =GB: (19)

Here V, B are vectors consisting of the velocity and stress vector respectively.

Note that every element of Hand G is a 22 submatrix, and every element of

V and B is a 2N 1 subvector. If there are N collocation points, and probably

include M corners, then H and G are actually 2N 2N and 2N 2(N +M)

matrices respectively, whereasV and B are 2N 1 and 2(N +M)1 vectors.

In order to build H and G the elliptic integrals (12) that occur in the

co-ecients are approximated by using a series representation(cf. Becker1] and Cody5]), i.e. K(k)log(4) + P ni=1 aii+ log( 1 )1 2 + P ni=1 bii] E(k)1 + P ni=1 cii+ log( 1 )P ni=1 dii (20) where = 1;k 2 :

When the element integrals become singular, we use Telles' transformation (cf. Partridge at al. 13] and Telles17]) to remove the singularity. More precisely we consider the following integral

I = Z

1 ;1

f()d

in which f() is singular at a point : One can choose a transformation

=a 3+ b 2+ c +d such that 8 > > > > < > > > > : (;1) =;1 (1) = 1 d dj   = 0 d2 d2 j   = 0

by which the coecients abcd can be determined. ThenI becomes

I = Z 1 ;1 f() d d d :

If f() has only a logarithmic singularity at , then as a function of ,f()dd is

(8)

As for the diagonal submatrix of H, we rst use a rigid-body motion in the

z{direction to obtain the elements that apply in this particular direction, i.e.

 Hrz Hzz ! ii = ; N X j= 1 i6=j  Hrz Hzz ! ij  fori= 12N: (21)

To determine the other two coecients Hrr and Hzr, we employ a particular

Stokes ow from which the velocity and tension can be computed for any arbi-trarily shaped region. For example, we use the following axi{symmetric Stokes

ow 

vc = (r=6;z=3)T



bc = (nr0)T:

Substituting this solution into Eqn. (19) we obtain

 Hrr Hzr ! ii =  Hrz Hzz ! ii ;z 3 i ; PN j= 1 i6=j  Hrr Hrz Hzr Hzz ! ij  r 6 ;z 3 ! j +PM j=1  Grr Grz Gzr Gzz ! ij  nr 0 ! j 1 A  6 r  i fori= 12N: (22)

Note, however, that if the point xi is on the z{axis, r becomes zero and

so it is impossible to calculate these diagonal terms from the above equation. Fortunately, there is no need to calculate these diagonal terms when the load

point is on the z{axis, because the radial velocity and stress at the z{axis must

be zero for all axi-symmetric problems. That means that if the load point is on

the z-axis then its diagonal terms in the radial direction have no inuence on

the overall system of linear algebraic equations, and so we can give any non-zero values to them.

The nal thing we need to consider is the application of the boundary con-ditions. If the stress vector is given at a boundary point then nothing has to

be done if the velocity vector is prescribed at boundary point xi then we have

to interchange the i{th column of H and the i{th column of G after reversing

its sign. The mixed boundary condition, by which we mean a combination of the velocity and the stress vector is prescribed at a boundary point, can also be handled in a similar way. After this kind of rearrangement we have

H 

Zunknown =G 

Zknown (23)

whereZunknownandZknownstand for unknown boundary data and known

(9)

2.3 Mesh Redistribution

In this problem, the shape of the mould and the plunger is given by a set of

discrete coordinates in the (rz) plane. Furthermore, after one time step, the free

boundary is only known at a set of boundary nodes. An algorithm is needed to remesh the boundary of the glass domain.

We assume that the boundary ; of the uid region can be parameterized

with respect to the arc length s, i.e.

x(s)2; 0sL:

A grid can be described by

fx(s 0) x(sN ;1) g or equivalently by fs 0 sN ;1 g: The quantity hi =si;si

;1 is called the step-length.

Letxi

;1 = x(si

;1) and

xi=x(si) be two given successive nodal points. The

next node xi

+1 has to lie at a distance

hi

+1 from

xi. Firstly the quasi{uniform

condition(cf. Trevelyan18]) is required, i.e. , for some A>0

hi A

hi +1

Ahi fori= 0N ;2: (24)

Secondly we need a rened mesh on the free boundary, especially near the plunger Note that the upper part of the mould has larger curvatures. To ap-proximate this part of the mould boundary we divide the mould boundary into two parts: the upper part and the lower part. On the upper part smaller

val-ues for h are needed than on the lower part of the mould boundary. Of course

at interfaces{between the mould and the free boundaries, between the plunger and the free boundaries, between the upper and the lower part of the mould boundaries{the condition (24) should be satised. A schematic grid is shown in Figure 3.

3 Time Integration

In this section, we discuss how the geometry of the glass is updated. Suppose

we are at time t

0. From the Stokes problem we can determine the velocity eld

v

0 at time

t

0, which is used to determine the geometry of the glass in cylindrical

coordinates at timet

1 := t

0+

tby discretising the ordinary di erential equation

(6). We rewrite this initial value problem as follows

( dx dt =v(x(t)) x(t 0) = x 0 : (25)

(10)

Figure 3: A schematic grid

We denote the area between the mould and the plunger at time tby At:Of

course the glass body at time tshould be in At for all t. If the initial geometry

of the glass, represented by its boundary pointsx

0 satises x 0 i 2At 0  for each i

the new geometry of the glassx

1 at time

t 1 :=

t 0+

thas to satisfy the constraint

x 1

i 2At 1

 for each i: (26)

Suppose we employ, say the Euler forward scheme, to discretize the problem (25)

y 1 = x 0+ tv 0 :

In general, constraint (26) will not be satised. That means we need a strategy

to reposition those points which are not in At

1: if y

1

i 2 At

1 then no further

action is needed otherwise they are redened to be at the intersection of the

boundary of At

1 and the line which is from

x 0

i to y 1

i (see Figure 4). We refer to

this repositioning step as theclipping algorithm. Applying the clipping algorithm

to y

1 we obtain the new geometry of the glass at time

t 1.

We prefer the mid{point rule as it is symplectic(cf. Stuart 16]). For the

initial value problem (25) the midpoint rule reads

x 1= x 0+ tv(x(t1 2 )) (27) where t1 2 :=t 0+ t

2 . But the information

v(x(t1 2

)) is unknown. To approximate

v(x(t1 2

)), we rst use the Euler forward scheme with time step t

(11)

x y 0 1 1 inside outside x free boundary i i i

Figure 4: The clipping algorithm

x(t1 2 ), i.e. x(t 1 2 ) =x 0+ t 2 v 0  (28) then v(x(t1 2

)) is achieved by solving a Stokes ow with a boundary represented

by x(t1

2

): In summary we have the following stages at each time step

Solve a Stokes problem on the domain represented byx

0 to get the velocity

v

0 at time

t 0.

Use the Euler forward scheme to predict the geometry of the glass x

1 2 at

time t1

2

.

Solve a Stokes problem on the domain represented by x

1 2 to predict the velocity v 1 2 at time t1 2 .

Use the midpoint rule to determine the new geometry of the glass x

1 at

time t

1.

Note that if necessary the clipping algorithm should be used before solving the Stokes problem.

4 Numerical Results

To show the eciency and the accuracy of the method above, we give some examples here.

Some initial shape of the glass drop has to be chosen when the plunger starts

(12)

t=0 (a) t=1.8 (b) t=3.6 (c)

t=4.6 (d) t=4.8 (e) t=5.1 (f)

Figure 5: Evolution of the glass domain

is assumed to have a given velocity the initial position of the plunger is a little bit in the glass drop(see Figure 5a).

The volume of the glass drop is determined by the volume of the parison since the glass is considered incompressible. The evolution of the glass domain is shown in Figure 5.

During the nal stage of the pressing phase, one can see several free bound-aries which are separated by the mould boundary, see also amplied version of Figure 5e, displayed in Figure 6.

The mass should be conserved. The mass as a function of time t is plotted

in Figure 8. In this example only 0:55% of the mass is lost due to the clipping

algorithm.

Intuitively, the velocity eld near the plunger has larger gradients if the no{ slip boundary condition is used. Figure 7 gives an example showing what happens

(13)

free boundaries

Figure 6: The free boundaries

mould plunger

Figure 7: The velocity eld on the free boundary

on the free boundary near the plunger. If a FEM is used then a very rened grid is needed. The BEM has the same problem on the free boundary near the plunger, but less serious. That means the BEM is cheaper although the FEM results in a sparse system. The BEM is more ecient in this situation.

References

1] Becker, A.A., The Boundary Element Method in Engineering,

McGRAW-HILL BOOK COMPANY, London, 1992.

2] Becker, A.A., The Boundary Integral Equation Method in Axisymmetric

Stress Analysis Problems, Springer{Verlag, 1986.

3] Brebbia, C.A., Telles, J.C.F. & Wrobel, L.C.,Boundary Element Techniques,

(14)

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 55 55.5 56 56.5 57 57.5 58 58.5 59 59.5 60

Figure 8: Mass as a function of time

4] Chandra, T.D. & Rienstra, S.W.,Analytical Approximations to the Viscous

Glass Flow Problem in the Mould{Plunger Pressing Process, RANA 97{08,

Eindhoven University of Technology, 1997.

5] Cody, W.J.,Chebyshev Approximation for the Complete Elliptic Integrals K

and E, Math. of Computation, 19, pp.105{112, 1965.

6] Fung, Y.C.,A First Course in Continuum Mechanics, Prentice{Hall, 1969.

7] Ladyzhenskaya, O.A., The Mathematical Theory of Viscous Incompressible

Flow, Gordon and Breach, New York{London, 1963.

8] Lu, Pin & Huang, Mao{Kuang,Some Problems in Boundary Element

Pro-gramming of Axisymmetric Elasticity, Applied Mathematics and Mechanics,

11, No.1, pp.45{52, 1990.

9] Knupp, P. & Steinbeg, S., Fundamentals of Grid Generation, CRC Press,

1993

10] Laevsky, C. & Mattheij, R.M.M., Development of a Simulation Model for

Glass Flow in Bottle and Jar Manufacturing, RANA 98{18, Eindhoven

Uni-versity of Technology, 1998.

11] Mattheij, R.M.M. & Molenaar, J.,Ordinary Dierential Equations in Theory

and Practice, Wiley, 1996.

12] Mattheij, R.M.M., Wang, K. & ter Morsche, H.G.,Modelling Glass Parisons,

in this book, 1999.

13] Partridge, P.W., Brebbia, C.A. & Wrobel, L.C., The Dual

Reci-procity Boundary Element Method, Computational Mechanics Publications,

(15)

14] Pozrikidis, C., Boundary Integral and Singularity Methods for Linearized

Viscous Flow, Cambridge University Press, 1992.

15] Simons, P., The Cooling of Molten Glass in a Mould, Technique Report,

Eindhoven University of Technology, 1996.

16] Stuart, A.M. & Humphries, A.R., Dynamical Systems and Numerical

Anal-ysis, Cambridge University Press, 1996.

17] Telles, J.C.F., A Self{Adaptive Co{ordinate Transformation for Ecient

Numerical Evaluation of General Boundary Element Integrals, International

Journal for Numerical Methods in Engineering,24, pp.959{973, 1987.

18] Trevelyan, J., Boundary Elements for Engineers{Theory and applications,

Computational Mechanics Publications, 1994.

19] van de Vorst, G.A.L. & Mattheij, R.M.M., A BEM{BDF scheme for

Cur-vature Driven Moving Stokes Flows, J. Comput. Phys. No. 1, 1995.

20] van de Vorst, G.A.L., Modelling and Numerical Simulation of Viscous,

Sin-tering, PH.D. Thesis, Eindhoven University of Technology, 1994.

21] van de Vorst, G.A.L. & Mattheij, R.M.M., Numerical analysis of a 2{D

vis-cous sintering problem with non{smooth boundaries, Computing,49, pp.239{

263, 1992.

22] van de Vorst, G.A.L. & Mattheij, R.M.M., A Boundary Element Solution

for Two{dimensional Viscous Sintering, J. Comput. Phys., 100, pp.50{63,

Referenties

GERELATEERDE DOCUMENTEN

3.1 Temperature The bio-char percentage yield is shown in figure 1 Figure 2 shows the effect of temperature on the biochar yield for both feedstocks with nitrogen as atmosphere...

de Wit onderzoeksschool PE&RC, Wageningen • Plant Research International, Wageningen • PPO, Naaldwijk • Wageningen Universiteit, Wageningen In dit rapport wordt een overzicht

Hierdie studie het ten doel om versorgers by 'n fasiliteit vir volwassenes met 'n intellektuele gestremdheid se persepsie van hulle lewensgehalte te verken en te beskryf vanuit

Ook deze kuil leverde enkele opmerkelijke vondsten op, waaronder fragmenten van een gasmasker (type 'small box respirator') en een vrij groot pakket menselijk haar. Het gaat om haar

- Voor waardevolle archeologische vindplaatsen die bedreigd worden door de geplande ruimtelijke ontwikkeling en die niet in situ bewaard kunnen blijven:. Wat is de

In order to move from a position of ‘knowing’ that we were better than the client to a position where we ‘knew’ that we had to work to- gether with the client’s experts required

This research seeks to establish the political role that the City Press defined for its black journalists in post-apartheid South Africa, and the role played by

Structural Health Monitoring of a helicopter tail boom using Lamb waves – Advanced data analysis of results obtained with integrated1. optical fibre