Modeling and Simulation of the Cooling Process of
Borosilicate Glass

Nicolas Barth, Daniel George, Saïd Ahzi, Yves Rémond, Véronique Doquet, Frédéric Bouyer, Sophie Bétremieux

► To cite this version:

Nicolas Barth, Daniel George, Saïd Ahzi, Yves Rémond, Véronique Doquet, et al. Modeling and Simulation of the Cooling Process of Borosilicate Glass. Journal of Engineering Materials and Technology, 2012, 134 (4), pp.041001. 10.1115/1.4006132. hal-00761649

HAL Id: hal-00761649

https://hal.science/hal-00761649v1

Submitted on 3 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY-SA 4.0 - Attribution - ShareAlike - International
License
# Modeling and Simulation of the Cooling Process of Borosilicate Glass

Nicolas Barth, Daniel George, Saïd Ahzi, Yves Rémond

Institut de Mécanique des Fluides et des Solides, Université de Strasbourg, CNRS, 2 rue, Boussingault, 67000 Strasbourg, France

## Véronique Doquet

Laboratoire de Mécanique des Solides, UMR CNRS 7649, Ecole Polytechnique, 91128, Palaiseau Cedex, France

## Frédéric Bouyer

CEA Centre de Marcoule DTCD/SECM, 30207 Bagnols-Sur-Cèze Cedex, France

## Sophie Bétremieux

AREVA, Direction Innovation Recherche et Projets, 1 place Jean Millier, 92084 Paris-La Défense 6, France

For a better understanding of the thermomechanical behavior of glasses used for nuclear waste vitrification, the cooling process of a bulk borosilicate glass is modeled using the finite element code Abaqus. During this process, the thermal gradients may have an impact on the solidification process. To evaluate this impact, the simulation was based on thermal experimental data from an inactive nuclear waste package. The thermal calculations were made within a parametric window using different boundary conditions to evaluate the variations of temperature distributions for each case. The temperature differences throughout the thickness of solidified glass were found to be significantly non-uniform throughout the package. The temperature evolution in the bulk glass was highly responsive to the external cooling rates applied; thus emphasizing the role of the thermal inertia for this bulky glass cast.

Keywords: cooling process, glass casting, thermal gradient, nuclear glass, FEM simulation

## 1 Introduction

Vitrification of high level radioactive waste has been recognized as a secure available technology to dispose of nuclear fission products on time scales going from thousands up to millions of years. With respect to cost and safety reasons, the manufacured nuclear glasses have as a key feature the nuclear fuel post products content mixed within the glassy matrix. The molten glass is mixed with calcined nuclear waste in a crucible melter at temperatures of around 1373 K since fission products have a good solubility at this temperature range. Once the waste content is in this optimum melted condition, the nuclear glass is poured into steel canisters and cooled down.

As the glass cooling process and solidification occur, the temperature gradients induce transient and residual stresses up to a point at which crack nucleation is likely to happen, as observed by Kamizono and Niwa [1] and Faletti and Ethridge [2]. The stress buildup results principally from thermal strains either directly in the solid glass or as a consequence of the mechanical equilibrium between the internal strained liquid and external solidified glass. Viscoelasticity acts as a dissipation effect of these thermal strains during the liquid solid phase transition, at about 800 K, thus reducing cracking during the cooling process. Nevertheless, if crack nucleation occurs, it may modify fission products release rates as soon as the high level waste (HLW) package becomes soaked in its final disposal conditions. The repercussions of the presence of fracture surfaces within the solidified structure have been demonstrated by Perez and Westsik [3] in a nuclear glass leaching process.

The borosilicate glass SON68 is an inert (non radioactive or inac tive) material developed in France as an equivalent reference for the nuclear glasses used in industry. Its chemical behavior and leach ability were studied by Frugier et al. [4]. Thermomechanical crack modeling studies were also undertaken on this material by thermal shocks on small samples [5]. We used this borosilicate glass to eval uate the cracking process only through the analysis of the thermal gradients during the cooling process, thus neglecting viscoelasticity and damage. In the context of the industrial methodology, it would also be necessary to take into account the thermal energy due to the radioactive decay which acts as a reducer of the thermal stresses buildup, and of cracking, as stated by Faletti and Ethridge [2] and observed by Kahl et al. [6]. In this study the radioactive decay energy was not introduced in the proposed modeling.

In many industrial processes concerning inert glasses, thermal simulations of cooling and solidifying glass are common for their applications in industries of tempered glass [7 9] or of float glass [10]. For these cases, thermal gradients are well controlled due to their small fabrication dimensions so that the glass is manufac tured without cracking and with appropriate residual stresses. On the contrary, the HLW vitrification is concerned with large vol umes which generate more complex cooling processes. They induce a behavior of the package mechanical response that is also complex. For the sake of simplicity and representativeness of the HLW package cooling, we neglect the pressure of the liquid glass in its contribution to the thermal stresses (to our knowledge of the literature, this has not been formalized before).

It is therefore crucial to study the coupled damage viscoelastic behavior since it allows for correctly predicting the cracking behavior of the disposed HLW packages over long periods. How ever, studying these behaviors requires a grasp of the temperature gradients evolution within HLW packages which can bring tangi ble interpretation of the package behavior. Hence, purely thermal analyses using the Abaqus finite element code [11] were under taken to study such evolution.

<footnote n="1">¹ Corresponding author.

Contributed by the Materials Division of ASME for publication in the JOURNAL OF ENGINEERING MATERIALS AND TECHNOLOGY. Manuscript received June 30, 2011; final manuscript received February 2, 2012; published online June 27, 2012. Assoc. Editor: Irene Beyerlein.</footnote>
All simulations are based on experimental results related to a full scale experiment where an inactive waste package was made, resulting in the vitrification of the SON68 glass. Then, a parametric study is undertaken for this particular experiment by applying different thermal boundary conditions, depending on the experimental process simulated, or through the modification of the cooling rate after the filling process.

## 2 Cooling Process and Materials

Necessary details of the cooling process of the vitrified waste will be stated in this section as well as the experimental data on which the current study relies.

**2.1 Cooling Process.** Following an experimental process at CEA, a full scale package specimen containing the SON68 glass is manufactured in order to collect a set of temperature boundary conditions results throughout this non radioactive process.

The cylindrical canisters used in this study are 1.3 m high with a diameter of 0.43 m including a wall thickness of 0.005 m. The containers are made of stainless steel and have a sole opening of 0.14 m in diameter at the top. The glass at a temperature of about 1373 K is poured into the preheated canister (at about 773 K) from the crucible melter located above it.

