Mathematical modelling of glass forming processes
Citation for published version (APA):Groot, J. A. W. M., Mattheij, R. M. M., & Laevsky, K. Y. (2009). Mathematical modelling of glass forming processes. (CASA-report; Vol. 0907). Technische Universiteit Eindhoven.
Document status and date: Published: 01/01/2009
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
providing details and we will investigate your claim.
Mathematical Modelling of Glass Forming Processes
J. A. W. M. Groot
R. M. M. Mattheij
K. Y. Laevsky
January 29, 2009
Contents
1 Introduction 4 1.1 Glass Forming . . . 4 1.2 Process Simulation . . . 6 1.3 Outline . . . 8 2 Mathematical Model 8 2.1 Geometry, Problem Domains and Boundaries . . . 82.2 Balance Laws . . . 10
2.2.1 Thermodynamics . . . 10
2.2.2 Mechanics . . . 13
3 Parison Press Model 15 3.1 Mathematical Model . . . 16
3.2 Slender-Geometry Approximation . . . 17
3.3 Motion of the Plunger . . . 23
3.4 Simulation model . . . 26 3.5 Results . . . 29 4 Blow Model 32 4.1 Mathematical Model . . . 32 4.2 Glass-Air Interfaces . . . 34 4.3 Variational Formulation . . . 36 4.4 Simulation model . . . 39 4.5 Results . . . 41
5 Direct Press Model 44 5.1 Mathematical Model . . . 44
5.2 Simulation model . . . 46
Abstract
An important process in glass manufacture is the forming of the product. The forming process takes place at high rate, involves extreme temperatures and is characterised by large deformations. The process can be modelled as a coupled thermodynamical/mechanical problem with corresponding interaction between glass, air and equipment. In this paper a general mathematical model for glass forming is derived, which is specified for different forming processes, in particular pressing and blowing. The model should be able to correctly represent the flow of the glass and the energy exchange during the process. Various modelling aspects are discussed for each process, while several key issues, such as the motion of the plunger and the evolution of the glass-air interfaces, are examined thoroughly. Finally, some examples of process simulations for existing simulation tools are provided.
Nomenclature
Br [-] Brinkman number Fr [-] Froude number Nu [-] Nusselt number Pe [-] P´eclet number Re [-] Reynolds number Fe [N] external force on plungerFg [N] force of glass on plunger
L [m] characteristic length
T [K] temperature
Tg [K] glass temperature
Tm [K] mould temperature
V [m s−1] characteristic flow velocity
Vp [m s−1] plunger velocity cp [J kg−1K−1] specific heat g [m s−2] gravitational acceleration p [Pa] pressure rp [m] radius of plunger t [s] time
¯zp [m] vertical plunger position
α [W m−2K−1] heat transfer coefficient
β [N m−3s] friction coefficient λ [W m−1K−1] effective conductivity µ [kg m−1s−1] (dynamic) viscosity ρ [kg m−3] density e [-] unit vector n [-] unit normal t [-] unit tangent u [m s−1] flow velocity uw [m s−1] wall velocity x [m] position I [-] unit tensor ˙
E [s−1] strain rate tensor
1
Introduction
Nowadays glass has a wide range of uses. By nature glass has some special characteristics, including shock-resistance, soundproofing, transparency and reflecting properties, which makes it particularly suitable for a wide range of applications, such as windows, television screens, bottles, drinking glasses, lenses, fibre optic cables, sound barriers and many other applications. It is therefore not surprising that glass manufacture is an extensive branch of industry.
1.1
Glass Forming
The manufacture of a glass product can be subdivided into several processes. Below a common series of glass manufacturing processes are described in the order of their application.
Melting: In industry the vast majority of glass products is manufactured by melting raw materials and recycled glass in tank furnaces at an elevated temperature [35, 61]. Examples of raw materials include silica, boric oxide, phosphoric oxide, soda and lead oxide. The temperature of the molten glass in the furnace usually ranges from 1200 to 1600◦C. A slow formation of the liquid is required
to avoid bubble forming [61].
Forming: The glass melt is cut into uniform gobs, which are gathered in a forming machine. In the form-ing machine the molten gobs are successively forced into the desired shape. The formform-ing technique used depends on the type of product. Forming techniques include pressing, blowing and combina-tions of both, and are discussed further on. During the formation the glass slightly cools down to below 1200◦C. After the formation the glass objects are rapidly cooled down as to take a solid form.
Annealing: Development of stresses during the formation of glass may lead to static fatigue of the product, or even to dimensional changes due to relaxation or optical refraction. The process of reduction and removal of stresses due to relaxation is called annealing [61].
In an annealing process the glass objects are positioned in a so-called Annealing Lehr, where they are reheated to a uniform temperature region, and then again gradually cooled down. The rate of cooling is determined by the allowable final permanent stresses and property variations throughout the glass [61].
Surface treatment: An exterior surface treatment is applied to reduce surface defects. Flaws in the glass surface are removed by chemical etching or polishing. Subsequent flaw formation may be prevented by applying a lubricating coating to the glass surface. Crack growth is prevented by chemical tem-pering (ion exchange strengthening), thermal temtem-pering or formation of a compressive coating. For more information about flaw removal and strengthening of the glass surface the reader is referred to [61].
The process step of interest in this paper is the actual glass formation. Below three widely used forming techniques are discussed. See [24, 61, 73] for further details on glass forming.
Press process: Commercial glass pressing is a continuous process, where relatively flat products (e.g. lenses, TV screens) are manufactured by pressing a gob that comes directly from the melt [61]. This
process is usually referred to as the press process. In this paper also the term direct press is used to distinguish the process from the parison press, which is explained further on.
The direct press is depicted in Fig. 1. Initially, the glass gob is positioned in the centre of a mould. Over the mould a plunger is situated. In order to enclose the space between the mould and the plunger, so that the glass cannot flow out during the process, a ring is positioned on top of the mould. During the direct press the plunger moves down through the ring so that the gob is pressed into the desired shape. ring glass mould plunger F 1: P
Press-blow process: A hollow glass object is formed by inflating a glass preform with pressurised air. This is called the final blow. The preform is also called parison.
In a press-blow process first a preform is constructed by a press stage, to avoid confusion with the direct press here referred to as the parison press. Figure 2 shows a schematic drawing of a press-blow process. In the press stage a glass gob is dropped down into a mould, called the blank mould, and then pressed from below by a plunger (see Fig. 2(a)). Once the gob is inside the blank mould, the baffle (upper part of the mould) closes and the plunger moves gradually up. When the plunger is at its highest position, the ring closes itself around the plunger, so that the mould-plunger construction is closed from below. Finally, when the plunger is lowered, the ring is decoupled from the blank mould and the glass preform is carried by means of a robotic arm to another mould for the blow stage (final blow). In the blow stage the preform is usually first left to sag due to gravity for a short period. Then pressurised air is blown inside to force the glass in a mould shape (see Fig. 2(b). It is important to know the right shape of the preform beforehand for an appropriate distribution of glass over the mould wall.
Blow-blow process: A blow-blow process is based on the same principle as a press-blow process, but here the preform is produced by a blow stage. By means of a blow stage a hollow preform can be formed, which is required for the production of narrow-mouth containers. In practice the glass gob is blown twice to create the preform. Figure 3 shows a schematic drawing of the first blow stage of a blow-blow process. First in the settle blow the glass gob is blown from above to form the neck of the container (see Fig. 3(a)), then in the counter blow from below to form the preform (see
gob
mould
plunger
baffle
(a) parison press
preform (b) final blow F 2: S - ring glass mould
(a) settle blow
preform
(b) counter blow
F 3: S -
Fig. 3(b)). After the counter blow the preform is carried to the mould for the second blow stage, the final blow. The final blow is basically the same as in the press-blow process, but the preform is typically different.
The temperature of the material of the forming machine is typically around 500◦C. Because of the high temperature of the gob, the surface temperature of the material will increase. To keep the temperature of the material within acceptable bounds the mould and plunger are heat insulated by means of water-cooled channels.
1.2
Process Simulation
In the recent past glass forming techniques were still based on empirical knowledge and hand on expe-rience. It was difficult to gain a clear insight into forming processes. Experiments with glass forming
machines were in general rather expensive and time consuming, whereas the majority had to be performed in closed constructions under complicated circumstances, such as high temperatures.
Over the last few decades numerical process simulation models have become increasingly important in understanding, controlling and optimising the process [35, 71]. Since measurements are often complicated, simulation models may offer a good alternative or can be used for comparison of results.
The growing interest of glass industry for computer simulation models has been a motivation for a fair number of publications on this subject. The earliest papers dealt with computer simulation of glass con-tainer blowing [7, 14]. The first publication in which both stages of the press-blow process were modelled was presented in [74], although the model did not include aspects such as the drop of the gob into the blank mould and the transfer of the preform to the mould for the final blow [30].
Shortly afterward also the first PhD theses on the modelling of glass forming processes appeared. A numerical model for blowing was published in [9]. Not much later a simplified mathematical model for pressing was presented in [29]. Both models assumed axial symmetry of the forming process.
Subsequent publications reported the development of more advanced models. A fully three-dimensional model for pressing TV panels was addressed in [37]. A complete model for the three-dimensional simu-lation of TV panel forming and conditioning, including gob forming, pressing cooling and annealing was developed by TNO [49]. A simulation model for the complete press-blow process, from gob forming until the final blow, was presented in [30]. The model also included effects of viscoelasticity and surface ten-sion. Finally, an extensive work on the mathematical modelling of glass manufacturing was published in the book entitled ‘Mathematical Simulation in Glass Technology’ [35].
More recent papers focus on optimisation of glass forming processes [27]. For instance a numerical optimisation method to find the optimum tool geometry in a model for glass pressing was introduced in [44]. Optimisation methods have also been developed to estimate heat transfer coefficients [11] or the initial temperature distribution in glass forming simulation models [43]. An engineering approach to find the optimum preform shape for glass blowing was addressed in [41]. Their algorithm attempts to optimise the geometry of the blank mould in the blow-blow process, given the mould shape at the end of the second blow stage. More recently, an optimisation algorithm for predictive control over a class of rheological forming processes was presented in [4].
By far most papers on modelling glass forming processes use FEM (Finite Element Methods) for the numerical simulation. Exceptions are [75], in which a Finite Volume Method is used, and [16], in which Boundary Element Methods are used. FEM in models for forming processes usually go together with re-meshing techniques, sometimes combined with an Eulerian formulation or a Lagrangian method. In [49] remeshing was completely avoided by using an arbitrary Euler-Lagrangian approach to compute the mesh for the changing computational domain due to the motion of the plunger and by using a Pseudo-Concentration Method to track the glass-air interfaces. In [22] a Level Set Method was used to track the glass-air interfaces in a FEM based model for glass blowing.
In this paper a general mathematical model for the aforementioned forming processes is derived. Sub-sequently, the model is specified for different forming processes, thereby discussing diverse modelling aspects. The paper focusses on the most relevant aspects of the forming process, rather than supplying a model that is as complete as possible. As discussed previously, a considerable amount of work on mod-elling glass forming processes has been done in the recent past. For completeness, comparison or reviewing
various references to different simulation models are included. Finally, some examples of process simula-tions for existing simulation tools are provided. Most of the work in this paper has been done by CASA1 (Centre for Analysis, Scientific computing and Applications).
1.3
Outline
The paper is structured as follows. First in § 2 the physical aspects of glass forming are described and a general mathematical model for glass forming is derived. Then the mathematical model is specified for the parison press in § 3, for both the counter blow and the final blow in § 4 and for the direct press in § 5. In each of these sections also a computer simulation model for the forming process in question is described and some results are presented.
2
Mathematical Model
This section presents a mathematical model for glass forming in general. The section has the following structure. First Section 2.1 defines the physical domains into which a glass forming machine can be subdi-vided and in which the boundary value problems for glass forming are defined. Then Section 2.2 describes the physical aspects of glass forming, sets up the resulting balance laws and derives the corresponding boundary value problems. Subsection 2.2.1 is concerned with the thermodynamics and Subsection 2.2.2 with the mechanics.
2.1
Geometry, Problem Domains and Boundaries
In order to formulate a mathematical model for glass forming, the forming machine is subdivided into subdomains. First the space enclosed by the equipment is subdivided into a glass domain and an air domain. Then separate domains for components of the equipment (e.g. mould, plunger) are considered. Subdomains of the equipment are of interest when modelling the heat exchange between glass, air and equipment. On the other hand, for less advanced heat transfer modelling it can be assumed that the equipment has constant temperature, so that the mathematical model can be restricted to the glass and air domain. In this case the equipment domains are disregarded.
Figure 4 illustrates the domain decomposition of 2D axi-symmetrical forming machines for the direct press, parison press, counter blow and final blow, respectively. The entire open domain of the forming machine, consisting of equipment, glass melt and air, is denoted byΣ. The ‘flow’ domain Ω consists of the open glass domainΩg, the open air domainΩa and the glass-air interface(s)Γi. For blowingΩ :=
Ωg∪Ωa∪Γiis fixed, while for pressingΩ changes in time. Furthermore, ΩgandΩaare variable in time
for any forming process. The boundaries of the domains are:
Γb : baffle boundary Γi : glass-air interface Γm : mould boundary
Γo : outer boundary Γp : plunger boundary Γr : ring boundary
Γs : symmetry axis
Γ Γ Γ Γi r z mould plunger ring glass air Ω Ω Ω Ωg ΩΩΩΩa t n ΓΓΓΓp Γ ΓΓ Γr Γ ΓΓ Γm Γ ΓΓ Γo Γ Γ Γ Γs Γ Γ Γ Γo Γ ΓΓ Γo (a) plunger mould baffle Ωg glass Ωaair Γs Γb Γo Γo Γm Γi Γp Γg,o t n z r (b) Γa,o Γg,o Γs Ωg glass Ωa air Γi Γm baffle Γb Γo Γo n t Γi mould r z (c) n t Ωa air Ωa air Ωg glass Γo Γm Γi Γi Γs Γo,a Γo,g mould r z (d) F 4: 2D -
Note that not necessarily all boundaries exist for each forming machine. DomainΣ is enclosed by Γo∪Γs.
DomainΩ is enclosed by ∂Ω := Ω ∩Γe∪Γo∪Γs
, with
Γe:= Γb∪Γm∪Γp∪Γr (equipment boundary). (2.1)
In addition, define ∂Ωa:= Ωa∩∂Ω and ∂Ωg:= Ωg∩∂Ω. Finally, the boundaries for the glass domain and
the air domain are distinguished:
Remark 2.1 In pressing sometimes only the forces acting on the glass are of interest. If the density and viscosity of air are negligible compared to the density and the viscosity of glass, the mathematical model may be restricted to the glass domainΩg. In this case the interfaceΓiis replaced byΓo.
2.2
Balance Laws
2.2.1 Thermodynamics
Glass forming involves high temperatures within a typical range of 800◦C-1400◦C. Temperature variations
within this range may cause significant changes in the mechanical properties of the glass.
• The range of viscosity for varying temperature is relatively large for glass: it amounts from 10 Pa s at the melting temperature (about 1500◦C) to 1020Pa s at room temperature. The viscosity increases rapidly as a glass melt is cooled, so that the glass will retain its shape after the forming process. Typical values for the viscosity in glass forming processes lie between 102and 105Pa s [61, 73]. The temperature dependence for the viscosity of glass within the forming temperature range is given by the VFT-relation, due to Vogel, Fulcher and Tamman [61, 73]:
µ(T )\[Pa s] = exp − A + B/(T − TL). (2.2)
Quantities A [-], B [◦C] and T
L[◦C] represent the Lakatos coefficients, which depend on the
com-position of the glass melt. Figure 5 shows how strongly the viscosity depends on temperature for soda-lime-silica glass with Lakatos coefficients A = 3.551, B = 8575◦C and TL= 259◦C [52].
800 850 900 950 1000 1050 1100 1150 1200 2 2.5 3 3.5 4 4.5 5 5.5 Temperature (oC) log 10 viscosity F 5: VFT-
• The following density-temperature relation can be deduced
ρ(T ) = ρ0
1 − αV T − Tref, (2.3)
– αV [◦C−1] is the volumetric thermal expansion coefficient,
– Tref[◦C] is a reference temperature,
– ρ0[kg m−3] is the density at the reference temperature.
The volumetric thermal expansion coefficient is often assumed constant; for molten glass it is typ-ically ranged from 5
·
10−5 to 8·
10−5 ◦C−1. The density of molten glass is of the order 2300 to2500 kg m−3and is 5 to 8% lower than at room temperature. Thus it is quite reasonable to assume incompressibility for molten glass [73].
Clearly, the mechanics of glass forming is related to the heat transfer in the glass. The heat transfer in glass, air and equipment is described by the heat equation for incompressible continua:
ρcp ∂T ∂t +u
·
∇T | {z } advection = ∇·
q |{z} conduction radiation + 2µ ˙E:
∇ ⊗u | {z } dissipation , in ◦Σ
, (2.4)where the temperature distribution T [K] is unknown. Here Σ denotes domain Σ minus the interfaces◦ between the continua. On the interfaces between different continua a steady state temperature transition is imposed,
hhλn
·
∇Tii = 0. (2.5)The heat flux q [W m−2] is the result of the contribution of both conduction and radiation,
q = −λ
·
∇T, (2.6)where λ is the effective conductivity [W m−1K−1], given by
λ = λc+ λr. (2.7)
Here λcis the thermal conductivity and λris the radiative conductivity. The thermal conductivity measures
1.0 W m−1K−1at room temperature for soda-lime glass and increases with approximately 0.1 W m−1K−1
per 100 K. In this paper it is assumed to be constant. The calculation of the radiative conductivity λr is
often a complicated process. However, for non-transparent glasses it can be simplified by the Rosseland approximation [18, 54, 73] λr(T )= 16 3 n2σT3 α , (2.8) where
• σ is the Stefan Boltzmann radiation constant [W m−2K−4],
• n is the average refractive index [-],
The radiative conductivity λrin the sense of (2.8) is called the Rosseland parameter. Relation (2.8) cannot
be applied for highly transparent glasses, since in this case not all radiation is absorbed by the glass melt [73]. A more simple approach is to omit the radiative term, which is often reasonable for clear glass [36, 39]. For more information on heat transfer in glass by radiation the reader is referred to [35, 39, 40, 45, 73]. The specific heat cpis in general slightly temperature dependent. In [3] an increase in the specific heat of
less than ten percent in a temperature range of 900 K to 1300 K for soda-lime-silica glass is reported. For various specific heat capacity models the reader is referred to [3, 60].
In order to analyse the energy exchange problem quantitatively the heat equation is written in dimen-sionless form. Define a typical: velocity V, length scale L, viscosity ¯µ, specific heat ¯cp, effective
conduc-tivity ¯λ, glass temperature Tgand mould temperature Tm. Then introduce the dimensionless variables
t∗:=Vt L, x ∗ := x L, u ∗ := u V, T ∗ := T − Tm Tg− Tm , c∗p:= cp ¯cp , λ∗ := λ ¯ λ, µ ∗= µ ¯µ,(2.9) For convenience all dimensionless variables, spaces and operators with respect to the dimensionless vari-ables are denoted with superscript∗. Substitution of the dimensionless variables (2.9) in the heat equation
and splitting up the effective conductivity into (2.7) lead to the dimensionless form, Pe ∂T ∗ ∂t∗ + u ∗
·
∇∗ T∗= ∇∗·
λ∗∇∗T∗ + 2µ∗Br ˙E∗:
∇∗⊗u∗, in ◦Σ
∗, (2.10) where Pe= ρ¯cpV L ¯ λ , (2.11) Br= ¯µV 2 ¯ λ Tg− Tm (2.12)are the P´eclet number and the Brinkman number, respectively. The P´eclet number represents the ratio of the advection rate to the diffusion rate. On the other hand, the Brinkman number relates the dissipation rate to the conduction rate. The dimensionless numbers are useful for assessing the order of magnitude of the different terms in (2.10).
The energy BCs follow from symmetry and heat exchange with the surroundings: λ∇T
·
n = 0, on Γsλ∇T
·
n = α T − T∞, on Γo,where T∞is the temperature of the surroundings. The heat transfer coefficient α [W m−2K−1] can differ
for separate equipment domains, such as the mould and the plunger. Let ¯α be a typical value for the heat transfer coefficient, then the Nusselt number is defined by
Nu=αL¯ ¯
λ . (2.13)
The dimensionless boundary conditions become λ∗∇T∗
·
n= 0, on Γs
λ∗∇T∗
·
n= Nu α∗ T∗− T∗∞, on Γo,
(2.14)
with α∗ = α/ ¯α. On the other hand, if the heat transfer in the equipment domain is not of interest, the
2.2.2 Mechanics
A balance law for the mechanics of the glass melt is formulated. In general, glass can be treated as an isotropic viscoelastic Maxwell material [3, 8, 9, 61], that is the strain rate tensor can be split up into an elastic and a viscous part:
˙
E= ˙Ee+ ˙Ev, (2.15)
where the elastic and viscous strain rate tensors, ˙Eeand ˙Evrespectively, are given by [9]
˙ Ee= 1 − 2ν E αV ∂T ∂t + 1 3tr( ˙T ) +1+ ν E dev( ˙T ) (2.16) ˙ Ev= 1 2µdev(T ). (2.17)
Here αV [◦C−1] is the volumetric thermal expansion coefficient, E [Pa] is the Young’s modulus, ν [-] is
the Poisson’s ratio and ˙T [Pa s−1] denotes the stress rate tensor. However, at relatively low viscosities the
relation between shear stress and viscosity becomes approximately linear. For example, for soda lime silica glasses the viscosity as a function of the strain rate and the temperature becomes [6, 63]
µ( ˙E, T ) = µ0(T ) 1+ 3.5 · 10−6E˙µ0.76
0 (T )
, (2.18)
where µ0 is the Newtonian viscosity. Consequently, the motion of glass is dominated by viscous flow
and the influence of elastic effects can be neglected [3, 8, 61]. Moreover, as verified in § 2.2.1, glass is practically incompressible in the forming temperature range, from which it follows that
tr( ˙E)= 0. (2.19)
It can be concluded that glass in the forming temperature range behaves as an incompressible Newtonian fluid [3, 8, 73], i.e. T = −pI + 2µ ˙E, (2.20) with ˙ E= 1 2 ∇ ⊗u+ (∇ ⊗ u) T. (2.21)
For simplicity the pressurised air is considered as an incompressible, viscous fluid with uniform viscos-ity. Thus the motion of glass melt and pressurised air is described by the Navier-Stokes equations for incompressible fluids. These involve the momentum equations,
ρ∂u ∂t +u
·
∇ ⊗u | {z } inertia = −∇p |{z} pressure + ρg |{z} gravity + 2∇·
µ ˙E | {z } viscosity , in Ω \ Γi, (2.22)and the continuity equation, which follows directly from (2.19),
The unknowns are the flow velocity u [m s−1] and the pressure p [Pa].
Flow problem (2.22)-(2.23) is coupled to the energy problem (2.10) in two ways: firstly the viscosity is temperature dependent and secondly the heat transfer is partly described by convection and diffusion.
In order to apply a quantitative analysis the Navier-Stokes equations are written in dimensionless form. First a dimensionless pressure is defined by
p∗:= Lp
¯µV. (2.24)
The gravity force can be written as −ρgez, where ez is the unit vector in z-direction. Substitution of the
dimensionless variables (2.9) and (2.24) into the Navier-Stokes equations (2.22)-(2.23) and division by the order of magnitude of the diffusion term, ¯µVL2, lead to the dimensionless Navier-Stokes equations
Re ∂u∗ ∂t∗ + u ∗
·
∇∗⊗u∗ = −∇∗p∗−Re Frez+ 2∇ ∗·
µ∗E˙∗, in Ω∗\Γ∗ i, ∇∗·
u∗= 0, in Ω∗\Γ∗i, (2.25) where Re= ρVL ¯µ , (2.26) Fr= V 2 gL (2.27)are the Reynolds number and the Froude number, respectively. The Reynolds number measures the ra-tio of inertial forces to viscous forces, while the froude number measures the rara-tio of inertial forces to gravitational forces.
The jump conditions between two immiscible viscous fluids are the continuity of the flow velocity,
[[u]]= 0, (2.28)
as well as the continuity of its tangential derivative,
[[ ∇ ⊗ u
·
t]]= 0, (2.29)
and a dynamic jump condition stating the balance of stress across the fluid interface [33, 47],
[[T n]]= 0. (2.30)
Remark 2.2 If the influence of surface tension is taken into account, the dynamic jump condition is
[[T n]]= −γκn, (2.31)
where n points in the air domain. Hereγ[N m−1] denotes the surface tension and κ[m−1] denotes the curvature. The influence of surface tension highly depends on the glass composition [61]. In this paper the influence of surface tension is simply neglected.
Boundary conditions for the flow problem can be determined as follows. OnΓssymmetry conditions
are imposed. OnΓg,ea suitable slip condition for the glass should be adopted. The air can escape through small cavities in (part of) the mould wall. This aspect can be modelled by allowing air to flow freely through the mould wall. Thus free-stress conditions are proposed on this part of the equipment boundary, which is referred to asΓ(1)a,e. A free-slip condition is prescribed on the remaining part,Γ(2)a,e. OnΓothe normal
stress should be equal to the external pressure.
A commonly used boundary condition to describe fluid flow at an impenetrable wall [15, 21, 31, 53] is Navier’s slip condition:
T n+ β(u − uw)
·
t= 0, (2.32)where β is the friction coefficient [N m−3s] and uwis the velocity of the wall [m s−1]. A similar condition
can be obtained by using the Tresca model [18]. Introduce a dimensionless friction coefficient, β∗:= Lβ
¯µ , (2.33)
then the dimensionless Navier’s slip condition reads:
T∗n+ β∗(u∗−u∗w)
·
t= 0. (2.34)The order of magnitude of the friction coefficient depends on many parameters, such as the type of glass, temperature, pressure or presence of a lubricant [17, 20, 53].
Remark 2.3 Forβ∗→ ∞ Navier’s slip condition together with the boundary condition for an impenetrable wall,(u∗−u∗
w)
·
n= 0, can be reformulated as a no slip condition:u∗= u∗w. (2.35)
In summary, the boundary conditions for the flow problem can be formulated as: u∗
·
n= 0, T∗n·
t= 0, on Γs, (u∗−u∗ w)·
n= 0, T ∗n+ β∗(u∗−u∗ w)·
t= 0, on Γ g,e, T∗n·
n= 0, T∗n·
t= 0, on Γ(1)a,e, (u∗−u∗ w)·
n= 0, T ∗n·
t= 0, on Γ(2) a,e, T∗n·
n= p0, T∗n·
t= 0, on Γo, (2.36)where p0is the external pressure.
3
Parison Press Model
This section presents a mathematical model for the parison press. Section 3.1 specifies the mathematical model described in § 2 for the parison press. By restricting the analysis to a narrow channel between the plunger and the mould, an analytical approximation of the flow can be derived. Section 3.2 explains this concept, known as the slender geometry approximation [53]. Section 3.3 describes the motion of the plunger by an ordinary differential equation. It appears that the motion of the plunger is coupled to the flow, which considerably complicates the parison press model. Section 3.4 presents a numerical simulation model for the parison press. The motion of the free boundaries is emphasised. Finally, Section 3.5 shows some examples of parison press simulations. The simulation tool used for these results is presented in [36].
3.1
Mathematical Model
In Section 2.2 the balance laws for glass forming were formulated. In this section these balance laws are further specified for the parison press process. Typical values for the parison press are:
glass density : ρg = 2.5 · 103kg m−3
glass viscosity : ¯µg = 104kg m−1s−1
gravitational acceleration : g = 9.8 m s−2 flow velocity : V = 10−1m s−1
length scale of the parison : L = 10−2m glass temperature : Tg = 1000◦C
mould temperature : Tm = 500◦C
specific heat of glass : ¯cp = 1.5 · 103J kg−1K−1
effective conductivity of glass : λ¯ = 5 W m−1K−1,
As a result the following dimensionless numbers are found:
Peglass≈ 7.6 · 102, Brglass≈ 4.0 · 10−3, Reglass≈ 2.5 · 10−4, Fr ≈ 1.0 · 10−1. (3.1)
Apparently, the heat transport in glass is dominated by thermal advection. As a result the heat equa-tion (2.10) simplifies to dT∗ dt∗ = ∂T∗ ∂t∗ + u ∗
·
∇∗T∗= 0, inΩ∗ g, (3.2)Thus the temperature remains constant along streamlines. Consequently, if the initial temperature distribu-tion in the glass is (approximately) uniform, the glass viscosity can be considered constant. From the small Reynolds number for glass it can be concluded that the inertia forces can be neglected with respect to the viscous forces. Furthermore,
Reglass
Fr ≈ 2.5 · 10
−3,
which means that also the contribution of gravitational forces is rather small. In conclusion, the glass flow can be described by the Stokes flow equations:
∇∗
·
T∗ = 0, ∇∗·
u∗ = 0, inΩ∗g, (3.3)where T∗is the dimensionless stress tensor, which satisfies (2.20) in terms of the dimensionless variables. The air domain is ignored for the following reasons. Firstly, the force of the plunger acts directly on the glass domain. Secondly, the density and viscosity of air are negligible compared to those of glass, so that air hardly forms any obstacle for the glass flow. Thirdly, the simplification of the heat equation to convection equation (3.2) gives reason to restrict the energy exchange problem to the glass domain, or all together ignore the heat transfer in case of an uniform initial glass temperature. Subsequently, the glass-air interfaces are treated as an outer boundary of the flow domain, Γo, on which free-stress conditions are
imposed. Note that although the problem can be considered as a free-boundary problem, the geometry is constrained by the mould and the plunger.
Remark 3.1 Close to the equipment wall extreme temperature variations occur over a small length scale. Therefore, the conductive heat flux close to the equipment wall should, strictly speaking, not be disregarded [53, 64]. Moreover, the viscosity in this region may increase by several orders, so that the fluid friction may not be negligible and the influence of heat generation by dissipation should be taken into account for optimal accuracy [53]. Although the reader should take notice of these boundary layer effects, in this paper it is simply assumed that they are small enough to be ignored. For more advanced heat modelling during the parison press the reader is referred to [25, 40, 64].
The boundary conditions for the flow problem (2.36) can be specified for the parison press:
u∗
·
n= 0, T∗n·
t= 0, on Γ s, (u∗− Vpez)·
n= 0, T∗n+ β∗p(u ∗− V pez)·
t= 0, on Γg,p, u∗·
n= 0, T∗n+ β∗ mu∗·
t= 0, on Γ g,m, T∗n·
n= 0, T∗n·
t= 0, on Γo. (3.4)In the remainder of this section it is assumed that the glass gob initially has an uniform temperature distribution, so that with (3.2) it follows that µ= ¯µ is constant.
3.2
Slender-Geometry Approximation
In the model for glass pressing the analysis can be restricted to the flow in a narrow channel between plunger and mould. In other words, the analysis focusses on the flow in the slender geometry around the plunger, while the flow between the plunger top and the baffle (Fig. 4(b)) is considered practically stagnant. For a more complete analysis the reader is referred to [53]. In the slender-geometry approximation of the flow two typical length scales ` and L are considered, with ` L, where ` is the length scale for the width of the channel and L is the length scale for the length of the channel. Thus variations in r-direction are scaled by ` and variations in z-direction are scaled by L. By means of this scaling the following dimensionless variables can be defined:
t∗:= Vt L, ε := ` L, r ∗ := r εL, z ∗ := z L, u ∗ r:= ur εV, u ∗ z := uz V, p ∗ := ε 2Lp µV . (3.5) In the remainder of this section all variables, spaces and operators are dimensionless and the superscript∗ is ignored. Substitution of the dimensionless variables (3.5) into the Navier-Stokes equations (2.22)-(2.23) leads to the dimensionless 2D axi-symmetrical Navier-Stokes equations
ε3Re `∂u∂t +r ur ∂ur ∂r +uz ∂ur ∂z = − ∂p ∂r +ε 2∂ ∂r 1 r ∂ ∂r rur + ε4 ∂2u r ∂z2 , inΩg, εRe`∂u∂t +z ur ∂uz ∂r +uz ∂uz ∂z = − ∂p ∂z + ∂ ∂r 1 r ∂ ∂r ruz + ε2 ∂2u z ∂z2 + ε Re` Fr , in Ωg, (3.6) 1 r ∂ ∂r rur+ ∂uz ∂z =0, inΩg,
where Re`is the Reynolds number with respect to length scale `, i.e.
Re`= ρV`
Thus, for small ε, (3.6) can be simplified to ∂p ∂r =O(ε2), inΩg, ∂p ∂z = ∂ ∂r 1 r ∂ ∂r ruz + O(ε2), inΩg, (3.8) 1 r ∂ ∂r rur+ ∂uz ∂z =0, inΩg.
If the O(ε2) terms are neglected, system of
equa-tions (3.8) can be recognised as Reynolds’ 2D axial-symmetrical lubrication flow [21, 53].
The equipment boundary conditions in the slender-geometry approximation can be simplified accordingly [53]. Define the plunger surface and the mould surface by:
r= rp z −¯zp(t), r= rm(z), (3.9)
respectively, where z= ¯zp(t) is the top of the plunger
at time t (see Fig. 6). Consider the plunger position ¯zp ≡ ¯zp(t) at a fixed time t. Then the outward unit
normal npand the counterclockwise unit tangent tp
on the plunger wall are given by
glass plunger mould rp rm zp r z
Γ
pΓ
oΓ
sΓ
m F 6: G np = εr0 pez−er q 1+ ε2r02 p , tp= −ez−εr0per q 1+ ε2r02 p . (3.10)Analogously, the outward unit normal nmand the counterclockwise unit tangent tmon the plunger wall on
the mould wall are given by
nm= −εr0 mez+ er p 1+ ε2r02 m , tm= ez+ εr0mer p 1+ ε2r02 m . (3.11)
By substituting (3.10)-(3.11) into Navier’s slip condition (2.32) and scaling the friction coefficient β by a factor ε the following boundary conditions are obtained
βp uz− Vp+ ε2urr0p = 1 − ε2r02 p ∂uz ∂r + ε2 ∂u∂zr − 2ε2r0 p ∂uz ∂z −∂u∂rr q 1+ ε2r02 p , on r= rp(z − ¯zp) (3.12) −βm uz+ ε2urr0m = 1 − ε2r02 m ∂uz ∂r + ε2 ∂u∂zr − 2ε2r0 m ∂uz ∂z −∂u∂rr p 1+ ε2r02 m , on r= rm(z), (3.13)
where Vpis the velocity of the plunger and βp, βmare the friction coefficients corresponding to the plunger
be written as βp uz− Vp= ∂uz ∂r +O(ε2), on r= rp(z − ¯zp), (3.14) −βmuz = ∂uz ∂r +O(ε2), on r= rm(z). (3.15)
The componentwise boundary condition for the impenetrable wall is
ur= (uz− Vp)rp0, on r= rp(z − ¯zp), (3.16)
ur= uzrm0, on r= rm(z). (3.17)
Note that for βm,p→ ∞ the error in boundary condition (3.14)-(3.15) due to the slender-geometry
approxi-mation vanishes.
Following [42, 53] the analytical solution to system of equations (3.8) with set of boundary conditions (3.14)-(3.17) on the equipment boundary and free-stress conditions on the other boundaries can be obtained. Neglecting the O(ε2) terms system of equations (3.8) becomes
∂p ∂r =0, inΩg, (3.18) ∂p ∂z = ∂ ∂r 1 r ∂ ∂r ruz, in Ωg, (3.19) 1 r ∂ ∂r rur+ ∂uz ∂z =0, inΩg. (3.20)
Firstly, from (3.18) it follows that p is a function of z only. Secondly, from (3.19) and the boundary conditions it follows that
uz(r, z)= 1 4 dp dz(z) r 2+ψm(z)χp(r, z) − ψp(z)χm(r, z) χm(r, z) − χp(r, z) ! + Vp χm(r, z) χm(r, z) − χp(r, z) , (3.21) with χp(r, z)= log rrp(z − ¯zp) −βprp(z−¯z1 p), ψp(z)= rp(z − ¯zp) rp(z − ¯zp) −β2p , χm(r, z)= log rrm(z)+βmr1m(z), ψm(z)= rm(z) rm(z)+β2m . (3.22)
Thirdly, from (3.20) and the boundary conditions it follows that
ur(r, z)= 1 r rmrm 0 uz(rm, z) + Z rm r s∂uz ∂z(s, z)ds ! =1 r d dz Z rm r s uz(s, z)ds. (3.23)
Finally, to find the pressure gradient, Gauss’ divergence theorem is applied to the continuity equation,
0= Z Ωg ∇
·
udΩ = Z Γg,o u·
ndΓ + Z Γp u·
ndΓ + Z Ωg∩Γs u·
ndΓ + Z Γm u·
ndΓ. (3.24)The fluxes through the symmetry axis and the mould wall are zero. Since the continuity equation also holds in the plunger domainΣp, it holds that
Z Γp u
·
ndΓ = − Z Σp∩Γo u·
ndΓ = πVpr2p. (3.25)As a result, Z Γg,o u
·
ndΓ = 2π Z rm rp ruzdr= −πVpr2p. (3.26)Substitution of (3.21) into (3.26) yields
r 2 1 8 dp dz r2−ψm(2χp+ 1) − ψp(2χm+ 1) χm−χp + Vp χm+12 χm−χp ! r=rm r=rp = −Vpr2p. (3.27)
By solving (3.27) the following solution for the pressure gradient is found: dp dz(z)= 4Vp ψp(z) − ψm(z) ψp(z) − ψm(z) 2 − χp−χm(z) ωp(z) − ωm(z) , (3.28) with ωp(z)= ψp(z)2− 4rp(z − ¯zp)2 β2 p , ωm(z)= ψm(z)2− 4rm(z)2 β2 m . (3.29)
Figure 7-8 plot the solution at ¯zp = 0.05m and ¯zp = 0.01m, respectively, for ε = 0.1 and Vp = −0.1ms−1.
The (dimensionless) geometries of the plunger and mould have been taken from [53]:
rp(z)= −0.1
√
5z, rm(z)= −0.1
√ 5z.
For simplicity no-slip boundary conditions at the mould and plunger wall have been used. In [53] it is reported that the results are in good agreement with the numerical solution obtained by using FEM. For more results the reader is referred to [42, 53].
0.005 0.01 0.015 0.02 0.06 0.08 0.1 0.12 0.14 r [m] z [m] (a) 0.046 0.06 0.08 0.1 0.12 0.14 0.16 7 8 9 10 11x 10 4 z [m] −dp/dz [Pa/m] (b) [m] (c)r- [m/s] (d)z- [m/s] F 7: A - ¯zp= 0.5 mould plunger
0.005 0.01 0.015 0.02 0.02 0.04 0.06 0.08 0.1 r [m] z [m] (a) 0 0.02 0.04 0.06 0.08 0.1 0.12 0 1 2 3 4 5 6x 10 5 z [m] −dp/dz [Pa/m] (b) [m] (c)r- [m/s] (d)z- [m/s] F 8: A - ¯zp= 0.1 mould plunger
3.3
Motion of the Plunger
In the previous section it was assumed that the plunger moves with a constant flow velocity. However, in practice the plunger is pushed by a piston. This means that the flow velocity is the result of an external force applied to the plunger. In this section it can be seen that the plunger velocity is coupled to the glass flow, which considerably complicates the parison press model.
The press process is initiated by applying an external force Feto the plunger. This causes the plunger
to move with velocity Vp(t). This plunger velocity is the result of the total force F on the plunger, which is
the sum of the external force Feand the force of the glass on the plunger Fg:
dVp dt (t)= F(t) mp = Fe+ Fg(t) mp , (3.30)
where mpis the mass of the plunger [36, 53]. The force Fgis determined by the mechanical forces of the
glass acting on the plunger and hence depends on the glass flow. The glass flow at its turn is caused by the plunger motion. Therefore, differential equation (3.30) is fully coupled to the flow problem. Clearly, if a constant external force is applied to the plunger, the plunger moves until the force of the glass on the plunger is equal in magnitude to the external force.
Below the equation for the motion of the plunger (3.30) is examined more thoroughly. For a complete analysis the reader is referred to [36, 53]. First the force of the glass on the plunger is analysed. At every time t the force can be fully described by the stress tensor (2.20) over the plunger surfaceΓp,
Fg=
Z
Γp
T n
·
ezdΓ. (3.31)
Using the definition of the plunger surface (3.9) the surface element dΓ can be written as dΓ = 2π
q 1+ r02
prpdz, on Γp. (3.32)
By means of the definition of the stress tensor (2.20) and the expressions for the unit normal and tangent (3.10) and the surface element (3.32), the expression for the force of the glass on the plunger becomes
Fg= 2π Z z0 ¯zp p −2µ∂uz ∂z r0p+ µ∂ur ∂z + ∂uz ∂r rpdz. (3.33)
where z = z0(t) is the bottom glass level at time t. Substitution of the dimensionless variables (3.5) into
(3.33) yields Fg= 2πµVL Z z∗ 0 ¯z∗ p p∗− 2ε2∂u ∗ z ∂z∗ r∗p0+ε2∂ur ∂z∗ + ∂u∗ z ∂r∗ r∗pdz∗. (3.34)
Considering expression (3.34) it makes sense to define the dimensionless force of the glass on the plunger as
Fg∗:= Fg
For small ε the dimensionless force can be simplified to Fg∗=: Z z∗ 0 ¯z∗ p p∗r∗p0+∂u ∗ z ∂r∗ r∗pdz∗+ O(ε2). (3.36)
The flow velocity and pressure in (3.36) are implicit functionals of the plunger velocity Vp. However,
the slender-geometry approximation in § 3.2 can be used to find approximate solutions for u∗
z and p∗. By
substituting (3.21) and (3.29) into (3.34)-(3.35) the dimensionless force can be written as
Fg∗(t∗)= Vp∗(t∗)I∗(t∗), (3.37)
with an assumed error of O(ε2), where the dimensionless function I∗(t∗) only depends on the geometry and
the friction coefficient β∗[36, 53]. Then the equation for the motion of the plunger (3.30) in dimensionless
form becomes dV∗ p dt∗(t ∗)=2πµL 2 mpV V∗p(t∗)I∗(t∗)+ F∗e+ O(ε2), (3.38)
where the dimensionless external force F∗eis defined in the same way as Fg∗in (3.35). Typical values in the slender-geometry approximation are:
viscosity : ¯µg = 104kg m−1s−1
flow velocity : V = 10−1m s−1 length scale of the parison : L = 10−1m
mass of the plunger : mp = 1 kg
As a result the dimensionless coefficient in (3.38) is typically 2πµL2
mpV
≈ 104.
This value is rather large, which indicates that (3.38) is a stiffness equation. According to [36] this phe-nomenon can also be observed if (3.38) is solved numerically using the Euler forward scheme for time integration. This means that one would have to resort to an implicit time integration scheme in order to solve (3.38). Unfortunately, a fully implicit scheme is practically impossible, since the plunger velocity is not known explicitly. See [36] for further details on this stiffness phenomenon.
The stiffness phenomenon previously described indicates that the coupling of the equation for the plunger motion (3.30) to the boundary conditions on the plunger for the Stokes flow problem is unde-sirable. Therefore, in the following the plunger velocity Vp(t) is decoupled from the parameter Vpin the
boundary conditions for the flow problem. For a more detailed analysis the reader is referred to [36]. The following lemma is used:
Lemma 3.1 Let (uυ, pυ) be the family of solutions of Stokes flow problem (3.3)-(3.4) with plunger velocity Vp= υ. Define (u, p)υ:= (uυ, pυ). Then
u, pk
1υ1+k2υ2= k1uυ1+ k2uυ2, p0+ k1(pυ1− p0)+ k2(pυ2− p0)
The lemma can be proven by direct substitution of the solution k1uυ1+k2uυ2, p0+k1(pυ1− p0)+k2(pυ2− p0)
into the Stokes flow problem with plunger velocity Vp= k1υ1+ k2υ2[36]. In the remainder of this section
all variables, spaces and operators are dimensionless and the superscript∗is ignored. From lemma 3.1 it immediately follows that
u, pυ= υu1, p0+ υ(p1− p0), (3.39)
and as a result,
Fg,υ= Fg,0+ υ Fg,1− Fg,0. (3.40)
Substitution of (3.40) with υ= Vp(t) into the dimensionless form of (3.30) yields
dVp dt (t)= 2πµL2 mpV Vp(t) Fg,1(t)+ Fg,0(t) + Fg,0(t)+ Fe . (3.41)
The force Fg,1(t) can be calculated by solving the Stokes flow problem in Ωg(t) with plunger velocity
Vp= 1. Thus the equation for the motion of the plunger is decoupled from the boundary conditions for the
Stokes flow problem. Note that the Stokes flow problem is still coupled to differential equation (3.41) as the geometry of the glass domain is determined by the plunger velocity.
The time dependency of the plunger velocity seems a bit awkward as the Stokes flow problem is not explicitly time dependent; the flow merely changes in time through the changing geometry. It seems more convenient to define the force of the glass on the plunger and the plunger velocity as functions of the plunger position z := ¯zp:
Fg:= Fg(z), Vp:= Vp(z). (3.42)
As a result also the motion of the plunger can be described by the plunger position. By the chain rule of differentiation,
dVp
dt (t)= dVp
dz (z)Vp(z). (3.43)
As for the initial condition, let z = 0 and Vp = V0 at t = 0. Then the motion of the plunger involves the
following initial value problem:
1 2 dV2p dz (z)= 2πµL2 mpV Vp(z) Fg,1(z)+ Fg,0(z) + Fg,0(z)+ Fe Vp(0)= V0 (3.44)
Since the plunger velocity does not need to be known to determine the glass domain at given plunger position z, an implicit time integration scheme can be used to solve (3.44), thus overcoming the stiffness problem [36].
3.4
Simulation model
A simulation tool for the parison press process was designed [36]. The tool is able to compute and visualise the velocity field, the pressure and also the temperature in the parison during the press stage. The input parameters for the simulations include the 2D axi-symmetrical parison geometry, i.e. the description of the mould and plunger surfaces, the initial positions of the plunger and the glass domain, as well as the physical properties of the glass. The geometry of the initial computational domain is defined by that of the glass gob (see § 2.1). At the beginning of each time step the computational domain is discretised by means of FEM. A typical mesh for the initial glass domain in a press simulation of the parison for a jar is depicted in Figure 9. The mesh distribution depends on the geometries of the mould and the plunger. For example, the mesh in the ring domain requires a relatively small scale in comparison with the mould domain and the plunger domain.
F 9: M [36].
The implicit Euler method is used to solve the initial value problem for the motion of the plunger (3.44): 1 2 Vkp+1 2 − Vk p 2 zk+1− zk = 2πµL2 mpV Vkp+1Fg,1(zk+1)+ Fg,0(zk+1) + Fg,0(zk+1)+ Fe V0p = V0 (3.45)
The plunger position and hence the geometry of the glass domain can be updated using the approximation
zk+1= zk+ ∆tkVp(zk), tk+1= tk+ ∆tk, (3.46)
where∆tkdenotes the kthtime step.
The motion of the free boundariesΓois described by the ordinary differential equation
dx
Let xi, i = 1, . . . , N be the nodes on a free boundary Γf ⊂Γoand let xki be an approximation of x(tk) (see
Fig. 10). Then the new position of the ithnode can be obtained by the explicit scheme
xki+1= xki + ∆t uki, (3.48)
with uk i := u(x
k
i). A particular question is how to deal with this moving boundary. Depending on the
velocity a situation may be encountered where the obtained position xki+1lies outsideΩg(see Fig. 10). In
t @ @thhht ((( t" " t P P P P PP t !! !! !! t A A A U t Γp Γm Γf Ωg(tk) H H H H H H Y Ωg(˜tk+1) ixki i xk+1 i vk i t @ @ @ R ˜xk+1 i XX XXXt!! !!! t F 10: T [36].
this situation one of the strategies described below may be used. For simplicity only explicit integration is considered.
One approach to deal with the moving boundary is to decrease the time step:
˜xki+1= xki + αik∆t uki, αki ∈ (0, 1], (3.49) where ˜xk+1
i is the new position of the i
thnode obtained with the decreased time step. The kthtime step can be
defined by∆˜tki := miniαki∆t, such that the nodes are situated inside the glass domain at time ˜tk+1:= ˜tk+∆˜tki,
as depicted in Fig. 10. Unfortunately, this algorithm introduces a variable time step that turns out to be too irregular in practice and can be excessively small. In order to have consistency in the topology of the computational domain the time step should be constant [36].
An alternative is illustrated in Fig. 11; the so-called clip algorithm leads the nodes along a discrete ‘solution curve’ until they reach the boundary, thereby clipping the trajectories on the boundary. Thus node i, which would originally leave the physical domain at time tk+1, is clipped on the boundary by (3.49), but with αk
i = 1 for a non-clipped node; the time step only differs for clipped nodes. Regrettably, also this
algorithm has an evident drawback: the clipping influences the mass conservation property of the glass domain [36].
The clip algorithm can be modified to enforce better mass conservation [36]. To this end alter the veloc-ities at the nodes that would otherwise end up outside the glass domain, such that their normal component stays the same, i.e.
˜uki
·
n= uki ·n, (3.50)
while the tangential component uk
i
·
t is obtained by rotating the velocity vector, such that x ki ends up on
the equipment boundary (see Fig. 12). This can be formulated as
˜uki := αkiRk iu
k
t @ @t hhht((( t" " t Q Q Q Q !! !! !! t A A A U t Γp Γm Γf Ωtk H H H H H Y Ωtk+1 i xk i ti ˜xk+1 i F 11: C . where αk
i is the scaling parameter and R k
i is the 2 × 2 rotation matrix,
Rk i := cos γki − sin γki sin γki cos γki . (3.52)
The net outflow for the modified velocity field remains zero and hence the algorithm should give better mass conservation, i.e.
Z Γf u
·
n dΓ = Z Γf ˜u·
n dΓ = 0. (3.53)An approximate position of a node i that would be outside the glass domain at time tk+1is obtained using the modified velocity field,
˜xki+1= xki + ∆t ˜uki. (3.54) t J J J J J J Jt Q Q Q Q Q Q ix k i vk i ˜vk i ˜vk i· n @ @ @ @ R ˜vk i· t @ @ @ @ @ @ @ @ ti F 12: M .
Instead of (3.48) also an implicit time integration scheme can be used:
xki+1= xki + ∆t uki+1. (3.55)
In order to use (3.55) it is necessary to compute the flow velocity inΩg(tk+1). However, even though the
plunger position is known in the sense of (3.46), the free boundary ofΩg(tk+1) is still unknown. This
difficulty can be overcome by employing an algorithm that iterates on xk+1
i . Unfortunately, this
straight-forward approach requires each iteration the solution of a Stokes flow problem, which is computationally too costly. In [36] a numerical tool is introduced that overcomes the essential difficulty of the implicitness of the scheme by using the fact that the flow velocity is autonomous. In [69] this matter is examined into more detail.
3.5
Results
In this section simulations of the pressing of a jar and a bottle parison are shown. This includes the visualisation of the velocity and pressure fields, the tracking of the free boundaries of the glass flow domain and the computation of the motion of the plunger. For more results of the simulation tool used the reader is referred to [36]. Furthermore, see [30] for results of a different parison press model.
First consider the simulation of the pressing of a jar parison. The simulation starts at the moment the gob of glass has entered the mould and the baffle has closed. Figure 13 visualises the flow velocity field in the glass domain during pressing. Here full slip of glass at the equipment boundary is assumed, that is β = 0. The plunger moves upward forcing the glass to fill the space between the mould and the plunger. When the glass hits the mould the pressure increases. Figure 14 depicts the flow velocity of the glass in the neck part ring of the jar in the final stage of the parison press. Figure 15 shows the pressure field. It is important to know the pressure of the glass onto the mould during the process, so that a similar pressure can be applied from the outside in order to keep the separate parts of the mould together [36].
Next consider the simulation of the pressing of a bottle parison. The initial position of the top of the plunger in the simulation is close to the attachment of the ring to the mould. When the glass gob is dropped into the mould it almost reaches the ring before the plunger starts moving (see Fig. 16(a)). Figure 16 visualises the flow velocity of the glass during pressing. Figure 17 depicts the flow velocity in the bottle neck in the final stage of the parison press. Figure 18 shows the pressure field. The results are similar to the jar parison simulation. Again full slip of glass at the equipment boundary is assumed.
0 0.02 0.04 0.06 0.08 0.1 0.12 (a) 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 (b) 0 0.02 0.04 0.06 0.08 0.1 (c) 0 0.02 0.04 0.06 0.08 (d) F 13: P : [m/s] [36] (a) (b) (c) (d) F 14: P ( ): [36] 1 1.05 1.1 1.15 1.2 (a) 1 1.05 1.1 1.15 1.2 1.25 1.3 1.35 (b) 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 (c) 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 (d) F 15: P : [bar] [36]
0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 (a) 0 0.05 0.1 0.15 (b) 0 0.02 0.04 0.06 0.08 (c) 0 0.02 0.04 0.06 0.08 (d) F 16: P : [m/s] [36] (a) (b) (c) (d) F 17: P ( ): [36] 1 1.05 1.1 1.15 (a) 1 1.05 1.1 1.15 1.2 (b) 1 1.05 1.1 1.15 1.2 1.25 1.3 1.35 1.4 (c) 1 1.1 1.2 1.3 1.4 1.5 (d) F 18: P : [bar] [36]
4
Blow Model
This section presents a mathematical model for blowing. Section 4.1 specifies the mathematical model described in § 2 for the counter blow and the final blow. In addition to the flow of glass, also the flow of air is modelled. A particular question is how to deal with the glass-air interfaces. Section 4.2 discusses several techniques. Subsequently, Section 4.3 explains how a variational formulation is used to combine the physical problems for glass and air as well as the jump conditions on the interfaces into one statement. Section 4.4 presents a numerical simulation model for the blow-blow process. Finally, Section 4.5 shows some examples of process simulations. The simulation tool used for the results is presented in [23].
4.1
Mathematical Model
In Section 2.2 the balance laws for glass forming were formulated. In this section these balance laws are further specified for glass blowing. Since the final blow stage starts with the preform obtained in either the parison press stage of the counter blow stage, the orders of magnitude of most physical parameters for these forming processes are typically the same. The temperature of the glass melt in the final blow is usually slightly lower than in the preceding stage, but this does not lead to a significant difference in the order of magnitude of the physical parameters. Therefore, the same typical values for both the counter blow and the final blow are considered:
glass density : ρg = 2.5 · 103kg m−3
glass viscosity : ¯µg = 104kg m−1s−1
gravitational acceleration : g = 9.8 m s−2 flow velocity : V = 10−2m s−1 length scale of the parison : L = 10−2m
glass temperature : Tg = 1000◦C
mould temperature : Tm = 500◦C
specific heat of glass : ¯cp = 1.5 · 103J kg−1K−1
effective conductivity of glass : λ¯ = 5 W m−1K−1,
The main difference from the physical parameters in the parison press is that the time duration of a blow stage is typically much larger, hence the flow velocity is much smaller. The dimensionless numbers corre-sponding to either blow stage are:
Peglass≈ 76, Brglass≈ 4.0 · 10−5, Reglass≈ 2.5 · 10−5, Fr ≈ 1.0 · 10−3. (4.1)
The P´eclet number for glass is moderately large, while the Brinkman number is negligibly small. Thus the heat transfer is dominated by advection, convection and radiation:
∂T∗ ∂t∗ + u ∗
·
∇∗T∗= ∇∗·
λ∗∇∗T∗, (4.2) with λ∗:= λ ρ¯cpV L , (4.3)rather than its definition (2.9) in § 2.2.1. The Reynolds number for glass is sufficiently small to neglect the inertia terms in the Navier-Stokes equations (2.25). The order of magnitude of the gravity term is given by
Reglass
Fr ≈ 2.5 · 10
−2,
which is not extremely small. Therefore, the glass flow is described by the Stokes flow equations:
∇∗
·
T∗ = g∗, ∇∗·
u∗ = 0, (4.4)where T∗is the dimensionless stress tensor, which satisfies (2.20) in terms of the dimensionless variables, and g∗is the dimensionless gravity force, given by
g∗ = −Re
Frez. (4.5)
Subsequently, the air domain is considered. Note that the arguments to ignore the air domain in the parison press model do not apply to the blow model. Firstly, the inflow pressure is applied at the mould entrance and not directly on the glass, so that in this case the transport phenomena in air are of higher interest for the blow model. Secondly, because the influence of the heat flux cannot be ignored, the energy exchange between the glass and its surroundings should be taken into account. For hot pressurised air the following typical values are considered2:
initial air temperature : T0 = 750◦C,
specific heat of air : cp = 103J kg−1K−1,
thermal conductivity of air : λ = 10−1W m−1K−1.
air density : ρ = 1 kg m−3, air viscosity : µ = 10−4kg m−1s−1.
The resulting dimensionless numbers are:
Peair≈ 102, Brair≈ 4 · 10−8, Reair≈ 1. (4.6)
It can be concluded that the heat transfer in air is described by (4.2), while the flow of air is described by the full Navier-Stokes equations (2.25). This means that the model for the air flow is more complicated model than for the glass flow, while the motion of glass is most interesting. Therefore, air is replaced by a fictitious fluid with the same physical properties as air, but with a much higher viscosity, e.g. µa= 1. Then
the Reynolds number of the fictitious fluid Re ≈ 10−4is small enough to reasonably neglect the influence
of the inertia forces. On the other hand, the viscosity of the fictitious fluid is still much smaller than the viscosity of glass, so that the pressure drop in the air domain is negligible compared to the pressure drop in the glass domain [2]. Thus the flow of the ficitious fluid can be described by the Stokes flow equations (4.4).
Remark 4.1 If the glass temperature and the pressure in air can be reasonably assumed to be uniform, the calculations can be restricted to the glass domain and the glass-air interfaces can be treated as free boundaries with a prescribed pressure. In this case it can be recommended to use Boundary Element Methods to solve the flow problem [15, 16].
2The true orders of magnitude may be slightly different from their rough estimates in the table, but this will not affect the final results.
The boundary conditions for the flow problem are given in § 2.2, but can be further specified for glass blowing. Free-stress conditions are imposed for air at the mould wall. If no lubricate is used, usually a no-slip condition is imposed for glass at the mould wall. For modelling of the slip condition in the presence of a lubricate the reader is referred to [17, 18, 38]. As a result the boundary conditions for the flow problem become: u∗
·
n= 0, T∗n·
t = 0, on Γs, u∗·
n= 0, u∗·
t = 0, on Γ g,e, T∗n·
n= 0, T∗n·
t = 0, on Γa,e, T∗n·
n= p0, T∗n·
t = 0, on Γo, (4.7) with p0= pin on Γa,o 0 on Γg,o, (4.8)where pinis the pressure at which air is blown into the mould. Alternatively, one may prefer to introduce
the boundary condition
u∗
·
n= 0, T∗n·
t= 0, on Γg,o. (4.9)This boundary condition avoids outflow of glass through the mould entrance during blowing. However, this involves the definition of separate boundariesΓa,oandΓg,o, rather than the single boundaryΓo. Moreover,
these boundaries can change in time. In the next section it can be seen that this is not always convenient, particularly if a single fixed mesh is used for the discretisation of the flow domain. InsteadΓg,o can be
conceived as the boundary between the glass and the ring, thereby imposing a no-slip condition,
u∗= 0, on Γg,o. (4.10)
Note that in this caseΓa,oandΓg,oremain fixed.
The boundary conditions for the energy exchange problem are defined on the boundary of the flow domain: λ∗∇T∗
·
n= 0, on Γ s∪Γe, λ∗∇T∗·
n= Nu α∗ T∗− T∗ ∞, on Γo, (4.11)4.2
Glass-Air Interfaces
A two-phase fluid flow problem is considered, involving the flow of both glass and air. The flow domainΩ is described by the geometry of the mould and hence fixed. On the other hand, the air domainΩaand the
glass domainΩgare separated by moving interfacesΓi, as depicted in Fig. 4(c)-4(d), and therefore change
in time. In order to model the two-phase fluid flow problem, the glass-air interfaces have to captured. There are different numerical techniques to deal with the moving interfaces in two-phase fluid flow problems. They can be classified in two main categories [22]: tracking techniques (ITT) and interface-capturing techniques (ICT).
Interface-tracking techniques (ITT) attempt to find the moving interfaces explicitly. ITT involve sep-arate discretisations of domainsΩaandΩg; the meshes of both domains are updated as the flow evolves
