Abstract

This report documents evaluation of simultaneous estimation of multiple interfacial heat transfer coefficients (HTCs) using transient measurements from an experiment designed for steady-state operation. The design of a mirror system for directing X-rays under cryogenic conditions requires knowledge of the interfacial HTC (contact conductance) between silicon and indium. An experimental apparatus was constructed to measure temperatures in a stack of five 7.62-mm thick pucks of silicon separated by 0.1-mm thick sheets of indium, which is operated under cryogenic temperatures in vacuum. Multiple pucks and interfaces are incorporated into the apparatus to allow evaluation of HTCs for surfaces of different roughness from a single experiment. Analysis of the sensitivity of each of the measured temperatures to each of the unknown HTCs reveals lack of linear independence of these sensitivities and suggests the recovery of the HTCs will be challenging. Artificially noised “data” were created from two different computational models by solving for temperatures and adding random Gaussian noise with a specified standard deviation. These data are subsequently analyzed using two different iterative parameter estimation methods: a Levenberg scheme and a Tikhonov iterative scheme. The required sensitivity matrix is computed using forward finite difference approximations. The results for the heat transfer coefficients for this model problem suggest that coefficients cannot be estimated independently, but the ratios relative to one of the unknowns can be recovered.

Introduction

Synchrotron light sources are used for research in a variety of fields including biology, medicine, condensed matter physics, and materials science. Many experiments using synchrotron light involve studying the structure of matter on the subnanometer to millimeter scale, for example, determining the structure of a protein molecule [1,2]. In a synchrotron light source, electrons traveling at relativistic speeds emit light. This light, typically of X-ray wavelengths, is directed to experimental stations or “end stations” by beamlines. A typical beamline may have one or more mirrors to direct and focus the X-rays.

A continuing challenge in the design of beamline mirrors is thermo-mechanical distortion. This distortion is caused by the energy absorbed by the mirror from the X-ray beam. Typical power loads are of the order of 400 W with a peak power density of 1 W/mm2. One strategy to limit this distortion is to cool a silicon mirror substrate to cryogenic temperatures [35]. This strategy is based on the material properties of silicon: near 125 K, the coefficient of thermal expansion of silicon is approximately zero, and the thermal conductivity is approximately four times greater than at room temperature. In practice, this type of mirror consists of a block of silicon clamped to a cooling block through which the cryogenic coolant flows. A collaboration of researchers from five United States Department of Energy light sources is currently developing a mirror cooled by liquid argon.

To design a mirror which operates near this target temperature of 125 K, accurate thermal models are required. One of the key parameters in these models is the thermal contact conductance of the interface between the mirror and its cooling block. In the liquid-argon-cooled design currently being developed, this interface consists of a silicon mirror substrate, a layer of indium foil, and a silicon cooling block (Fig. 1).

Fig. 1
The silicon mirror substrate is cooled by liquid argon via a silicon cooling block and a layer of indium foil. (Image courtesy Andrew Ringwall, Plamen Velikov, Thomas Rabedeau, SLAC nation accelerator laboratory).
Fig. 1
The silicon mirror substrate is cooled by liquid argon via a silicon cooling block and a layer of indium foil. (Image courtesy Andrew Ringwall, Plamen Velikov, Thomas Rabedeau, SLAC nation accelerator laboratory).
Close modal

Other authors have measured the thermal contact conductance across interfaces for high heat load X-ray optics using steady-state measurements made at room temperature. Khounsary and Takeuchi found thermal contact conductance values as a function of applied interfacial pressure for a silicon-indium-copper interface [6]. Khounsary also measured contact conductance of a silicon-silver-copper interface [7]. Asano measured thermal contact conductance values across a silicon carbide-gallium-copper interface and a silicon carbide-indium-gallium-copper interface [8]. Marion measured thermal contact conductance of copper–tungsten, copper–copper, and copper–silicon under different contact pressures and with and without interstitial materials gold and indium foil [9].

While calculating the thermal contact conductance from steady-state measurements is straightforward, the time required to reach steady-state conditions can be significant, and therefore in practice can limit the number of measurements that can be taken in a given amount of time. A potentially more efficient method to determine the thermal contact conductance is to use measurements made during the transient response.

Using transient data, however, poses other challenges. When using transient data, the heat flux is not directly measurable and has to be calculated using inverse methods [10]. In Ref. [10], Fieberg and Kneer use thermographic imaging to capture the transient response of two heated specimens that are brought into contact under pressure. Using Beck's sequential estimation method [11], they solved the inverse problem to determine the heat flux across the interface as a function of time. With knowledge of the temperature distribution from the thermographic imagery, they computed the HTC at the interface.

Analysis to find the HTCs in the experiment under consideration in the current investigation is a parameter estimation problem, as compared to a function estimation problem. In parameter estimation a “small” number of unknowns are sought, whereas in function estimation, the history or variation over time of an unknown stimulus needs to be found.

A nonlinear estimation method was developed by Beck to calculate contact conductance of thermally thick specimens as a constant or function of time; he also presented criteria to determine the optimal experimental approach for accurately estimating contact conductance from transient data [12]. Woodbury et al. describe the use of nonlinear estimation methods in solving transient inverse heat conduction problems with multiple unknown heat fluxes [11]. As mentioned previously, nonlinear estimation, along with infrared thermography, was used by Fieberg and Kneer to find contact conductance values of steel and steel–aluminum alloy combinations at high contact pressures [10]. Bi et al. utilized a different thermal imaging method, Laser Photothermal Method, to find contact conductance from transient data at cryogenic temperatures as a function of contact pressure for stainless steel [13].

While thermal contact conductance has been measured from (1) joints that include silicon and indium using steady-state data at room temperature and (2) other materials using transient data, this study addresses the unique case of using transient temperature data to find the thermal contact conductance across a silicon-indium-silicon interface at cryogenic temperatures as a function of interface pressure and surface roughness.

Experimental Setup