The canister is filled in two phases, each casting being half of the nominal volume. The total volume of filled material is about 400 kg of SON68 whose composition is given by Frugier et al. [4]. After a waiting time of about 2 h below the melter, the package is taken out and placed into a room where the cooling process progresses further. The convection at the container surface is considered as natural (no forced air flows).

**2.2 Experimental Results.** Two thermocouples have been fixed and used to monitor the temperature at the container wall, at

the bottom and top, in order to measure experimental temperature evolution during the whole process. The results from this test are reported in Fig. 1 and show the temperature evolution during casting and cooling process up to point B (2000 min. process time). In this study, these results will be used as boundary conditions to simulate this particular experimental process.

The SON68 material thermal conductivity, heat capacity and density are given by Dubé et al. [5] for temperatures between 293 and 773 K. Thermal conductivity values above 773 K were also calculated from the following formula given by Barthélémy [12]: $k(T) = 3.66 \times 10^{-3} T$ 0.269 for $T > 773$ K with $T$ in [°C]. Heat capacity values above 773 K for the liquid glass were taken constant to the value at 773 K. The thermal properties are referenced in Table 1 for a density of 2845 kg/m³.

## 3 Modeling of the Process

The finite element code Abaqus FEA [11] was used to compute a 3D model of a half canister totally filled (with 21,500 hexahedral elements, 24,000 nodes and 24,000 dof). The half canister is designed as a half cylinder, longitudinally cut, defined by a vertical axis of 1.1 m and a diameter of 0.4 m. It models a package whose boundary is delimited by the canister external surface, except for the upper face of the cylinder which is identified with the top surface of the glass ("free surface") of a nominally filled canister. The geometrical approximations used in this analysis are shown in Fig. 2.

**3.1 Governing Thermal Equations.** For the numerical analysis, energy balance was computed through heat transfer equations. We consider a volume $V$ of isotropic material with the external surface $S$ and without material flow. The classical heat transfer equation within this material is given by:

<figure><img src="image_2.png" /><figcaption>Fig. 1 External temperatures during the filling and cooling of inert glass in the experimental two casting process (positions are near the bottom and the top of the nominal cast, outside on the container). The "heterogeneous" temperature simulations begin from point $S$ onwards – the heterogeneity in the glass being induced mainly by this particular process in two castings.</figcaption></figure>

<table><caption>Table 1 SON68 glass thermal properties from Refs. [5,12]</caption><tr><td>Temperature [K]</td><td>293</td><td>373</td><td>473</td><td>573</td><td>673</td><td>773</td><td>873</td><td>973</td><td>1073</td><td>1173</td><td>1273</td><td>1373</td></tr><tr><td>Heat capacity [J/kg/K] (±7% &lt;873 K)</td><td>800</td><td>800</td><td>875</td><td>920</td><td>950</td><td>1000</td><td>1000</td><td>1000</td><td>1000</td><td>1000</td><td>1000</td><td>1000</td></tr><tr><td>Thermal conductivity [W/m/K] (±12% &lt;873 K)</td><td>1.03</td><td>1.03</td><td>1.08</td><td>1.11</td><td>1.11</td><td>1.15</td><td>1.93</td><td>2.29</td><td>2.66</td><td>3.02</td><td>3.39</td><td>3.76</td></tr></table>
<figure><img src="image_3.png" /><figcaption>Fig. 2 Geometry approximations made for the finite element design of the canister, rendered as a half-cylinder, and symmetry considered in the model (which could have been reduced to an axisymmetry). The arrow points to the cylinder center whose height delimits the two glass castings. The gray shades are the 16 layers of stacked plates used to simulate the glass casting process.</figcaption></figure>

$$ \rho \left[ \frac{\partial H(T(x, y, z, t))}{\partial t} \right]_p = \rho c_p (T(x, y, z, t)) \frac{\partial T(x, y, z, t)}{\partial t} \\ = \nabla \cdot [k(T(x, y, z, t)) \nabla T(x, y, z, t)] \\ + \dot{Q}(x, y, z, t) \qquad (1) $$

where $H$ is the specific enthalpy, $\rho$ the density, $(x,y,z)$ the position in the body, $t$ the time, $T(x,y,z,t)$ the temperature, $c_p(T)$ the specific heat capacity, $k(T)$ the isotropic thermal conductivity, and $\dot{Q}$ the internal heat source.

Apart from initial temperatures, thermal boundary conditions applied to the body of interest are needed for Eq. (1) to be resolved in time through the method described in this section. These boundary conditions can be basically specified for the same area of interest either as prescribed temperatures $T=T(x,y,z,t)$, or as heat fluxes: “q” for a surface flux or “$\dot{Q}$” for a bulk heat flux the heat fluxes flow into/from the body of interest through the concerned surface $S_q$ ($S_q \subset S$) or through the concerned volume $V_Q$ of the body ($V_Q \subset V$).

Classically, the heat flux $q$ throughout the surface $S_q$ could be of the surface convection type:

$$ q = h_c(T - T_0) \qquad (2) $$

and/or of the surface radiation type:

$$ q = a(T^4 - T_0^4) \qquad (3) $$

where $h_c=h_c(x,y,z,t)$ is the film coefficient, $T_0=T_0(x,y,z,t)$ the sink temperatures, eventually different from each other, and $a$ the radiation constant.

As the material is condensed matter for which there is slow dilatation in relation with $T$, we assume:

$$ dU \approx dH \approx c \cdot dT(x, y, z, t) \qquad (4) $$

through the uniqueness of the heat capacities: $c_p=c_v=c$.

$\dot{U}$ is the time rate of the specific internal energy, which is used in the equations from this point on. Yet it has to be converted from this continuous time description to a form which is thereafter integrated over time using the backward difference algorithm in order to solve incrementally the nonlinear system by a modified Newton method [11]:

$$ \dot{U}_{t+\Delta t} = \frac{U_{t+\Delta t} - U_t}{\Delta t} \qquad (5) $$

If the body is geometrically approximated with finite elements, the temperature is interpolated with the shape functions and we obtain the space discretization in the form:

$$ T(x, y, z, t) = N^N(x, y, z)T^N(t), \quad N = 1, 2, \dots \qquad (6) $$

where the array of interpolating functions is $N^N$, and $T^N$ the nodal temperatures.

The standard Galerkin approach assumes that $\delta T$, the variational field satisfying the essential boundary conditions applied to the body, is interpolated by the same functions:

