Algorithms for the estimation of transient surface heat flux during ultra-fast surface cooling

International Journal of Heat and Mass Transfer 100 (2016) 1–10

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt

Algorithms for the estimation of transient surface heat flux during ultra-fast surface cooling Zhi-Fu Zhou, Teng-Yu Xu, Bin Chen ⇑ State Key Laboratory of Multiphase Flow in Power Engineering, Xi’an Jiaotong University, Xi’an 710049, China

a r t i c l e

i n f o

Article history: Received 9 January 2016 Received in revised form 17 April 2016 Accepted 19 April 2016 Available online 4 May 2016 Keywords: Surface heat flux Inverse heat conduction problem Duhamel’s theorem Sequential function specification method Transfer function method Cryogen spray cooling

a b s t r a c t Surface heat flux is an important parameter in various industrial applications. This paper compared several algorithms for estimating surface heat flux using both simulated heat flux and measured temperature in response to ultra-fast surface cooling. Two different strategies were employed to measure the surface temperature, using the thin film thermocouple with a 2 lm depth deposited directly on the cooling substrate surface (direct surface temperature measurement) and the fine thermocouple with a 100 lm bead diameter placed on the surface of the cooling substrate and covered with aluminum foil (indirect surface temperature measurement). The algorithms of Duhamel’s theorem, sequential function specification method (SFSM) and transfer function method are briefly analyzed. A new method is proposed that the surface temperature is first calculated based on Duhamel’s theorem and then the surface heat flux can be estimated, in order to improve the accuracy of results in the case of the indirect surface temperature measurement method. The transfer function is obtained by solving an auxiliary problem and thus the transfer function method can be implemented to solve multilayer, semi-infinite inverse heat conduction problems. Accuracy and sensitivity to noise are examined using both the simulated triangular pulse heat flux and the measured temperature data. Duhamel’s theorem is insensitive to noise, but is unsuitable for predicting surface heat flux in the indirect measurement method of surface temperature. The SFS method can provide acceptable results of surface heat flux using both measurement methods. However, a noticeable discrepancy exists as the heat flux changes extremely quickly and it suffers considerably from noise. The transfer function and newly proposed methods are both effective in inhibiting noise, and produce very similar results, which match the simulated heat flux exactly and have a negligible standard deviation and residual temperature. Ó 2016 Elsevier Ltd. All rights reserved.

1. Introduction Spray cooling is an efficient and powerful thermal management technique for achieving ultra-fast surface cooling and high heat flux removal. It is widely applied in the steel industry, electronics devices, power plant and laser dermatology for vascular skin lesions [1–7]. Generally, in such fields of thermal management engineering, surface heat flux is a key parameter for assessing whether the equipment or body works under the normal heat load condition. However, it is relatively difficult to measure the timevarying heat flux directly at the solid surface for the initial cooling period or the pulse spray cooling. Alternatively, it can usually be estimated from the temperature measurement made at accessible locations, which is termed an inverse heat conduction problem (IHCP). ⇑ Corresponding author. E-mail address: [email protected] (B. Chen). http://dx.doi.org/10.1016/j.ijheatmasstransfer.2016.04.058 0017-9310/Ó 2016 Elsevier Ltd. All rights reserved.

The IHCP is mathematically ill-posed, and far more difficult than the direct problem. It is extremely sensitive to measurement noise and suffers considerably from the lag and damping of the measurements. Several analytical and numerical methods have been proposed for the solution of IHCPs, such as specified sequential function method (SFSM), regularization method, and transfer function method. The SFS method is one of the most widely used algorithms for solving IHCPs, as proposed by Beck in 1985 [8]. This method minimizes the effect of random errors by using future temperature data based on the least square method. The regulation method estimates all of the heat flux components simultaneously for all time and are usually presented as whole domain methods. It can be applied generally, whereas it is relatively complicated in mathematical terms [9,10]. The transfer function method analogizes the heat conduction problems to dynamic systems, where heat flux is treated as the input of the system and the temperature profile as the response [11]. By using the Laplace transform, the relation between input and output can be given by the transfer

2

Z.-F. Zhou et al. / International Journal of Heat and Mass Transfer 100 (2016) 1–10

Nomenclature c f G h H k K n q Q r S SM t T x x(t) X y(t) Y

specific heat capacity (J/(kg K)) surface temperature (K) Green’s function transfer function Laplace transform of transfer function thermal conductivity (W/(m K)) whole time step of measured data whole time step of simulated data heat flux (W/m2) Laplace transform of heat flux future time steps standard deviation of temperature (K) error between measured and estimated temperature (K2) time (s) temperature (K) spatial coordinate (m) input function Laplace transform of input function output function Laplace transform of output function, or the measure temperature (K)

function as determined by the Green’s function [12]. This method is simple in concept and one of the most accurate ways of estimating surface heat flux. However, it is relatively difficult to determine the analytic solution of the transfer function for the complex geometry problem. Although IHCPs have been extensively investigated with regard to various other applications, little work has been conducted related to heat flux on skin surface during cryogen spray cooling in laser dermatology. Cryogen spray cooling with several tens of milliseconds provides an effective way to cool the epidermis and therefore reduce the thermal injuries induced by the absorption of laser energy by melanin, when lasers are radiated on the skin of port wine stain (PWS) patients [13,14]. The heat flux on skin surface during CSC is strongly time-dynamic, which is a critical parameter for judging the cooling effect and characterizing the interaction between skin and spray cooling. However, it is relatively difficult to measure surface heat flux directly. Therefore, the indirect technique employing internal temperature measurement is often used to estimate surface heat flux. Two algorithms have been employed to estimate heat flux based on measured temperature data due to CSC. Tunnell et al. predicted surface heat flux during and following cryogen spurt using the SFS method based on internal temperature measurement, in which a wire-like thermocouple with a bead diameter of 30 lm was imbedded in epoxy resin [15–18]. Aguilar et al. also employed the SFS method to estimate surface heat flux based on internal temperature data using a faster response temperature measurement method [19–21]. In their work, a fine thermocouple (50 lm bead diameter) was placed on the surface of cooling substrate (Plexiglas) and covered by a thin (20 lm) layer of aluminum. Using the same temperature measurement method, Franco et al. calculated the surface heat flux by solving a direct problem through Duhamel’s theorem, where the measured temperature was treated as the real surface temperature [22,23]. More recently, Zhou et al. used a thin film thermocouple with a depth of 2 lm deposited on an epoxy resin surface through the Magnetron spurting technique in order to measure surface temperature directly during CSC, and then predicted the surface heat flux based on Duhamel’s theorem [6,13,24]. From the above literature review, various temperature measurement methods

