Modeling of thermal shock-induced damage in a borosilicate glass

Mechanics of Materials 42 (2010) 863–872

Contents lists available at ScienceDirect

Mechanics of Materials journal homepage: www.elsevier.com/locate/mechmat

Modeling of thermal shock-induced damage in a borosilicate glass Martine Dubé a, Véronique Doquet a,*, Andrei Constantinescu a, Daniel George b, Yves Rémond b, Saïd Ahzi b a b

Laboratoire de Mécanique des Solides, UMR CNRS 7649, Ecole Polytechnique, 91128 Palaiseau cedex, France Institut de Mécanique des Fluides et des Solides, Université de Strasbourg, CNRS, 2 rue Boussingault, 67000 Strasbourg, France

a r t i c l e

i n f o

Article history: Received 17 December 2009 Received in revised form 7 July 2010

Keywords: Thermal shock Glass Cracks Continuum damage mechanics

a b s t r a c t Continuum damage mechanics is used to model the damage induced by a thermal shock to the R7T 7 glass, the French borosilicate glass used for nuclear waste vitrification. A finite element model of the thermal shock is developed in which the elastic constitutive equations are coupled with an anisotropic stress-based damage evolution law. The Weibull distributions of strength measured at various temperatures by biaxial flexural tests are used to identify the parameters of this damage evolution law. Vibration tests are conducted to identify the elastic properties of the glass and to determine the effect of thermal shockinduced damage on the glass residual stiffness. The residual stiffness predicted by the damage model agrees with that measured experimentally. The fractured surface in the glass after a thermal shock is estimated, assuming that all the elastic energy associated with the stresses that are higher than a pre-determined threshold is dissipated in the creation of new surfaces. Comparison between the predicted and experimentally-measured fractured surface is performed. The model is shown to capture the saturation of the crack network density for severe thermal shocks; whereas this is not the case if damage is not accounted for in the constitutive equations. Ó 2010 Elsevier Ltd. All rights reserved.

1. Introduction Vitrification of the nuclear wastes consists in solidifying the liquid waste into a glass matrix, and pouring the mixture into stainless steel canisters. In France, the R7T 7 glass is the reference borosilicate glass for this application. During the natural cooling of a waste block, thermal gradients induce triaxial thermal stresses. If these stresses were large enough, cracks, could appear in the glass (Faletti and Ethridge, 1988; Kamizono and Niwa, 1984). This multiple cracking generation is a potential issue as it would increase the area accessible to underground water during the disposal phase of the vitrified wastes, once the carbon steel over-pack has lost its water tightness. A higher accessible surface in turn means an increased radionuclide release rate in the environment (Kamizono and Senoo, 1983; Perez

* Corresponding author. Tel.: +33 1 69 33 57 65; fax: +33 1 69 33 57 06. E-mail address: [email protected] (V. Doquet). 0167-6636/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.mechmat.2010.07.002

and Westsik, 1981). That is why it is important to know the leachable surface of a glass package. Many experimental and analytical or numerical studies of thermal shock damage in brittle materials can be found in the literature (Hasselman, 1969; Tomba et al., 2000; Maensiri and Roberts, 2002; Jeong et al., 2005; Jin and Feng, 2008). These studies generally focus on the residual stiffness and strength of the materials after a thermal shock and on the parameters (microstructure, thermomechanical properties or surface condition) that control the resistance of a material to thermal shocks. The aim of the present work is not to improve the resistance of the glass or to compare its performances with challenger materials, but rather to identify, from simple experiments on small samples of SON 68 glass (inactive R7T 7 type glass), a thermo-mechanical model able to describe thermal shock-induced damage (especially in terms of cracked surface) for a future application in the simulation of the industrial process of melted waste solidification and potential cracking.

864

M. Dubé et al. / Mechanics of Materials 42 (2010) 863–872