$$ \delta T(x, y, z, t) = N^N(x, y, z) \delta T^N(t) \qquad (7) $$

Equations (1) and (4) lead to a spatial discretization, in the continuous time description form:

$$ \int_V \rho \dot{U} \delta T dV + \int_V \nabla \delta T \cdot \left( j_F(x, y, z, t) \right) dV \\ = \int_{S_q} q \delta T dS + \int_{V_Q} \dot{Q} \delta T dV \qquad (8) $$

where the Fourier heat flux $j_F$ is given by:

$$ j_F(x, y, z, t) = k(T(x, y, z, t)) \nabla T(x, y, z, t) \qquad (9) $$

The Fourier flux will be used to interpret the gradients of the thermal evolution for each studied case.

Using Eq. (7), Eq. (8) becomes the following set of equations for the geometric approximation (as $\delta T^N$ are arbitrarily chosen):

$$ \int_V N^N \rho \dot{U} dV \qquad \int_V \nabla N^N \cdot j_F(x, y, z, t) dV \\ = \int_{S_q} q N^N dS + \int_{V_Q} \dot{Q} N^N dV \qquad (10) $$

Using Eq. (5), and neglecting the slow dependence of conductivity with temperature, the time integration leads to:

$$ \frac{1}{\Delta t} \int_V \rho (U_{t+\Delta t} - U_t) N^N dV + \int_V \nabla N^N \cdot k(T(x, y, z, t + \Delta t)) \\ \times \nabla T(x, y, z, t + \Delta t) dV = \int_{S_q} q(x, y, z, t + \Delta t) N^N dS \\ + \int_{V_Q} \dot{Q}(x, y, z, t + \Delta t) N^N dV \qquad (11) $$

## 3.2 Model's Phenomenology.
The radiation and convection phenomena at the canister's edges were not taken into account since the experimental temperatures (measured from the thermo couples) are used as prescribed boundary conditions ("BC") on the edges. These experimental temperatures include the combined effects of conduction, convection and radiation at the stainless steel wall. Moreover, the influence of the stainless steel envelope was neglected due to its very small influence on the conduction of the prescribed temperatures to the glass interface, so that the SON68 glass is the only material modeled. All thermal calculations are based on the thermal properties presented in Table 1.
In addition to the finite element design, the double cast was simulated using eight layers per casting as seen by the gray shades in Fig. 2 (the total volume filled is modeled by a stack of piled plates). For each layer, the application of the BC is uniform over the layer's surface or volume. In the simulations, the 16 layers were included from the start even if, at first, the upper layers were not filled with glass in the experimental process. The pouring of the glass, and its rise in the container, were modeled by means of boundary conditions control. The BC were synchronized with the experimental process and appropriately applied for each filled layer. The modeling of the glass flow from the melter down to the bottom of the container and the continuous rise of the glass free surface along the whole height of a single layer were only approximated by this protocol.

Thermal studies of glass include several exchange mechanisms such as scattering of thermal radiation or internal emission [13]. In this study, the internal radiative transfers are assumed to be included in the thermal conductivity properties given in Table 1, given the fact that SON68 is a rather opaque borosilicate glass as seen by the sample photographs in Ref. [5]. As the stainless steel is opaque, the influence of internal emission or radiative scattering is assumed to be insignificant at the glass steel interface.

Finally, we need to discuss the thermal exchanges at the glass free surface at the top of the cast (i.e. between the glass and the steel cavity above it). Depending on the number of filled layers in the simulations, from 1 to 16, this surface represents respectively 37% to 8% of the total glass interface. Hence, the thermal responses at the 16 successive free surfaces inside the canister are approximated in the current model. Moreover, the values of tem- perature applied at these internal surfaces are those of the stainless steel envelope (external temperatures).

**3.3 Boundary Conditions.** Five different simulations were undertaken in this study, which differ from their boundary condi- tions. The first one is meant to be as close as possible to the exper- imental setup for which external temperature results are available on the canister's edges. Following are four kinds of BC variations for which the cooling process will be slower or faster than the ref- erence case (experimentally based case). The main differences between the BC are summed up in Table 2, and are also shown in Fig. 3.

The four BC variations are namely “Homogenous low T” where the model simulates the annealing at a uniform temperature at which the glass is scarcely liquid, “Homogeneous high T” where the modeling of an instantaneous single cast is simulated, “linear BC” where the cooling process is towards a better thermal insula- tion of the package during solidification, and finally “quench BC” where the cooling process is towards a quench of the package at ambient temperature during solidification.

**3.3.1 The Reference Simulation.** The boundary conditions described here concern the simulation named “nonlinear BC” or cited as “reference simulation” or “ref. model”. Unless otherwise

specified in the following sections, they are applied to all under taken simulations. These boundary conditions are axisymmetri- cally applied on the 3D model.

The application of the prescribed temperatures on the model is split in two cases depending on the temperature range between the top and bottom of the canister. If the experimental temperatures have a difference lower than 30 K between the top and bottom, the mean temperature is uniformly applied on the whole surface of the cylindrical model (the 16 layer temperatures are the same). This happens after the filling process, as soon as the package is placed in the room where the cooling progresses further (at point A in Fig. 1 or Fig. 3). Otherwise, the 16 temperatures applied on each layer's surface are linearly extracted between the top and bottom temperature values.

The temperatures applied at the beginning of the analysis are those at the end of the preheating of the canister. The external temperatures applied closely follow the experimental temperature data during the process which are represented in Fig. 1. Fifty hours are simulated from the beginning of the first casting so that ambient temperature is obtained at the end.

During the cast, internal temperatures of unfilled layers are set at the corresponding value of the external experimental tempera- tures. As the pouring proceeds, the 16 layers are in due turn instantaneously heated at 1373 K, with the external temperature maintained as a priority as previously described. The internal tem- peratures thereafter only evolve with time according to the ther- mal calculations in the filled layers.

**3.3.2 Internal Temperature Homogenization.** In the two sim- ulations presented in this section, this casting process is not mod- eled. The simulations start with a totally filled canister whose internal state is at the same uniform (or “homogeneous”) tempera- ture. The external boundary temperatures are applied to closely follow the experimental cooling rates of the reference model, which are assumed unaltered by this homogenization (this point is further discussed in Sec. 4.1.1).