Z

sensitivity coefficient (K m2/W)

Greek symbols a thermal diffusivity (m2/s) f standard deviation of heat flux (W/m2) h temperature difference (K) H Laplace transform of temperature difference q density (kg/m3) s current time step u non-dimensional stream function Subscripts c position of temperature measurement erf error function i time index or layer index k chosen time index L Laplace transform r residual 0 initial state 1 first layer

and algorithms for estimating surface heat flux have been employed. However, few of them provide a comparison of the algorithms for the estimation of surface heat flux under different surface temperature measurement methods. The accuracy and applicability of the algorithms should be further explored to promote the accurate estimation of surface heat flux due to CSC. In this paper, we use epoxy resin as the cooling substrate and measure the surface temperature during and following a pulse spray cooling using R404A with two different measurement methods: the thin film thermocouple with a 2 lm depth deposited directly on the substrate surface and the fine thermocouple with a 100 lm bead diameter placed on the surface of the substrate and covered with aluminum foil. Duhamel’s theorem and the SFS method are implemented to predict time-varying surface heat flux. Then, a new method is proposed based on Duhamel’s theorem to improve the accuracy of surface heat flux when using the indirect surface temperature measurement method. An auxiliary problem is established to obtain the transfer function numerically, and thus this method can be implemented to solve multilayer, semi-infinite inverse heat conduction problems. The accuracy and applicability of each algorithm in different measuring methods are evaluated based on both the simulated heat flux and the experimental temperature data. 2. Experiment system 2.1. Spray cooling system Fig. 1 shows a schematic of the cryogen spray cooling system. It consists of a cryogen container for the storage of liquid cryogen, a 3D positioner for controlling the distance between the nozzle tip and cooling substrate, a fast response solenoid valve (Parker) and a straight tube nozzle. The nozzle is made of stainless steel with an inner diameter of 0.38 mm and length of 40 mm, which resembles those used in clinical surgeries. Cryogen R404A (Dupont) is used to form the pulse spray due to its superior cooling capacity compared with R134a to protect the epidermis [6,25]. The spurt duration of the pulse spray (usually less than 100 ms in clinical surgeries) can be accurately controlled since the response time of

Z.-F. Zhou et al. / International Journal of Heat and Mass Transfer 100 (2016) 1–10

3

2.2. Temperature measurements

Fig. 1. Schematic of the cryogen spray cooling system.

the solenoid valve is less than 1 ms. An epoxy resin substrate (50 mm
50 mm
5 mm) is employed as the skin phantom, since its thermal properties are similar to those of human skin and so it has been extensively used previously by many researchers [17,18,26]. A DAQ board (NI: M6251) is used to acquire the temperature signal and generate the pulse voltage to control the solenoid valve. In this study, the spray distance and spurt duration are 30 mm and 20 ms, respectively.

Two different methods are employed to measure the surface temperature during and following the pulse spray cooling. Fig. 2 shows a schematic of these two methods: the thin film thermocouple (TFTC) and the fine thermocouple covered by aluminum foil. The TFTC of type T with a depth of 2 lm is deposited directly onto the surface of the epoxy resin using the Magnetron technique. It has perfect contact with the underlying substrate and a fast response time (1 ls). More information about TFTC can be found in our previous work [13]. The other surface temperature measurement method was first introduced by Aguilar et al., with a response time in the order of a few milliseconds [19]. As can be seen from Fig. 2(b), it consists of a fine thermocouple of type T (100 lm bead diameter), located beneath a thin layer of aluminum foil (10 lm), which is positioned on the surface of the epoxy resin substrate. A thermal paste is placed between the aluminum foil and the epoxy resin, ensuring good thermal contact and also providing mechanical support. Thus, there are totally three layers in this method, among which the first is the aluminum layer (L1) with depth of 10 lm, the second is the thermocouple layer (L2) with depth of 100 lm and the last is the epoxy resin layer (L3) with depth of 5 mm. The thermal properties of each layer are listed in Table 1. Since the TFTC has an extremely thin depth of 2 lm and is in direct contact with the liquid film on the substrate, this may be regarded as the direct surface temperature measurement method. On the contrary, the latter method for measuring surface temperature is far more complicated than the TFTC measurement method, with a larger thermocouple and an auxiliary aluminumfoil layer. The liquid film first comes into contact with the aluminum foil before the heat is conducted to the thermocouple

(a) Direct measurement of surface temperature with the thin film thermocouple

(b) Indirect measurement of surface temperature with the fine thermocouple covered by aluminum foil Fig. 2. Schematic of the two different surface temperature strategies.

4

Z.-F. Zhou et al. / International Journal of Heat and Mass Transfer 100 (2016) 1–10 Table 1 Thermal properties of each layer. Property

Aluminum

T thermocouple (average of copper and constantan)

Epoxy resin

k (W m1 K1) q (kg m3) cp (J kg1 K1) a (m2 s1)

236 2710 902 9.65
105

210 8910 400 5.90
105

0.14 1019 1631 0.843
107

