Co-development of experimental and simulation methods for the laser flash heating and melting technique: The thermoelastic effects of UO2

International Journal of Thermal Sciences 132 (2018) 174–185

Contents lists available at ScienceDirect

International Journal of Thermal Sciences journal homepage: www.elsevier.com/locate/ijts

Co-development of experimental and simulation methods for the laser flash heating and melting technique: The thermoelastic effects of UO2

T

M.J. Wellanda,b,∗, R. Böhlera, L. Vlahovica, K. Boboridisa, D. Manaraa,∗∗ a b

European Commission, Joint Research Centre - Karlsruhe, Postfach 2340, Karlsruhe, 76125, Germany Canadian Nuclear Laboratories, Chalk River, ON K0J 1J0, Canada

A R T I C LE I N FO

A B S T R A C T

Keywords: UO2 Laser flash melting Laser flash Thermoelasticity

The thermoelastic effects on the behaviour of subsecond laser heating experiments on UO2 to temperatures approaching 3000 K are explored using a thermoelastic model coupling heat transport and large deformation hyperelasticity, with consideration of thermoelastic effects. A series of steady state experiments at temperatures up to ∼2600 K are conducted and used to calibrate heat loss models. A series of subsecond transient laser flash melting experiments are then conducted and simulated. Simulation results reproduce the observations well without tuning of the transient model. The importance of considering thermoelastic effects in modelling the behaviour of solid materials with temperature and stress gradients typical of laser flash heating experimentation is discussed.

1. Introduction Subsecond laser heating coupled with fast optical pyrometry provides an experimental method of observing the material behaviour of refractory materials at temperatures up to, and in excess of, 3000 K while limiting chemical modification such as noncongruent vaporization and reaction with the sample crucible and atmosphere. In particular, this technique is being used to observe the high temperature properties up to and beyond the melting temperatures of some nuclear fuels and has led to significant changes in the phase diagram of materials such as mixed oxide (MOX) fuels [1–5]. In subsecond laser flash experiments, samples of oxide, carbide, or nitride fuels are subjected to a laser pulse which heats the surface of the sample to temperatures potentially beyond 3000 K in times ranging from tens to hundreds of milliseconds. Cooling is then allowed to occur naturally via radiation and convective/conductive exchange with the apparatus and atmosphere. The surface temperature is observed with fast optical pyrometry during heating and cooling cycles. Material properties such as the melting temperature [1,2,6], thermal conductivity [7,8] and heat capacity [9,10] may be determined from the resulting thermograms, sometimes through the aid of simulation. Due to the short high-temperature dwelling time and localized heating area, the subsecond laser-flash heating technique is able to circumvent many of the difficulties associated with very high temperatures such as reaction with experimental apparatus or buffer gas

∗

and non-congruent vaporization leading to compositional shift [3]. However, the large temperature gradients incurred can produce large stresses via uneven thermal expansion, which can result in mechanical failure leading to sample destruction or non-reproducible heating curves due to cracks acting as thermal barriers. Additionally, the heating of material in the laser spot which is constrained by the cooler surroundings may generate sufficient stresses to encounter thermoelastic effects including modified heat capacities, melting temperature, and enthalpy of fusion. These effects complicate simulation of the experiment, which otherwise aid in analysis. Models for this experiment have previously been developed to help elucidate the important phenomena leading to the observed thermogram [11,12]. For example, such models capture non-congruent phasechange for samples of binary systems. The transient phase transition, and possibly thermomigration, leads to compositional redistribution within the sample during the pulse and thus it is not immediately clear what composition is being observed to solidify by the pyrometer. While informative, previous simulations left room for improvement in the prediction of the thermogram and freezing plateau length and trend. The current work describes the co-development of the laser flash heating method to increase the quality of the simulation predictions. A revised experimental procedure is investigated, and the significance of thermoelastic effects is explored. The current model is derived in a selfconsistent fashion from the thermodynamic description including hyperelasticity. The large range in temperatures during the experiment,

Corresponding author. European Commission, Joint Research Centre - Karlsruhe, Postfach 2340, Karlsruhe 76125, Germany. Corresponding author. E-mail addresses: [email protected] (M.J. Welland), [email protected] (D. Manara).

∗∗

https://doi.org/10.1016/j.ijthermalsci.2018.05.035 Received 26 February 2018; Received in revised form 16 May 2018; Accepted 22 May 2018 1290-0729/ © 2018 The Authors. Published by Elsevier Masson SAS. This is an open access article under the CC BY license (http://creativecommons.org/licenses/BY/4.0/).

International Journal of Thermal Sciences 132 (2018) 174–185

M.J. Welland et al.

approximately 1400–3500 K, precludes an infinitesimal strain approximation and so the model is treated with a large deformation formulation, and the heat balance calculated on the deformed domain. Derivation from the thermodynamic basis ensures thermoelastic effects, as described partially by the Grüneisen parameter, are captured while treating coupled heat transport and elasticity. The current analysis is limited to the solid regime as a matter of simplicity. It is the long-term goal of this experimental and simulation development to investigate the non-congruent melting/freezing transition, however this requires understanding of the state of the material before melting. Therefore, it was decided to first reproduce the solid heating characteristics reliably. The extension of the work to melting will be investigated separately. This work is divided as follows: In Section 2, the model including the governing equations, hyperelastic thermodynamic model, and material properties is derived. In Section 3 an experimental procedure involving preheating in order to increase the reproducibility of the observations is discussed. In Section 4 a set of steady state experiments to calibrate heat loss terms as a function of temperature is discussed. In Section 5, sets of transient experiments compared to simulations using the heat loss terms from Section 4 are shown. Finally, the importance of the results and interpretation with regards to the importance of thermoelastic effects, and some applications of the present investigation to the laser-heating analysis of particularly challenging materials that has been ongoing while this document was prepared is discussed in Section 6.

Fig. 1. Photograph of a sample mounted in the apparatus, with the four screws suspending it in the centre of the graphite ring. This sample was shot 30 times for a total of approximately 10 s, with significant time in the molten state which results in the clearly visible refrozen area in the centre of the sample. The sample was 8.35 mm in diameter. The laser beam was aimed at the centre of the front face with a nominal spot size of 8 mm. Vapour has clearly been transported vertically by convective currents and condensed in a yellow material on the screws and the graphite ring. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

2.1. Governing equations 2. Model development