* Homogeneous low T (“H low T”): This model simulates the homogeneous annealing of the HLW package at a tempera- ture closely above the solidification temperature, after which this model undergoes the same cooling treatment as in the experimental process. Consequently, the whole glass volume is set at a uniform temperature of 856 K, in this model, at the moment the package is taken into the room where the cooling progresses further (at point A in Fig. 1 or Fig. 3). This value is very close to the mean experimental temperature (<5 K) at this point.

* Homogeneous high T (“H high T”): This model simulates the instantaneous casting of the total glass volume in the HLW package, then undergoing the experimental cooling down process (it would also be the case of an annealing of the pack- age at the uniform temperature of the melter). The whole glass volume is set at a uniform temperature of 1373 K at the end of the second casting. During the waiting time of

<table><caption>Table 2 Main features of the boundary conditions</caption><thead><tr><th></th><th>Ref. model</th><th>H low T</th><th>H high T</th><th>Linear BC</th><th>Quench BC</th></tr></thead><tbody><tr><td>Simulation start, in Fig. 1 timescale [min]</td><td>0 (point S, first casting start)</td><td>777 (point A)</td><td>627 (second casting end)</td><td>0 (point S)</td><td>0 (point S)</td></tr><tr><td>Number of casting (or as thermal initial states of the glass: count of significantly different volumes)</td><td>2</td><td>1</td><td>1</td><td>2</td><td>2</td></tr><tr><td>External cooling rate applied</td><td>Experimental (“exp.”)</td><td>Exp.</td><td>Exp. (top thermocouple)</td><td>Lower than exp. at solidification start</td><td>Higher than exp. at solidification start</td></tr><tr><td>End of the simulation in Fig. 1 timescale [min]</td><td>4938</td><td>4938</td><td>4938</td><td>4938</td><td>4938</td></tr></tbody></table>
<figure><img src="image_4.png" /><figcaption>Fig. 3 External temperatures applied after the filling process for the different cooling down periods. Ref., quench and linear BC simulations begin at point S in Fig. 1 (here this Fig. is at the end of Fig. 1). The ref. simulation follows closely experimental thermocouple values between point A and point B. (I) all simulations are represented in this particular time window, where different temperatures are prescribed. (II) the internal temperature homogenization moment is given by the filled symbols for the corresponding simulations, which then follow the experimental data. (III) the external temperature alteration simulations (linear or quench BC) follow the linear paths marked by an open symbol (in the case of linear BC: chord A–B; in the case of quench BC: the initially vertical path from point A).</figcaption></figure>

150 min which then follows, only the experimental data from the top thermocouple is applied on the whole surface ($t=0$ 150 min in Fig. 3). After this period of time, the surface cooldown rates are prescribed identically to the previous simulations.

### 3.3.3 External Temperature Alteration.
In the last two simulations (linear and quench BC), the thermal history of the filling process before cooling is conserved (point S to point A in Fig. 1) but the surface cooling down process is altered. This procedure is started at the moment the package is taken into the cooling room (point A in Fig. 1 or Fig. 3).

* "linear BC": This model simulates a modification of the process to improve the thermal insulation of the HLW package, thus avoiding an initial air quench. The temperature variation of nonlinear BC is replaced by a 21 h long linear cooling down period defined by section A B on Fig. 3, thus cooling the package slower (as if it was in a thermally controlled oven).

* "quench BC": This model simulates an acceleration of the cooling down period by means of quench methods. For this purpose, the package surface is instantaneously set at 300 K.

## 4 Results and Discussion

The variations obtained within the material during the cooling process depend all on the specified boundary conditions used. These BC are completely dependent on the methodology and assumptions made at the beginning of the cooling. They also represent different cases of manufacturing processes and cannot be assumed as negligible since a strong influence was observed on the cooling evolution of the glass. The manufactured dimensions of the material volume have a strong influence in the process since the thermal gradients, as a function of time, are dependent on the cooling rate of the whole material volume. The results obtained with the different boundary conditions applied will be presented hereafter.

### 4.1 Spatial Characterization of Temperature Distributions.
All simulations show various temperature gradients and maps during the cooling process at all points within the model geometry. In order to evaluate precisely the differences between each case,

temperature profiles are plotted along a diameter through the canister as will be shown in Fig. 4.

#### 4.1.1 Solidification Characterization.
Four reference times and temperature profiles are chosen to represent each solidification process simulated. At each of these reference times, the diameter's positions of the corresponding profile along the cylindrical axis are where the radial temperature difference between center and edge is maximum (which equals the vertical position of the maximum centerline temperature as a uniform external temperature is applied at the solidification stage). First of all, the glass solidification temperature is assumed to be fixed at the representative temperature of 773 K which is shown over the temperature pro files in Fig. 4 (see next section). The four reference times used in our simulations are: (i) as the external solidification starts at external edge (radial position $R = \pm 0.2$ m), the temperature profile curve being fully above the dashed solidification line, (ii) the temperature profile crosses the dashed solidification temperature line at 1/4 length of the solidified canister radius ($R = \pm 0.15$ m), (iii) the temperature profile crosses the dashed solidification line at half the length of the radius solidified ($R = \pm 0.1$ m), and (iv) the temperature profile is fully under the solidification temperature the glass has just completely solidified inside the whole canister. The above cited reference times are presented for each simulation in Table 3.

The results and interpretation of the reference time evolution will be detailed further.

The computed radial thermal fluxes (on the edge surface) are presented in Table 4 at the reference times for the different simulations undertaken. The heat fluxes show the variability of thermal behavior between the different cases. In practice, at a given temperature, all fluxes should be identical if the external boundary conditions are assumed to extract the same radiative and convective flux as the reference case, but we see in Table 4, that it is not the case. Let us assume that the numerical reference case corre lates well with the experimental case. The heat flux at solidification start time (773 K external temperature) is 9848 W/m². For the "H low T" and "H high T" cases, we follow the same external cooldown rates boundary conditions (prescribed temperatures) but fluxes are different although we wanted to be close to the
<figure><img src="image_5.png" /><figcaption>Fig. 4 Series of temperature profiles at reference times (i) to (iv) during solidification for each BC case. The experimental temperatures on the edges are taken for each reference time as the mean value of the experimental curves from the two thermocouples at the same moment on the time axis (see Fig. 3).</figcaption></figure>

