Dr Norbert Siedow
email: norbert.siedow@itwm.fraunhofer.de

Dominique Lochegnies

Fabien Béchet

Philippe Moreau

Hiroshi Wakatsuki

Nobuhiro Inoue

Axisymmetric modeling of the thermal cooling, including radiation, of a circular glass disk

Keywords:

teaching and research institutions in France or abroad, or from public or private research centers.

Introduction

Proper glass cooling is very important to achieve desired product quality. If cooling does not occur properly, undesired stresses may occur inside the glass, which could lead to glass breakage either during the cooling itself or during subsequent product manipulations. Numerical simulations can be used to study the physical behavior of glass (temperature, stresses) during cooling to improve cooling process design. Modeling glass cooling is a complex, non-linear thermo-mechanical problem. In the last decades, glass-cooling models have been widely developed [START_REF]Modelling of Glass Forming and Tempering[END_REF][START_REF]Modelling and Control of Glass Forming and Tempering[END_REF][START_REF] Siedow | Application of a new method for radiative heat transfer to flat glass tempering[END_REF][START_REF] Nielsen | Finite element implementation of a glass tempering model in three dimensions[END_REF] using commercial software packages or homemade codes.

Glass is a semi-transparent material. Consequently, in addition to heat conduction and heat convection, radiation plays an important role in thermal exchanges, especially at high temperatures where it is the dominant heat transfer process. Thorough assessments of radiative heat transfer can be found in [START_REF] Modest | Radiative heat transfer[END_REF], [START_REF] Siegel | Thermal radiation heat transfer[END_REF], and [START_REF] Viskanta | Heat transfer in semitransparent solids[END_REF], while the application of radiative heat transfer to the glass industry is assessed in [START_REF] Krause | Mathematical simulation in glass technology, Schott series on glass and glass ceramics[END_REF] and [START_REF] Farina | Mathematical models in the manufacturing of glass[END_REF].

In the literature, different solutions were proposed to model glass cooling and account for radiation effects. The simplest one is to completely ignore radiation [START_REF] Nielsen | Finite element implementation of a glass tempering model in three dimensions[END_REF][START_REF] Chen | Numerical Simulation and Experimental Stud of Residual Stresses in Compression Molding of Precision Glass Optical Components[END_REF]. Another solution is to consider surface radiation only by applying Stefan-Boltzmann's law [START_REF] Sellier | An iterative algorithm for optimal mold design in high-precision compression moulding[END_REF], which is appropriate for opaque materials like metals. However, since glass is a semi-transparent material, internal radiation cannot be ignored. A widely used solution involves treating radiation as a correction of heat conduction by using an equivalent conductivity (such as the active thermal conductivity method [START_REF] Fotheringham | Active thermal conductivity of hot glass[END_REF][START_REF] Mann | Determination of specific heat and true thermal conductivity of glass from dynamic temperature data[END_REF]) or the Rosseland approximation [START_REF] Dusserre | Thermo-mechanical stresses in cast steel dies during glass pressing process[END_REF]. These methods are fast and simple to integrate into commercial software packages. Originally derived by S. Rosseland in 1924 for stellar radiation [START_REF] Rosseland | Note on the absorption of radiation within a star[END_REF], the Rosseland approximation is only valid, however, for optically thick glass, i.e. 𝑑 • 𝜅(𝜆) ≫ 1, where 𝜅(𝜆) denotes the wavelength depending absorption coefficient and 𝑑 the distance to the boundary. In [START_REF] Siedow | Application of a new method for radiative heat transfer to flat glass tempering[END_REF] and [START_REF] Béchet | Two-dimensional finite element modeling of glass forming and tempering processes, including radiative effects[END_REF] it was shown that the use of the Rosseland approximation for glass cooling could lead to significant stress calculation error.

The right way to model thermal radiation is to use the radiative transfer equation (RTE):

Ω ̅ • ∇𝐼(𝑥̅ , Ω ̅ , 𝜆, 𝑡) + 𝜅(𝜆)𝐼(𝑥̅ , Ω ̅ , 𝜆, 𝑡) = 𝜅(𝜆)𝐵(𝑇(𝑥̅ , 𝑡), 𝜆).

This equation is non-linear and high-dimensional regarding the spectral radiative intensity 𝐼(𝑥̅ , Ω, ̅ 𝜆, 𝑡), and therefore very time consuming to solve. Lee and Viskanta [START_REF] Lee | Transient conductive-radiative cooling of an optical quality glass disk[END_REF] used the Discrete Ordinate Method (DOM) in axisymmetric cylindrical coordinates to solve this equation in the case of the cooling of an optical-quality glass disk by natural cooling. These researchers suggested surrounding the disk with air at a constant temperature. To validate the solution method, they modified the boundary condition to obtain a one-dimensional solution from the two-dimensional formulation and compared it with experimental data from Field and Viskanta [START_REF] Field | Measurement and Prediction of the Dynamic Temperature Distributions in Soda-Lime Glass Plates[END_REF]. The difference between the simulated temperature and the experimental data is quite small. As an alternative to the time consuming DOM a fast and sufficiently accurate method based on the formal solution of the radiative transfer equation was developed in [START_REF] Siedow | Application of a new method for radiative heat transfer to flat glass tempering[END_REF]. This method is used here to simulate the cooling of a glass disk supported on its edge by a metal support mold. Due to the contact between glass and metal, a boundary condition that describes the exchange of radiative energy not only in the opaque wavelength region but also in the semitransparent region must be taken into account. The selected numerical solution method [START_REF] Siedow | Application of a new method for radiative heat transfer to flat glass tempering[END_REF] must be modified to incorporate this kind of boundary condition.

The paper is organized as follows. In section 2, the axisymmetric model considered in this paper is defined. The geometries of the disk and the mold, the mechanical and thermal equations and the radiative heat transfer model are discussed. To validate this model, the selected numerical solution is compared with values from [START_REF] Lee | Transient conductive-radiative cooling of an optical quality glass disk[END_REF] and [START_REF] Field | Measurement and Prediction of the Dynamic Temperature Distributions in Soda-Lime Glass Plates[END_REF] in a specific one-dimensional solution. It turns out that the proposed numerical method for radiative transfer is very close to the experimental data. In section 3, the method described in section 2 is applied to the axisymmetric problem of the cooling of the glass disk supported on a metal mold. The results obtained in terms of temperatures and stresses will be discussed. It will be shown that considering internal radiation is very important for simulating glass cooling.

Definition of the two-dimensional glass-cooling model 2.1 The geometric model for glass cooling

This paper examines the cooling of a circular disk supported on its edge by a metal mold as described in Figure 1. The glass domain is defined by 𝐷 𝑔 = {0 ≤ 𝑟 ≤ 𝑅, 0 ≤ 𝑧 ≤ 𝑒}. The mold domain is defined by 𝐷 𝑚 = {𝑅 1 ≤ 𝑟 ≤ 𝑅 2 , -𝐸 ≤ 𝑧 ≤ 0}. Since cooling occurs between time 0 and time 𝑡 𝑚𝑎𝑥 , 𝐷 𝑡 𝑔 = 𝐷 𝑔 × {0 ≤ 𝑡 ≤ 𝑡 𝑚𝑎𝑥 } and 𝐷 𝑡 𝑚 = 𝐷 𝑚 × {0 ≤ 𝑡 ≤ 𝑡 𝑚𝑎𝑥 } respectively denote the domains occupied by the glass and the mold over the time. The boundaries of these two domains are denoted by 𝜕𝐷 𝑡 𝑔 and 𝜕𝐷 𝑡 𝑚 respectively. The contact area between the glass and the mold is denoted by 𝜕𝐷 𝑡 𝑔𝑚 .

Figure 1: Axisymmetric model of the circular glass disk supported on its edge by the metal mold

In the present modeling study, it is assumed that: A1: At the bottom of the mold, there is no displacement in the z-direction: 𝑢 𝑧 (𝑟̅ , 𝑡) = 0, (𝑟̅ , 𝑡) ∈ 𝜕𝐷 𝑡 𝑚 , 𝑧 = -𝐸.

A2:

The mold is considered as a thermoelastic body that dilates when heated. A3: At time 𝑡 = 0 𝑠, when tempering begins, temperatures in the glass and in the mold are homogeneous. A4: For time 𝑡 > 0 𝑠, the glass and the mold are cooled by air and the cooling is uniform throughout the glass disk and mold surfaces. Even the bottom of the mold (𝑧 = -𝐸, 𝑅 1 ≤ 𝑟 ≤ 𝑅 2 ) undergoes cooling. A5: Mechanically, sliding contact is considered at the interface between the glass and the mold.