To assess the fractured surface, thermo-mechanical finite element simulations of the cooling process will be performed. However, if the thermo-mechanical model assumes an elastic, and not visco-elastic, material behavior and if it is not coupled with any damage evolution law, large overestimations of the stresses can be expected, so that the calculations could not be relied onto predict the damage in the waste canisters. There is a clear need for a model coupled with a damage evolution law. In the first part of the paper, a damage model, based on continuum damage mechanics, is presented and the parameters that need to be identified are highlighted. The second part presents the experimental identification of these parameters. Finally, in the third part, thermomechanical finite element simulations of thermal shock that make use of the damage model are conducted. A method to predict the total fractured surface is proposed and a comparison with experimental results is performed. 2. Glass constitutive equations At temperatures approaching the glass transition temperature of 502 °C, the SON 68 glass exhibits a visco-elastic behavior. Relaxation tests were conducted to characterize this behavior; however, since the present work focuses on low temperatures (T < 200 °C) and short term experiments, the viscosity can be neglected and the visco-elastic behavior of the glass will be reported in a subsequent paper. Since the pioneering work of Kachanov (1980), Lemaitre (1971) and Chaboche (1974), continuum damage mechanics was extensively used to describe the stiffness degradation of a material undergoing damage, especially for brittle materials such as concrete and rocks. Existing isotropic damage models were modified to take into account the damage-induced anisotropy in brittle materials and the tension/compression asymmetry due to the crack closure effect (Desmorat et al., 2007). A simple anisotropic damage model was developed by Sun and Khaleel (2004) to investigate the damage caused by an indenter on a soda-lime glass. The model predicted the cracking damage pattern and zone size adequately and was in agreement with experimental results (Ismail et al., 2008; Zhao et al., 2006). This model is simple enough for the modeling of a complex vitrification process that includes phase transition from liquid to visco-elastic solid and finally elastic brittle solid at low temperatures. It was shown to be well adapted for glass and, therefore was selected as a proper damage model for this study. In modeling the damage caused to the SON 68 glass, the initial undamaged material was considered to be isotropic

8 > > > > > > > > <

r1 9 > > > > r2 > > > > r3 = > r12 > > > > > > > > > > > r > > 23 > > : ; r31

and stress-free. The stiffness degradation arising from the damage accumulation was simulated by adding a damage tensor into the glass constitutive equation (Sun and Khaleel, 2004):

n

o 

rij ¼ K eijkl ðTÞ þ K dijkl ðTÞ  ekl  ethkl



ð1Þ

where rij and ekl are the stress and strain tensor components, respectively, and K eijkl ðTÞ and K dijkl ðTÞ are the temperature-dependent fourth order stiffness tensors representing the undamaged isotropic material and the added influence of damage, respectively. eth kl represent the thermal strain tensor components. The components of the stiffness tensors are given by:

K eijkl ¼ kðTÞdij dkl þ lðTÞðdik djl þ dil dkj Þ

ð2Þ

K dijkl ¼ C 1 ðTÞðdij Dkl þ dkl Dij Þ þ C 2 ðTÞðdjk Dil þ dil Djk Þ

ð3Þ

where k(T) and l(T) are Lame’s constants, Dij are the damage parameters and C1(T) and C2(T) are two constants of the material. Since we are interested in the thermomechanical behavior of the glass, Lame’s constants and C1 and C2 are considered temperature-dependent. The terms of the damage tensor, Dij, are function of the stress state and take values between zero and one, zero meaning that the material is virgin (undamaged) and one meaning that the material is fully damaged, i.e., can no longer sustain stresses. Though the original damage model of Sun and Khaleel (2004) was meant to take into account the effect of the shear stresses on the damage evolution, only the diagonal terms of the damage tensor, i.e., the damage parameters due to tensile principal stresses (opening mode), were accounted for in this study. Their values are assumed to follow a linear evolution law such that:

Dii ¼

8 > <0

ri rth ðTÞ > rc ðTÞrth ðTÞ

:

1

ri 6 rth rth < ri < rc ri P rc

ð4Þ

where i = 1, 2, 3, ri is the ith principal tensile stress, rth is the temperature-dependent threshold stress under which no damage occurs and rc is the temperature-dependent critical stress above which the material is fully damaged. Note that in the application considered below (thin disk specimens subjected to a thermal shock) loading is proportional (rrr = rhh), so that the principal axes do not change with time and correspond to the radial, tangential and normal directions. The damage tensor will thus be naturally defined in this coordinate system. From Eqs. (1)–(3) and neglecting the non-diagonal terms of the damage tensor, the glass constitutive equation is:

8 k þ 2l þ 2D11 ðC 1 þ C 2 Þ > > > > > k þ C 1 ðD11 þ D22 Þ k þ 2l þ 2D22 ðC 1 þ C 2 Þ SYM > > > < k þ C ðD þ D Þ k þ C 1 ðD22 þ D33 Þ k þ 2l þ 2D33 ðC 1 þ C 2 Þ 1 11 33 ¼ > 0 0 0 l > > > > > 0 0 0 0 > > : 0 0 0 0

98 > > > > > > >> > > > > > > > < =>

9

e1 > > > > e2 > > > > e3 = > > 0 > e12 > > > > > > >> > > > > > > l > e > > > 23 > > > ;: ; 0 l e31

865

M. Dubé et al. / Mechanics of Materials 42 (2010) 863–872

The values of C1 and C2 are calculated by stating that the stresses r11 and r22 drop to zero when D11 = D22 = 1.0, for equi-biaxial tension in directions 1 and 2 i.e., for D33 = 0. For such an equi-biaxial tension test, e11 = e22 = e:

r11 ¼ r22 ¼ 0 ¼ f2ðk þ lÞ þ 2ð2C 1 þ C 2 Þge þ fk þ C 1 ge33 ð5Þ

r33 ¼ 0 ¼ f2ðk þ C 1 Þge þ fk þ 2lge33

ð6Þ

Eqs. (5) and (6) must be equivalent. Following this logic, we find C1 = 2l and C2 = 3l. This choice of C1 and C2 is different from that of Sun and Khaleel (2004) for soda-lime glass. They chose the two parameters (C1 = l and C2 = 1.5l) so that the stress r11 drops to zero when D11 = 1.0, for uniaxial tension in direction 1. In the application presented here (glass disks subjected to thermal shocks) loading is equi-biaxial and identification of the parameters in equi-biaxial loading was found to allow a better description of thermal shock damage. 3. Identification of the thermo-mechanical model parameters 3.1. Thermal properties The thermal properties of the SON 68 glass (thermal conductivity, k, thermal expansion coefficient, aa, and specific heat capacity, cp) must be identified in order to perform the thermo-mechanical finite element simulation of the thermal shock. The heat capacity was determined by differential scanning calorimetry and the diffusivity, j, by the laser flash method in which one face of a cylindrical specimen is suddenly exposed to a sharp temperature rise and the temperature of the other face is monitored. The rate at which the temperature becomes uniform within the specimen volume is function of the specimen geometry and material diffusivity. The thermal expansion coefficient was determined by dilatometry and the thermal conductivity was deduced from to the following equation:

k ¼ qcp j

ð7Þ

where q is the density. The measured values are reported in Table 1. 3.2. Mechanical properties The mechanical properties that are needed are Lame’s constants that are related to Young’s modulus (E(T)) and Poisson’s ratio (t(T)) by the classical formulae:

kðTÞ ¼

EðTÞtðTÞ ; ð1 þ tðTÞÞð1  2tðTÞÞ

lðTÞ ¼

EðTÞ 2ð1 þ tðTÞÞ

ð8Þ

and, for the damage model, the threshold and critical stresses (rth(T), rc(T)). Those parameters were identified using biaxial flexural tests and vibration tests conducted at temperatures varying between 20 and 506 °C. 3.2.1. Biaxial flexural tests The biaxial flexural tests were carried out according to the ASTM C 1499 – 05 standard test method, on 2.04 mmthick disk specimens with a diameter of 40.0 mm. The face of the specimens subjected to tensile stresses was polished to a roughness of 1.0 lm. The disks were placed between a 30 mm diameter support ring and a 10 mm diameter concentric loading ring. Displacement-controlled monotonic loading was applied on the loading ring until a clear force drop was observed, due to the failure of the specimens. The tests were conducted at temperatures of 20, 200 and 506 °C. A resistive oven, in which the specimens were heated and kept at the desired temperature for 60 min, was used for the high temperatures tests. No less than 17–20 specimens were tested at each temperature in order to have a minimum of 13 valid tests. Tests were considered valid when the fracture origin was located inside the diameter of the smaller loading ring, i.e., at a radial position r < 5 mm. Examples of valid and non-valid tests are presented on Fig. 1. The force–displacement data were recorded for determining the glass elastic properties and its biaxial flexural strength. The broken specimens were then observed using optical microscopy to identify the likely failure origin. The cumulative strength distributions for T = 20, 200 and 506 °C are shown in Fig. 2. The biaxial flexural strengths were calculated from (ASTM C1499 – 05):

rf ¼

"

3F 2

2p h

ð1  tÞ

2

2

ds  d1 2

2dd

þ ð1 þ tÞ ln

ds dl

# ð9Þ

where rf is the biaxial flexural strength, F is the load at failure, h is the plate thickness and ds, dl and dd are the diameters of the support ring, loading ring and disk specimen, respectively. The Weibull modulus, which is a measure of the material microstructural variability, was found to have a value of 6.0, for all investigated temperatures. This is in agreement with the typical Weibull modulus of other brittle materials such as ceramics, which generally varies between 5.0 and 10.0 (Ovri, 2000; Forquin et al., 2003; Danzer et al., 2008). The threshold stress for damage was determined from the measured strength distributions as the stress above which damage can possibly be created. The lowest strength obtained during the tests, i.e.,

Table 1 SON 68 thermal and mechanical properties. Temperature (°C)

Young’s modulus (GPa) (%)

Poisson’s ratio

Coefficient of thermal expansion (°C1) (%)

Heat capacity (J/kg/°C) (%)

Thermal conductivity (W/m/°C) (%)

20 100 200 300 400 500

92.8 ± 1.5 89.3 ± 2.9 88.3 ± 2.2 87.2 ± 2.4 86.1 ± 3.0 39.8 ± 2.3