<table><caption>Table 3 Reference times [minutes right from the end of the cast] for each case of cooling process (simulations “H low T” and “H high T” actually begins at t = 150 and t = 0 min, respectively, in Fig. 3)</caption><thead><tr><th>Ref. times</th><th>Ref. model</th><th>H low T</th><th>H high T</th><th>Linear BC</th><th>Quench BC</th></tr></thead><tbody><tr><td>(i)</td><td>159</td><td>161</td><td>159</td><td>386</td><td>150</td></tr><tr><td>(ii)</td><td>273</td><td>215</td><td>294</td><td>516</td><td>198</td></tr><tr><td>(iii)</td><td>374</td><td>308</td><td>400</td><td>593</td><td>295</td></tr><tr><td>(iv)</td><td>478</td><td>408</td><td>508</td><td>754</td><td>367</td></tr></tbody></table>

<table><caption>Table 4 Maximum radial heat flux component, in [W/m²], of the surface elements at the reference times for each cooling case; the corresponding external temperatures, prescribed, are on the edge of the temperature profiles curves in Fig. 4</caption><thead><tr><th>Ref. times</th><th>Ref. model</th><th>H low T</th><th>H high T</th><th>Linear BC</th><th>Quench BC</th></tr></thead><tbody><tr><td>(i)</td><td>9848</td><td>3813</td><td>11400</td><td>1958</td><td>11680</td></tr><tr><td>(ii)</td><td>5343</td><td>4753</td><td>5212</td><td>1759</td><td>11550</td></tr><tr><td>(iii)</td><td>3573</td><td>3209</td><td>3559</td><td>1696</td><td>5922</td></tr><tr><td>(iv)</td><td>2785</td><td>2520</td><td>2789</td><td>1694</td><td>4203</td></tr></tbody></table>

experimental cooling behavior. This phenomenon is due to the ho mogenized temperature inside the glass which modifies the ther mal equilibrium and generates different heat fluxes at the edge for this reference time (i) cited as an example. Identical observations can be made for the two other heterogeneous cases although the applied boundary conditions are different. For example, the low cooling rate defining the linear BC case generates a much lower heat flux than the other cases and the corresponding reference time (i) is clearly later than the other cases.

In all these cases, the fluxes present an intensity decrease. This is particularly expected for the reference case where the experi mental measurements show a loss of internal energy in the pack age under radiative and/or air cooling processes (see the monotonically decrease of experimental temperature evolution in Fig. 3). These fluxes are all directly dependent on the applied external BC. However, these results enable us only to partly repre

<table><caption>Table 5 Maximum radial heat flux component, in [W/m²], of the solidification front's interfacial elements at the reference times for each cooling case; the corresponding temperature is 773 K</caption><thead><tr><th>Ref. times</th><th>Ref. model</th><th>H low T</th><th>H high T</th><th>Linear BC</th><th>Quench BC</th></tr></thead><tbody><tr><td>(i)</td><td>9848</td><td>3813</td><td>11400</td><td>1958</td><td>11680</td></tr><tr><td>(ii)</td><td>4138</td><td>2137</td><td>4287</td><td>1248</td><td>8029</td></tr><tr><td>(iii)</td><td>1981</td><td>1431</td><td>2076</td><td>760</td><td>3540</td></tr><tr><td>(iv)</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr></tbody></table>

sent the thermal evolution of the cooling process. In Table 5, we present the fluxes inside the package for a radial position at the solidification front interface (which is not the outer edge except for reference time (i)). Incidentally, the intensity decrease of these fluxes during solidification is mainly due to the temperature gradient of the Fourier flux Eq. (9). Both values from Table 4 and Table 5 help us shape the thermal gradient defined in the solid phase.

### 4.1.2 Temperature Profiles at Solidification Reference Times.
Since the thermal gradients strain the package from a stress free reference time assumed to be the reference time (i), the evolution of the thermal gradients in time at the beginning of the solidification induce small total strain. Later on, the shape of thermal gradients spatially defined over the radius may evolve during the cooling period inducing various mechanical responses of the solid liquid equilibrium or of the solid phase during the cooling process.

The temperature profiles associated to the times in Table 3 (and the above mentioned fluxes represented by the slope of the temper ature profiles times the conductivity) are presented in Fig. 4. These profiles are plotted for the different boundary conditions and in the center they correspond to the BC of the reference case (experimen tal). Each graph presents the temperature profile through the diam eter of the package (from R = 0.2 m to + 0.2 m).

The evolution of temperature profiles for each case is very spe cific and shows various thermal gradients as a function of time (from left to right in Fig. 4):

* For the linear BC cooling rate case, the temperature profiles through the canister diameter are almost identical at all
specified times (slow cooling and nearly identical thermal fluxes). The cooling rate at all points of the model is identical. There is no evident temperature gradient evolution between the center and the edge of the package as a function of time during the solidification.

* The “H low T” case shows an evolution different from all other cases because the thermal history of the cast process is completely ignored. In this case, the thermal internal state for the reference time (i) is close to the homogeneous start temperature. Here, the simulated thermal flux is lower than the reference flux at (i) (see Table 4), which proves the fact that the cooling of “H low T” is underestimated with the pre scribed temperature boundary conditions.

* The reference case shows the narrowest gap in our simulations between physical cooling and numerical results. The model takes into account the thermal history of the cooling process with thermal gradients and edge cooling due to the external experimental boundary conditions applied.

* The “H high T” case shows results close to the reference case. This is expected since these cooling processes have some similarities in their prescribed BC. In practice, the cooling should involve the same radiative and natural convection thermal fluxes at time (i) as the experimental simulation. Here, the simulated thermal flux is greater than the reference flux at time (i), which proves the fact that the cooling of “H high T” is overestimated.

* Finally, the quenched case shows temperature distributions rather similar to the two previous cases but with a much faster cooling rate; thus extracting most rapidly the internal thermal energy through the highest thermal fluxes as reported in Table 4.

In addition, the experimental temperatures (mean measure of the thermocouples at the same time on the edge of the container) are plotted in Fig. 4 for comparison with our numerical results. As expected for the reference case, “H low T” and “H high T” cases, there is a direct correlation between the applied temperature boundary conditions and the one observed on the numerical model since the model closely follows the experimental data. However, this is not the case for the two other simulations (linear and quench BC) since the cooling processes are modeled with different external thermal behavior, as represented in Fig. 3. For the linear cooling, the edge experimental temperatures are lower than the numerical results because of the slower cooling rates. On the contrary, for the quench case, the corresponding temperatures for each reference times are higher than those of the experimental temperatures.