and the epoxy resin substrate. Thus it belongs to the indirect measurement method of the surface temperature. It is noted that there exists large difference of thermocouple sizes between two measurement methods. However, the previous study indicated that the thermocouple with 40 lm only introduced a maximum temperature disturbance of 0.5 °C. Errors in the prediction of surface heat flux were minimal by using such small thermocouple [18]. Therefore, it can be inferred that the errors is minimal and the effect of our large thermocouple size of 100 lm (2.5 times as large as 40 lm) on the accuracy is negligible for such linear heat conduction problem once the appropriate algorithm is chosen. In the experiments, only the center temperature of the impingement area on the substrate surface is measured under both the direct and indirect measurement methods. 3. Mathematical description of the IHCP The depth of interest (500 lm) is far smaller than the dispersion of the cryogen spray (in the order of 20 mm) [20,23,24]. In addition, the previous studies on radial and temporal variation in surface heat transfer during cryogen spray cooling by Franco et al. [23] and Wang et al. [24] indicated that there always existed a subregion of uniform cooling with a radius of 2 mm in the cooling central area where the temperature and heat flux did not changes in the lateral direction. Therefore, the heat transfer in the single or multilayer slab (shown in Fig. 2) can be treated as a onedimensional heat conduction problem along the depth direction (x direction). Actually, most of the previous studies treated this heat transfer problem as 1D heat conduction problem [6,13,17– 21,25,26] during pulsed cryogen spray cooling in laser dermatology application. The dynamic temperature distribution can be described as: 2

@ T i ðx; tÞ 1 @T i ðx; tÞ ¼ @x2 ai @t

ð1Þ

where a is the thermal diffusivity, a = k/qc, and i denote the parameters in the i layer. The initial temperature is T0 in each region. The epoxy resin can be treated as a semi-infinite thick body and all of the layers are in perfect thermal contact. The thermal properties of each material remain unchanged during and following cryogen spray cooling. Then, the boundary condition and initial condition are:


i ðx;tÞ
ki @T @x


¼ kiþ1

x¼xi




@T iþ1 ðx;tÞ

@x x¼xiþ1

Tðx; tÞjx¼xi ¼ T iþ1 ðx; tÞjx¼xiþ1
@T 1 ðx; tÞ

k1 @x


¼ qðt Þ

9 = ;

ð2aÞ

ð2bÞ

x¼0

Tðx; 0Þ ¼ T 0

ð2cÞ

In IHCPs, the boundary condition q(t) is unknown. The objective is to estimate the transient surface heat flux q(t), based on the temperature measurement data Y(t) at x = c location

T c ðtÞ ¼ YðtÞ

ð3Þ

According to above instruction of two temperature measurement methods, the values of c are 0 and 10 lm, since the measured temperature are the real surface temperature and the one measured beneath the aluminum layer with thickness of 10 lm, respectively. 4. Algorithms for estimation of heat flux 4.1. Duhamel’s theorem In linear heat conduction problems, Duhamel’s theorem provides a convenient way of obtaining a solution with timedependent boundary conditions. When applying Duhamel’s theorem, the measured internal temperature Y(t) is introduced into the calculation of heat flux directly as the surface temperature of the epoxy resin substrate. Then, a one-dimensional, direct heat conduction problem for the two surface temperature measurement strategies may be solved to obtain the temperature distribution within the epoxy resin substrate [23,27]

Tðx; tÞ ¼ f ð0Þuðx; tÞ þ

Z

t

0

uðx; t  sÞ

dTðsÞ ds ds

ð4Þ

where u(x,t) is the temperature distribution response function of the substrate to a unit step (unit step function) in the surface temperature; f(0) and T(s) are the initial and time-varying surface temperature. The unit step function for a semi-infinite planar solid is [27]

uðx; tÞ ¼ 1  erf



x pffiffiffiffiffi 2 at



ð5Þ

Using Eq. (4), the temperature gradient at the surface can be solved and, subsequently, the surface heat flux from Fourier’s law is

qðtÞ ¼ k1



Z t @T

@ uðx; t  sÞ

dT ¼ k 1
ds ds @x
x¼0 @x 0 x¼0

ð6Þ

Substituting Eqs. (5) into (6), the surface heat flux, q(t), can be written as

qðtÞ ¼

rffiffiffiffiffiffiffiffi Z kqc

p

t 0

1 dT pffiffiffiffiffiffiffiffiffiffiffi ds t  s ds

ð7Þ

4.2. New method based on Duhamel’s theorem When applying Duhamel’s theorem to estimate surface heat flux from temperature data measured indirectly, the simplification that the experimental data equals the surface temperature somewhat changes the physics of the problem, where neither the heat dissipating capacity of the materials above the thermocouple nor the thermocouple itself is taken into account. This may cause significant errors in cases where the surface temperature changes rapidly. Therefore, a new method is proposed to deal with this problem based on Duhamel’s theorem, rather than directly

Z.-F. Zhou et al. / International Journal of Heat and Mass Transfer 100 (2016) 1–10

estimate surface heat flux from the measured temperature data with utilization of Duhamel’s theorem. Firstly, the real surface temperature should be calculated from the experimental data measured inside the substrate at x = c location, which can be expressed by a piecewise constant function of time as

f ¼ f ði  DtÞ; i ¼ 1; 2; . . . ; n

ð8Þ

where f represents the surface temperature at the time of iDt. Then, the temperature at x = c location can be represented from Duhamel’s theorem as [27]

T c ðn; DtÞ ¼ f ð0Þuðx; n  Dt Þ þ ½f ð1  Dt Þ  f ð0Þuðx; ðn  1Þ  Dt Þ ð9Þ

Eq. (9) may be rearranged and written in expanded matrix form as

2

DuðDtÞ 0 6 Duð2DtÞ D u ðDtÞ 6 6 .. .. 6 4 . . DuðnDtÞ Duððn  1ÞDtÞ 3 2 T c ð0  DtÞ 6 T c ð1  DtÞ 7 7 6 7 ¼6 .. 7 6 5 4 . T c ðn  DtÞ



0

... .. .

0 0

32 76 76 76 76 54

   DuðDtÞ

