Abstract

This work reports a method to measure thermal diffusivity of thin disk samples at high temperatures (approximately 900–1150 K) using a modified Ångström's method. Conventionally, samples are heated indirectly from the surroundings to reach high temperatures for such measurements, and this process is time-consuming, typically requiring hours to reach stable temperatures. In this work, samples are heated directly in a custom instrument by a concentrated light source and are able to reach high steady-periodic temperatures in approximately 10 min, thus enabling rapid thermal diffusivity characterization. Further, existing Ångström's methods for high temperature characterization use thermocouples for temperature detection that are commonly attached to samples via drilling and welding, which are destructive to samples and introduce thermal anomalies. We use an infrared camera calibrated to 2000 °C for noncontact, nondestructive, and data-rich temperature measurements and present an image analysis approach to process the infrared (IR) data that significantly reduces random noise in temperature measurements. We extract amplitude and phase from processed temperature profiles and demonstrate that these metrics are insensitive to uncertainty in emissivity. Previous studies commonly use regression approaches for parameter estimation that are ill-posed (i.e., nonunique solutions) and lack rigorous characterization of parameter uncertainties. Here, we employ a surrogate-accelerated Bayesian framework and a “no-U-turn” sampler for uncertainty quantification. The reported results are validated using graphite and copper disks and exhibit excellent agreement within 5% as compared to reference values obtained by other methods.

1 Introduction

Accurate thermal diffusivity characterization at high temperatures is crucial in the aerospace and energy industries, among others. Currently, the laser-flash [15], transient plane source [6], hot wire [7], and Ångström's [8,9] methods are commonly used for such characterization. Compared to other methods, Ångström's method is robust and simple to realize. Typically, a portion of a sample is periodically heated, and thermal diffusivity is extracted by analyzing amplitude and phase differences between temperature profiles along the heat conduction direction. The method typically does not require knowledge or precise control of (1) contact resistance between the sample and the heater, (2) heat input to the sample, and (3) sample heat loss to the surroundings.

Ångström's method has been used for thermal diffusivity characterization at high temperatures in a relatively small number of previous studies that can be categorized into in-plane and cross-plane measurements. For in-plane measurements, a sample is typically circular in shape and heated to high temperatures by a furnace. Upon reaching steady temperatures, the curved outer surface is heated periodically by the furnace, whereas flat faces are often thermally insulated. Thermal waves are propagated along the radial direction, two thermocouples are placed at different radii to measure oscillatory temperature profiles, and thermal diffusivity is computed from temperature amplitude decay, phase shift, or both. Katsura [10] used such a technique on silica glass up to 9 GPa and 1200 K with approximately 3% accuracy. Xu et al. [11] used a similar apparatus and measured olivine, wadsleyite, and ringwoodite up to 1373 K with 12% accuracy. For thin sheets or foils, Hatta et al. [12] used a mask to cover part of a strip sample and irradiated both the sample and the mask using a lamp with periodic intensity. Thermal diffusivity was calculated by measuring the decay of temperature oscillation at multiple locations within the masked region. This method achieved 5% accuracy for in-plane thermal diffusivity of nickel foil up to 500 °C.

For cross-plane measurements, a cylindrical sample can be periodically heated at either a top or bottom surface. Typically, two or three thermocouples are placed at different axial locations to measure temperature oscillations, and thermal diffusivity is extracted based on amplitude decay and phase shift. Sidles and Danielson [9] developed an apparatus to measure copper and nickel rods (>50 cm) up to 500 °C. Abeles et al. [13] reported an apparatus that measures thermal diffusivity of short cylindrical solids (Armco iron and germanium) up to 1000 °C with an accuracy of approximately 2% using relatively high heating frequencies. Vandersande and Pohl [14] reported an instrument that is more convenient to assemble and is capable of measuring thermal diffusivity between 80 and 500 K with 3%–7% accuracy. An alternative technique for cross-plane thermal diffusivity measurements was developed by Cowan [15]. The method is similar to laser flash and is suitable for electrically conductive thin disks. In Cowan's method, a high-intensity electron gun with uniform intensity distribution bombards the front side of the disk periodically. The temperature phase difference between the front and the back surfaces of the disk is used to determine thermal diffusivity. Wheeler [16] used this approach and measured refractory metals up to 3000 K with 5% accuracy. Tanaka and Suzuki [17] employed a similar technique on pyrolytic graphite up to 1900 K with 5% accuracy.

Most of the foregoing methods require a heat source external to the sample to heat it and the surroundings to high temperature prior to periodic heating. This process is typically lengthy, requiring hours to reach steady, elevated temperature. To mitigate heat loss, radiation shields are often used that preclude the use of noncontact, nonintrusive temperature detection. Thermocouples are commonly used, and drilling or welding is typically involved to attach thermocouples to the sample, introducing heat losses and potentially altering the sample's structure. In addition, temperature measurements are limited by the number of thermocouples, leading to unreliable results if defects are present between the thermocouples [18].

In this work and in contrast to methods that heat samples from surroundings, we employ a concentrated light source to directly heat samples (3.5″ outer diameter, OD) thin disks) to high temperatures to measure in-plane thermal diffusivity. Upon reaching a steady target temperature (typically under 10 min), we modulate the light source's heat flux periodically. This approach significantly decreases the total time for characterization and increases testing throughput. Similar to Hatta et al.'s approach [12], we partially block the light source such that only the center of the sample is heated. However, several key differences exist. First, in Hatta's work, samples were heated to high temperatures (up to 500 °C) from the surroundings. Second, in Hatta et al.'s work, the radiation losses were linearized, whereas in our case, temperature gradients along the direction of heat conduction exist even under steady heating, and nonlinear radiation losses are included in the model. Third, Hatta's work requires thermocouples for temperature detection. In our study, an infrared (IR) camera calibrated up to 2000 °C is employed for noncontact and nonintrusive temperature detection. Moreover, IR thermography offers data-rich temperature profiles that increase measurement reliability [18] and decrease uncertainties for parameters of interest [19]. Here, we employ a Bayesian framework that benefits from the data-rich IR thermography to quantify uncertainty rigorously for the parameters of interest.

2 Physical System Modeling

Existing Ångström's methods for high temperature characterization employ linearized radiation loss models [9,12,13,15]. These methods can be accurate if the sample's temperature T remains uniform along the direction of heat conduction and if the temperature oscillation amplitude is negligible compare to T. In this work, a thin disk sample (88.50 mm OD) is partially heated using a concentrated light source and relatively significant temperature gradient exists along the heat conduction direction. Therefore, radiation losses cannot be linearized, and we instead develop a numerical solution to account for nonlinear radiation loss.