### 4.1.3 Results and Discussions from the Literature.
Results from the literature are scarce for similar processes. Faletti and Ethridge [2] and Kahl et al. [6] provide typical temperature pro files (“T profiles”) for canister diameters of 0.3 to 0.6 m comparable with Fig. 4. Around their solidification start, Faletti and Ethridge [2] report a 150 K difference between centerline and edge temperatures, whereas it is near to 350 K in our reference simulation (see Fig. 4 and it will also be displayed in Fig. 5). As the edge temperature drops to 473 K the authors’ profiles display an 85 K difference, while it is still high at 300 K in our reference case. Faletti and Ethridge [2] also reviewed various experiments concerning T profiles that have a 773 K wall temperature, which corresponds well to our solidification start (i). Numerous temperature differences as high as 300 K are reviewed, which are closer to the value of our reference simulation (345 K).

Kahl et al. [6] reported thermal gradients over 0.15 m radius up to a maximum value of 767 and 1000 K/m for active packages (one being thermally insulated). Initially, at the thermal equilibrium in their annealing oven (at a high temperature around 1200 K uniform external temperature), initial radial gradients existed in the package, due to the internal source thermal term, which underwent a difference value of 433 and 667 K/m respectively. These thermal gradients are lower in comparison with our

<figure><img src="image_6.png" /><figcaption>Fig. 5 Temperature differences evolution over constant radius positions; the solidification front in the three simulations modeled is progressing from 0 up to its end between 251–349 min</figcaption></figure>

simulations, as will be illustrated later in Figs. 6 and 7, even if Riley et al. [14] mentioned that these gradients increase as the canister diameter increases. Their results also showed a decrease of thermal gradients with the addition of radioactive material into the canister. This is in good agreement with the assumption made in this work where our large thermal gradients and fluxes are considered as extreme cases whereas the linear cases tend to be close to these particular cooling rates.

As for the tempered glass process, Daudeville and Carré [7] specify thermal gradients in thin plates up to 150 K between the center and edge, which corresponds to 49,180 K/m. This value, which is much higher than in our simulations, is reported for plates over a half thickness of 3.05 mm, whereas the maximum gradients recorded in our simulations (between 1300 and 9500 K/m) are obtained at (ii) (see Figs. 6 and 7) for a solid glass thickness of 50 mm.

The amplitude of the thermal stresses in the experimental cooling process is strongly dependent on the definition of a stress free reference, inherited from the liquid state. This concept is inti mately linked to the quantity of internal liquid and its thermal strain potential as a trapped material phase enclosed by the solid glass. This strain behavior can even be time dependent, as for the case of specific volume evolution during the liquid solid transition that was illustrated for example by Kong et al. [9].

A realistic approach resulting from thermomechanical models is reported by Daudeville and Carré [7] for tempered soda lime silicate glass plates. At the beginning of the cooling, they report that the temperatures decrease more rapidly at the surface than in

<figure><img src="image_7.png" /><figcaption>Fig. 6 Maximum radial temperature differences averaged in the solid glass (over the relevant radial distance $d$) for the reference and homogeneous simulations. The three consecutive first points are respectively at reference times (ii), (iii) and (iv); the continuous lines being used for the totally solidified package only ($d = 0.2$ m).</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Fig. 7** Maximum radial temperature differences averaged in the solid glass (over the relevant radial distance *d*) for the reference and altered boundary conditions simulations. The three consecutive first points are respectively at reference times (ii), (iii) and (iv); the continuous lines being used for the totally solidified package only.</figcaption></figure>

the inner part, which leads to a compression of the core and tran-
sient tensile stresses at the surface (these stresses could be criti-
cal). As the solidification proceeds, the core continues to cool and
contracts whereas the surface is undergoing much less thermal
strain in the solid state. This phenomenon ultimately leads to re-
sidual stresses for which the surface is in a compression state in
equilibrium with the core that has built up tensile stresses.

Since a depression could arise in the liquid as a compressive
stress builds up in the solid, we prefer to adopt in our study
another simplifying stress free reference where the liquid does not
contribute to the thermal strains. Thermal gradients in the current
study were considered in the solid phase, i.e., they are only spa-
tially defined, contrary to the definitions mentioned above which
are dependent on the relative evolution of the core and surface
temperatures. Thus, it is assumed that the liquid phase is always
in a pressure free reference; for example by preserving a liquid
free surface at the top of the glass castings during solidification, or
a gas cavity at standard pressure in the liquid core.

**4.1.4 Spatial Interpretation of the Thermal Gradients.** The internal temperature maps of the studied homogeneous simulations are symmetrical in relation to the cylinder center and axis because this is also the symmetry of all their BC applied. Inciden tally, all of the profiles for the two homogeneous simulations are then positioned at the center of the cylinder, as particular care was taken to only apply uniform temperature at the whole cylinder sur face after the temperature homogenizations.

It is not the case for the reference simulation, which is consid-
ered as “heterogeneous.” At an early stage of the solidification,
the second top casting of the reference simulation has higher
temperatures then more elevated radial thermal gradients than
the first casting. This phenomenon illustrates the unsymmetrical
heterogeneity (far from the bottom and top edge) of T profiles
along the cylinder axis, whereas it is symmetrical (in relation to
the cylinder center) for the homogeneous simulations.

Thermal gradients evolution are assumed to be maximum at the beginning of solidification since at some point the cooling down rate of the centerline is higher than the cooldown rate on the edge. This relative variation happens in the reference simulation as soon as the curves decrease in Fig. 5. It would well correspond to the case where the residual internal tensile stresses begin to build up in the inner part in the Daudeville and Carré [7] approach discussed previously, except that the liquid is assumed to be in a stress free state in our hypotheses. In all our simulations under taken, the tensile stress transient buildup is assumed to be maximal in the ascending part of the curves in Fig. 5.

The space scale is also the major difference between our pro-
cess involving bulky glass and other industrial processes where

one of the glass product dimensions is reduced to the centimeter
scale. Their glass thermal inertia should be then a lesser problem
compared with the studied vitrification process.

**4.2 Time Evolution of Temperature Distributions.** From the temperature profiles (T profiles) introduced above and plotted versus radial distances, we define temperature differences along a section or the entire radius of the canister. These differences are generally plotted as a function of time and called “ΔT profile.” These temperature differences are given at unspecified times of the cooling process, not necessarily at the reference times.

Regarding the validity of the assumptions described for this
model, and as the envelope temperature is the same at the differ-
ent solidification starts (i), the differences of thermal fluxes in the
first line of Table 4 between the altered external cooldown rate
simulations (linear BC or quench BC) and the reference simula-
tion are expected to be the differences in the convection fluxes.

