Lagrangian finite element model for the 3D simulation of glass forming processes

==1 Introduction==

In spite of its very long history, bottle manufacturing remains a challenging process and requires further improvements. The continuously growing competition calls for the optimization of the existing processes, diminishing the risks of producing deficient bottles. Thus, such optimization must be based upon a detailed knowledge of the process variables (such as the final topology, wall thickness, stress and temperature distribution) and their dependence upon the input parameters (inlet pressure, cooling conditions, etc). The typical questions that need to be answered are:  How can the container weight be reduced without compromising on its strength? What are the optimal operation conditions (air pressure, mould temperature) and duration of the different forming stages?    Up-to-date the answers to these and similar questions are predominantly based upon the experience and craftsmanship rather than scientific knowledge.

Numerical modeling and simulation can serve as an efficient tool for answering many questions arising when facing unexpected  effects in the real products. Apart from being considerably cheaper than conducting expensive and time-consuming trial-and-error procedures common to factories, only numerical simulation can  provide such (otherwise impossible or difficult to obtain) results as: stress distributions within the solidifying melt and temperature distribution.

Up-to-date, there exist several computational tools used by industries for bottle manufacturing simulation.  Usually these software model the glass manufacturing process using axis-symmetric formulations.  This approximation greatly reduces the associated computational costs.  However, it over-simplifies the process: even though many bottle molds are purely axis-symmetric, nearly all containers produced do have non-symmetrical thickness distributions. Additionally, axisymmetric formulations cannot be applied to modeling bottles with non-circular cross-sections, such as e.g. fragrance containers. Thus, 3D simulations appear to be obligatory for obtaining reliable predictions. However, 3D simulations are typically characterized by excessive computational times.

Generally, two classes of methods can be applied to the problem at hand: the fixed mesh (Eulerian) and the mesh-moving (Lagrangian or Arbitrary Lagrangian Eulerian (ALE)) ones. Eulerian formulations require excessively fine meshes for the correct representation of the domain evolution and typically introduce errors in mass conservation (whenever the Level Set method is used for representing the glass-air interface evolution).  On the other hand, Lagrangian approaches lead to strongly non-linear systems of equations and large mesh deformations. Thus, a robust and computationally efficient 3D model still presents a major challenge.