0.19–0.25 0.19–0.30 0.16–0.26 0.16–0.27 0.18–0.30 0.15–0.25

8.3E6 ± 5.0 8.3E6 ± 5.0 8.3E6 ± 5.0 8.3E6 ± 5.0 9.1E6 ± 5.0 9.1E6 ± 5.0

800 ± 7.0 800 ± 7.0 875 ± 7.0 920 ± 7.0 950 ± 7.0 1000 ± 7.0

1.03 ± 12 1.03 ± 12 1.08 ± 12 1.11 ± 12 1.11 ± 12 1.15 ± 12

866

M. Dubé et al. / Mechanics of Materials 42 (2010) 863–872

Fig. 1. Examples of valid (a) and non-valid (b) tests in biaxial flexural tests, based on the likely failure origin.

Optical microscopy revealed that near-surface bubbles in the glass acted as initiation sites. Those bubbles had diameters as large as 1.0 mm and were uniformly-distributed over the specimens lot. They were formed during the glass manufacturing process and resulted in stress concentration which was responsible for the formation of cracks during the tests. The force–displacement data recorded during the tests were used to obtain the elastic properties of the glass. The Young modulus for a circular plate under biaxial flexural loading can be estimated from (ASTM C 1499 – 05, 2005):

1 T = 20ºC

Failure Probability

0.8

T = 200ºC T = 506ºC

0.6

0.4 Critical Stress (σc)

Threshold Stress (σ th) 0.2

0 40

2

50

60

70

80

Biaxial Flexural Strength (MPa) Fig. 2. Biaxial flexural strength distributions of SON 68 glass at various temperatures.

41 MPa, was selected as the threshold stress. On the other hand, the critical stress is the stress at which the material is fully damaged. The damage mechanics model is a phenomenological model in which the various parameters are adjusted to represent at best an experimental behavior. The choice of the mean value of the biaxial bending strength distribution, 64 MPa, (in agreement with previous studies from Sun and Khaleel (2004) or Ismail et al. (2008)) was found to be more effective than that of the maximum value, even though it might seem less intuitive. Besides the damage mechanics model is deterministic and cannot capture the scatter inherent to brittle fracture. Had the maximum value of the distribution been chosen for the critical stress, the fact that most specimens fail below this stress would not be predicted. Choosing the mean value, the ‘‘average behavior” is captured. Fig. 2 shows that the strength of the glass does not vary within the investigated temperature range. Therefore, the threshold and critical stresses were considered to be the same, for all temperatures.

Effi

3Fð1  t2 Þdl 8pdh

3

2

ds

2

dl

" 1þ

# 2 2 ð1  tÞðds  dl Þ 2

2ð1 þ tÞdd

 ! ds  1 þ ln dl ð10Þ

where d is the disk deflection. A first relationship between Young’s modulus and Poisson’s ratio is obtained with this equation but a second test is needed to identify the two properties separately. 3.2.2. Impulse excitation of vibration Specimens made of linear elastic materials possess modal frequencies that are determined by their Young modulus, Poisson ratio, mass and geometry. The elastic properties of a material can therefore be deduced from the measured modal frequencies (Fletcher and Rossing, 1998). In this experiment, SON 68 glass disks with diameter of 40 mm and thickness of 2.04 mm were used. The specimens were supported by four rubber bands that were bonded to the sides of the specimens and attached to a metallic structure. The specimens were excited using an impact hammer and their acoustic response was captured by a microphone that was connected to an amplifier and an oscillator (Fig. 3). A Fourier transform was applied on the signal in order to obtain the first and second modal frequencies of the specimens. According to (Fletcher and Rossing, 1998), the

M. Dubé et al. / Mechanics of Materials 42 (2010) 863–872

Uz

Vibrating Mode I

Microphone

Vibrations Test Schematic

4. Predicting thermal shock damage 4.1. Thermal shock experiments

Impact

Vibrating Mode II

Fig. 3. Schematic of the vibration test and iso-contours of normal displacements corresponding to the first two vibrating modes of a freely supported disk, as computed by F.E.M.

equation relating the modal frequency of the first vibrating mode to the elastic properties of a thin disk is:

f0 ¼

accounts for the experimental error. A large range of possible Poisson’s ratio values exists. By cross-referencing the values obtained here with the measurements performed during the glass qualification tests, a constant value of 0.24 was taken as the Poisson ratio in all calculations made throughout this study.

Amplifier

Oscillator

Disk Specimen

867

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi E h b qð1  t2 Þ R2

ð11Þ

