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
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.
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
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)
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 !
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)
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
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
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)
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
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
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
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,
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,
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,