Two coordinate systems are defined for this problem. Coordinates ⎯→ ⎯ X are defined on the undeformed reference domain corresponding to the original configuration of the material at a reference temperature and pressure. The coordinates → x describe the current, deformed configuration which changes during the simulation as the shape of the sample does, and describes the modelled object as would be observed in real space. Transport phenomena such as the heat equation should therefore be described in terms of the → x coordinates (e.g. the gradient is taken with respect to these coordinates) which are then mapped back to ⎯→ ⎯ the stationary X domain for calculation. This mapping is described by → → ⎯→ ⎯ the vector-valued displacement d = x − X which, under certain conditions such as no tearing, can be considered smoothly differentiable [14]. The deformation gradient F, in which bold font denotes tensors, can then be defined which transforms between coordinate systems:

Laser flash experiments typically employ localized laser heating coupled with fast pyrometry. Composition change is minimized through the high experiment speed, and containment of the heated sample region within the cooler sample periphery; i.e.: the observed region of interest is self-crucibled. Disk shaped samples of industrial grade UO2 1–2 mm thick were mounted in a pressurized cell under a controlled atmosphere of Ar + 2%H2 at 3 atm with optical glass ports to allow for laser and pyrometer transmission. The disks were wrapped in a Teflon ring to prevent broken fragments from falling appart. The wrapped disks were suspended by four radially arranged graphite screws, with the flat face oriented vertically, perpendicular to the horizontal laser. The laser was a TRUMPF® Nd:YAG continuous wave laser radiating at 1064 nm, programmable with a complex power/time profile between 45 and 4500 W in increments of 45 W. Both the laser and the pyrometer were aimed at the centre of the sample. The details of the experimental apparatus are described elsewhere [13]. A photograph of a sample is shown in Fig. 1 in which the seemingly polished area is a result of having melted and resolidifed, and thus shows the position of the laser beam relative to the sample. The important physical phenomena are: 1. 2. 3. 4. 5.

→ ∂→ ∂d x = F = ⎯→ ⎯ ⎯→ ⎯ +I ∂X ∂X

(1)

where I is the identity matrix. The deformation gradient is related to the material strain, which is calculated from the stress balance on the reference domain. In a thermoelastic material, the material stress response is the derivative of the Helmholtz energy by the strain. The principle of frame indifference requires that the free energy be independent of spatial rotations. Therefore, the Lagrangian strain tensor, E = .5(F T F − I ) is an appropriate choice of thermodynamic state parameters. The volume ratio associated with the deformation, J, is calculated by:

Heat transport in the sample, Elasticity in the sample, Heat transport in the buffer gas, Mass transport in the buffer gas of vapour species, and Compressible fluid dynamics in the buffer gas.

J= This section proceeds as follows: Firstly, the governing equations which describe the heat transport and elasticity are derived. Secondly, the hyperelastic thermodynamic model of the system is described with particular emphasis on the heat capacity at constant strain. Finally, the material properties used in the model are given. The resulting system of partial differential equations were solved numerically by the finite element method, in the commercial software package COMSOL Multiphysics, v4.2 on a desktop computer.

ρ0 = det(F ) = ρ

2det(E ) + 1

(2)

where ρ and ρ0 are the current and initial densities respectively. The specific internal energy and Helmholtz differentials may be written in terms of their natural variables, the specific entropy, s, the Lagrangian strain tensor, E , and the temperature, T [14]:

du (s, E ) = Tds +

175

1 S: dE ρ0

(3)

International Journal of Thermal Sciences 132 (2018) 174–185

M.J. Welland et al.

3. Conduction into the screws, 4. Conduction/convection into the buffer gas, and 5. Enthalpy of vaporization of the sample material.

Table 1 Definitions of thermodynamic properties from Helmholtz energy. Expression

Symbol

Name

∂a ∂T ∂a ∂E ∂2a

−s

Entropy

1 S ρ0

Second Piola stress

∂T2 ∂2a

∂T ∂E

−

∂s ∂T

=−

−

∂s ∂E

=

da (T , E ) = −sdT +

Heat capacity at constant strain

CE T

1 ∂S ρ0 ∂T

1 S: dE ρ0

The driver of this experiment was the laser, and the result is sensitive to the radial power distribution and the transmissivity of the optical setup. The total laser power is measured at the source, but the laser must then traverse the optical fiber, focusing apparatus, and window before reaching the sample. Characterization of the laser power was done separately to obtain the correlation between the recorded laser power from the source, and the total power which reached the sample. The laser power profile was measured using a different optical apparatus which was necessary to limit the beam power. Therefore the radial power profile was stretched by a factor corresponding to the different beam spot diameters, which was determined by the simulation of the steady state temperatures. The elasticity boundary conditions were the pressure of the buffer gas and the fixed displacement of the (assumed) incompressible screws imposed on the faces of the screw's contact. The particular boundary conditions depended on the symmetry and scope of the model, as discussed in sections 4 and 5. The system evolution is therefore governed by equations (7b) and (10b), which solve for variables T and d subject to the boundary conditions described above as applicable.

=

1 Γ ρ0

Grüneisen tensor

(4)

A summary of derivation and relationships between thermodynamic variables is reviewed in Table 1. The time derivative of the specific internal energy is:

ρu˙ = ρTs˙ + J −1S: E˙

(5a)

ρu˙ = ρTs˙ + σ : D

(5b)

˙ . The where the Lagrangian time derivative is denoted by superscript □ Cauchy stress is defined as σ = J −1FSF T , and assumed to be symmetric due to conservation of angular momentum [15]. The deformation-rate tensor, D is calculated by E˙ = F T DF [16]. Moving on to the conservation laws for momentum and energy [15,17]: ρ→ v˙ = ∇⋅σ

2.2. Hyperelastic thermodynamic model The heat capacity at constant strain, cE is not a commonly quoted material property but may be related to the heat capacity at constant volume, cV and the shear modulus as shown below. An isotropic, homogeneous, thermoelastic model is considered. Furthermore, multiplicative decomposition is considered in order to separate thermal expansion and elastic deformation [14,18]. The specific Helmholtz energy can be written as:

(6a)

ρu˙ = −∇⋅→ q + σ: D (6b) → where q = −k∇T is the reduced heat flux, and k is the thermal conductivity. Eq. (6a) is the origin of the elastic balance equation which, assuming quasistatic (zero acceleration) conditions becomes: 0 = ∇⋅σ