Thermally, in the contact zone, in addition to accounting for radiative transfer, heat exchange by conduction is considered and modeled by a constant heat transfer. A6: Gravity is not considered as only the residual stresses field is studied and not on the final shape. In fact, the bending stresses due to gravity are very small and have a negligible impact on the residual stresses.

Mathematical Model for the mechanical behavior

Since the problem is axisymmetric, a cylindrical coordinate system (𝑟, 𝜃, 𝑧) is used. The problem does not depend on 𝜃 and there is no tangential displacement. The displacement vector in the glass and in the mold is:

𝑢 ̅(𝑟̅ , 𝑡) = { 𝑢 𝑟 (𝑟̅ , 𝑡) 0 𝑢 𝑧 (𝑟̅ , 𝑡) }. (1)
In the tempering operation, displacements are very small and the linearized strain tensor can be used.

In absence of gravity and inertial forces, mechanical equilibrium is:

z r R2 R1 R e E ∇ 𝑟̅ ⋅ 𝝈(𝑟̅ , 𝑡) = 0 ̅ , (𝑟̅ , 𝑡) ∈ 𝐷 𝑡 𝑔 , (3)
where 𝜎 is the Cauchy stress tensor and ∇ 𝑟̅ the divergence operator in cylindrical coordinates.

Time is not explicitly present in equation ( 4), but will come into play through the temperature dependence of the material properties and thermal dilatation.

The boundary conditions are the following:  Due to the axisymmetry, the radial displacement vanishes on axis 𝑟 = 0: 𝑢 𝑟 (𝑟̅ , 𝑡) = 0, (𝑟̅ , 𝑡) ∈ 𝜕𝐷 𝑡 𝑔 , 𝑟 = 0.

(5)  Due to the presence of the mold, there is a unilateral contact condition on boundary 𝜕𝐷 𝑡 𝑔𝑚 . The displacements of both bodies must satisfy the Signorini condition stating that the bodies cannot interpenetrate and that a contact force only exists when the distance between both bodies vanishes [START_REF] Wriggers | Computational Contact Mechanics[END_REF].  The rest of the boundary is free of forces and the natural boundary condition:

𝝈 ⋅ 𝑛 ̅ 𝑔 = 0, (𝑟̅ , 𝑡) ∈ 𝜕𝐷 𝑡 𝑔 \𝜕𝐷 𝑡 𝑔𝑚 , (6)
where 𝑛 ̅ 𝑔 is the outer unit normal vector to boundary 𝜕𝐷 𝑡 𝑔 \𝜕𝐷 𝑡 𝑔𝑚 .

Stress tensor 𝝈(𝑟̅ , 𝑡) is linked to deformation tensor 𝜺(𝑟̅ , 𝑡) (and consequently to the displacements) through the constitutive law of glass.

Applying the constitutive law of glass in a cooling process is very complex. At a given temperature, glass is viscoelastic [START_REF] Tool | Concerning the Annealing and Characteristics of Glass[END_REF]. The stress and strain tensors are decomposed into a deviatoric tensor and a hydrostatic tensor using following relations:

𝝈(𝑟̅ , 𝑡) = 𝒔(𝑟̅ , 𝑡) + 𝑡𝑟𝑎𝑐𝑒(𝝈(𝑟̅ , 𝑡)) 3 𝑰 = 𝒔(𝑟̅ , 𝑡) + 𝜎 ℎ (𝑟̅ , 𝑡) 3 𝑰, (7)
𝜺(𝑟̅ , 𝑡) = 𝒆(𝑟̅ , 𝑡) + 𝑡𝑟𝑎𝑐𝑒(𝜺(𝑟̅ , 𝑡))

3 𝑰 = 𝒆(𝑟̅ , 𝑡) + 𝜀 ℎ (𝑟̅ , 𝑡) 3 𝑰, (8)
where 𝒆(𝑟̅ , 𝑡) is the deviatoric strain tensor, 𝒔(𝑟̅ , 𝑡) is the deviatoric stress tensor, 𝑰 is the unit tensor, 𝜀 ℎ (𝑟̅ , 𝑡) is the first strain tensor invariant, and 𝜎 ℎ (𝑟̅ , 𝑡) is the stress tensor invariant. The constitutive law then states:

𝒔(𝑟̅ , 𝑡) = ∫ 𝐺(𝜉 -𝜉 ′ ) 𝜕𝒆 𝜕𝜉 ′ (𝜉 ′ )𝑑𝜉 𝜉(𝑟̅ ,𝑡) 0 (9)
and 𝜎 ℎ (𝑟̅ , 𝑡) = ∫ 𝐾(𝜉(𝑟̅ , 𝑡) -𝜉 ′ ) 𝜕 𝜕𝜉 ′ [𝜀 ℎ (𝜉 ′ ) -𝜀 𝑡ℎ (𝜉 ′ )]𝑑𝜉 ′ , 𝜉(𝑟̅ ,𝑡) 0 [START_REF] Chen | Numerical Simulation and Experimental Stud of Residual Stresses in Compression Molding of Precision Glass Optical Components[END_REF] where 𝐺(𝑡) and 𝐾(𝑡) are respectively the shear and bulk modulus and 𝜀 𝑡ℎ the thermal strain. Variable 𝜉 is the "reduced time" used to take the temperature dependence into account through the thermorheological simplicity assumption. It is determined through the following equation:

𝜉(𝑟̅ , 𝑡) = ∫ 𝜙[𝑇(𝑟̅ , 𝑡′)]𝑑𝑡′ 𝑡 0 , ( 11
)
where 𝜙 is the "shift function" [START_REF] Schwarzl | Time-Temperature Dependence of Linear Viscoelastic Behavior[END_REF] that will be defined later. The behavior of glass during cooling is very complex since structural relaxation needs to be considered, generally by applying the fictive temperature concept [START_REF] Tool | Relation between Inelastic Deformability and Thermal Expansion of Glass in Its Annealing Range[END_REF][START_REF] Gardon | Stress and Volume Relaxation in Annealing Flat Glass[END_REF]. Roughly speaking, fictive temperature 𝑇 𝑓 represents the deviation of the structure of glass from its equilibrium state. Fictive temperature is determined through the following equation:

𝑇 𝑓 (𝑟̅ , 𝑡) = 𝑇(𝑟̅ , 𝑡) -∫ 𝑀(𝜉 -𝜉 ′ ) 𝑑𝑇 𝑑𝜉 ′ 𝑑𝜉 ′ 𝜉(𝑟̅ ,𝑡) 0 , (12)
where 𝑀(𝑡) is the fictive temperature's relaxation modulus, which depends only on the glass being considered. The shift function [START_REF] Gardon | Stress and Volume Relaxation in Annealing Flat Glass[END_REF] is defined by:

𝜙(𝑟̅ , 𝑡) = 𝑒𝑥 𝑝 (- 𝐻 𝑅 𝑔 [ 𝑥 𝑇(𝑟̅ , 𝑡) + 1 -𝑥 𝑇 𝑓 (𝑟̅ , 𝑡) - 1 𝑇 𝑟 ]), ( 13
)
where 𝐻 is an activation energy, 𝑅 𝑔 = 8.314 𝐽/(𝑚𝑜𝑙 𝐾) the universal gas constant, 𝑥 a material parameter and 𝑇 𝑟 a reference temperature at which 𝐺(𝑡) is measured. Moreover, the fictive temperature 𝑇 𝑓 is also used when determining thermal strain 𝜀 𝑡ℎ (r̅ , 𝑡) = 𝛽 𝑙 (𝑇 𝑓 (r̅ , 𝑡) -𝑇 0 (r̅ , 𝑡)) + 𝛽 𝑔 (𝑇(r̅ , 𝑡) -𝑇 𝑓 (r̅ , 𝑡)),

where 𝛽 𝑙 and 𝛽 g are the dilatation coefficients in the liquid and solid state respectively.

Remark: Material properties 𝐺(𝑡), 𝐾(𝑡) and 𝑀(𝑡) are usually provided as prony series, which is convenient for numerical computations. In the present work, they will be of the form1

𝐺(𝑡) = 𝐺 ∑ 𝐴 𝑖 𝑒 - 𝑡 𝜏 𝑖 𝑛 𝑖=1 , (15a)
𝐾(𝑡) = 𝐾, (15b)
𝑀(𝑡) = ∑ 𝑤 𝑖 𝑒 - 𝑡 𝜏 𝑠𝑖 𝑛 𝑠 𝑖=1 (15c)
This particular behavior needed to be implemented in Abaqus ® using specific subroutines for the shift function [START_REF] Mann | Determination of specific heat and true thermal conductivity of glass from dynamic temperature data[END_REF], for the computation of the fictive temperature [START_REF] Fotheringham | Active thermal conductivity of hot glass[END_REF] and for thermal dilatation [START_REF] Dusserre | Thermo-mechanical stresses in cast steel dies during glass pressing process[END_REF]. The thermo-viscoelastic glass properties used for the tempering process are listed in section 3.