where R is the disk radius and b is a coefficient that depends on the boundary conditions. For free edges boundary conditions, b takes a value of 0.2413 (Fletcher and Rossing, 1998). The modal frequency corresponding to the second vibrating mode (see Fig. 3 for first and second mode shapes) is 1.73fo. The described experiment was performed at temperatures varying between room temperature and 506 °C thanks to a radiant oven in which the experimental set-up could be adapted. Eqs. (10) and (11) provide two independent relations from which the Young modulus and Poisson ratio can be deduced. The obtained values are summarized in Table 1. The range of presented values at each temperature

Thermal shock experiments were carried out on SON 68 glass disk specimens of thickness 2.04 mm and diameter 40 mm. Both faces of the specimens were polished to 1.0 lm prior to the thermal shock in order to make sure that cracks would not initiate on machining scratches. The disks were heated to a pre-determined temperature in a radiant oven and kept at this temperature for 60 min in order to insure a uniform temperature within their volume. The disks were then quenched into water kept at room temperature. Following quenching and cooling to room temperature, the specimens were observed by optical microscopy in order to assess the thermal shockinduced damage. Approximately 20 micrographs were taken for each specimen. The observations revealed many curved and branched cracks (Fig. 4). This curvature can be explained by (i) the equi-biaxial tensile stress field which arises during the thermal shock and (ii) the mode-mixity induced by the interactions between cracks (Kane and Doquet, 2006). The cumulated developed length of cracks on both the upper and lower surfaces, denoted by lc, was determined from these micrographs, taking the curvature into account. In order to measure the depth of the surface cracks in the thickness direction, several specimens were cut along different diameters. The cracks were curved in that plane as well, due to mutual interactions which increase with the severity of the thermal shock and the density of the crack network (Fig. 5). The average developed ‘‘depth” of cracks on those transverse sections, denoted by dc, was

Fig. 4. Cracks formed at the surface of a specimen after a thermal shock of amplitude DT = 157 °C. The large black dots are bubbles and the white spots are second phase Ru and Pd particles.

868

M. Dubé et al. / Mechanics of Materials 42 (2010) 863–872

Cracks

Damaged layer

z Undamaged layer

500 μm

during the thermal shocks is 180 °C and stress relaxation tests performed in biaxial bending at temperatures up to 500 °C, showed that viscous effects below 180 °C can be neglected for such brief tests. A convection heat transfer coefficient h = 10,000 W/m2/ °C was selected, in accordance with literature for convection in water (Incropera and DeWitt, 1996; Aydıner and Üstündag, 2005). The time scale of the order of 0.05 s permitted to assume that the specimen passes directly from the oven to a complete immersion into water. The timedependent temperature field was defined by the classical heat transfer equation:

Damaged layer

qcp ðTÞDT

Fig. 5. Transverse section of a disk showing thermal shock-induced cracks (DT = 157 °C).

measured. The cumulated fractured surface was finally approximated as:

Stot ¼ lc dc

ð12Þ

As seen on Fig. 5, following the thermal shock, the glass specimens can be thought of as a composite material in which two anisotropic damaged layers exist at the surfaces with an undamaged, and therefore isotropic, core.

dT ¼ r½kðTÞrT þ Q_ dt

ð13Þ

in which the internal heat source Q_ is zero. A mechanical computation was then performed in which the loading was defined by the preceding temperature field. The thermal stresses were computed based on the anisotropic damage evolution law presented in Section 2. The constitutive law was directly implemented in Cast3m. We considered a forward integration scheme with respect to the damage evolution. Since the virgin material is assumed to be free of any defect, the damage parameters were set to zero at the beginning of the analysis and due to the irreversibility of the damage process, Dii was considered a monotonically increasing function of time:

Dii ¼ maxðDnii ; Dn1 ii Þ 4.2. Finite element model of the thermal shock An axial-symmetric thermo-mechanical finite element model of the thermal shock was developed using the Cast3m finite element program (www-cast3m.cea.fr). The cross-section of the disks was meshed with eight-node quadratic elements that were used to represent both the temperature and displacement fields. Eighteen and 40 elements were used across the thickness and radius of the disk, respectively. Convection boundary conditions were applied on the free surfaces of the geometry (Fig. 6). The thermal calculations were based on the thermal properties of the glass (Table 1) and density q = 2845 kg/m3. The elastic properties were assumed to vary linearly with temperature between the measured values. It takes around 0.2 s during the thermal shock for the maximum stresses to be developed at the surfaces of the specimens and 10 s in the biaxial flexural tests. This is a significant difference. However, the maximum temperature reached

Axial-symmetric axis

ð14Þ

where n and (n  1) are the nth and (n  1)th time step of the analysis, respectively. The analysis was conducted incrementally, i.e., the stresses were computed at each step as:

n

o n

rnij ¼ K eijkl þ K d;ðn1Þ  enkl  eth;n ijkl kl

o

ð15Þ

