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 [3–5]. 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).
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.
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.
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.
Parameter Estimation
Estimation Framework.
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).
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.
and, therefore, should decrease with iterations with improved estimates for the unknowns.
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.
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.
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.
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%.
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.
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.”
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.
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 (T–TTC) 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.
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).
hest W/m2-K | hact W/m2-K | Rel. err. (%) | |
---|---|---|---|
h1 | 4664 | 5000 | 6.7 |
h2 | 17,572 | 19,000 | 7.5 |
h3 | 19,937 | 22,000 | 9.4 |
h4 | 22,956 | 25,000 | 8.2 |
h5 | 12,722 | 14,000 | 9.1 |
h6 | 3613 | 4000 | 9.7 |
hest W/m2-K | hact W/m2-K | Rel. err. (%) | |
---|---|---|---|
h1 | 4664 | 5000 | 6.7 |
h2 | 17,572 | 19,000 | 7.5 |
h3 | 19,937 | 22,000 | 9.4 |
h4 | 22,956 | 25,000 | 8.2 |
h5 | 12,722 | 14,000 | 9.1 |
h6 | 3613 | 4000 | 9.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.
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.
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.
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.
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.
hest W/m2-K | hact W/m2-K | Rel. err. (%) | |
---|---|---|---|
h1 | 3,352 | 5,000 | −33.0 |
h2 | 13,722 | 19,000 | −27.8 |
h3 | 13,972 | 22,000 | −36.5 |
h4 | 15,735 | 25,000 | −37.1 |
h5 | 10,451 | 14,000 | −25.4 |
h6 | 2,434 | 4,000 | −39.1 |
hest W/m2-K | hact W/m2-K | Rel. err. (%) | |
---|---|---|---|
h1 | 3,352 | 5,000 | −33.0 |
h2 | 13,722 | 19,000 | −27.8 |
h3 | 13,972 | 22,000 | −36.5 |
h4 | 15,735 | 25,000 | −37.1 |
h5 | 10,451 | 14,000 | −25.4 |
h6 | 2,434 | 4,000 | −39.1 |
hguess = 5000 | hguess = 10,000 | hguess = 15,000 | hguess = 20,000 | |||||
---|---|---|---|---|---|---|---|---|
hest | Rel. err. | hest | Rel. err. | hest | Rel. err. | hest | Rel. err. | |
h1 | 1,753 | −64.9% | 3,352 | −33.0% | 4,791 | −4.2% | 6,485 | 29.7% |
h2 | 7,174 | −62.2% | 13,722 | −27.8% | 19,599 | 3.2% | 26,577 | 39.9% |
h3 | 7,287 | −66.9% | 13,972 | −36.5% | 20,018 | −9.0% | 27,026 | 22.8% |
h4 | 8,201 | −67.2% | 15,735 | −37.1% | 22,429 | −10.3% | 30,467 | 21.9% |
h5 | 5,452 | −61.1% | 10,451 | −25.4% | 14,961 | 6.9% | 20,241 | 44.6% |
h6 | 1,269 | −68.3% | 2,434 | −39.1% | 3,481 | −13.0% | 4,716 | 17.9% |
Smin | 9.1 | 2.6 | 2.8 | 2.5 | ||||
heff | 515.3 W/m2-K | 987.3 W/m2-K | 1,411.6 W/m2-K | 1,911.7 W/m2-K | ||||
qss | 0.0384 MW/m2 | 0.0736 MW/m2 | 0.1053 MW/m2 | 0.1425 MW/m2 |
hguess = 5000 | hguess = 10,000 | hguess = 15,000 | hguess = 20,000 | |||||
---|---|---|---|---|---|---|---|---|
hest | Rel. err. | hest | Rel. err. | hest | Rel. err. | hest | Rel. err. | |
h1 | 1,753 | −64.9% | 3,352 | −33.0% | 4,791 | −4.2% | 6,485 | 29.7% |
h2 | 7,174 | −62.2% | 13,722 | −27.8% | 19,599 | 3.2% | 26,577 | 39.9% |
h3 | 7,287 | −66.9% | 13,972 | −36.5% | 20,018 | −9.0% | 27,026 | 22.8% |
h4 | 8,201 | −67.2% | 15,735 | −37.1% | 22,429 | −10.3% | 30,467 | 21.9% |
h5 | 5,452 | −61.1% | 10,451 | −25.4% | 14,961 | 6.9% | 20,241 | 44.6% |
h6 | 1,269 | −68.3% | 2,434 | −39.1% | 3,481 | −13.0% | 4,716 | 17.9% |
Smin | 9.1 | 2.6 | 2.8 | 2.5 | ||||
heff | 515.3 W/m2-K | 987.3 W/m2-K | 1,411.6 W/m2-K | 1,911.7 W/m2-K | ||||
qss | 0.0384 MW/m2 | 0.0736 MW/m2 | 0.1053 MW/m2 | 0.1425 MW/m2 |
A regularized estimation procedure.
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%.
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.
hguess = 5000 | hguess = 10,000 | hguess = 15,000 | hguess = 20,000 | |||||
---|---|---|---|---|---|---|---|---|
hest | Rel. err. | hest | Rel. err. | hest | Rel. err. | hest | Rel. err. | |
h1 | 4,058 | −18.8% | 4,087 | −18.3% | 4,659 | −6.8% | 4,894 | −2.1% |
h2 | 16,598 | −12.9% | 16,756 | −11.8% | 19,086 | 0.5% | 20,038 | 5.5% |
h3 | 16,893 | −23.2% | 17,039 | −22.6% | 19,405 | −11.8% | 20,392 | −7.3% |
h4 | 19,030 | −23.9% | 19,180 | −23.3% | 21,880 | −12.5% | 22,977 | −8.1% |
h5 | 12,645 | −9.7% | 12,753 | −8.9% | 14,533 | 3.8% | 15,268 | 9.1% |
h6 | 2,948 | −26.3% | 2,971 | −25.7% | 3,389 | −15.3% | 3,560 | −11.0% |
Smin | 3.0 | 3.0 | 2.8 | 2.8 | ||||
heff | 1195.1 W/m2-K | 1204.6 W/m2-K | 1373.4 W/m2-K | 1442.5 W/m2-K | ||||
qss | 0.0891 MW/m2 | 0.0898 MW/m2 | 0.1024 MW/m2 | 0.1076 MW/m2 |
hguess = 5000 | hguess = 10,000 | hguess = 15,000 | hguess = 20,000 | |||||
---|---|---|---|---|---|---|---|---|
hest | Rel. err. | hest | Rel. err. | hest | Rel. err. | hest | Rel. err. | |
h1 | 4,058 | −18.8% | 4,087 | −18.3% | 4,659 | −6.8% | 4,894 | −2.1% |
h2 | 16,598 | −12.9% | 16,756 | −11.8% | 19,086 | 0.5% | 20,038 | 5.5% |
h3 | 16,893 | −23.2% | 17,039 | −22.6% | 19,405 | −11.8% | 20,392 | −7.3% |
h4 | 19,030 | −23.9% | 19,180 | −23.3% | 21,880 | −12.5% | 22,977 | −8.1% |
h5 | 12,645 | −9.7% | 12,753 | −8.9% | 14,533 | 3.8% | 15,268 | 9.1% |
h6 | 2,948 | −26.3% | 2,971 | −25.7% | 3,389 | −15.3% | 3,560 | −11.0% |
Smin | 3.0 | 3.0 | 2.8 | 2.8 | ||||
heff | 1195.1 W/m2-K | 1204.6 W/m2-K | 1373.4 W/m2-K | 1442.5 W/m2-K | ||||
qss | 0.0891 MW/m2 | 0.0898 MW/m2 | 0.1024 MW/m2 | 0.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.
hguess = 5000 | hguess = 10,000 | hguess = 15,000 | hguess = 20,000 | |
---|---|---|---|---|
h1/h1 | 1 | 1 | 1 | 1 |
h2/h1 | 4.08 | 4.12 | 4.09 | 4.11 |
h3/h1 | 4.16 | 4.20 | 4.18 | 4.20 |
h4/h1 | 4.66 | 4.66 | 4.68 | 4.70 |
h5/h1 | 3.11 | 3.12 | 3.12 | 3.12 |
h6/h1 | 0.72 | 0.72 | 0.73 | 0.72 |
hguess = 5000 | hguess = 10,000 | hguess = 15,000 | hguess = 20,000 | |
---|---|---|---|---|
h1/h1 | 1 | 1 | 1 | 1 |
h2/h1 | 4.08 | 4.12 | 4.09 | 4.11 |
h3/h1 | 4.16 | 4.20 | 4.18 | 4.20 |
h4/h1 | 4.66 | 4.66 | 4.68 | 4.70 |
h5/h1 | 3.11 | 3.12 | 3.12 | 3.12 |
h6/h1 | 0.72 | 0.72 | 0.73 | 0.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.
hguess = 5000 | hguess = 10,000 | hguess = 15,000 | hguess = 20,000 | |
---|---|---|---|---|
h1/h1 | 1 | 1 | 1 | 1 |
h2/h1 | 4.09 | 4.10 | 4.10 | 4.09 |
h3/h1 | 4.16 | 4.17 | 4.17 | 4.17 |
h4/h1 | 4.69 | 4.69 | 4.70 | 4.69 |
h5/h1 | 3.12 | 3.12 | 3.12 | 3.12 |
h6/h1 | 0.73 | 0.73 | 0.73 | 0.73 |
hguess = 5000 | hguess = 10,000 | hguess = 15,000 | hguess = 20,000 | |
---|---|---|---|---|
h1/h1 | 1 | 1 | 1 | 1 |
h2/h1 | 4.09 | 4.10 | 4.10 | 4.09 |
h3/h1 | 4.16 | 4.17 | 4.17 | 4.17 |
h4/h1 | 4.69 | 4.69 | 4.70 | 4.69 |
h5/h1 | 3.12 | 3.12 | 3.12 | 3.12 |
h6/h1 | 0.73 | 0.73 | 0.73 | 0.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