The schematic of our system is shown in Fig. 1(a). The concentrated light source is a 10 kW Xenon short arc lamp with Lorentzian intensity distribution [20,21]. The power of the concentrated light source is controlled using a 0–10 V external control voltage. The heat flux of the light source is given as
(1)

where A(V) is nominal peak heat flux and is determined by the control voltage V, and σsolar indicates the width of the distribution and is determined by the reflector of the light source and is relatively independent of V; r is the radius on the sample. We use a function generator to produce sinusoidal voltage signals to control the light source to output periodic heat flux. Details for characterizing and modeling the light source are provided in the Supplemental Material S1 on the ASME Digital Collection.

Fig. 1
(a) Schematic of the experimental system for thermal diffusivity characterization, (b) ray tracing study for intensity distribution of the concentrated light source with and without light blockers, and (c) a window function to model the truncated Lorentzian intensity distribution
Fig. 1
(a) Schematic of the experimental system for thermal diffusivity characterization, (b) ray tracing study for intensity distribution of the concentrated light source with and without light blockers, and (c) a window function to model the truncated Lorentzian intensity distribution
Close modal

In our study, the concentrated light source is partially blocked by three annular tungsten rings (referred to as light blockers) so that only the center of the sample is exposed to the light source. A custom tray tracing code [21] was employed to evaluate the intensity distribution of the concentrated light source as truncated by the light blockers. The simulated intensities with and without truncation are shown in Fig. 1(a). A window function, defined by the ratio of intensities with and without light blocker, is shown in Fig. 1(c). As shown by the window function, the light blockers only affect the intensity distribution near the inner edge of the light blocker. The window function in Fig. 1(c) shows the simulated truncated intensity distribution with the light blockers.

Having introduced the concentrated light source, we present a finite difference solver for the transient heating process, under the following assumptions to simplify the analysis:

  1. The fin approximation is made for these thin disk samples, because the thermal penetration depth is much larger than the sample thickness [12].

  2. All surfaces are gray and diffuse.

  3. Thermal expansion of the sample is not considered.

  4. Heat loss via conduction and convection is neglected.

  5. To facilitate IR thermography, a thin layer of high-emissivity coating was applied to the back side of the sample. This layer contributes approximately 2% to the total sample mass and is not considered in the model.

  6. We use the following expression [2224] for temperature-dependent thermal diffusivity
    (2)

This empirical form has been shown to fit a broad range of conditions, largely because it captures the main effects of increased thermal carrier scattering with increasing temperature.

A schematic of the finite difference network is shown in Fig. 2(a) along with a control volume (dashed box) for a node m at time-step p. Applying conservation of energy to the control volume shown in Fig. 2(a) we obtain
(3)
where qr and qr+dr indicate heat transfer via conduction, expressed as
(4)
(5)
qsolar indicates energy absorbed from the concentrate light source
(6)

where η indicates sample absorptivity from the concentrated light source, and qs is given in Eq. (1).

Fig. 2
(a) Finite difference scheme for the transient heating process with the radiation network. The dashed rectangular box indicates a control volume for node m on the sample at time-step p. The shaded area indicates the concentrated light and the solid curves indicate the intensity distribution before and after the light blockers. (b) Temperature distributions on the light blocker. For simplicity, we use three discrete temperatures for modeling.
Fig. 2
(a) Finite difference scheme for the transient heating process with the radiation network. The dashed rectangular box indicates a control volume for node m on the sample at time-step p. The shaded area indicates the concentrated light and the solid curves indicate the intensity distribution before and after the light blockers. (b) Temperature distributions on the light blocker. For simplicity, we use three discrete temperatures for modeling.
Close modal
The net heat flow between the light blocker and the sample, qlb,net, is calculated using a radiation network as shown in Fig. 2(a). Based on experimental observations (see Supplemental Material S3 on the ASME Digital Collection for light blocker temperature measurements), the temperature gradient of the light blocker along the radial direction is typically low (approximately 50 K). Therefore, we use three discrete constant temperatures to approximate the light blocker temperature distribution in different regions as shown in Fig. 2(b). We assign Tlb,1 as the temperature between radius Rlb,H and Rlb/3, where Rlb,H is the radius of the light blocker's hole and Rlb is the radius of the light blocker. Similarly we assign Tlb,2,Tlb,3, and Tlb,H for the region bounded by Rlb/3 to 2Rlb/3,2Rlb/3 to Rlb and 0 to Rlb,H, respectively. qlb,net is expressed as
(7)
where σ is the Stefan-Boltzmann constant; εlb is the emissivity of the light blocker; and εfront is the emissivity of the front side of the sample. AlbH, Alb1,Alb2, and Alb3 are the light blocker areas corresponding to discrete temperatures TlbH, Tlb1,Tlb2, and Tlb3, respectively. FlbHm,Flb1m,Flb2m,Flb3m are view factors from AlbH, Alb1,Alb2, and Alb3 to the differential ring element at node m. Cm,lb and qm,lb are expressed as
(8)
(9)
qe,back indicates radiation loss at the back side of the sample
(10)

where εback is the emissivity of the back side of the sample.

The term qsurr,back indicates the rate of radiation absorbed by the back side of the sample. We assume the surrounding temperature in the vacuum chamber at the back side of the sample is constant. Therefore, qsurr,back is also a constant. More details for calculating qsurr,back are provided in the Supplemental Material S1 on the ASME Digital Collection.

The energy accumulation term E˙ is expressed as
(11)
Substitution of Eqs. (4)(11) into the governing equation Eq. (3) yields
(12)

The initial (IC) and the boundary (BC) conditions are

(13a)

where TW indicates the temperature of the vacuum chamber wall. Equation (12) is discretized using a central-difference implicit finite difference scheme. A detailed discussion of the finite difference scheme and the implicit Newton–Raphson solver are given in the Supplemental Material S2 on the ASME Digital Collection. A mesh independence analysis has been performed to ensure the convergence of the numerical solution (see Supplemental Material S2). The total number of nodes along the radial direction is 220.

Having obtained the simulated temperature profiles, we next extract amplitude and phase from the simulation results. First, we choose a reference radius R0 that is close to the inner diameter of the light blocker but not irradiated by the concentrated light source. We also choose a region of analysis bounded by the reference radius R0 and an outer bound RN, in which significant amplitude decay is observed. Our code automatically detects quasi-steady-state oscillation when the maximum mean temperature rise within one cycle is less than 0.1% compared to the oscillation amplitude, as shown in Fig. 3(a). Amplitude decay and phase shift for all nodes are calculated with respect to the temperature profile at the reference radius. Representative amplitude decay and phase shift plots are shown in Figs. 3(b) and 3(c), respectively.