(7a)

0 = ∇X ⋅S

(7b)

a (T , E ) = aσ = 0 (T ) + ael (T , Eel )

where aσ = 0 is the specific Helmholtz energy at σ = 0 , i.e.: when the material is free to expand without applied stress, and ael (T , E ) follows a Neo-Hookean hyperelastic model:

where the symbol ∇X is the gradient taken with respect to the reference ⎯→ ⎯ coordinate system X . In order to arrive at a heat balance equation, Eqs. (5) and (6b) are compared to eliminate u˙ :

q ρTs˙ = −∇⋅→

ael (T , Eel ) =

(8)

cE ˙ T + J −1Γ : E˙ T

where Jth =

∂σ ρcE T˙ = −∇⋅→ q −T :D ∂T

(10b)

ρ0 ρσ = 0

is the thermal dilation in the absence of stress, a known

material property based on the density as a function of temperature. ρ The elastic dilation, Jel = σρ= 0 therefore describes the change in density between the current deformed shape and what it would be if no stress were applied. The term I1 in equation (12) is the first modified invariant of the elastic Lagrangian strain tensor:

where cE is the heat capacity at constant strain and Γ is the Grüneisen tensor. Eqs. (8) and (9) may be combined to yield an expression for the rate of change of temperature in terms of heat flux and the rate of deformation, (10a)

(12)

(13)

J = Jel Jth

(9)

ρcE T˙ = −∇⋅→ q − TJ −1Γ : E˙

1 [κT (Jel − 1)2 + GT (I1 − 3)] 2

where κT and GT are the isothermal bulk and shear moduli respectively. Jel is the elastic dilation and obeys:

where the elastic deformation has vanished from the expression of entropy production, a consequence of the reversibility of ideal elasticity. The entropy may also be expressed in terms of the variables T and E , such that the partial derivative becomes:

ρs˙ (T , E ) = ρ

(11)

1 Eel = Jth−2/3 ⎡E − (Jth2/3 − 1) I ⎤ 2 ⎦ ⎣

(14)

The first modified invariant is given by:

I1, el = [2⋅tr(Eel ) + 3]⋅Jel−2/3

(15)

Alternately, equation (11) may be decomposed uniquely into terms which depend on the dilation (pure volume change) and another which depends on shear stresses [14]:

The rightmost term, which contains the rate of change of deformation, describes the heat effect of changing the state of strain and lies at the heart of the thermoelastic effects on heat transport. The relevant thermal boundary conditions on the sample were:

a (T , E ) = a v (T , J ) + as (T , E )

(16)

where a v (T , J ) is the Helmholty energy component only dependent on the dilation, and as (T , E ) is the component which depends on shear stresses and is constructed such that as (T , E = 0) = 0 [14].

1. Laser radiation, 2. Radiative heat loss,

176

International Journal of Thermal Sciences 132 (2018) 174–185

M.J. Welland et al.

3. Refinement of the experimental technique through preheating

The heat capacity at constant strain can then likewise be decomposed:

cE = −T

∂2a v (T , J ) ∂2as (T , E ) −T ∂T 2 ∂T 2

In order to reasonably expect accurate simulation results, the experimental process must be well understood and controlled. This motivated a study of the reproducibility of the thermograms and led to the introduction of a step to preheat the sample before the pulse as opposed to beginning the experiment with the sample at room temperature. Data presented in this section at temperatures below approximately 1800 K may not be accurate due to the calibration and sensitivity of the pyrometer, however this does not change the relevance of the cracking. The observed temperatures in Section 4 are from a different pyrometer which is calibrated in the reported temperature range.

(17)

where it is clear that under condition of persevering zero shear stress, cE = c v . The heat capacity at constant volume is readily obtained by the well known relation:

c v = cp − Tα 2κT

(18)

where α is the volumetric thermal expansion. This equation is however strictly valid only when σ = 0 , i.e.: when the current volume is equal to the unstrained, thermally expanded volume, as opposed to an arbitrary, potentially compressed state. The expression for heat capacity at constant volume, for a material which is at non-zero pressure (and therefore elastically dilated) is, to a first expansion in 1 − Jel [19,20]:

c v (T , J ) = c v (T , σ = 0) + + α 2κT

3.1. Starting at room temperature Four sets of experiments with laser powers 135, 180, 270, and 360 W and a 5 mm nominal beam spot1 were performed. The samples were badly fractured after each set and so replaced before the next. The set of pulses at 360 W is representative of the other powers and shown in Fig. 3 with pulse durations 0, 50, and 100 ms beyond the initial power ramp. Each power/pulse duration pair was repeated, and the corresponding thermograms are labeled in chronological order. Immediately obvious is the lack of reproducibility between shots with identical conditions, and furthermore that subsequent shots at the same pulse duration result in increasing heating rates and maximum temperatures. After a couple shots, this trend saturates and reproducibility is achieved as indicated by shots 3 & 4, and 6 & 7. As may be expected, the first shot at a new pulse duration, shots 5 & 8, follow the previous curve, corresponding to the saturation of the previous duration. Of particular note are shots 10 & 11 in which the duration was returned to 0 ms, following the completion of all the other durations. Shot 10 now lies well above the previously saturated curves 3 & 4, instead following the curves corresponding to the 100 ms pulses, indicating an irreversible sample modification has occurred, consistent with a crack. There is a small blip in curve 1 at approximately 1625 K however it is not reproduced in any of the other shots, nor was anything similar systematically observed in other samples not shown here. The origin may be a piece of the sample falling during the first shot, a surface modification on initial heating or thermomechanical effects associated with sudden stress relaxation upon initial cracking. All thermograms show a difference in heating rate at approximately 1625 K, which may be attributed to surface modification, sample movement, or a reflection of the laser as discussed further below and examined in Fig. 6. Often, the sample would crumble when removed from the holder or if the Teflon ring was removed. Of the debris, flakes 80–200 μm in thickness and several larger fragments were typical, indicating that both parallel and perpendicular cracks occurred. Fig. 4a shows an image of typical sample which has fractured substantially, with the flakes shown in the inset. Due to the stochastic nature of the cracking, simulation of these thermograms have little value and is therefore not performed here.

∂κ (T , p) ∂α (T , P ) 1 (1 − Jel ) T ⎡κT + 2α T ⎢ T ∂T ρ ∂ ⎣