where the stiffness tensor accounting for damage is calculated based on the stresses of the previous time increment, and the thermal strains are subtracted from the total strains. Convergence in terms of mesh refinement and time increment was verified. 4.3. Results 4.3.1. Thermal stresses and damage parameters Upon cooling, tensile stresses are generated on the outer surfaces of the disk specimens. The state of stress is equi-biaxial with rrr = rhh and rzz = 0. As mentioned in

Twater = 23ºC

Uniform Temperature

Applied convection on free surfaces (h = 10,000 W/m 2/°C) Fig. 6. Finite element mesh of the disk specimen.

M. Dubé et al. / Mechanics of Materials 42 (2010) 863–872

Section 2, rrr and rhh correspond to the principal stresses responsible for damage. The evolution of the stresses with time is depicted in Fig. 7 for a point located at mid-radius, on the outer surface of the specimen, for a thermal shock amplitude DT = 157 °C. The curve obtained using a simple elastic thermomechanical model, without any damage taken into account, is also plotted for comparison. When accounting for damage, the stresses are cut off at the critical stress rc = 64 MPa, whereas they reach 88 MPa, which is above the maximum of the strength distribution plotted in Fig. 2, if damage is not taken into account. The peak stress is reached at 0.11 s. Fig. 8 shows the profile of rrr in the thickness direction, at t = 0.11 s, at the same radial position. Only the upper half of the disk is presented, as the profile is symmetrical with respect to z. As expected, rrr is positive on the outer surfaces of the disk and becomes negative at a certain depth. Fig. 9 shows the iso-contours of the damage parameter D11 after a thermal shock amplitude DT = 157 °C. Since rrr = rhh, D22 presents the same distribution as D11. As in the experiments, two superficial damaged layers appear. Finally, Fig. 10 compares the profile of the damage parameter D11 along the z-axis of the specimen after ther-

100

σrr or σθθ (MPa)

80

σc

Thermo-mechanical model without damage

60 40

Point considered

z r

Thermo-mechanical model coupled with damage law

20 Axial-symmetric axis 0

0

0.05

0.1

0.15

0.2

0.25

Time (s) Fig. 7. Evolution of radial and tangential stresses with time during a thermal shock of amplitude DT = 157 °C for a point located at the surface of the specimen.

70 60 σc

50

σ th

σrr (MPa)

40 30 20 10 0 0

0.2

0.4

0.6

0.8

1

-10 -20 -30

Position along the z-axis (mm) Fig. 8. rrr profile across the specimen depth, at t = 0.11 s, for a thermal shock amplitude DT = 157 °C. Z = 0 corresponds to the specimen midthickness.

869

mal shocks of various amplitudes. The larger the thermal shock amplitude, the deeper the damage penetrates inside the specimens. 4.3.2. Fractured surface One of the objectives of the study is to estimate the fractured surface in a glass specimen subjected to a thermal shock. To do this, a conservative assumption was made stating that all the elastic energy associated with the stresses that are higher than the threshold stress is dissipated in the creation of new surfaces. The peak value of the elastic energy over the time was taken as:

Maxtime

1 2

Z V rPrth

rij eij dV  2cs Stot

ð16Þ

in which ½ rijeij (J/m3) is the elastic energy density, V rPrth (m3) is the volume over which the threshold stress is reached by at least one principal stress, Stot (m2) is the total fractured surface and cs (J/m2) is the surface energy which is determined by:

cs ¼

K 2Ic 2E

ð17Þ

p in which KIc is the material toughness of 0.85 MPa m, which is assumed to be constant over the investigated temperature range (20–180 °C), based on measurements done at 20 and 424 °C that yielded similar KIc, allowance made for experimental scatter. The factor two in Eq. (16) accounts for the two faces created for each crack. The measured and predicted fractured surfaces are compared in Fig. 11 for various thermal shock amplitudes. The fractured surface predicted by a standard thermo-mechanical model that does not account for damage is also shown as a reference. At low thermal shock amplitudes, a fairly accurate prediction is obtained for both models. By contrast, at high thermal shock amplitudes, i.e., DT P 157 °C, the model that does not account for damage does not predict the saturation of the fractured surface. This result was expected as the stresses, and therefore the computed elastic energy, should keep increasing with the thermal shock amplitude. The model accounting for damage succeeds in predicting the saturation of the fractured area at DT P 157 °C. Though the model somewhat overestimates the fractured surface, it provides an upper bound of the fractured area. It is assumed that the excess elastic energy, i.e., the part that is not dissipated into the creation of new surfaces, is dissipated into other forms such as vibrations, acoustic energy, etc. Some error on the convection heat transfer coefficient, h, might also partly explain the discrepancy. 4.3.3. Residual stiffness The damage model states that the stiffness of the glass decreases during the thermal shock when the tensile principal thermal stresses are higher than the threshold stress. This stiffness reduction was assessed by subjecting predamaged disk specimens to the vibration test described in Section 3.2.2. All the tests were done at room temperature. As stated before, the damage model includes an anisotropic damage evolution law. Due to the damageinduced anisotropy and heterogeneity of the specimens