Fig. 3
(a) Simulated temperature profile at the reference radius and the boundary of the region of analysis, (b) amplitude decay of temperature profiles within the region of analysis, and (c) phase shift of temperature profiles within the region of analysis
Fig. 3
(a) Simulated temperature profile at the reference radius and the boundary of the region of analysis, (b) amplitude decay of temperature profiles within the region of analysis, and (c) phase shift of temperature profiles within the region of analysis
Close modal

3 Experimental Techniques

3.1 Experimental Setup.

In this work, the samples under test are thin disks (graphite: 0.92 mm, copper: 0.79 mm) with 88.50 mm outer diameter. Measurements were conducted under vacuum, and the vacuum level was approximately 0.05 Torr. The heat source is a 10 kW Xenon short arc lamp (Superior: SQP-SX10000FXT, Bethlehem, PA), and an electroformed parabolic reflector (Optiforms, Ag coated) concentrates the light. The inner diameter of the light blockers is 20.0 mm, and the distance between neighboring light blockers is 0.90 mm. A sample holder is designed to hold the light blockers and the sample concentrically as shown in Fig. 4(a). The sample is fixed to the sample holder using four pairs of small thermally insulating ceramic washers to minimize heat losses via conduction. The distance between the sample and the neighboring light blocker is 1.60 mm. We use three light blockers so that temperature oscillation on the last light blocker is minimal. The last light blocker facing the sample is coated with a thin layer of high-emissivity coating (Aremco 840CMX) to reduce radiative reflection from the sample and enable temperature measurements using IR thermography. We use an IR camera (Flir A655sc, calibrated to 2000 °C) for temperature measurements. A schematic of the system is shown in Fig. 4(b), and a photograph of the system is shown in Fig. 4(c).

Fig. 4
(a) Tungsten light blockers and the sample holder for measurements, (b) schematic of the system, and (c) a photograph of the system
Fig. 4
(a) Tungsten light blockers and the sample holder for measurements, (b) schematic of the system, and (c) a photograph of the system
Close modal

The concentrated light source is powered by a three-phase power supply with fixed 110 V with programmable current ranging from 100 to 200 A. The operating current is controlled externally using a 0–10 V function generator. Within the range of operation, the concentrated light source exhibits a linear relationship between external control voltage and the operating current. To heat the sample to high temperatures, we first control the solar simulator to produce steady heat flux. Upon reaching steady temperature (typically in 10 min), we modulate the function generator to output biased sinusoidal control voltage signals. Because of the linear relationship between external control voltage and the operating current, the concentrated light source also outputs biased sinusoidal heat flux. When the system reaches quasi-steady-state oscillation (approximately 1 min after introducing the oscillating heat flux), we start recording the temperature profiles using the IR camera.

3.2 Infrared Thermography and Processing.

To facilitate IR thermography, the back side of the sample is coated with a thin layer of high-emissivity coating (Aremco, 840MX). Our IR camera uses uncooled microbolometer detector with spectral range between 7.5 and 14 μm. From the spectral emissivity of the coating between 7.5 and 14 μm [25], we estimated the total emissivity in the IR camera's range is between 0.88 and 0.92. Because the coating's spectral emissivity is relatively independent of temperature, we measured its emissivity using a thermocouple (type K, ±0.75% accuracy) and the IR camera at relatively lower temperatures (370 °C), and the emissivity was 0.92. Therefore, we use 0.92 as the sample's emissivity for IR measurement. The temperature measurement error for the IR camera is approximately ±2% [26].

Because of symmetry, the isothermal lines in the IR image are concentric circles, as shown in Fig. 5(a). At steady-state, temperature measured for a pixel located at R0 exhibits much stronger noise in a high temperature environment as compared to a room temperature condition, as illustrated in Fig. 5(b). To reduce high random noise, we developed a computationally efficient method using bi-linear interpolation to obtain the average temperature at a given radius. The averaged temperature profiles at two representative locations R0 and RN are shown in Fig. 5(c). Details of the IR image processing are available in the Supplemental Material S4 on the ASME Digital Collection. We employ a parallel approach to average all the IR images in a recorded video, and the time required to fully process an IR recording consists of 500 images is less than one minute.

Fig. 5
(a) Temperature contour plot of a sample. The area exposed to the concentrated light source is indicated by the shaded region. (b) Steady temperature profile at pixel R0 at room and high temperatures. (c) averaged temperature profile at two representative radii R0 and RN.
Fig. 5
(a) Temperature contour plot of a sample. The area exposed to the concentrated light source is indicated by the shaded region. (b) Steady temperature profile at pixel R0 at room and high temperatures. (c) averaged temperature profile at two representative radii R0 and RN.
Close modal

Having obtained averaged temperature profiles at different radii, we extract amplitude decay and phase shift with respect to the reference radius R0. We employ a Fourier transform to calculate amplitude decay and phase shift for two sinusoidal temperature signals [27]. The main benefit of working with amplitude decay and phase shift instead of absolute temperature is the robustness against emissivity [18,19,28]. For the high-emissivity coating, a slight variation of the emissivity setting in the IR camera from 0.92 to 0.88 causes approximately 20 K difference in absolute temperature, as shown in Fig. 6(a). Because emissivity affects only absolute temperature and does not change the shape of the oscillating temperature profiles, intuitively, the phase of a sinusoidal temperature profile is not affected by the emissivity setting in the IR camera. Therefore, phase difference can be accurately measured even if emissivity is not known precisely, as shown in Fig. 6(b).

Fig. 6
(a) Temperature profiles measured at R0 and RN using different emissivity settings. (b) Phase difference calculated using temperature profiles between R0 and RN under different emissivity settings. (c) Amplitude decay calculated using temperature profiles between R0 and RN under different emissivity settings.
Fig. 6
(a) Temperature profiles measured at R0 and RN using different emissivity settings. (b) Phase difference calculated using temperature profiles between R0 and RN under different emissivity settings. (c) Amplitude decay calculated using temperature profiles between R0 and RN under different emissivity settings.
Close modal
For amplitude decay, we consider arbitrary locations Rx on the sample. If we neglect the effect of the surroundings, the radiation measured by the IR camera is proportional to the radiation emission at Rx. We denote the true emissivity at Rx as ε and true temperature as Tx, whereas the emissivity setting in the IR camera is ε and the corresponding measured temperature is Tx. The following equation assumes that the IR camera measures radiance correctly:
(14)
Because of oscillating heat flux, the true peak-to-valley temperature change is ΔT, and the temperature change observed in the IR camera is ΔT. Considering two locations on the sample R0 and Rx, and based on Eq. (14), we obtain the following relationships:
(15a)
The amplitude decay can be expressed as
(16)
where γ=(ε/ε)1/4. Based on Eq. (14), we also have
(17)
As a result, γ cancels, and the ratio of oscillation amplitude shows no dependence on emissivity. Fig. 6(c) shows that amplitude decay extracted from temperature profiles does not depend on emissivity
(18)