The apparatus shown in Fig. 2 was constructed to measure the thermal contact conductance of silicon-indium-silicon interfaces at cryogenic temperatures. This apparatus allows application of a known interface pressure and heat load to a sample while measuring the sample temperature. The sample is compressed between a cold plate and a hot plate. The cold plate is made of copper and is cooled by forced convection by a flow of liquid nitrogen. The hot plate is also made of copper and contains an electric resistance cartridge heater connected to a variable power supply. The hot plate is mounted on a spring-loaded lead screw mechanism through a load cell, which allows variation and measurement of the sample interface pressure. To minimize heat transfer between the sample and the lab environment, the sample, lead screw, load cell, cold plate, and heater are all enclosed in a vacuum chamber. Additional thermal isolation is provided by PEEK spacers between the hot plate and load cell and the cold plate and vacuum chamber base. The sample temperature is measured using calibrated silicon diode temperature sensors.

Fig. 2
The experimental apparatus consists of a sample compressed between a hot and cold plate, and enclosed in a vacuum chamber. The sample is a stack of five silicon pucks separated by layers of indium foil.
Fig. 2
The experimental apparatus consists of a sample compressed between a hot and cold plate, and enclosed in a vacuum chamber. The sample is a stack of five silicon pucks separated by layers of indium foil.
Close modal

The sample is a stack of five approximately cylindrical silicon “pucks” in which each puck is separated from its neighbor by a layer of indium foil as shown in Fig. 3(a). Each puck has a nominal diameter of 25.4 mm and thickness of 7.62 mm. At each of the four silicon-indium-silicon interfaces the indium foil is nominally 0.1-mm thick. Additionally, there are also interfaces between the silicon pucks and the copper hot and cold surfaces. This stacked sample configuration is selected to allow simultaneous measurement and comparison of the thermal contact conductance of multiple nominally identical interfaces while under the same contact pressure.

Fig. 3
Experimental apparatus (a) and nomenclature for analysis (b)
Fig. 3
Experimental apparatus (a) and nomenclature for analysis (b)
Close modal

Preliminary Discussion.

There are six unknown interfacial HTCs (four silicon-indium-silicon and two silicon-copper) and seven temperature measurement histories. Because two of the measurements must be used as constraints (boundary conditions) in the analysis, there are only five temperatures to match through proper determination of the HTCs. A system of equations to determine the six HTCs will be underdetermined. This point will be explored further through analysis of the sensitivity coefficients.

Analytical Modeling

Lumped Capacitance Justification.

Any parameter estimation or inverse problem solution requires a suitable direct or forward model for the process. An order of magnitude analysis shows that a lumped capacitance model of the pucks in the stack is appropriate.

Consider one of the silicon pucks which have thermal communication through the interfacial HTCs to neighboring pucks above and below. The lateral surfaces of the pucks are assumed to be adiabatic due to the vacuum chamber environment (this neglects any radiation heat transfer to the walls of the vacuum chamber). The volume to heat transfer surface area ratio becomes L/2 = 0.0045 m. The thermal conductivity of the silicon is a strong function of temperature but for the range of temperatures of interest is of the order of 1000 W/m-K. The interfacial HTCs are to be determined in the experiment but are of the order of 10,000 W/m2-K. The corresponding Biot number is then Bi = hL/k =0.05. The criterion for valid assumption of lumped capacitance requires that the Biot number is very small relative to 1.0, and so this value suggests that lumped capacitance is a reasonable assumption.

Lumped Capacitance Model.

The nomenclature for the model is shown in Fig. 3(b). An energy balance on each of the pucks in the stack result in the following set of coupled nonlinear ordinary differential equations:
(1a)
(1b)
(1c)
(1d)
(1e)
The temperatures of the hot (Th) and cold (Tc) plates are used as boundary conditions in these equations. For the forward model, these equations can be cast into the following matrix form:
(2a)
where
(2b)
and the coefficient for the conductance hi associated with puck “k” is
(2c)
The vector b on the right side of Eq. (2a) contains the boundary condition information from the temperature measurements on the hot and cold plates. Specifically, the vector b is
(2d)

Parameter Estimation

Estimation Framework.

To estimate the interfacial conductances between the pucks, the left side of Eq. (2a) can be reformulated to consider the HTCs as unknowns as follows:
(3a)
Here, the ΔT matrix is
(3b)
Now, the history of conductances can be estimated globally over the whole time domain by minimizing the sum of the squared errors
(4a)
Here the vectors T and Y are concatenated vectors of all the temperatures at each time
(4b)
(4c)

T is the vector of temperatures computed from solution of the model equation, Eq. (3a), which depend on the unknown HTCs. Y is a vector of corresponding measured temperatures from the experiment to determine the HTCs. If there are n time steps for measurements and calculations, then vectors T and Y will be of dimension 5n ×1. The goal is to determine the values of the HTCs that best match these two vectors by minimizing the metric in Eq. (4a).

By linearizing the solution about an assumed vector of HTCs, a Gaussian iterative method using a Levenberg-like stabilizing term [1417] can be used to minimize the sum of squares function
(5a)
(5b)

Here, k is in iteration index, and Δh(k) is a correction to the current value of h(k) to improve the vector for the next iteration.

The idea behind the damping term in Eq. (5b) [14] is that the coefficient, λ(k), is proportional to the error in the estimation from the current iteration, e(k)
(5c)

and, therefore, should decrease with iterations with improved estimates for the unknowns.

The sensitivity matrix, X, in Eq. (5a) is of utmost importance. This matrix quantifies the dependence of the solution for T on the conductance parameters, h
(6a)
Since the solution for T depends on h (Eqs. (2a)(2c)), the X matrix is approximated at each iteration using a forward finite difference scheme
(6b)

Here, ε is a small fraction used to scale the current values of hi for the finite difference approximation.

If there are p parameters (HTCs), m thermocouple measurement locations, and n measurement sampling times, X is a matrix of size m × n rows and p columns. Since there are m =5 measurement locations, and p =6 parameters, the dimensions of the X matrix are 5n ×6. X is constructed from the vertical concatenation of five n ×6 sensitivity matrices, one for each measurement location. These five matrices result from computing the derivative of the five temperatures in Eq. (4b) with respect to each of the six HTCs. Each of these five resulting matrices are full; that is, not block diagonal or any other special structure. This is because a variation in any one of the h components causes a change in every measured temperature.