∂κT (T , p) ⎤ ⎥ ∂p ⎦

(19)

where the coefficient, α and κT are evaluated at zero pressure. The second term in equation (17) requires the component of the Helmholtz energy which depends on E . Considering the simple Neohookean hyperelastic model in equation (12), it is clear this corresponds to the second term. Therefore the total heat capacity at constant strain can be calculated as:

cE = c v (T , J ) −

T ∂2GT (I1 − 3) 2 ∂T 2

(20)

2.3. Material properties Material properties required by the model are itemized in Table 2. Material property measurements at temperatures approaching 3000 K are scarce, and so extrapolation from low temperature models is common. The enthalpy of vaporization is calculated from the vapour pressure according to the one component Clausius-Clapeyron equation,

h vap = RT 2

d ln(Pvap) dt

(21)

which in fact results in an effective enthalpy of vaporization since the vapour pressure measured and reported in Table 2 is the total vapour pressure of all uranium bearing species [29]. Of particular importance to this work, the elastic moduli for the full temperature range of interest, particularly at higher temperatures as dilation related effects scale roughly linearly with temperature. The neutron scattering results of Clausen et al. [30] were selected as several data points are above 2000 K. The published results are for the elastic stiffness constants of the cubic UO2 structure, c11, c12 and c44 , not for the polycrystaline material. The elastically isotropic elastic parameters, κT and GT are calculated:

κT =

c11 − 2c12 3

(22)

GT =

c11 − 2c12 ⎛ 2 + 3A A ⎞ +5 4 3 + 2A ⎠ ⎝ 5

(23)

3.2. With preheating The series of experiments with a preheating of 45 W for at least 30 s was performed with powers: 225, 360, 450, 540, 675, and 810 W and with an 8 mm nominal beam spot size, where the increase of spot size helps avoid radial temperature gradients. The preheating drastically improved the structural stability of the sample as can be seen in Fig. 4b. This stability eliminated the need for the teflon ring and permitted

inwhich the shear modulus is the arithmetic average of the Voigt and Reuss approximations of effective elastic constants, and the anisotropy factor A = 2c44/(c11 − c12) [31]. A quadratic was fit to each data set as shown in Fig. 2 and given in Table 2.

1 The beam spot size is the nominal size corresponding to the optics used. Auxiliary experiments revealed that the actual spot size can be up to 20% larger due to difficulties focusing the beam.

177

International Journal of Thermal Sciences 132 (2018) 174–185

M.J. Welland et al.

Table 2 Table of material properties using the short notation Tk = T/1000 . Property UO2 cp (T , p = 0)

Jth (T , p = 0)

Equation

2.4570 ⋅ 107e548.68/ T 2 T2 (e548.68/ T − 1)

+ 4.570⋅10−3T +

4.373 ⋅ 1011e−18531.7/ T T2

9.973⋅10−1 + 9.082⋅10−3Tk − 2.705⋅10−4Tk2 + 4.391⋅10−4Tk3

εh kUO2 (T )

Source

J mol−1 K−1

[21] [21]

273 > T > 923

9.9672⋅10−1 + 1.179⋅10−2Tk − 2.429⋅10−3Tk2 + 1.219⋅10−3Tk3

ε1064

Units

[21]

923 > T > 3120

[22]

0.836 + 4.321·10−6 (T − 3120) 0.85 100 7.5408 + 17.692Tk + 3.614Tk2

+ 6400Tk5/2 e−16.35/ Tk⋅

1 − (2.6 − .5Tk ) por 1 − (2.6 − .5Tk )0.05

W m−1 K−1

[21,23] [21]

κT (T , p = 0)

− 1.507⋅10 4T2 + 4.657T + 2.245⋅1011

Pa

Current work

GT (T )

− 8054T2 + 5.223⋅106T + 8.69⋅1010 6.1

Pa

Current work

Pvap

10−31284/ T + 7.616

MPa

[25]

Ar + 2%H2 kAr

[1.11 + 3.66Tk − 0.385Tk2]⋅10−2

W m−1 K−1

[26]

[5.987 + 5.259(T − 1200)]⋅105

W m−1 K−1

[27]

p RT 5 R 2

mol m−3

Ideal gas

J mol−1 K−1

Ideal gas

23⋅10−6

Pa s−1

[28]

∂κT ∂p

k H2 ρAr cp, Ar νmixture

[24]

Fig. 2. Calculated polycrystaline isothermal bulk modulus (left) and shear modulus (right) for UO2 from data in Ref. [30]. Quadratic polynomial fits are shown, and given in Table 2.

reusing the same sample for the sets at 225, 360, and 450 W. Two sets of experiments are described in detail here, at powers of 675 and 810 W. Fig. 5 shows the thermograms with laser power 675 W, and increasing pulse durations. Fig. 6 shows the thermograms with maximum laser power 810 W. All pulses were repeated 3 times. At the end of each series, a shot was made which melted enough surface material to produce a small temperature plateau upon freezing at 3120 K, in excellent agreement with literature values [21] which helps verify correct pyrometer calibration. The reproducibility of the shots in these figures is excellent, although some deviations occur at longer durations. As the temperature increases, the rate of heating curves over as the temperature pulse propagates into the sample. These curves gradually separate with subsequent shots reaching a maximum deviation at, in both figures, pulses 16, 17, and 18, and then appear to recover slightly for shot 19. This behaviour suggests that some parallel cracks are occurring despite the preheating but the melt is filling them in, therefore restoring the original solid upon freezing. Close examination of the beginning of these thermograms, enlarged

Fig. 3. The temperature and laser power profiles for a set of shots at 360 W and a 5 mm beam spot size, with pulse durations (excluding the power ramp) of 0, 50, and 100 ms. Since the laser outputs in discrete levels, the power ramp appears as steps. Thermograms are labeled chronologically. 178

International Journal of Thermal Sciences 132 (2018) 174–185

M.J. Welland et al.

Fig. 4. Microscope image of the front face of the samples following experimental series: a) without preheating corresponding to 3, showing the front face of sample while still in Teflon ring and sample of flakes and small debris in the inset, and b) with preheating after series of shots at 225, 360, and 450 W without the teflon ring since it was unnecessary.