To summarize, the averaging of the IR image over the radius significantly reduces random noise in the high-temperature environment. The averaging approach is computationally efficient and requires less than 1 min to process an IR recording typically consisting of 500 IR images. We demonstrate that amplitude decay and phase shift extracted from measured temperature profiles are not affected by errors in emissivity, and therefore we work with amplitude decay and phase shift instead of absolute temperature profiles.

4 Model Sensitivity Analysis

In this study, we employ a full factorial design of experiments at two levels [29] to investigate the sensitivity of the physical model output (amplitude decay and phase shift) to the parameter of interest (thermal diffusivity). Here we use graphite's properties for demonstration. In addition to the parameter of interest, we are uncertain about several others. First, the intensity distribution of the concentrated light source is difficult to characterize precisely. Based on a previous study for a similar concentrated light source [21] and our preliminary investigation using a calorimetry method, the estimated intensity distribution width (σsolar) is 0.01 m. We assume a ±30% variation associated with σsolar. Second, in our study the light blocker temperature T0 is difficult to measure precisely. We assume a ±5% variation associated with T0 and assume the light blocker temperature is uniform and constant for sensitivity analysis. Finally, we treat thermal diffusivity as a temperature-independent parameter for sensitivity analysis. For a typical graphite sample the estimated thermal diffusivity α at 1000 K is 2.2×105m2/s. We assigned a relatively low range of variation (±5%) for α. Table 1 shows a summary of these parameters and the corresponding range of variations. We employ a main effect approach [29] to quantify model output sensitivity to an input parameter θ. For a two-level design, the main effect ME(θ) is expressed as
(19)

where f indicates the model, and θ+ and θ indicate the model at “high” and “low” parameter values. We evaluate the main effect of parameters defined in Table 1 to understand how model outputs (amplitude decay and phase shift) are influenced by input parameter changes. We choose an area (R0=60 pixels with 12 pixels span) that is not irradiated by the concentrated light for amplitude and phase calculations. This area is chosen because of high signal-to-noise ratio based on experimental observation. A prior study showed that increased heating frequency is beneficial to enhance sensitivity in thermal diffusivity [19].

Table 1

Parameters for a two level full factorial design of experiment

ParameterDefinitionRange
αThermal diffusivity(2.09, 2.31) ×105m2/s, 5% variation
σsolarLight source intensity distribution(0.7, 1.3)×102 m, 30% variation
T0Light blocker temperature(950, 1050) K, 5% variation
ParameterDefinitionRange
αThermal diffusivity(2.09, 2.31) ×105m2/s, 5% variation
σsolarLight source intensity distribution(0.7, 1.3)×102 m, 30% variation
T0Light blocker temperature(950, 1050) K, 5% variation

In this study, the concentrated light source is modulated sinusoidally to apply periodic heating to samples. We limit the periodic heating frequency to 0.15 Hz due to safety and operational stability of the Xenon light bulb [30]. The lower limit of the heating frequency is 0.02 Hz because of the lengthy measurement time at low heating frequencies. The main effect results for amplitude decay of selected input parameters are shown in Fig. 7 (top row). Similar to a prior study [19], sensitivity in thermal diffusivity increases with the heating frequency. For relatively fast heating frequencies (fheating0.08 Hz), emissivity and light blocker temperature exhibit low sensitivity. We also observe that the intensity distribution width σsolar contributes little to the main effect. This finding is advantageous compared to Hatta et al.'s approach [12], Cowan's approach [15] and the laser flash method [1], each of which are sensitive to the intensity distribution of the light source. The main effect results calculated using phase shift are shown in Fig. 7 (bottom row). Similarly, faster heating frequencies are advantageous for thermal diffusivity detection, and σsolar makes an insignificant contribution to the main effect.

Fig. 7
Top row: percentage main effect for amplitude ratio of different parameters at several heating frequencies. Bottom row: percentage main effect for phase difference of different parameters at several heating frequencies.
Fig. 7
Top row: percentage main effect for amplitude ratio of different parameters at several heating frequencies. Bottom row: percentage main effect for phase difference of different parameters at several heating frequencies.
Close modal

The foregoing results indicate that high sensitivity for thermal diffusivity can be achieved with fheating0.08 Hz for a graphite sample. The results also indicate that amplitude and phase are insensitive to the intensity distribution of the light source σsolar.

5 Uncertainty Quantification and Results

5.1 Uncertainty Quantification Using a Bayesian Framework.

In this study, we estimate thermal diffusivity using amplitude decay and phase shift measurements. Conventional parameter estimation methods convert an inverse problem to an optimization problem by minimizing residuals between measurement results and the physical model. In this case, the inverse problem is ill-posed, and the solution is not unique. Also such methods lack rigorous uncertainty estimates for parameters. This study employs a Bayesian framework that incorporates prior knowledge of parameters and measurement processes to produce a complete statistical description of the parameters. Growing interest exists in applying the Bayesian framework to inverse heat transfer problems, and several studies have demonstrated its benefits [19,3133].

The Bayesian framework is based on the Bayes' theorem. Consider a set of parameters θ and measurement results D; the posterior distribution for θ after observing D is expressed using Bayes' theorem
(20)

In this work, D consists of amplitude, phase measurements within a region of analysis, and average temperature profiles at inner (R0) and outer bounds (RN) of the region. θ consists of four parameters: Aα,Bα, σsolar, and Tbias. Aα and Bα are the coefficients for temperature-dependent thermal diffusivity and are defined in Eq. (2). In our study, the intensity distribution of the concentrated light source is difficult to characterize precisely. Therefore, we treat σsolar as a unknown parameter and update our state of knowledge for σsolar using the Bayesian framework. Finally, the light blocker temperature measured using the IR camera (see Supplemental Material S3 on the ASME Digital Collection for details) also suffers slight system error because the measurement process differs from real experiments. Therefore, we introduce another parameter Tbias to augment the measured light blocker temperature to model the light blocker temperature under real measurement conditions.