Discussion: Levenberg iterative technique

The Levenberg iterative technique is not a regularizing method, in the sense of inverse problems, since it does not modify the underlying objective function (the sum of squared errors). The method is best interpreted as a technique that returns a solution that is reasonably close to the starting point of iteration, but it does not address the existence of other solutions that are “equally good.” Aster et al. [18] offer this explanation:

The λI term in the LM method … looks a lot like Tikhonov regularization…. It is important to understand that this is not actually a case of Tikhonov regularization. The λI term is used to stabilize the solution of the linear system of equations which determines the search direction to be used. Because the λI term is only used as a way to improve the convergence of the algorithm and does not enter into the objective function that is being minimized, it does not regularize the nonlinear least squares problem.

Numerical Experiment.

Insight into the ability of the parameter estimation algorithm to recover the contact HTCs from an actual experiment can be gained by classic numerical experimentation. Values are assumed for the “unknown” HTCs, and, to test the ability of the algorithm to resolve these, a range of different values is assumed. Arbitrary, but reasonable, variations in the hot- and cold-plate temperatures are also assumed, and the numerical model in Eqs. (2a)(2c) is used to calculate ideal temperature measurements for the experiment. These ideal temperatures are corrupted by adding random noise of an expected magnitude to the values to create simulated data for the parameter estimation.

Figure 4 shows the variation assumed for the numerical experiment. This variation is simulating a possible response if the entire apparatus is initially at room temperature (295 K) when the cryogenic liquid suddenly begins circulating at the cold plate and lowers Tc to 80 K. No external heating is applied at the hot plate but an artificial lag is imposed on Th before the temperature decreases linearly to a hypothetical steady-state level.

Fig. 4
Assumed variation in temperature of the heater (top line) and the cold plate (bottom line) for the numerical experiment
Fig. 4
Assumed variation in temperature of the heater (top line) and the cold plate (bottom line) for the numerical experiment
Close modal

When the values of the HTCs are relatively large (of the order of 10,000 W/m2-K), the system of equations in Eqs. (2a)(2c) is very stiff. Accordingly, a lower-order solver in MATLAB, ode23s, is used to integrate the equations. Actual values for the six “unknown” HTCs used in the numerical experiments are h1 = 5000 W/m2-K, h2 = 19,000 W/m2-K, h3 = 22,000 W/m2-K, h4 = 25,000 W/m2-K, h5 = 14,000 W/m2-K, and h6 = 4000 W/m2-K.

Figure 5 shows the computed ideal temperatures responses at the five sensor locations. Gaussian random noise with a standard deviation of 0.05 K is added to the values seen in Fig. 5.

Fig. 5
Ideal temperatures resulting from the numerical experiment simulation
Fig. 5
Ideal temperatures resulting from the numerical experiment simulation
Close modal

Sensitivity Coefficients.

Some insights and expectations regarding the parameter estimation procedure can be obtained by examining the sensitivity coefficients. Figure 6 shows the scaled sensitivity coefficients for each of the six HTCs. There are five distinct sections to each of the curves corresponding to the five different sensor measurements, which are concatenated together for the estimation (see Eqs. (4b) and (4c)). Figure 6(a) shows the coefficients for the HTCs near the boundaries, which are much larger in magnitude than the others. Figures 6(b) and 6(c) plot the sensitivities for the sensors in the interior of the stack.

Fig. 6
Scaled sensitivity coefficients for the six conductances (a) conductances 1 and 6 near the boundaries, (b) interior conductances 2, and 3, and (c) interior conductances 4 and 5. ε = 0.02 is used in Eq. (6b) to approximate the derivatives. Note differences in scales on the subfigures.
Fig. 6
Scaled sensitivity coefficients for the six conductances (a) conductances 1 and 6 near the boundaries, (b) interior conductances 2, and 3, and (c) interior conductances 4 and 5. ε = 0.02 is used in Eq. (6b) to approximate the derivatives. Note differences in scales on the subfigures.
Close modal

There are two desirable characteristics of these scaled sensitivity coefficients. First, the absolute value of the scaled sensitivity coefficients should be “large” in relation to the magnitude of the sensor signal. Secondly, the sensitivity coefficients should be linearly independent [14]. Linear independence can be observed qualitatively by examining plots of the sensitivity coefficients. If any curve for a sensitivity coefficient has the same, or nearly the same, shape as another curve, then the coefficients are linearly dependent. It means that one curve can be obtained by merely scaling the other one. The mathematical consequence of linearly dependent sensitivities is that the inverse of the XTX matrix in Eq. (5a) may be singular, or nearly singular, near the minimum of S in Eq. (4a).

Figure 6 displays the sensitivity coefficients for the six HTCs. These coefficients are computed using the exact values of the assumed HTCs and so are the values at the minimum of S. By experimentation, the value of ε = 0.02 in Eq. (6b) was found adequate for the finite difference approximation. Each of the six sensitivity histories has five distinct sections corresponding to the five different measurements.

The values in Fig. 6 are the scaled sensitivity coefficients, which means that each of the coefficients computed using Eq. (6b) are multiplied by the corresponding value of the contact conductance. The resulting scaled sensitivity values have physical units of measured temperature. This is a useful way to compare the sensitivity of a particular sensor's signal to the magnitude of the response for that sensor. The magnitude of the values of the coefficients in Fig. 6 indicate that the sensitivities are of the order of at least 5–10 degrees Kelvin with those near the boundary of the order of tens of degrees Kelvin, which is sufficiently large in comparison to the 80 K to 300 K range of temperatures expected.

