The entire glass sheet forming process consists of heating and forming a glass sheet and cooling and tempering it afterwards. For the first step, the glass sheet is heated using a local radiative source and deforms by sagging. In the thermo-mechanical calculations, temperature dependent glass viscosity, heat exchange with the ambient air and radiative source effects should be considered. A two-dimensional finite element model with plane deformation assumptions is developed. Using the P1-Approximation, the formulation and numerical resolution of the Radiative Transfer Equation (RTE) are performed on the glass domain as it changes over time to estimate the flux of the radiative body at each position in the glass. In the next step, the sheet is cooled. Narayanaswamy’s model is used to describe the temperature dependent stress relaxation and the structural relaxation. The RTE is again solved using the P1-Approximation to consider the internal radiative effects during the cooling. There is a discussion using the P1-Approximation and comparing the results to other existing methods for the temperature changes of the glass throughout the forming process, for the deformed shape at the end of the forming step and for the residual stresses after tempering.
Glass is a semi-transparent material with highly temperature dependent mechanical and thermal behaviors. In reality as well as in modeling, deforming glass to achieve the final desired geometry requires mastering the coupled changes of the temperature inside the glass and the product shape. Modeling a glass forming process, where large deformations occur, means computing the solution of one complex and non-linear thermo-mechanical problem. On one hand, it is necessary to determine the heat conduction of the glass taking the conditions of the boundaries between the glass and the ambient air and the glass and the forming tools into account. On the other hand, the model would also need to solve the mechanical problem with the temperature dependent viscosity of the glass and the changing contact conditions imposed by the forming tools. In the last two decades, the modeling of glass forming has been widely developed [1,2] using commercial software packages or homemade codes. For specific applications, glass must be tempered to make it more resistant and safer. During the tempering phase, the deformations are limited but there is a coupling between temperature and stress relaxation.
Besides heat conduction and heat convection, radiation plays an important role and for high temperatures, thermal radiation is the dominant heat transfer process. Very complete assessments of radiative heat transfer can be found in [3–5] and its application to the glass industry can be found in [6,7].
Since the geometry of glass changes during forming, different glass forming modeling solutions were proposed to account for radiation effects. The first and simplest one involves totally ignoring the radiation effects [8]. Another solution consists of using the Stefan-Boltzmann’s law and considering only surface radiation [9], This law is normally applied to opaque bodies, which is not the case for glass. Only a certain part of the radiant energy corresponding to the opaque spectrum of the glass is directly absorbed at the surface [10]. From a numerical point of view, the surface radiation is often approximated and taken into account by modifying the convection coefficient. This provides a linear relation instead of the T4 non-linearity in temperature [9].
Glass is a semi-transparent material, and internal radiative effects also occur inside the glass. A widely used solution involves using equivalent conductivity (such as the active thermal conductivity method [11,12]), which is temperature dependent, or the Rosseland’s approximation [13]. The Rosseland’s approximation treats thermal radiation as a correction of heat conductivity, which is computed before the FEM computation using absorption coefficients. This is why it is so quick and easy to integrate into commercial software packages. It is used not only in glass forming modeling, where the domain changes overtime, but also extensively used for fixed domains with negligible deformations, such as glass tempering modeling. Originally, the method was derived in 1924 by S. Rosseland [14] to investigate stellar radiation. This is why the method is valid only for optically thick glass. Furthermore, it was shown in [15] that, for glass tempering, using the Rosseland’s approximation produces vast errors in transient stress calculations.
The exact method for taking radiation in glass into account is to solve the radiative transfer equation. From a numerical point of view, this is a challenge because of the high-dimensionality and the non-linearity of RTE. A detailed discussion about different numerical methods for solving the radiative transfer equation and many more references can be found in [3] and [7]. Due to that fact that, during glass forming, one must deal with changes in the geometry of the glass plate being formed, the P1-Approximation for numerically solving the radiative transfer equation is used for the research presented here. Using the method of moments [7], one can obtain a system of two-dimensional diffusion equations instead of the high-dimensional radiative transfer equation.
In the present glass sagging modeling case, a laser is used to create local heating of the glass, as examined for the thermal computations of the cutting [16], drilling [17] or scribing [18] processes. In these studies, radiative heating was usually performed in a simple way. A first method for taking the radiative source into account is to consider it as a surface flux [16,19,20] whose intensity is directly related to the power of the radiative source. This approximation is good if the source wavelength is in the opaque zone of the glass under consideration. If this is not the case, a more precise method must be used. Another method consists of applying Beer’s law, which involves nothing more than solving the one-dimensional radiative transfer equation [17,20]. In [21], Li et al. compare these two methods. The comparison reveals that both methods are very similar for optically thick glass. However, for optically thin glass, the Beer’s law model provides much better results with respect to experimental data.
The research discussed here focuses on the two-dimensional (2-D) modeling of the extent of gravity sagging in a glass sheet being heated with a laser and of the tempering after sagging. This thermo-mechanical problem, for which both deformations and temperatures of the glass sheet must be computed, is a complex incremental and iterative problem that can be solved using commercial software. In the case of a deformable body discussed here, the P1-Approximation is used, for the forming and the tempering steps, to solve the RTE instead of using more simplified solutions found in the literature. An initial discussion on the temperature changes in the glass sheet during forming and on the deformed shape of the glass at the end of the forming step is proposed to compare the solution obtained by the P1-Approximation for radiation to other simplified ones. The discussion continues with the temperatures during the cooling phase and the residual stresses after tempering.
In this paper, the glass was first exposed to a local heat using a laser source, and subsequently deformed through gravity sagging. In a second step, the deformed shape was tempered. The glass sheet was clamped on the left side (y = 0) and exposed to uniform convection with ambient air on all its surfaces. Considering uniform thermal conditions in the x-direction and assuming that the dimension of the sheet in the x-direction was much larger than in the(y, z)-directions, the problem can be reduced to two dimensions with generalized plane strain conditions. Under this assumption, there is no heat transfer in the (x, z) plane but dilatation effects are allowed in the x-direction. Finally, the problem was solved on the following domain:
w denotes the sheet thickness, l the length (Figure 1) and tmax the heating duration.
Figure 1: Description of 2-D glass sagging under radiative laser heating.
During the heating phase, the glass plate was surrounded by hot air and a laser was applied at point (d, w) with constant power. The glass is only formed due to the gravity. For the tempering step, the laser was switched off and cool air was blown all over the glass sheet.
The static equilibrium of the deformable glass sheet in the presence of gravitational effects without inertial effects is described using:
where σ is the Cauchy stress tensor in the deformed glass sheet, ρ the density of the glass and 
the gravitational force. Note that the derivative is taken with respect to position 
in the actual deformed glass sheet. The boundary conditions of the glass surface ∂D are affected by a non-displacement condition imposed by clamped side ∂Lu of the glass sheet and the fact that there was no external force acting on the other glass sheet boundary ∂D\∂Lu. They are described using:
where 
is the displacement vector at position and 
is the normal unit vector for the glass surface.
With a very low strain rate during sagging, the glass behavior was assumed to be viscoelastic around transition temperature Tg. The elastic part is characterized by the instantaneous Young modulus E and Poisson’s ratio v. At a given temperature, the glass is viscoelastic. The stress and strain tensors were split into a deviatoric tensor and a hydrostatic part using following relationships:
(3) 
(4) 
where e(
, t)is the deviatoric strain tensor, s(, t) the deviatoric stress tensor, I the unit tensor, εh(, t) the first strain tensor invariant and σh(, t) the first stress tensor invariant. In the following, a generalized Maxwell model was considered for the shear part. This leads to:
(5) 
with shear modulus 
is the weight at relaxation time τi and n the number of relaxation times used to describe the behavior of the glass (n is generally equal to 6). Bulk modulus K was used as a constant. Variable ξ was the so-called “reduced time”, which was used to take temperature dependence into account through the thermo-rheological simplicity assumption. It is defined by:
(6) 
where φ is the “shift function” defined by (8). The behavior of glass during a cooling process is very complex since structural relaxation must be taken into account. This is usually done by using the concept of fictive temperature. Roughly speaking, fictive temperature Tf represents the deviation of the structure of the glass from its equilibrium state. The fictive temperature is determined as follows:
(7) 
M(t) is the relaxation modulus of the fictive temperature, which depends only on the material. Shift function [24] is defined by:
H is an activation energy, Rg = 8.314 J · mol−1 · K−1 the universal gas constant, x a material parameter and Tr a reference temperature at which G(t) is measured.
Moreover, the fictive temperature Tf also contributes to thermal strain:
βl and βg are the dilatation coefficients in the liquid and solid states respectively.
Since gravity sagging will create large deformations of the glass sheet, the highly nonlinear system (1,2) must be solved to get the displacements. The high temperature dependency of the behavior of the glass (8) means that, in the modeling, we must consider the heat transfer in the glass during gravity sagging. The heat transfer in the 2-D glass sheet is described by the well-known heat transfer equation:
T denotes the temperature depending on position and time t, cp is the specific heat capacity, and kh the heat conductivity. T0() denotes the initial temperature of the glass. On the right hand side of (10), 
(, T) denotes the radiative flux vector, which is defined as the first moment of radiative intensity 
) with respect to direction vector 
by:
(12) 
S2 denotes the unit sphere. At the boundary, it is proposed to describe heat flux using:
This means that heat flux is composed of two terms. The first one represents convection with the surrounding air at a temperature T∞. α is the convection coefficient. The second represents the difference between the radiation of the glass and the irradiation of the surroundings in the opaque wavelength region, ε denotes hemispherical emissivity and λ the wavelength in the glass, λmin ≤ λ ≤ λmax. B(T(, t),λ) denotes Planck’s function given as
k is Boltzmann’s constant, h Planck’s constant, c0 the speed of light in a vacuum, and ng the refractive index of the glass. Due to the glass sheet forming, the heat transfer problem (10–14) must be solved on a deformable body. There is coupling between the equilibrium equations (1–9) and the heat equations (10–14).
The glass sheet is locally heated ( 
∂Lq) using a laser with a power denoted by qlaser(t) acting on a surface area denoted by Alaser. The radiation of the hot glass in the semi-transparent wavelength region is described by the Radiative Transfer Equation (RTE) [22]. If a band model for the absorption coefficient is considered:
(15) 
Then the radiative intensity must satisfy the following equation in each band k[23]:
with the boundary condition for 
on domain ∂Lqt affected by the laser heating
and the boundary condition elsewhere
In equations (16) to (18), the following are used:
Due to the RTE (16–18), the whole system (10–18) used to compute the temperatures in the glass sheet during glass forming is high-dimensional, non-linear, and therefore, challenging to solve numerically. This is why we use the P1-Approximation for the radiative part (16–18). Instead of (16) with boundary conditions (17) and (18), this approximation leads to
where
In Equations (10–17), ∂D and ∂Lq describe the surface of the deformable glass sheet, regardless of the state of deformation being considered during gravity sagging.
The divergence of the radiative flux in (19) can then be computed using
The derivation of the fundamental equations and the P1-Approximation can be found in [7]. For the mathematical analysis of the coupled system (10, 11, 16), (19, 20), (21, 22), refer to [23].
To solve the RTE using P1-Approximation, one can note that (19) is very close in form to the steady-state heat transfer equation obtained by deleting the left-hand side in (10). Not only does the equation have a similar form, but the boundary conditions do as well (comparing (19) to (13)). Only term κkGk(), present in the P1-Approximation (19) has no equivalent in the steady-state heat transfer equation (10). Consequently, the decision was made to solve the P1-Approximation using ABAQUS® finite element software because of its ability to solve the steady-state heat equation. For this reason, using a DC2D8 thermal finite element, which is an 8-node quadrangle element with bilinear interpolation in ABAQUS®, the RTE can be solved using the P1-Approximation (19) if, at each finite element level:
) is considered as a temperature to be computed for each of the nodes of the finite element,
represents material conductivity,
and the film temperature equal to Gka(),By acting at each finite element level according to the aforementioned considerations to solve the RTE using P1-Approximation, the UEL (User ELement) user-subroutine in ABAQUS® was employed to add extra term κkGk() present in (19) and not present in the steady-state heat transfer equation (10) to the thermal stiffness matrix. With the DFLUX user-subroutine, the temperature dependent right-hand side in (19) was introduced.
The P1-Approximation (19) must be solved for each band using the aforementioned procedure in ABAQUS®. Afterwards, the divergence of the radiative flux ∇ · (, t) (23) is computed to determine the heat transfer in the glass sheet (10–13).
To summarize, for a given temperature map in the glass sheet (meaning for a given time t and for given temperature values at the nodes of the 2-D finite element mesh of the glass sheet), the body caused by radiation is computed as follows:
, t)) (14) to be used in (13) are computed. The surface radiation effects appearing in the second term of the boundary conditions (13) used for the heat equation (10) can be directly computed (a first FORTRAN program was developed for this purpose).) (21–22) for the Mk bands is computed using functions Bk(T(, t)) (this is done using ABAQUS® according to the aforementioned procedure)(, t) is computed using results Gk() of the Mk bands (23) (a second FORTRAN program was developed for this purpose).The method to solve the RTE with the P1-Approximation (19–22) using the ABAQUS® finite element software as described above was validated with a one-dimensional (1-D) problem by using just one band κk = 10. m−1 to get an analytical solution. The 1-D problem was solved using 2-D rectangular geometry with insulation on the two horizontal boundaries and boundary conditions on the vertical edges described by (20). The analytical solution for G(x) and the ABAQUS® solution were consistent with each other (difference of less than 0.01% at each node) using a mesh of 50 uniform elements in the 1-D solution. Several κk values were tested. The conclusion was that the mesh must be refined near the boundaries as κk increases (i.e., for the more opaque bands). This is due to the appearance of boundary layers when κk becomes large. The coupling of the RTE with the heat transfer equation (10–13) was also validated using a 1-D thermal problem (10–13) considering a 2-band model for radiation. The value of the surface temperature was successfully compared with the solution of a 1-D finite difference program developed by Siedow et Al. [15]. In this program, the P1-Approximation was also implemented (difference less than 0.25% for the surface temperature at each time step).
In this section, different methods to take radiation into account are compared and their influences on the temperature changes during the process, on the deformed shape of the glass sheet and on the residual stresses.
The geometry described in Figure 1 was considered with thickness w equal to 6.· 10−3m and length l to 150.· 10−3 m. The initial uniform temperature of the glass was T0() = 873.15 K. All the mechanical, thermal and radiative data are given in Appendix. The modeling is divided into two parts. During the first 25 s of forming, the laser heating occurred at a distance d (Figure 1) of 5.· 10−3 m. The characteristics of the laser were: width δ = 3.· 10−3 m, surface flux 
= 250. kW · m−2, and wavelength 2.75 μm ≤ λlaser ≤ 4.50 μm. Note that the laser wavelength was chosen in the k = 2 band of the band model for the glass absorption ((15) and Table 5). Natural convection around the sheet was considered with the temperature of the surrounding air T∞ equal to 873.15 K and convection coefficient α equal to 20. W · m−2 · K−1. For the next 250 s of tempering, the laser was switched off and the glass sheet was cooled down by forced convection described by T∞ = 273.15 K and α = 300. W · m−2 · K−1.
Three different approaches were considered in the modeling. They differ in the way the laser source, internal radiation and surface radiation are taken into account.
*Case 1, denoted “Surface”: the laser heating was modeled using a surface flux q” applied on a zone of boundary ∂Lqt defined by 
and centered on the laser entry point located at (d, w). Term q” was added to the right-hand side of (13) only on the area concerned by the laser flux. The surface radiation was taken into account with Stefan-Boltzmann’s law [9]. The radiation effects inside the glass were modeled with the Rosseland’s approximation [14] using an equivalent conductivity ke computed from the band model used and the corresponding absorption coefficients (Table 7 in Appendix). In Case 1, body flux ∇ · (, t) in equation (10) was not computed and taken as equal to zero.
*Case 2, denoted “Beer”: the laser was modeled with Beer’s law [20] considering body flux Q(z) = 
obtained by solving the 1-D RTE [20]. In the present 2-D modeling, it was assumed that the glass region concerned by the heat flux was defined as 0 < z < w and – 
< y < d + ; Alaser = δ · 1 (1 is the dimension in the x –direction). Constant κ2 is the absorption coefficient of the glass at laser wavelength λlaser; κ2 is equal to 330. m−1 considering Table 5 in Appendix. The body flux Q(z) is added to the right-hand side of (10). In case 2, the surface radiation was also taken into account with the Stefan-Boltzmann’s law and, the radiation effects were taken into consideration using the Rosseland’s approximation. As in Case 1, body flux ∇ · (, t) in equation (10) was not computed and taken as equal to zero.
* Case 3, denoted “P1”: the laser and the radiation effects inside glass are modeled with the P1-Approximation using a three-band-model for the absorption coefficient (Table 5 in Appendix). The resolution of the thermal problem and RTE were performed made using the equations [15–23]. Surface radiation was only considered in the opaque region of the spectrum using the second expression in the right-side in (13). Since the laser wavelength belonged to [2.75 μm; 4.50 μm], the laser heat was also considered in equation (16) when k = 2. Consequently, for k = 1 and k = 3, the term qklaser vanished in (16).
The thermo-mechanical problem of the entire glass sheet forming process was incrementally and iteratively solved, using ABAQUS® finite element software. Both mechanical and thermal equations were solved using the mesh in Figure 2. The mesh is composed of 8987 nodes and 2920 CPEG8T elements with generalized plane strain conditions, biquadratic interpolation for displacements and bilinear interpolation for temperatures. Before each time step, the RTE equation was solved for each band using the P1-Approximation (18–23) and the same mesh as in Figure 2, but with DC2D8 thermal elements to solve the RTE with the modifications described in the procedure above.
Forty elements were used in the thickness with refinement near the upper and lower glass surfaces, under the laser location and on the glass edges. The mesh was refined to get a correct estimation of the temperature gradients during the heating (forming) and cooling (tempering) steps. These gradients correspond respectively to the laser heating and to the air convection cooling. Moreover, the refinement provides to get an accurate computation of the radiative energy (19–22), of the bending effects during the forming and of the residual stresses after tempering.
The first analysis concerns the forming step (25 s). Figure 3 shows the temperature changes for P1, Surface and Beer cases in three locations: in the upper surface of the glass sheet throughout the laser entry point touched by the laser (Figure 3-a), in the mid-plane of the glass sheet under the laser (Figure 3-b), and on the lower surface under the laser (Figure 3-c). These three locations are respectively denoted A(d, w), B
and C(d, 0).
Figure 3: Temperature changes during the forming step for all cases (P1 
, Surface 
and Beer 
cases : (a) at Point A (d, w), (b) at Point B 
, (c) at Point C (d, 0).
Because of the position of the laser over the sheet, the highest temperatures are obtained at Point A (Figure 3-a) regardless of the case, and the lowest temperatures are obtained at Point C (Figure 3-c).
Regardless of the point, the temperature values obtained during the forming step in the Surface case are higher than in the Beer and P1 cases. In fact, in the modeling for the Surface case, the laser energy is totally absorbed by the glass and in contrast, in the formulation of the two other models, one part passes through the glass.
At the end of the forming step, at point A (Figure 3-a) close to the laser, the temperature obtained in the P1 case is 25. K higher than in the Beer case. On the other side of the glass sheet, Point C in Figure 3-c, the temperature obtained in P1 is 20. K lower than in Beer. Inside the glass (Points B and C), the temperatures are lower in P1 than in Surface and Beer, As found in [3,5,22], the Rosseland’s approximation is too diffusive.
The difference between Beer and Surface is greater than in Tian et Li works [20]. In fact, they used a laser with a larger laser wavelength (10.6 μm) and at this wavelength, a larger absorption coefficient for the glass (≈ 30000. m−1, i.e. quite opaque glass). In consequence, less of radiation entered in the glass in their studies than in the present study.
At the end of the forming step, the temperature fields plotted on the initial glass geometry for the three cases are presented in Figure 4. The analysis zone is limited to d − 3 · w ≤ y ≤ d + 3 · w (it corresponds to 32 × 10−3 m ≤ y ≤ 68 × 10−3 m). Outside of this zone, the temperature remains equal to the initial temperature with no laser heating effect.
Comparing temperature fields in Figure 4, the zone concerned by the temperature changes at the end of forming is larger for Surface and Beer than for P1. In y –direction, the temperature remains equal to the initial temperature (873. K) for P1 at a distance 1.5 · w under the laser entry point, for Surface and Beer, the distance is over 2.5 · w. There is also more homogeneous temperature distribution in z –direction for Surface and Beer compared with P1. Once again, the Rosseland’s approximation is too diffusive.
Figure 5 gives the displacement in z –direction of the glass sheet middle line at the end of the forming step (25 s). With the differences in the temperature fields in Figures 3 and 4, the computed displacements on the midline are different for the P1, Surface and Beer cases.
Figure 4: Temperature fields at the end of the forming step in the glass zone defined by d − 3 · w ≤ y ≤ d + 3 · w (32 × 10−3 m ≤ y ≤ 68 × 10−3 m): (a) P1, (b) Surface and (c) Beer cases.
Figure 5: Representation of z-displacements along the glass sheet middle line defined by z = w/2, 0 ≤ y ≤ 1 (0. ≤ y ≤ 150. · 10−3m, z = 3 × 10−3m) for P1 
, Surface 
and Beer 
cases.
Regardless of the case, displacement is a combination of the temperature level in the glass and of the zone concerned by the temperature change. Since glass viscosity depends on the temperature (3–8), the viscosity level is more or less reduced in a large zone of slightly varying size and glass sagging occurs. Based on the computations, one may observe that the displacements along the vertical z –axis at glass sheet edge y = 150 × 10−3 m are respectively 8.9 × 10−3 m, 11.6 × 10−3 m and 17.2 × 10−3 m for the P1, Beer and Surface cases. This corresponds to z –displacement with thickness ratios equal to 1.48, 1.93 and 2.87. The Surface case produces the largest displacements and P1 the smallest. This can be explained by the higher heat energy absorbed by the glass in the Surface case. When comparing Beer and P1, the modeling shows that the glass region affected by the highest temperature changes is larger in case of Beer.
The second analysis concerns the tempering step (250 s) when the laser is shut off and the glass is cooled by forced convection. Since no laser heating is present, the term related to qlaser in equations (17) and (22) vanishes. The P1-Approximation (19–22) to solve the RTE remains the same as the forming step. In this tempering step, the P1-Approximation results (denoted “P1”) compared with the results of the Stefan-Boltzmann’s law and Rosseland’s approximation are used (denoted “BS-Rosseland”).
For the tempering, Figure 4 shows the modeling results outside of the zone affected by the laser forming, where the temperatures after the forming step remain uniform and equal to initial temperature 873.15 K. The analysis is made on the line defined by y = 2 · d and 0. ≤ z ≤ w (y = 100 × 10−3 m, 0. ≤ z ≤ 6 × 10−3 m). Due to tempering conditions, there is a symmetry in the temperature distribution with respect to the glass thickness in this location.
The temperature gradient between the core and the surface of the glass sheet is one major parameter for judging the quality of the tempering [15]. Figure 6 shows the change in temperature gradient between surface and core during tempering for the SB-Rosseland and P1 cases are presented. The maximal value of the gradient is 170. K for P1, whereas it is only 80. K for SB-Rosseland. This is once again due to the overestimated diffusion of the Rosseland’s approximation which leads to a smaller difference between the surface temperature and the core temperature. Consequently, the computed residual stresses in the tempered glass are different, as shown in Figure 7, The residual stress distributions along the thickness of the glass present the well-known parabolic shape with tension in the core and compression on the surface. Like the consequences of the different temperature gradients of the P1 and BS-Rosseland cases (Figure 6), the magnitudes of core and surface stresses are different. In the core, tension stress is +13.8 MPa for SB-Rosseland and, +47.7 MPa for P1. On the surface, compression stress is -34.9 MPa for SB-Rosseland and, −89.1 MPa for P1. In [24], for a convection coefficient equal to 320. W · m−2 · K−1 and an initial glass temperature equal to 873.15 K, the experiments performed by R. Gardon and O.S. Narayanaswamy showed +50. MPa tensile stress in the core. With a convection coefficient of 300. W · m−2 · K−1 in the present modeling and the same initial temperature used in [24], we may conclude that:
Figure 6: Temperature gradients between surface (z = 6 × 10−3 m) and core (z = 3 × 10−3 m) during the tempering step in location y = 2 · d = 100.· 10 −3m) for P1 
and SB-Rosseland 
cases.
Figure 7: Distribution of residual stresses along the glass thickness in location y = 2 · d = 100 × 10−3 m at the end of the entire glass sheet forming process for P1 and SB-Rosseland 
cases.
In the present study, 2-D mathematical formulations of the entire glass sheet forming process, including radiation, are proposed. In comparison with the literature, where approximated methods have been developed to account for internal and surface radiation in glass, the Radiative Transfer Equation (RTE) was solved using the P1-Approximation even through the glass domain is not a fixed domain during the forming phase. By using a band-model for the absorption coefficient of glass, the radiative body flux needed to solve the heat transfer equation was computed using ABAQUS® finite element software by developing specific user-subroutines in ABAQUS® and two complementary subroutines in FORTRAN. Comparing the results with approximated methods, it was proven that significant errors exist when the RTE is not correctly solved. These errors not only concern the temperature changes in the glass but also the deformed shape after forming and the residual stresses after tempering.
This research was supported by 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. The authors are most grateful to Prof. G.W. Scherer (Department of Civil Engineering, Princeton Materials Institute, Princeton University) and to Prof. S.M. Rekhson (Cleveland State University) for assisting and guiding us in our work.
[1] D. Lochegnies and R.M.M. Mattheij (Eds.). Modelling of Glass Forming and Tempering. International Journal of Forming Processes. Paris: Hermes Science Publications. ISBN 2-7462-0081-3, 1999.
[2] D. Lochegnies and M. Cable (Eds.). Modelling and Control of Glass Forming and Tempering. International Journal of Forming Processes. Paris: Hermes Science Publications. Vol. 7, N°4. ISBN 2-7462-1041-X, 2004.
[3] M.F. Modest. Radiative heat transfer. Academic Press, 2003.
[4] R. Siegel and J.R. Howell. Thermal radiation heat transfer. Taylor & Francis Inc., USA, 1992.
[5] R. Viskanta and E.E. Anderson. Heat transfer in semitransparent solids. In: Irvine, T.F. Jr., Harnett, J.P. (eds.) Advances in Heat Transfer, 11,317–441, 1975.
[6] D. Krause and H. Loch. Mathematical simulation in glass technology. Schott series on glass and glass ceramics. Springer-Verlag Berlin Heidelberg, 2002.
[7] A. Farina, A. Klar, R.M.M. Mattheij, A. Mikelic, and N. Siedow. Mathematical models in the manufacturing of glass. Lecture Notes in Mathematics, Springer, 2011.
[8] S. Gregoire, J.M.A. Cesar de Sa, P. Moreau, and D. Lochegnies. Modelling of heat transfer at glass/mould interface in press and blow forming processes. Computers and Structures, 85,15–16,1194–1205, 2007.
[9] M. Sellier, M. Breitbach, H. Loch and N, Siedow. An iterative algorithm for optimal mould design in high-precision compression moulding, Proceedings of The Institution of Mechanical Engineers Part B-journal of Engineering Manufacture – Proc. Inst. Mech. Eng. B-J Eng. Ma., 221,1,25–33, 2007.
[10] A. Sarhadi, J. Hattel, H. Hansen, C. Tutum, L. Lorenzen, and P. Skovgaard. Thermal modelling of the multi-stage heating system with variable boundary conditions in the wafer based precision glass moulding process. Journal of Materials Processing Technology, 212,8,1771–1779,2012.
[11 ] U. Fotheringham and F.T. Lentes. Active thermal conductivity of hot glass. Glass Science and Technology, 67,12,335–342, 1994.
[12] D. Mann, R.E. Field and R. Viskanta. Determination of specific heat and true thermal conductivity of glass from dynamic temperature data. Heat and Mass Transfer, 27,225–231, 1992.
[13] G. Dusserre, F. Schmidt, G. Dour, and G. Bernhart. Thermo-mechanical stresses in cast steel dies during glass pressing process. Journal of Materials Processing Technology, 162–164,484–491,2005.
[14] S. Rosseland. Note on the absorption of radiation within a star. M.N.R.A.S., 84,525–545, 1924.
[15] N. Siedow, T. Grosan, D. Lochegnies, and E. Romero. Application of a new method for radiative heat transfer to flat glass tempering. Journal of the American Ceramic Society, 88,8,2181–2187,2005.
[16] J. Jiao and X. Wang. A numerical simulation of machining glass by dual CO2-laser beams. Optics & Laser Technology,40,297–301,2008.
[17] A. Tseng, Y. Chen, C. Chao, K. Ma and T.P. Chen. Recent developments on microablation of glass materials using excimer lasers. Optics and Lasers in Engineering, 45:975–992, 2007.
[18] K. Yamamotoa, N. Hasaka, H. Morita and E. Ohmuraa. Three-dimensional thermal stress analysis on laser scribing of glass. Precision Engineering, 32,301–308,2008.
[19] A.J.C. Grellier, N.K. Zayer, and C.N. Pannell. Heat transfer modelling in CO2 laser processing of optical fibres. Optics Communications, 152,324–328, 1998.
[20] W. Tian and W. K. S. Chiu. Temperature prediction for CO2 laser heating of moving rods. Optics & Laser Technology, 36,131–137,2004.
[21] J.F. Li, L. Li and F.H. Stott. Comparison of volumetric and surface heating sources in the modeling of laser melting of ceramic materials. International Journal of Heat and Mass Transfer, 47,1159–1174,2004.
[22] F.T. Lentes and N. Siedow. Three-dimensional radiative heat transfer in glass cooling processes. Glasstechnische
Berichte, Glass Science and Technology, 72,188–196,1999.
[23] R. Pinnau. Analysis of optimal boundary control for radiative heat transfer modeled by the sp1-system. Communication in Mathematical Sciences, 5,4,951–969,2007.
[24] R. Gardon and O. S. Narayanaswamy, Stress and volume relaxation in annealing flat glass. Journal of the American Ceramic Society, 53,7,380–385, 1970.
Table 1: Elastic properties of glass
| Density | ρ = 2500. kg · m−3 |
| Elastic part | E = 71. GPa |
| v = 0.22 | |
| Glass dilatation coefficient | αg = 9.2 10−6 |
| Liquid dilatation coefficient | αl = 18.4 10−6 |
Table 2: Relaxation properties of glass (5)
Table 3: Data for shift function φ(, t) (8)
| Tref = 746.15 K |
| H/R = 7.65 104 |
| x = 0.5 |
Table 4: Thermal properties of glass
| Conductivity | kh = 1. W · m−1 · K−1 |
| Specific heat | Cp = 1250. J · kg−1 · K−1 |
| Density | ρ = 2500. kg · m−3 |
Table 5: Absorption coefficients
Table 6: Other radiation properties
| refractive index | ng = 1.46 |
| emissivity | ε = 0.92 |
Table 7: Equivalent conductivity
| T(K) | ke (W · m−2 · K−1) |
| 873. | 7.19 |
| 923. | 9.38 |
| 973. | 12.04 |
| 1023. | 15.19 |
| 1073. | 18.86 |
| 1123. | 23.07 |
| 1173. | 27.83 |
