- Email: [email protected]
3D simulation of glass forming process
Journal of Materials Processing Technology 164–165 (2005) 1197–1203
3D simulation of glass forming process E. Feulvarch a,b , N. Moulin a,∗ , P. Saillard b , T. Lornage c , J.-M. Bergheau a a
Laboratoire de Tribologie et Dynamique des Syst`emes, UMR 5513 CNRS/ECL/ENISE, 58 rue Jean Parot, 42100 Saint Etienne, France b ESI Group, Le Discover, 84 Bd Vivier Merle, 69485 Lyon Cedex 03, France c BSN Glasspack, European Technical Center, St Romain en Gier, 69702 Givors, France
Abstract The simulation of glass forming is a complex thermo-mechanical problem. For high temperatures, the mechanical glass behaviour is assumed to be as a viscous incompressible material. Its viscosity is strongly temperature dependent and the process involves large displacements and large strains. Therefore, the finite element modelling of such a process needs to take account of large geometrical transformations. From the numerical point of view, this problem can be studied by means of a solid approach where the degrees of freedom are displacements. This modelling is highly non-linear. Moreover, the numerical convergence can be affected on account of the mechanical contact between the mould and the glass. The aim of this paper is to propose a linear finite element formulation for 3D glass forming where the degrees of freedom are velocity and hydrostatic pressure. The mechanical contact is taken into account by means of an explicit algorithm related to the time integration scheme. The thermal coupling is not considered. The results are discussed and compared to a non-linear formulation where the degrees of freedom are displacements. © 2005 Elsevier B.V. All rights reserved. Keywords: Finite element method; Glass forming; Mechanical contact
1. Introduction The simulation of glass forming is a complex thermomechanical problem. For high temperatures, the mechanical glass behaviour is assumed to be as a viscous incompressible material. Its viscosity is strongly temperature dependent and the process involves large displacements and large strains. Therefore, the modelling of such a process needs to take large transformations into account. From the numerical point of view, this problem can be studied by means of solid mechanics or fluid dynamics. As far as the finite element method is considered, the main difference comes from the degrees of freedom which are displacements in a solid approach and velocity in a fluid model. The thermal coupling is not considered. The aim of this paper is to present an explicit contact model related to a fluid approach in a linear modelling. ∗
Corresponding author. E-mail address: [email protected] (N. Moulin).
0924-0136/$ – see front matter © 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.jmatprotec.2005.02.135
Indeed, in a classical solid approach, the finite element method leads to a non-linear equation system. A 3D model cannot be studied easily because of convergence problems due to the iterative solver. Therefore, the modelling of glass forming for complex shape seems to be difficult. For this reason, a 3D linear model has been developed. It consists in treating implicitly the equilibrium equation and the conservation of mass in a velocity approach. At each time step, the displacements are calculated using an explicit method. The glass is assumed to stick on the mould and a linear time step calculation algorithm has been developed for this purpose. In the first part of this article, the industrial glass maker process is described and the context of the forming process is detailed. In the second part, the physical problem is presented. In the third section, the velocity finite element formulation is detailed with a special attention to the mechanical contact. The results are discussed and compared to a displacement finite element formulation implemented in the software SYSTUS® .
1198
E. Feulvarch et al. / Journal of Materials Processing Technology 164–165 (2005) 1197–1203
2. Glass maker process Before describing the simulation of glass making, it is useful to describe the properties of the raw material as well as the general process of manufacture from the fusion to its forming. Glass is a very attractive material and it is appreciated in particular for its visual qualities (transparency, hardness, etc.) and even in spite of its brittleness, glass is present in various forms: • Glass of packing for agroalimentary industry, the industry of the luxury (perfumery and cosmetics) and pharmaceutical industry. • Flat glass for the car and construction trade. • Glass of domestic use, the glassmaking of art. • Special glasses like glass TV, optical fibres and optical glasses. In opposition with crystallized solids – where the atoms are arranged according to a precise and periodic structure – glass is a non-crystalline solid with random distributed atoms in space. The mechanical properties of glass are isotropic. During the material annealing, the mechanical behaviour shows a transition from a liquid state into a solid state through a pasty fusion. Nowadays, 10–90% of the raw materials are recycled, the remainder being sand, soda and lime. First, this mixture is molten in the furnace at 1450 ◦ C. During fusion, it is necessary to homogenize the mixture and to remove it from its impurities. Subsequently, glass go through a siphon where it is maintained in rest in a basin before it is brought at the right temperature to obtain the desired viscosity. The molten glass is then cut out to form drops of glass of a predetermined weight as shown in Fig. 1. The drops of glass are brought to a first mould to form the preform, glass is then either puffed up using compressed air, or pressed by a punch against the preform mould. The preform is then transferred by inversion in the finishing mould where it will be puffed in order to obtain the final form. The
Fig. 2. Geometry of the glass forming process (section of 90◦ ).
aim of this work is to model this final glass forming step (Fig. 2).
3. Physical modelling For the glass forming process, we are interested in the glass shape evolution. The equilibrium of a continuous medium Ω is governed by two global conservation principles which are the stress equilibrium and the mass conservation. 3.1. Stress equilibrium The stress equilibrium equation must be satisfied for each material element. Including dynamic and static effects internally, and gravity as an external body force, the equation of motion in a homogenous medium Ω can be written in an Eulerian frame as follows: ¯ + ρ · g = ρ · → div(σ)
Dv Dt
(1)
where σ¯ is the Cauchy stress tensor, ρ is the density, g is the acceleration vector due to the gravity, v the velocity vector ∂(·) D and Dt is the material derivative operator D(·) · Dt = ∂t + v → grad(·).In this study, the inertial and gravity forces can be considered as much smaller than the static effects to be neglected: ¯ = 0 → div(σ)
(2)
At typically glass forming temperatures, it is reasonable to assume that the large strains are perfectly viscoplastic. For such a mechanical behaviour, the total stress tensor σ¯ can be decomposed into hydrostatic and viscous parts as for a viscous incompressible fluid: Fig. 1. Industrial process.
σ¯ = S¯ − p · I¯
(3)
E. Feulvarch et al. / Journal of Materials Processing Technology 164–165 (2005) 1197–1203
1199
where p denotes the hydrostatic pressure, I¯ represents the unit tensor and S¯ is the viscous stress tensor. The stress equilibrium equation can therefore be written as: ¯ − → grad(p) = 0 → div(S)
(4)
The viscous stress tensor S¯ is related to the strain-rate tensor ¯ through the dynamic viscosity µ: D ¯ S¯ = 2 · µ · D
(5) Fig. 3. Regularly shaped tetrahedral element.
¯ is the strain-rate tensor defined by means of the where D velocity field v ¯ = D
1 2
T
· (grad(v) + grad (v))
(6)
The boundary conditions on the surface ∂Ω of a medium Ω can be written as: • Dirichlet boundary conditions: v = vd on surface ∂Ωv . • Neumann boundary conditions: σ¯ · n = T on surface ∂ΩT with ∂Ω = ∂Ωv ∪ ∂ΩT and ∂Ωv ∩ ∂ΩT = ∅.
a specific tetrahedral element has been implemented. The velocity field is discretized with continuous piecewise linear functions N1 , N2 , N3 and N4 associated to nodes 1, 2, 3 and 4 enriched by a bubble function Nb associated to a fifth node b which is the centroid of the element [3,4]. These functions can be developed for a regularly shaped element as shown in Fig. 3. The velocity, at each time step t, is expressed in the following way at each point (ξ, η, ζ) in the element
3.2. Conservation of mass The continuity equation states that mass cannot be lost or gained. It implies that velocity field must be well behaved: ∂ρ + div(ρ · v) = 0 (7) ∂t Considering a constant viscosity, the continuity equation becomes: div(v) = 0
(8)
This equation can be referred to an incompressibility condition.
vt (ξ, η, ζ) =
4 i=1
Ni (ξ, η, ζ) · vti + Nb (ξ, η, ζ) · vtb
(9)
where vi is the velocity at each node i and vb is the velocity at node b. The pressure is approximated by means of the linear functions N1 , N2 , N3 and N4 : pt (ξ, η, ζ) =
4 i=1
Ni (ξ, η, ζ) · pti
(10)
where pti is the pressure at each node i = 1, 2, 3, 4.
4. Finite element model
4.2. Finite element formulation
In this section, a linear finite element formulation is detailed. In this model, the degrees of freedom are velocity and hydrostatic pressure. The stress equilibrium equation and the mass conservation equation are also called the Stokes equations. To approximate the solution of this problem, the velocity space discretization and the pressure space discretization must be properly chosen to satisfy the inf-sup condition of the mixed method theory [10].
The finite element method is based on the variational formulation of the problem obtained from the weighted residual method and the Green–Gauss theorem (integration by parts) [5,6]. When the finite element method is applied, at each time step t, to Eqs. (4) and (8), it leads to the solution of a linear equation system (residual vector):
4.1. Discretization
m A
e (v, p) R v
A
ep (v, p) R
e=1 m e=1
The inf-sup condition of the mixed method theory imposes the use of specific finite elements. One possibility considers a classical velocity discretization and a constant pressure in the finite element [1]. Other possibility takes account of the incompressibility thanks to a penalty method [2]. In this study,
=
v (v, p) R p (v, p) R
= {0}
(11)
m
where A denotes the element matrices assembling symbol. e=1
e (v, p)k and Rep (v, p)k represent the kth components of R v e (v, p) and R ep (v, p) corresponding to the element vectors R v
1200
E. Feulvarch et al. / Journal of Materials Processing Technology 164–165 (2005) 1197–1203
surface (Fig. 4) [14]. The explicit Euler forward integration scheme allows us to calculate easily a corrected time step t* considering the positive lowest value of ηi as follows: t ∗ = MIN(ηi ) ·
t
(17)
The displacements can therefore be calculated by means of this new times step t* : xt+
t∗
= xt + u t = xt + vt ·
t∗
(18)
t∗
From the new geometry xt+ , each node of the glass is considered as being into contact when the distance between the node and its orthogonal projection on the mould surface
Fig. 4. Contact modelling.
node k:
e ∂Nk · D11 + ∂x Ωe ∂N e k e (v, p)k = R · D12 + v Ωe ∂x ∂Nke · D13 + ∂x Ωe e Rp (v, p)k = Nke · div(v) dv Ωe
∂Nke · D12 + ∂y ∂Nke · D22 + ∂y ∂Nke · D23 + ∂y
∂Nke · D13 dv − ∂z
Ωe ∂Nke · D23 dv − ∂z
Ωe ∂Nke · D33 dv − ∂z Ωe
t
(14)
where t is the time step which needs to be relatively small [6]. u t corresponds to the displacements between the geometry at time t and the geometry at time t + t. Therefore, the geometry is updated by the following way: xt+
t
= xt + u t = xt + vt ·
t
is lower than a value ε. In this case, the velocity is assumed to be equal to zero. Several nodes can therefore be detected into contact at the same time t + t* avoiding a drastic decrease of the time step t* when the number of nodes becoming into contact is large. This procedure is used at each time step.
5. Glass forming simulation The objective of the simulation for the glass forming process is the study of the final geometry in term
(15)
For the treatment of the mechanical contact, the literature exhibits a large range of numerical methods as penalty methods [7–9], Lagrange multipliers [7,10] or augmented Lagrangians [11–13]. In this study, the glass is supposed to stick on the mould which is assumed to be as a rigid body. A coefficient ηi is calculated for each node becoming into contact with the mould between t and t + t as shown in Fig. 4 ηi =
t+ηi · t − X t || ||X i i t+ t t || ||X −X i
(12)
(13)
The integrals are evaluated numerically by means of the Gauss quadrature [6]. The solution of this linear equation system gives the velocity vt and the pressure pt at each node. One can note that the solution is obtained in an Eulerian frame where the large deformations are taken into account by means ¯ of the strain rates tensor D. t The displacements u are calculated in an Euler explicit way [5] as follows: u t = vt ·
∂Nke · p dv − Nke · T · x ds ∂x ∂Ωe∩∂ΩT e ∂Nk e · p dv − Nk · T · y ds ∂y ∂Ωe∩∂ΩT ∂Nke e Nk · T · z ds · p dv − ∂z ∂Ωe∩∂ΩT
(16)
i
t corresponds to node i at t, X t+ t represents its position at X i i t+ηi · t is its orthogonal projection on the mould t+ t and X i
Fig. 5. 2D mesh for the SYSTUS® model.
E. Feulvarch et al. / Journal of Materials Processing Technology 164–165 (2005) 1197–1203
1201
of glass thickness. In this section, the final geometry and its evolution during the process are compared. In this section, a forming simulation is presented for a small bottle. At first, a 2D simulation has been done using the software SYSTUS® [15] in an axisymetric option (Fig. 5). The degrees of freedom are the displacements. In this case, the mechanical behaviour is perfectly viscoplastic. A 3D simulation has been also studied for a section of 10◦ with the velocity/pressure formulation presented in the previous section (Fig. 6). For the boundary conditions, the mould is assumed to be rigid (velocity equal to zero). The nodes of the glass which are located on the axisymetric axis are allowed to move only along this axis. Fig. 6. 3D mesh for the velocity approach.
Fig. 7. 3D simulation based on the velocity/pressure formulation.
1202
E. Feulvarch et al. / Journal of Materials Processing Technology 164–165 (2005) 1197–1203
The nodes which are on the top part of the glass are not allowed to move to maintain the perform as in the industrial process. Dynamic viscosity values can be found in the literature [16]. In this example, the dynamic viscosity is taken equal to 5 × 10−3 MPa s. For the 3D model, the nodes which are on the two lateral faces are not allowed to move along the orthoradial directions to model the interaction with the other part of the bottle which is not modelled. In this example, a depression is imposed between the mould and the glass. The geometry evolution is drawn in Fig. 7 for the 3D model. The final shapes are superimposed in Fig. 8. The results obtained from the velocity/pressure formulation compare well with the SYSTUS® model. In Fig. 9, the norm of the displacement for two nodes is plotted and superimposed to compare the evolution during the forming process. The displacement results are quite similar. Fig. 8. Final shape for both calculations.
Fig. 9. Comparison of the displacements.
E. Feulvarch et al. / Journal of Materials Processing Technology 164–165 (2005) 1197–1203
6. Conclusion A model has been studied with the software SYSTUS® for a classical solid approach where the degrees of freedom are displacements. This modelling takes large deformations into account. The solution is obtained by means of an iterative solver. In this paper, a linear finite element formulation has been presented without thermal coupling. The results show a good correlation in term of geometry evolution. Moreover, the final shapes obtained by means of each model are quite similar as shown in Fig. 9. The model presented in this article could model complex 3D problems. In terms of numerical treatment, the presented model only requires a linear solver removing convergence problems. This formulation is therefore much faster than the solid model (displacements) which is based on an iterative solver. Anyway, there are still some physical considerations which are not taken into account. For example, the thermal coupling could be implemented so as to consider the density gradient in the glass. One can note that a low thermomechanical coupling model will not be affected the linear characteristic of the velocity/pressure formulation.
References [1] J.N. Reddy, D.K. Gartling, The Finite Element Method in Heat Transfer and Fluid Dynamics, CRC Press LLC, 2001.
1203
[2] R.H. Wagoner, J.-L. Chenot, Metal Forming Analysis, Cambridge University Press, 2001. [3] D.N. Arnold, F. Brezzi, M. Fortin, A stable element for the Stokes equations, Clacolo 21 (1984) 337–344. [4] A. Ern, J.-L. Guermond, El´ements Finis: th´eorie, applications, mise en oeuvre, Springer, 2002. [5] T. Belytschko, W.K. Liu, B. Moran, Nonlinear Finite Elements for Continua and Structures, Wiley, 2001. [6] K.J. Bathe, Finite Element Procedure, Prentice-Hall, 1996. [7] T. Belytschko, M.O. Neal, Contact-impact by the pinball algorithm with penalty and Lagrangian methods, Int. J. Numer. Meth. Eng. 31 (1991) 547–572. [8] D. Peric, D.R.J. Owen, Computational model for 3D contact problems with friction based on the penalty method, Int. J. Numer. Meth. Eng. 35 (1992) 1289–1309. [9] P. Papadopoulos, R.L. Taylor, A simple algorithm for threedimensional finite element analysis of contact problems, Comput. Struct. 46 (1993) 1107–1118. [10] A.B. Chaudhary, K.J. Bathe, A solution method for static and dynamic analysis for three-dimensional contact problems with friction, Comput. Struct. 24 (6) (1986) 855–873. [11] T.A. Laursen, B.N. Maker, An augmented Lagrangian quasi-Newton solver for constrained non linear finite element applications, Int. J. Numer. Meth. Eng. 38 (1995) 3571–3590. [12] T.A. Laursen, B.N. Maker, Algorithmic symetrisation of coulomb frictional problems using augmented Lagrangians, Comput. Meth. Appl. Mech. Eng. 108 (1993) 133–146. [13] J.C. Simo, T.A. Laursen, An augmented Lagrangian treatment of contact involving friction, Comput. Struct. 42 (1) (1992) 97–116. [14] M. Oldenburg, L. Nilsson, The position code algorithm for contact searching, Int. J. Meth. Eng. 37 (1994) 359–386. [15] SYSTUS® , User’s Manual, ESI Group, 2004. [16] G.S. Fulcher, Analysis of recent measurements of the viscosity of glasses, J. Am. Ceram. Soc. 8 (6) (1925) 339–355.
3D simulation of glass forming process E. Feulvarch a,b , N. Moulin a,∗ , P. Saillard b , T. Lornage c , J.-M. Bergheau a a
Laboratoire de Tribologie et Dynamique des Syst`emes, UMR 5513 CNRS/ECL/ENISE, 58 rue Jean Parot, 42100 Saint Etienne, France b ESI Group, Le Discover, 84 Bd Vivier Merle, 69485 Lyon Cedex 03, France c BSN Glasspack, European Technical Center, St Romain en Gier, 69702 Givors, France
Abstract The simulation of glass forming is a complex thermo-mechanical problem. For high temperatures, the mechanical glass behaviour is assumed to be as a viscous incompressible material. Its viscosity is strongly temperature dependent and the process involves large displacements and large strains. Therefore, the finite element modelling of such a process needs to take account of large geometrical transformations. From the numerical point of view, this problem can be studied by means of a solid approach where the degrees of freedom are displacements. This modelling is highly non-linear. Moreover, the numerical convergence can be affected on account of the mechanical contact between the mould and the glass. The aim of this paper is to propose a linear finite element formulation for 3D glass forming where the degrees of freedom are velocity and hydrostatic pressure. The mechanical contact is taken into account by means of an explicit algorithm related to the time integration scheme. The thermal coupling is not considered. The results are discussed and compared to a non-linear formulation where the degrees of freedom are displacements. © 2005 Elsevier B.V. All rights reserved. Keywords: Finite element method; Glass forming; Mechanical contact
1. Introduction The simulation of glass forming is a complex thermomechanical problem. For high temperatures, the mechanical glass behaviour is assumed to be as a viscous incompressible material. Its viscosity is strongly temperature dependent and the process involves large displacements and large strains. Therefore, the modelling of such a process needs to take large transformations into account. From the numerical point of view, this problem can be studied by means of solid mechanics or fluid dynamics. As far as the finite element method is considered, the main difference comes from the degrees of freedom which are displacements in a solid approach and velocity in a fluid model. The thermal coupling is not considered. The aim of this paper is to present an explicit contact model related to a fluid approach in a linear modelling. ∗
Corresponding author. E-mail address: [email protected] (N. Moulin).
0924-0136/$ – see front matter © 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.jmatprotec.2005.02.135
Indeed, in a classical solid approach, the finite element method leads to a non-linear equation system. A 3D model cannot be studied easily because of convergence problems due to the iterative solver. Therefore, the modelling of glass forming for complex shape seems to be difficult. For this reason, a 3D linear model has been developed. It consists in treating implicitly the equilibrium equation and the conservation of mass in a velocity approach. At each time step, the displacements are calculated using an explicit method. The glass is assumed to stick on the mould and a linear time step calculation algorithm has been developed for this purpose. In the first part of this article, the industrial glass maker process is described and the context of the forming process is detailed. In the second part, the physical problem is presented. In the third section, the velocity finite element formulation is detailed with a special attention to the mechanical contact. The results are discussed and compared to a displacement finite element formulation implemented in the software SYSTUS® .
1198
E. Feulvarch et al. / Journal of Materials Processing Technology 164–165 (2005) 1197–1203
2. Glass maker process Before describing the simulation of glass making, it is useful to describe the properties of the raw material as well as the general process of manufacture from the fusion to its forming. Glass is a very attractive material and it is appreciated in particular for its visual qualities (transparency, hardness, etc.) and even in spite of its brittleness, glass is present in various forms: • Glass of packing for agroalimentary industry, the industry of the luxury (perfumery and cosmetics) and pharmaceutical industry. • Flat glass for the car and construction trade. • Glass of domestic use, the glassmaking of art. • Special glasses like glass TV, optical fibres and optical glasses. In opposition with crystallized solids – where the atoms are arranged according to a precise and periodic structure – glass is a non-crystalline solid with random distributed atoms in space. The mechanical properties of glass are isotropic. During the material annealing, the mechanical behaviour shows a transition from a liquid state into a solid state through a pasty fusion. Nowadays, 10–90% of the raw materials are recycled, the remainder being sand, soda and lime. First, this mixture is molten in the furnace at 1450 ◦ C. During fusion, it is necessary to homogenize the mixture and to remove it from its impurities. Subsequently, glass go through a siphon where it is maintained in rest in a basin before it is brought at the right temperature to obtain the desired viscosity. The molten glass is then cut out to form drops of glass of a predetermined weight as shown in Fig. 1. The drops of glass are brought to a first mould to form the preform, glass is then either puffed up using compressed air, or pressed by a punch against the preform mould. The preform is then transferred by inversion in the finishing mould where it will be puffed in order to obtain the final form. The
Fig. 2. Geometry of the glass forming process (section of 90◦ ).
aim of this work is to model this final glass forming step (Fig. 2).
3. Physical modelling For the glass forming process, we are interested in the glass shape evolution. The equilibrium of a continuous medium Ω is governed by two global conservation principles which are the stress equilibrium and the mass conservation. 3.1. Stress equilibrium The stress equilibrium equation must be satisfied for each material element. Including dynamic and static effects internally, and gravity as an external body force, the equation of motion in a homogenous medium Ω can be written in an Eulerian frame as follows: ¯ + ρ · g = ρ · → div(σ)
Dv Dt
(1)
where σ¯ is the Cauchy stress tensor, ρ is the density, g is the acceleration vector due to the gravity, v the velocity vector ∂(·) D and Dt is the material derivative operator D(·) · Dt = ∂t + v → grad(·).In this study, the inertial and gravity forces can be considered as much smaller than the static effects to be neglected: ¯ = 0 → div(σ)
(2)
At typically glass forming temperatures, it is reasonable to assume that the large strains are perfectly viscoplastic. For such a mechanical behaviour, the total stress tensor σ¯ can be decomposed into hydrostatic and viscous parts as for a viscous incompressible fluid: Fig. 1. Industrial process.
σ¯ = S¯ − p · I¯
(3)
E. Feulvarch et al. / Journal of Materials Processing Technology 164–165 (2005) 1197–1203
1199
where p denotes the hydrostatic pressure, I¯ represents the unit tensor and S¯ is the viscous stress tensor. The stress equilibrium equation can therefore be written as: ¯ − → grad(p) = 0 → div(S)
(4)
The viscous stress tensor S¯ is related to the strain-rate tensor ¯ through the dynamic viscosity µ: D ¯ S¯ = 2 · µ · D
(5) Fig. 3. Regularly shaped tetrahedral element.
¯ is the strain-rate tensor defined by means of the where D velocity field v ¯ = D
1 2
T
· (grad(v) + grad (v))
(6)
The boundary conditions on the surface ∂Ω of a medium Ω can be written as: • Dirichlet boundary conditions: v = vd on surface ∂Ωv . • Neumann boundary conditions: σ¯ · n = T on surface ∂ΩT with ∂Ω = ∂Ωv ∪ ∂ΩT and ∂Ωv ∩ ∂ΩT = ∅.
a specific tetrahedral element has been implemented. The velocity field is discretized with continuous piecewise linear functions N1 , N2 , N3 and N4 associated to nodes 1, 2, 3 and 4 enriched by a bubble function Nb associated to a fifth node b which is the centroid of the element [3,4]. These functions can be developed for a regularly shaped element as shown in Fig. 3. The velocity, at each time step t, is expressed in the following way at each point (ξ, η, ζ) in the element
3.2. Conservation of mass The continuity equation states that mass cannot be lost or gained. It implies that velocity field must be well behaved: ∂ρ + div(ρ · v) = 0 (7) ∂t Considering a constant viscosity, the continuity equation becomes: div(v) = 0
(8)
This equation can be referred to an incompressibility condition.
vt (ξ, η, ζ) =
4 i=1
Ni (ξ, η, ζ) · vti + Nb (ξ, η, ζ) · vtb
(9)
where vi is the velocity at each node i and vb is the velocity at node b. The pressure is approximated by means of the linear functions N1 , N2 , N3 and N4 : pt (ξ, η, ζ) =
4 i=1
Ni (ξ, η, ζ) · pti
(10)
where pti is the pressure at each node i = 1, 2, 3, 4.
4. Finite element model
4.2. Finite element formulation
In this section, a linear finite element formulation is detailed. In this model, the degrees of freedom are velocity and hydrostatic pressure. The stress equilibrium equation and the mass conservation equation are also called the Stokes equations. To approximate the solution of this problem, the velocity space discretization and the pressure space discretization must be properly chosen to satisfy the inf-sup condition of the mixed method theory [10].
The finite element method is based on the variational formulation of the problem obtained from the weighted residual method and the Green–Gauss theorem (integration by parts) [5,6]. When the finite element method is applied, at each time step t, to Eqs. (4) and (8), it leads to the solution of a linear equation system (residual vector):
4.1. Discretization
m A
e (v, p) R v
A
ep (v, p) R
e=1 m e=1
The inf-sup condition of the mixed method theory imposes the use of specific finite elements. One possibility considers a classical velocity discretization and a constant pressure in the finite element [1]. Other possibility takes account of the incompressibility thanks to a penalty method [2]. In this study,
=
v (v, p) R p (v, p) R
= {0}
(11)
m
where A denotes the element matrices assembling symbol. e=1
e (v, p)k and Rep (v, p)k represent the kth components of R v e (v, p) and R ep (v, p) corresponding to the element vectors R v
1200
E. Feulvarch et al. / Journal of Materials Processing Technology 164–165 (2005) 1197–1203
surface (Fig. 4) [14]. The explicit Euler forward integration scheme allows us to calculate easily a corrected time step t* considering the positive lowest value of ηi as follows: t ∗ = MIN(ηi ) ·
t
(17)
The displacements can therefore be calculated by means of this new times step t* : xt+
t∗
= xt + u t = xt + vt ·
t∗
(18)
t∗
From the new geometry xt+ , each node of the glass is considered as being into contact when the distance between the node and its orthogonal projection on the mould surface
Fig. 4. Contact modelling.
node k:
e ∂Nk · D11 + ∂x Ωe ∂N e k e (v, p)k = R · D12 + v Ωe ∂x ∂Nke · D13 + ∂x Ωe e Rp (v, p)k = Nke · div(v) dv Ωe
∂Nke · D12 + ∂y ∂Nke · D22 + ∂y ∂Nke · D23 + ∂y
∂Nke · D13 dv − ∂z
Ωe ∂Nke · D23 dv − ∂z
Ωe ∂Nke · D33 dv − ∂z Ωe
t
(14)
where t is the time step which needs to be relatively small [6]. u t corresponds to the displacements between the geometry at time t and the geometry at time t + t. Therefore, the geometry is updated by the following way: xt+
t
= xt + u t = xt + vt ·
t
is lower than a value ε. In this case, the velocity is assumed to be equal to zero. Several nodes can therefore be detected into contact at the same time t + t* avoiding a drastic decrease of the time step t* when the number of nodes becoming into contact is large. This procedure is used at each time step.
5. Glass forming simulation The objective of the simulation for the glass forming process is the study of the final geometry in term
(15)
For the treatment of the mechanical contact, the literature exhibits a large range of numerical methods as penalty methods [7–9], Lagrange multipliers [7,10] or augmented Lagrangians [11–13]. In this study, the glass is supposed to stick on the mould which is assumed to be as a rigid body. A coefficient ηi is calculated for each node becoming into contact with the mould between t and t + t as shown in Fig. 4 ηi =
t+ηi · t − X t || ||X i i t+ t t || ||X −X i
(12)
(13)
The integrals are evaluated numerically by means of the Gauss quadrature [6]. The solution of this linear equation system gives the velocity vt and the pressure pt at each node. One can note that the solution is obtained in an Eulerian frame where the large deformations are taken into account by means ¯ of the strain rates tensor D. t The displacements u are calculated in an Euler explicit way [5] as follows: u t = vt ·
∂Nke · p dv − Nke · T · x ds ∂x ∂Ωe∩∂ΩT e ∂Nk e · p dv − Nk · T · y ds ∂y ∂Ωe∩∂ΩT ∂Nke e Nk · T · z ds · p dv − ∂z ∂Ωe∩∂ΩT
(16)
i
t corresponds to node i at t, X t+ t represents its position at X i i t+ηi · t is its orthogonal projection on the mould t+ t and X i
Fig. 5. 2D mesh for the SYSTUS® model.
E. Feulvarch et al. / Journal of Materials Processing Technology 164–165 (2005) 1197–1203
1201
of glass thickness. In this section, the final geometry and its evolution during the process are compared. In this section, a forming simulation is presented for a small bottle. At first, a 2D simulation has been done using the software SYSTUS® [15] in an axisymetric option (Fig. 5). The degrees of freedom are the displacements. In this case, the mechanical behaviour is perfectly viscoplastic. A 3D simulation has been also studied for a section of 10◦ with the velocity/pressure formulation presented in the previous section (Fig. 6). For the boundary conditions, the mould is assumed to be rigid (velocity equal to zero). The nodes of the glass which are located on the axisymetric axis are allowed to move only along this axis. Fig. 6. 3D mesh for the velocity approach.
Fig. 7. 3D simulation based on the velocity/pressure formulation.
1202
E. Feulvarch et al. / Journal of Materials Processing Technology 164–165 (2005) 1197–1203
The nodes which are on the top part of the glass are not allowed to move to maintain the perform as in the industrial process. Dynamic viscosity values can be found in the literature [16]. In this example, the dynamic viscosity is taken equal to 5 × 10−3 MPa s. For the 3D model, the nodes which are on the two lateral faces are not allowed to move along the orthoradial directions to model the interaction with the other part of the bottle which is not modelled. In this example, a depression is imposed between the mould and the glass. The geometry evolution is drawn in Fig. 7 for the 3D model. The final shapes are superimposed in Fig. 8. The results obtained from the velocity/pressure formulation compare well with the SYSTUS® model. In Fig. 9, the norm of the displacement for two nodes is plotted and superimposed to compare the evolution during the forming process. The displacement results are quite similar. Fig. 8. Final shape for both calculations.
Fig. 9. Comparison of the displacements.
E. Feulvarch et al. / Journal of Materials Processing Technology 164–165 (2005) 1197–1203
6. Conclusion A model has been studied with the software SYSTUS® for a classical solid approach where the degrees of freedom are displacements. This modelling takes large deformations into account. The solution is obtained by means of an iterative solver. In this paper, a linear finite element formulation has been presented without thermal coupling. The results show a good correlation in term of geometry evolution. Moreover, the final shapes obtained by means of each model are quite similar as shown in Fig. 9. The model presented in this article could model complex 3D problems. In terms of numerical treatment, the presented model only requires a linear solver removing convergence problems. This formulation is therefore much faster than the solid model (displacements) which is based on an iterative solver. Anyway, there are still some physical considerations which are not taken into account. For example, the thermal coupling could be implemented so as to consider the density gradient in the glass. One can note that a low thermomechanical coupling model will not be affected the linear characteristic of the velocity/pressure formulation.
References [1] J.N. Reddy, D.K. Gartling, The Finite Element Method in Heat Transfer and Fluid Dynamics, CRC Press LLC, 2001.
1203
[2] R.H. Wagoner, J.-L. Chenot, Metal Forming Analysis, Cambridge University Press, 2001. [3] D.N. Arnold, F. Brezzi, M. Fortin, A stable element for the Stokes equations, Clacolo 21 (1984) 337–344. [4] A. Ern, J.-L. Guermond, El´ements Finis: th´eorie, applications, mise en oeuvre, Springer, 2002. [5] T. Belytschko, W.K. Liu, B. Moran, Nonlinear Finite Elements for Continua and Structures, Wiley, 2001. [6] K.J. Bathe, Finite Element Procedure, Prentice-Hall, 1996. [7] T. Belytschko, M.O. Neal, Contact-impact by the pinball algorithm with penalty and Lagrangian methods, Int. J. Numer. Meth. Eng. 31 (1991) 547–572. [8] D. Peric, D.R.J. Owen, Computational model for 3D contact problems with friction based on the penalty method, Int. J. Numer. Meth. Eng. 35 (1992) 1289–1309. [9] P. Papadopoulos, R.L. Taylor, A simple algorithm for threedimensional finite element analysis of contact problems, Comput. Struct. 46 (1993) 1107–1118. [10] A.B. Chaudhary, K.J. Bathe, A solution method for static and dynamic analysis for three-dimensional contact problems with friction, Comput. Struct. 24 (6) (1986) 855–873. [11] T.A. Laursen, B.N. Maker, An augmented Lagrangian quasi-Newton solver for constrained non linear finite element applications, Int. J. Numer. Meth. Eng. 38 (1995) 3571–3590. [12] T.A. Laursen, B.N. Maker, Algorithmic symetrisation of coulomb frictional problems using augmented Lagrangians, Comput. Meth. Appl. Mech. Eng. 108 (1993) 133–146. [13] J.C. Simo, T.A. Laursen, An augmented Lagrangian treatment of contact involving friction, Comput. Struct. 42 (1) (1992) 97–116. [14] M. Oldenburg, L. Nilsson, The position code algorithm for contact searching, Int. J. Meth. Eng. 37 (1994) 359–386. [15] SYSTUS® , User’s Manual, ESI Group, 2004. [16] G.S. Fulcher, Analysis of recent measurements of the viscosity of glasses, J. Am. Ceram. Soc. 8 (6) (1925) 339–355.