In the present work a 3D viscous incompressible fluid formulation using an updated Lagrangian framework is proposed, where the current configuration is the reference one.  It adopts the basic features of the Particle Finite Element Method (PFEM) <span id='citeF-1'></span>[[#cite-1|[1]]]. The key idea of the PFEM is that the variables of interest are stored at the nodes instead of the Gauss points. This results in a hybrid between a standard FE and a mesh-free method. A finite element mesh is created at every time step of the transient problem and the solution is then stored at the nodes. At every time step the governing equations are solved in the standard Finite Element (FE) fashion. The discrete operators are updated at every non-linear iteration step according to the newly obtained domain configuration, ensuring excellent convergence of the iterative procedure. The nodes obtain their new positions and the mesh is re-generated using an unconstrained Delaunay technique. The approach is adapted to the problem at hand introducing a simple but efficient contact algorithm, boundary tracking and a re-meshing strategy. In terms of the method for solving the governing equations, we use a modified fractional step approach, combining the classical technique with a quasi-incompressible prediction <span id='citeF-2'></span>[[#cite-2|[2]]], <span id='citeF-3'></span>[[#cite-3|[3]]]. The approach on the one hand allows for a highly computationally efficient solution and, on the other hand, leads to an accurate mass conservation, which is essential for the problems of interest.

The paper concludes with two numerical examples. The first one is used for the validation of the model. The second one shows the potential of the method. Moreover, it can be used as a reference for the comparison and validation  of the future models for the bottle forming simulation. Even though several glass forming simulation results are published in literature (<span id='citeF-4'></span>[[#cite-4|[4]]], <span id='citeF-5'></span>[[#cite-5|[5]]], <span id='citeF-6'></span>[[#cite-6|[6]]], <span id='citeF-7'></span>[[#cite-7|[7]]]), there exist no well-established benchmark up-to-date. The example we propose focuses on the modeling of the ''final blow'' stage of a glass manufacturing process. The material and geometrical data as well as all the boundary conditions necessary for reproducing the example are specified.

==2 Glass forming==

===2.1 Industrial process===

Prior to presenting the model for the glass  (in particular, bottle) forming, let us review the industrial process and introduce the corresponding terminology. The standard bottle manufacturing process is sketched on Fig. [[#img-1|1]].

<div id='img-1'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_García-Espinosa_370675429-STEPS_GLASS_FORMING.png|600px|Different stages of bottle manufacturing process]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 1:''' Different stages of bottle manufacturing process
|}

Typically, high speed machines are fed a stream of molten glass<span id="fnc-1"></span>[[#fn-1|<sup>1</sup>]] that is cut with a shearing blade into “gobs” of predetermined weight.  These gobs fall into the first blank mold as shown in Fig. [[#img-1|1]] a), where the temperature drops to the so-called “working temperature range” (some 1150<math>\quad ^{\circ }</math>C for soda-lime glasses).  At the base of the mold a cylindrical plunger for shaping the bottle neck is inserted. Air pressure or a plunger is applied in order to push the gob to the bottom of the mold (Fig. [[#img-1|1]] b)). Afterwards, air compressed to <math display="inline">\approx </math> 0,15 MPa is blown from the bottom of the mold forcing the gob to rise and take the shape of the mold (Fig. [[#img-1|1]] c)). This stage is known as a ''counter-blow'' process, suggesting that the “blow” is performed against gravity. The intermediate product of the counter-blow is known as “parison”.

Afterwards, the parison is removed from the first mold, turned upside down (Fig. [[#img-1|1]] d)) and transferred into the second mold, where it is hung in order to stretch due to gravity  (Fig. [[#img-1|1]] e)), usually until the contact with the mold bottom is established. Finally, air pressure (slightly lower than the one used in the counter blow mold) is applied leading to the final shape of the bottle (Fig. [[#img-1|1]] f)). This stage is known as ''final blow'' process.  The bottle is then removed from the mold and is transferred to the annealing oven where it is reheated to remove the residual stresses produced during forming and finally it is cooled to the ambient temperature. The forming process, from the time when the gob is dropped into the first mold until the final product is removed from the second mold takes around 6 seconds.

<span id="fn-1"></span>
<span style="text-align: center; font-size: 75%;">([[#fnc-1|<sup>1</sup>]]) The most prevalent glass used for glass containers is soda-lime <math>Na_2O-CaO-SiO_2</math> glass.</span>

===2.2 Material properties===

Glass is a visco-elastic material. At low temperatures (roughly, below 400 <math>\quad ^{\circ }</math>C) elastic effects dominate, while above 550 <math>\quad ^{\circ }</math>C elastic effects are negligible. One can also consider the transition zone  where both effects are important (see (Fig. fig:glass_viscosity_a).   In the following we shall consider the temperature dependent properties of a typical soda-lime glass.

Mechanical properties: viscosity and density

The viscosity is the most important property in the glass forming process.  For example, in a typical temperature range encountered in glass forming processes (between approximately 700 and 1200<math>\quad ^{\circ }</math>C) glass viscosity varies from <math display="inline">140 \cdot 10^6</math> Pa s at 700<math>\quad ^{\circ }</math>C to some 400 Pa s. The dependence of glass viscosity upon temperature is typically given by Fulcher expression <span id='citeF-8'></span>[[#cite-8|[8]]]:

<span id="eq-1"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: right;" | 
| <math>log_{10} \mu = A+\frac{B}{T-T_0}  </math>
|}
| style="width: 5px;text-align: right;" | (1)
|}

where <math display="inline">T_0</math>, A and B are constants from experiments. In the present work, the following parameter values are considered (except for example 2): <math display="inline">T_0</math>=220<math>\quad ^{\circ }</math>C, B=4700, A=-2.8. The glass viscosity as a function of temperature is shown in a logarithmic scale on Fig. [[#img-2a|2a]]. Due to very large variations of viscosity during the forming process it is mandatory to  include thermo-mechanical coupling in the model in  order to obtain realistic predictions.

<div id='img-2a'></div>
<div id='img-2'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_García-Espinosa_370675429-visc_total.png|400px|]]
|[[Image:draft_García-Espinosa_370675429-visc_partial.png|400px|Viscosity-temperature curve of soda-lime glass used for beverage containers <span id='citeF-8'></span>[[#cite-8|[8]]].]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 2:''' Viscosity-temperature curve of soda-lime glass used for beverage containers <span id='citeF-8'></span>[[#cite-8|[8]]].
|}

The temperature distribution in the glass is non-uniform at any stage of the glass manufacturing process.  Thus, the viscosity is also non-homogeneous both in space and time. In the context of numerical modeling, Lagrangian formulations are very advantageous when dealing with non-constant properties:  the property (for example, viscosity) is automatically transported being “attached” to the moving nodes. In the Eulerian framework, the representation of each non-constant property evolution would require solving the corresponding transport equation.

Glass density does not undergo considerable changes (2438 kg/m3 at 700<math>\quad ^{\circ }</math>C and 2367 kg/m3 at 1100 <math>\quad ^{\circ }</math>C), thus constant density is an acceptable approximation.

Thermal properties

Heat transfer in the glass is governed not exclusively by conduction, but also by the radiation, which may even be predominant. For strongly absorbing semi-transparent materials this radiation can be  modeled as a diffusion process, thus an effective conductivity taking into account both processes is often defined <span id='citeF-4'></span>[[#cite-4|[4]]]. Real radiation models are complex and computationally expensive. In the scope of this work they are not discussed.

Variation of the specific heat in the temperature range of interest is negligible (1400-1420 J/kg*K between 700 and 1200<math>\quad ^{\circ }</math>C).  The value of the diffusivity <math display="inline">D</math> changes from 0.0000015 <math display="inline">m^2</math>/s at 700<math>\quad ^{\circ }</math>C to 0.0000065 <math display="inline">m^2</math>/s at 1100 <math>\quad ^{\circ }</math>C.

==3 Numerical model==

The numerical model for the complex phenomenon of glass forming developed here consists of a mechanical and thermal modules. These models are coupled in order to take into account the temperature-dependent material properties' evolution.

Kinematic framework

Lagrangian description for modeling the glass forming processes appears attractive as it allows to track the evolution of the deforming domain naturally. The position of the evolving domain coincides with the position of the mesh nodes, defined by the solution of the flow problem. An accurate interface capturing using Eulerian approaches would require much higher mesh resolution to achieve similar precision.  The additional cost associated to the Lagrangian description is due to the necessity of re-meshing the computational domain in order to  avoid excessive element distortion.

The Particle Finite Element Method (PFEM) applied in the present work for modeling the glass forming process is based on the updated Lagrangian description of the governing equations. The key idea of the PFEM is that the variables of interest are stored at the nodes instead of the Gauss points. A finite element mesh is created at every  time step of the transient problem and the solution is then stored at the nodes <span id='citeF-1'></span>[[#cite-1|[1]]]. The nodes are generally maintained (unless adaptive refinement or erasal is performed), thus the mesh re-generation consists in reconnecting the existing nodes. The nodes move to their new position according to their velocity and then the finite element mesh is re-generated using an unconstrained Delaunay triangulation/tesselation <span id='citeF-9'></span>[[#cite-9|[9]]].

In the following the governing equations for the glass are specified. Air is neglected, thus the glass-air interface becomes simply the glass domain boundary. The fully sticking contact with the mold is considered. In Section [[#3.2 Contact modeling and boundary mesh|3.2]] the issues related to the boundary identification and contact treatment are explained in detail.

===3.1 Mechanical model===

At forming temperatures, elastic effects in the glass are negligible and the behavior is nearly iscochoric. Thus, hot glass can be modeled as a viscous incompressible fluid. The total stress tensor can be decomposed into hydrostatic and viscous parts and, therefore, the motion of the hot glass can be described by the Navier-Stokes equations for viscous incompressible fluid.  These can be written for the glass domain <math display="inline">\Omega </math> in Cartesian coordiantes as

<span id="eq-2"></span>
<span id="eq-3"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: right;" | 
| <math>\rho \frac{D \hbox{v} }{D t}+\nabla p-\nabla \cdot (2\mu \epsilon ( \hbox{v} ))=\rho \hbox{g} </math>
|}
| style="width: 5px;text-align: right;" | (2)
|}
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: right;" | 
| <math> \nabla \cdot  \hbox{v} =0  </math>
|}
| style="width: 5px;text-align: right;" | (3)
|}

where <math display="inline"> \hbox{v} </math> is the velocity vector, <math display="inline">p</math> the pressure, <math display="inline">t</math> - time, <math display="inline">\frac{D \hbox{v} }{D t}</math> is the material derivative, <math display="inline">\hbox{g}</math> the '''gravity''', <math display="inline">\rho </math> the density, <math display="inline">\mu </math> the dynamic viscosity  and <math display="inline">\epsilon =\frac{\nabla  \hbox{v} +\nabla ^T \hbox{v} }{2}</math> - the deviatoric strain rate.

''Boundary conditions''. At the mold walls <math display="inline">\Gamma _d</math> and other fixed domain parts (e.g. bottle neck during the gravity stretching), homogeneous Dirichlet boundary conditions are prescribed, i.e.

<span id="eq-4"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: right;" | 
| <math>\hbox{v} = 0 </math>
| style="text-align: right;" | <math> \hbox{ } \hbox{at} \hbox{ } </math>
| <math> \Gamma _d </math>
|}
| style="width: 5px;text-align: right;" | (4)
|}

Two type of Neumann boundary conditions will be distinguished. The first one is the so-called “free-surface” condition that can be approximated for vanishing velocity gradients as: <math display="inline">p=0</math> at <math display="inline">\Gamma _n</math>. This condition will be prescribed at the evolving outer surface of the glass prior to the contact with the mold and at the inner surface prior to the application of the compressed air. A Neumann boundary condition will be used in order to account for the  air pressure <math display="inline">p_a</math> at <math display="inline">\Gamma _n</math>:

<span id="eq-5"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: right;" | 
| <math>\hbox{t}_p = p_a \hbox{n} </math>
| style="text-align: right;" | <math>  \hbox{ } \hbox{at} \hbox{ } </math>
| <math> \Gamma _n </math>
|}
| style="width: 5px;text-align: right;" | (5)
|}

This condition is prescribed at the inner surface of the glass.

====3.1.1 Finite Element formulation====

An equal order linear interpolations for the velocity and pressure variables over 4-noded tetrahedra (3D) is used here for the space discretization of the governing equations Eqs. ([[#eq-2|2]]), ([[#eq-3|3]]). We use Backward Euler time discretization scheme exclusively for the sake of simplicity. All the arguments presented in the paper are valid for any implicit time integration scheme.   Being standard, the space and time discretization are not discussed here (see e.g. <span id='citeF-10'></span>[[#cite-10|[10]]], [Lohner]).  A pressure stabilization term is added due to the use of the equal order velocity-pressure interpolation (see e.g. <span id='citeF-11'></span>[[#cite-11|[11]]] or <span id='citeF-12'></span>[[#cite-12|[12]]]).

Given <math display="inline"> \bar{\hbox{v}} _n</math> and <math display="inline"> \bar{\hbox{p}} _n</math> at <math display="inline">t_n</math>, the time discrete problem consists in finding <math display="inline"> \bar{\hbox{v}} _{n+1}</math>  and <math display="inline"> \bar{\hbox{p}} _{n+1}</math> at <math display="inline">t_{n+1}</math> as the solution of

<span id="eq-6"></span>
<span id="eq-7"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: right;" | 
| <math> \hbox{M}\frac{ \bar{\hbox{v}} _{n+1}- \bar{\hbox{v}} _{n}}{\Delta t}+\mu \hbox{L} \bar{\hbox{v}} _{n+1}+\hbox{G} \bar{\hbox{p}} _{n+1}=\bar{\hbox{F}} </math>
|}
| style="width: 5px;text-align: right;" | (6)
|}
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: right;" | 
| <math> \hbox{D} \bar{\hbox{v}} _{n+1}+ \hbox{S} \bar{\hbox{p}} _{n+1}=0 </math>
|}
| style="width: 5px;text-align: right;" | (7)
|}

where <math display="inline">\hbox{M}</math> is the mass matrix, <math display="inline">\hbox{L}</math> is the Laplacian matrix, <math display="inline">\hbox{G}</math> is the gradient matrix, <math display="inline"> \bar{\hbox{v}} </math> and <math display="inline"> \bar{\hbox{p}} </math> are the velocity and pressure respectively and <math display="inline">\bar{\hbox{F}}</math> is the body force vector.  Note the absence of the convective term due to Lagrangian kinematic framework.

The matrices and vectors are assembled from the element contributions defined as

<span id="eq-8"></span>
<span id="eq-13"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: right;" | 
| <math> \hbox{M}=\rho \int _{\Omega _e}\hbox{N}\hbox{N}^Td\Omega _e</math>
|}
| style="width: 5px;text-align: right;" | (8)
|}
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: right;" | 
| <math> \hbox{L}=\int _{\Omega _e}\nabla \hbox{N} \nabla \hbox{N}^Td\Omega _e</math>
|}
| style="width: 5px;text-align: right;" | (9)
|}
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: right;" | 
| <math> \hbox{G}=-\int _{\Omega _e}\nabla \hbox{N}\hbox{N}d\Omega _e</math>
|}
| style="width: 5px;text-align: right;" | (10)
|}
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: right;" | 
| <math> \bar{\hbox{F}}=\int _{\Omega _e}\hbox{N}\rho \hbox{g}d\Omega _e +\int _{\Gamma _n}p_{a}\hbox{n}d\Gamma _n</math>
|}
| style="width: 5px;text-align: right;" | (11)
|}
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: right;" | 
| <math>\hbox{D}=-\hbox{G}^T</math>
|}
| style="width: 5px;text-align: right;" | (12)
|}
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: right;" | 
| <math> \hbox{S} = \tau \hbox{L}</math>
|}
| style="width: 5px;text-align: right;" | (13)
|}