f ð0  DtÞ f ð1  DtÞ .. .

3 7 7 7 7 5

f ððn  1Þ  DtÞ

where Du(iDt) represents u (x, iDt)  u (x, (i  1)Dt). The first matrix on the left-hand side of the equation is the temperature transformation matrix, which converts the surface temperature into the internal temperature. Then, the surface temperature history can be obtained by multiplying the inverse temperature transformation matrix on both sides of the equation from the measured internal temperature data. Subsequently, the surface heat flux can be solved by Eq. (7). From the above discussion, the new method extends the application of Duhamel’s theorem which is originally applicable to heat conduction problems with known time-dependent boundary conditions. 4.3. Sequential function specification method

SM ¼

Y kþr1  T i;kþr1

ð11Þ

Several functional forms for q(t) from tk to tk+r1 has been proposed. The simplest one is that q(t) is a constant

ð12Þ

Then, the temperature distribution is represented as a function of surface heat flux qk and the temperature field is expanded in a ^k as Taylor series about a known value of surface heat flux q

^k Þ þ Z i;k ðqk  q ^k Þ þ    T i;k ðqk Þ ¼ T i;k ðq

^k Þ Z ki;r1 Y kþr1  T i;kþr1 ðq PR k 2 i¼1 X i;r1

The dynamic system and heat conduction problem share a key feature: both are described by a set of differential equations [11]. Therefore, the transfer functions, which are used to characterize the relationship between the input and output of a dynamic system, can also be used to analyze the linear heat conduction problems. Thus, the heat flux is treated as the input of the system, and the temperature field as the response. The Laplace transform is frequently used when analyzing the dynamic system to provide the mathematical relationship between the input and output of the dynamic system. The relation between input and output, in a linear dynamic system (Fig. 3), is given by the multiplication in the complex variable s domain,

ð16Þ

Applying the convolution theorem to (16), the inverse integral y (t) in terms of Y(s) is given by the convolution between h(t) and x (t), represented symbolically by the operator (⁄). The relationship between the input and output in the time domain is provided by the convolution integral

yðt Þ ¼ L1 ½Y ðsÞ ¼ L1 ½HðsÞ  X ðsÞ ¼ hðtÞ  xðt Þ Z t ¼ hðt  sÞxðsÞds

@T i;k ðqk Þ @qk

ð17Þ

0

The transfer function of a system is defined as the Laplace transform of the impulse respond of the system H(s) = L[h(t)]. The integral solution to the problem given by Eqs. (1) and (2) is obtained by Green’s function [12]

a k

Z 0

t

Gc;1 ðx; tj0; sÞqðsÞds

ð18Þ

where G is the green function. By subtracting the initial temperature T0 on both sides of the Eq. (18), the temperature distribution can then be represented as

Z hc ðt Þ ¼

t

hc ðt  sÞqðsÞds

ð19Þ

where hc(t) = Tc  T0, hc(t–s) = (a/k)Gc,1(x,t|0,s). As expressed in Eq. (19), the relation between surface heat flux (input) and temperature distribution (respond) in the complex variable s domain is given by the multiplication expressed in Eq. (16). Thus, in terms of heat flux/ temperature pairs, we have yields

L½hc ðtÞ ¼ L½hc ðtÞ  qðtÞ ! Hc ðsÞ ¼ Hc ðsÞ  Q ðsÞ Or, we can write

ð13Þ

where Zi,k is the sensitivity coefficient defined by

Z i;k 

ð15Þ

In this paper, the value of r is selected based on the residual principle [28], which are 4 and 3 for the direct and indirect temperature measurement methods, respectively. The detailed information about the SFS method can be found in references [8,18].

0

i¼1

qk ¼ qkþ1 ¼    ¼ qkþr

i¼1

T c ðtÞ ¼ T 0 þ

The SFS method is one of the most widely-used technologies for solving inverse heat conduction problems. The least-square method is used to minimize the error between the measured temperature Yk and estimated temperature Tk for the current time and r future time steps [8]

2

PR  ^k þ qk ¼ q

Y ðsÞ ¼ HðsÞ  X ðsÞ ð10Þ

R X 

The solution for the estimated surface heat flux at time tk can be obtained by minimizing Eq. (11) with respect to qk

4.4. Transfer function method

þ ½f ð2  Dt Þ  f ð1  DtÞuðx; ðn  2Þ  Dt Þ þ    þ ½f ðn  DtÞ  f ððn  1Þ  Dt Þuðx; 0Þ

5

ð14Þ Fig. 3. Dynamic system with one input and one output.

ð20Þ

6

Z.-F. Zhou et al. / International Journal of Heat and Mass Transfer 100 (2016) 1–10

Q ðsÞ ¼

1 Hc ðsÞ Hc ðsÞ

ð21Þ

Eq. (21) in the time domain is equivalent to the deconvolution

qðt Þ ¼ L1





1 1  Hc ðsÞ ¼ L1  hc ðt Þ Hc ðsÞ Hc ðsÞ

ð22Þ

Therefore, the solution to the problem is the estimation of the surface heat flux, of which the input is the measured temperature and the transfer function is given by 1/Hc(s). Using Eq. (22), we can estimate the heat flux simultaneously for all time. The inverse Laplace transform of 1/Hc(s) can be determined analytically in finite and single-layer heat conduction problems. In other cases, such as multiply-layers problems, however, this will be difficult to calculate using the major inverse Laplace transform numerical approaches, such as the fast Fourier transform [29], Dubner method [30], Crump method [31] and Stehfest method [32], the first three of which require limited finite singularity in 1/Hc(s) and the last of which is incapable of providing stable inverse solutions near the origin. Thus, in order to apply this method to solve the multiply-layers IHCP, we proposed a method to broaden the application of this transfer function method. Firstly, we rewrite Eq. (22) as

Z

t

qðt Þ ¼ 0

hc ðt  sÞf ðsÞds