surface and into the pyrometer. The magnitude of this effect changes on each experiment due to modification of the sample surface. During simulation, this factor was fit by trial and error to each thermogram individually in such a way as to eliminate these steps. Since this effect is present during preheating, the initial temperature in Figs. 5 and 6 is incorrect and, when the correction is applied, drops to approximately 1300 K. 4. Calibration of heat losses Heat losses outlined in 2.1 were calibrated by a series of experiments in which the laser was held at a particular power for 30 s such that the sample reached steady state. The same laser optics and configuration used for the preheating tests in Section 3.2 were used. The power was then increased and the process repeated for a series of powers, providing laser power vs. surface temperature data. The steady state problem does not show such high axial gradients in temperatures and stresses and is therefore more computationally tractable. Additionally, and perhaps more importantly, in removing transients from equation (10), the heat capacity effects and the Grüneisen tensor are eliminated, allowing the fitting of various parameters associated with heat loss, discussed in section 2.1 without introducing the uncertainty of the material properties associated with the thermoelastic effects or the complexity of its calculation on heat transport. In examining Fig. 1, it can be seen that due to the screws and the vertical movement of the buffer gas over the sample, this problem is three dimensional and exhibits mirror symmetry across the plane dividing the vessel in two, and the vertical axis. A 3-D model which included all the physical phenomena outlined in 2.1 was therefore developed for the steady state experiments. The fluid flow in the buffer gas was simulated concurrently in order to capture natural convection which can be seen to be relevant from the vapour deposit in Fig. 1. In particular, coupled heat transport, mass diffusion and Navier-Stokes fluid dynamics were simulated in the adjacent domain which was limited to the apparatus geometry. The concentration of the vapour in the buffer gas domain was set to the vapour pressure in Table 2 as a boundary condition, and therefore the vaporization rate was determined by convection/diffusion into the buffer gas. The coupling between the sample and this domain implies two types of convection-augmented heat loss were captured mechanistically: heat conduction into the moving gas, and the diffusion of the vapour species, which drives further vaporization and associated enthalpy loss. The heat flux at the boundary was calculated as:

Fig. 5. The temperature and laser power profiles for a set of shots at 675 W and an 8 mm beam spot size and multiple pulse durations labeled chronologically. Each pulse was repeated three times to verify reproducibility. Pulse 19 is included which shows a thermal inflection from freezing at 3120 K to verify the calibration of the pyrometer calibration.

Fig. 6. The temperature and laser power profiles for a set of shots at 810 W and an 8 mm beam spot size, with increasing pulse durations labeled chronologically. Each pulse is repeated three times to verify reproducibility. Pulse 19 shows a thermal inflection from the freezing at 3120 K to verify the pyrometer calibration. The inset shows an enlargement of the beginning of the pulse and the reflection of the laser exaggerating the thermogram signal.

→ → q ⋅nˆ = −plaser (t ) plaser profile (r ) tsapphire ε1064 + εh σSB T 4 + h vap Jvap⋅nˆ

in the inset of Fig. 6, shows a pyrometer signal which exactly reproduces the laser profile down to a scalar factor. As such a high rate of temperature increase is not realistic, therefore this signal must arise from a fraction of the laser light which reflects (diffusely) off the sample

(24)

q is the heat flux, nˆ is the normal vector, plaser (t ) is the time where → dependent power emitted by the laser, plaser profile (r ) is the fractional laser power density as a function of radius, r, tsapphire is the 179

International Journal of Thermal Sciences 132 (2018) 174–185

M.J. Welland et al.

transmissivity of the sapphire window, ε1064 is the emissivity at 1064 nm which was assumed to be equal to the absorptivity by Kirchhoff's law. εh is the total-hemispherical emissivity, σSB is the Stephan-Boltzmann constant. h vap is the enthalpy of vaporization and Jvap was the total flux of vapour (diffusive and convective) away from the surface, calculated simultaneously by a simulation of the natural convection in the buffer gas. By coupling heat conduction into naturally-circulating buffer gas convection simultaneously, the heat loss by conduction and convection into the buffer was also accounted for implicitly. The radiative heat transport term neglects absorbed radiation from the apparatus as a simplification since the temperature of the vessel, which was measured using a thermocouple during the experiment and never deviated from room temperature by a significant degree, is much smaller than the surface temperature. In analysis of the solution, this term is also insignificant compared to heat loss to the buffer gas. The screws were assumed to make full cross-sectional contact, although in the experiments this was seldom the case. The screws were not modelled explicitly, instead being accounted for by the equation:

→ q ⋅nˆ = kscrews (T − Tamb)

Fig. 8. The steady state temperature in the centre of the sample's front surface plotted as a function of the laser power, and simulation results.

(25)

at the contact areas. The value of kscrews was fit to reproduce the experimental data. For the steady state elastic calculations, the screws were assumed to be completely rigid such that they did not deform. The remaining surfaces were simply subjected to the external pressure, 3 atm applied by the buffer gas in the deformed coordinate system. The two parameters which were fit in the steady state simulations were kscrews and the beam spot radius of the defocused laser beam. The thermal conductivity of the screws, kscrews , was found to be 2 × 102 Wm−1K−1, comparable to literature values of graphite [32]. The fractional power density of the defocussed beam, plaser profile (r ) was measured in auxiliary experiments using a 2D profilometer and a different beam spot size. The profile was then scaled radially, with a corresponding reduction in power density such that the integral remained constant, by a parameter fit to reproduce the experimental observations which also matched exposure of photographic paper. The laser power profile is shown in Fig. 7, and is not changed throughout the remainder of this work. The results of the fitting are shown in Fig. 8 and show excellent agreement with observations. Simulation results of the temperature and vapour concentrations corresponding to 45 W, the power used for preheating, are shown in Fig. 9. Streamlines in the buffer gas are shown as white lines and qualitatively match the yellow deposit shown in Fig. 1.

5. Transient experiment The final series of experiments were of transient pulses, which are representative of the laser flash melting experiments the model is being developed to help analyze. As a simplification, the convective heat loss terms (conduction and enthalpy of vaporization) were simplified. The heat losses due to heat conduction and mass transport (with associated enthalpy of vaporization) in the buffer gas were examined as a function of the surface temperature, and a simplified heat loss equations were replaced by a term of the form hbuffer (T − Tamb) where hbuffer = 300 kAr , where kAr is the thermal conductivity of argon, which was found to reproduce the results shown in Fig. 8 to a high degree of accuracy. Once the steady state powers were well simulated and the generalized convective heat loss parameter fit, the transient simulations could proceed without further fitting. Since the transient model only considered heat transport and elasticity in the sample on an axi-symmetric 2-D domain, it could be solved with much higher resolution which is required to accurately represent the large temperature and strain gradients present in the transient experiments. The simplified heat loss boundary conditions for the 2-D antisymmetric model were:

→ q ⋅nˆ = −plaser (t ) plaser profile (r ) twindows ε1064 + εh σSB T 4 + hbuffer (T − Tamb) (26) Additionally, the heat loss to the screws was averaged to obtain as axisymmetric equivalent term. Specifically, the kscrew value found from the steady state calculation was divided by the ratio of areas between the approximate screw contact area, and a cylinder of the same height and added to the sides of the geometry. It was noted in the 3-D model that, even with the screws being considered completely unyielding, the remainder of the free surface of the sample expanded unabated. Therefore, in the 2-D approximation, only the external buffer gas pressure was applied as the elastic boundary conditions. The results of the simulation corresponding to thermograms at 675 W and 810 W, first shown in Figs. 5 and 6 respectively, are shown in Figs. 10 and 11, along with the first shot of each pulse duration. As an example of the temperature and stress profiles in the sample, plots of the temperature and pressure at 290 ms in Fig. 11 are given in Figs. 12 and 13 respectively. Also visible is the extent of predicted deformation of the sample. The calculated Cauchy stress components for the 810 W series in Fig. 11 are shown in Fig. 14 along with the hydrostatic pressure over time. The calculated Lagrangian strain components are shown in Fig. 15, and the calculated density compared to the density at zero pressure is shown in Fig. 16.

Fig. 7. The fractional laser power density as a function of radius, determined by scaling a recorded profile for a smaller spot diameter in auxiliary experiments by a factor determined to match photographic paper exposure and for the steady state simulation observations. 180

International Journal of Thermal Sciences 132 (2018) 174–185

M.J. Welland et al.

Fig. 9. Simulation results for steady state shot at 45 W, corresponding to the preheating power. Plots (a) and (b) show the temperature and vapour concentration on the plane of the irradiated sample surface, and plots (c) and (d) on the plane defined by the vertical and the surface normal with the laser on the right. Streamlines are plotted in white.

181

International Journal of Thermal Sciences 132 (2018) 174–185

M.J. Welland et al.

6. Discussion 6.1. The preheating technique The samples heated without preheating are largely affected by cracking, both structurally and in terms of their temperature response. Cracking results in the fragmentation of the sample by perpendicular cracks, and formation of flakes at the surface due to parallel cracks. Poor heat transport across the parallel cracks are primarily responsible for retarding heat transport from the surface into the depth of the sample, where this resistance is considered to be dependent on the depth of the crack from the surface, the area of solid-solid contact and the distance of solid-solid separation. Repeated shots at the same pulse power and duration appear to form the crack, then shift the pieces further apart until a certain point is reached at which point subsequent shots yield the same thermogram, effectively saturating the effect. Upon a new pulse duration, the thermal behaviour follows that of the previous duration and extends to a higher temperature. Subsequent shots at the longer duration follow new paths, indicating the formation of new cracks, or the separation of existing cracks occurs on cooling. The simulation results shown in Fig. 12 show the penetration depth of the heat pulse exceeds the thickness of the flakes shown in Fig. 4a, approximately 100 μm which supports the notion that the formation of these flakes were presenting a thermal barrier in the experiments without preheating.

Fig. 10. Simulation results compared to observed thermograms for pulses of various durations at 675 W, showing experimental data in broken lines and simulation results as full lines.

6.2. With preheating The effect of preheating on avoiding significant cracking is obvious, and reproducibility of subsequent shots is significantly enhanced. It is unclear from the current analysis if this is due to the reduced temperature gradient incurred during the pulse or the material properties at the preheating temperature. Some deviations occur with longer pulse durations indicating that a higher preheating temperature might be required to avoid cracking at these conditions although such a proposal must be weighed against the risk of sample reaction or compositional shift at elevated temperatures. This point is emphasized by the appearance of a ring of orange on the periphery of Fig. 1 indicative of vaporization from the hot centre and condensation on the cooler periphery. An interesting feature shown in Fig. 5 is the appearance of two slight inflection points in the ascending flank, located at approximately 2600 and 2900 K. In other experiments at different powers, similar features persist although corresponding neither in temperature nor time. The temperatures are not consistent enough to be attributed to the Bredig transition in UO2, known to occur at 2670 K [33], which has too

Fig. 11. Simulation results compared to observed thermograms for pulses of various durations at 810 W, showing experimental data in broken lines and simulation results as full lines.

Fig. 12. Simulated temperature corresponding to 290 ms on the 810 W pulse shown in Fig. 11.

182

International Journal of Thermal Sciences 132 (2018) 174–185

M.J. Welland et al.

Fig. 13. Simulated pressure corresponding to 290 ms on the 810 W pulse shown in Fig. 11. The negative pressure in the middle of the sample implies at least one of the stress components to be significantly tensile, therefore perhaps causing fracture in the sample at these points.

small a heat effect to be seen on heating or cooling in these experiments. Also noted is that the anomaly only occurs on heating, not cooling. Finally, none of the samples shot without preheating showed this phenomenon. The simulation results show that, for these conditions, the temperature perturbation only penetrates on the order of 100 μm into the sample, roughly the width of the flakes formed when not preheating, thus supporting the hypothesis that cracks, which may be microscopic, are not present to the extent that they form thermal barriers in the preheated samples. 6.3. Thermomechanical effects Fig. 14. Calculated Cauchy stress components, and the hydrostatic pressure at the observed point corresponding to the 810 W series in Fig. 11.