<math>\hbox{N}</math> stands for standard linear FE shape functions vector, <math display="inline">\Omega _e</math> is the element integration domain, <math display="inline">\tau </math> is an algorithmic stabilization coefficient defined as <math display="inline">\tau =\left(\frac{2|| \bar{\hbox{v}} ||}{h}+\frac{4\nu }{h^{2}}\right)^{-1}</math>, where <math display="inline">h</math> is the element size  <span id='citeF-11'></span>[[#cite-11|[11]]].  Note also that the discrete operators given by Eqs. [[#eq-8|8]]-[[#eq-13|13]] correspond to the unknown current configuration <math display="inline">\hbox{X}_{n+1}</math> according to the updated Lagrangian approach <span id='citeF-1'></span>[[#cite-1|[1]]], <span id='citeF-13'></span>[[#cite-13|[13]]].  Thus, the system is non-linear and must be solved in an iterative manner by updating the discrete operators at every non-linear iteration.

=====3.1.1.1 Modified fractional step technique for the efficient solution of governing equations=====

Solving the governing equation system monolithically (i.e. for velocities and pressure simultaneously) has a considerable computational cost <span id='citeF-14'></span>[[#cite-14|[14]]].  Pressure segregation or fractional step methods are known for their high computational efficiency, however they often lead to mass conservation problems <span id='citeF-15'></span>[[#cite-15|[15]]].  We propose to use the modified version of the fractional step approach presented in <span id='citeF-2'></span>[[#cite-2|[2]]] and further developed in <span id='citeF-3'></span>[[#cite-3|[3]]].  On the one hand, it inherits the high computational efficiency of  the original fractional step procedure ( <span id='citeF-16'></span>[[#cite-16|[16]]],  <span id='citeF-17'></span>[[#cite-17|[17]]], <span id='citeF-18'></span>[[#cite-18|[18]]]) due to the decoupling of the velocity and the pressure. On the other hand, it has much better mass conservation properties.

The fractional step split is applied here at purely algebraic level to the governing equations system defined by Eqs. ([[#eq-6|6]])-([[#eq-7|7]]). The momentum equation is split into two parts by introducing the intermediate velocity <math display="inline">\tilde{ \hbox{v} }</math> and the original monolithic system  is replaced by

<span id="eq-14"></span>
<span id="eq-15"></span>
<span id="eq-16"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: right;" | 
| <math> \hbox{M}\frac{\tilde{ \hbox{v} }- \bar{\hbox{v}} _{n}}{\Delta t}+\mu \hbox{L} \bar{\hbox{v}} _{n+1}+\hbox{G} \bar{\hbox{p}} _{n+1}^{g}=\bar{\hbox{F}}</math>
|}
| style="width: 5px;text-align: right;" | (14)
|}
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: right;" | 
| <math> \hbox{M}\frac{ \bar{\hbox{v}} _{n+1}-\tilde{ \hbox{v} }}{\Delta t} +\hbox{G}\left( \bar{\hbox{p}} _{n+1}- \bar{\hbox{p}} _{n+1}^{g}\right)=0</math>
|}
| style="width: 5px;text-align: right;" | (15)
|}
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: right;" | 
| <math> \hbox{D} \bar{\hbox{v}} _{n+1}+ \hbox{S} \bar{\hbox{p}} _{n+1}=0</math>
|}
| style="width: 5px;text-align: right;" | (16)
|}

where <math display="inline">\tilde{ \hbox{v} }</math> is an auxiliary vector, representing the intermediate or “fractional” nodal velocities and <math display="inline"> \bar{\hbox{p}} _{n+1}^g</math>  is the guess of the end-of-step nodal pressure vector.  Eq. ([[#eq-14|14]]) is known as “fractional momentum” and  Eq. ([[#eq-15|15]]) as “end-of-step momentum” equations. The novelty of <span id='citeF-2'></span>[[#cite-2|[2]]] in contrast to the classical fractional step approach consists in using <math display="inline">  \bar{\hbox{p}} _{n+1}^g</math> instead of <math display="inline"> \bar{\hbox{p}} _{n}</math> in the fractional momentum equation.

The pressure Poisson's equation is obtained by enforcing the incompressibility condition (Eq. ([[#eq-16|16]])) with the end-of-step momentum equation (Eq. ([[#eq-15|15]])), leading to

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\hbox{D}\tilde{ \hbox{v} }= \Delta t\hbox{D}\hbox{M}^{-1}\hbox{G}\left( \bar{\hbox{p}} _{n+1}- \bar{\hbox{p}} _{n+1}^g\right)+ \hbox{S} \bar{\hbox{p}} _{n+1} </math>
|}
| style="width: 5px;text-align: right;" | (17)
|}

Using the typical approximation <math display="inline">\hbox{D}\hbox{M}^{-1}\hbox{G}\approx \hbox{L}</math>, we arrive at the final system of discretized equations to be solved:

<span id="eq-18"></span>
<span id="eq-19"></span>
<span id="eq-20"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\left(\hbox{M}\frac{\tilde{ \hbox{v} }- \bar{\hbox{v}} _{n}}{\Delta t}+\mu \hbox{L}\tilde{ \hbox{v} } +\hbox{G} \bar{\hbox{p}} _{n+1}^g\right)=\bar{\hbox{F}}  </math>
|}
| style="width: 5px;text-align: right;" | (18)
|}
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math> \hbox{D}\tilde{ \hbox{v} }= \Delta t\hbox{L}\left( \bar{\hbox{p}} _{n+1}- \bar{\hbox{p}} _{n+1}^g\right)+\hbox{S} \bar{\hbox{p}} _{n+1} </math>
|}
| style="width: 5px;text-align: right;" | (19)
|}
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math> \hbox{M}\frac{ \bar{\hbox{v}} _{n+1}-\tilde{ \hbox{v} }}{\Delta t} +\hbox{G}\left( \bar{\hbox{p}} _{n+1}- \bar{\hbox{p}} _{n+1}^g\right)=0  </math>
|}
| style="width: 5px;text-align: right;" | (20)
|}

=====3.1.1.2 Fractional momentum equation solution and pressure prediction=====

The momentum equation (Eq. ([[#eq-18|18]])) is non-linear due to the dependence of the discrete operators on the unknown current configuration <math display="inline">\hbox{X}_{n+1}</math>. Thus, it must be solved iteratively. For this reason, let us define the residual of the fractional momentum equation as

<span id="eq-21"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: right;" | <math></math>
| <math>\tilde{\hbox{r}}_m=\bar{\hbox{F}}-\left(\hbox{M}\frac{\tilde{ \hbox{v} }- \bar{\hbox{v}} _{n}}{\Delta t}+\mu \hbox{L} \bar{\hbox{v}} _{n+1}+\hbox{G} \bar{\hbox{p}} _{n+1}^{g}\right) </math>
|}
| style="width: 5px;text-align: right;" | (21)
|}

The modified fractional step method proposed in  <span id='citeF-2'></span>[[#cite-2|[2]]] consists in considering the pressure variation in the fractional momentum equation. Thus, <math display="inline"> \bar{\hbox{p}} _{n+1}^g</math> cannot be considered constant in the residual expression. In order to obtain the fractional momentum equation dependent on the velocity exclusively (and thus maintain the decoupling of the velocities and the pressure), a prediction for the pressure <math display="inline">p_{n+1}^g</math> can be computed using an assumption of slight compressibility. Following this assumption the current-step pressure is obtained by adding the term proportional to the divergence of velocity  to the pressure of the previous step (see <span id='citeF-13'></span>[[#cite-13|[13]]], <span id='citeF-2'></span>[[#cite-2|[2]]] for details):

<span id="eq-22"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: right;" | 
| <math>p^g_{n+1}=p_{n}+\kappa \int _{t_n}^{t_{n+1}}\nabla \cdot  \hbox{v} dt  </math>
|}
| style="width: 5px;text-align: right;" | (22)
|}

where <math display="inline">\kappa </math> is the compressibility parameter of the fluid. The discrete form of Eq. ([[#eq-22|22]]) using linear velocity-pressure finite elements reads

<span id="eq-23"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: right;" | 
| <math>\hbox{M}_p \bar{\hbox{p}} _{n+1}^g=\hbox{M}_p \bar{\hbox{p}} _{n}+\kappa \int _{t_n}^{t_{n+1}}\hbox{D} \bar{\hbox{v}} dt  </math>
|}
| style="width: 5px;text-align: right;" | (23)
|}

where <math display="inline">\hbox{M}_p</math> is the pressure mass matrix.

In order to avoid matrix inversion for obtaining the current step pressure, the pressure mass matrix <math display="inline">\hbox{M}_p</math> will be taken in the lumped form.  Performing the integration, Eq. ([[#eq-23|23]]) can be re-written as

<span id="eq-24"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: right;" | 
| <math>\bar{\hbox{p}} _{n+1}^g= \bar{\hbox{p}} _{n}+\kappa \Delta t\hbox{M}_{p}^{-1}\hbox{D} \bar{\hbox{v}} _{n+1}  </math>
|}
| style="width: 5px;text-align: right;" | (24)
|}

Expressing the pressure in terms of velocity according to Eq. ([[#eq-24|24]]) allows to define the iterative solution of the non-linear equation <math display="inline">\tilde{\hbox{r}}_m=0</math> (with <math display="inline">\tilde{\hbox{r}}_m</math>  defined by Eq. ([[#eq-21|21]])) exclusively in terms of the nodal velocities:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: right;" | 
| <math>\hbox{H}\delta  \bar{\hbox{v}} =\tilde{\hbox{r}}_m( \bar{\hbox{v}} _i,  \bar{\hbox{p}} _i)</math>
|}
| style="width: 5px;text-align: right;" | (25)
|}
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: right;" | 
| <math> \bar{\hbox{v}} ^{i+1}= \bar{\hbox{v}} ^{i}+\delta  \bar{\hbox{v}} </math>
|}
| style="width: 5px;text-align: right;" | (26)
|}
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: right;" | 
| <math>\hbox{X}^{i+1}=\hbox{X}^{i}+\Delta t \bar{\hbox{v}} ^{i+1} </math>
|}
| style="width: 5px;text-align: right;" | (27)
|}

where “i” stands for the non-linear iteration index at time <math display="inline">t_{n+1}</math> and <math display="inline">\hbox{H}</math> is the tangent matrix defined as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: right;" | 
| <math>\hbox{H}=-\frac{\partial \tilde{\hbox{r}}_{m}}{\partial  \bar{\hbox{v}} } </math>
|}
| style="width: 5px;text-align: right;" | (28)
|}

As the non-constant pressure term is now included in the residual, it must be accounted for in the linearization.

According to <span id='citeF-3'></span>[[#cite-3|[3]]] the dynamic tangent matrix can be approximated as:

<span id="eq-29"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: right;" | 
| <math>\hbox{H}=\frac{\hbox{M}}{\Delta t}+\mu \hbox{L}+\Delta t \int _{\Omega _e}\hbox{B}^T\hbox{C}_{\kappa } \hbox{B} d \Omega _e  </math>
|}
| style="width: 5px;text-align: right;" | (29)
|}

