- Email: [email protected]
mould interface in press and blow forming processes
Computers and Structures 85 (2007) 1194–1205 www.elsevier.com/locate/compstruc
Modelling of heat transfer at glass/mould interface in press and blow forming processes Se´bastien Gre´goire a, Jose´ M.A. Ce´sar de Sa´ a
b,*
, Philippe Moreau a, Dominique Lochegnies
a
Laboratoire d’Automatique, de Me´canique et d’Informatique Industrielles et Humaines, UMR CNRS 8530, Universite´ de Valenciennes, Le Mont-Houy, Jonas 2, 59313 Valenciennes Cedex 9, France b Departamento de Engenharia Mecaˆnica e Gesta˜o Industrial, Faculdade de Engenharia, Universidade do Porto, Rua Dr. Roberto Frias, s/n, 4200-465 Porto, Portugal Received 3 March 2006; accepted 21 November 2006 Available online 23 January 2007
Abstract Numerical models may play an important role in the optimization of the quality of hollow-ware glass articles in glass industry. Due to the complexity of the phenomena involved a coupling between thermal and mechanical aspects is crucial. One of the key aspects required is a deep knowledge of the heat transfer at the glass/mould interface, a factor that may influence dramatically the final product. In this work an interface element developed to deal with heat transfer between glass and moulds in glass forming processes was introduced in a commercial code as a user element. The initial purpose was to make it possible to consider the dependence of the interface heat transfer resistance on process parameters like glass temperature and viscosity, mould temperature and contact pressure. In addition its implementation also proved to be an efficient tool to eliminate some perturbations introduced by the penalised form adopted in the commercial code to deal with non-matching meshes. A press and blow forming simulation of a Champagne flute was performed to assess the element performance on a real case study. 2006 Elsevier Ltd. All rights reserved. Keywords: Heat conduction; Finite elements; Heat transfer coefficient; Interface element; Glass forming; Press/blow process
1. Introduction The forming processes of hollow-ware glass articles is a complex coupled thermal/mechanical problem with interaction between the heat transfer analysis and the viscous flow of molten glass. The glass viscosity, which governs the glass flow during the forming process, is strongly dependent on the temperature [1,2] and concomitantly the contact with the moulds and the geometry changes, due to gravity and pressure, alter significantly the heat transfer process. The final product thickness distribution is one of the aspects which is more affected by the changing
*
Corresponding author. Tel.: +351 225081738; fax: +351 225081445. E-mail address: [email protected] (J.M.A. Ce´sar de Sa´).
0045-7949/$ - see front matter 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruc.2006.11.023
process conditions [2,4–7]. When the glass hits the walls of the moulds the local rapid cooling results in a considerable increase of the viscosity in theses zones and, consequently, in a more significant stretching and thinning of the zones where the temperature is higher. A perfect balance of the process conditions, meaning a perfect control of the highly dependence of viscosity on temperature, must therefore be achieved to obtain a good final product [1–8]. Manufacturers are facing an evolving and increasing competition in the development of new products and therefore software tools may play an important role in their goal to obtain a right – first time product design. Software tools are available which are capable of dealing with the most important aspects of these forming processes [9]. Nevertheless some important particular aspects of glass forming processes need a dedicated software or, at least, some user intervention in order to adapt commercial codes to this particular
S. Gre´goire et al. / Computers and Structures 85 (2007) 1194–1205
field, which may not be an easy task. One of these aspects is the modelling of the heat resistance between glass and mould at their interface which will be addressed here. The heat resistance at the interface is an average parameter thermally characterising the interface between the glass and the mould, which may be composed of air, lubricant or combustion products from the lubrication [10–14]. That parameter is a function of the interface pressure, of the glass temperature and therefore its viscosity, and of the asperities in the mould and its temperature. For the same pressure, the higher the temperature of the glass the smaller is the viscosity and therefore the physical contact is more effective as the glass fills better the micro cavities in the mould surface. When the glass cools against the mould walls then rapidly the viscosity rises near to a solid state and the glass may not even attain or may lose the contact with the mould as the pressures involved are usually very low [15,16]. When a simulation code, based on the finite element method, is used to deal with the processes described, typically glass and mould regions are modelled with nonmatching meshes [16,17]. Perfect mechanical and thermal contact is usually treated by a Lagrange multiplier technique or approximately by penalty formulations [17,18]. In particular this last approach may be used to model non-perfect heat contact by varying the penalty parameter, which, in fact, acts as an equivalent heat transfer coefficient between the two materials at the interface. Unfortunately for the user of these codes it is almost impossible to manipulate those parameters in order to make them directly or indirectly a function of the pressure, viscosity, temperature of the glass and mould or other local defects. These facts triggered the study presented here in which an interface finite element is implemented in the commercial code Abaqus in order to deal with heat transfer between glass and mould for non-matching meshes and where the heat transfer coefficient can be made, locally, a function of pressure, temperature or viscosity. The layout of this paper is the following: in Section 2 the governing equations in heat transfer problems are briefly introduced and the inclusion of the interface element justified; the definition of the interface element and its implementation in a commercial code are presented in Section 3 and a test to appraise its performance is carried out in Section 4. In the last section, the influence of the heat transfer coefficient on a hollow-ware glass article is studied through a simulation of a Champagne flute.
‘‘attached’’ to the material. In this case, for each time step, both in glass, Xg, and mould, Xm, Fig. 1, the governing differential equation is kr2 T þ Q ¼ 0
ð1Þ
where T is the temperature, Q is a heat source term and k the thermal conductivity of the medium. Regarding only the glass domain the natural boundary conditions may be of heat flux at the interface with the mould, in C1, and in the zones in contact with the surrounding air, in C2, i.e., oT ¼ h1 ðT T m Þ in C1 ð2Þ kg on1 oT kg ¼ h2 ðT T a Þ in C2 ð3Þ on2 where h1 and h2 are the corresponding heat transfer coefficients and Tm and Ta are, respectively, the mould temperature at the interface and the ambient temperature. Essential boundary conditions, i.e., prescribed temperatures, may exist in other part of the boundary, in C3, ð4Þ
T ¼ T in C3
It is possible to obtain an equivalent weak form for the problem by writing Z Z oT ðk g r2 T þ Qg ÞdT dXg þ k g h1 ðT T m Þ dT dC on1 Xg C1 Z oT þ k g h2 ðT T a Þ dT dC ¼ 0 ð5Þ on 2 C2 Performing an integrating by parts in the first term we obtain Z Z oT dT dC ðk g rT :rðdT Þ þ Qg dT Þ dXg þ kg on Xg C Z oT þ k g h1 ðT T m Þ dT dC on 1 C1 Z oT k g h2 ðT T a Þ dT dC ¼ 0 ð6Þ þ on2 C2 Noting that some flux terms cancel in boundaries C1 and C2, and that in boundary C3 the variation on T is zero then finally Eq. (6) may be written as
Ta
Ωm
n2 Γ2
2. Heat conduction on the glass: weak and strong forms. Interface element for natural boundary conditions In order to introduce the definition of the interface element let us state briefly the heat transfer problem in the process addressed. For simplicity let us assume, in the glass, that radiation effects are neglected and that, in our model, for each time step only diffusion takes place, i.e., the heat is convected by the moving Lagrangian mesh
1195
Ωg Γ3 T
Γ1
n1
Tm
T Γ2
n2 Ta Fig. 1. Glass and mould domains. Boundary conditions in the glass.
S. Gre´goire et al. / Computers and Structures 85 (2007) 1194–1205
1196
Z
ðk g rT rðdT Þ þ Qg dT Þ dXg Xg Z ðh2 ðT T a ÞÞdT dC ¼ 0
Z
ðh1 ðT T m ÞÞdT dC C1
ð7Þ
C2
This last equation may be viewed as a variational principle dPg ¼ 0
ð8Þ
associated with the functional Z 1 1 2 2 k g krT k Qg T dXg þ h1 ðT T m Þ dC Xg 2 C1 2 Z 1 h2 ðT T a Þ2 dC þ ð9Þ 2 C2
Pg ¼
Z
transfer coefficient assumes the value of h2. When contact is established the boundary condition is represented by Eq. (2). The heat coefficient h1 may then be defined, at the sampling integrating points of the interface element, as a function of glass viscosity, temperature or pressure which is known at the nodes of the underlying element. Likewise the mould temperature Tm at each contacting node may be obtained by determining its position, in terms of natural coordinates, on the side of element on the mould and subsequently using the shape functions to interpolate its value from the mould node temperatures. The details of the implementation of this interface element are discussed in the following sections. 3. Thermal interface element 3.1. FEM formulation
If the thermal contact between glass and mould would be perfect, corresponding to an infinite value of the heat transfer coefficient at the interface, then the temperatures at the two sides of the interface would be the same at each point. The second term of the right hand side of the functional of Eq. (9) would then be viewed as a penalty term imposing equality of temperatures on both sides. The parameter h1 would then act as a penalty parameter assuming a large number corresponding to a big heat transfer coefficient. If the contact is not perfect and some interface thermal resistance is present then the parameter h1 would be smaller, corresponding to a local equivalent heat transfer coefficient, and resulting in different temperatures in the glass and the mould at the interface. Based on this evidence an interface element is introduced to specially deal with that term on the functional that corresponds, on the strong form, to Eq. (2). On the mesh representing the glass, Fig. 2, an element boundary that is or may enter in contact with the mould is assumed to be insulated. An interface element is then superposed on that boundary and takes care of the heat transfer to the air and to the moulds. If there is no contact the heat
The finite element formulation may be obtained from the weak form in Eq. (7). In particular in what respects the interface element the discretisation is obtained from the term Z ... ðh1 ðT T m ÞÞdT dC . . . ð10Þ C1
The contribution of each interface element to the discretised form of (10) is . . . dT e H e T e dT e f e . . . e
ð11Þ e
in which T are the element nodal temperatures, H is the element thermal ‘‘stiffness’’ matrix and f e the element contribution for the overall thermal ‘‘force’’ vector. Theses matrices are obtained from the element shape functions using the standard procedure for isoparametric axisymmetric elements (see Fig. 3 for details) by performing the integrations in the natural coordinates as Z 1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi e T H ¼ 2p h1 ðnÞrðnÞN N ðox=onÞ2 þ ðoy=onÞ2 dn ð12Þ 1 Z 1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 e T ðox=onÞ þ ðoy=onÞ dn ð13Þ f ¼ 2p h1 ðnÞrðnÞT m N 1
3.2. Implementation of the interface element on a FEM code
Axisymetric axis
Glass contact surface
Interface elements
Mould contact surface
Fig. 2. Positioning of the interface elements.
3.2.1. Non-matching meshes As referred above typically non-matching meshes will be used for the glass and the mould. When a node of the mesh of the glass enters in contact with an element of the mould, the temperature at the mould in the contact point must be determined to establish the thermal contact. Since the mechanical contact is not perfect, as it is dealt with penalty functions, the contact nodes may be slightly outside or inside the element of the mould mesh. A safe procedure must be then performed to assure that the thermal contact is made. In fact the same procedure is also adopted whenever a remeshing operation is per-
S. Gre´goire et al. / Computers and Structures 85 (2007) 1194–1205
Nj : (1+ξ)/2
Ni : (1-ξ)/2 i
ξ1(-1/ √3)
ξ
j
ξ2(1/ √3) Y
i
Xi, Yi, hi, Ti
r(ξ)
dL
∂y ∂ξ
∂x ∂ξ
2
+
∂y ∂ξ
2
∂x ∂ξ
Lij
j
dL =
X
Xj, Yj, hj, Tj dS=2πr(ξ)dL
Fig. 3. Interface element: (a) in natural coordinates, (b) in global coordinates.
formed, as some nodes of the new mesh may slightly fall outside the mesh of the previous domain. A gap function is defined as 2
/ðn; gÞ ¼ ðxem ðn; gÞ xg Þ þ ðy em ðn; gÞ y g Þ
2
ð14Þ
The coordinates xem ; y em to where the temperatures of the mould element are interpolated (extrapolated) should then minimize the function /ðn; gÞ. This implies that the gradient of /ðn; gÞ should be zero, i.e., r/ðn; gÞ ¼ O
ð15Þ
The solution is obtained in an iterative manner using the Newton–Raphson method. Departing from ðn0 ; g0 Þ ¼ ð0; 0Þ at iteration n the following system is solved to obtain ðDnn ; Dgn Þ ( Dðo/=onÞnþ1 ¼ ðo2 /=on2 ÞDnn þ ðo2 /=onogÞDgn ð16Þ Dðo/=ogÞnþ1 ¼ ðo2 /=ogonÞDnn þ ðo2 /=og2 ÞDgn The natural coordinates are then updated and convergence criterion is checked by defining a very small tolerance (toler ’ 104) ðnnþ1 ; gnþ1 Þ ¼ ðnn ; gn Þ þ ðDnn ; Dgn Þ kðnnþ1 ; gnþ1 Þ ðnn ; gn Þk=kðnn ; gn Þk < toler
ð17Þ ð18Þ
The convergence is typically reached in two or three iterations. The temperature at the contact point, is then determined by using the shape functions of the element Tm ¼
nX odes
N i ðnc ; gc ÞT i
ð19Þ
i¼1
where the subscript c stands for converged value. This definition of Tm should then be included in Eq. (13) coupling the degrees of freedom in the glass and the mould and therefore allowing the heat transfer between the two materials to take place.
1197
3.2.2. Implementation in a commercial code One would like to include the features described in a commercial code where many facilities are available like different constitutive equations, element type choice, contact handling, remeshing operations and the usual pre and post processing tools. In this work the Abaqus commercial code was selected as the host code for the implementation of the proposed element, via the user subroutine UEL. Nevertheless part of the procedure indicated in the previous section, which should be pursued when dealing with a home code where the access is free, will be a difficult (impossible?) task to achieve when a commercial code is used. This is the case with the definition of Tm in Eq. (13). Whilst the node temperatures in the interface elements, as they are permanently ‘‘attached’’ to the same elements of the glass mesh, can be made equal to the corresponding glass temperatures using multiple point constraints, the definition of Tm as a function of the temperatures of the mould elements is a hard and cumbersome task in this case. This is due not only to the fact that, specially in the press operation, the contacting nodes are moving from one element to another in the contact surface but also that it is not possible to obtain internally the number (label) of the element of the mould where contact has been established. This can only be done indirectly using the procedure referred in the last section. Therefore a simplified staggered algorithm was established. The value of Tm in Eq. (13) is evaluated, in each iteration, according to Eq. (19). The glass cooling against the mould is, by this process, guaranteed. However, because the temperatures in the nodes of the mould mesh are not constrained internally with glass temperature values, an equivalent heat transfer coefficient is assumed at each contacting side of an element of the mould. Its value is a weighted average of the values of the heat transfer coefficients of the surrounding contacting nodes of the glass. 4. Assessment of the performance of the interface element A set of simple tests was carried out in order to validate the performance of the interface element and its implementation in the Abaqus code. In those tests, the solutions obtained with standard thermal and thermomechanical models in Abaqus are compared with those obtained with a thermomechanical model in Abaqus that includes the interface element. The main objective is to assess the heat transfer process between two non-matching meshes, which is present in glass forming simulation due to the evolving contact between glass and mould meshes. This is an important feature because the glass viscosity is strongly dependent on the temperature. 4.1. Tests description The axisymetric model in Fig. 4 represents two cylinders of two different materials, a crystal glass and a steel (XC38), initially put into contact. Each cylinder has a
S. Gre´goire et al. / Computers and Structures 85 (2007) 1194–1205
1198 Z 30mm Glass (1150˚C)
10mm r
nodes at the steel interface (Fig. 5b). The mesh for the steel cylinder is identical for both simulations (7 elements in Zdirection and 16 in R-direction). The glass has 7 elements in the Z-direction and is meshed by 16 elements for the matching meshes simulation and 21 elements for the nonmatching meshes in the R-direction (Fig. 5).
10mm Steel (485˚C)
4.2. Simulation with matching meshes
Fig. 4. Axisymetric model representing a glass and a steel cylinder into contact.
radius of 30 mm and a height of 10 mm. The initial temperatures are 1150 C on the glass and 485 C on the steel. The thermal properties of the two materials at these temperatures are given in Table 1. By using the interface element it would be possible to include at the glass/steel interface a heat transfer coefficient, simulating the heat resistance, dependant on the glass and the mould contact temperatures, on the contact pressure or in the local viscosity. However when Abaqus code is used only contact pressure or gap thickness at the interface can influence the heat transfer coefficient. Therefore, in order to compare the results with and without the interface element, a constant heat transfer coefficient of 5000 W/m2/ C was assumed at the interface. There is only heat transfer between the cylinders as they are insulated at the other faces. Only the self weight of the cylinders is considered as mechanical force. The friction coefficient between the two materials is set equal to the unity. The heat transfer process is analysed during an initial short period of 10 s. Two different meshes for the glass cylinder are used: one that matches the mould mesh at the interface (Fig. 5a) and a finer mesh with non-matching Table 1 Thermal properties of glass and steel at, respectively, 1150 C and 485 C
Thermal conductivity (W/m/C) Specific heat (J/kg/C) Density (kg/m3)
Glass (crystal) at 1150 C
Steel (XC38) at 485 C
2.77 1282 2500
37 620 7800
Fig. 5. Different used meshes: (a) matching meshes, (b) Non-matching meshes.
In the thermomechanical model with the interface element, Fig. 6a and in the standard thermal simulation Fig. 6c the same temperatures are obtained, with average differences of less than 0.8 C for all the length of the analysis. The temperature isocurves are straight lines parallel to the interface. But in the standard thermomechanical model, i.e., Abaqus code with no interface element, the temperatures are not identical along the steel contact surface (Fig. 6b) as the temperature isocurves depart slightly from parallel lines to the interface, specially closer to the axisymmetry axis where the temperatures are higher than the previous simulations. The evolution of the temperatures at the extreme points of the interface for both the glass and the steel sides are represented in Figs. 7 and 8. The temperature evolutions between the contact element model and the thermomechanical model are different on the left sampling point located at the steel contact surface. After 0.4 s, a difference of 12 C is noticed (Fig. 8a). At the same time, the same temperature difference is read at the left glass point (Fig. 7a). On the right comparison points, the temperatures evolutions of these three models are similar (Figs. 7b and 8b). 4.3. Results for non-matching meshes simulation In this case, both in the standard Abaqus thermal and thermomechanical models, the solution is perturbed due to the fact that the two meshes do not match. In Fig. 9b and c, it is clear that in these cases the temperature isocurves depart from parallel lines to the interface, especially near the axis. When the interface element is included in the analysis, these perturbations are eliminated and the correct solution is recovered (Figs. 9a and 10). The comparison of the temperature profiles at the interface are shown in Fig. 10 where its clearly seen how the perturbations on the temperature are eliminated with the inclusion of the proposed element. The evolution of the temperatures at the extreme points of the interface for both the steel and the glass sides are represented in Figs. 11 and 12. The temperature evolutions are similar for the heat transfer model and the thermomechanical (maximal difference of 4 C after 0.6 s) but they differ from the solution with the interface element, specially on the steel side (Fig. 11a) where differences as high as 24 C on the steel side are obtained in this simple example. It is known in glass forming how important is the temperature distribution in the mould for the final product [1–6].
S. Gre´goire et al. / Computers and Structures 85 (2007) 1194–1205
1199
Fig. 6. Temperature isocurves inside the steel at time 10 s: (a) contact element model, (b) thermomechanical model, (c) heat transfer model.
1150
1150
1100
Contact element model Thermomechanical model Heat transfer model
Temperature (°C)
1050 1000 950 900 850 800
Glass (1150˚C)
1050 Temperature (°C)
1100
1000 950 900 850 800 750
750 700
Glass node
650
Steel (485˚C)
700
Glass node
650 600
600 0
2
4
6
8
0
10
2
4
6
8
10
Contact time (s)
Contact time (s)
Fig. 7. Temperature evolution at two points on the glass interface for the three models (matching meshes).
620
620
Contact element model Thermomechanical model Heat transfer model
580 560 540
Glass (1150˚C)
520 Steel node
500
Steel (485˚C)
600 Temperature (°C)
Temperature (°C)
600
Steel node
580 560 540 520 500 480
480 0
2
4 6 8 Contact time (s)
10
0
2
4
6
8
10
Contact time (s)
Fig. 8. Temperature evolution at two points on the steel interface for the three models (matching meshes).
Fig. 9. Temperature isocurves inside the steel after 10 s of simulation: (a) contact element model, (b) thermomechanical model, (c) heat transfer model.
In the glass side also differences of 14 C are detected after 0.6 s (Fig. 12a) what may correspond, in the interval of
forming temperatures, to changes in the glass viscosity that may range from 20% to 500%.
S. Gre´goire et al. / Computers and Structures 85 (2007) 1194–1205
1200 680
Contact element model Thermomechanical model Heat transfer model
Glass interface
Temperature (˚C)
660
640 Glass
620
D
Steel Mould interface
600 0
5
10
15
20
25
30
Distance D from the axisymetric axis (mm)
Fig. 10. Glass and steel nodal temperature along the glass/steel interface after a run of 10 s.
620
620
Contact element model Thermomechanical model Heat transfer model
580
600 Temperature (°C)
Temperature (°C)
600
560 540 Glass (1150°C)
520
580 560 540 520
Steel node
500
Steel (485°C)
Steel node
500 480
480 2
0
8 6 4 Contact time (s)
0
10
2
8 6 4 Contact time (s)
10
Fig. 11. Temperature evolution at two points on the steel interface for the three models (non-matching meshes).
1150
1150
Contact element model Thermomechanical model Heat transfer model
1050 Temperature (°C)
1000 950 900 850 800
1100 1050 Temperature (°C)
1100
Glass (1150°C)
1000 950 900 850 800 750
750 700
Glass node
650
Steel (485°C)
Glass node
700 650 600
600 0
2
4 6 8 Contact time (s)
10
0
2
4 6 8 Contact time (s)
10
Fig. 12. Temperature evolution at two points on the glass interface for the three models (non-matching meshes).
4.4. Comments The results of the simple test described in the previous sections showed that the Abaqus standard thermal an ther-
mal–mechanics models may introduce some perturbations in the temperature distribution when non-matching meshes are used, which is probably due to the penalised form used to impose thermal and mechanical contact. This is a crucial
S. Gre´goire et al. / Computers and Structures 85 (2007) 1194–1205
1201
issue for glass forming process simulation due to the fact that glass viscosity is highly dependent on the temperature (see Eq. (20)). The inclusion of the interface element, with the initial purpose of making it possible to introduce the dependence of heat resistance on process parameters, proved to be also an efficient tool to eliminate those perturbations making the numerical model more robust and reliable.
with the thermomechanical axisymetric model on the Abaqus commercial code, in which the interface thermal element is included. The total forming time for the Champagne flute is divided in three phases:
5. Application to a press and blow process
The initial temperature of the plunger, the blank mould and the blow mould is 485 C and the initial temperature of the glass is 1150 C. The thermal conductivity and the specific heat of the XC38 steel used for the moulds and the punch at 600 C are, respectively, 33.9 W/m/C and 708 J/kg/C. At 1000 C the thermal conductivity and the specific heat of the crystal are, respectively, 2.5 W/m/C and 1265 J/ kg/C. The temperature dependences of the thermal conductivity and specific heat of the crystal and of the XC38 steel are presented in Figs. 15 and 16. Glass viscosity is highly dependent on the temperature and is commonly represented, in the forming range temperatures, by the Vogel, Fulcher and Thalmann law [21]. For the particular case of the crystal used in the Champagne flute this law takes the following form:
The forming of an upper part of a Champagne flute, requiring three sequential forming phases – press, gravity stretch and blow – is studied in this section. This forming process, and particularly some of its forming parameters, has already been studied by Noiret [19]. The study proposed here concerns the analysis of the influence of the heat transfer coefficient at the interface during the full process. Final temperature and thickness distributions will be the important parameters to look at carefully. 5.1. Experimental data In order to obtain representative heat transfer coefficients at the interface, an experimental device was designed [14,20] which puts into contact a crystal glass at an initial temperature of 1150 C and a XC38 steel plunger at an initial temperature of 485 C. A set of heat transfer coefficients were obtained for different increasing contact pressures and temperatures. Two of these cases are shown in Fig. 13 where it is clear how the pressure may influence the heat transfer value.
• a pressing time: 0.3 s (Fig. 14a and b); • a gravity stretching time: 4 s (Fig. 14c and d); • a blowing time: 7.5 s (Fig. 14d and e).
logðgÞ ¼ 1:9458 þ
5407 Pa s ðT 174:5Þ
ð20Þ
in which T is the temperature and g is the viscosity. In Fig. 17 this law is represented in the temperature range of [500 C, 1200 C].
5.2. Modelling of the forming processes In order assess the pressure influence two simulations of the full forming process, using the two heat transfer coefficient evolutions represented in Fig. 13, were performed
Heat transfer coefficient (W/m²/°C)
8000
6000
h1I ( t ) 4000
h1II ( t )
2000 0
1
2 3 4 Contact time (s)
5
6
Fig. 13. Heat transfer coefficient evolution obtained by experimentation. hI1 ðtÞ: increasing pressure from 0 to 3.0 MPa, hII 1 ðtÞ: increasing pressure from 0 to 1.5 MPa.
Fig. 14. Champagne flute forming phases (thermal isocurves). Scale 1:2. (a) Initial position; (b) end of pressing; (c) before stretching; (d) end of stretching; (e) end of blowing.
S. Gre´goire et al. / Computers and Structures 85 (2007) 1194–1205
1202
0
0.05
0.1
0.15
0.2
0.25
0.3
0
45 40
-0.01
35 moulds and plunger : XC38 steel
30 25
glass : cristal
20 15 10
Displacement (m)
Thermal conductivity (W.m-1.°C-1)
50
-0.02 -0.03 -0.04
5
-0.05
0 400
600 800 1000 Temperature (˚C)
1200
-0.06
Fig. 15. Thermal conductivity of glass, moulds and plunger.
Time (s)
Fig. 18. Plunger vertical displacement versus time.
• a heat transfer coefficient (as given in Fig. 13) in case of contact between glass and mould or glass and plunger.
1400
-1
Specific heat (J.kg-1.°C )
1200 1000 800
moulds and plunger : XC38 steel
600
glass : cristal
400 200 0 400
600 800 1000 Temperature (°C)
1200
Fig. 16. Specific heat of glass, moulds and plunger.
16
For this stage four remeshing operations were needed to obtain the final parison. The subsequent stage is a stretching operation of the glass part due to gravity only in which a certain homogenisation of the temperatures in the glass is obtained. Prior to begin this stage a new remeshing operation is performed due to the mesh distortion in the final period of the pressing stage. The thermal boundary conditions for the stretching stage are: • a convection with the ambient air (h2 = 20 W/m/C and Ta = 50 C) where there is no contact between glass and mould contact; • a heat transfer coefficient (as given in Fig. 13) in case of contact between glass and mould.
14
Log(η) (Pa.s)
12 10 8 6 4 2 0 500
600
700 800 900 1000 1100 1200 Temperature (°C)
Fig. 17. Temperature dependence of the crystal viscosity.
In the pressing stage, the parison (a hollow-shaped blank) is obtained after a displacement of 51 mm of the plunger (Fig. 18). The thermal boundary conditions are: • a convection with the ambient air (h2 = 20 W/m/C and Ta = 50 C) in the case of non-contact between the glass and mould or glass and plunger;
The mesh of the blowing phase is identical to the mesh of the stretching phase. In the previous study of Noiret [19] the blowing pressure was set to 1 kPa for a NBS 710 standard glass. The crystal glass adopted here has a much higher viscosity at the forming temperatures and therefore to perform the final blowing another pressure law, represented in Fig. 19, was used. The flute forming is completed after a blowing time between 2 s and 4 s (depending on the heat transfer coefficient, i.e., on the heat resistance at the interface). After full contact some extra time is needed to cool the glass in order to remove it out of the mould without damaging the final shape. The boundary conditions for the stretching phase are: • a convection with the pressurised air (h2 = 40 W/m/C and Ta = 200 C); • a convection with the ambient air (h2 = 20 W/m/C and Ta = 50 C) if contact has not been established; • a heat transfer at the glass and mould interface after contact has been established coefficient with the appropriate heat transfer coefficient h1(t).
S. Gre´goire et al. / Computers and Structures 85 (2007) 1194–1205
But, according to the same figure, all the temperatures in the neck have dropped significantly, which is explained by the thinning effect in this zone caused by the pressing phase. During the stretching phase, in both simulations, a homogenization of the temperatures in glass is noticed (Fig. 21b) but, as expected, at the end of this phase the final temperatures are lower for the case of the simulation with hI1 ðtÞ (Fig. 20). For the same time operations when the blowing phase begins the stiffness of the glass in the two cases will be different due to high dependence of the viscosity on the temperature. Therefore the contact between glass and mould begins 1.2 s later for the case of the simulation with hI1 ðtÞ due to the glass viscosity which is higher owing to lower temperatures. In both simulations, the temperatures level reached, in the end, at the mould’s sampling point have the same values but with an offset time of 1.2 s.
25
Air pressure (kPa)
20 15 10 5 0 0
1
2
3 4 5 Blowing time (s)
6
1203
7
Fig. 19. Evolution of the blowing pressure versus time.
5.3. Results and discussion In this section, by utilising the two different heat transfer coefficients described above, hI1 ðtÞ and hII 1 ðtÞ, the influence of the heat transfer coefficient on the temperature evolution inside the glass and the moulds and on the final thickness of the glass part is analysed. 5.3.1. Temperature field and stretching During the pressing, as expected, the glass surface cools more in the simulation with hI1 ðtÞ and, correspondingly, in the mould the temperature increases more. This is illustrated in Fig. 20 where the temperature evolutions, for two particular points in the glass and mould that will contact during the blow phase, are represented. In both simulations, no significant temperature variation in the heart of the glass in the lower zone of parison is noticed, particularly in the simulation with hI1 ðtÞ, as shown in Fig. 21a.
5.3.2. Thickness distribution The glass temperatures after the pressing phase are lower in the simulation with hI1 ðtÞ. Therefore the glass deformation during the gravity stretching is less important. In this case, the sampling node in glass bottom has a displacement of 8 mm after 4 s of stretching compared to a displacement of 13.5 mm at the same point in the simulation with hII 1 ðtÞ (Fig. 22). In Figs. 23 and 24 the thickness of the Champagne flute is represented and measured from the internal side of the flute along its height. In the simulation with hI1 ðtÞ, the final thickness is relatively homogeneous (an average of 2.4 mm). This can be explained by the fact that the relatively low temperatures in the neck restrict the deformations in that zone during the stretching phase. Indeed, at the end of the pressing,
1150 Stretching
Blowing
1050
Temperature (˚C)
950 Glass 850 750 650 550
Mould
450 0
2
4
6
8
10
12
punch + glass + blank mould
glass + blow mould
Process time (s)
Fig. 20. Temperature evolution on a glass and on a mould node on the contact surfaces during the full process. Thick curve: simulation with hI1 ðtÞ. Thin curve: simulation with hII 1 ðtÞ.
S. Gre´goire et al. / Computers and Structures 85 (2007) 1194–1205
1204
Fig. 21. Temperature field in the glass during the forming process in the simulation with hI1 ðtÞ. (a) End of pressing, (b) end of stretching, (c) end of blowing.
0
Displacement (m)
-0.002
h (t ) I 1
-0.004 -0.006
h1II ( t )
-0.008 -0.01 -0.012 -0.014 0
1
2 Creeptime (s)
3
4
Fig. 22. Displacement of a node in the glass part during stretching phase in the two simulations.
120mm 100
80
h1I ( t )
60
40
h1II ( t ) 20
0mm 0 1
2
3 4 Thickness (mm)
5
Fig. 24. Thickness of the Champagne flute at the end of the press–blow process.
120 Distance from the bottom of the glass (mm)
Distance from the bottom of the glass (mm)
120
100
100mm
80
60
40
0mm
20
0 2
4 6 Thickness (mm)
8
Fig. 23. Thickness of the Champagne flute at the end of the pressing.
that zone has a thickness that varies between 3 mm and 4 mm (Fig. 23) and at the end of the full process, the same zone has an average thickness of 2.5 mm. Most of the deformation is localized in the lower zone of parison. In that zone the thicknesses are varying between 5 mm and 7.5 mm after the pressing and are dropping to an average of 2.3 mm. For the simulation with the hII 1 ðtÞ, the glass is thinner in the neck of the flute (value of 2.1 mm) than in the lower zone of parison (4.5 mm). Here, contrary to the previous simulation, the glass is softer (low viscosity) in the neck and allows a bigger stretching. Then, the bottom zone keeps more or less the same thicknesses which decrease from values comprised between 5 mm and 7.5 mm to values comprised between 3 and 4.5 mm.
S. Gre´goire et al. / Computers and Structures 85 (2007) 1194–1205
6. Conclusions An interface finite element was implemented, at the University of Porto, in the commercial code Abaqus in order to deal with heat transfer between glass and mould for nonmatching meshes. The associated heat transfer coefficient can be, locally, set as a function of pressure, temperature or viscosity. The interface finite element was validated on a simple benchmark. It proved to perform well and to regularise some perturbations introduced by the penalised form adopted in the commercial code to deal with non-matching meshes. This is an important feature as the glass viscosity and therefore the final shape are strongly dependent on the temperature distribution. An application on the forming of a Champagne flute by a press and blow operation illustrates the influence of the heat transfer coefficient throughout the whole process. Two simulations are performed with two heat transfer coefficients obtained experimentally at the University of Valenciennes showing a great influence in the final product in terms of glass temperature values and distribution. These differences have a great impact on the glass viscosity and therefore on the forming operation times and conditions. The time for each phase (pressing, stretching and blowing), the value of the pressure at any instant and when it should be activated play a decisive role in the final shape and thicknesses distribution. To succeed the optimization of the quality of the manufactured glass products, the numerical thermomechanical models require a perfect knowledge of the heat exchange at the glass/mould interface. The procedure described here may play an important role in the robustness of the modelling of glass forming processes. Acknowledgements The present work is supported by the Department of Mechanical Engineering and Industrial Management of the University of Porto (Portugal), the Laboratoire d’Automatique de Me´canique et d’Informatique Industrielles et Humaines de l’Universite´ de Valenciennes et du Hainaut-Cambre´sis, the CNRS and the Ministe`re de l’Education Nationale et de la Recherche. The authors gratefully acknowledge the support of these institutions. References [1] Williams JH, Owen DRJ, Cesar de Sa JMA. The numerical modelling of glass forming process. In: Proceedings of the XIVth international congress on glass, New Delhi, India; 1986.
1205
[2] Ce´sar de Sa´ JMA. Numerical modelling of glass forming processes. Eng Comput 1986;3(6):266–75. [3] Storck K. Thermal system analysis. Heat transfer in glass forming and fluid temperature-control systems. Institute of Technology, Linko¨pings Universitet. Linko¨pings Studies in Science and Technology. Dissertation no. 542. Linko¨ping; 1998. [4] Lochegnies D, Mare´chal C. How glass properties can affect the glass flow? In: Proceedings of the 2nd international colloquium modelling of glass forming and tempering. Valenciennes, France: Presses Universitaires de Valenciennes; 2002. [5] Ce´sar de Sa´ JMA, Natal-Jorge RM, Cardoso RPR. Simulation of forming of hollow glass containers. In: Proceedings of the Euromech 388 colloqium on modelling of glass forming processes. Presses Universitaires de Valenciennes. Lez Valenciennes no. 27. Valenciennes, France; 1999. [6] Ce´sar de Sa´ JMA, Natal-Jorge RM, Silva CMC, Cardoso RPR. A computacional model for glass container forming processes. In: Wunderlich W, editor. Proceedings of ECCM’99 – European conference on computational mechanics, Munich, Germany; 1999. [7] Lochegnies D, Moreau P, Guilbaut R. A reverse engineering approach to the design of the blank mould for the glass and blow process. Glass Technol 2005;46(2):116–20. [8] Mare´chal C, Moreau P, Lochegnies D. Numerical optimization of a new robotized glass blowing process. Eng Comput 2004;19:233–40. [9] Hyre M, Rubin Y. Comprehensive modelling of the glass container forming process. In: Ken-ichiro Mori, editor. Simulation of materials processing: theory, methods and applications. In: Proceedings of the 7th international conference on numerical methods in industrial forming processes-Numiform 2001, Toyohashi, Japan; 2001. [10] Fellows CJ, Shaw F. A laboratory investigation of glass to mould heat transfer during pressing. Glass Technol 1978;19(1):4–9. [11] Pchelyakov SK, Guloyan YA. Heat transfer at the glass mould interface. Glass Ceram 1985;42:400–3. [12] Storck K, Loyd D, Augustsson A. Heat transfer modelling of the parison forming in glass manufacturing. Glass Technol 1998;39(6): 210–6. [13] Keijts G, van der Werff K. Heat transfer in glass container production during the final blow. Glass Technol 2001;42(3):104–8. [14] Moreau P, Gre´goire S, Lochegnies D. Me´thodologie expe´rimentale pour la de´termination du coefficient de transfert thermique a` l’interface verre/outil. Verre 2003;9(2):24–7. [15] McGraw DA. Transfer of heat in glass during forming. J Am Ceram Soc 1961;44(7):353–63. [16] Mori M, Shimazaki T, Naraki T, Wakatsuki H. Numerical simulation system for mould designer. In: Proceedings of the 2nd international colloquium modelling of glass forming and tempering. Valenciennes, France: Presses Universitaires de Valenciennes; 2002. [17] Ce´sar de Sa´ JMA, Natal-Jorge RM, Cardoso RPR. Improved finite element analysis on blow forming. In: Proceedings of the 2nd international colloquium modelling of glass forming and tempering. Valenciennes, France: Presses Universitaires de Valenciennes; 2002. [18] Wriggers P. Computational contact mechanics. England: John Wiley & Sons; 2002. [19] Noiret C, Thibaud C, Lochegnies D, Oudin J. Le formage du verre creux: l’approche nume´rique sur deux exemples. Verre 1996;2:19–22. [20] Gre´goire S, Moreau P, Cesar de Sa JMA. Integration of the heat transfer coefficient h(t) in glass forming modelling. In: Proceedings of GlassDays 2005, Kaiserslautern, Germany; 2005. [21] Fulcher GS. Analysis of recent measurements of the viscosity of glasses. J Am Ceram Soc 1925;8(6):339–55.
Modelling of heat transfer at glass/mould interface in press and blow forming processes Se´bastien Gre´goire a, Jose´ M.A. Ce´sar de Sa´ a
b,*
, Philippe Moreau a, Dominique Lochegnies
a
Laboratoire d’Automatique, de Me´canique et d’Informatique Industrielles et Humaines, UMR CNRS 8530, Universite´ de Valenciennes, Le Mont-Houy, Jonas 2, 59313 Valenciennes Cedex 9, France b Departamento de Engenharia Mecaˆnica e Gesta˜o Industrial, Faculdade de Engenharia, Universidade do Porto, Rua Dr. Roberto Frias, s/n, 4200-465 Porto, Portugal Received 3 March 2006; accepted 21 November 2006 Available online 23 January 2007
Abstract Numerical models may play an important role in the optimization of the quality of hollow-ware glass articles in glass industry. Due to the complexity of the phenomena involved a coupling between thermal and mechanical aspects is crucial. One of the key aspects required is a deep knowledge of the heat transfer at the glass/mould interface, a factor that may influence dramatically the final product. In this work an interface element developed to deal with heat transfer between glass and moulds in glass forming processes was introduced in a commercial code as a user element. The initial purpose was to make it possible to consider the dependence of the interface heat transfer resistance on process parameters like glass temperature and viscosity, mould temperature and contact pressure. In addition its implementation also proved to be an efficient tool to eliminate some perturbations introduced by the penalised form adopted in the commercial code to deal with non-matching meshes. A press and blow forming simulation of a Champagne flute was performed to assess the element performance on a real case study. 2006 Elsevier Ltd. All rights reserved. Keywords: Heat conduction; Finite elements; Heat transfer coefficient; Interface element; Glass forming; Press/blow process
1. Introduction The forming processes of hollow-ware glass articles is a complex coupled thermal/mechanical problem with interaction between the heat transfer analysis and the viscous flow of molten glass. The glass viscosity, which governs the glass flow during the forming process, is strongly dependent on the temperature [1,2] and concomitantly the contact with the moulds and the geometry changes, due to gravity and pressure, alter significantly the heat transfer process. The final product thickness distribution is one of the aspects which is more affected by the changing
*
Corresponding author. Tel.: +351 225081738; fax: +351 225081445. E-mail address: [email protected] (J.M.A. Ce´sar de Sa´).
0045-7949/$ - see front matter 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruc.2006.11.023
process conditions [2,4–7]. When the glass hits the walls of the moulds the local rapid cooling results in a considerable increase of the viscosity in theses zones and, consequently, in a more significant stretching and thinning of the zones where the temperature is higher. A perfect balance of the process conditions, meaning a perfect control of the highly dependence of viscosity on temperature, must therefore be achieved to obtain a good final product [1–8]. Manufacturers are facing an evolving and increasing competition in the development of new products and therefore software tools may play an important role in their goal to obtain a right – first time product design. Software tools are available which are capable of dealing with the most important aspects of these forming processes [9]. Nevertheless some important particular aspects of glass forming processes need a dedicated software or, at least, some user intervention in order to adapt commercial codes to this particular
S. Gre´goire et al. / Computers and Structures 85 (2007) 1194–1205
field, which may not be an easy task. One of these aspects is the modelling of the heat resistance between glass and mould at their interface which will be addressed here. The heat resistance at the interface is an average parameter thermally characterising the interface between the glass and the mould, which may be composed of air, lubricant or combustion products from the lubrication [10–14]. That parameter is a function of the interface pressure, of the glass temperature and therefore its viscosity, and of the asperities in the mould and its temperature. For the same pressure, the higher the temperature of the glass the smaller is the viscosity and therefore the physical contact is more effective as the glass fills better the micro cavities in the mould surface. When the glass cools against the mould walls then rapidly the viscosity rises near to a solid state and the glass may not even attain or may lose the contact with the mould as the pressures involved are usually very low [15,16]. When a simulation code, based on the finite element method, is used to deal with the processes described, typically glass and mould regions are modelled with nonmatching meshes [16,17]. Perfect mechanical and thermal contact is usually treated by a Lagrange multiplier technique or approximately by penalty formulations [17,18]. In particular this last approach may be used to model non-perfect heat contact by varying the penalty parameter, which, in fact, acts as an equivalent heat transfer coefficient between the two materials at the interface. Unfortunately for the user of these codes it is almost impossible to manipulate those parameters in order to make them directly or indirectly a function of the pressure, viscosity, temperature of the glass and mould or other local defects. These facts triggered the study presented here in which an interface finite element is implemented in the commercial code Abaqus in order to deal with heat transfer between glass and mould for non-matching meshes and where the heat transfer coefficient can be made, locally, a function of pressure, temperature or viscosity. The layout of this paper is the following: in Section 2 the governing equations in heat transfer problems are briefly introduced and the inclusion of the interface element justified; the definition of the interface element and its implementation in a commercial code are presented in Section 3 and a test to appraise its performance is carried out in Section 4. In the last section, the influence of the heat transfer coefficient on a hollow-ware glass article is studied through a simulation of a Champagne flute.
‘‘attached’’ to the material. In this case, for each time step, both in glass, Xg, and mould, Xm, Fig. 1, the governing differential equation is kr2 T þ Q ¼ 0
ð1Þ
where T is the temperature, Q is a heat source term and k the thermal conductivity of the medium. Regarding only the glass domain the natural boundary conditions may be of heat flux at the interface with the mould, in C1, and in the zones in contact with the surrounding air, in C2, i.e., oT ¼ h1 ðT T m Þ in C1 ð2Þ kg on1 oT kg ¼ h2 ðT T a Þ in C2 ð3Þ on2 where h1 and h2 are the corresponding heat transfer coefficients and Tm and Ta are, respectively, the mould temperature at the interface and the ambient temperature. Essential boundary conditions, i.e., prescribed temperatures, may exist in other part of the boundary, in C3, ð4Þ
T ¼ T in C3
It is possible to obtain an equivalent weak form for the problem by writing Z Z oT ðk g r2 T þ Qg ÞdT dXg þ k g h1 ðT T m Þ dT dC on1 Xg C1 Z oT þ k g h2 ðT T a Þ dT dC ¼ 0 ð5Þ on 2 C2 Performing an integrating by parts in the first term we obtain Z Z oT dT dC ðk g rT :rðdT Þ þ Qg dT Þ dXg þ kg on Xg C Z oT þ k g h1 ðT T m Þ dT dC on 1 C1 Z oT k g h2 ðT T a Þ dT dC ¼ 0 ð6Þ þ on2 C2 Noting that some flux terms cancel in boundaries C1 and C2, and that in boundary C3 the variation on T is zero then finally Eq. (6) may be written as
Ta
Ωm
n2 Γ2
2. Heat conduction on the glass: weak and strong forms. Interface element for natural boundary conditions In order to introduce the definition of the interface element let us state briefly the heat transfer problem in the process addressed. For simplicity let us assume, in the glass, that radiation effects are neglected and that, in our model, for each time step only diffusion takes place, i.e., the heat is convected by the moving Lagrangian mesh
1195
Ωg Γ3 T
Γ1
n1
Tm
T Γ2
n2 Ta Fig. 1. Glass and mould domains. Boundary conditions in the glass.
S. Gre´goire et al. / Computers and Structures 85 (2007) 1194–1205
1196
Z
ðk g rT rðdT Þ þ Qg dT Þ dXg Xg Z ðh2 ðT T a ÞÞdT dC ¼ 0
Z
ðh1 ðT T m ÞÞdT dC C1
ð7Þ
C2
This last equation may be viewed as a variational principle dPg ¼ 0
ð8Þ
associated with the functional Z 1 1 2 2 k g krT k Qg T dXg þ h1 ðT T m Þ dC Xg 2 C1 2 Z 1 h2 ðT T a Þ2 dC þ ð9Þ 2 C2
Pg ¼
Z
transfer coefficient assumes the value of h2. When contact is established the boundary condition is represented by Eq. (2). The heat coefficient h1 may then be defined, at the sampling integrating points of the interface element, as a function of glass viscosity, temperature or pressure which is known at the nodes of the underlying element. Likewise the mould temperature Tm at each contacting node may be obtained by determining its position, in terms of natural coordinates, on the side of element on the mould and subsequently using the shape functions to interpolate its value from the mould node temperatures. The details of the implementation of this interface element are discussed in the following sections. 3. Thermal interface element 3.1. FEM formulation
If the thermal contact between glass and mould would be perfect, corresponding to an infinite value of the heat transfer coefficient at the interface, then the temperatures at the two sides of the interface would be the same at each point. The second term of the right hand side of the functional of Eq. (9) would then be viewed as a penalty term imposing equality of temperatures on both sides. The parameter h1 would then act as a penalty parameter assuming a large number corresponding to a big heat transfer coefficient. If the contact is not perfect and some interface thermal resistance is present then the parameter h1 would be smaller, corresponding to a local equivalent heat transfer coefficient, and resulting in different temperatures in the glass and the mould at the interface. Based on this evidence an interface element is introduced to specially deal with that term on the functional that corresponds, on the strong form, to Eq. (2). On the mesh representing the glass, Fig. 2, an element boundary that is or may enter in contact with the mould is assumed to be insulated. An interface element is then superposed on that boundary and takes care of the heat transfer to the air and to the moulds. If there is no contact the heat
The finite element formulation may be obtained from the weak form in Eq. (7). In particular in what respects the interface element the discretisation is obtained from the term Z ... ðh1 ðT T m ÞÞdT dC . . . ð10Þ C1
The contribution of each interface element to the discretised form of (10) is . . . dT e H e T e dT e f e . . . e
ð11Þ e
in which T are the element nodal temperatures, H is the element thermal ‘‘stiffness’’ matrix and f e the element contribution for the overall thermal ‘‘force’’ vector. Theses matrices are obtained from the element shape functions using the standard procedure for isoparametric axisymmetric elements (see Fig. 3 for details) by performing the integrations in the natural coordinates as Z 1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi e T H ¼ 2p h1 ðnÞrðnÞN N ðox=onÞ2 þ ðoy=onÞ2 dn ð12Þ 1 Z 1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 e T ðox=onÞ þ ðoy=onÞ dn ð13Þ f ¼ 2p h1 ðnÞrðnÞT m N 1
3.2. Implementation of the interface element on a FEM code
Axisymetric axis
Glass contact surface
Interface elements
Mould contact surface
Fig. 2. Positioning of the interface elements.
3.2.1. Non-matching meshes As referred above typically non-matching meshes will be used for the glass and the mould. When a node of the mesh of the glass enters in contact with an element of the mould, the temperature at the mould in the contact point must be determined to establish the thermal contact. Since the mechanical contact is not perfect, as it is dealt with penalty functions, the contact nodes may be slightly outside or inside the element of the mould mesh. A safe procedure must be then performed to assure that the thermal contact is made. In fact the same procedure is also adopted whenever a remeshing operation is per-
S. Gre´goire et al. / Computers and Structures 85 (2007) 1194–1205
Nj : (1+ξ)/2
Ni : (1-ξ)/2 i
ξ1(-1/ √3)
ξ
j
ξ2(1/ √3) Y
i
Xi, Yi, hi, Ti
r(ξ)
dL
∂y ∂ξ
∂x ∂ξ
2
+
∂y ∂ξ
2
∂x ∂ξ
Lij
j
dL =
X
Xj, Yj, hj, Tj dS=2πr(ξ)dL
Fig. 3. Interface element: (a) in natural coordinates, (b) in global coordinates.
formed, as some nodes of the new mesh may slightly fall outside the mesh of the previous domain. A gap function is defined as 2
/ðn; gÞ ¼ ðxem ðn; gÞ xg Þ þ ðy em ðn; gÞ y g Þ
2
ð14Þ
The coordinates xem ; y em to where the temperatures of the mould element are interpolated (extrapolated) should then minimize the function /ðn; gÞ. This implies that the gradient of /ðn; gÞ should be zero, i.e., r/ðn; gÞ ¼ O
ð15Þ
The solution is obtained in an iterative manner using the Newton–Raphson method. Departing from ðn0 ; g0 Þ ¼ ð0; 0Þ at iteration n the following system is solved to obtain ðDnn ; Dgn Þ ( Dðo/=onÞnþ1 ¼ ðo2 /=on2 ÞDnn þ ðo2 /=onogÞDgn ð16Þ Dðo/=ogÞnþ1 ¼ ðo2 /=ogonÞDnn þ ðo2 /=og2 ÞDgn The natural coordinates are then updated and convergence criterion is checked by defining a very small tolerance (toler ’ 104) ðnnþ1 ; gnþ1 Þ ¼ ðnn ; gn Þ þ ðDnn ; Dgn Þ kðnnþ1 ; gnþ1 Þ ðnn ; gn Þk=kðnn ; gn Þk < toler
ð17Þ ð18Þ
The convergence is typically reached in two or three iterations. The temperature at the contact point, is then determined by using the shape functions of the element Tm ¼
nX odes
N i ðnc ; gc ÞT i
ð19Þ
i¼1
where the subscript c stands for converged value. This definition of Tm should then be included in Eq. (13) coupling the degrees of freedom in the glass and the mould and therefore allowing the heat transfer between the two materials to take place.
1197
3.2.2. Implementation in a commercial code One would like to include the features described in a commercial code where many facilities are available like different constitutive equations, element type choice, contact handling, remeshing operations and the usual pre and post processing tools. In this work the Abaqus commercial code was selected as the host code for the implementation of the proposed element, via the user subroutine UEL. Nevertheless part of the procedure indicated in the previous section, which should be pursued when dealing with a home code where the access is free, will be a difficult (impossible?) task to achieve when a commercial code is used. This is the case with the definition of Tm in Eq. (13). Whilst the node temperatures in the interface elements, as they are permanently ‘‘attached’’ to the same elements of the glass mesh, can be made equal to the corresponding glass temperatures using multiple point constraints, the definition of Tm as a function of the temperatures of the mould elements is a hard and cumbersome task in this case. This is due not only to the fact that, specially in the press operation, the contacting nodes are moving from one element to another in the contact surface but also that it is not possible to obtain internally the number (label) of the element of the mould where contact has been established. This can only be done indirectly using the procedure referred in the last section. Therefore a simplified staggered algorithm was established. The value of Tm in Eq. (13) is evaluated, in each iteration, according to Eq. (19). The glass cooling against the mould is, by this process, guaranteed. However, because the temperatures in the nodes of the mould mesh are not constrained internally with glass temperature values, an equivalent heat transfer coefficient is assumed at each contacting side of an element of the mould. Its value is a weighted average of the values of the heat transfer coefficients of the surrounding contacting nodes of the glass. 4. Assessment of the performance of the interface element A set of simple tests was carried out in order to validate the performance of the interface element and its implementation in the Abaqus code. In those tests, the solutions obtained with standard thermal and thermomechanical models in Abaqus are compared with those obtained with a thermomechanical model in Abaqus that includes the interface element. The main objective is to assess the heat transfer process between two non-matching meshes, which is present in glass forming simulation due to the evolving contact between glass and mould meshes. This is an important feature because the glass viscosity is strongly dependent on the temperature. 4.1. Tests description The axisymetric model in Fig. 4 represents two cylinders of two different materials, a crystal glass and a steel (XC38), initially put into contact. Each cylinder has a
S. Gre´goire et al. / Computers and Structures 85 (2007) 1194–1205
1198 Z 30mm Glass (1150˚C)
10mm r
nodes at the steel interface (Fig. 5b). The mesh for the steel cylinder is identical for both simulations (7 elements in Zdirection and 16 in R-direction). The glass has 7 elements in the Z-direction and is meshed by 16 elements for the matching meshes simulation and 21 elements for the nonmatching meshes in the R-direction (Fig. 5).
10mm Steel (485˚C)
4.2. Simulation with matching meshes
Fig. 4. Axisymetric model representing a glass and a steel cylinder into contact.
radius of 30 mm and a height of 10 mm. The initial temperatures are 1150 C on the glass and 485 C on the steel. The thermal properties of the two materials at these temperatures are given in Table 1. By using the interface element it would be possible to include at the glass/steel interface a heat transfer coefficient, simulating the heat resistance, dependant on the glass and the mould contact temperatures, on the contact pressure or in the local viscosity. However when Abaqus code is used only contact pressure or gap thickness at the interface can influence the heat transfer coefficient. Therefore, in order to compare the results with and without the interface element, a constant heat transfer coefficient of 5000 W/m2/ C was assumed at the interface. There is only heat transfer between the cylinders as they are insulated at the other faces. Only the self weight of the cylinders is considered as mechanical force. The friction coefficient between the two materials is set equal to the unity. The heat transfer process is analysed during an initial short period of 10 s. Two different meshes for the glass cylinder are used: one that matches the mould mesh at the interface (Fig. 5a) and a finer mesh with non-matching Table 1 Thermal properties of glass and steel at, respectively, 1150 C and 485 C
Thermal conductivity (W/m/C) Specific heat (J/kg/C) Density (kg/m3)
Glass (crystal) at 1150 C
Steel (XC38) at 485 C
2.77 1282 2500
37 620 7800
Fig. 5. Different used meshes: (a) matching meshes, (b) Non-matching meshes.
In the thermomechanical model with the interface element, Fig. 6a and in the standard thermal simulation Fig. 6c the same temperatures are obtained, with average differences of less than 0.8 C for all the length of the analysis. The temperature isocurves are straight lines parallel to the interface. But in the standard thermomechanical model, i.e., Abaqus code with no interface element, the temperatures are not identical along the steel contact surface (Fig. 6b) as the temperature isocurves depart slightly from parallel lines to the interface, specially closer to the axisymmetry axis where the temperatures are higher than the previous simulations. The evolution of the temperatures at the extreme points of the interface for both the glass and the steel sides are represented in Figs. 7 and 8. The temperature evolutions between the contact element model and the thermomechanical model are different on the left sampling point located at the steel contact surface. After 0.4 s, a difference of 12 C is noticed (Fig. 8a). At the same time, the same temperature difference is read at the left glass point (Fig. 7a). On the right comparison points, the temperatures evolutions of these three models are similar (Figs. 7b and 8b). 4.3. Results for non-matching meshes simulation In this case, both in the standard Abaqus thermal and thermomechanical models, the solution is perturbed due to the fact that the two meshes do not match. In Fig. 9b and c, it is clear that in these cases the temperature isocurves depart from parallel lines to the interface, especially near the axis. When the interface element is included in the analysis, these perturbations are eliminated and the correct solution is recovered (Figs. 9a and 10). The comparison of the temperature profiles at the interface are shown in Fig. 10 where its clearly seen how the perturbations on the temperature are eliminated with the inclusion of the proposed element. The evolution of the temperatures at the extreme points of the interface for both the steel and the glass sides are represented in Figs. 11 and 12. The temperature evolutions are similar for the heat transfer model and the thermomechanical (maximal difference of 4 C after 0.6 s) but they differ from the solution with the interface element, specially on the steel side (Fig. 11a) where differences as high as 24 C on the steel side are obtained in this simple example. It is known in glass forming how important is the temperature distribution in the mould for the final product [1–6].
S. Gre´goire et al. / Computers and Structures 85 (2007) 1194–1205
1199
Fig. 6. Temperature isocurves inside the steel at time 10 s: (a) contact element model, (b) thermomechanical model, (c) heat transfer model.
1150
1150
1100
Contact element model Thermomechanical model Heat transfer model
Temperature (°C)
1050 1000 950 900 850 800
Glass (1150˚C)
1050 Temperature (°C)
1100
1000 950 900 850 800 750
750 700
Glass node
650
Steel (485˚C)
700
Glass node
650 600
600 0
2
4
6
8
0
10
2
4
6
8
10
Contact time (s)
Contact time (s)
Fig. 7. Temperature evolution at two points on the glass interface for the three models (matching meshes).
620
620
Contact element model Thermomechanical model Heat transfer model
580 560 540
Glass (1150˚C)
520 Steel node
500
Steel (485˚C)
600 Temperature (°C)
Temperature (°C)
600
Steel node
580 560 540 520 500 480
480 0
2
4 6 8 Contact time (s)
10
0
2
4
6
8
10
Contact time (s)
Fig. 8. Temperature evolution at two points on the steel interface for the three models (matching meshes).
Fig. 9. Temperature isocurves inside the steel after 10 s of simulation: (a) contact element model, (b) thermomechanical model, (c) heat transfer model.
In the glass side also differences of 14 C are detected after 0.6 s (Fig. 12a) what may correspond, in the interval of
forming temperatures, to changes in the glass viscosity that may range from 20% to 500%.
S. Gre´goire et al. / Computers and Structures 85 (2007) 1194–1205
1200 680
Contact element model Thermomechanical model Heat transfer model
Glass interface
Temperature (˚C)
660
640 Glass
620
D
Steel Mould interface
600 0
5
10
15
20
25
30
Distance D from the axisymetric axis (mm)
Fig. 10. Glass and steel nodal temperature along the glass/steel interface after a run of 10 s.
620
620
Contact element model Thermomechanical model Heat transfer model
580
600 Temperature (°C)
Temperature (°C)
600
560 540 Glass (1150°C)
520
580 560 540 520
Steel node
500
Steel (485°C)
Steel node
500 480
480 2
0
8 6 4 Contact time (s)
0
10
2
8 6 4 Contact time (s)
10
Fig. 11. Temperature evolution at two points on the steel interface for the three models (non-matching meshes).
1150
1150
Contact element model Thermomechanical model Heat transfer model
1050 Temperature (°C)
1000 950 900 850 800
1100 1050 Temperature (°C)
1100
Glass (1150°C)
1000 950 900 850 800 750
750 700
Glass node
650
Steel (485°C)
Glass node
700 650 600
600 0
2
4 6 8 Contact time (s)
10
0
2
4 6 8 Contact time (s)
10
Fig. 12. Temperature evolution at two points on the glass interface for the three models (non-matching meshes).
4.4. Comments The results of the simple test described in the previous sections showed that the Abaqus standard thermal an ther-
mal–mechanics models may introduce some perturbations in the temperature distribution when non-matching meshes are used, which is probably due to the penalised form used to impose thermal and mechanical contact. This is a crucial
S. Gre´goire et al. / Computers and Structures 85 (2007) 1194–1205
1201
issue for glass forming process simulation due to the fact that glass viscosity is highly dependent on the temperature (see Eq. (20)). The inclusion of the interface element, with the initial purpose of making it possible to introduce the dependence of heat resistance on process parameters, proved to be also an efficient tool to eliminate those perturbations making the numerical model more robust and reliable.
with the thermomechanical axisymetric model on the Abaqus commercial code, in which the interface thermal element is included. The total forming time for the Champagne flute is divided in three phases:
5. Application to a press and blow process
The initial temperature of the plunger, the blank mould and the blow mould is 485 C and the initial temperature of the glass is 1150 C. The thermal conductivity and the specific heat of the XC38 steel used for the moulds and the punch at 600 C are, respectively, 33.9 W/m/C and 708 J/kg/C. At 1000 C the thermal conductivity and the specific heat of the crystal are, respectively, 2.5 W/m/C and 1265 J/ kg/C. The temperature dependences of the thermal conductivity and specific heat of the crystal and of the XC38 steel are presented in Figs. 15 and 16. Glass viscosity is highly dependent on the temperature and is commonly represented, in the forming range temperatures, by the Vogel, Fulcher and Thalmann law [21]. For the particular case of the crystal used in the Champagne flute this law takes the following form:
The forming of an upper part of a Champagne flute, requiring three sequential forming phases – press, gravity stretch and blow – is studied in this section. This forming process, and particularly some of its forming parameters, has already been studied by Noiret [19]. The study proposed here concerns the analysis of the influence of the heat transfer coefficient at the interface during the full process. Final temperature and thickness distributions will be the important parameters to look at carefully. 5.1. Experimental data In order to obtain representative heat transfer coefficients at the interface, an experimental device was designed [14,20] which puts into contact a crystal glass at an initial temperature of 1150 C and a XC38 steel plunger at an initial temperature of 485 C. A set of heat transfer coefficients were obtained for different increasing contact pressures and temperatures. Two of these cases are shown in Fig. 13 where it is clear how the pressure may influence the heat transfer value.
• a pressing time: 0.3 s (Fig. 14a and b); • a gravity stretching time: 4 s (Fig. 14c and d); • a blowing time: 7.5 s (Fig. 14d and e).
logðgÞ ¼ 1:9458 þ
5407 Pa s ðT 174:5Þ
ð20Þ
in which T is the temperature and g is the viscosity. In Fig. 17 this law is represented in the temperature range of [500 C, 1200 C].
5.2. Modelling of the forming processes In order assess the pressure influence two simulations of the full forming process, using the two heat transfer coefficient evolutions represented in Fig. 13, were performed
Heat transfer coefficient (W/m²/°C)
8000
6000
h1I ( t ) 4000
h1II ( t )
2000 0
1
2 3 4 Contact time (s)
5
6
Fig. 13. Heat transfer coefficient evolution obtained by experimentation. hI1 ðtÞ: increasing pressure from 0 to 3.0 MPa, hII 1 ðtÞ: increasing pressure from 0 to 1.5 MPa.
Fig. 14. Champagne flute forming phases (thermal isocurves). Scale 1:2. (a) Initial position; (b) end of pressing; (c) before stretching; (d) end of stretching; (e) end of blowing.
S. Gre´goire et al. / Computers and Structures 85 (2007) 1194–1205
1202
0
0.05
0.1
0.15
0.2
0.25
0.3
0
45 40
-0.01
35 moulds and plunger : XC38 steel
30 25
glass : cristal
20 15 10
Displacement (m)
Thermal conductivity (W.m-1.°C-1)
50
-0.02 -0.03 -0.04
5
-0.05
0 400
600 800 1000 Temperature (˚C)
1200
-0.06
Fig. 15. Thermal conductivity of glass, moulds and plunger.
Time (s)
Fig. 18. Plunger vertical displacement versus time.
• a heat transfer coefficient (as given in Fig. 13) in case of contact between glass and mould or glass and plunger.
1400
-1
Specific heat (J.kg-1.°C )
1200 1000 800
moulds and plunger : XC38 steel
600
glass : cristal
400 200 0 400
600 800 1000 Temperature (°C)
1200
Fig. 16. Specific heat of glass, moulds and plunger.
16
For this stage four remeshing operations were needed to obtain the final parison. The subsequent stage is a stretching operation of the glass part due to gravity only in which a certain homogenisation of the temperatures in the glass is obtained. Prior to begin this stage a new remeshing operation is performed due to the mesh distortion in the final period of the pressing stage. The thermal boundary conditions for the stretching stage are: • a convection with the ambient air (h2 = 20 W/m/C and Ta = 50 C) where there is no contact between glass and mould contact; • a heat transfer coefficient (as given in Fig. 13) in case of contact between glass and mould.
14
Log(η) (Pa.s)
12 10 8 6 4 2 0 500
600
700 800 900 1000 1100 1200 Temperature (°C)
Fig. 17. Temperature dependence of the crystal viscosity.
In the pressing stage, the parison (a hollow-shaped blank) is obtained after a displacement of 51 mm of the plunger (Fig. 18). The thermal boundary conditions are: • a convection with the ambient air (h2 = 20 W/m/C and Ta = 50 C) in the case of non-contact between the glass and mould or glass and plunger;
The mesh of the blowing phase is identical to the mesh of the stretching phase. In the previous study of Noiret [19] the blowing pressure was set to 1 kPa for a NBS 710 standard glass. The crystal glass adopted here has a much higher viscosity at the forming temperatures and therefore to perform the final blowing another pressure law, represented in Fig. 19, was used. The flute forming is completed after a blowing time between 2 s and 4 s (depending on the heat transfer coefficient, i.e., on the heat resistance at the interface). After full contact some extra time is needed to cool the glass in order to remove it out of the mould without damaging the final shape. The boundary conditions for the stretching phase are: • a convection with the pressurised air (h2 = 40 W/m/C and Ta = 200 C); • a convection with the ambient air (h2 = 20 W/m/C and Ta = 50 C) if contact has not been established; • a heat transfer at the glass and mould interface after contact has been established coefficient with the appropriate heat transfer coefficient h1(t).
S. Gre´goire et al. / Computers and Structures 85 (2007) 1194–1205
But, according to the same figure, all the temperatures in the neck have dropped significantly, which is explained by the thinning effect in this zone caused by the pressing phase. During the stretching phase, in both simulations, a homogenization of the temperatures in glass is noticed (Fig. 21b) but, as expected, at the end of this phase the final temperatures are lower for the case of the simulation with hI1 ðtÞ (Fig. 20). For the same time operations when the blowing phase begins the stiffness of the glass in the two cases will be different due to high dependence of the viscosity on the temperature. Therefore the contact between glass and mould begins 1.2 s later for the case of the simulation with hI1 ðtÞ due to the glass viscosity which is higher owing to lower temperatures. In both simulations, the temperatures level reached, in the end, at the mould’s sampling point have the same values but with an offset time of 1.2 s.
25
Air pressure (kPa)
20 15 10 5 0 0
1
2
3 4 5 Blowing time (s)
6
1203
7
Fig. 19. Evolution of the blowing pressure versus time.
5.3. Results and discussion In this section, by utilising the two different heat transfer coefficients described above, hI1 ðtÞ and hII 1 ðtÞ, the influence of the heat transfer coefficient on the temperature evolution inside the glass and the moulds and on the final thickness of the glass part is analysed. 5.3.1. Temperature field and stretching During the pressing, as expected, the glass surface cools more in the simulation with hI1 ðtÞ and, correspondingly, in the mould the temperature increases more. This is illustrated in Fig. 20 where the temperature evolutions, for two particular points in the glass and mould that will contact during the blow phase, are represented. In both simulations, no significant temperature variation in the heart of the glass in the lower zone of parison is noticed, particularly in the simulation with hI1 ðtÞ, as shown in Fig. 21a.
5.3.2. Thickness distribution The glass temperatures after the pressing phase are lower in the simulation with hI1 ðtÞ. Therefore the glass deformation during the gravity stretching is less important. In this case, the sampling node in glass bottom has a displacement of 8 mm after 4 s of stretching compared to a displacement of 13.5 mm at the same point in the simulation with hII 1 ðtÞ (Fig. 22). In Figs. 23 and 24 the thickness of the Champagne flute is represented and measured from the internal side of the flute along its height. In the simulation with hI1 ðtÞ, the final thickness is relatively homogeneous (an average of 2.4 mm). This can be explained by the fact that the relatively low temperatures in the neck restrict the deformations in that zone during the stretching phase. Indeed, at the end of the pressing,
1150 Stretching
Blowing
1050
Temperature (˚C)
950 Glass 850 750 650 550
Mould
450 0
2
4
6
8
10
12
punch + glass + blank mould
glass + blow mould
Process time (s)
Fig. 20. Temperature evolution on a glass and on a mould node on the contact surfaces during the full process. Thick curve: simulation with hI1 ðtÞ. Thin curve: simulation with hII 1 ðtÞ.
S. Gre´goire et al. / Computers and Structures 85 (2007) 1194–1205
1204
Fig. 21. Temperature field in the glass during the forming process in the simulation with hI1 ðtÞ. (a) End of pressing, (b) end of stretching, (c) end of blowing.
0
Displacement (m)
-0.002
h (t ) I 1
-0.004 -0.006
h1II ( t )
-0.008 -0.01 -0.012 -0.014 0
1
2 Creeptime (s)
3
4
Fig. 22. Displacement of a node in the glass part during stretching phase in the two simulations.
120mm 100
80
h1I ( t )
60
40
h1II ( t ) 20
0mm 0 1
2
3 4 Thickness (mm)
5
Fig. 24. Thickness of the Champagne flute at the end of the press–blow process.
120 Distance from the bottom of the glass (mm)
Distance from the bottom of the glass (mm)
120
100
100mm
80
60
40
0mm
20
0 2
4 6 Thickness (mm)
8
Fig. 23. Thickness of the Champagne flute at the end of the pressing.
that zone has a thickness that varies between 3 mm and 4 mm (Fig. 23) and at the end of the full process, the same zone has an average thickness of 2.5 mm. Most of the deformation is localized in the lower zone of parison. In that zone the thicknesses are varying between 5 mm and 7.5 mm after the pressing and are dropping to an average of 2.3 mm. For the simulation with the hII 1 ðtÞ, the glass is thinner in the neck of the flute (value of 2.1 mm) than in the lower zone of parison (4.5 mm). Here, contrary to the previous simulation, the glass is softer (low viscosity) in the neck and allows a bigger stretching. Then, the bottom zone keeps more or less the same thicknesses which decrease from values comprised between 5 mm and 7.5 mm to values comprised between 3 and 4.5 mm.
S. Gre´goire et al. / Computers and Structures 85 (2007) 1194–1205
6. Conclusions An interface finite element was implemented, at the University of Porto, in the commercial code Abaqus in order to deal with heat transfer between glass and mould for nonmatching meshes. The associated heat transfer coefficient can be, locally, set as a function of pressure, temperature or viscosity. The interface finite element was validated on a simple benchmark. It proved to perform well and to regularise some perturbations introduced by the penalised form adopted in the commercial code to deal with non-matching meshes. This is an important feature as the glass viscosity and therefore the final shape are strongly dependent on the temperature distribution. An application on the forming of a Champagne flute by a press and blow operation illustrates the influence of the heat transfer coefficient throughout the whole process. Two simulations are performed with two heat transfer coefficients obtained experimentally at the University of Valenciennes showing a great influence in the final product in terms of glass temperature values and distribution. These differences have a great impact on the glass viscosity and therefore on the forming operation times and conditions. The time for each phase (pressing, stretching and blowing), the value of the pressure at any instant and when it should be activated play a decisive role in the final shape and thicknesses distribution. To succeed the optimization of the quality of the manufactured glass products, the numerical thermomechanical models require a perfect knowledge of the heat exchange at the glass/mould interface. The procedure described here may play an important role in the robustness of the modelling of glass forming processes. Acknowledgements The present work is supported by the Department of Mechanical Engineering and Industrial Management of the University of Porto (Portugal), the Laboratoire d’Automatique de Me´canique et d’Informatique Industrielles et Humaines de l’Universite´ de Valenciennes et du Hainaut-Cambre´sis, the CNRS and the Ministe`re de l’Education Nationale et de la Recherche. The authors gratefully acknowledge the support of these institutions. References [1] Williams JH, Owen DRJ, Cesar de Sa JMA. The numerical modelling of glass forming process. In: Proceedings of the XIVth international congress on glass, New Delhi, India; 1986.
1205
[2] Ce´sar de Sa´ JMA. Numerical modelling of glass forming processes. Eng Comput 1986;3(6):266–75. [3] Storck K. Thermal system analysis. Heat transfer in glass forming and fluid temperature-control systems. Institute of Technology, Linko¨pings Universitet. Linko¨pings Studies in Science and Technology. Dissertation no. 542. Linko¨ping; 1998. [4] Lochegnies D, Mare´chal C. How glass properties can affect the glass flow? In: Proceedings of the 2nd international colloquium modelling of glass forming and tempering. Valenciennes, France: Presses Universitaires de Valenciennes; 2002. [5] Ce´sar de Sa´ JMA, Natal-Jorge RM, Cardoso RPR. Simulation of forming of hollow glass containers. In: Proceedings of the Euromech 388 colloqium on modelling of glass forming processes. Presses Universitaires de Valenciennes. Lez Valenciennes no. 27. Valenciennes, France; 1999. [6] Ce´sar de Sa´ JMA, Natal-Jorge RM, Silva CMC, Cardoso RPR. A computacional model for glass container forming processes. In: Wunderlich W, editor. Proceedings of ECCM’99 – European conference on computational mechanics, Munich, Germany; 1999. [7] Lochegnies D, Moreau P, Guilbaut R. A reverse engineering approach to the design of the blank mould for the glass and blow process. Glass Technol 2005;46(2):116–20. [8] Mare´chal C, Moreau P, Lochegnies D. Numerical optimization of a new robotized glass blowing process. Eng Comput 2004;19:233–40. [9] Hyre M, Rubin Y. Comprehensive modelling of the glass container forming process. In: Ken-ichiro Mori, editor. Simulation of materials processing: theory, methods and applications. In: Proceedings of the 7th international conference on numerical methods in industrial forming processes-Numiform 2001, Toyohashi, Japan; 2001. [10] Fellows CJ, Shaw F. A laboratory investigation of glass to mould heat transfer during pressing. Glass Technol 1978;19(1):4–9. [11] Pchelyakov SK, Guloyan YA. Heat transfer at the glass mould interface. Glass Ceram 1985;42:400–3. [12] Storck K, Loyd D, Augustsson A. Heat transfer modelling of the parison forming in glass manufacturing. Glass Technol 1998;39(6): 210–6. [13] Keijts G, van der Werff K. Heat transfer in glass container production during the final blow. Glass Technol 2001;42(3):104–8. [14] Moreau P, Gre´goire S, Lochegnies D. Me´thodologie expe´rimentale pour la de´termination du coefficient de transfert thermique a` l’interface verre/outil. Verre 2003;9(2):24–7. [15] McGraw DA. Transfer of heat in glass during forming. J Am Ceram Soc 1961;44(7):353–63. [16] Mori M, Shimazaki T, Naraki T, Wakatsuki H. Numerical simulation system for mould designer. In: Proceedings of the 2nd international colloquium modelling of glass forming and tempering. Valenciennes, France: Presses Universitaires de Valenciennes; 2002. [17] Ce´sar de Sa´ JMA, Natal-Jorge RM, Cardoso RPR. Improved finite element analysis on blow forming. In: Proceedings of the 2nd international colloquium modelling of glass forming and tempering. Valenciennes, France: Presses Universitaires de Valenciennes; 2002. [18] Wriggers P. Computational contact mechanics. England: John Wiley & Sons; 2002. [19] Noiret C, Thibaud C, Lochegnies D, Oudin J. Le formage du verre creux: l’approche nume´rique sur deux exemples. Verre 1996;2:19–22. [20] Gre´goire S, Moreau P, Cesar de Sa JMA. Integration of the heat transfer coefficient h(t) in glass forming modelling. In: Proceedings of GlassDays 2005, Kaiserslautern, Germany; 2005. [21] Fulcher GS. Analysis of recent measurements of the viscosity of glasses. J Am Ceram Soc 1925;8(6):339–55.