The Radiative Heat Transfer Model

During glass tempering, thermal convection and radiation (in addition to thermal conduction) are important thermal processes. It is well known [START_REF] Fotheringham | Active thermal conductivity of hot glass[END_REF] that for a thick glass disk, radiation is the most important one. In the following, the importance of thermal radiation for tempering thin circular glass disks will be studied. The rotational symmetric heat transfer in the circular glass disk (Figure 1) is calculated using the heat transfer equation in cylindrical coordinates:

𝑐 𝑝 𝑔 𝜌 𝑔 𝜕𝑇 𝑔 𝜕𝑡 (𝑟̅ , 𝑡) = ∇ ̅ 𝑟̅ • (𝑘 ℎ 𝑔 ∇ ̅ 𝑟̅ 𝑇 𝑔 (𝑟̅ , 𝑡)) -∇ ̅ 𝑟̅ • 𝑞 ̅ 𝑟𝑎𝑑 (𝑟̅ , 𝑇 𝑔 ) , (𝑟̅ , 𝑡) ∈ 𝐷 𝑡 𝑔 , (16)
𝑇 𝑔 (𝑟̅ , 0) = 𝑇 0 𝑔 (𝑟̅ ), 𝑟̅ ∈ 𝐷 𝑔 , (17)
where the quantities with the subscript "𝑔" denote quantities related to the glass. Here 𝑇 𝑔 denotes the temperature in the glass disk depending on position 𝑟̅ = (𝑟, 𝑧) and time 𝑡, 𝑐 𝑝 𝑔 the specific heat capacity, 𝑘 ℎ 𝑔 the heat conductivity and 𝑇 0 𝑔 (𝑟̅ ) the initial temperature of the glass. Furthermore, the following notation is used:

∇ ̅ 𝑟̅ = ( 𝜕 𝜕𝑟 ⁄ , 𝜕 𝜕𝑧 ⁄ ) 𝑇 , ∇ ̅ 𝑟̅ • 𝐴 ̅ = 1 𝑟 𝜕 𝜕𝑟 (𝑟𝐴 𝑟 ) + 𝜕 𝜕𝑧 𝐴 𝑧 , 𝐴 ̅ = (𝐴 𝑟 , 𝐴 𝑧 ) 𝑇 .
The last term of the right hand side in [START_REF] Béchet | Two-dimensional finite element modeling of glass forming and tempering processes, including radiative effects[END_REF] represents the divergence of the radiative flux 𝑞 ̅ 𝑟𝑎𝑑 in the semitransparent wavelength region [START_REF] Modest | Radiative heat transfer[END_REF][START_REF] Siegel | Thermal radiation heat transfer[END_REF][START_REF] Viskanta | Heat transfer in semitransparent solids[END_REF] depending on the glass temperature 𝑇 𝑔 . At the boundary, heat flux is calculated using:

-𝑘 ℎ 𝑔 𝜕𝑇 𝑔 𝜕𝑛 𝑔 (𝑟̅ , 𝑡) = 𝛼 𝑔 (𝑇 𝑔 (𝑟̅ , 𝑡) -𝑇 ∞ (𝑡)) + 𝑛 ̅ 𝑔 • 𝑞 ̅ 𝑟𝑎𝑑 𝑔 , (𝑟̅ , 𝑡)𝜖𝜕𝐷 𝑡 𝑔 ∖ 𝜕𝐷 𝑡 𝑔𝑚 , (18)
-𝑘 ℎ 𝑔 𝜕𝑇 𝑔 𝜕𝑛 𝑔 (𝑟̅ , 𝑡) = 0, 𝑟 = 0, (19)
where 𝑛 ̅ 𝑔 is the outer normal vector of the glass disk. Equation ( 18) means that heat flux is composed on two terms. The first one represents convection with the ambient air at a temperature of 𝑇 ∞ (𝑡), and 𝛼 𝑔 is the convection coefficient between the glass and the surroundings. The second term represents the difference between the glass radiation and the irradiation of the surroundings at the glass surface in the opaque wavelength region. Equation ( 19) is a symmetry condition with no flux in the middle of the glass disk for 𝑟 = 0.

In addition to the heat transfer inside the glass disk, heat transfers in the mold and the interaction between these two domains must be considered.

For the metal mold, the following equations must be used:

𝑐 𝑝 𝑚 𝜌 𝑚 𝜕𝑇 𝑚 𝜕𝑡 (𝑟̅ , 𝑡) = ∇ 𝑟̅ • (𝑘 ℎ 𝑚 ∇ 𝑟̅ 𝑇 𝑚 (𝑟̅ , 𝑡)), (𝑟̅ , 𝑡) ∈ 𝐷 𝑡 𝑚 , ( 20
) 𝑇 𝑚 (𝑟̅ , 0) = 𝑇 0 𝑚 (𝑟̅ ), 𝑟̅ ∈ 𝐷 𝑚 , (21)
-𝑘 ℎ 𝑚 𝜕𝑇 𝑚 𝜕𝑛 𝑚 (𝑟̅ , 𝑡) = 𝛼 𝑚 (𝑇 𝑚 (𝑟̅ , 𝑡) -𝑇 ∞ (𝑡)) + 𝑛 ̅ 𝑚 • 𝑞 ̅ 𝑟𝑎𝑑 𝑚 , (𝑟̅ , 𝑡)𝜖𝜕𝐷 𝑡 𝑚 ∖ 𝜕𝐷 𝑡 𝑔𝑚 , (22)
-𝑘 ℎ 𝑚 𝜕𝑇 𝑚 𝜕𝑛 𝑚 (𝑟̅ , 𝑡) = 0, 𝑟 = 0, (23)
where the quantities with the subscript "𝑚" are related to the mold. 𝑇 𝑚 is the mold temperature and 𝑛 ̅ 𝑚 the outer normal vector of the mold.

As can be seen from ( 18) and ( 22), the interface condition of the contact area must be defined. It is assumed that thermal flux is proportional to the temperature difference between the glass and the mold:

-𝑘 ℎ 𝑔 𝜕𝑇 𝑔 𝜕𝑛 𝑔 (𝑟̅ , 𝑡) = 𝛼 𝑔𝑚 (𝑇 𝑔 (𝑟̅ , 𝑡) -𝑇 𝑚 (𝑟̅ , 𝑡)) + 𝑛 ̅ 𝑔 • 𝑞 ̅ 𝑟𝑎𝑑 𝑔𝑚 , (𝑟̅ , 𝑡)𝜖𝜕𝐷 𝑡 𝑔𝑚 , (24)
-𝑘 ℎ 𝑚 𝜕𝑇 𝑚 𝜕𝑛 𝑚 (𝑟̅ , 𝑡) = 𝛼 𝑔𝑚 (𝑇 𝑚 (𝑟̅ , 𝑡) -𝑇 𝑔 (𝑟̅ , 𝑡)) + 𝑛 ̅ 𝑚 • 𝑞 ̅ 𝑟𝑎𝑑 𝑚𝑔 , (𝑟̅ , 𝑡)𝜖𝜕𝐷 𝑡 𝑔𝑚 . (25)
Once more, it is important to mention that, in addition to thermal contact exchange, with coefficient 𝛼 𝑔𝑚 , it is necessary to consider radiative flux as well. To determine radiative flux, the radiative transfer equation (RTE) must be considered. Thermal radiation is a three-dimensional process. In the three-dimensional Cartesian coordinate system, rays travel on straight lines. This is why for radiative transfer, the three-dimensional Cartesian coordinate system is considered. The glass disk is defined as follows:

𝑅 𝑔 = {(𝑥, 𝑦): 𝑥 2 + 𝑦 2 ≤ 𝑅 2 , 0 ≤ 𝑧 ≤ 𝑒}, 𝑅 𝑡 𝑔 = 𝑅 𝑔 × {0 ≤ 𝑡 ≤ 𝑡 𝑚𝑎𝑥 }.
The radiative transfer inside the glass disk is described by the radiative transfer equation (RTE):