Visual analysis of the curves in Fig. 6 suggests that all of the sensors have highest sensitivity for the HTCs h1 and h6 in Fig. 6(a). Sensor 1 in Fig. 6(a) has the highest sensitivity to h1 of all the sensors, and sensor 5 has highest sensitivity (in magnitude) to h6. This makes sense physically since these sensors are nearest the corresponding HTC. Figures 6(b) and 6(c) show that several of the HTCs 2 through 5 have overlapping sensitivities for some sensors. For example, sensor 1 (measurement numbers less than 181) has very similar sensitivity for all four interior HTCs, and sensor 2 (measurement numbers between 181 and 362) has very similar sensitivity for HTCs h3 and h4. Similar observations can be made for sensors 3, 4, and 5 regarding HTCs h2, h3, and h4.

The numerical consequence of linear dependence of the sensitivity coefficients is seen by examination of the coefficient matrix XTX in Eq. (5b). The determinant of this matrix is 1.184 × 10−22, which is small and suggests the matrix is near singular. Another metric is the condition number of the matrix, which has the value 6.92 × 103.

The condition number is the ratio of the largest to the smallest eigenvalue of the matrix resulting from its singular value decomposition. As the condition number becomes larger, small perturbations (caused by measurement noise or model error) are amplified more in the solution (inversion of XTX). Hansen [19] gives good discussion of large condition number on solution of ill-posed problems.

In summary, there are five sensors and six unknown parameters. All sensors have strong sensitivity to the HTCs at the boundaries, but sensitivities to HTCs h2, h3, and h4 are generally confounded (linearly dependent). The lack of linear independence of the sensitivities means that unique estimation of the HTCs may not be possible [14].

Estimated Parameters

Procedure.

The noised “data” based on the precise values depicted in Fig. 5 were used in the estimation procedure described by Eq. (5) to recover the HTCs used in the numerical experiment. Convergence of the iterative procedure was gauged by monitoring both the average relative change in the sum of squares objective function (Scnvrg) and the maximum relative change in the values of the HTCs (hcnvrg). Iterations ceased when |Scnvrg| < 1% and |hcnvrg| < 0.5%.

Figures 7 and 8 display these aspects of the convergence history. As seen in these figures, the final values of the HTCs are largely realized by iteration 12 and that very small additional changes are made until the final convergence criteria are satisfied.

Fig. 7
Convergence history: the sum of squares error at each iteration
Fig. 7
Convergence history: the sum of squares error at each iteration
Close modal
Fig. 8
Convergence history: values of the conductances at each iteration
Fig. 8
Convergence history: values of the conductances at each iteration
Close modal
Results(1).

The estimated HTCs are listed in the first column of Table 1, and the true values are listed under hact with the corresponding relative errors listed in the final column. The estimates are in excellent agreement with the actual values.

Table 1

Estimated conductances, confidence intervals, actual values, and relative errors when using the original model to generate data

hest W/m2-Khact W/m2-KRel. err. (%)
h150095000−0.2
h219006190000.0
h32205822000−0.3
h42522025000−0.9
h513983140000.1
h640044000−0.1
hest W/m2-Khact W/m2-KRel. err. (%)
h150095000−0.2
h219006190000.0
h32205822000−0.3
h42522025000−0.9
h513983140000.1
h640044000−0.1

Extended Lumped Capacitance Model.

When additional physics are present in the “experiment,” which are not incorporated in the parameter estimation forward model, characteristic signatures in the residuals persist. This is demonstrated by extending the model used to generate “data.”

Two types of additional physics are incorporated: radiation heat loss from the pucks, and a simple thermocouple dynamics model. Radiation heat loss from the pucks is modeled by adding an additional term to the energy balance equations in Eqs. (1a)(1e). For puck “i
(7)

The value for emissivity (ε) is taken as 1.0 as the worst-case scenario, and the walls of the chamber (T) are assumed as 298 K.

The second phenomenon incorporated into the extended model is the effect of thermocouple sensor dynamics modeled as a simple resistance-capacitance effect [20]
(8)

Here, T is the “true” temperature computed from the solution of Eq. (7), 1/(hTCA) is the resistance to external heat flow to the thermocouple, and (ρcpV) is the capacitance of the thermocouple bead. The thermocouple is modeled as a pure chrome bead (ρcp = 3.44 × 106 J/m3-K) of diameter 0.559 mm (roughly equivalent to 30 gauge wire), and the conductance to the surrounding medium is assumed to be hTC = 1000 W/m2-K.

Equations (7) and (8) are solved and the difference (TTTC) plotted as the “thermocouple error.” The results in Fig. 9 indicate that the sensor signal lags behind the “true” temperature by as much as 18 K at the cold plate when the temperature drops suddenly. The discrepancies diminish until the temperature of the hot plate begins to decrease, then differences of about 1 K are seen.

Fig. 9
Thermocouple errors
Fig. 9
Thermocouple errors
Close modal

The TTC values are also noised with random Gaussian errors with standard deviation of 0.05 K. The resulting “data” are processed with the original algorithm to estimate the HTCs.

Results(2).

As before, the procedure converged quickly. The results of the estimation are listed in Table 2. Of course, the results are now biased due to missing physics in the estimation model. However, the results are within about 10% of the true values. Notably h6 has the highest relative error, which is due to its strong sensitivity to the data (Fig. 6) and its proximity to the T5 sensor which has the largest measurement error (Fig. 9).

Table 2

Estimated conductances, confidence intervals, actual values, and relative errors when using the extended lumped capacitance model to generate data

hest W/m2-Khact W/m2-KRel. err. (%)
h1466450006.7
h217,57219,0007.5
h319,93722,0009.4
h422,95625,0008.2
h512,72214,0009.1
h6361340009.7
hest W/m2-Khact W/m2-KRel. err. (%)
h1466450006.7
h217,57219,0007.5
h319,93722,0009.4
h422,95625,0008.2
h512,72214,0009.1
h6361340009.7