Radiative and convective fluxes (calculated with sink tempera-
tures of T₀ = 300 K) are hereby compared to the corresponding
computed thermal fluxes at reference times (i) to (iv) of our exper-
imental simulation. For the radiative flux, the emissivity taken is
that of the stainless steel which is of 0.22 to 0.45, respectively
between 373 to 698 K [15]. For natural convection, the film coef-
ficient used is $h_c = 1.3[(T \ T_0)/0.4]^{0.25}$ [16].

For (i) and (ii) the difference between reference simulation and temperature extrapolated fluxes are within a 2000 W range. For (iii) and (iv) this difference is within 1000 W. This is a relatively good match (<40%) given the fact that the chosen thermal radia- tive and convective phenomena are not parametrically tuned to match the exact case of the experimental thermal flux behavior.

In the following two sections, we show specifically the effect of
the temperature gradient evolution versus time in the propagation
of the solid phase.

**4.2.1 Internal Temperature Homogenization.** In Fig. 5, the thermal responses of the two castings in the reference simulation are compared with both “H low T” and “H high T” homogenized simulations. Five centerline to edge ΔT profiles are plotted versus time regardless of the glass state (glassy or liquid). We use the same position for each homogenized case (the center’s height). Moreover, to illustrate the thermal response in the reference simulation, three particular ΔT profiles are recorded at three vertical positions on the cylinder axis: at the first quarter, half and last quarter of the glass height. In the container, these positions correspond to the height of the first casting center, the package center height (also corresponding to the glass level between castings) and the height of the second casting center, respectively.

In Fig. 5, all the ΔT profiles present an increase in the begin-
ning of glass solidification, signifying a cooldown rate that is
higher at the surface than at the corresponding point on the center
line. Towards the end of solidification, and well after it in the sol
identified package, the ΔT profiles monotonically decrease towards 0
(the cooldown rate is higher at the core than at the surface during
this period). This evolution happens quite simultaneously between
the three positions recorded along the two castings of the refer-
ence model and also between the reference and homogenized
simulations.

In Fig. 5, the centerline to edge plots enable us to draw compar-
isons between some specific positions in the models during the
solidification, which is the period where these global thermal gra-
dients are highest. As expected, the reference case is located in
between the two other cases “H low T” and “H high T” with the
second casting closer to the “H high T” and the first casting closer
to the “H low T”. Interestingly, the temperature gap between each
ΔT profile is about the same.

To evaluate the maximum impact of temperature gradients on
the physical behavior of the glass, maximum ΔT profiles are com-
pared after averaging over the solidified length *d* in the models.
The reference and the homogenized simulations are represented in
Fig. 6.
A high decrease of solid glass thermal gradients takes place in the beginning of solidification. This is opposite to the Fig. 5 where the ΔT profiles started to increase before decreasing.

Moreover, the unsymmetrical heterogeneity along the height for the nonlinear BC simulation tends to disappear after 10 h from solidification start due to the uniform temperatures imposed on the edge. It becomes rather symmetric in relation with the cylinder center (not represented here). Ultimately, the glass initial internal state tends to be insignificant because the boundary conditions are applied the same way over sufficient time, as illustrated by the superposition of the homogeneous and nonlinear BC models' curves in Figs. 5 or 6.

Concerning the gradients only taken in the solid phase, Fig. 6 clearly shows how nonlinear BC, more precisely its top casting during the solidification, and homogeneous high T models dis play similar maximum thermal gradients. The real cases of an annealing experiment and a single cast experiment would prob ably draw the maximum thermal gradients even closer to the case of the maximum thermal gradients in our reference simula tion. This is because the homogenized simulations yield thermal fluxes which are underestimated and overestimated for annealing and single cast models respectively. The maximum thermal gra dient buildup between (i) and (ii) proves the major roles of local effects and also of the boundary conditions, during cooling, as there is a relatively small influence of the initial extent of the in ternal thermal energy.

### 4.2.2 External Temperature Alteration.
In Fig. 7, maximum temperature differences are represented for the two externally altered cooling down rate simulations, and compared with the ref erence model.

The maximum gradient observed is due to the quench simula tion at reference time (ii) for which it is about 9500 K/m. After the reference time (iii) for the quench simulation, we observe a behavior similar to the reference model after refer ence time (ii).

Concerning the linear BC model, in comparison with the other simulations, the thermal gradients are at a very low level (around 1000 K/m) but for a longer period of time. Although the gradient is small at the beginning since the cooling rate is lower it shows thermal gradients higher than the other heterogeneous sim ulations after 600 min from solidification start.

The evolution of the maximum thermal gradients is also very different in linear BC from the other cases. Due to an edge effect in the maximum local slopes of the T profile linked to these gra dients in the solid phase only, the highest thermal gradients of this simulation are not necessarily arising near the stainless steel envelop, contrary to all the other simulations undertaken.

The internal temperature maps for the last two heterogeneous simulations, not represented here, are qualitatively similar to those of nonlinear BC, although the dynamics are very different. Within the framework of the current applied boundary conditions, a par tial modification of the final external cooling rates as observed in Fig. 7 modifies greatly the thermal gradients of the heterogeneous simulations in comparison with the modification of the internal state for the homogeneous cases.

The experimental nonlinear BC model, with inert glass, is near the conditions of a rapid quench except for the initial time of sol idification. As for most of the values quoted from literature in Sec. 4.1.3, they are mainly closest to the linear BC model where particular care was taken to simulate a thermal insulation of the package during the most critical cooling.

Concerning this linear BC model, in comparison with the other simulations, the thermal gradients are at a very low level but for a much longer period. Indeed, modeling an increase in thermal insulation induces thermal gradient reduction. In this case, the thermal insulation simulation during the cooling down period presents the lowest gradients in the current study, which proves the point of controlling the thermal insulation during the solidification.

# 5 Conclusion

A thermal finite elements simulation was developed to model the filling with glass in a stainless steel container with a scale of 1:1 by taking into account the cooling process conditions. The aim is to better understand the impact of the thermal gradient on the whole final glass structure. The internal thermal state through out the process has been calculated using the experimental top and bottom external temperature measures and process parame ters. Radial thermal gradients were then extracted to investigate the stress buildup and to study the influence of boundary condi tions of the inactive waste package. Homogeneous internal tem peratures simulations were observed to be incapable of accurately representing the intrinsic internal temperatures map of the pack age during our vitrification process.