The simulation results shown in Figs. 10 and 11 reproduce the observed thermograms well considering the uncertainty in material property data at these temperatures. The preheating technique was successful in controlling the cracks, which enables the elastic treatment. By using the steady state experiments to fix the heat loss conditions, the uncertainty associated with the thermoelastic material properties, such as the effective heat capacity, are avoided. By adapting the heat losses calculated and fit in the 3D steady state model, a 2D model is possible which is much more computationally tractable at the very high mesh resolution required to accurately resolve the large temperature and stress gradients during the transient experiments. In the centre of Fig. 13, the model predicts a negative hydrostatic pressure, which means at least one of the stress components must be highly tensile. This may cause fracture in the material, producing a parallel crack that is deep inside the sample. Indeed, some thicker samples did show such cracks but they were too far from the surface to affect the heat transport directly, although they would affect the thermoelastic behaviour. The remaining discrepancy between the simulation results and observed temperature may be explained as resulting from the extrapolated material properties, unaccounted for cracks affecting the thermomechanics if not the heat transport directly, or errors associated with averaging the boundary conditions from the 3D model to the 2D one. No fitting of parameters was performed during the transient simulations. The good correspondence between the simulation results and the observations is a direct reflection of the reasonable quality of the high temperature material properties and the ability of the model to correctly account for the relevant phenomena, in this case the impact of thermoelastics on heat transport and storage. In this work, the thermoelastic effect modified the effective heat capacity, and similarly the enthalpy, depending on the local state of strain. As shown in Fig. 17 for the point observed by the pyrometer, this effect can account for a reduction in heat capacity by as much as 25%. It is important to note that the effect on the heat capacity/enthalpy depends on the local strain state, and therefore this factor must be

Fig. 15. Calculated Lagrangian strain components at the observed point corresponding to the 810 W series in Fig. 11.

Fig. 16. Calculated density compared to the zero pressure density at the observed point corresponding to the 810 W series in Fig. 11.

183

International Journal of Thermal Sciences 132 (2018) 174–185

M.J. Welland et al.

paramount in accurately benchmarking such experimental conditions. Further and even more ambitious applications are currently foreseen, such as the investigation of irradiated and spent nuclear fuels, corium and lava core samples coming from Chernobyl and Fukushima, and tests on new materials (typically hafnium carbo-nitrides) which are predicted to melt at even higher temperatures than the aforementioned carbides [35]. 7. Conclusions Cracking is demonstrated to have a large impact on the reproducibility of fast laser-flash heating experiments in which the surface temperature is observed. Preheating to approximately 1300 K prior to the pulse is shown to drastically reduce cracking, increasing the reproducibility of the thermogram and the structural stability of the samples even after twenty shots, an important consideration when investigating materials of limited availability (e.g.: actinide compounds), and increased fragility due to, e.g., irradiation damage. The separation of experiments into a series of steady state conditions, useful for calibration of model parameters, and a transient series for validation of the model produced good results. The thermoelastic effect is apparently an important consideration in modelling the heat transport in such laser flash experiments where cracking is avoided. The analysis performed in this work has paved the way to a significant improvement of the fast laser heating technique, permitting the observation of the melting behaviour in particularly challenging materials, such as heavily radiation damaged transmutation fuels of very high-melting ultra-refractory ceramics.

Fig. 17. Comparison of the effective heat capacity under conditions of constant pressure, constant volume (at atmospheric pressure), constant volume at the current volume, and constant strain at the current strain. Note that neither the heat capacity at current volume nor current strain is a single path, since the volume/strain at a given temperature during the simulation varied with time.

considered locally. Upon melting, the strain in the solid would be released by fluid flow, and the effective enthalpy of fusion would increase beyond the value typically measured (at constant pressure) to include that which was removed by the strain state [12]. Similarly, the local state of strain would also modify the fusion temperature, akin to the Clausius-Clapeyron relation, although this analysis would be complicated by the condition relevant to a strained solid and free flowing liquid. Solidification would be more complicated, as the liquid would not likely solidify into the same strain-state as the underlying solid. Analysis of this effect is beyond the scope of this work.

Acknowledgements The European Commission and Atomic Energy of Canada Ltd. is thanked for the funding provided to support the current work. References

6.4. Applications [1] F. De Bruycker, K. Boboridis, R.J.M. Konings, M. Rini, R. Eloirdi, C. Guéneau, N. Dupin, D. Manara, On the melting behaviour of uranium/plutonium mixed dioxides with high-Pu content: a laser heating study, J. Nucl. Mater. 419 (1–3) (December 2011) 186–193. [2] R. Böhler, M.J. Welland, F. De Bruycker, K. Boboridis, A. Janssen, R. Eloirdi, R.J.M. Konings, D. Manara, Revisiting the melting temperature of NpO2 and the challenges associated with high temperature actinide compound measurements, J. Appl. Phys. 111 (2012) 113501. [3] R. Böhler, M.J. Welland, D. Prieur, P. Cakir, T. Vitova, T. Pruessmann, I. Pidchenko, C. Hennig, C. Guéneau, R.J.M. Konings, D. Manara, Recent advances in the study of the UO2 - PuO2 phase diagram at high temperatures, J. Nucl. Mater. 448 (1) (2014) 330–339. [4] O. Cedillos-Barraza, D. Manara, K. Boboridis, T. Watkins, S. Grasso, D.D. Jayaseelan, R.J.M. Konings, M.J. Reece, W.E. Lee, Investigating the highest melting temperature materials: a laser melting study of the TaC-HfC system, Sci. Rep. 6 (2016) 37962. [5] D. Manara, L. Soldi, S. Mastromarino, K. Boboridis, D. Robba, L. Vlahovic, R. Konings, Laser-heating and radiance spectrometry for the study of nuclear materials in conditions simulating a nuclear power plant accident, JoVE 130 (2017) e54807. [6] D. Manara, C. Ronchi, M. Sheindlin, M. Lewis, M. Brykin, Melting of stoichiometric and hyperstoichiometric uranium dioxide, J. Nucl. Mater. 342 (1–3) (June 2005) 148–163. [7] M. Sheindlin, D. Staicu, C. Ronchi, L. Game-Arnaud, B. Remy, A. Degiovanni, Experimental determination of the thermal conductivity of liquid UO2 near the melting point, J. Appl. Phys. 101 (9) (2007) 093508. [8] L. Vlahovic, D. Staicu, A. Kuest, R.J.M. Konings, Thermal diffusivity of UO2 up to the melting point, J. Nucl. Mater. 499 (2018) 504–511. [9] C. Ronchi, M. Sheindlin, M. Musella, G.J. Hyland, Thermal conductivity of uranium dioxide up to 2900 K from simultaneous measurement of the heat capacity and thermal diffusivity, J. Appl. Phys. 85 (2) (1999) 776. [10] T.R. Pavlov, M.R. Wenman, L. Vlahovic, D. Robba, R.J.M. Konings, P. Van Uffelen, R.W. Grimes, Measurement and interpretation of the thermo-physical properties of UO2 at high temperatures: the viral effect of oxygen defects, Acta Mater. 139 (2017) 138–154. [11] M.J. Welland, B.J. Lewis, W.T. Thompson, A comparison of Stefan and Phase Field modeling techniques for the simulation of melting nuclear fuel, J. Nucl. Mater. 376 (2) (May 2008) 229–239. [12] M.J. Welland, The effect of structural mechanics on the simulation of melting