above, it is clear that the transfer function f(t) is determined directly by the Green’s function Gc,1(x,t|0,s), which depends only on the physical property and dimension parameters of the composite medium [12,27]. Therefore, ~ F is the same for a given heat conduction problem even when the boundary condition differs from the practical situation. Then, a simple auxiliary problem is formulated where q(t) in Eq. (22) is now set as a constant of 1. The hc,i in the lower triangular matrix can be calculated numerically using an implicit central finite difference technique. Therefore, the appropriate vector ~ F can be obtained by simply multiplying the inverse matrix on both sides of Eq. (26). Then, the surface heat flux can be determined by substituting the numerical solution of the inverse Laplace transform of L1[1/[Hc(s)] in Eq. (22). 5. Computational results using simulated data A hypothetical triangular pulse heat flux is commonly used to examine the accuracy and sensitivity to random noises of the algorithms [8,11], as shown in Fig. 4. This given heat flux is first set as

600

ð23Þ

500

where f(s) = L1[1/[H(s)]. Then, the discrete convolution of Eq. (23) can be represented as

þ hc ð0Þ  f ðiDt Þ  Dt

2

ð24Þ

The expanded matrix form of Eq. (24) can be written as

2

q1

3

2

hc;1

6 q 7 6 hc;2 6 27 6 6 7¼6 6 .. 7 6 .. 4 . 5 4 . hc;n

qn

0

...

hc;1 .. .

... .. .

hc;n1

32

0

f1

3

6 7 0 7 76 f 2 7 76 7 76 . 7 0 54 .. 5

   hc;1

400

q (kW/m )

qðiDt Þ ¼ hc ðiDtÞ  f ð0Þ  Dt þ hc ðði  1ÞDt Þ  f ð1Dt Þ  Dt þ   

300 200

Hypothetical Duhamel SFSM Transfer function

100 ð25Þ

0

fn

0

Note that qi = q(iDt), hc,i = hc(iDt), fi = f(iDt)Dt. In a more compact form, Eq. (25) is

~¼H ~c  ~ F Q

80

t (ms)

120

160

200

(a) Under the direct surface temperature measurement method with TFTC

ð26Þ

hc,i can be calculated from the measuring data in the lower triangular matrix. In order to obtain the history of the surface heat flux, we still need to identify the fi for the given problem. From the analysis

600 500 400

2

q (kW/m )

600

40

500

Hypothetical Duhamel SFSM Transfer function New method

200

q (kW/m

2

)

400

300

300

100

200

0

100 0 0

40

80

t (ms)

120

160

Fig. 4. Hypothetical triangular pulse heat flux as a function of time.

200

0

40

80

t (ms)

120

160

200

(b) Under the indirect surface temperature measurement method with fine thermocouple covered by aluminum foil Fig. 5. Comparison between the hypothetical triangular pulse heat flux and estimated heat fluxes through different algorithms without random errors.

7

Z.-F. Zhou et al. / International Journal of Heat and Mass Transfer 100 (2016) 1–10

triangular pulse heat flux. This indicates that all of the algorithms can predict heat flux accurately under the direct surface temperature method. When using the indirect measuring method (Fig. 5(b)), the estimated heat flux obtained from Duhamel’s theorem changes nonlinearly over time. Furthermore, it can hardly match the exact heat flux well. Thus, it is unsuitable for estimating the surface heat flux in the indirect surface temperature measurement. The result predicted by the SFS method only matches the given heat flux well at the increasing period, whereas the deviation starts to emerge after the maximum heat flux. This indicates that the SFS method can only provide accurate predictive results during the initial period. In comparison, both the transfer function and new methods can provide linear and accurate results that match the given heat flux exactly.

600 500

2

q (kW/m )

400 300 200 Hypothetical Duhamel SFSM Transfer function

100 0 0

40

80

120

t (ms)

160

200 5.2. Predicted temperature with noise

(a) Under the direct surface temperature measurement method with TFTC

The addition of random noise is a procedure for analyzing the ability of an algorithm to deal with temperature data that are suffering disturbance due to experimental errors. In most cases, a small perturbation in the experimental signals can be amplified largely in the IHCPs. In order to analyze the effect of measurement noise on the estimated heat flux quantitatively, we define the standard deviation of the estimated surface heat flux due to the random noise fq as [8]

600 500

2

q (kW/m )

400



300 fq ðr Þ ¼

200

Hypothetical Duhamel SFSM Transfer function New method

100 0 0

40

80

t (ms)

120

160

200

(b) Under the indirect surface temperature measurement method with the fine thermocouple covered by aluminum foil Fig. 6. Comparison between the hypothetical triangular pulse heat flux and the estimated heat fluxes by different algorithms with random errors.

the surface boundary condition, then the temperature at the measured point corresponding to the two different measurement methods is obtained by solving the direct heat conduction problem (Eqs. (1) and (2)). New surface heat flux can be predicted based on the calculated temperature through different algorithms. Two test cases will be analyzed. The first considers the calculated temperature without the addition of noise. The second adds random noise to the calculated temperature and investigates the influence of this on the new predicted surface heat flux. The noise used here follows a normal distribution with the zero mean, and has a standard deviation of 0.08 °C. The results of the two test cases are shown in Figs. 5 and 6.

5.1. Predicted temperature without noise Fig. 5 presents the estimated surface heat flux from the temperature without noise compared with the given heat flux under the two temperature measurement methods. In Fig. 5(a), it can be seen that the heat fluxes calculated by Duhamel’s theorem, SFS and the transfer function method all agree perfectly with the hypothetical

1=2 1 Xn 0 ðq  q Þ e;i e;i i¼1 n1

ð27Þ

where the qe,i and q0 e,i denote the estimated surface heat flux both with and without the addition of noise at time ti. Fig. 6 shows the new estimated heat fluxes using different algorithms with random noise under the direct and indirect surface temperature measurement methods. It is clear that the random noise exerts little influence on the estimated surface heat flux in the direct surface temperature measurement with TFTC. However, for the indirect surface temperature measurement with the fine thermocouple covered by aluminum foil, the random noise plays a different role in affecting the estimated surface heat flux. The influence of the random noise on the predicted heat flux remains negligible for Duhamel’s theorem because of its nature of solving the direct heat conduction problem. The estimated heat flux using the SFS method undergoes considerable fluctuation during the whole time domain, which indicates that the SFS method is highly sensitive to random noise. The predicted results by the transfer function and new methods display a virtually identical response to the random noise. As deconvolution is performed in both algorithms, the random noise is not amplified but linearly transmitted to the estimated heat flux. The quantitative perturbation of the estimated heat flux can be seen more clearly in Table 2. They are all quite small for different algorithms in the case of directly measuring the surface temperature. However, the deviations in the estimated heat flux are far larger under the indirect measurement method than those before, except when Duhamel’s theorem is applied. Furthermore, the fq of the transfer function and new methods are far less than that

Table 2 Standard deviation of estimated heat flux in the test case with random errors.

TFTC fq (kW/m2) Fine thermocouple fq (kW/m2)

Duhamel

SFSM

Transfer function

New method

4.2145 4.2106

1.6700 61.9621

3.1283 21.5485

21.4632

8

Z.-F. Zhou et al. / International Journal of Heat and Mass Transfer 100 (2016) 1–10

where Yk is the measured temperature, and Tkc is the temperature calculated from the estimated heat flux qk through solving the direct heat conduction problem with Eqs. (1) and (2). The standard deviation ST provides an indirect criterion for judging the accuracy of the algorithm to predict the surface heat flux. The smaller the value of ST, the more able the algorithm to estimate heat flux. The residual of the temperature is the difference between the measured temperature and the calculated temperature at each time point, given as Residual = Yk  Tkc .

6.1. Heat flux based on the direct surface temperature measurement method

Fig. 7. Variation in the surface temperature as measured by the TFTC and the fine thermocouple covered by aluminum foil.

of the SFS method, which implies that the former two approaches are far more capable of inhibiting random errors in temperature data. 6. Results and discussion using measured temperature data Fig. 7 shows the variation in surface temperature as a function of time, as measured by the direct measurement method with the TFTC and the indirect measurement method with the fine thermocouple covered by aluminum foil. The temperature first decreases rapidly as the droplets of R404A impinge on the substrate surface, following a relatively slow change when the liquid film forms on the surface. It begins to resume the ambient temperature as the liquid film completely evaporates. It is notable that the surface temperature measured by the TFTC decreases far faster than that by the indirect measurement method, which shows a definite delay. Since the exact heat flux on the substrate surface can hardly be known in this CSC experiment, an alternative expression of the standard deviation between the estimated and measured temperature is employed to assess the quality of different algorithms, given as [15]

ST ðrÞ ¼

2 1 XK ðY k  T kc Þ k¼1 K 1

1=2 ð28Þ

800

15 z = 30 mm t = 20 ms Duhamel SFSM Transfer function

400

200

0

10

20

t (ms)

30

40

z = 30 mm = 20 ms Duhamel SFSM Transfer function

10

Residual (oC)

2

q (kW/m )

600

0

When the surface temperature is measured by TFTC, Fig. 8 presents the variation in surface heat flux as a function of time predicted by different algorithms of Duhamel’s theorem, the SFS method and transfer function method. It is clear that the three algorithms predict very similar results for heat flux. They first increase rapidly to a maximum value qmax, then drop quickly to a certain value (350 kW/m2), before gradually decreasing to zero. This variation is consistent with the history of measured surface temperature, shown in Fig. 7. Note that the estimated result by the SFS method increases faster and has a slightly greater maximum heat flux than that obtained through employing the other two algorithms. Fig. 9 presents the residual temperature between the measured and predicted temperature under the direct surface temperature measurement method with TFTC. It can also be seen that, for each algorithm, the maximum temperature residual rmax occurs approximately at the same time as qmax happens (tmax). Then, the residual drops quickly to zero and remains virtually unchanged as t increases. One can see that the largest residual temperature exists during the short period around tmax, indicating that the predicted heat flux contains non-negligible errors when the SFS method is used. Table 3 shows the quantitative standard deviation of the surface temperature during the whole time domain. Corresponding with Fig. 9, the SFS method has the largest ST, whereas the transfer function method has the smallest ST. The middle one is that predicted by Duhamel’s theorem. Therefore, the SFS method is less accurate than Duhamel’s theorem and the transfer function method by comparing the standard deviation and residual temperature for the direct surface temperature measurement.

50

Fig. 8. Surface heat flux predicted by different algorithms under the direct surface temperature measurement method with TFTC.

5

0

-5

0

10

20

t (ms)

30

40

50

Fig. 9. Residual temperature between the predicted and measured temperature under the direct surface temperature measurement method with TFTC.

Z.-F. Zhou et al. / International Journal of Heat and Mass Transfer 100 (2016) 1–10 Table 3 Standard deviation between the estimated and measured temperature by the TFTC.

ST(r) (°C)

SFSM

Duhamel’s theorem

Transfer function

1.47829

0.479494

0.119881

800 z = 30 mm t = 20ms Duhamel SFSM Transfer function New method

2

q (kW/m )

600

400

200

0

0

10

20

t (ms)

30

40

50

Fig. 10. Surface heat flux predicted by different algorithms under the indirect surface temperature measurement method with a fine thermocouple.

10 5

Residual (oC)

0 -5 z = 30 mm = 20ms Duhamel SFSM Transfer function New method

-10 -15 -20 -25

0

10

20

t (ms)

30

40

50

Fig. 11. Residual of temperature between the predicted and measured temperature under the indirect surface temperature measurement method with a fine thermocouple.

6.2. Heat flux based on the indirect surface temperature measurement method Fig. 10 presents the heat flux predicted by different algorithms under the indirect surface temperature measurement method with the fine thermocouple. As mentioned above, Duhamel’s theorem may change the physics of the problem and is inappropriate for predicting surface heat flux directly when the temperature is measured inside the substrate. Thus, the value of qmax and the

Table 4 Standard deviation between the estimated and measured temperature by the fine thermocouple.

ST(r) (°C)

SFSM

Duhamel’s theorem

Transfer function

New method

2.8809

7.4666

6.9201e5

2.5570e2

9

magnitude of its increase rate by Duhamel’s theorem are both far lower than is the case with other algorithms, as shown in Fig. 10. It seems that the heat flux is significantly delayed and reduced. For other results of surface heat flux, the two predicted by the transfer function and new methods have almost the same varying history, with a larger qmax than that predicted by the SFS method. Nevertheless, these three results are relatively similar, compared with that predicted by Duhamel’s theorem. Correspondingly, Fig. 11 presents the residual temperature between the measured and predicted temperature. It is clear that Duhamel’s theorem generates the largest residual temperature, and the same phenomenon is also observed for the SFS theorem. Quite different from Fig. 9, however, the residual temperature always exists in the whole time domain. For the transfer function and new methods, there exists virtually no residual temperature. Furthermore, from Table 4, which shows the standard deviation in temperature for different algorithms, Duhamel’s theorem and the SFS method generate far larger deviations than is the case with the other two algorithms, which are negligible. This implies that both the transfer function and new methods are more appropriate for predicting surface heat flux under the indirect surface temperature measurement method. For all cases, Duhamel’s theorem can only calculate surface heat flux accurately when the surface temperature is measured directly. Otherwise, a wide deviation might be caused, so it should be not suitable for predicting the surface heat flux accurately. The estimated results under the SFS method of around tmax suffer partial distortion. The future information used to obtain the surface heat flux is effective in minimizing the effects of random noise, under the hypothesis that the surface heat flux remains constant over r future time steps. This is proved to be reasonable when the heat flux changes slowly. During a pulsed cryogen spurt, however, the surface heat flux increases dramatically once the cryogen droplets impinge on the surface, so that even a few time steps can produce substantial changes, which makes the assumption quite inappropriate and causes significant deviation in the estimated results. Therefore, the SFS method is slightly less accurate than the transfer function and new methods in our cases. These two algorithms produce virtually identical results, which have far smaller deviation of temperature compared to others. In addition, they have the ability to inhibit the measured noise.

7. Conclusions Several algorithms are analyzed and compared in predicting time-varying surface heat flux due to ultra-fast surface cooling. A new method is proposed based on Duhamel’s theorem to get the accurate results through the transformation of internal temperature into surface temperature, when the indirect surface temperature measurement method is used. A convenient method is proposed that involves solving an auxiliary problem to obtain the transfer function. As a result, the transfer function method can be implemented to solve the multilayer, semi-infinite IHCPs. A hypothetical triangular pulse heat flux is first employed to analyze the accuracy and sensitivity to noise of the algorithms under different measurement methods. The estimated result of Duhamel’s theorem widely deviates from the given heat flux under the indirect temperature measurement method, whereas the transfer function and new methods all provide the exact estimated heat flux under both measuring conditions. The random noise has little effect on the solution when Duhamel’s theorem is applied in all cases. The result of the SFS method suffers from random errors, being heavily based on the indirect surface temperature measurements. The noises are linearly transmitted to the results without being amplified and the estimated

10

Z.-F. Zhou et al. / International Journal of Heat and Mass Transfer 100 (2016) 1–10

heat fluxes are stable in all measuring approaches for the transfer function and new methods. Then, the algorithms are implemented to calculate the surface heat flux based on the measured surface temperature of the epoxy resin that is exposed to pulse cryogen spray cooling using two different measurement methods. In the condition that the surface temperature can be accurately measured by the TFTC, excellent estimated resulted can be obtained through employing Duhamel’s theorem. However, it is inappropriate to predict surface heat flux in cases where the surface temperature is measured inside the substrate. The SFS method can produce acceptable results under both measurement conditions, except in the time region where the heat flux changes extremely rapidly. The transfer function and new methods produce highly similar and excellent results with a negligible standard deviation of temperature for both the direct and indirect surface temperature measurement methods. The transfer function method is more rigorous with regard to mathematics and may be extended to other inverse heat conduction problems conveniently. Acknowledgements This work was supported by the Natural Science Foundation of China (51336006, 51406151), the International Science & Technology Cooperation Plan of Shaanxi Province (2013KW30-05), China Postdoctoral Science Foundation (2013M540749) and the Fundamental Research Funds for the Central Universities. References [1] S.-S. Hsieh, S.-Y. Luo, Droplet impact dynamics and transient heat transfer of a micro spray system for power electronics devices, Int. J. Heat Mass Transfer 92 (2016) 190–205. [2] R. Dou, Z. Wen, G. Zhou, Heat transfer characteristics of water spray impinging on high temperature stainless steel plate with finite thickness, Int. J. Heat Mass Transfer 90 (2015) 376–387. [3] P. Bhattacharya, A. Samanta, S. Chakraborty, Spray evaporative cooling to achieve ultra fast cooling in runout table, Int. J. Therm. Sci. 48 (9) (2009) 1741– 1747. [4] S. Chen, J. Liu, X. Liu, Y. Hou, An experimental comparison of heat transfer characteristic between R134-a and R22 in spray cooling, Exp. Thermal Fluid Sci. 66 (2015) 206–212. [5] J.S. Nelson, T.E. Milner, B. Anvari, B.S. Tanenbaum, S. Kimel, L.O. Svaasand, S.L. Jacques, Dynamic epidermal cooling during pulsed laser treatment of portwine stain. A new methodology with preliminary clinical evaluation, Arch. Dermatol. 131 (6) (1995) 695–700. [6] Z. Zhou, B. Chen, Y. Wang, L. Guo, G. Wang, An experimental study on pulsed spray cooling with refrigerant R-404a in laser surgery, Appl. Therm. Eng. 39 (2012) 29–36. [7] D. Li, B. Chen, W.J. Wu, G.X. Wang, Y.L. He, Multi-scale modeling of tissue freezing during cryogen spray cooling with R134a, R407c and R404a, Appl. Therm. Eng. 73 (2) (2014) 1489–1500. [8] J.V. Beck, B. Blackwell, C.R.S. Clair Jr., Inverse Heat Conduction: Ill-Posed Problems, James Beck, 1985. [9] P.C. Hansen, Regularization tools: a Matlab package for analysis and solution of discrete ill-posed problems, Numer. Algorithms 6 (1) (1994) 1–35.

[10] O. Alifanov, Solution of an inverse problem of heat conduction by iteration methods, J. Eng. Phys. Thermophys. 26 (4) (1974) 471–476. [11] A.P. Fernandes, M.B. dos Santos, G. Guimarães, An analytical transfer function method to solve inverse heat conduction problems, Appl. Math. Model. 39 (22) (2015) 6897–6914. [12] K.D. Cole, J.V. Beck, A. Haji-Sheikh, B. Litkouhi, Heat Conduction Using Green’s Functions, Taylor & Francis, 2010. [13] Z.-f. Zhou, B. Chen, R. Wang, F.-l. Bai, G.-x. Wang, Coupling effect of hypobaric pressure and spray distance on heat transfer dynamics of R134a pulsed flashing spray cooling, Exp. Thermal Fluid Sci. 70 (2016) 96–104. [14] G. Aguilar, B. Choi, M. Broekgaarden, O. Yang, B. Yang, P. Ghasri, J.K. Chen, R. Bezemer, J.S. Nelson, A.M. van Drooge, A. Wolkerstorfer, K.M. Kelly, M. Heger, An overview of three promising mechanical, optical, and biochemical engineering approaches to improve selective photothermolysis of refractory port wine stains, Ann. Biomed. Eng. 40 (2) (2012) 486–506. [15] B.M. Pikkula, J.H. Torres, J.W. Tunnell, B. Anvari, Cryogen spray cooling: effects of droplet size and spray density on heat removal, Lasers Surg. Med. 28 (2) (2001) 103–112. [16] B.M. Pikkula, J.W. Tunnell, B. Anvari, Methodology for characterizing heat removal mechanism in human skin during cryogen spray cooling, Ann. Biomed. Eng. 31 (5) (2003) 493–504. [17] J.H. Torres, J.W. Tunnell, B.M. Pikkula, B. Anvari, An analysis of heat removal during cryogen spray cooling and effects of simultaneous airflow application, Lasers Surg. Med. 28 (5) (2001) 477–486. [18] J.W. Tunnell, J.H. Torres, B. Anvari, Methodology for estimation of timedependent surface heat flux due to cryogen spray cooling, Ann. Biomed. Eng. 30 (1) (2002) 19–33. [19] G. Aguilar, G.X. Wang, J.S. Nelson, Effect of spurt duration on the heat transfer dynamics during cryogen spray cooling, Phys. Med. Biol. 48 (14) (2003) 2169– 2181. [20] W. Jia, G. Aguilar, G.X. Wang, J.S. Nelson, Heat-transfer dynamics during cryogen spray cooling of substrate at different initial temperatures, Phys. Med. Biol. 49 (23) (2004) 5295–5308. [21] G. Aguilar, H. Vu, J.S. Nelson, Influence of angle between the nozzle and skin surface on the heat flux and overall heat extraction during cryogen spray cooling, Phys. Med. Biol. 49 (10) (2004) N147–N153. [22] W. Franco, J. Liu, R. Romero-Mendez, W. Jia, J.S. Nelson, G. Aguilar, Extent of lateral epidermal protection afforded by a cryogen spray against laser irradiation, Lasers Surg. Med. 39 (5) (2007) 414–421. [23] W. Franco, J. Liu, G.X. Wang, J.S. Nelson, G. Aguilar, Radial and temporal variations in surface heat transfer during cryogen spray cooling, Phys. Med. Biol. 50 (2) (2005) 387–397. [24] R. Wang, Z. Zhou, B. Chen, F. Bai, G. Wang, Surface heat transfer characteristics of R404A pulsed spray cooling with an expansion-chambered nozzle for laser dermatology, Int. J. Refrig 60 (2015) 206–216. [25] T. Dai, M.A. Yaseen, P. Diagaradjane, D.W. Chang, B. Anvari, Comparative study of cryogen spray cooling with R-134a and R-404a: implications for laser treatment of dark human skin, J. Biomed. Opt. 11 (4) (2006) 041116. [26] J.H. Torres, J.S. Nelson, B.S. Tanenbaum, T.E. Milner, D.M. Goodman, B. Anvari, Estimation of internal skin temperatures in response to cryogen spray cooling: implications for laser therapy of port wine stains, IEEE J. Sel. Top. Quantum Electron. 5 (4) (1999) 1058–1065. [27] D.W. Hahn, M.N. Ozisik, Heat Conduction, John Wiley & Sons, 2012. [28] J. Beck, B. Blackwell, A. Haji-Sheikh, Comparison of some inverse heat conduction methods using experimental data, Int. J. Heat Mass Transfer 39 (17) (1996) 3649–3657. [29] H. Baher, The fast Fourier transform and its applications, Signal Process. Integr. Circuits (1990) 149–191. [30] H. Dubner, J. Abate, Numerical inversion of Laplace transforms by relating them to the finite Fourier cosine transform, J. ACM (JACM) 15 (1) (1968) 115– 123. [31] K.S. Crump, Numerical inversion of Laplace transforms using a Fourier series approximation, J. ACM (JACM) 23 (1) (1976) 89–96. [32] H. Stehfest, Algorithm 368: numerical inversion of Laplace transforms [D5], Commun. ACM 13 (1) (1970) 47–49.