Two kinds of parametric investigations were made using the reference simulation: assuming firstly different uniform internal temperatures and secondly different external cooling rates. This parametric window allowed us to understand the influence of the cooling process parameters on the thermal gradients in the bulk glass. The package showed a significant reduction of the thermal gradients through the simulation of the annealing above the sol idification temperature, and an even higher decrease of thermal gradients through the simulation of the thermal insulation control.

A complementary work of the presented thermal study would be to model the thermal energy increase due to the radioactive decay, which acts as a reducer of the thermal stresses buildup. This would make it possible to be closer to the real conditions of an active waste.

# Acknowledgment

The authors wish to acknowledge the CEA (LCLT) for the financial support as well as providing the experimental data. The ANDRA (French National Radioactive Waste Management Agency) and AREVA NC are also gratefully acknowledged for supporting this study.

# Nomenclature

## SI Units

$H$ = specific enthalpy [J/kg]

$N$ = $N^{th}$ node in the geometrical approximation of the body

$N^N$ = array of interpolating functions

$\dot{Q}$ = internal heat source [W/m$^3$]

$S$ = surface of the body [m$^2$]

$S_q$ = surface concerned with thermal flux $q$ [m$^2$]

$T$ = temperature [K]

$T_0$ = sink temperature [K]

$T^N$ = nodal temperatures [K]

$U$ = specific internal energy [J/kg]

$\dot{U}$ = time rate of the specific internal energy [W/kg]

$V$ = volume of the body [m$^3$]

$V_Q$ = volume concerned with thermal flux $\dot{Q}$ [m$^3$]

$a$ = radiation constant (Stefan Boltzmann constant times the emissivity) [W/m$^2$/K$^4$]

$c$ = specific heat capacity [J/kg/K]

$c_p$ = specific heat capacity at constant pressure [J/kg/K]

$c_v$ = specific heat capacity at constant volume [J/kg/K]

$k$ = thermal conductivity [W/m/K]

$h_c$ = film coefficient [W/m$^2$/K]

$j_F$ = Fourier heat flux [W/m$^2$]

$\dot{q}$ = surface heat flux [W/m$^2$]

$t$ = time [s]

$(x, y, z)$ = 3D position vector in the body

$\delta T^{(N)}$ = variational field in the Galerkin approach (at nodes)

$\Delta t$ = time increment [s]

$\rho$ = density [kg/m$^3$]
References

[1] Kamizono, H., and Niwa, K., 1984, "An Estimation of the Thermal Shock Resistance of Simulated Nuclear Waste Glass Under Water Quenching Conditions," *J. Mater. Sci. Lett.*, **3**(7), pp. 588-590.

[2] Faletti, D., and Ethridge, L., 1988, "A Method for Predicting Cracking in Waste Glass Canisters," *Nucl. Chem. Waste Manage.*, **8**(2), pp. 123-133.

[3] Perez, J. M., Jr., and Westsik, J. H., Jr., 1981, "Effects of Cracks on Glass Leaching," *Nucl. Chem. Waste Manage.*, **2**(2), pp. 165-168.

[4] Frugier, P., Gin, S., Minet, Y., Chave, T., Bonin, B., Godon, N., Lartigue, J.-E., Jollivet, P., Ayral, A., De Windt, L., and Santarini, G., 2008, "SON68 Nuclear Glass Dissolution Kinetics: Current State of Knowledge and Basis of the New GRAAL Model," *J. Nucl. Mater.*, **380**(1-3), pp. 8-21.

[5] Dubé, M., Doquet, V., Constantinescu, A., George, D., Rémond, Y., and Ahzi, S., 2010, "Modeling of Thermal Shock-Induced Damage in a Borosilicate Glass," *Mech. Mater.*, **42**(9), pp. 863-872.

[6] Kahl, L., Kroebel, R., Storch, W., Holton, L., and Dierks, R., 1991, "Investigation of the Thermal Behaviour of High Level Waste Glass in Steel Canisters During Cooling and at Equilibrium," *Nucl. Eng. Des.*, **130**(1), pp. 77-88.

[7] Daudeville, L., and Carré, H., 1998, "Thermal Tempering Simulation of Glass Plates: Inner and Edge Residual Stresses," *J. Therm. Stresses*, **21**(6), pp. 667-689.

[8] Bernard, F., Daudeville, L., and Gy, R., 2004, "Load Bearing Capacity of Connections in Tempered Glass Structures," *Struct. Eng. Int.*, **14**(2), pp. 107-110.

[9] Kong, J., Kim, J., and Chung, K., 2007, "Residual Stress Analysis With Improved Numerical Methods for Tempered Plate Glasses Based on Structural Relaxation Model," *Met. Mater. Int.*, **13**(1), pp. 67-75.

[10] Prieto, M., Díaz, J., and Egusquiza, E., 2002, "Analysis of the Fluid-Dynamic and Thermal Behaviour of a Tin Bath in Float Glass Manufacturing," *Int. J. Therm. Sci.*, **41**(4), pp. 348-359.

[11] Abaqus FEA (Abaqus/Standard Version 6.9-1), 2009, "Theory Manual Section 2.11.1," SIMULIA/Dassault Systèmes.

[12] Barthélémy, B., 2003, "(Radioactive Waste Combustion-Vitrification Under Arc Plasma: Thermal and Dynamic Modeling) Combustion-Vitrification de déchets radioactifs par plasma d'arc: Modélisation de la thermique et de la dynamique," Ph.D. Thesis, Université de Limoges, France (in French).

[13] Siegel, R., 1998, "Transient Effects of Radiative Transfer in Semitransparent Materials," *Int. J. Eng. Sci.* (Oxford, U. K.), **36**(12-14), pp. 1701-1739.

[14] Riley, W. P., Jr., O'Holleran, T., Johnson, S. G., DiSanto, T., Bateman, K. J., and Goff, K. M., 1999, "Scale Up Issues Involved With the Ceramic Waste Form: Ceramic-Container Interactions and Ceramic Cracking Quantification," *Proceedings of American Ceramic Society 101st Annual Meeting, Symposium: Waste Management Science and Technology in the Ceramic and Nuclear Industries*, Indianapolis, IN.

[15] Mattei, S., 2005, "(Heat Radiation of Opaque Materials) Rayonnement Thermique des Matériaux Opaques," *Techniques de l'ingénieur, Génie Energétique*, BE2(BE8210), pp. BE8210.1-BE8210.24 (in French).

[16] Cereser, W., 1979, "(Pipings: Heat Transfer) Tuyauteries : Transmission de la Chaleur," *Techniques de l'ingénieur, Génie mécanique*, BL4(A820), pp. A820.1-A820.17 (in French).