The analysis performed in this work has paved the way for a significant improvement of the fast laser heating technique on projects that have been performed while this document was prepared. The present investigation has led to the introduction of a laser pre-heating stage and multiple shots between which the sample is maintained at temperatures well above 298 K (typically not lower than 1500 K). This novel approach has drastically reduced thermal shocks, preventing the investigated material from cracking during the experiments. This has permitted the mechanical stabilization of refractory materials for several seconds at temperatures up to and beyond the melting point. It has thus become possible to use this approach to obtain repeated measurements of the solid/liquid transition of particularly challenging material samples, some of which could hardly be performed previously. For example, the pre-heating and multiple-shot laser heating technique has been employed, in recent years, for the melting behaviour analysis of highly radioactive nuclear transmutation fuel samples, containing large amounts of 241Am (up to 20 mol %), and therefore extremely brittle due to the self-irradiation damage [34]. Another significant example concerns the application of the laser heating technique, improved thanks to the results of the current analysis, to socalled ”ultra refractory ceramics”, such as tantalum and hafnium carbides [4]. Such compounds display the highest melting temperatures ever observed at ambient pressure (well above 4000 K), and are therefore important candidates for aerospace applications and other extreme conditions. It is clear, that the huge thermal gradients and shocks induced by laser heating would have quickly resulted in the destruction of such samples before completion of any significant measurement, unless the experimental conditions had been carefully studied and set-up before the tests. The present investigation has been 184

International Journal of Thermal Sciences 132 (2018) 174–185

M.J. Welland et al.

[13]

[14]

[15] [16] [17] [18]

[19] [20] [21] [22] [23]

[24]

nuclear fuel, Proceedings of the International Youth in Nuclear Congress 2010 spring Meeting, 2010, pp. 12–18 number 109, pages. D. Manara, M. Sheindlin, W. Heinz, C. Ronchi, New techniques for high-temperature melting measurements in volatile refractory materials via laser surface heating, Rev. Sci. Instrum. 79 (11) (November 2008) 113901. B.J. Plohr, J.N. Plohr, Large Deformation Constitutive Laws for Isotropic Thermoelastic Materials, Technical report Los Alamos National Laboratory, 2005 Report No. LA-UR-05–5471; July; 2005. D. Jou, J. Casas-Vazquez, G. Lebnon, Extended Irreversible Thermodynamics, Springer, 2001. J. Lubliner, Plasticity Theory, revised edition Pearson Education, Inc, 2006, p. 471. H.J. Kreuzer, Nonequilibrium Thermodynamics and its Statistical Foundations, Oxford University Press, New York, 1981. V.A. Lubarda, Constitutive theories based on the multiplicative decomposition of deformation gradient: thermoelasticity, elastoplasticity, and biomechanics, Appl. Mech. Rev. 57 (2) (2004) 95. D.C. Wallace, Thermodynamics of Crystals, Wiley, 1972. G. Grimvall, Thermophysical Properties of Materials, Elsevier Science Pub. Co. Inc, New York, NY, 1999. J.K. Fink, Thermophysical properties of uranium dioxide, J. Nucl. Mater. 279 (1) (March 2000) 1–18. M. Bober, H.U. Karow, K. Muller, Study of the spectral reflectivity and emissivity of liquid ceramics, High. Temp. - High. Press. 12 (1980) 161. J.H. Harding, D.G. Martin, P.E. Potter, Thermophysical and Thermochemical Properties of Fast Reactor Materials, Technical report Harwell Laboratory, 1989, p. 107. M.C. Pujol, M. Idiri, L. Havela, S. Heathman, J. Spino, Bulk and Youngs modulus of

[25] [26] [27] [28] [29] [30]

[31] [32] [33] [34]

[35]

185

doped UO2 by synchrotron diffraction under high pressure and Knoop indentation, J. Nucl. Mater. 324 (2–3) (January 2004) 189–197. M. Tetenbaum, P.D. Hunt, Total pressure of uranium-bearing species over oxygendeficient urania, J. Nucl. Mater. 34 (January 1970) 86–91. F.M. Faubert, Measurement of the thermal conductivity of argon, krypton, and nitrogen in the range 800 to 2000 K, J. Chem. Phys. 57 (6) (1972) 2333–2340. N.C. Blais, J.B. Mann, Thermal conductivity of helium and hydrogen at high temperatures, J. Chem. Phys. 32 (5) (1960) 1459. B.A. Younglove, H.J.M. Hanley, The viscocity and thermal conductivity coefficients of gaseous and liquid argon, J. Phys. Chem. Ref. Data 15 (4) (1986) 1323. W. Breitung, K.O. Reil, Oxygen self and chemical diffusion coefficients in UO2 ±x , Nucl. Sci. Eng. 101 (1) (1989) 26. K. Clausen, W. Hayes, M.T. Hutchings, J.K. Kjems, J.E. Macdonald, R. Osborn, Lattice dynamics and elastic constants of uranium dioxide at high temperatures investigated by neutron scattering, High Temp. Sci. 19 (1985) 189–197. H.M. Ledbetter, E.R. Naimon, Relationship between single-crystal and polycrystal elastic constants, J. Appl. Phys. 45 (1) (1974) 66. C. Uher, 4.3.2 Temperature Dependence of thermal Conductivity of Graphite, Pages 430–439, Springer Berlin Heidelberg, Berlin, Heidelberg, 1991. J.P. Hiernaut, G.J. Hyland, C. Ronchi, Premelting transition in uranium dioxide, Int. J. Thermophys. 14 (2) (1993) 259–283. D. Prieur, F. Lebreton, M. Caisso, P.M. Martin, A.C. Scheinost, T. Delahaye, D. Manara, Melting behaviour of americium-doped uranium dioxide, J. Chem. Thermodyn. 97 (2016) 244–252. Q.J. Hong, A. van de Walle, Prediction of the material with highest known melting point from ab initio molecular dynamics calculations, Phys. Rev. B 92 (020104) (Jul 2015).