The posterior distribution is proportional to the likelihood function p(D|θ) and the prior distribution p(θ). The likelihood function p(D|θ) models the measurement process and indicates the probability of observing data D given parameter θ and measurement noise σD. First we consider amplitude decay measurements. We denote measured amplitude decay at radius Ri as ΔAi and assume that it follows a Gaussian distribution with mean A[f(θ)] and standard deviation σΔA,i. A[f(θ)] indicates the amplitude decay from the physical model f evaluated using parameter θ. For amplitude decay at Ri, the likelihood function is given by
(21)
Assuming amplitude measurements at different radii are independent, then for all amplitude decay measurements ΔA within the region of analysis, the likelihood function is
(22)
The likelihood functions for all phase shift measurements Δϕ within the region of analysis are derived similarly. ϕ[f(θ)] indicates the phase shift from the physical model f evaluated using parameter θ
(23)
In addition, the likelihood functions for the mean temperature profiles at the inner (Tmean,R0) and outer (Tmean,RN) bounds of the region of analysis are treated as follows. Tm([f(θ)]) indicates the mean temperature obtained from the physical model f evaluated using parameter θ. σT indicates the standard deviation in mean temperature measurements. Unlike amplitude and phase measurements, we manually choose σT = 2 K to ensure the mean simulated temperatures within the region of analysis match experimental measurements. Because we use temperature-dependent properties in this study, matching simulated and experimentally measured mean temperatures is important. The likelihood function becomes
(24)
By further assuming that ΔA,Δϕ,Tmean,R0 and Tmean,RN are statistically independent, the simplified likelihood function becomes
(25)
The prior distribution p(θ) represents the prior knowledge for θ. We assume Aα,Bα, σsolar, and Tbias are statistically independent. Therefore, p(θ) is expressed as
(26)
Substituting the likelihood function (Eq. (25)) and the prior distribution (Eq. (26)) into the Bayes posterior defined in Eq. (20), the posterior distribution is expressed as
(27)

For an arbitrary distribution with unknown normalization constant, the Markov chain Monte Carlo (MCMC) method is typically used to obtain random samples to reconstruct such a distribution. The most commonly used MCMC algorithm is random walk Metropolis [34] because of its simplicity [19,31,32]. However, this algorithm suffers poor convergence and low efficiency for problems with large parameter spaces [35]. To expand the exploration of the parameter space and to improve convergence of the Markov chain, our work employs a no-U-turn (NUT) sampler [36]. In addition, MCMC requires frequent evaluations of the physical model to calculate the likelihood function. This becomes infeasible for numerical models that require lengthy execution times. A previous study [32] demonstrated significant time efficiency improvement using polynomial chaos for an inverse heat transfer problem. In this study, we develop a nonintrusive polynomial chaos (order 4) surrogate model using a python module “chaospy” [37] to accelerate the original finite difference model. The surrogate model exhibits high accuracy (maximum 0.5% error) for steady temperature profiles, amplitude decay, and phase shift predictions, yet the execution time is reduced to a few milliseconds from several minutes. We employ the python module pymc3 [38] with custom modifications [39] to implement the NUT sampler.

5.2 Nut Sampler's Results Using a Graphite Sample.

To demonstrate the Bayesian framework and NUT sampler, we first tested an isotropic graphite sample (manufacturer:Entegris, grade:TM). The sample's thickness is 0.92 mm with outer diameter 88.50 mm. The density of the graphite sample is 1740 kg/m3. We used a polynomial fit for specific heat [40] as
(28)

The gray emissivity of the graphite is 0.77, and we coated a layer of thin layer of diffuse high-emissivity (0.88) material on the back side of the sample to facilitate IR thermography. We modulated the light source using a 0.1 Hz voltage signal to achieve relatively high sensitivity in thermal diffusivity. We heated the sample from room temperature to high temperature with 119 A solar simulator current. Steady temperatures were obtained after 10 min, and then a 0.1 Hz oscillation current signal with amplitude 22 A was superimposed on the 110 A bias current. We started the IR recording after approximately one to two minutes when the sample exhibited quasi-steady-state oscillations.

Here, we analyze a region on the IR image bounded by R0=61 and RN = 73 pixels to obtain temperature profiles and corresponding amplitude decay and phase shift. This region is chosen for high signal-to-noise ratio in both amplitude ratio and phase shift. In addition, the region is chosen to be relatively narrow to keep the temperature difference low (approximately 40 K), such that Eq. (2) remains valid to model thermal diffusivity as a function of temperature.

The results obtained using the NUT sampler are shown in Fig. 8. The top row of subfigures contains the trace plots for four unknown parameters. A trace plot shows the sequential value of a parameter obtained from the NUT sampler. The trace plots exhibit good mixing, and different chains for each parameter oscillate around nearly the same constant value, indicating convergence of the sampling process. The middle row contains histograms obtained from the trace plots of corresponding parameters. The normalized histograms represent the posterior distribution of parameters and indicate the updated state of knowledge given the prior knowledge and the measurement process. The bottom row contains the autocorrelation of the corresponding trace plots in the top row. The autocorrelation decreases rapidly to zero for all parameters, indicating high sampling efficiency.

Fig. 8
Top row: trace plots of parameters. (a1) Intensity distribution width for the concentrated light source σsolar. (b1) and (c1) Parameters for temperature-dependent thermal diffusivity Aα and Bα. (d1) Bias temperature of the light blocker Tbias. Middle row: histograms of the trace plots in the top row. Bottom row: autocorrelation of the trace plots in the top row.
Fig. 8
Top row: trace plots of parameters. (a1) Intensity distribution width for the concentrated light source σsolar. (b1) and (c1) Parameters for temperature-dependent thermal diffusivity Aα and Bα. (d1) Bias temperature of the light blocker Tbias. Middle row: histograms of the trace plots in the top row. Bottom row: autocorrelation of the trace plots in the top row.
Close modal

Given the posterior distributions of Aα and Bα, the temperature-dependent thermal diffusivity within the region of analysis (ROA) is shown in Fig. 9. The mean temperatures of the inner and the outer boundary of the ROA are approximately 960 K and 1000 K, respectively. The solid line indicates the posterior mean for thermal diffusivity, and the blue band indicates corresponding uncertainty. The same grade graphite samples have been measured using the laser flash method by the manufacturer [40] and are shown as green lines. Good agreement between the posterior results and reference values [40] is observed for the entire temperature range within the ROA.

Fig. 9
Thermal diffusivity as a function of temperature evaluated using the posterior mean of Aα and Bα (solid line) and corresponding uncertainties (shade region). The reference values are shown as the solid dots with error bars [40].
Fig. 9
Thermal diffusivity as a function of temperature evaluated using the posterior mean of Aα and Bα (solid line) and corresponding uncertainties (shade region). The reference values are shown as the solid dots with error bars [40].
Close modal

The fitted amplitude (Fig. 10(a)), phase (Fig. 10(b)), and temperatures (Fig. 10(c)) are plotted using the parameters' posterior means, and they match well with experimental measurements.

Fig. 10
Fitted (solid black lines) amplitude decay (a), phase shift (b) and temperature profiles at inner and outer boundary of the ROA (c) using the posterior means of parameters from Fig. 8(a2)–8(d2)
Fig. 10
Fitted (solid black lines) amplitude decay (a), phase shift (b) and temperature profiles at inner and outer boundary of the ROA (c) using the posterior means of parameters from Fig. 8(a2)–8(d2)
Close modal

5.3 Results for Graphite at Different Frequencies and Bias Voltages.