870

M. Dubé et al. / Mechanics of Materials 42 (2010) 863–872

Fig. 9. Iso-contours of damage parameter D11 after a thermal shock of amplitude DT = 157 °C.

1

ΔT = 157ºC ΔT = 137ºC ΔT = 117ºC ΔT = 97ºC ΔT = 77ºC ΔT = 57ºC

0.8

D11

0.6

0.4

0.2

0 0.0

0.2

0.4

0.6

0.8

1.0

Position along the z-axis (mm) Fig. 10. Profiles of damage parameter D11 along the z-axis computed after thermal shocks of amplitudes DT. Z = 0 corresponds to the specimen midthickness.

after the thermal shocks, their elasticity cannot be expressed solely in terms of Young’s modulus and Poisson’s ratio. The ratio of modal frequencies measured after and before the thermal shock provides a global and more pertinent evaluation of the stiffness reduction. Following the thermal shock simulation, an axialsymmetric finite element model of the vibration test was developed in which the final stiffness tensor obtained at the last time increment of the thermal shock simulation, N i.e., K eijkl þ K dijkl , was implemented. The finite element model used the same mesh as that used for the thermal shock simulation (Fig. 6). The modal frequency of the damaged material was obtained from this simulation and comparison was performed between the stiffness reduction predicted by finite element and that obtained experimentally (Fig. 12). As can be seen, the predicted stiffness degradation follows the same trend as that measured experimentally with the predicted values falling within the experimental error bars.

Experimental Thermo-mechanical model without damage Thermo-mechanical model coupled with damage law Experimental Thermo-mechanical model coupled with damage law

Fig. 11. Predicted and measured fractured surfaces in glass specimens subjected to thermal shocks of amplitudes DT.

Fig. 12. Residual stiffness after thermal shocks of amplitudes DT, illustrated by the ratio of the natural frequencies of the damaged and undamaged glass.

M. Dubé et al. / Mechanics of Materials 42 (2010) 863–872

Biaxial Flexural Strength (MPa)

4.3.4. Residual strength Following the thermal shocks, specimens were tested to failure using the biaxial flexural test described in Section 3.2.1 in order to assess the effect of the thermal shock-induced damage on their residual strength. The results (Fig. 13) are in agreement with Hasselman’s conclusions (1969) that a sharp strength reduction is obtained at a critical thermal shock amplitude DTc. In our case, the critical thermal shock amplitude for which the strength reduction is of the order of 75%, is DTc = 77 °C. The residual strength then remains constant at larger thermal shock amplitudes. This behavior is typical of brittle materials and was observed before for ceramics (Maensiri and Roberts, 2002; Jin and Feng 2008). Observations of the broken specimens revealed a different failure behavior from that observed for the undamaged specimens. As shown in Fig. 14, the cracks followed an irregular path due to further extension and coalescence of several pre-existing cracks generated by the thermal shock. This contrasts with the straighter crack paths observed for the undamaged material (Fig. 1) which failed by extension of a single crack that initiated near the center of the disk and then branched and propagated along more or less radial directions.

80

60

40

20

871

5. Conclusions The thermo-mechanical behavior of a borosilicate glass used for nuclear waste vitrification was investigated by experiments and numerical simulations. A continuum damage mechanics model was identified and used to predict the damage induced by a thermal shock. Comparison between predicted and measured damage based on the residual stiffness of the glass was performed, and a satisfactory agreement was obtained. An upper bound of the total fractured surface induced during a thermal shock was estimated, assuming that the elastic energy associated to the transient stresses that are higher than a threshold value is dissipated in the creation of new surfaces. The damage model captured the saturation of the crack network density for severe thermal shocks. The residual strength of the glass after a thermal shock was measured and a critical thermal shock amplitude of DT = 77 °C was determined. Moreover, the cracking pattern of the thermally-shocked disks was different from that of the undamaged specimens. A finite element simulation using the presented damage model cannot be straightforwardly applied to the biaxial flexural test to predict the residual strength. Indeed, numerical difficulties (mesh dependency, non-convergence close to final fracture, etc.) that are quite common with such a local continuum damage model would be encountered (Lemaitre et al., 2009). A non-local formulation of the model should be developed in the future, allowing for these difficulties to be overcome and the residual strength to be predicted. Acknowledgements

0 0

50

100

150

200

ΔT (ºC) Fig. 13. Residual biaxial flexural strength after thermal shocks of amplitudes DT.