Figure 10 compares the recovered estimated temperatures from the procedure and the corresponding error residuals in relation to the measured “data” for the baseline case using the original model (Fig. 10(a)) and the extended model (Fig. 10(b)). The residuals in the figure are scaled by a factor of 50 for better visualization and the two subfigures are plotted to the same scale for easy comparison. When the model used in the estimation is appropriate for the “experiment” (in this case, the numerical simulation), residuals should be small and approximately randomly distributed, and this behavior is seen in Fig. 10(a). If the model used in the estimation is inadequate to describe the physics of the experiment, the residuals may be relatively large and have a signature of some sort that deviates from a random pattern. This behavior is seen in Fig. 10(b). In a real experiment, where the underlying model for the process is not known, this analysis of the residuals may help determine the appropriateness of the modeling in the estimation. However, this nonrandom behavior may be necessary, but not sufficient, to conclude the adequacy of the model. For ill-posed problems, it is possible that the model error will manifest as a bias with very small residual.

Fig. 10
Recovered temperature estimates and corresponding residuals relative to the data (note residuals are scaled 50x) using (a) original model and (b) extended model with unaccounted physics
Fig. 10
Recovered temperature estimates and corresponding residuals relative to the data (note residuals are scaled 50x) using (a) original model and (b) extended model with unaccounted physics
Close modal

Impact of Modeling Errors—Inverse Crime?

The excellent results found through this numerical experiment in spite of the warnings presented by the sensitivity coefficient analysis exemplify the danger of the “inverse crime” [21]:

Using the same model to generate and subsequently analyze a synthetic dataset is called an inverse crime. Inverse crimes can be avoided [by]…generat[ing] and analyz[ing] the data using two completely different types of models (p. 859).

To better assess the ability of the procedure to recover the HTCs, a different type of model is implemented to generate “data.” A finite element model is used to generate temperature histories, which are subsequently noised and used as “data” for the estimation algorithm.

Finite Element Model.

The finite element model of the experiment includes the hot plate, puck stack, and cold plate. An image of the model is seen in Fig. 11. The indium layers and temperature sensors are not meshed; instead, the indium layers are replaced with contact elements to facilitate parameterization of the thermal contact conductance at each interface. Temperature-dependent material properties are used for both copper and silicon. Rather than using temperature changes at the hot and cold plate to simulate heating, boundary conditions more consistent with the experiment are used. Specifically, a heat flux of 30,417 W/m2 is applied at the interface between the cartridge heater inside the copper hot plate assembly and the walls of the enclosing cylindrical cavity. A convection coefficient of 10,000 W/m2-K is applied to the cooling fins of the cold plate interacting with a coolant temperature of 80 K. Initial conditions consist of a uniform temperature distribution of 87 K. Assumed HTCs for the interfaces are h1 = 5000, h2 = 19,000, h3 = 22,000, h4 = 25,000, h5 = 14,000, and h6 = 4000 W/m2-K.

Fig. 11
The finite element model used to calculate responses of the experiment includes the hot plate (top), puck stack (center section), and cold plate (bottom)
Fig. 11
The finite element model used to calculate responses of the experiment includes the hot plate (top), puck stack (center section), and cold plate (bottom)
Close modal
Data for the estimation.

Temperature measurements are extracted from the finite element model by tagging the temperature on the surface at the vertical center of each puck. The values extracted from the simulation are temperatures at irregular time intervals due to the variable time stepping routine used in the solver. These values are interpolated using cubic splines to obtain temperature values on uniform time intervals, and these values are noised with random values with 0.05 K standard deviation. The resulting “data” are seen in Fig. 12.

Fig. 12
Simulated data from the finite element calculation
Fig. 12
Simulated data from the finite element calculation
Close modal

The estimation model retains the lumped capacitance formulation of Eq. (2) and is used to compute the sensitivity coefficients for the FEM “data” and the exact HTCs. These sensitivities are computed using the value of ε = 0.10 in Eq. (6b) are shown in Fig. 13. The values are lower in magnitude and noisier than previously but display similar interdependence discussed earlier in reference to Fig. 6. The lack of linear independence indicates estimation of the parameters may not be possible.

Fig. 13
Scaled sensitivity coefficients for the FEM model for the six conductances: (a) conductances 1 and 6 near the boundaries, (b) interior conductances 2, and 3, and (c) interior conductances 4 and 5. ε = 0.20 is used in Eq. (6b) to approximate the derivatives.
Fig. 13
Scaled sensitivity coefficients for the FEM model for the six conductances: (a) conductances 1 and 6 near the boundaries, (b) interior conductances 2, and 3, and (c) interior conductances 4 and 5. ε = 0.20 is used in Eq. (6b) to approximate the derivatives.
Close modal
Results(3).

The algorithm in Eq. (5a) is used to process the data from Fig. 12. A 60 s time-step was used to extract simulated measurements from the FEM output and a uniform initial guess of 10,000 W/m2-K is assumed for the unknown HTCs. The iterative procedure converged smoothly in 17 steps.

The values found for the six HTCs are listed in Table 3 along with the actual values used in the FEM simulation and the relative errors. The values found are in poor agreement with the actual values. Additional estimations with different initial guesses produce the results displayed in Table 4. Note that a fortunate initial guess, such as hguess = 15,000 W/m2-K, produces reasonable results, but general arbitrary guesses produce poor results.

Table 3

Estimated conductances, confidence intervals, actual values, and relative errors when using the FEM model to generate data

hest W/m2-Khact W/m2-KRel. err. (%)
h13,3525,000−33.0
h213,72219,000−27.8
h313,97222,000−36.5
h415,73525,000−37.1
h510,45114,000−25.4
h62,4344,000−39.1
hest W/m2-Khact W/m2-KRel. err. (%)
h13,3525,000−33.0
h213,72219,000−27.8
h313,97222,000−36.5
h415,73525,000−37.1
h510,45114,000−25.4
h62,4344,000−39.1
Table 4

Estimated conductances and relative errors when using Levenberg iterations and the FEM model to generate data for differing initial guesses for h