Ω ̅ • ∇𝐼(𝑥̅ , Ω ̅ , 𝜆, 𝑡) + 𝜅(𝜆)𝐼(𝑥̅ , Ω ̅ , 𝜆, 𝑡) = 𝜅(𝜆)𝐵(𝑇 𝑔 (𝑥̅ , 𝑡), 𝜆), (𝑥̅ , 𝑡) ∈ 𝑅 𝑡 𝑔 , (26)
where radiative intensity 𝐼(𝑥̅ , Ω, ̅ 𝜆, 𝑡) depends on position 𝑥̅ = (𝑥, 𝑦, 𝑧), direction Ω ̅ , wavelength λ, and time 𝑡 [START_REF] Modest | Radiative heat transfer[END_REF]. 𝜅(𝜆) is the absorption coefficient depending on wavelength λ. Planck's function, 𝐵(𝑇 𝑔 , 𝜆), is given by: Furthermore, a band model for the absorption coefficient [START_REF] Farina | Mathematical models in the manufacturing of glass[END_REF][START_REF] Lentes | Three-dimensional radiative heat transfer in glass cooling processes[END_REF] with 𝑀 𝑘 bands is assumed:

𝐵(𝑇 𝑔 , 𝜆) = 2ℎ𝑐 0 2 𝑛 𝑔 2 𝜆
𝜅(𝜆) = 𝜅 𝑘 = 𝑐𝑜𝑛𝑠𝑡., 𝜆 𝑘-1 ≤ 𝜆 < 𝜆 𝑘 , 𝑘 = 1,2, … , 𝑀 𝑘 . Using the following notation 𝐼 𝑘 (•) = ∫ 𝜆 𝑘-1 𝜆 𝑘 𝐼(•, 𝜆)𝑑𝜆, 𝐵 𝑘 (•) = ∫ 𝜆 𝑘-1 𝜆 𝑘 𝐵(•, 𝜆)𝑑𝜆,
instead of ( 26), the RTE is obtained as follows:

Ω ̅ • ∇𝐼 𝑘 (𝑥̅ , Ω ̅ , 𝑡) + 𝜅 𝑘 𝐼 𝑘 (𝑥̅ , Ω ̅ , 𝑡) = 𝜅 𝑘 𝐵 𝑘 (𝑇 𝑔 (𝑥̅ , 𝑡)), (𝑥̅ , 𝑡) ∈ 𝑅 𝑡 𝑔 , 𝑘 = 1,2, … , 𝑀 𝑘 . (27)
Now, it is necessary to add the boundary conditions. Emission and specular reflection are assumed. The boundary condition is given in the usual form:

for 𝑥̅ ∈ 𝜕𝐷 𝑡 𝑔 ∖ 𝜕𝐷 𝑡 𝑔𝑚 : 𝐼 𝑘 (𝑥̅ , Ω ̅ ) = 𝜌(Ω ̅ )𝐼 𝑘 (𝑥̅ , Ω ̅ ′) + (1 -𝜌(Ω ̅ ))𝐵 𝑔 𝑘 (𝑇 ∞ ), 𝑛 ̅ 𝑔 • Ω ̅ < 0, (28)
for 𝑥̅ ∈ 𝜕𝐷 𝑡 𝑔𝑚 :

𝐼 𝑘 (𝑥̅ , Ω ̅ ) = 𝜌(Ω ̅ )𝐼 𝑘 (𝑥̅ , Ω ̅ ′) + (1 -𝜌(Ω ̅ ))𝐵 𝑔 𝑘 (𝑇 𝑟 ), 𝑛 ̅ 𝑔 • Ω ̅ < 0, (29)
where Ω′ ̅ denotes reflected direction

Ω′ ̅ = Ω ̅ -2(Ω ̅ • 𝑛 ̅ 𝑔 )𝑛 ̅ 𝑔 .
The 𝑀 𝑘 systems of partial differential equations [START_REF] Kurkjian | Relaxation of Torsional Stress in the Transformation Range of a Soda-Lime-Silica Glass[END_REF][START_REF] Gardon | Calculation of Temperature Distributions in Glass Plates Undergoing Heat-Treatment[END_REF][START_REF] Micheletti | Space-Time Adaptation for Purely Diffusive Problems in an Anisotropic Framework[END_REF] are high dimensional and can only be solved numerically. Due to the rotational symmetry of the present problem, it is appropriate to consider only points in plane 𝑅 𝑦 𝑔 = { 0 ≤ 𝑥 ≤ 𝑅, 𝑦 = 0, 0 ≤ 𝑧 ≤ 𝑒}.

Since the radiative intensity in domain 𝑅 𝑦 𝑔 is known, it is easy to calculate all previously discussed radiative fluxes. The divergence of the radiative flux for the band model is defined as:

∇ ̅ 𝑟̅ • 𝑞 ̅ 𝑟𝑎𝑑 (𝑟̅ , 𝑇 𝑔 ) = ∇ 𝑥̅ • 𝑞 ̅ 𝑟𝑎𝑑 (𝑥̅ , 𝑇 𝑔 ) = ∑ 𝜅 𝑘 (4𝜋𝐵 𝑘 (𝑇 𝑔 ) -∫ 𝐼 𝑘 (𝑥̅ , Ω ̅ )𝑑Ω 𝑆 2 ) 𝑀 𝑘 𝑘=1 , (30)
where 𝑆 2 denotes the unit sphere. Equation (30) represents the flux in the semitransparent wavelength region.

For the surface flow, the radiative part leaving the glass disk at 𝜕𝐷 𝑡 𝑔 and the irradiation from the surroundings are considered:

𝑛 ̅ 𝑔 • 𝑞 ̅ 𝑟𝑎𝑑 𝑔 = 𝛾 𝑔 ∫ (𝐵 𝑎 (𝑇 𝑔 , 𝜆) -𝐵 𝑎 (𝑇 ∞ , 𝜆))𝑑𝜆 𝑜𝑝𝑎𝑞𝑢𝑒 , (31)
where 𝐵 𝑎 (. , 𝜆) denotes Planck's function for the wavelength in a vacuum (≈ air). Factor 𝛾 𝑔 is a measure of the "current emissivity" of the glass surface in the opaque wavelength rangethe hemispherical emissivitydefined as [START_REF] Krause | Mathematical simulation in glass technology, Schott series on glass and glass ceramics[END_REF]]

𝛾 𝑔 = ∫ (1 -𝜌(Ω ̅ ))(𝑛 ̅ 𝑔 • Ω ̅ ) 𝑛 ̅ 𝑔 •Ω ̅ >0 𝑑Ω.
Since the metal mold is an opaque material, the surface flow of the mold is obtained by:

𝑛 ̅ 𝑚 • 𝑞 ̅ 𝑟𝑎𝑑 𝑚 = 𝜎[(𝑇 𝑚 ) 4 -(𝑇 ∞ ) 4 ], (32)
where 𝜎 is Stefan-Boltzmann's constant, calculated as 𝜎 = 5.67 • 10 -8 𝑊/(𝑚 2 𝐾 4 ).

Finally, the interface condition at the glass-mold contact surface 𝜕𝐷 𝑡 𝑔𝑚 must be investigated. Figure 2 illustrates the heat exchange by radiation between the glass disk and the mold [START_REF] Krause | Mathematical simulation in glass technology, Schott series on glass and glass ceramics[END_REF]. For glass, only the surface emission in the opaque wavelength region needs to be considered:

𝑛 ̅ 𝑔 • 𝑞 ̅ 𝑟𝑎𝑑 𝑔𝑚 = 𝛾 𝑔 ∫ (𝐵 𝑎 (𝑇 𝑔 , 𝜆) -𝐵 𝑎 (𝑇 𝑚 , 𝜆)) 𝑜𝑝𝑎𝑞𝑢𝑒 𝑑𝜆 . (33)
However, this is not the case for the metal mold. Besides surface emission, radiation from the inside of the glass must be considered. A part of the radiation in the semitransparent wavelength region, On the other hand, part of