The authors are indebted to Frédéric Bouyer from the CEA for numerous helpful discussions as well as providing the SON 68 glass and Xavier Boutillon from the Laboratoire de Mécanique des Solides of the Ecole Polytechnique for valuable advices regarding the stiffness measurements. The CEA, Andra and AREVA are also gratefully acknowledged for supporting this study. References

Fig. 14. Aspect of a thermally-shocked specimen after being tested to failure by biaxial flexural tests.

Aydıner, C.C., Üstündag, E., 2005. Residual stresses in a bulk metallic glass cylinder induced by thermal tempering. Mechanics of Materials 37 (1), 201–212. ASTM C 1499 – 05, 2005. Standard test method for monotonic equibiaxial flexural strength of advanced ceramics at ambiant temperature. ASTM International. Chaboche, J.-L., 1974. Une Loi Différentielle d’Endommagement de Fatigue avec Cumulation Non-Linéaire. Revue Française de Mécanique, 50–51 (in French). Danzer, R., Lube, T., Supancic, P., Damani, R., 2008. Fracture of ceramics. Advanced Engineering Materials 10 (4), 275–298. Desmorat, R., Gatuingt, F., Ragueneau, F., 2007. Nonlocal anisotropic damage model and related computational aspects for quasi-brittle materials. Engineering Fracture Mechanics 74, 1539–1560. Faletti, D.W., Ethridge, L.J., 1988. A method for predicting cracking in waste glass canisters. Nuclear and Chemical Waste Management 8, 123–133. Fletcher, N.H., Rossing, T.D., 1998. The Physics of Musical Instruments. Springer Science + Business Media Inc.

872

M. Dubé et al. / Mechanics of Materials 42 (2010) 863–872

Forquin, P., Denoual, C., Cottenot, C.E., Hild, F., 2003. Experiments and modeling of the compressive behavior of two SiC ceramics. Mechanics of Materials 35 (10), 987–1002. Hasselman, D.P.H., 1969. Unified theory of thermal shock fracture initiation and crack propagation in brittle ceramics. Journal of American Ceramic Society 52 (11), 600–604. Incropera, F.P., DeWitt, D.P., 1996. Introduction to Heat Transfer. John Wiley. Ismail, J., Zaïri, F., Naït-Abdelaziz, M., Azari, Z., 2008. Computational modeling of static indentation-induced damage in glass. Computational Materials Science 42, 407–415. Jeong, J., Adib, H., Pluvinage, G., 2005. Proposal of new damage model for thermal shock based on dynamic fracture on the brittle materials. Journal of Non-Crystalline Solids 351, 2065–2075. Jin, Z.-H., Feng, Y.Z., 2008. Effects of multiple cracking on the residual strength behavior of thermally shocked functionally graded ceramics. International Journal of Solids and Structures 45, 5973–5986. Kachanov, M., 1980. Continuum model of medium with cracks. Journal of the Engineering Mechanics Division, ASCE 106 (5), 1039–1051. Kamizono, H., Niwa, K., 1984. An estimation of the thermal shock resistance of simulated nuclear waste glass under water quenching conditions. Journal of Materials Science Letters 3, 588–590. Kamizono, H., Senoo, M., 1983. Thermal shock resistance of a simulated high-level waste glass. Nuclear and Chemical Waste Management 4, 329–333.

Kane, A., Doquet, V., 2006. Surface crack and cracks networks in biaxial fatigue. Engineering Fracture Mechanics 73, 233–251. Lemaitre, J., 1971. Evaluation of Dissipation and Damage in Metals Submitted to Dynamic Loading, ICM 1. Kyoto, Japan. Lemaitre, J., Chaboche, J.-L., Benallal, A., Desmorat, R., 2009. Mécanique des Matériaux Solides, 3rd ed. Dunod, Paris. Maensiri, S., Roberts, S.G., 2002. Thermal shock of ground and polished alumina and Al2O3/SiC nanocomposites. Journal of European Ceramic Society 22, 2945–2956. Ovri, J.E.O., 2000. A parametric study of the biaxial strength test for brittle materials. Materials Chemistry and Physics 66, 1–5. Perez Jr., J.M., Westsik Jr., J.H., 1981. Effects of cracks on glass leaching. Nuclear and Chemical Waste Management 2, 165–168. Sun, X., Khaleel, M.A., 2004. Modeling of glass fracture damage using continuum damage mechanics–static spherical indentation. International Journal of Damage Mechanics 13, 263–284. Tomba, M., Cavalieri, A.L., 2000. Ceramic surface finish: its influence on the heat transfer in thermal shock tests. Materials Letters 42, 240– 245. Zhao, S., Dharani, L.R., Chai, L., Barbat, S., 2006. Analysis of damage in laminated automotive glazing subjected to simulated head impact. Engineering Failure Analysis 13, 582–597.