hguess = 5000hguess = 10,000hguess = 15,000hguess = 20,000
hestRel. err.hestRel. err.hestRel. err.hestRel. err.
h11,753−64.9%3,352−33.0%4,791−4.2%6,48529.7%
h27,174−62.2%13,722−27.8%19,5993.2%26,57739.9%
h37,287−66.9%13,972−36.5%20,018−9.0%27,02622.8%
h48,201−67.2%15,735−37.1%22,429−10.3%30,46721.9%
h55,452−61.1%10,451−25.4%14,9616.9%20,24144.6%
h61,269−68.3%2,434−39.1%3,481−13.0%4,71617.9%
Smin9.12.62.82.5
heff515.3 W/m2-K987.3 W/m2-K1,411.6 W/m2-K1,911.7 W/m2-K
qss0.0384 MW/m20.0736 MW/m20.1053 MW/m20.1425 MW/m2
hguess = 5000hguess = 10,000hguess = 15,000hguess = 20,000
hestRel. err.hestRel. err.hestRel. err.hestRel. err.
h11,753−64.9%3,352−33.0%4,791−4.2%6,48529.7%
h27,174−62.2%13,722−27.8%19,5993.2%26,57739.9%
h37,287−66.9%13,972−36.5%20,018−9.0%27,02622.8%
h48,201−67.2%15,735−37.1%22,429−10.3%30,46721.9%
h55,452−61.1%10,451−25.4%14,9616.9%20,24144.6%
h61,269−68.3%2,434−39.1%3,481−13.0%4,71617.9%
Smin9.12.62.82.5
heff515.3 W/m2-K987.3 W/m2-K1,411.6 W/m2-K1,911.7 W/m2-K
qss0.0384 MW/m20.0736 MW/m20.1053 MW/m20.1425 MW/m2
A regularized estimation procedure.
Better convergence behavior can be obtained by modifying the objective function of Eq. (4a) as
(9a)
Here, the last term is a Tikhonov-like regularizing term intended to constrain the behavior of h in some way ([11, Chapter 4]. Let
(9b)
and substitute for both T(k+ 1) and h(k+ 1) (from Eq. (5a)) into Eq. (9a). The result is
(9c)
After minimizing Eq. (9c) with respect to Δh(k), the regularized iteration correction equation is
(9d)

The superscript “(k)” has been dropped in Eq. (9d) for clarity in the matrix notations.

Equation (9d) is very similar to Eq. (5b) but differs significantly because of the last term in the parenthesis on the right side. This term helps stabilize the solution by decreasing the error differential (Y − T) as h increases. The regularizing parameter αTik must be chosen by the analyst to achieve reasonable performance.

Results(4).

The algorithm of Eq. (9d) is used to process the “data” from the FEM in the same manner as the Levenberg algorithm of Eq. (5b). Zeroth-order regularization is used (meaning W = I, the identity matrix) in Eq. (9d), and a value of the Tikhonov regularizing parameter of αTik = 5 × 10−10 was found heuristically to give reasonable results. It was not possible to realize convergence to the same level as unregularized Levenberg method, but convergence was declared when S changes less than 1% and the maximum change in h is less than 2%.

For starting values of h less than 10,000 W/m2-K, convergence was readily achieved. However, some larger starting values of h resulted in nonconverging solutions. This is attributable to the regularizing term (last term on the right side) of Eq. (9d): when h becomes very large, it dominates the correction procedure. This behavior was modified by performing a line search at each iteration to find the optimal fraction of the new Δh(k) to apply to minimize S for that iteration. That is, instead if Eq. (5a), a modified update equation is used
(10)

The value of Δh is found using Eq. (9d), and the value of f that minimizes the sum squared error of Eq. (4a) is used at each iteration. This line search made the iterative procedure more robust and convergent from different starting points.

Table 5 shows the results obtained for the HTCs using Tikhonov iterations. The accuracy of the results depends on the initial guess value, but the results from different initial guesses are much more consistent than those obtained with Levenberg iterations.

Table 5

Estimated conductances and relative errors when using Tikhonov iterations and the FEM model to generate data for differing initial guesses for h

hguess = 5000hguess = 10,000hguess = 15,000hguess = 20,000
hestRel. err.hestRel. err.hestRel. err.hestRel. err.
h14,058−18.8%4,087−18.3%4,659−6.8%4,894−2.1%
h216,598−12.9%16,756−11.8%19,0860.5%20,0385.5%
h316,893−23.2%17,039−22.6%19,405−11.8%20,392−7.3%
h419,030−23.9%19,180−23.3%21,880−12.5%22,977−8.1%
h512,645−9.7%12,753−8.9%14,5333.8%15,2689.1%
h62,948−26.3%2,971−25.7%3,389−15.3%3,560−11.0%
Smin3.03.02.82.8
heff1195.1 W/m2-K1204.6 W/m2-K1373.4 W/m2-K1442.5 W/m2-K
qss0.0891 MW/m20.0898 MW/m20.1024 MW/m20.1076 MW/m2
hguess = 5000hguess = 10,000hguess = 15,000hguess = 20,000
hestRel. err.hestRel. err.hestRel. err.hestRel. err.
h14,058−18.8%4,087−18.3%4,659−6.8%4,894−2.1%
h216,598−12.9%16,756−11.8%19,0860.5%20,0385.5%
h316,893−23.2%17,039−22.6%19,405−11.8%20,392−7.3%
h419,030−23.9%19,180−23.3%21,880−12.5%22,977−8.1%
h512,645−9.7%12,753−8.9%14,5333.8%15,2689.1%
h62,948−26.3%2,971−25.7%3,389−15.3%3,560−11.0%
Smin3.03.02.82.8
heff1195.1 W/m2-K1204.6 W/m2-K1373.4 W/m2-K1442.5 W/m2-K
qss0.0891 MW/m20.0898 MW/m20.1024 MW/m20.1076 MW/m2

Discussion

The inability of the Levenberg estimation algorithm to converge to a unique minimum is not surprising in consideration of the sensitivity analysis. There simply is not enough information conveyed in the five sensor signals to uniquely identify the six free parameters in the process.

The nonconvergence of the Levenberg iteration is a consequence of the shallowness of the objective function. This means that there is one “global minimum” surrounded by a very long valley in m-dimension space, where m is the number of unknowns. This also reinforces the point that the Levenberg–Marquardt method is not a regularizing technique. Regularizing techniques explicitly change the curvature of the objective function.

Results at the bottom of Table 4 further inform the nonconvergence and data insufficiency. For every value of h1 in Table 4, there are complementary values of the remaining HTCs that result in a similar Smin. This confirms that the S(h) surface has a long flat valley, so that the final value of Smin for each search is similar in magnitude but corresponds to a different collection of HTCs.

The value of heff in Table 4 is the effective value for the converged HTCs for each search case which results from the series combination of these HTCs. As seen in the table, the effective values increase with increasing magnitude of initial guess. The last line of the table contains the result of a calculation of the steady-state heat flux through the puck stack based on the average final values of Th and Tc (see Fig. 12): qss = heff(Th − Tc). This means that the same temperature history at the measurement locations will result for different heat fluxes applied to the stack but for different HTCs.

Table 6 shows the results from Table 4 when they are normalized relative to the values of h1 found from each initial guess solution. The values are highly insensitive to the initial guess value, indicating that the Levenberg algorithm can determine the relationship between the HTCs but not the precise values. This is a result of the linear dependence of the sensitivity coefficients: unique estimation of all six parameters is not possible. Table 6 suggests that, if one of the HTCs can be found independently, then the estimation procedure can be used to find the remaining values.

Table 6

Results from Table 4 normalized as h/h1

hguess = 5000hguess = 10,000hguess = 15,000hguess = 20,000
h1/h11111
h2/h14.084.124.094.11
h3/h14.164.204.184.20
h4/h14.664.664.684.70
h5/h13.113.123.123.12
h6/h10.720.720.730.72
hguess = 5000hguess = 10,000hguess = 15,000hguess = 20,000
h1/h11111
h2/h14.084.124.094.11
h3/h14.164.204.184.20
h4/h14.664.664.684.70
h5/h13.113.123.123.12
h6/h10.720.720.730.72

Table 5 illustrates how regularization can help extract results from an ill-conditioned system of equations. The results in Table 5 are much more consistent for differing starting values. The ratios of these HTCs relative to h1 are shown in Table 7 and are very consistent for the different initial guesses.

Table 7

Results from Table 5 normalized as h/h1

hguess = 5000hguess = 10,000hguess = 15,000hguess = 20,000
h1/h11111
h2/h14.094.104.104.09
h3/h14.164.174.174.17
h4/h14.694.694.704.69
h5/h13.123.123.123.12
h6/h10.730.730.730.73
hguess = 5000hguess = 10,000hguess = 15,000hguess = 20,000
h1/h11111
h2/h14.094.104.104.09
h3/h14.164.174.174.17
h4/h14.694.694.704.69
h5/h13.123.123.123.12
h6/h10.730.730.730.73

At least two approaches can be used to close the problem to find the unique set of parameters from the experiment. The first approach is to use a special interface at the location of h1 which has a known resistance. If this value is known, then the five thermocouple readings will be sufficient to uniquely estimate the remaining five parameters. The second approach is to always run the experiment to steady-state (or approximately steady-state) so that the heat flux through puck stack will be known. With a known value for the steady heat flow through the stack, the unique collection of HTCs that cause this heat flux can be determined (see Table 4).

Conclusion

This investigation underscores the importance of examination of sensitivity coefficients in a parameter estimation procedure and the pitfalls of the “Inverse Crime.” An existing experiment is evaluated for application to determine interfacial heat transfer coefficients (conductances) in silicon-indium-silicon contact. A simple lumped capacitance model and Levenberg-like iterative procedure is used to estimate the conductances. Examination of the sensitivity coefficients provides insight that the estimation of all the unknown parameters may be difficult or impossible due to the linear dependence of these coefficients. Initial numerical experiments using the same lumped capacitance model to generate “data” suggest that the estimation may be possible. Subsequent numerical experiments using “data” computed from a FEM clearly show that only poor estimates of the values are possible. Modification of the iterative algorithm to include some regularization, such as Tikhonov regularization, improves the ability of the iterations to provide a consistent set of results for different initial starting vectors. Because the parameters are not linearly independent, they cannot all be determined simultaneously. However, if one of the conductances can be found independently, the numerical experiment suggests the other values can be recovered from the estimation procedure.

Acknowledgment

The authors are grateful to the reviewers and the associate editor for their time evaluating the paper and for many helpful comments. The paper is substantially strengthened through the review process.

Funding Data

  • Lawrence Berkeley National Laboratory (Award ID: (internal); Funder ID: 10.13039/100006235).

Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

Nomenclature

A =

coefficient matrix

b =

right hand side vector

cp =

specific heat, J/kg-K

f =

fraction for update, Eq. (10)

h =

conductance (heat transfer coefficient) W/m2-K

h =

vector of all the heat transfer coefficients

H =

coefficient (Eq. (2c))

HTC =

heat transfer coefficient

i =

integer index

k =

iteration index

m =

number of measurements

m =

mass, kg

n =

number of samples

p =

number of parameters

S =

sum of squared error (Eq. (4a))

t =

time, s

t =

Student's t-statistic

T =

vector of discrete temperatures

Tl =

temperature, K

X =

sensitivity matrix (Eq. (6a))

Y =

vector of temperature measurements

αTik =

Tikhonov regularization parameter

ε =

Fraction used in finite difference approximation (Eq. (6b))

λ =

Levenberg iteration parameter (Eq. (5b))

ρ =

density, kg/m3

References

1.
BERKELEY LAB
,
2023
, “
About the ALS
,” accessed Aug. 4, 2023, https://als.lbl.gov/about/about-the-als/
2.
BERKELEY LAB
,
2023
, “
ALS Work Using Protein Crystallography
,” accessed Aug. 9, 2023, https://als.lbl.gov/tag/protein-crystallography/
3.
Cutler
,
G.
,
Cocco
,
D.
,
DiMasi
,
E.
,
Morton
,
S.
,
Sanchez del Rio
,
M.
, and
Padmore
,
H.
,
2020
, “
A Cantilevered Liquid-Nitrogen-Cooled Silicon Mirror for the Advanced Light Source Upgrade
,”
J. Synchrotron Radiat.
,
27
(
5
), pp.
1131
1140
.10.1107/S1600577520008930
4.
Cutler
,
G.
,
Bryant
,
D.
,
Cocco
,
D.
,
DiMasi
,
E.
,
Mccombs
,
K.
,
Morton
,
S.
,
del Rio
,
M. S.
,
Smith
,
N.
, and
Padmore
,
H.
,
2020
, “
An Update on Development of a Cryogenically Cooled-Silicon Mirror for the Advanced Light Source Upgrade Project
,”
Advances in X-Ray/EUV Optics and Components XV
, Aug. 24–Sept. 4, pp.
50
56
.10.1117/12.2568101
5.
Cutler
,
G.
,
Cocco
,
D.
,
Bentley
,
B.
,
Cervantes
,
M.
,
Chavez
,
P.
,
Chrzan
,
J.
,
DiMaggio
,
S.
, et al.,
2023
, “
Experimental Testing of a Prototype Cantilevered Liquid-Nitrogen-Cooled Silicon Mirror
,”
J. Synchrotron Radiat.
,
30
(
1
), pp.
76
83
.10.1107/S1600577522010700
6.
Khounsary
,
A. M.
,
Chojnowski
,
D.
,
Assoufid
,
L.
, and
Worek
,
W. M.
,
1997
, “
Thermal Contact Resistance Across a Copper-Silicon Interface
,” High Heat Flux Synchrotron Radiation Beamlines, San Diego, CA, July 27–Aug. 1, pp.
45
51
.10.1117/12.294497
7.
Takeuchi
,
T.
,
Tanaka
,
M.
,
Senba
,
Y.
,
Ohashi
,
H.
, and
Goto
,
S.
,
2011
, “
Thermal-Contact-Conductance Measurement for High-Heat-Load Optics Components at SPring-8
,”
Advances in X-Ray/EUV Optics and Components VI
, San Diego, CA, Aug. 21–25, pp.
303
308
.10.1117/12.894612
8.
Asano
,
M.
,
Ogata
,
J.
, and
Yosinaga
,
Y.
,
1993
, “
Quantitative Evaluation of Contact Thermal Conductance in a Vacuum as a Result of Simulating the Effect of Cooling
,”
High Heat Flux Engineering
, San Diego, CA, July 19–24, pp.
652
656
.10.1117/12.140520
9.
Marion
,
P.
,
Zhang
,
L.
,
Vallet
,
L.
, and
Lesourd
,
M.
,
2004
, “
Thermal Contact Resistance Study
,”
Proceedings of MEDSI
, Grenoble, France, May 24–27.
10.
Fieberg
,
C.
, and
Kneer
,
R.
,
2008
, “
Determination of Thermal Contact Resistance From Transient Temperature Measurements
,”
Int. J. Heat Mass Transfer
,
51
(
5–6
), pp.
1017
1023
.10.1016/j.ijheatmasstransfer.2007.05.004
11.
Woodbury
,
K. A.
,
Najafi
,
H.
,
de Monte
,
F.
, and
Beck
,
J. V.
,
2023
,
Inverse Heat Conduction: Ill-Posed Problems
, 2nd ed.,
Wiley, Inc.
,
Hoboken, NJ
.
12.
Beck
,
J. V.
,
1969
, “
Determination of Optimum, Transient Experiments for Thermal Contact Conductance
,”
Int. J. Heat Mass Transfer
,
12
(
5
), pp.
621
633
.10.1016/0017-9310(69)90043-X
13.
Bi
,
D.
,
Chen
,
H.
, and
Ye
,
T.
,
2012
, “
Influences of Temperature and Contact Pressure on Thermal Contact Resistance at Interfaces at Cryogenic Temperatures
,”
Cryogenics
,
52
(
7–9
), pp.
403
409
.10.1016/j.cryogenics.2012.03.006
14.
Beck
,
J. V.
, and
Arnold
,
K. J.
,
1977
,
Parameter Estimation in Engineering and Science
,
Wiley
,
New York
.
15.
Levenberg
,
K.
,
1944
, “
A Method for the Solution of Certain Non-Linear Problems in Least Squares
,”
Q. Appl. Math.
,
2
(
2
), pp.
164
168
.10.1090/qam/10666
16.
Marquardt
,
D. W.
,
1963
, “
An Algorithm for Least-Squares Estimation of Nonlinear Parameters
,”
J. Soc. Ind. Appl. Math.
,
11
(
2
), pp.
431
441
.10.1137/0111030
17.
Marquaridt
,
D. W.
,
1970
, “
Generalized Inverses, Ridge Regression, Biased Linear Estimation, and Nonlinear Estimation
,”
Technometrics
,
12
(
3
), pp.
591
612
.10.2307/1267205
18.
Aster
,
R. C.
,
Orchers
,
B.
, and
Thurber
,
C. H.
,
2019
,
Parameter Estimation and Inverse Problems
, 3rd ed.,
Elsevier
,
Amsterdam, The Netherlands
.
19.
Hansen
,
P. C.
,
1998
,
Rank-Deficient and Discrete Ill-Posed Problems
,
Siam
,
Philadelphia, PA
.
20.
Woodbury
,
K. A.
,
1990
, “
Effect of Thermocouple Sensor Dynamics on Surface Heat Flux Predictions Obtained Via Inverse Heat Transfer Analysis
,”
Int. J. Heat Mass Transfer
,
33
(
12
), pp.
2641
2649
.10.1016/0017-9310(90)90200-E
21.
Howell
,
J. R.
,
Menguc
,
M. P.
,
Daun
,
K.
, and
Siegel
,
R.
,
2021
,
Thermal Radiation Heat Transfer, Seventh
,
CRC Press
,
Boca Raton, FL, p. 859
.