∫ (1 -𝜌(Ω ̅ ))(𝑛 ̅ 𝑔 • Ω ̅ ) ∑ 𝐼 𝑘 (𝑥̅ , Ω ̅ )𝑑Ω , 𝑀 𝑘 𝑘=1 𝑛 ̅ 𝑔 •Ω ̅ >0 (
𝛾 𝑚 ∑ 𝐵 𝑎 𝑘 (𝑇 𝑚 ) 𝑀 𝑘 𝑘=1 (35)
is emitted by the mold and absorbed by the glass disk. 𝛾 𝑚 is the hemispherical emissivity of the mold

𝛾 𝑚 = ∫ (1 -𝜌(Ω ̅ ))𝑛 ̅ 𝑚 • Ω ̅ 𝑛 ̅ 𝑚 •Ω ̅ >0
𝑑Ω.

The difference between (34) and ( 35) is the net flow absorbed by the mold in the semitransparent wavelength region. Together with the net flow in the opaque wavelength region, the following flow term is obtained for the mold:

𝑛 ̅ 𝑚 • 𝑞 ̅ 𝑟𝑎𝑑 𝑚𝑔 = 𝛾 𝑚 ∫ (𝐵 𝑎 (𝑇 𝑚 , 𝜆) -𝐵 𝑎 (𝑇 𝑔 , 𝜆))𝑑𝜆 𝑜𝑝𝑎𝑞𝑢𝑒 -∫ (1 -𝜌(Ω ̅ ))𝑛 ̅ 𝑔 • Ω ̅ ∑ 𝐼 𝑘 (𝑥̅ , Ω ̅ ) 𝑀 𝑘 𝑘=1 𝑑Ω + 𝛾 𝑚 𝑛 ̅ 𝑔 •Ω ̅ >0 ∑ 𝐵 𝑎 𝑘 (𝑇 𝑚 ) 𝑀 𝑘 𝑘=1 . (36)

Numerical Solution of the Glass Tempering Model

The numerical solution of the glass tempering model consists of two parts. The heat transfer and the mechanical calculations can be done using the Abaqus ® commercial software package. For thermal radiation, proprietary software code was developed. In the following paragraph, the numerical solution of the radiation model is described and the relationship to Abaqus ® will be discussed.

The numerical solution of the radiative transfer equation ( 27) is a challenging problem due to the non-linearity and high-dimensionality. The so-called Formal Solution Approximation (FSA) presented in [START_REF] Siedow | Application of a new method for radiative heat transfer to flat glass tempering[END_REF] is used. Since specular reflecting boundary conditions [START_REF] Gardon | Calculation of Temperature Distributions in Glass Plates Undergoing Heat-Treatment[END_REF][START_REF] Micheletti | Space-Time Adaptation for Purely Diffusive Problems in an Anisotropic Framework[END_REF] are present, the FSA method must be modified. For a given point inside the glass plate, the radiative intensity is calculated by tracing all rays back to a faraway point which has no influence on the calculated intensity.

The formal solution of the radiative transfer equation is provided through the following equation:

𝐼 𝑘 (𝑥̅ , Ω ̅ ) = 𝐼 𝑘 (𝑥̅ 𝑏 0 , Ω ̅ )𝑒 -𝜅 𝑘 𝑑(𝑥̅ ,Ω ̅ ) + 𝜅 𝑘 ∫ 𝐵 𝑘 (𝑇 𝑔 (𝑥̅ -𝑠Ω ̅ , 𝑡))𝑒 -𝜅 𝑘 𝑠 𝑑(𝑥̅ ,Ω ̅ ) 0 𝑑𝑠.
Here, 𝑑(𝑥̅ , Ω ̅ ) denotes the distance from point 𝑥̅ to boundary point 𝑥̅ 𝑏 0 in the (-Ω ̅ ) direction (see Figure 3).

𝐼 𝑘 (𝑥̅ , Ω ̅ ) ≈ 𝐼 𝑘 (𝑥̅ 𝑏 0 , Ω ̅ )𝑒 -𝜅 𝑘 𝑑 0 + (1 -𝑒 -𝜅 𝑘 𝑑 0 )𝐵 𝑘 (𝑥̅ ) - 1 𝜅 𝑘 [1 -(1 + 𝜅 𝑘 𝑑 0 )𝑒 -𝜅 𝑘 𝑑 0 ] 𝑑𝐵 𝑘 𝑑𝑇 (𝑥̅ )Ω ̅ • ∇𝑇 𝑔 (𝑥̅ , 𝑡), (37)
where 𝑑 0 = 𝑑(𝑥̅ , Ω ̅ ) and 𝐵 𝑘 (𝑥̅ ) = 𝐵 𝑘 (𝑇 𝑔 (𝑥̅ , 𝑡)). 𝐼 𝑘 (𝑥̅ 𝑏 0 , Ω ̅ ) denotes the spectral intensity in direction Ω ̅ at boundary point 𝑥̅ 𝑏 0 . For this point, boundary condition [START_REF] Gardon | Calculation of Temperature Distributions in Glass Plates Undergoing Heat-Treatment[END_REF][START_REF] Micheletti | Space-Time Adaptation for Purely Diffusive Problems in an Anisotropic Framework[END_REF] can be applied: 𝐼 𝑘 (𝑥̅ 𝑏 0 , Ω ̅ ) = 𝜌 0 𝐼 𝑘 (𝑥̅ 𝑏 0 , Ω ̅ 0 ) + (1 -𝜌 0 )𝐵 𝑔 𝑘 (𝑇 𝑏 0 𝑎 ), where 𝑇 𝑏 0 𝑎 = 𝑇 𝑏 0 ∞ for 𝑥̅ 𝑏 0 ∈ 𝜕𝐷 𝑔 and 𝑇 𝑏 0 𝑎 = 𝑇 𝑏 0 𝑚 for 𝑥̅ 𝑏 0 ∈ 𝜕𝐷 𝑔𝑚 . Ω ̅ 0 is the new direction according to Snell's law Ω ̅ 0 = Ω ̅ -2(Ω ̅ , 𝑛 ̅ 𝑔 )𝑛 ̅ 𝑔 . 𝑛 ̅ 𝑔 denotes the outer normal vector at glass boundary point 𝑥̅ 𝑏 0 .

Point 𝑥̅ 𝑏 0 is an inner point with respect to new direction Ω ̅ 0 . Thus, formula (37) can be used once again to approximate 𝐼 𝑘 (𝑥̅ 𝑏 0 , Ω ̅ 0 ).

This process can continue to be used. The following approximation for spectral intensity 𝐼 𝑘 (Ω ̅ ) in point 𝑥̅ and for direction Ω ̅ is then obtained:

𝐼 𝑘 = 𝐼 𝑘 (𝑥̅ , Ω ̅ ) ≈ (1 -𝑒 -𝜅 𝑘 𝑑 0 )𝐵 𝑘 (𝑥̅ ) + ∑ 𝐺 𝑙 (1 -𝑒 -𝜅 𝑘 𝑑 𝑙+1 ) 𝑛𝑏 𝑙=0 𝐵 𝑘 (𝑥̅ 𝑏 𝑙 ) + ∑ 𝐹 𝑙 (1 -𝜌 𝑙 )𝐵 𝑘 (𝑇 𝑏 𝑙 𝑎 ) 𝑛𝑏 𝑙=0 - - 1 𝜅 𝑘 [1 -(1 + 𝜅 𝑘 𝑑 0 )𝑒 -𝜅 𝑘 𝑑 0 ] 𝑑𝐵 𝑘 𝑑𝑇 (𝑥̅ )Ω ̅ • ∇𝑇 𝑔 (𝑥̅ ) - (38) - 1 𝜅 𝑘 ∑ 𝐺 𝑙 [1 -(1 + 𝜅 𝑘 𝑑 𝑙+1 )𝑒 -𝜅 𝑘 𝑑 𝑙+1 ] 𝑑𝐵 𝑘 𝑑𝑇 (𝑥̅ 𝑏 𝑙 )Ω ̅ 𝑙 • ∇𝑇 𝑔 (𝑥̅ 𝑏 𝑙 ). 𝑛𝑏 𝑙=0
The following denotations were used: 𝑑 0 = 𝑑(𝑥̅ , Ω ̅ ), 𝑑 𝑙+1 = 𝑑(𝑥̅ 𝑏 𝑙 , Ω ̅ 𝑙 ), 𝜌 𝑙 : reflection coefficient at point 𝑥̅ 𝑏 𝑙 , 𝑛𝑏: number of backtrackings, 𝐺 𝑙 = 𝜌 0 • … • 𝜌 𝑙 𝑒 -κ 𝑘 (𝑑 0 +...+𝑑 𝑙 ) , 𝐹 0 = 𝑒 -κ 𝑘 𝑑 0 , 𝐹 𝑙 = 𝜌 𝑙-1 𝑒 -κ 𝑘 𝑑 𝑙 𝐹 𝑙-1 .

In practice, not all backtrackings will be performed. The process is stopped if 𝐺 𝑙 is smaller than a predefined small number (𝜀).

Knowing the spectral radiative intensity for all spatial points and all directions, it is easy to calculate the divergence of the radiative flux vector (30):

∇ 𝑥̅ • 𝑞 ̅ 𝑟𝑎𝑑 (𝑥̅ 𝑖 ) ≈ ∑ 𝜅 𝑘 (4𝜋𝐵 𝑘 (𝑇 𝑔 (𝑥̅ 𝑖 )) -∑ 𝐼 𝑖𝑗 𝑘 𝑤 𝑗 𝑚-1 𝑗=0 ) 𝑀 𝑘 𝑘=1
, where 𝐼 𝑖𝑗 𝑘 is the spectral intensity of the 𝑘-th band in point 𝑥̅ 𝑖 , and direction Ω ̅ 𝑗 calculated according to (38).

To discretize the direction Ω ̅ , the LSH-8 method proposed by Fiveland [START_REF] Fiveland | The selection of discrete ordinate quadrature sets for anisotropic scattering[END_REF] was used. Here, 𝑤 𝑗 denotes the weight for direction Ω ̅ 𝑗 . In a similar way, the boundary flux (36) can be calculated:

𝑛 ̅ 𝑚 • 𝑞 ̅ 𝑟𝑎𝑑 𝑚𝑔 ≈ 𝛾 𝑚 (𝐵 𝑎 𝑜𝑝 (𝑇 𝑚 ) -𝐵 𝑎 𝑜𝑝 (𝑇 𝑔 )) - ∑ (1 -𝜌(Ω ̅ 𝑗 )) 𝑛 ̅ 𝑔 • Ω ̅ 𝑗 𝑛 ̅ 𝑔 •Ω ̅ 𝑗 >0 ∑ 𝐼 𝑖𝑗 𝑘 𝑤 𝑗 𝑀 𝑘 𝑘=1 + 𝛾 𝑚 ∑ 𝐵 𝑎 𝑘 (𝑇 𝑚 ), 𝑀 𝑘 𝑘=1 (39)
where

𝐵 𝑎 𝑜𝑝 (𝑇) = ∫ 𝐵 𝑎 (𝑇, 𝜆)𝑑𝜆.

𝑜𝑝𝑎𝑞𝑢𝑒

The calculation of the other boundary fluxes is straightforward.

The Abaqus® software package was used to solve the mechanical and thermal equations to simulate glass disk tempering. The thermal calculation is first performed at one time step ∆𝑡 between time 𝑡 and time 𝑡 + ∆𝑡. The mechanical calculation is performed after. For the thermal calculation, Abaqus ® provides information on glass geometry and temperatures to the radiation model via specific subroutines developed for that purpose. Then, radiation terms are computed as previously described and resubmitted to Abaqus® to complete the thermal computations. Finally, at each time step, mechanical computations are performed to determine the stresses in the glass. The whole process ends when tempering is complete.

Validation of the Radiative Heat Transfer Model

The radiative heat transfer model (Section 2.3) and its numerical solution (Section 2.4) were validated by studying only the natural cooling of a glass disk. This means that no support mold was considered, and consequently, boundary 𝜕𝐷 𝑡 𝑔𝑚 does not exist. Equations [START_REF] Tool | Concerning the Annealing and Characteristics of Glass[END_REF][START_REF] Schwarzl | Time-Temperature Dependence of Linear Viscoelastic Behavior[END_REF][START_REF] Tool | Relation between Inelastic Deformability and Thermal Expansion of Glass in Its Annealing Range[END_REF][START_REF] Gardon | Stress and Volume Relaxation in Annealing Flat Glass[END_REF] pertaining to the mold and the boundary conditions [START_REF] Lentes | Three-dimensional radiative heat transfer in glass cooling processes[END_REF][START_REF] Fiveland | The selection of discrete ordinate quadrature sets for anisotropic scattering[END_REF] pertaining to the thermal exchanges between the mold and the glass were not considered in the present section.

This example is pertinent to validate the radiative heat transfer model because, since there was natural cooling, there were fewer convection effects and the cooling mainly occurred through radiation. In 1990, Field and Viskanta developed a specific experiment for the cooling of a square plate (width 0.197 𝑚, thickness 11.68 • 10 -3 𝑚) [START_REF] Field | Measurement and Prediction of the Dynamic Temperature Distributions in Soda-Lime Glass Plates[END_REF]. They succeeded in measuring the temperature changes in the plate, in the half thickness and on the surface for the first 450𝑠 of cooling. In 1998, Lee and Viskanta developed an axisymmetric simulation model for glass disk cooling by convection [START_REF] Lee | Transient conductive-radiative cooling of an optical quality glass disk[END_REF]. By using the Finite Volume Method to solve the heat transfer equation and the Discrete Ordinate Method to solve the radiative transfer equation, these researchers validated their model by comparing their solution in the disk's center to the measurements in the plate's center eight years earlier [START_REF] Field | Measurement and Prediction of the Dynamic Temperature Distributions in Soda-Lime Glass Plates[END_REF]. The input data used by Lee and Viskanta [START_REF] Lee | Transient conductive-radiative cooling of an optical quality glass disk[END_REF] were considered when validating the present model.

The computation was made using a structured mesh of 2000 DCAX8 elements in Abaqus ® (50 along the thickness, 40 along the radius). The mesh was refined near the surface and the edge because of the high temperature gradient change in these locations. Like for Lee and Viskanta [START_REF] Lee | Transient conductive-radiative cooling of an optical quality glass disk[END_REF], insulation was applied to the axisymmetry axis to consider the axisymmetry; at the circumference (𝑟 = 𝑅), the reflectivity was assumed to be unity and the boundary condition at that surface for conductive transfer was assumed to be adiabatic; on the upper and lower surfaces of the glass disk, the natural convection in the surrounding air was considered: air temperature 𝑇 ∞ (𝑡) = 293.15 𝐾 and convection coefficient between the glass and the surroundings 𝛼 𝑔 = 4.25 𝑊/ (𝑚 2 𝐾). At time 𝑡 = 0 𝑠, when the cooling began, the temperature in the glass disk was assumed to be uniform and equal to 798.15 𝐾.

The temperature changes in the disk's center obtained by the method presented in Sections 2.3 and 2.4 were compared with the numerical results obtained by Lee and Viskanta [START_REF] Lee | Transient conductive-radiative cooling of an optical quality glass disk[END_REF] and the measurements of Field and Viskanta [START_REF] Field | Measurement and Prediction of the Dynamic Temperature Distributions in Soda-Lime Glass Plates[END_REF]. In Figure 4(a), the comparison was made at the half-thickness and in Figure 4(b), on the surface. Using the radiative model (Section 2.3) and the numerical solution (Section 2.4), the results obtained for the core (Figure 4(a)) and surface (Figure 4(b) temperatures were consistent with the experimental data. After 200 s, the temperature difference was about 5K and decreased 3K after 400s of cooling. In comparison, the difference between experimentation and simulation for Lee and Viskanta [START_REF] Lee | Transient conductive-radiative cooling of an optical quality glass disk[END_REF] was 11K after 200 s and 15K after 400 s. Since the numerical results were closer to actual measurements than for Lee and Viskanta, the radiative heat transfer model and numerical solution presented in this paper were validated. 3

Modeling results and discussion

The finite element modeling including volume radiation presented in Section 2 was applied to the cooling of a glass disk with dimensions 𝑅 = 0.150 𝑚, 𝑒 = 0.006 𝑚 (Figure 1) and initial temperature 𝑇 0 𝑔 (𝑟̅ ) = 923 𝐾. The dimensions of the metal mold were: 𝑅 1 = 0.140 𝑚, 𝑅 2 = 0.170 𝑚 and 𝐸 = 0.005 𝑚. The initial temperature 𝑇 0 𝑚 (𝑟̅ ) of the mold was 293 𝐾. Natural cooling occurred with ambient air at temperature 𝑇 ∞ (𝑡) equal to 293 𝐾 and a convection coefficient on the glass 𝛼 𝑔 and on mold surfaces 𝛼 𝑚 equal to 15 𝑊/(𝑚 2 𝐾). The value of the heat transfer coefficient between glass and mold, 𝛼 𝑔𝑚 = 100 𝑊/(𝑚 2 𝐾), was given by AGC. Due to the derivation in 2.3 the coefficient is free of any radiative contribution. Following the discussion and laboratory experiments described in [START_REF] Krause | Mathematical simulation in glass technology, Schott series on glass and glass ceramics[END_REF] it is reasonable to take such a value.

All thermal and mechanical data needed for the simulations are shown in Table 1 -Table 6. All the thermal and elastic properties were taken as constant on this study like in [START_REF] Soules | Finite-Element Calculation of Stresses in Glass Parts Undergoing Viscous Relaxation[END_REF]. 1: Elastic properties of glass [START_REF] Soules | Finite-Element Calculation of Stresses in Glass Parts Undergoing Viscous Relaxation[END_REF] Table 2: Relaxation properties of glass [START_REF] Soules | Finite-Element Calculation of Stresses in Glass Parts Undergoing Viscous Relaxation[END_REF][START_REF] Kurkjian | Relaxation of Torsional Stress in the Transformation Range of a Soda-Lime-Silica Glass[END_REF] The thermomechanical problem presented in section 2.1 was incrementally and iteratively solved using Abaqus ® finite element software. Both mechanical and thermal equations were solved using the mesh in Figure 5(a). The meshes of the disk and the mold were made of CAX8T elements with biquadratic interpolation for displacement and bilinear interpolation for temperatures. The glass is composed of 5007 nodes and 1608 elements whereas the mold is composed of 2663 nodes and 840 elements.

Along the thickness, 25 elements were used with refinement near surfaces to get a correct estimation of the temperature gradients due to convection and surface radiation during tempering (Figure 5(b)). The mesh was also refined in the area in contact with the mold. Elements near the center of the plate are anisotropic, this is adapted to the anisotropic nature of the physical phenomenon: gradients of temperature and gradients of radiative intensity are high in the z-direction and nearly vanish in the r-direction. An example of anisotropic (and adaptive) mesh for the heat equation with convection can be seen in [START_REF] Micheletti | Space-Time Adaptation for Purely Diffusive Problems in an Anisotropic Framework[END_REF].

(a) (b) To see the differences between the model with volume and surface radiation as presented in Section 2 and a model ignoring volume radiation, two numerical simulations were performed: one simulation denoted "Vol" as described in Section 2, in which volume and surface radiation are considered (resulting in (𝑇, 𝜎) 𝑉𝑜𝑙 ) and a second simulation denoted "Surf" considering only surface radiation (resulting in (𝑇, 𝜎) 𝑆𝑢𝑟𝑓 ). The temperatures and stresses during the tempering process were compared at two locations: at the center of the glass (𝑟 = 0. 𝑚𝑚) and in the glass/mold contact area (𝑟 = 143.7 𝑚𝑚), with special attention paid to the upper and lower surfaces of the glass disk at these positions.

Temperature results

Figure 6(a) shows the temperature changes on the central upper and lower surfaces of the glass disk computed by the "Vol" model. There was almost no difference between the lower surface and the upper surface due to the presence of the same convection cooling on the surfaces and the long distance from the support mold. After 1000 s, the glass cooled to a temperature of about 310 K. At this temperature, the residual stresses were no longer changing. The temperature differences between the "Vol" and "Surf" models ∆𝑇 𝑔 = (𝑇 𝑔 ) 𝑉𝑜𝑙 -(𝑇 𝑔 ) 𝑆𝑢𝑟𝑓 are presented in Figure 6(b). For the first 50 s, the temperatures calculated by the "Vol" model were higher than those of the "Surf" model in which volume radiation was ignored. The maximum difference of +12 𝐾 was obtained after 5 s of cooling. In this first period, the "Surf" model calculated a transfer of too much energy from the glass surface to the surroundings. The thermal flux was too large due to the boundary condition 4 ) of the "Surf" model. The entire wavelength region was assumed to be opaque. In contrast, in the "Vol" model, there was a transfer of energy from the boundary to the surroundings in the opaque wavelength region and from the middle of the glass to the boundaries in the semitransparent wavelength region. After 50 s of cooling, the surface temperatures calculated with the "Vol" model were lower than those obtained by "Surf" model (Figure 6(b)). Here, the maximal temperature difference of -13 𝐾 occurred after about 255 s of cooling. The cooling of the glass surfaces at the glass/mold contact area (Figure 7(a)) was very similar to what was computed at the center (Figure 6(a)). The lower surface of the glass disk cooled a little bit faster than the upper surface due to the presence of the mold, which had been at ambient temperature when cooling started (the average difference was 57 𝐾 during the 100 s of cooling). The trends in differences in temperatures between "Vol" and "Surf" models for both surfaces (Figure 7(b)) were also very similar to those at the disk's center (Figure 6(b)). The maximal positive temperature difference between the "Vol" and "Surf" models was obtained after the same cooling duration, 5 s, and was equal to +12 𝐾. The greatest negative value -13 𝐾 was the same but appeared earlier than in the center, and at different moments for each of the surfaces, i.e., 140 s for the lower surface and 165 s for the upper one (Figure 7(b)). In the disk's center (Figure 6(a)) and in the glass/mold contact area (Figure 7(a)), during the first 50 s of cooling, the temperature on both glass surfaces remained over 700 𝐾 (i.e., over 𝑇 𝑔 ). In that period, the cooling on the glass surfaces using the "Surf" model was faster than with the "Vol" model in which volume radiation was considered (Figures 6(b) and 7(b)). It is known that the temperature distribution in the glass over and near 𝑇 𝑔 is important for the residual stress of the glass. In the following section, the computational results for transient and residual stresses are discussed.