We fixed the bias current of the solar simulator at 119 A and measured multiple graphite samples at different frequencies. Here, the tested frequency range is from 0.08 Hz to 0.14 Hz to achieve relatively high sensitivity in thermal diffusivity. We analyze amplitude decay and phase shift within R0=61 and RN = 73 pixels, and the results are shown in Fig. 11 in which thermal diffusivity values exhibit good consistency among the different frequencies and multiple samples.

Fig. 11
Thermal diffusivity measured at different frequencies for multiple graphite samples compared to reference values [40]. Note sample 1 exhibits a different temperature range because it was tested with a different light bulb and slightly different reflector condition.
Fig. 11
Thermal diffusivity measured at different frequencies for multiple graphite samples compared to reference values [40]. Note sample 1 exhibits a different temperature range because it was tested with a different light bulb and slightly different reflector condition.
Close modal

Next, we fixed heating frequency at 0.1 Hz and changed the bias current. Here, we choose 119, 152, and 179 A as bias current for the solar simulator and 22 A as the oscillation amplitude. We also explore different region of analysis on the sample. In addition to the region bounded by R0=61 and RN = 73 pixels used in Secs. 5.2 and 5.3), we analyze another region bounded by R0=73 and RN = 85 pixels shown as region 2 in Fig. 12(a). This region exhibits lower signal-to-noise ratio compared to region 1, but it also experiences lower temperatures and can extend the temperature measurement range for each test. The results for different currents and different regions are shown in Fig. 12(b). We achieve thermal diffusivity measurements over a broader temperature range (220 K) using a wider range of bias currents. The results obtained under different bias currents exhibit good agreement with the reference values [40] over 900–1120 K.

Fig. 12
(a) Region of analysis in an IR image. Region 1 (inner shaded area) bounded by R0=61 and RN = 73 pixels. Region 2 (outer shaded area) bounded by R0=73 and RN = 85 pixels. (b) Thermal diffusivity measured under different bias currents and in different regions, compared to reference values [40].
Fig. 12
(a) Region of analysis in an IR image. Region 1 (inner shaded area) bounded by R0=61 and RN = 73 pixels. Region 2 (outer shaded area) bounded by R0=73 and RN = 85 pixels. (b) Thermal diffusivity measured under different bias currents and in different regions, compared to reference values [40].
Close modal

5.4 Results for a Copper Sample.

We also tested a copper sample (certified copper 110, 99.9% purity). The sample's thickness is 0.79 mm with outer diameter 88.50 mm. We use a polynomial fit for temperature-dependent specific heat [41]. The copper sample's front and back surfaces were coated with a thin layer of diffuse high-emissivity (0.88) coating (Aremeco, 840MX). The weight contribution of the coating is less than 2% of the sample's total weight and is neglected. The copper sample was measured using 119 A bias current under several different frequencies. The results are shown in Fig. 13, and good agreement with reference values [42] is observed.

Fig. 13
Thermal diffusivity of copper measured at different heating frequencies with 119 A bias current
Fig. 13
Thermal diffusivity of copper measured at different heating frequencies with 119 A bias current
Close modal

5.5 Limitations of the Custom Instrument.

This work reports a custom Ångström's method for rapid and accurate high-temperature thermal diffusivity characterization. In this section, several limitations are presented and should be addressed in future work.

First, the custom instrument and methods reported here are suitable for materials with high diffusivity (e.g., above approximately 1×105 m2/s) and thermal properties that are isotropic in the in-plane directions. The concentrated light intensity causes significant temperature gradients and thermal stresses on materials with low diffusivity. As a result, such materials tend to bend or crack during testing. An image for a titanium sample after testing and corresponding measurement results are shown in the Supplemental Material S6 on the ASME Digital Collection.

Second, compared to past Ångström's methods, the custom instrument requires the quantification of energy input to the sample and energy lost to the surroundings because a temperature gradient exists along the heat conduction direction, and temperature-dependent properties must be used in the model. In addition, an analytical solution does not exist, and the modeling process requires numerical computation, which can be time consuming. However, we developed an accurate and rapid surrogate model to mitigate this issue under the fin approximation, and future work could develop two-dimensional models for thicker samples.

Third, samples under test typically require a high-emissivity coating layer to obtain reliable amplitude decay and phase shift results. Such coatings can deteriorate when temperature increases beyond 1400 °C [25], thus limiting the temperature range for some materials. Damaged coatings can become difficult to remove and make repeated measurement on a single sample difficult.

Fourth, because the focal spot of the concentrated light source is similar in size to the sample's radius, we employed several light blockers so that only the center of the sample is exposed. This approach increases the complexity of the model and experimental setup. A potential improvement is to use an energy source with smaller focal spot, such as a powerful laser source to obviate the need for light blockers.

Fifth, the concentrated light source in this study is not designed to sustain prolonged power fluctuations [30]. We observed relatively fast degradation of bulb life span under oscillating conditions (approximately 100 h.) compared to its nominal life span (approximately 800 h.). Such performance degradation may lead to variations in testing results.

6 Conclusions

In this work, we present a custom instrument for characterizing thermal diffusivity at high temperatures using a modified Ångström's method. The sample under test is heated directly using a concentrated light source and reaches steady high temperatures rapidly (approximately 10 min). This is advantageous compared to existing methods that heat samples from surroundings and typically require hours to reach steady high temperatures. We employ noncontact, nondestructive and data-rich IR thermography to measure quasi-steady-state oscillating temperature profiles, and an efficient averaging method to process the noisy IR images and obtain temperature profiles with low noise. Amplitude and phase results extracted from such temperature profiles are robust against errors in emissivity. For uncertainty quantification, a surrogate-accelerated Bayesian framework and a NUT sampler are developed to obtain statistical distributions for parameters of interest. For several graphite and copper samples under test, we observe good agreement with reference values under different heating frequencies. We also measure different temperature ranges by adjusting solar simulator current. In general, the results obtained using the custom instrument exhibit approximately 5% difference compared to reference values obtained using the laser flash method. The custom instrument's time efficiency, accuracy, inherent ability to quantify uncertainty rigorously, and relatively low overall cost are expected to benefit researchers for characterizing thermal diffusivity at high temperatures.

Acknowledgment

The authors are grateful for Professor Ilias Bilionis at Purdue University for his free online uncertainty quantification course. The authors thank Professor Jeff Eldredge at UCLA for improving the computational efficiency of the numerical model. The authors thank Dr. Michael Barako from Northrop Grumman and Dr. James Wilson from Raytheon Technologies for fruitful feedback on the project. The authors appreciate the help from Ben Alpert of Aremco Products Inc. for high-temperature coating. The authors are grateful for the support from MAE staff Ben Tam, Miguel Lozano, Amanda Gordillo, and Marla Cooper during the challenging covid period. This work was supported by the Center for Integrated Thermal Management of Aerospace Vehicles and a UCLA Dissertation Year Fellowship.