where the operator <math display="inline">\hbox{B}</math> and the volumetric constitutive matrix <math display="inline">\hbox{C}_{\kappa }</math> are defined (in 2D) as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\hbox{B}=\begin{pmatrix}\frac{\partial N_1}{\partial x}       &  0             & \frac{\partial N_2}{\partial x}    &      0        &   \frac{\partial N_3}{\partial x}   &           0\\       0                   & \frac{\partial N_1}{\partial y}  & 0                & \frac{\partial N_2}{\partial y} &            0      & \frac{\partial N_3}{\partial y}\\       \frac{\partial N_1}{\partial y}       & \frac{\partial N_1}{\partial x}  & \frac{\partial N_2}{\partial y}    & \frac{\partial N_2}{\partial x} & \frac{\partial N_3}{\partial y}    & \frac{\partial N_3}{\partial x}  \end{pmatrix}       </math>
|}
| style="width: 5px;text-align: right;" | (30)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\hbox{C}_{\kappa }=\frac{1}{2}\begin{pmatrix}\kappa & \kappa  & 0\\      \kappa & \kappa  & 0\\       0 & 0 & 0    \end{pmatrix}  </math>
|}
| style="width: 5px;text-align: right;" | (31)
|}

Elimination of the unknown pressure variables from the left-hand-side of the modified fractional momentum equation enables us to solve it for the velocities only (as it is done in the standard fractional step), thereby preserving the decoupling of  the nodal velocities and pressures. However at every non-linear iteration the nodal pressure needs to be updated using Eq. ([[#eq-24|24]]).

=====3.1.1.3 Pressure Poisson's equation and the correction step=====

The next step to be carried out is the correction of the pressure, i.e. obtaining the end-of-step incompressible pressure field using Eq. ([[#eq-19|19]]).  Solution of Eq. ([[#eq-19|19]]) requires to impose the pressure boundary conditions due to the presence of the Laplacian <math display="inline">\hbox{L}</math>. According to the methodology presented in <span id='citeF-2'></span>[[#cite-2|[2]]], <math display="inline"> \bar{\hbox{p}} _{n+1}= \bar{\hbox{p}} _{n+1}^g</math>  can be used as an essential boundary condition for the pressure necessary for solving the Poisson's equation for the pressure.  The quality of this approximation is related exclusively to the value of the compressibility constant <math display="inline">\kappa </math> used in the prediction step.  Having the pressure fixed to the predicted value <math display="inline"> \bar{\hbox{p}} _{n+1}^g</math> at the free surface, the pressure Poisson's equation is solved elsewhere in the domain to yield the end-of-step pressure vector <math display="inline"> \bar{\hbox{p}} _{n+1}</math>.

This step can be thus viewed as a correction of the predicted pressure <math display="inline"> \bar{\hbox{p}} _{n+1}^g</math> to the correct end-of-step one everywhere except for the free surface, where the “slightly compressible” pressure is maintained.  Consequently, the projection step is carried out according to Eq. ([[#eq-20|20]]) and returns the end-of-step divergence-free velocity everywhere in the domain except for the pressure boundary, where the divergence-free velocity is approximated.

The implementation of the modified fractional step scheme is very similar to that of the classical fractional step method. It is summarized next.

<ol>

<li>Solve the fractional momentum, Eq. ([[#eq-18|18]])</li>
* Update the nodal pressure vector according to Eq. ([[#eq-24|24]]) and add it to the fractional momentum residual
* Update the discrete operators according to the new nodal positions <math display="inline">\hbox{X}_{n+1}</math>
* Repeat two previous steps until convergence for the velocity field is achieved

<li>Solve the pressure Poisson's Eq. ([[#eq-19|19]]), using the predicted pressure as the boundary condition</li> 
<li>Solve the end-of-step momentum equation (Eq. ([[#eq-20|20]]))</li>

</ol>

===3.2 Contact modeling and boundary mesh===

====3.2.1 Contact modeling====

Accurate modeling of contact between the glass and the wall of the mold is essential for predicting the wall thickness of the glass object correctly. Literature exhibits a large number of numerical methods. In Eulerian non-matching grid formulations a large variety of methods has been proposed ranging from penalty approaches (<span id='citeF-19'></span>[[#cite-19|[19]]], <span id='citeF-20'></span>[[#cite-20|[20]]]) to augmented Lagrange multipliers <span id='citeF-21'></span>[[#cite-21|[21]]].  For the matching grid methods the contact modeling is usually done using simple geometric considerations <span id='citeF-4'></span>[[#cite-4|[4]]].  In this study we assume that the mold is a rigid body and the glass sticks to the mold. This is a reasonable and commonly accepted approximation (see e.g. <span id='citeF-6'></span>[[#cite-6|[6]]], <span id='citeF-4'></span>[[#cite-4|[4]]]).

We propose here a simple algorithm that allows preserving the overall strategy of the PFEM, but at the same time does not lead to the incorrect prediction of the wall thickness. For the sake of minimizing the computational cost associated with re-meshing, PFEM utilizes unconstrained Delaunay triangulation/tetrahedrization ([tetgen]) equipped with the alpha-shape technique for detecting boundaries.

====3.2.2 Modified PFEM contact algorithm====

Considering that the nodes follow a variable distribution <math display="inline">h(x)</math>, which is the minimum distance between two nodes, the ''alpha-shape'' technique (see <span id='citeF-22'></span>[[#cite-22|[22]]], <span id='citeF-1'></span>[[#cite-1|[1]]]) applied to the  unconstrained Delaunay mesh allows to  distinguish whether an element defined by the four nodes must be created or not. The radius <math display="inline">r</math> of a sphere defined by the examined nodes is compared to the  corresponding size <math display="inline">h</math> (average size of the inter-nodal distances). If the ratio <math display="inline">\frac{r}{h}>\alpha </math>, where <math display="inline">\alpha </math> is the alpha-shape parameter typically taken as 1.5 in 3D, the element is not created and the nodes are marked as belonging to the boundary (for details on the alpha-shape techniques one can see <span id='citeF-1'></span>[[#cite-1|[1]]]).  For the correct contact representation we propose to combine the alpha-shape technique with the boundary markers.

Let us distinguish the outer boundary of the glass domain and the mold wall with an ''interface''  and ''wall'' flags, respectively.  As the ''interface'' nodes approach the ''wall'' nodes, the distance between them diminishes and an element is created.  According to the standard PFEM algorithm such an element immediately obtains the properties of the fluid, thus bringiung the fluid and the  solid into contact (see Fig. [[#img-3|3]]).

<div id='img-3a'></div>
<div id='img-3'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_García-Espinosa_370675429-contact_0_new.png|400px|]]
|[[Image:draft_García-Espinosa_370675429-contact_1_new_PFEM.png|400px|Schematic representation of the standard PFEM contact algorithm]]
|- style="text-align: center; font-size: 100%;"
| (a) Glass and wall are separated
| (b) Fluid element is created
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 3:''' Schematic representation of the standard PFEM contact algorithm
|}

However, this newly created element (see Fig. [[#img-3a|3a]]) defines a “premature” contact (both mechanical and thermal) between the glass and the wall, as the distance between the original fluid boundary nodes and the wall are still of order <math display="inline">h</math>. This implies introducing an error of order <math display="inline">h</math> in the determination of the thickness of the fluid layer (glass wall in our case). In this work we propose a simple change in the contact algorithm that leads to a considerable improvement in the accuracy of thickness determination.

We propose to consider such the newly created contact element a purely geometrical entity, that does not contribute to the governing equation system and thus does not result in any resistance to the fluid motion.  This can be formally described as: if the mesher creates an element that contains both ''interface'' and ''wall'' nodes, mark it as a ''fictitious element'' and multiply its  contribution (residual and tangent matrix) by zero during the finite element assembly process. Nevertheless, fictitious elements can be used for tracking the distance between the ''interface'' nodes and the wall.  Moreover, creating fictitious elements as geometric entities enables maintaining the original PFEM unconstrained meshing strategy.

<div id='img-4'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_García-Espinosa_370675429-contact_0_new.png|400px|]]
|[[Image:draft_García-Espinosa_370675429-contact1_new.png|400px|]]
|- style="text-align: center; font-size: 100%;"
| (a) Glass and wall are separated
| (b) Fictitious contact elements are created
|-
|[[Image:draft_García-Espinosa_370675429-contact2_new.png|400px|]]
|[[Image:draft_García-Espinosa_370675429-contact3_new.png|400px|Schematic representation of the modified PFEM contact algorithm]]
|- style="text-align: center; font-size: 100%;"
| (c) Removal of interface nodes close to the wall
| (d) Stick contact between the glass and the wall
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 4:''' Schematic representation of the modified PFEM contact algorithm
|}

The distance <math display="inline">d</math> between the ''interface'' node and the wall is computed within each fictitious element. As long as <math display="inline">d</math> is larger than a given tolerance <math display="inline">\epsilon </math>, the fictitious element is maintained and  no actual contact between the fluid and the solid is established. As soon as <math display="inline">d</math> reached a value lower than a given tolerance <math display="inline">\epsilon </math>, the ''interface'' node is erased, thus the glass domain gets connected to the wall via a real element in a fully sticking contact (see Fig. [[#img-4|4]] c) and d)).  This way, the error in the determination of the thickness of the bottle becomes controlled by the user-defined tolerance <math display="inline">\epsilon </math>. In the simulations used in the present work, <math display="inline">\epsilon </math> was taken as 10 <math display="inline">%</math> of the element size.

====3.2.3 Preserving the boundary mesh quality====

The domain deformation occurring in a glass forming process can be viewed as a stretch in different directions. This leads to the danger of decreasing the mesh resolution and possibility of creating non-physical “holes” at the glass domain boundary due to alpha-shape (see stretching of a domain fragment fixed at the top on Fig. fig:mesh_ref_a). In order to preserve the mesh quality and avoid the creation of holes, an adaptive refinement of glass domain boundary is used. A new node is introduced at the center of the boundary elements edge as the corresponding length reaches of 1.5 times the original value (see Fig. fig:mesh_ref_b).  <div id='img-5'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_García-Espinosa_370675429-ref_wrong.png|400px|]]
|[[Image:draft_García-Espinosa_370675429-ref_right.png|400px|Mesh refinement at the boundary]]
|- style="text-align: center; font-size: 100%;"
| (a) Mesh deterioration
| (b) Mesh-preserving refinement
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 5:''' Mesh refinement at the boundary
|}

===3.3 Thermal model===

During glass forming process the temperature distribution within the glass is highly non-uniform. First of all, the temperature varies in the thickness direction of the glass domain due to cooling of the free surfaces. Second, when the glass  reaches the mold wall a local rapid cooling occurs, resulting in a strong increment of viscosity, thus altering the shape evolution: zones with the higher temperature will  undergo stronger stretching and thinning <span id='citeF-23'></span>[[#cite-23|[23]]]. Thus, inclusion of the heat equation is mandatory and purely mechanical models ignoring the heat transfer must be discarded when modeling the problem of interest.

Assuming that the conduction obeys the Fourier law which relates the heat flux with the temperature and the thermal conductivity, the heat transfer equation can be written in the Lagrangian reference frame as:

<span id="eq-32"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\rho c \frac{D  T }{D t}  = k\nabla ^2 T + q  </math>
|}
| style="width: 5px;text-align: right;" | (32)
|}

where <math display="inline">c</math> is the specific heat, <math display="inline">T</math> is the temperature, <math display="inline">k</math> is the thermal conductivity and <math display="inline">q</math> is the internally generated heat flux. Using linear temperature finite elements and the Backward Euler scheme for space and time discretization, respectively, leads to the following equation:

<span id="eq-33"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\rho c\frac{\bar{\hbox{T}}_{n+1}-\bar{\hbox{T}}_n}{\Delta t} \hbox{M}  = k\hbox{L} \bar{\hbox{T}}_{n+1} + \hbox{Q}  </math>
|}
| style="width: 5px;text-align: right;" | (33)
|}

where <math display="inline">\hbox{M}</math> and <math display="inline">\hbox{L}</math> are the mass and Laplacian matrices and <math display="inline">\hbox{Q}=\int _{\Omega _e}\hbox{N}q d\Omega _e </math>. Note the absence of the convective term in the Lagrangian framework. Heat convection is automatic in the Lagrangian model, as the heat is convected by the moving Lagrangian grid “attached to the material”.

The problem must be equipped with the Dirichlet boundary conditions (fixed temperature <math display="inline">T=T_d</math> at mold in the locations to be specified) and Neumann boundary condition (adiabatic condition will be considered in this work. Heat conduction prior to the contact with the mold (occurring due to convection of the air in contact with the glass) will not be considered here. Once the contact with the mold has taken place, the direct heat conduction between the glass and the mold takes place. Temperature field <math display="inline">T_{n+1}</math> obtained as a solution of Eq. ([[#eq-33|33]]) is used for computing the viscosity distribution (Eq. ([[#eq-1|1]])) for the mechanical model.  More details on the thermo-mechanical coupling can be found in <span id='citeF-24'></span>[[#cite-24|[24]]], <span id='citeF-25'></span>[[#cite-25|[25]]].

===3.4 Overall solution algorithm===

To this end all the “ingredients” of the thermo-mechanical model are defined. The overall solution algorithm for the bottle forming process is presented next.

<ol>

<li>Discretize the glass domain with a finite element mesh.    </li>
<li>Identify the external boundaries for the fluid and glass-mold wall contact elements (using the alpha-shape technique + boundary markers)       </li>
<li>Solve the Lagrangian equations for the glass. Obtain nodal velocities, pressure and displacements.       </li>
<li>Move the mesh nodes to a new position (according to the computed displacements).    </li>
<li>Solve the heat transfer equation. Obtain the temperature distribution of the glass.    </li>
<li>Update viscosity according to the obtained temperature distribution using Eq. [[#eq-1|1]]    </li>
<li>Re-generate the mesh for the glass domain (nodes are generally maintained)    </li>
<li>Go back to the next time step. Start the solution from Step 2          </li>

</ol>

==4 Examples==

In this section two numerical examples are solved.  We mentioned previously that there exist practically no established benchmark for the glass forming models. In the majority of literature on the numerical modeling of the glass manufacturing process the data provided is insufficient for reproducing the tests. Thus, in this section we first validate the proposed formulation by solving  one of the very few examples suitable for the validation of the glass forming simulation model on a simple rectangular geometry. Next we propose a numerical test dealing with the bottle-forming that could serve as a reference in the future. The main aim is to accurately define all the data necessery for reporducibility of this example.

===4.1 TV panel pressing===

This example reported by Berndhauser in <span id='citeF-26'></span>[[#cite-26|[26]]] models the manufacturing of a TV glass panel. Even though it does not deal explicitly with the bottle manufacturing process, the characteristics of the example are very similar (pressing the glass gob into a mold, accounting for thermo-mechanical nature of the process). Moreover, it provides simple nearly rectangular geometry which facilitates the validation. A portion of molten glass in the shape of cylindrical gob, is pressed into the mold by the plunger (see Fig [[#img-6|6]]).  The plunger is moving downwards with a prescribed time-dependent velocity. During pressing, the glass is flowing while it cools down because of the contact with the colder mold and plunger. The problem settings and material properties presented below are taken from <span id='citeF-26'></span>[[#cite-26|[26]]].

<div id='img-6'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_García-Espinosa_370675429-TV_panel_model.png|600px|Cross section of the glass pressing process: initial and final states]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 6:''' Cross section of the glass pressing process: initial and final states
|}

The problem to be simulated accounts for the thermo-mechanical process involving the plunger, the glass and the mold. The mold and the plunger are considered to be rigid. The overall time of 3 seconds is simulated. The plunger descends with the vertical velocity <math display="inline">v_z=(-0.0842 e^{-1.535 t}+0.00842)</math> m/s for 0.0s <math display="inline"><</math> t <math display="inline"><</math> 1.5 s. For 1.5s <math display="inline"><</math> t <math display="inline"><</math> 3.0s it remains still. The no-slip condition at the glass-wall interface is considered.

The thermal conditions are as follows: the initial glass temperature <math display="inline">T_g=1000</math> <math>\quad ^{\circ }</math>C. The initial temperature of the plunger and the mold is <math display="inline">T_p=T_m</math>=500 <math>\quad ^{\circ }</math>C. Dirichlet boundary condition for the temperature is prescribed at the outer surface of the mold and the outer surface of the plunger.  The rest of the mold and the plunger surfaces are considered to be adiabatic (heat flux <math display="inline">Q=0</math>).

The geometry of the model is shown in Fig. [[#img-7|7]]. Geometrical details can be found in Fig. [[#img-8|8]]. The initial distance between the mold and the plunger <math display="inline">d_{pl}^{0}=4.67376</math> cm is equal to the original height of the glass gob. The gob has a cylindrical shape with a radius of <math display="inline">r=13.243</math> cm.

Material properties of the plunger, the mold and the glass are summarized in Table [[#table-3|3]]. Gravity of 9.8 <math display="inline">m/s^2</math> is considered. <div id='img-7'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_García-Espinosa_370675429-screen_model.png|600px|TV screen forming set-up (one quarter): mold, gob and plunger]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 7:''' TV screen forming set-up (one quarter): mold, gob and plunger
|}

<div id='img-8a'></div>
<div id='img-8'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_García-Espinosa_370675429-plunger.png|400px|]]
|- (a) Cross section of the plunger; front and side view
|-
|[[Image:draft_García-Espinosa_370675429-mould.png|400px|Dimensions of the mould and the plunger.]]
|- (a) Cross section of the mold; front and side view
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 8:''' Dimensions of the mould and the plunger.
|}


{| class="wikitable" style="text-align: center; margin: 1em auto;"
|+ <span id='table-3'></span>Table. 3 Material properties of the model components

|- style="border-top: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" colspan="2" | Glass
|- style="border-top: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" | Viscosity <math display="inline">\mu </math> 
| style="border-left: 2px solid;border-right: 2px solid;" | <math>10^{-2.8+\frac{4700}{T-220}}</math>Pa s 
|-
| style="border-left: 2px solid;border-right: 2px solid;" | Conductivity <math display="inline">\lambda _g</math> 
| style="border-left: 2px solid;border-right: 2px solid;" | 5 W/m K 
|-
| style="border-left: 2px solid;border-right: 2px solid;" | Specific heat <math display="inline">c_g</math> 
| style="border-left: 2px solid;border-right: 2px solid;" | 1400 J/kg K 
|-
| style="border-left: 2px solid;border-right: 2px solid;" | Density <math display="inline">\rho _g</math> 
| style="border-left: 2px solid;border-right: 2px solid;" | 2500 kg/<math display="inline">m^3</math> 
|- style="border-top: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" colspan="2" | Plunger and mold
|- style="border-top: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" |   Conductivity <math display="inline">\lambda _m</math> 
| style="border-left: 2px solid;border-right: 2px solid;" | 20 W/m K 
|-
| style="border-left: 2px solid;border-right: 2px solid;" | Specific heat <math display="inline">c_m</math> 
| style="border-left: 2px solid;border-right: 2px solid;" | 500 J/kg K 
|- style="border-bottom: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" | Density <math display="inline">\rho _m</math> 
| style="border-left: 2px solid;border-right: 2px solid;" | 8000 kg/<math display="inline">m^3</math> 

|}

====4.1.1 Results====

According to the benchmark proposal, the temperature, pressure and the glass front evolution are analyzed. The values of the variables of interest are inspected at several locations indicated in Fig. [[#img-9|9]]. “Sensor” lines are located at four corners of the plain part of the mold bottom as L1: x=0, y=0; L2: x=0.22 m y=0; L3: x=0, y=0.16 m; L4: x=0.22, y=0.16 m. <div id='img-9'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_García-Espinosa_370675429-Monitor_Lines_Pavel.png|600px|Location of the inspection lines (“sensors”) <span id='citeF-26'></span>[[#cite-26|[26]]]]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 9:''' Location of the inspection lines (“sensors”) <span id='citeF-26'></span>[[#cite-26|[26]]]
|}

Fig. [[#img-10|10]] displays the evolution of the temperature at L1 reported in <span id='citeF-26'></span>[[#cite-26|[26]]]. One can see the compressing gob due to the plunger descend. In the beginning the temperature field exhibits a sharp  discontinuity at the glass/mold and the glass/plunger interface. As time evolves, the transition region appears. At t=1.5 s the plunger stops and the process is governed by diffusion only.

<div id='img-10'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 70%;max-width: 100%;"
|-
|[[Image:draft_García-Espinosa_370675429-reference_res_clean_temperature.png|600px|Temperature evolution at L1 according to <span id='citeF-26'></span>[[#cite-26|[26]]]]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 10:''' Temperature evolution at L1 according to <span id='citeF-26'></span>[[#cite-26|[26]]]
|}

A comparison of the results obtained in the present work with the reference data is displayed in Fig. [[#img-11|11]]. One can see a good agreement between the compared values.  In the beginning the temperature is distributed almost uniformly across the glass thickness (z-direction). Around 0.75 s the maximum temperature concentrates in the middle plain of the gob. At the end of  the simulation the maximum temperature in the middle of the glass thickness reduces from 1000 <math>\quad ^{\circ }</math>C to some 950 <math>\quad ^{\circ }</math>C.

<div id='img-11a'></div>
<div id='img-11'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_García-Espinosa_370675429-temper_comparison_00s.png|400px|]]
|[[Image:draft_García-Espinosa_370675429-temper_comparison_0-25s.png|400px|]]
|-
|[[Image:draft_García-Espinosa_370675429-temper_comparison_0-5s.png|400px|]]
|[[Image:draft_García-Espinosa_370675429-temper_comparison_0-75s.png|400px|]]
|-
|[[Image:draft_García-Espinosa_370675429-temper_comparison_1-25s.png|400px|]]
|[[Image:draft_García-Espinosa_370675429-temper_comparison_3s.png|400px|Temperature along L1 at different time steps. Comparison with <span id='citeF-26'></span>[[#cite-26|[26]]]]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 11:''' Temperature along L1 at different time steps. Comparison with <span id='citeF-26'></span>[[#cite-26|[26]]]
|}

The slightly smaller temperature observed at t=0.01 s in the lowest edge of the moving plunger (z=0.075) is a spurious “overshoot”, likely to occur due to the steep temperature gradient  (this phenomenon is addressed, for example, <span id='citeF-27'></span>[[#cite-27|[27]]]). Steep temperature gradient is present at both contact interfaces, namely at z=0.03 and z=0.075.  However, since mold and the adjacent glass at the bottom (z=0.03) are not moving, no convection occurs there and thus the overshoot is observed only at the contact above (z=0.075), where glass is  in contact with the moving plunger. We note that the overshoot is only of some 3 ºC (which is less than 1 <math display="inline">%</math> of the temperature difference between the mold and the glass).

Evolution of glass front in the vertical direction of the xz-cross section (long side, x=0.23 m, y=0) is shown in Fig. [[#img-13|13]]. In the reference <span id='citeF-26'></span>[[#cite-26|[26]]] the location of the glass front  is reported only for t<math display="inline">></math>1 s. We included the results spanning over the total simulation time. Since the plunger reaches zero velocity at t=1.5 s, the glass front position remains constant afterwards. One can see good agreement between the results of the present work and the ones obtained in <span id='citeF-26'></span>[[#cite-26|[26]]].  Slightly slower increment in height may be attributed to either different mesh resolutions or presence of the gravity in the present work (gravity was neglected in <span id='citeF-26'></span>[[#cite-26|[26]]]). One other factor that may influence this result is the volume conservation of the method. Spurious volume gain would lead to the faster domain evolution, while volume losses would lead to a slower growth. In order to exclude this source of inaccuracy of our method, the glass volume evolution had been recorded throughout the simulation (see Fig. [[#img-12|12]]). One can see excellent volume conservation features of the present method. Maximum volume variation observed was below 4 <math display="inline">%</math>. Volume conservation features of the method are further analyzed in the next example.              <div id='img-12'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_García-Espinosa_370675429-pressure_0-5_clean.png|400px|]]
|[[Image:draft_García-Espinosa_370675429-volume_evol_screen.png|400px|TV screen pressing: pressure and volume data.]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 12:''' TV screen pressing: pressure and volume data.
|}

The pressure evolution at four inspection lines is shown in Fig. fig:pressure_sensors.  It is important to note that pressure varies insignificantly in the thickness direction. Thus a single value is representative. Originally, glass does not reach the position of sensors L2, L3 and L4 (thus, zero values are observed). First, the glass front reaches L2, then L3 and, finally, the further most corner L4. The maximum pressure of appx. 1.1 MPa  is reached at L1 (t=1.25 s). An example of the pressure distribution at t=0.50 s is shown in Fig. fig:pressure_example.

<div id='img-13'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_García-Espinosa_370675429-4_press_sensors.png|400px|]]
|[[Image:draft_García-Espinosa_370675429-z_coord_evol.png|400px|Pressure and glass front evolution in the gob]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 13:''' Pressure and glass front evolution in the gob
|}

===4.2 Bottle manufacturing: final blow process===

In this example  we concentrate our attention on the second stage of the bottle manufacturing process, namely the gravity stretching and the final blow.  The objective of this example is to test the capability of the presented numerical  approach to solving this challenging problem.  Moreover, we strive to establish an example that may be used as a first approximation to a benchmark for the thermo-mechanical simulations of the bottle forming process.  In order to facilitate the modeling and reduce computational time, the mold is modeled merely as a set of nodes discretizing its inner wall providing Dirichlet boundary for both the mechanical and the thermal problem. The thermal problem is modeled considering convective and diffusive heat transfer modes (radiation is neglected). We note again that adopting the Lagrangian framework, convection in the glass phase is automatically solved due to the mesh motion. No convection condition at the glass domain boundary is prescribed.

The initial geometry consists of the mold wall and the parison. The geometry data is shown in Fig. fig:final_blow_bench_model_c.


{| class="wikitable" style="text-align: center; margin: 1em auto;"
|+ <span id='table-4'></span>Table. 4 Material properties of glass in the final blow simulation

|- style="border-top: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" | Viscosity 
| style="border-left: 2px solid;border-right: 2px solid;" | <math>\mu=265677693762693 e^{-0.0233569026*T}</math>Pa s 
|-
| style="border-left: 2px solid;border-right: 2px solid;" | Conductivity <math display="inline">\lambda _g</math> 
| style="border-left: 2px solid;border-right: 2px solid;" | 1.5 W/m K 
|-
| style="border-left: 2px solid;border-right: 2px solid;" | Specific heat <math display="inline">c_g</math> 
| style="border-left: 2px solid;border-right: 2px solid;" | 1409 J/kg K 
|- style="border-bottom: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" | Density <math display="inline">\rho _g</math> 
| style="border-left: 2px solid;border-right: 2px solid;" | 2400 kg/<math display="inline">m^3</math> 

|}

.

<div id='img-14a'></div>
<div id='img-14'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 75%;max-width: 100%;"
|-
|[[Image:draft_García-Espinosa_370675429-bot_model_bench.png|400px|]]
|
|[[Image:draft_García-Espinosa_370675429-bot_model_bench_points.png|400px|]]
|
|[[Image:draft_García-Espinosa_370675429-bot_model_bench_3d.png|400px|]]
|
|[[Image:draft_García-Espinosa_370675429-bottle_initial_temperature1.png|400px|Model for the final blow simulation]]
|- style="text-align: center; font-size: 75%;" 
| colspan="7"  |'''Figure 14:''' Model for the final blow simulation
|}


{| class="wikitable" style="text-align: center; margin: 1em auto;"
|+ <span id='table-5'></span>Table. 5 Model for the final blow simulation
|- style="border-top: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" |    Points 
| style="border-left: 2px solid;border-right: 2px solid;" | x 
| style="border-left: 2px solid;border-right: 2px solid;" | y 
| style="border-left: 2px solid;border-right: 2px solid;" | z 
|- style="border-top: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" | A 
| style="border-left: 2px solid;border-right: 2px solid;" | <math>0</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>0</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>-0.2815</math>
|-
| style="border-left: 2px solid;border-right: 2px solid;" | B 
| style="border-left: 2px solid;border-right: 2px solid;" | <math>-0.0346</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>0.0</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>-0.30848</math>
|-
| style="border-left: 2px solid;border-right: 2px solid;" | C 
| style="border-left: 2px solid;border-right: 2px solid;" | <math>-0.042</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>0</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>-0.3020</math>
|-
| style="border-left: 2px solid;border-right: 2px solid;" | D 
| style="border-left: 2px solid;border-right: 2px solid;" | <math>-0.0429</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>0</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>-0.20632</math>
|-
| style="border-left: 2px solid;border-right: 2px solid;" | E 
| style="border-left: 2px solid;border-right: 2px solid;" | <math>-0.0275</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>0</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>-0.1127</math>
|-
| style="border-left: 2px solid;border-right: 2px solid;" | F 
| style="border-left: 2px solid;border-right: 2px solid;" | <math>-0.0163</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>0</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>-0.0662</math>
|-
| style="border-left: 2px solid;border-right: 2px solid;" | G 
| style="border-left: 2px solid;border-right: 2px solid;" | <math>-0.0144</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>0</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>-0.0185</math>
|-
| style="border-left: 2px solid;border-right: 2px solid;" | H 
| style="border-left: 2px solid;border-right: 2px solid;" | <math>-0.0145</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>0</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>-0.0014</math>
|-
| style="border-left: 2px solid;border-right: 2px solid;" | I 
| style="border-left: 2px solid;border-right: 2px solid;" | <math>-0.0100</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>0</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>0</math>
|-
| style="border-left: 2px solid;border-right: 2px solid;" | J 
| style="border-left: 2px solid;border-right: 2px solid;" | <math>-0.0086</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>0</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>-0.0185</math>
|-
| style="border-left: 2px solid;border-right: 2px solid;" | K 
| style="border-left: 2px solid;border-right: 2px solid;" | <math>-0.01004</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>0</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>-0.0803</math>
|-
| style="border-left: 2px solid;border-right: 2px solid;" | L 
| style="border-left: 2px solid;border-right: 2px solid;" | <math>-0.01139</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>0</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>-0.1420</math>
|-
| style="border-left: 2px solid;border-right: 2px solid;" | M 
| style="border-left: 2px solid;border-right: 2px solid;" | <math>-0.01874</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>0</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>-0.2241</math>
|-
| style="border-left: 2px solid;border-right: 2px solid;" | N 
| style="border-left: 2px solid;border-right: 2px solid;" | <math>-0.01963</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>0</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>-0.2440</math>
|-
| style="border-left: 2px solid;border-right: 2px solid;" | O 
| style="border-left: 2px solid;border-right: 2px solid;" | <math>0</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>0</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>-0.2447</math>
|-
| style="border-left: 2px solid;border-right: 2px solid;" | P 
| style="border-left: 2px solid;border-right: 2px solid;" | <math>0</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>0</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>-0.2815</math>
|-
| style="border-left: 2px solid;border-right: 2px solid;" | Q 
| style="border-left: 2px solid;border-right: 2px solid;" | <math>-0.028</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>0</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>-0.2276</math>
|-
| style="border-left: 2px solid;border-right: 2px solid;" | R 
| style="border-left: 2px solid;border-right: 2px solid;" | <math>-0.034</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>0</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>-0.27</math>
|-
| style="border-left: 2px solid;border-right: 2px solid;" | S 
| style="border-left: 2px solid;border-right: 2px solid;" | <math>-0.029</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>0</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>-0.22</math>
|- style="border-bottom: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" | T 
| style="border-left: 2px solid;border-right: 2px solid;" | <math>-0.0246</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>0</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>-0.169</math>

|}
.

The 3D model can be obtained by applying rotation to the surface (representing the parison) and the curve (mold wall) around the vertical axis (see Fig. fig:final_blow_bench_model_a).  The location of the points necessary for reproducing the geometry is shown in Fig. fig:final_blow_bench_model_b and the corresponding coordinates are summarized in Table  [[#table-5|5]]. The mold and the glass domain share the point <math display="inline">G</math>. The bottle neck and the mold are fixed. The domain boundaries are defined as follows: Dirichlet boundary <math display="inline">\Gamma _d</math> comprises of the bottle neck and the mold wall.  The mold wall becomes part of the computationally domain <math display="inline">\Omega </math> only once the contact with the glass is established (see Section  [[#3.2 Contact modeling and boundary mesh|3.2]] for details). The mold wall, the glass domain as well as Dirichlet and Neumann boundaries are shown in Fig. [[#img-15|15]].  <div id='img-15'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_García-Espinosa_370675429-domain_parts.png|600px|The initial domain and the location of Dirichlet and Neumann boundaries in the ''final blow'' mold.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 15:''' The initial domain and the location of Dirichlet and Neumann boundaries in the ''final blow'' mold.
|}

We shall also distinguish the free outer boundary of the glass domain <math display="inline">\Gamma _i</math>. Note, that as soon as the glass node belonging to <math display="inline">\Gamma _i</math> touches the wall, it becomes a part of <math display="inline">\Gamma _d</math>.

The temperature at mold walls is fixed to 800 <math>\quad ^{\circ }</math>C. At the bottle neck the temperature is fixed to 724 <math>\quad ^{\circ }</math>C.  At the inner and outer surfaces of the glass the initial temperature is set to 950 <math>\quad ^{\circ }</math>C. Elsewhere in the glass domain, the initial temperature distribution is prescribed according to <math display="inline">T_{init}=950+(1140-950)(1-\frac{z-z_{bot}}{z_{top}-z_{bot}})^3</math>, where <math display="inline">z</math> is the vertical coordinate of the given node, <math display="inline">z_{bot}</math> and <math display="inline">z_{top}</math> are the minimum and maximum vertical coordinates of the parison. The initial temperature distribution is shown in Fig. [[#img-14a|14a]]. During first two seconds the parison is exposed to gravity exclusively. At t=2 s air pressure of 0.14 MPa is applied at the inner cavity. At t=2.6 the process is completed. The glass density is <math display="inline">\rho =</math> 2400 <math display="inline">kg/m^3</math>. The viscosity is computed from the temperature field according to: <math display="inline">\mu=265677693762693 e^{-0.0233569026\cdot T}</math>. The industrial partner (see Acknowledgements) has provided the authors directly with the viscosity values for a  number of temperatures, rather than the values of Fulcher constants. Thus, simple exponential fitting to find a function approximating this data was used.

====4.2.1 Results====

Fig. [[#img-16|16]] shows the evolution of the glass object.  One can see that around 0.3 seconds glass reaches the mold and accommodates there until the air pressure is applied. It takes less than 0.5 seconds to press the glass into the mold.

<div id='img-16'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_García-Espinosa_370675429-temp00.png|400px|]]
|[[Image:draft_García-Espinosa_370675429-temp03.png|400px|]]
|-
|[[Image:draft_García-Espinosa_370675429-temp1.png|400px|]]
|[[Image:draft_García-Espinosa_370675429-temp2.png|400px|]]
|-
|[[Image:draft_García-Espinosa_370675429-temp2-01.png|400px|]]
|[[Image:draft_García-Espinosa_370675429-temp2-02.png|400px|]]
|-
|[[Image:draft_García-Espinosa_370675429-temp2-03.png|400px|]]
|[[Image:draft_García-Espinosa_370675429-temp2-5.png|400px|Gravity stretching and final blow]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 16:''' Gravity stretching and final blow
|}

Fig. [[#img-17|17]] shows the temperature distribution along the vertical coordinate at the inner bottle surface. In the beginning of  the simulation temperature of 950<math>\quad ^{\circ }</math>C is set at inner (and outer)  surfaces (green line). The initial temperature distribution in the middle of the wall (thickness) is  shown in blue. The temperature (and, consequently, viscosity) at the surface of the glass evolves due to convection and diffusion leading to the final temperature distribution at the inner surface (shown in red).

Fig. [[#img-18|18]] displays the evolution of the glass volume. In case of using the standard PFEM algorithm (not using the fictitious contact elements approach) the total volume increases by as much as 35 <math display="inline">%</math>. One can distinguish two stages: the one corresponding to the gravity stretch (the contact between  the parison on the mold becomes established at around 0.3 seconds) and the one of the air pressure application. Since the contact area during the gravity stretch is mainly restricted to the bottom of the bottle, the overall volume growth due to contact is minor. On the other hand, air presses the entire outer surface of the bottle towards the mold wall. This process starts at t=2 seconds and at around t=2.3 seconds the entire outer surface of the glass is in contact with the wall. A major spurious volume increment is observed at this stage. From this point onwards no considerable volume change takes place.

<div id='img-17'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 60%;max-width: 100%;"
|-
|[[Image:draft_García-Espinosa_370675429-temperetaure_height.png|600px|Temperature distribution along the bottle height]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 17:''' Temperature distribution along the bottle height
|}

<div id='img-18'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 60%;max-width: 100%;"
|-
|[[Image:draft_García-Espinosa_370675429-volume_evolution.png|600px|Temporal evolution of glass volume]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 18:''' Temporal evolution of glass volume
|}

<div id='img-19a'></div>
<div id='img-19'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 60%;max-width: 100%;"
|-
|[[Image:draft_García-Espinosa_370675429-bottle_thickness.png|400px|]]
|[[Image:draft_García-Espinosa_370675429-thickness_along_wall_long1.png|400px|Wall thickness distribution along the bottle height]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 19:''' Wall thickness distribution along the bottle height
|}

On the other hand, if the fictitious contact elements algorithm is used, the volume conservation considerably improves. One can observe the line corresponding to  the apparent volume (actual volume of the glass + volume of the fictitious contact elements): the apparent volume increases when the contact elements are created, but since these elements provide no resistance it reduces to the actual one when these contract and finally vanish. Thus the apparent and the actual glass volumes coincide at the end of  the simulations (and, practically coincide at the end of each sub-stage: the gravity stretch and the air blow). In case of using the proposed contact algorithm,  volume gain reduces to 5 <math display="inline">%</math>.

Fig. [[#img-19|19]] shows the wall thickness distribution at t=2.6 s (vertical axis being the distance from the top measured along the wall coordinate L). This distribution can be used as the main reference data for the future comparisons.

====4.2.2 Summary and conclusions====

In this paper we have presented a 3D Lagrangian thermo-mechanical model for glass forming simulation. The model follows the PFEM philosophy combining the features of a classical FEM and a Lagrangian mesh-free method. A modified fractional step methodology has been introduced for obtaining a computationally efficient solution without compromising on the mass conservation.

The model was validated using a TV-panel pressing example, a test characterized by features similar to bottle forming but with a simpler geometry. A numerical example dealing with the ''gravity stretch'' and ''final blow'' stage of the bottle manufacturing process has been proposed.  We intended to create an example that would fulfill the need of a standardized reference lacking up-to-date in the field of bottle manufacturing simulation. All the data necessary for reproducing this test (geometry, material's properties, boundary conditions) has provided. Even though several strong approximations are introduced in this example (such as simplified thermal contact, initial temperature distribution), the example defines a first approximation to a benchmark in the field.

The final blow process was simulated in a reasonable computational time (around 1 hour using a quad-core I7 PC) on a computational mesh containing around 300 000 four-noded tetrahedra elements. The facility of the free-surface evolution tracking together with the automatic transfer of the non-homogeneous and temperature-dependent material properties proves that Lagrangian approaches are  more advantageous for the problem at hand than Eulerian ones. Thermal contact did not require any special technique, being an intrinsic feature of a matching-mesh method. The proposed mechanical contact methodology is promising, allowing the user to control the error in the thickness direction.  In the present work the geometrical contact tolerance  (i.e. the critical distance leading to erasal of the interface node of a fictitious element) was set to 10 <math display="inline">%</math> of the element size. The simulations carried out have proven that the method possesses excellent mass conservation features ( <math display="inline">\approx </math> 5 <math display="inline">%</math> of volume variation for the given contact tolerance).

Keeping in mind all the advantages of the formulation, it is important to note that it also has some limitations.  For optimal functionality of the method the mesh size distribution must be as uniform as possible. Thus only the mesh quality can be easily controlled. Moreover, the time step size was restricted due to the danger of element inversion. Using a novel Jacobian-independent Lagrangian explicit stream-line temporal integration <span id='citeF-28'></span>[[#cite-28|[28]]] is a promising option that must be investigated further. While keeping the overall architecture of the approach proposed here, this alternative Lagrangian formulation may lead to considerable advantages in computational efficiency removing the time step restrictions faced by the present method. It is important to note that bottle forming simulations deal with a problem involving very large changes in viscosity due to temperature variation. Thus, purely mechanical simulations are meaningless for the practical purposes. Even though the heat equation was included in the present model and viscosity was computed according to the changing temperature field, the present model did not account for radiative heat transfer. It is essential to include this feature  in the next development step of this numerical method.

==Acknowledgements==

This work was supported under the auspices of the FPDI-2013-18471 grant of the Spanish ''Ministerio de Economia y Competitividad''. It was also partially funded by the NUMEXAS project of the European Comission (FP7) and the COMETAD project of the National RTD Plan (ref. MAT2014-60435-C2-1-R) of the Spanish ''Ministerio de Economia y Competitividad''. The authors of this work express their gratitude to Dr. J. Jimenez from Vidrio&Engineering COMM. V. for the consultations regarding the bottle forming technology and, in particualr, the data provided for Example 2.

===BIBLIOGRAPHY===

<div id="cite-1"></div>
'''[[#citeF-1|[1]]]''' Oñate, E. and Idelsohn, S. and Del Pin, F. and Aubry, R.. ''The Particle Finite Element Method: an overview.'', Volume 1. International Journal of Computational Methods 267-307, 2004

<div id="cite-2"></div>
'''[[#citeF-2|[2]]]''' Ryzhakov, P. and Oñate, E. and Rossi, R. and Idelsohn, S.. ''Improving mass conservation in simulation of incompressible flows'', Volume 90/12. Int. Jour. for Num.Methds. in Eng. 1435-1451, 2012

<div id="cite-3"></div>
'''[[#citeF-3|[3]]]''' Ryzhakov, P.. ''A modified fractional step method for fluid-structure interaction problems''. Revista Intern. Met. Num. Ing. (RIMNI), 2016

<div id="cite-4"></div>
'''[[#citeF-4|[4]]]''' Cesar de Sa, J.M.A.. ''Numerical modelling of glass forming processes'', Volume 3. Eng. Comput. 266&#8211;275, 1986

<div id="cite-5"></div>
'''[[#citeF-5|[5]]]''' Matthew, H.. ''Numerical simulation of glass forming and conditioning'', Volume 85. Wiley Online Library. Journal of the American Ceramic Society 5 1047&#8211;1056, 2002

<div id="cite-6"></div>
'''[[#citeF-6|[6]]]''' Feulvarch, E. and Moulin, N. and Saillard, P. and Lornage, T. and Bergheau, J-M. ''3D simulation of glass forming process'', Volume 164. Elsevier. Journal of materials processing technology 1197&#8211;1203, 2005

<div id="cite-7"></div>
'''[[#citeF-7|[7]]]''' Giannopapa, C.G. and Groot, J.. ''Modeling the blow-blow forming process in glass container manufacturing: a comparison between computations and experiments'', Volume 133. American Society of Mechanical Engineers. Journal of Fluids Engineering 2 021103, 2011

<div id="cite-8"></div>
'''[[#citeF-8|[8]]]''' Fulcher, G.S.. ''Analysis of recent measurements of the viscosity of glasses'', Volume 8. Wiley Online Library. Journal of the American Ceramic Society 6 339&#8211;355, 1925

<div id="cite-9"></div>
'''[[#citeF-9|[9]]]''' Delaunay, B.. ''Sur la sphere vide'', Volume 7. Izvestia Akademii Nauk SSSR, Otdelenie Matematicheskikh i Estestvennykh  Nauk 793-800, 1934

<div id="cite-10"></div>
'''[[#citeF-10|[10]]]''' Donea, J. and Huerta, A.. ''Finite element method for flow problems'', J. Wiley Edition, 2003

<div id="cite-11"></div>
'''[[#citeF-11|[11]]]''' Codina, R.. ''A stabilized finite element method for generalized stationary incompressible  flows'', Volume 190. Computer Methods in Applied Mechanics and Engineering 20-21 2681 - 2706, 2001

<div id="cite-12"></div>
'''[[#citeF-12|[12]]]''' Oñate E.. ''A stabilized finite element method for incompressible viscous flows  using a finite increment calculus formulation.'', Volume 182. Computer Methods in Applied Mechanics and Engineering 355-370, 2000

<div id="cite-13"></div>
'''[[#citeF-13|[13]]]''' Ryzhakov, P. and Rossi, R. and Idelsohn, S. and Oñate, E. . ''A monolithic Lagrangian approach for fluid-structure interaction  problems'', Volume 46/6. Journal of Computational Mechanics 883-399, 2010

<div id="cite-14"></div>
'''[[#citeF-14|[14]]]''' Ryzhakov, P. and Cotela, J. and Rossi, R. and Oñate, E.. ''A two-step monolithic method for the efficient simulation of incompressible flows'', Volume 74. Wiley Online Library. International Journal for Numerical Methods in Fluids 12 919&#8211;934, 2014

<div id="cite-15"></div>
'''[[#citeF-15|[15]]]''' Idelsohn, S. and Oñate, E. . ''The challenge of mass conservation in the solution of free-surface  flows with the fractional-step method: problems and solutions'', Volume 26. International Journal for Numerical Methods in Biomedical Engineering 1313-1330, 2010

<div id="cite-16"></div>
'''[[#citeF-16|[16]]]''' Chorin, A.J.. ''A numerical method for solving incompressible viscous problems'', Volume 2. Journal of Computational Physics 12-26, 1967

<div id="cite-17"></div>
'''[[#citeF-17|[17]]]''' Temam, R.. ''Sur l’approximation de la solution des equations de Navier-Stokes  par la methode des pase fractionaires'', Volume 32. Archives for Rational Mechanics and Analysis 135-153, 1969

<div id="cite-18"></div>
'''[[#citeF-18|[18]]]''' Codina, R.. ''Pressure stability in fractional step finite element method for incompressible  flows'', Volume 170. Journal of Computational Physics 112-140, 2001

<div id="cite-19"></div>
'''[[#citeF-19|[19]]]''' Peri, D. and Owen, D.. ''Computational model for 3-D contact problems with friction based on the penalty method'', Volume 35. Wiley Online Library. International journal for numerical methods in engineering 6 1289&#8211;1309, 1992

<div id="cite-20"></div>
'''[[#citeF-20|[20]]]''' Papadopoulos, P. and Taylor, R.. ''A simple algorithm for three-dimensional finite element analysis of contact problems'', Volume 46. Elsevier. Computers & Structures 6 1107&#8211;1118, 1993

<div id="cite-21"></div>
'''[[#citeF-21|[21]]]''' Simo, J. and Laursen, T.. ''An augmented Lagrangian treatment of contact problems involving friction'', Volume 42. Elsevier. Computers & Structures 1 97&#8211;116, 1992

<div id="cite-22"></div>
'''[[#citeF-22|[22]]]''' Akkiraju, N. and Edelsbrunner, H. and Facello, M. and Fu P. and Mucke,  E. P. and Varela C. ''Alpha shapes: definition and software''. Proceedings of International Computational Geometry Software Workshop, 1995

<div id="cite-23"></div>
'''[[#citeF-23|[23]]]''' Grégoire, S. and Cesar de Sá, J. and Moreau, P. and Lochegnies, D.. ''Modelling of heat transfer at glass/mould interface in press and blow forming processes'', Volume 85. Elsevier. Computers & structures 15 1194&#8211;1205, 2007

<div id="cite-24"></div>
'''[[#citeF-24|[24]]]''' Codina R. and Vazquez M. and Zienkiewizc O.C.. ''A general algorithm for compressible and incompressible flows'', Volume 27. International Journal for Numerical Methods in Fluids 13-32, 1998

<div id="cite-25"></div>
'''[[#citeF-25|[25]]]''' Ryzhakov, Pavel and Rossi, Riccardo and Oñate, Eugenio. ''An algorithm for the simulation of thermally coupled low speed flow problems'', Volume 70. Wiley Online Library. International Journal for Numerical Methods in Fluids 1 1&#8211;19, 2012

<div id="cite-26"></div>
'''[[#citeF-26|[26]]]''' Berndhauser, C.. ''TC25 - Modelling of glass forming processes Short review of activities 2000 - 2009''. International Comission on Glass, 2013

<div id="cite-27"></div>
'''[[#citeF-27|[27]]]''' Augustin M. and Caiazzo A. and Fiebach, A. and Fuhrmann, J. and Volker J. and Linke A. and Umla R.. ''An assessment of discretizations for convection-dominated convection&#8211;diffusion equations'', Volume 200. Elsevier. Computer Methods in Applied Mechanics and Engineering 47 3395&#8211;3409, 2011

<div id="cite-28"></div>
'''[[#citeF-28|[28]]]''' Idelsohn, S.R. and Nigro, N. and Gimenez, J. and Rossi, R. and Marti, J. ''A fast and accurate method to solve the incompressible Navier-Stokes equations'', Volume 30. Emerald Group Publishing Limited. Engineering Computations 2 197&#8211;222, 2013