-𝑘 ℎ 𝑔 𝜕𝑇 𝑔 𝜕𝑛 𝑔 (𝑟̅ , 𝑡)~((𝑇 𝑔 ) 4 -(𝑇 ∞ )

Stress results

Glass is annealed to reduce residual stresses at the end of the cooling. To avoid breakage during the annealing process or during the subsequent use of the annealed glass, one must avoid tensile stress on the glass surface. All the stresses discussed here are the circumferential stresses 𝜎 𝜃𝜃 .

In Figure 8, the final repartition of tensile and compressive stresses is represented. Classically, tensile stresses are obtained at the surfaces whereas tensile stresses are obtained inside the glass. The presence of the mold increases the compressive area at the edge of the plate which corresponds nearly to the contact area length. The time evolutions of the stresses at the 4 points corresponding to the 4 crosses in Figure 8 are given in Figure 9 and Figure 10. For the surface points at the disk's center, the change in transient stresses is shown in Figure 9 (a). Although there were no differences in temperature, the transient compressive stresses on the upper surface are greater than on the lower surface (the negative values indicate compressive stresses). This difference increased during the first 250 s and stabilized afterwards. After 1000 s of cooling, the residual stress level was reached and the compressive stresses were respectively -32 𝑀𝑃𝑎 for the lower surface and -37 𝑀𝑃𝑎 for the upper surface. The upper surface was in a more compressed stress state because of the influence of the support mold on the behavior of the glass, even in the center. Due to the annealing by natural cooling on the glass faces, the residual stresses were relatively small. The stress differences ∆𝜎 𝑔 = (𝜎 𝑔 ) 𝑉𝑜𝑙 -(𝜎 𝑔 ) 𝑆𝑢𝑟𝑓 in the glass calculated by "Vol" and "Surf" models werecompared with the absolute valuessignificant (Figure 9(b)). The difference for the two surfaces was over +10 𝑀𝑃𝑎 after 250 s of cooling and remained at this level until the end. The conclusion was that if volume radiation is ignored in the simulation and only surface radiation is considered, the calculated residual compressive stresses are overestimated. For the "Surf" model, after 1000 s of cooling the residual stresses were respectively -48 𝑀𝑃𝑎 on the upper surface and -42 𝑀𝑃𝑎 on the lower surface. The difference between the two radiation models is therefore about -11 𝑀𝑃𝑎, or more than 20% compared with the results with volume radiation. In the glass/mold contact area (Figure 10(a)), the time to obtain the final residual stresses was much longer than in the disk's center (Figure 9(a)). Although the temperature differences between the upper and lower surfaces were small in Figure 7(b), the differences in transient and residual stresses were very large. Whereas the residual stress of the upper surface was -63 𝑀𝑃𝑎, the residual stress of the lower surface was -136 𝑀𝑃𝑎.

This means that, due to the presence of the support mold, even if natural cooling is considered, the lower glass surface in contact with the mold has been tempered.

Regarding the differences between "Vol" and "Surf" models (Figure 10

Influence of the value of the heat transfer coefficient 𝜶 𝒈𝒎

In the previous section, the value of heat transfer coefficient 𝛼 𝑔𝑚 between the mold and the glass was equal to 100 𝑊/(𝑚 2 𝐾). For glass with uniform initial temperature of 1773 K in [START_REF] Krause | Mathematical simulation in glass technology, Schott series on glass and glass ceramics[END_REF] this coefficient 𝛼 𝑔𝑚 was calculated and varies from 300 to 1000 𝑊/(𝑚 2 𝐾). In the present study the glass is only at 923 K. In order to see the influence of this coefficient on the results (temperature and residual stresses), computations of temperature and stresses were performed with two different values of the heat transfer coefficient, respectively 𝛼 𝑔𝑚 = 500 𝑊/(𝑚 2 𝐾) and 𝛼 𝑔𝑚 = 1000 𝑊/(𝑚 2 𝐾).

The results show that changing the coefficient 𝛼 𝑔𝑚 has almost no influence on the temperature at the center of the plate. This is evident as the heat flux in this part is exclusively in the 𝑧-direction. Consequently, only the results in the glass/mold contact area are presented here (Figures 11 and12). (a) temperature changes in the glass during cooling ("Vol" model) (b) difference in the temperature changes of "Vol" and "Surf" models By increasing the value of 𝛼 𝑔𝑚 , and consequently the heat exchange between the glass and the mold, the temperature decreases faster in the contact area (Figures 11(a) and 12(a)). The temperature differences at both surfaces between "Vol" and "Surf" models are nearly constant, when 𝛼 𝑔𝑚 increases (Figures 11(b) and 12(b)). Finally, the results for residual stresses at both surfaces are compared for the 3 values of 𝛼 𝑔𝑚 and for "Vol" and "Surf" models in Table 8. In the contact area, the compressive residual stresses increase greatly when 𝛼 𝑔𝑚 increases. The glass is tempered by the mold especially at the lower surface. At the upper surface this effect is substantially lower. When 𝛼 𝑔𝑚 is multiplied by 10 (from 100 𝑊/(𝑚 2 𝐾) to 1000 𝑊/(𝑚 2 𝐾)), the residual stresses are multiplied by 2.5 at the lower surface and by 1.5 at the upper surface. At the lower surface, the difference between both model "Vol" and "Surf" is nearly constant around 5.5 MPa. At the upper surface, the difference is +6.1 MPa for 𝛼 𝑔𝑟 = 100 𝑊/(𝑚 2 𝐾) but becomes -3.1 MPa for 𝛼 𝑔𝑚 = 1000 𝑊/(𝑚 2 𝐾).

In the center of the disk (𝑟 = 0), despite identical temperature evolutions with the 3 values of 𝛼 𝑔𝑚 , the final residual stresses are modified. Compressive stresses are more and more asymmetrical when 𝛼 𝑔𝑚 increases: the residual stress at the upper surface is nearly constant (around 38 MPa) whereas the residual stress at the lower surface diminishes from 31.9 𝑀𝑃𝑎 to 20.9 𝑀𝑃𝑎.

With the "Surf" model, the differences are even more important. Residual stresses are globally much higher in the center of the plate. At the upper surface, the difference between both model reaches 33.2 𝑀𝑃𝑎.

Conclusions

To correctly model cooling, realistic physical models must be used. The radiative transfer equation was developed and solved to incorporate volume radiation into the finite element modeling of the annealing of a circular disk supported on the edge by a mold.

In comparison with considering surface radiation only, like for an opaque body or in many current glass tempering simulations, the volume radiation model developed in this study demonstrates an approximately 20% difference in the stress calculated at the center of the disk. It was also shown that the support mold changes the expected cooling behavior of the glass disk, producing late tempered stress state on the edge on the glass in the vicinity of the mold contact area.

Other computations with different conditions (higher convection or thinner glass) were also carried out. They show that the importance of radiation diminishes with higher cooling rates. When the glass is thinner, the residual stresses are globally lower but the full radiation model may be important in certain cases to predict dangerous tensile stresses at the surface. The computations with different values of heat transfer coefficient reveal that a too high exchange with the mold cause higher stresses on the whole plate which is not appropriate if annealed glass is expected. For industrial use, especially if the safety region for breakage is small, the application of the model should be as precise as possible. Since glass is a semitransparent material, the presented volume radiation model is closer to reality than a surface radiation model or a model that does not take radiation into consideration at all.

Figure 2 :

 2 Figure 2: Heat exchange by radiation between glass and mold

Figure 3 :

 3 Figure 3: Back tracing of a ray

Figure 4 :

 4 Figure 4: Comparison between numerical and experimental temperature changes in the disk's center:(a) Core temperature (b) Surface temperature (blue dashes: the method presented in Sections 2.3 and 2.4, red circles: measurements[START_REF] Field | Measurement and Prediction of the Dynamic Temperature Distributions in Soda-Lime Glass Plates[END_REF], black dots: numerical results[START_REF] Lee | Transient conductive-radiative cooling of an optical quality glass disk[END_REF]).

Figure 5 :

 5 Figure 5: Finite element mesh of the glass disk and the support mold: (a) general view (b) close-up of the contact area.

Figure 6 :

 6 Figure 6: Temperature results on the upper (dashed line) and lower surfaces (dotted line) at the center of the glass disk (𝑟 = 0. 𝑚𝑚): (a) temperature changes in the glass during cooling ("Vol" model) (b) difference between the temperature changes of the "Vol" and "Surf" models

Figure 7 :

 7 Figure 7: Temperature results on the upper (dashed line) and lower surfaces (dotted line) in the glass/mold contact area (𝑟 = 143.7 𝑚𝑚): (a) temperature changes in the glass during cooling ("Vol" model) (b) difference in the temperature changes of "Vol" and "Surf" models

Figure 8 :

 8 Figure 8: Final residual stresses with "Vol" model

Figure 9 :

 9 Figure 9: Stress results on the upper (dashed line) and lower surfaces (dotted line) in the center of the glass disk (𝑟 = 0. 𝑚𝑚): (a) changes in stress of the glass during cooling ("Vol" model) (b) difference of changes in stress between "Vol" and "Surf" models

  (b)), transient and residual stresses on the upper surface were overestimated by the "Surf" model (as observed at the center of surfaces, Figure9(b)). They are underestimated for lower surface in contact with the support mold. After 1000 s, the absolute difference in the residual stresses ∆𝜎 𝑔 = (𝜎 𝑔 ) 𝑉𝑜𝑙 -(𝜎 𝑔 ) 𝑆𝑢𝑟𝑓 was about 6 𝑀𝑃𝑎 for both surfaces. Comparing the residual stresses computed by "Vol" model, -63 𝑀𝑃𝑎 for the upper surface, -136 𝑀𝑃𝑎 for the lower surface, these stresses represent 4.4 % of the lower surface in contact with the mold and 9.5 % of the upper surface, which is important.

Figure 10 :

 10 Figure 10: Stress results on the upper (dashed line) and lower surfaces (dotted line) in the glass/mold contact area (𝑟 = 143.7 𝑚𝑚): (a) changes in glass stress during cooling ("Vol" model), (b) difference in changes in glass stress between "Vol" and "Surf" models

Figure 11 :Figure 12 :

 1112 Figure 11: Temperature results on the upper (dashed line) and lower surfaces (dotted line) in the glass/mold contact area (𝑟 = 143.7 𝑚𝑚) for 𝛼 𝑔𝑚 = 500 𝑊/(𝑚 2 𝐾): (a) temperature changes in the glass during cooling ("Vol" model) (b) difference in the temperature changes of "Vol" and "Surf" models

  values and at different positions (in the center of the glass disk and in the glass/mold area, on the upper and the lower surfaces)

  𝑘 = 1.3807 • 10 -23 𝐽/𝐾 -Boltzmann's constant, 𝑐 0 = 2.998 • 10 -8 𝑚/𝑠 -speed of light in vacuum). 𝑛 𝑔 is the refractive index of the glass.

	5 (𝑒	𝑛 𝑔 𝜆𝑘𝑇 𝑔 -1) ℎ𝑐 0	,
	(ℎ = 6.626 • 10 -34 𝐽𝑠 -Planck's constant,

Table 3 .

 3 Data for the shift function[START_REF] Soules | Finite-Element Calculation of Stresses in Glass Parts Undergoing Viscous Relaxation[END_REF]

Table 4 .

 4 Thermal properties[START_REF] Soules | Finite-Element Calculation of Stresses in Glass Parts Undergoing Viscous Relaxation[END_REF]

	Band number	𝜆 𝑘-1 in 𝜇𝑚	𝜆 𝑘 in 𝜇𝑚	𝜅 𝑘 in 𝑚 -1
	1	0.5	2.75	20
	2	2.75	4.5	330	Refractive index of glass	𝑛 𝑔 = 1.46
	3	4.5	6.0	5000	Emissivity of glass	𝛾 𝑔 = 0.92
		6.0	∞	opaque	Emissivity of the mold	𝛾 𝑚 = 0.8

Table 5 .

 5 Absorption coefficients[START_REF] Gardon | Calculation of Temperature Distributions in Glass Plates Undergoing Heat-Treatment[END_REF]

Table 6 .

 6 Other radiation properties (given by AGC)

Table 8 -

 8 Residual stresses (in MPa) for different 𝜶 𝒈𝒎

Sometimes, a long term shear modulus 𝐺 0 is added in the expression (15a). This long term modulus will be present in the sequel through a very long characteristic time 𝜏 𝑖 .

Acknowledgments

This research was supported by the International Campus on Safety and Intermodality in Transportation, the Nord/Pas-de-Calais Region, the European Community, the Regional Delegation for Research and Technology, the Ministry of Higher Education and Research and the National Center for Scientific Research. The authors gratefully acknowledge the support of these institutions.