References

1.
Parker
,
W.
,
Jenkins
,
R.
,
Butler
,
C.
, and
Abbott
,
G.
,
1961
, “
Flash Method of Determining Thermal Diffusivity, Heat Capacity, and Thermal Conductivity
,”
J. Appl. Phys.
,
32
(
9
), pp.
1679
1684
.10.1063/1.1728417
2.
Clark Iii
,
L.
, and
Taylor
,
R.
,
1975
, “
Radiation Loss in the Flash Method for Thermal Diffusivity
,”
J. Appl. Phys.
,
46
(
2
), pp.
714
719
.10.1063/1.321635
3.
Chu
,
F.
,
Taylor
,
R.
, and
Donaldson
,
A.
,
1980
, “
Thermal Diffusivity Measurements at High Temperatures by the Radial Flash Method
,”
J. Appl. Phys.
,
51
(
1
), pp.
336
341
.10.1063/1.327377
4.
Cezairliyan
,
A.
,
1984
, “
A Dynamic Technique for Measurements of Thermophysical Properties at High Temperatures
,”
Int. J. Thermophys.
,
5
(
2
), pp.
177
193
.10.1007/BF00505499
5.
Baba
,
T.
, and
Cezairliyan
,
A.
,
1994
, “
Thermal Diffusivity of POCO AXM-5Q1 Graphite in the Range 1500 to 2500 K Measured by a Laser-Pulse Technique
,”
Int. J. Thermophys.
,
15
(
2
), pp.
343
364
.10.1007/BF01441590
6.
Gustafsson
,
S.
,
1991
, “
Transient Plane Source Techniques for Thermal Conductivity and Thermal Diffusivity Measurements of Solid Materials
,”
Rev. Sci. Instrum.
,
62
(
3
), pp.
797
804
.10.1063/1.1142087
7.
Healy
,
J.
,
De Groot
,
J.
, and
Kestin
,
J.
,
1976
, “
The Theory of the Transient Hot-Wire Method for Measuring Thermal Conductivity
,”
Phys. B+C
,
82
(
2
), pp.
392
408
.10.1016/0378-4363(76)90203-5
8.
Angström
,
A.
,
1863
, “
XVII. New Method of Determining the Thermal Conductibility of Bodies
,”
Philos. Mag.
,
25
(
166
), pp.
130
142
.10.1080/14786446308643429
9.
Sidles
,
P.
, and
Danielson
,
G.
,
1954
, “
Thermal Diffusivity of Metals at High Temperatures
,”
J. Appl. Phys.
,
25
(
1
), pp.
58
66
.10.1063/1.1721521
10.
Katsura
,
T.
,
1993
, “
Thermal Diffusivity of Silica Glass at Pressures Up to 9 GPa
,”
Phys. Chem. Miner.
,
20
(
3
), pp.
201
208,
.10.1007/BF00200122
11.
Xu
,
Y.
,
Shankland
,
T.
,
Linhardt
,
A.
,
Rubie
,
D.
,
Langenhorst
,
F.
, and
Klasinski
,
K.
,
2004
, “
Thermal Diffusivity and Conductivity of Olivine, Wadsleyite and Ringwoodite to 20 GPa and 1373 K
,”
Phys. Earth Planet. Inter.
,
143–144
, pp.
321
336
.10.1016/j.pepi.2004.03.005
12.
Hatta
,
I.
,
Sasuga
,
Y.
,
Kato
,
R.
, and
Maesono
,
A.
,
1985
, “
Thermal Diffusivity Measurement of thin Films by Means of an AC Calorimetric Method
,”
Rev. Sci. Instrum.
,
56
(
8
), pp.
1643
1647
.10.1063/1.1138117
13.
Abeles
,
B.
,
Cody
,
G.
, and
Beers
,
D.
,
1960
, “
Apparatus for the Measurement of the Thermal Diffusivity of Solids at High Temperatures
,”
J. Appl. Phys.
,
31
(
9
), pp.
1585
1592
.10.1063/1.1735897
14.
Vandersande
,
J.
, and
Pohl
,
R.
,
1980
, “
Simple Apparatus for the Measurement of Thermal Diffusivity Between 80–500 K Using the Modified Ångström Method
,”
Rev. Sci. Instrum.
,
51
(
12
), pp.
1694
1699
.10.1063/1.1136158
15.
Cowan
,
R.
,
1961
, “
Proposed Method of Measuring Thermal Diffusivity at High Temperatures
,”
J. Appl. Phys.
,
32
(
7
), pp.
1363
1370
.10.1063/1.1736235
16.
Wheeler
,
M.
,
1965
, “
Thermal Diffusivity at Incandescent Temperatures by a Modulated Electron Beam Technique
,”
Br. J. Appl. Phys.
,
16
(
3
), pp.
365
376
.10.1088/0508-3443/16/3/311
17.
Tanaka
,
T.
, and
Suzuki
,
H.
,
1972
, “
The Thermal Diffusivity of Pyrolytic Graphite at High Temperatures
,”
Carbon
,
10
(
3
), pp.
253
257
.10.1016/0008-6223(72)90323-5
18.
Hahn
,
J.
,
Reid
,
T.
, and
Marconnet
,
A.
,
2019
, “
Infrared Microscopy Enhanced Angstrom's Method for Thermal Diffusivity of Polymer Monofilaments and Films
,”
ASME J. Heat Transfer-Trans. ASME
,
141
(
8
), p.
081601
.10.1115/1.4043619
19.
Hu
,
Y.
, and
Fisher
,
T.
,
2020
, “
Accurate Thermal Diffusivity Measurements Using a Modified Ångström's Method With Bayesian Statistics
,”
ASME J. Heat Transfer-Trans. ASME
,
142
(
7
), p.
071401
.10.1115/1.4047145
20.
Gomez-Garcia
,
F.
,
Santiago
,
S.
,
Luque
,
S.
,
Romero
,
M.
, and
Gonzalez-Aguilar
,
J.
,
2016
, “
A New Laboratory-Scale Experimental Facility for Detailed Aerothermal Characterizations of Volumetric Absorbers
,”
AIP Conf. Proc.,
1734
, p.
030018
.10.1063/1.4949070
21.
Abuseada
,
M.
,
Ophoff
,
C.
, and
Ozalp
,
N.
,
2019
, “
Characterization of a New 10 kWe High Flux Solar Simulator Via Indirect Radiation Mapping Technique
,”
ASME J. Sol. Energy Eng.
,
141
(
2),
p. 021005.10.1115/1.4042246
22.
Seipold
,
U.
,
1998
, “
Temperature Dependence of Thermal Transport Properties of Crystalline Rocks-a General Law
,”
Tectonophysics
,
291
(
1–4
), pp.
161
171
.10.1016/S0040-1951(98)00037-7
23.
Miao
,
S.
,
Li
,
H.
, and
Chen
,
G.
,
2014
, “
Temperature Dependence of Thermal Diffusivity, Specific Heat Capacity, and Thermal Conductivity for Several Types of Rocks
,”
J. Therm. Anal. Calorim.
,
115
(
2
), pp.
1057
1063
.10.1007/s10973-013-3427-2
24.
Abdulagatova
,
Z.
,
Kallaev
,
S.
,
Omarov
,
Z.
,
Bakmaev
,
A.
,
Grigor'ev
,
B.
, and
Abdulagatov
,
I.
,
2020
, “
Temperature Effect on Thermal-Diffusivity and Heat-Capacity and Derived Values of Thermal-Conductivity of Reservoir Rock Materials
,”
Geomech. Geophys. Geo-Energy
,
6
(
1
), pp.
1
23
.10.1007/s40948-019-00131-2
25.
Aremco Products,
2020
, High Temperature High Emissivity Coatings,
Aremco Product
Valley Cottage, NY, accessed Feb. 2, https://www.aremco.com/high-emissivity-coatings/
26.
FLIR Systems,
2017
, User's Manual FLIR A6xx Series,
FLIR Systems,
Wilsonville, OR, accessed Feb. 2, 2020, https://www.flir.com/products/a655sc/
27.
Hu
,
Y.
,
Brahim
,
S.
,
Maat
,
S.
,
Davies
,
P.
,
Kundu
,
A.
, and
Fisher
,
T.
,
2020
, “
Rapid Analytical Instrumentation for Electrochemical Impedance Spectroscopy Measurements
,”
J. Electrochem. Soc.
,
167
(
2
), p.
027545
.10.1149/1945-7111/ab6ee8
28.
Visser
,
E.
,
Versteegen
,
E.
, and
van Enckevort
,
W.
,
1992
, “
Measurement of Thermal Diffusion in Thin Films Using a Modulated Laser Technique: Application to Chemical-Vapor-Deposited Diamond Films
,”
J. Appl. Phys.
,
71
(
7
), pp.
3238
3248
.10.1063/1.350970
29.
Montgomery
,
D.
,
2017
,
Design and Analysis of Experiments
,
Wiley
, Hoboken, NJ.
30.
Superior Quartz Products,
2020
, Lamp Characteristics and Description,
Superior Quartz Products
, Bethlehem, PA, accessed Feb. 2, https://www.sqpuv.com/PDFs/TechnicalSpecificationGuide.pdf
31.
Wang
,
J.
, and
Zabaras
,
N.
,
2004
, “
A Bayesian Inference Approach to the Inverse Heat Conduction Problem
,”
Int. J. Heat Mass Transfer
,
47
(
17–18
), pp.
3927
3941
.10.1016/j.ijheatmasstransfer.2004.02.028
32.
Rynn
,
J.
,
Cotter
,
S.
,
Powell
,
C.
, and
Wright
,
L.
,
2019
, “
Surrogate Accelerated Bayesian Inversion for the Determination of the Thermal Diffusivity of a Material
,”
Metrologia
,
56
(
1
), p.
015018
.10.1088/1681-7575/aaf984
33.
Xu
,
X.
,
Zhang
,
Q.
,
Hao
,
M.
,
Hu
,
Y.
,
Lin
,
Z.
,
Peng
,
L.
,
Wang
,
T.
,
Ren
,
X.
,
Wang
,
C.
,
Zhao
,
Z.
,
Wan
,
C.
,
Fei
,
H.
,
Wang
,
L.
,
Zhu
,
J.
,
Sun
,
H.
,
Chen
,
W.
,
Du
,
T.
,
Deng
,
B.
,
Cheng
,
G. J.
,
Shakir
,
I.
,
Dames
,
C.
,
Fisher
,
T. S.
,
Zhang
,
X.
,
Li
,
H.
,
Huang
,
Y.
, and
Duan
,
X.
,
2019
, “
Double-Negative-Index Ceramic Aerogels for Thermal Superinsulation
,”
Science
,
363
(
6428
), pp.
723
727
.10.1126/science.aav7304
34.
Metropolis
,
N.
,
Rosenbluth
,
A.
,
Rosenbluth
,
M.
,
Teller
,
A.
, and
Teller
,
E.
,
1953
, “
Equation of State Calculations by Fast Computing Machines
,”
J. Chem. Phys.
,
21
(
6
), pp.
1087
1092
.10.1063/1.1699114
35.
Brooks
,
S.
,
Gelman
,
A.
,
Jones
,
G.
, and
Meng
,
X.
,
2011
,
Handbook of Markov Chain Monte Carlo
,
CRC Press
Boca Raton, FL.
36.
Hoffman
,
M.
, and
Gelman
,
A.
,
2014
, “
The No-U-Turn Sampler: Adaptively Setting Path Lengths in Hamiltonian Monte Carlo
,”
J. Mach. Learn. Res.
,
15
(
1
), pp.
1593
1623
.https://arxiv.org/abs/1111.4246
37.
Feinberg
,
J.
, and
Langtangen
,
P.
,
2015
, “
Chaospy: An Open Source Tool for Designing Methods of Uncertainty Quantification
,”
J. Comput. Sci.
,
11
, pp.
46
57
.10.1016/j.jocs.2015.08.008
38.
Salvatier
,
J.
,
Wiecki
,
T.
, and
Fonnesbeck
,
C.
,
2016
, “
Probabilistic Programming in Python Using PyMC3
,”
PeerJ Comput. Sci.
,
2
, p.
e55
.10.7717/peerj-cs.55
39.
Matt Pitkin
,
2018
, “
Using a Black Box Likelihood Function
,” Matt Pitkin, accessed Feb. 2, 2020, http://mattpitkin.github.io/samplers-demo/pages/pymc3-blackbox-likelihood/
40.
Entegris,
2020
, Graphite Properties and Characteristics,
Entegris
, Bellerica, MA, accessed June 19, https://poco.entegris.com/content/dam/poco/resources/reference-materials/brochures/brochure-graphite-properties-and-characteristics-11043.pdf
41.
White
,
G.
, and
Collocott
,
S.
,
1984
, “
Heat Capacity of Reference Materials: Cu and W
,”
J. Phys. Chem. Ref. Data
,
13
(
4
), pp.
1251
1257
.10.1063/1.555728
42.
Touloukian
,
Y. S.
,
Powell
,
R. W.
,
Ho
,
C. Y.
, and
Nicolaou
,
M. C.
,
1974
,
Thermophysical Properties of Matter, The TPRC Data Series. Volume 10. Thermal Diffusivity. Data Book
, West Lafayette, IN, Vol.
10
, p.
1
.

Supplementary data