The Fundamental and Application of Surface Heat Flux Estimation by Inverse Method in Cryogen Spray Cooling

Abstract

Surface heat flux is an important parameter in various industrial applications, which is often estimated based on measured temperature by solving inverse heat conduction problem (IHCP). In this chapter, the available IHCP methods including sequential function specification (SFS), transfer function (TF) and Duhamel’s theorem were compared, taking the example of surface heat flux estimation during spray cooling. The Duhamel’s theorem was improved to solve 1D multi-layer ICHP. Considering the significant nonuniformity of heat transfer, the 2D filter solution method was proposed to estimate surface heat flux for 2D multi-layer mediums. The maximum heat flux calculated by the 1D method was underestimated by 60% than that calculated by 2D filter solution, indicating that the lateral heat transfer cannot be ignored. The cooling performances based on 2D filter solution demonstrated that substituting the environment friendly R1234yf for R134a can remarkably reduce global warming potential to <1, but its cooling capacity is insufficient. The effective heat flux of R1234yf can be enhanced by 18.8% by reducing the nozzle diameter and decreasing the back pressure, providing the theoretical basis for the clinical potential substitution of R1234yf with low global warming potential (GWP) for commercial R134a with high GWP in laser dermatology.

1. Introduction

Heat flux is an important parameter to characterize heat transfer performance in many industrial applications, such as thermal protection of space shuttles [1], thermal management of electronic devices [2], metal heat treatment [3], maintenance of boilers [4] and nuclear reactors [5], spray cooling [6], geophysics [7], etc. Heat flux is often estimated by surface or internal temperature, which is also termed as inverse heat conduction problem (IHCP).

IHCPs are mathematically ill-posed, and a small error in temperature may significantly affect the accuracy of heat flux estimation [8, 9]. Several analytical and numerical methods have been proposed for the solution of IHCPs, such as sequential function specification (SFS) method, Tikhonov regularization (TR) method, transfer function (TF) method, Duhamel’s theory, etc. The SFS method is commonly used to solve IHCPs by minimizing the effect of random errors using temperature data at future time steps based on the least square method [10]. The TR method estimates all of the heat flux simultaneously for all time steps and is usually presented as whole time domain form, which often causes heavy computational load [11, 12]. The TF method analogizes the heat conduction problems to dynamic systems, in which heat flux is treated as the input of the system and the temperature profile as the response [13]. 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. Duhamel’s theorem is based on the principle of superposition and assumes that the substrate thermal response at t equals the total sum of what the substrate experienced in small steps prior to t [14]. This method is simple and widely used in the surface heat flux estimation with known surface temperature. However, its assumption that the internal temperature equals to the surface one for indirect temperature measurement often causes significant calculation errors. Thus, Duhamel’s theory needs to be improved.

Recently, a filter solution based on TR has attracted the interest of many researchers [9, 15, 16, 17, 18, 19], which minimizes the sum of the squares of the errors between estimated and measured temperatures and stabilized by Tikhonov regularization. This solution is expressed in a digital filter form, allowing for a real-time heat flux estimation, and has been used for heat flux estimation of directional flame thermometer [20]. This method demonstrates superiority when solving IHCPs with a complex geometry. However, it can only be used to solve 1D single-layer IHCPs. The multidimensional, multilayer IHCP has yet to be solved.

Although IHCPs have been extensively investigated with regard to various other applications, little work has been conducted related to surface heat flux estimation during cryogen spray cooling (CSC) in laser dermatology. Spray cooling is widely applied in metallurgy, electronics, power plant, and laser dermatology for vascular skin lesions [21, 22, 23, 24, 25, 26, 27], because of the superiorities of high power density, ultrafast cooling rate, uniformity of heat removal, and low fluid inventory. In the laser treatment of vascular skin lesions (e.g., port wine stain, PWS), CSC can be implemented to prevent unwanted thermal damage of the epidermis induced by high laser absorption of melanin. Different with traditional steady spray cooling, CSC is a highly transient process with several tens of milliseconds to avoid cold injury. The transient surface heat flux is crucial for cooling performance evaluation, which needs accurate heat flux estimation method and rapid temperature measurement with fast response and small lag, as well as damping of algorithm [28].

In transient CSC, two typical temperature measurements of fine thermocouple (FTC) and thin-film thermocouple (TFTC) are widely used to measure internal and surface temperature. Aguilar et al. [29] used the SFS method to estimate surface heat flux by internal temperature measured by a type-T FTC placed underneath a thin-layer aluminum foil, positioned on the top of epoxy resin surface to provide rapid heat transfer and mechanical support. Zhou et al. [28, 30, 31] measured time-dependent surface temperature by a type-T TFTC with thickness of 2 μm directly deposited onto the epoxy resin surface; this measurement accurately captured the temperature variation owing to its ultrafast thermal response (∼1.2 μs). Then, the surface heat flux was estimated by Duhamel’s theorem. However, TFTC cannot be used to measure the metal material temperature due to the electrical conductivity. Moreover, TFTC corrodes and oxidizes easily in high-temperature environments. Therefore, FTC measurement is widely used in many industries owing to its reliability and stability. Different with TFTC measurement with single-layer geometry, FTC measurement consists of three layers, namely, aluminum, thermal paste, and epoxy resin. Moreover, the radial and temporal surface temperature variations during CSC result in significant nonuniformity of the surface heat flux [14, 32]. Therefore, lateral heat transfer must be considered. For generality, 2D multilayer IHCPs need to be developed.

In this chapter, the SFS, TF, and Duhamel’s theorem methods for TFTC and FTC measurements were compared based on one hypothetical heat flux. Duhamel’s theorem was improved to increase the accuracy of surface heat flux when using the indirect three-layer FTC measurement. Afterwards, the 2D filter solution was proposed to calculate the surface heat flux for 2D multilayer mediums. Six hypothetical triangular pulse heat fluxes were used to examine the accuracy and sensitivity of the algorithm. Finally, the 2D filter solution was employed to calculate the surface heat flux to investigate the cooling performance. Using estimated heat flux as evaluation index, the possibility of substituting the commercial cryogen R134a with high GWP (1430) by environment friendly cryogen with low GWP (<1) was discussed.

2. 1D algorithms for heat flux estimation

2.1 Sequential function specification method

The SFS method is widely used to solve IHCPs by minimizing the error between the measured temperature Y_k and estimated temperature T_k for the current time and r future time steps based on the least square method [10]

$$S_M=\sum _{i=1}^R{\left(Y_{k+r-1}-T_{i,k+r-1}\right)}^2$$ E1

Several functional forms of q(t) from t_k to t _k + r−1 have been proposed. The simplest one is that q(t) is a constant

$$q_k=q_{k+1}=\cdots =q_{k+r}$$ E2

Then, the temperature distribution is represented as a function of surface heat flux q_k , and the temperature field is expanded in a Taylor series about a known value of surface heat flux ${\widehat{q}}_k$ as

$$T_{i,k}\left(q_k\right)=T_{i,k}\left({\widehat{q}}_k\right)+Z_{i,k}\left(q_k-{\widehat{q}}_k\right)+\cdots$$ E3

where Z_i,k is the sensitivity coefficient defined by

$$Z_{i,k}\equiv \frac{\partial T_{i,k}\left(q_k\right)}{\partial q_k}$$ E4

The solution for the estimated surface heat flux at time t_k can be obtained by minimizing Eq. (1) with respect to q_k :

$$q_k={\widehat{q}}_k+\frac{\sum _{i=1}^R\left[Y_{k+r-1}-T_{k+r-1}\left({\widehat{q}}_k\right)\right]Z_{i,r-1}^k}{\sum _{i=1}^R{\left(X_{i,r-1}^k\right)}^2}$$ E5

The values of r are selected based on the residual principle [33], which are Eqs. (3) and (4) for the direct and indirect temperature measurement methods, respectively. The detailed information about the SFS method can be found in previous publications [10, 34].

2.2 Transfer function method

The transfer functions establish the relationship between the input and output in a dynamic system, which can also be used to solve the linear heat conduction problems, where the heat flux is treated as the input of the system and the temperature is treated as the response [13].

The estimated surface heat flux using Laplace transform can be described as

$$q\left(t\right)=L^{-1}\left[\frac{1}{H_c\left(s\right)}\right]\ast {\theta }_c\left(t\right)$$ E6

where θ_c (t) is the temperature allowance (θ_c (t) = T_c (t) − T ₀), the subscript c denotes the measurement position, the superscript* is the convolution integral operator, and T ₀ is the initial temperature. L ⁻¹[1/H(s)] is the Laplace inverse transform of the transfer function, which can be written as the function of time t, as follows:

$$L^{-1}\left[\frac{1}{H\left(s\right)}\right]=f\left(t\right)$$ E7

Substituting Eq. (7) into Eq. (6), q(t) becomes

$$q\left(t\right)={\int }_0^t{\theta }_c\left(t-\tau \right)f\left(\tau \right)\mathrm{d}\tau$$ E8

After dispersing Eq. (8), it takes the following form:

$$\left[\begin{array}{c}{q}_1\\ q_2\\ ⋮\\ q_n\end{array}\right]=\left[\begin{array}{cccc}{\theta }_{c,1}& 0& \cdots & 0\\ {\theta }_{c,2}& {\theta }_{c,1}& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ {\theta }_{c,n}& {\theta }_{c,n-1}& \cdots & {\theta }_{c,1}\end{array}\right]\left[\begin{array}{c}{f}_0\\ f_1\\ ⋮\\ f_n\end{array}\right]$$ E9

An assist question can be established to solve the temperature allowance θ assuming the boundary condition q(t) = 1. Substituting the numerical solution into Eq. (8), the surface heat flux q(t) can be obtained.

2.3 Duhamel’s theorem method

Duhamel’s theorem directly treated the measured internal temperature as the surface temperature of the substrate when using indirect FTC temperature measurement. A one-dimensional, direct heat conduction problem for the two surface temperature measurements can be solved to calculate the temperature distribution in the substrate [14, 35]:

$$T\left(x,t\right)=f\left(0\right)\phi \left(x,t\right)+{\int }_0^t\phi \left(x,t-\tau \right)\frac{\mathit{dT}\left(\tau \right)}{\mathit{d\tau }}\mathit{d\tau }$$ E10

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

$$\phi \left(x,t\right)=1-\mathit{erf}\left(\frac{x}{2\sqrt{\mathit{\alpha t}}}\right)$$ E11

Solving Eq. (10), the temperature gradient at the surface can be obtained:

$$q\left(t\right)=\sqrt{\frac{\mathit{\lambda \rho c}}{\pi }}{\int }_0^t\frac{1}{\sqrt{t-\tau }}\frac{\mathit{dT}}{\mathit{d\tau }}\mathit{d\tau }$$ E12

The simplification that the internal temperature equals the surface ones ignoring the heat dissipating capacity of the materials above the thermocouple, thereby often causing significant heat flux errors. Therefore, Duhamel’s theorem needs to be improved to deal with this problem, rather than directly estimate surface heat flux from the measured temperature data. Firstly, the real surface temperature is needed to be calculated from the internal measured temperature at x = c location, which can be expressed by a piecewise constant function of time as

$$f=f\left(i\cdot \mathrm{\Delta }t\right),i=1,2,\dots ,n,$$ E13

where f is the surface temperature at the time of i·Δt. The temperature at x = c location can be solved by Duhamel’s theorem as [35].

$$\begin{array}{l}{T}_c\left(n\cdot \mathrm{\Delta }t\right)=f\left(0\right)\phi \left(x,n\cdot \mathrm{\Delta }t\right)+\left[,f\left(1\cdot \mathrm{\Delta }t\right)-f\left(0\right)\right]\phi \left(x,\left(n-1\right)\cdot \mathrm{\Delta }t\right)\\ + \left[,f\left(2\cdot \mathrm{\Delta }t\right)-f\left(1\cdot \mathrm{\Delta }t\right)\right]\phi \left(x,\left(n-2\right)\cdot \mathrm{\Delta }t\right)\cdots \\ + \left[,f\left(n\cdot \mathrm{\Delta }t\right)-f\left(\left(n-1\right)\cdot \mathrm{\Delta }t\right)\right]\phi \left(x,0\right)\end{array}$$ E14

Eq. (14) can be written in expanded matrix form as

$$\left[\begin{array}{cccc}\mathrm{\Delta }\phi \left(\mathrm{\Delta }t\right)& 0& \cdots & 0\\ \mathrm{\Delta }\phi \left(2\mathrm{\Delta }t\right)& \mathrm{\Delta }\phi \left(\mathrm{\Delta }t\right)& \cdots & 0\\ ⋮& ⋮& \ddots & 0\\ \mathrm{\Delta }\phi \left(n\mathrm{\Delta }t\right)& \mathrm{\Delta }\phi \left(\left(n-1\right)\mathrm{\Delta }t\right)& \cdots & \mathrm{\Delta }\phi \left(\mathrm{\Delta }t\right)\end{array}\right]\left[\begin{array}{c}f\left(0\cdot \mathrm{\Delta }t\right)\\ f\left(1\cdot \mathrm{\Delta }t\right)\\ ⋮\\ f\left(\left(n‐1\right)\cdot \mathrm{\Delta }t\right)\end{array}\right]=\left[\begin{array}{c}{T}_c\left(0\cdot \mathrm{\Delta }t\right)\\ T_c\left(1\cdot \mathrm{\Delta }t\right)\\ ⋮\\ T_c\left(n\cdot \mathrm{\Delta }t\right)\end{array}\right]$$ E15

where Δφ(iΔt) represents φ(x, i·Δt)−φ (x, (i−1)·Δt). The surface temperature can be obtained by multiplying the inverse temperature transformation matrix on both sides of the equation from the measured internal temperature data. Afterwards, the surface heat flux can be estimated using Eq. (12).

2.4 Validation of 1D filter solution

Two different methods were employed to measure the surface temperature during CSC. Figure 1 shows a schematic of the thin-film thermocouple (TFTC) and the fine thermocouple (FTC) measurements. The type-T TFTC with a thickness of 2 μm was directly deposited onto the epoxy resin surface using the magnetron technique. It has perfect contact with the underlying substrate and a fast response time (∼1 μs). The FTC measurement with a response time of 3.33 ms [36] consists of a fine type-T thermocouple (∼10 μm bead diameter) underneath a thin layer of aluminum foil (∼10 μm), which is positioned on the surface of the epoxy resin substrate (∼5 mm), as shown in Figure 1(b). The thermal paste with thickness of 100 μm is placed between the aluminum foil and the epoxy resin, ensuring good thermal contact and providing mechanical support.

Figure 1.

Schematic of the TFTC and FTC measurements with single temperature sensor [28]: (a) TFTC measurement and (b) FTC measurement.

A hypothetical triangular pulse heat flux is commonly used to examine the accuracy and sensitivity to random noises of the algorithms [10, 13]. This given heat flux is first set as the surface boundary condition, and then the temperature adding random noise (ε = 1°C) at the measured point corresponding to the two different measurement methods is obtained by solving the direct heat conduction problem. New surface heat flux can be predicted based on the calculated temperature through different algorithms.

Figure 2 shows the new estimated heat fluxes using different algorithms with random noise using TFTC and FTC measurements. The heat fluxes calculated by Duhamel’s theorem, SFS, and the transfer function method all agreed well with the hypothetical ones using the TFTC measurement. When using FTC measurement (Figure 2(b)), the estimated heat flux obtained from Duhamel’s theorem and SFS methods changed nonlinearly over time. Furthermore, they can hardly match the exact heat flux well. Thus, they were unsuitable for estimating the surface heat flux in the indirect surface temperature measurement with FTC. In comparison, both the transfer function and improved Duhamel’s theorem can provide linear and accurate results that match the given heat flux exactly.

Figure 2.

Results of the comparison between estimated heat fluxes (scatters) and hypothetical ones (solid lines) [28]. (a) single-layer IHCP with single TFTC and (b) three-layer IHCP with single FTC.

2.5 1D heat flux estimation

Figure 3 shows the variation in surface temperature as a function of time, measured by the direct measurement method with the TFTC and the indirect measurement method with FTC. The temperature first decreased rapidly when the R404A droplets impinge on the substrate surface, following a relatively slow change as the thin liquid film forms on the surface. It began to resume to the ambient temperature after the liquid film completely evaporated. It is notable that the surface temperature measured by the TFTC measurement decreased faster than that by the indirect FTC measurement method, which showed a definite delay.

Figure 3.

Variations of surface temperature measured by TFTC and FTC measurements [28].

Figure 4(a) depicts the time-varied surface heat flux predicted by Duhamel’s theorem, SFS method, and transfer function method with TFTC measurement. It can be seen that all the three algorithms predicted very similar results. They first increased rapidly to their maximum values, then dropped quickly to a certain value (∼350 kW/m²), and finally gradually decreased to zero. Note that the estimated result by the SFS method increased faster and had a slightly greater maximum heat flux than that obtained by the other two algorithms.

Figure 4.

Variations of surface heat flux predicted by different algorithms [28]. (a) TFTC measurement and (b) FTC measurement.

Figure 4(b) represents the heat flux predicted by different algorithms under the FTC measurement. As mentioned above, Duhamel’s theorem was inappropriate for predicting surface heat flux directly from internal measured temperature. Thus, the q _max and the magnitude of its increase rate by Duhamel’s theorem were both far lower than the case with other algorithms. It seems that the heat flux was significantly delayed and reduced. Other surface heat fluxes predicted by the transfer function and newly improved Duhamel’s theory had almost the same varying history, with a larger q _max than that predicted by the SFS method.

3. 2D algorithm for heat flux estimation

Although the TF and improved Duhamel’s theorem can accurately estimate the surface heat flux with multilayer mediums, it still failed to solve 2D IHCPs. Moreover, large nonuniformity radial and temporal surface temperature variations during CSC often occurred [14, 32]. Therefore, lateral heat transfer must be considered. The 2D algorithm for 2D multilayer medium is urgent to be developed. In this chapter, the 2D filter solution was proposed and validated based on the three-layer FTC temperature measurement.

3.1 2D filter solution and optimal comparison criterion

The 2D general model with K heat fluxes, J temperature sensors, and I layers is shown in Figure 5. The multiple surface heat fluxes are assumed to be uniform in a specific upper lateral surface area (e.g., 0 ≤ x ≤ x ₁). The temperature sensors are placed at the lateral midpoint of each uniform surface heat flux range in the ith layer, and the interval between two thermocouples is ∆x. All sensors are placed at an identical depth (y = H _T) from the upper surface. The adiabatic boundary condition is considered on the three other sides (x = 0, x = W, and y = H). The initial temperature of the substrates is T ₀.

Figure 5.

Geometry model of a 2D multilayer IHCP [36].

The discrete solution of temperature for direct heat conduction problems with K heat fluxes and J sensors can be presented in a matrix form as follows [37].

$$\mathrm{T}=\mathbf{X}\mathrm{q}$$ E16

where T is the real temperature matrix, X is the sensitivity matrix, and q denotes multiple heat fluxes. The sensitivity matrix X is [19]

$$\mathbf{X}=\left[\begin{array}{ccccc}\mathrm{a}\left(1\right)& 0& 0& \cdots & 0\\ \mathrm{a}\left(2\right)& \mathrm{a}\left(1\right)& 0& \cdots & 0\\ \mathrm{a}\left(3\right)& \mathrm{a}\left(2\right)& \mathrm{a}\left(1\right)& \cdots & ⋮\\ ⋮& ⋮& ⋮& \ddots & 0\\ \mathrm{a}\left(\mathrm{N}\right)& \mathrm{a}\left(\mathrm{N}-1\right)& \mathrm{a}\left(\mathrm{N}-2\right)& \cdots & \mathrm{a}\left(1\right)\end{array}\right]$$ E17

where

$$\mathrm{a}\left(n\right)=\left[\begin{array}{ccccc}{a}_{11}\left(n\right)& a_{12}\left(n\right)& a_{13}\left(n\right)& \cdots & a_{1\mathrm{K}}\left(n\right)\\ a_{21}\left(n\right)& a_{22}\left(n\right)& a_{23}\left(n\right)& \cdots & a_{2\mathrm{K}}\left(n\right)\\ a_{31}\left(n\right)& a_{32}\left(n\right)& a_{33}\left(n\right)& \cdots & a_{3\mathrm{K}}\left(n\right)\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ a_{\mathrm{J}1}\left(n\right)& a_{\mathrm{J}2}\left(n\right)& a_{\mathrm{J}3}\left(n\right)& \cdots & a_{\mathrm{JK}}\left(n\right)\end{array}\right]$$ E18

$$a_{\mathit{jk}}\left(n\right)=\frac{\partial \left(T\left(x_{{\mathrm{TC}}_j},{\mathrm{W}}_{\mathrm{T}},t\right)-T_0\right)}{\partial \left[q_k\left(1\right)\right]}={\phi }_{\mathit{jk}}\left(n\right)$$ E19

where N and n are total time steps and current time step. Excessive temperature ϕ can be presented as follows:

$${\varphi }_{\mathit{jk}}\left(n\right)=T\left(x_{{\mathrm{TC}}_j},H_{\mathrm{T}},{\mathit{nt}}_{\mathrm{d}}\right)\left|{}_{q_k=1,q_1=q_2=\cdots \cdots =q_{\mathrm{K}}=0}\right.-T_0$$ E20

where t _d is the time step. The sum of the squares of the errors between the estimated and measured temperatures plus one-order regularization, S, is

$$\mathrm{S}={\left(\mathrm{Y}‐\mathrm{T}\right)}^{\text{'}}\left(\mathrm{Y}‐\mathrm{T}\right)+{\mathrm{\alpha }}_{\mathrm{t}}{\left[{\mathrm{H}}_{\mathrm{t}}\mathrm{q}\right]}^{\text{'}}\left[{\mathrm{H}}_{\mathrm{t}}\mathrm{q}\right]+{\mathrm{\alpha }}_{\mathrm{s}}{\left[{\mathrm{H}}_{\mathrm{s}}\mathrm{q}\right]}^{\text{'}}\left[{\mathrm{H}}_{\mathrm{s}}\mathrm{q}\right]$$ E21

where Y is the measured temperature matrix. α_t and α_s are regularization parameters with respect to temporal and spatial terms. The superscript ‘denotes the transpose of a matrix. H_t and H_s are temporal and spatial regularization matrixes. Minimizing S, estimated heat flux matrix $\widehat{\mathbf{q}}$ using the entire domain data can be obtained [8]:

$$\widehat{\mathrm{q}}={\left[{\mathrm{X}}^{\text{'}}\mathrm{X}+{\mathrm{\alpha }}_{\mathrm{t}}{{\mathrm{H}}_{\mathrm{t}}}^{\text{'}}{\mathrm{H}}_{\mathrm{t}}+{\mathrm{\alpha }}_{\mathrm{s}}{{\mathrm{H}}_{\mathrm{s}}}^{\text{'}}{\mathrm{H}}_s\right]}^{-1}{\mathrm{X}}^{\text{'}}\mathrm{Y}=\mathrm{FY}$$ E22

Generally, the 2D multilayer IHCPs are solved layer by layer, starting from the ith layer with the known experimental measured temperature. The heat flux is estimated at the interface with the (i−1)th layer. Thereafter, the interface heat flux between the (i−1)th and (i−2)th layers can be calculated by using the known interface heat flux as the input. Finally, the surface heat flux is estimated [18]. The interface heat fluxes ${\widehat{\mathbf{q}}}_i$ between the ith and (i−1)th layers can be described as

$${\widehat{\mathrm{q}}}_{\mathrm{i}}={\left[{{\mathrm{X}}_{\mathrm{i}}}^{\text{'}}{\mathrm{X}}_{\mathrm{i}}+{\mathrm{\alpha }}_{\mathrm{t}}{{\mathrm{H}}_{\mathrm{t}}}^{\text{'}}{\mathrm{H}}_{\mathrm{t}}+{\mathrm{\alpha }}_{\mathrm{s}}{{\mathrm{H}}_{\mathrm{s}}}^{\text{'}}{\mathrm{H}}_{\mathrm{s}}\right]}^{‐1}{{\mathrm{X}}_{\mathrm{i}}}^{\text{'}}\mathrm{Y}={\mathrm{F}}_{\mathrm{i}}\mathrm{Y}$$ E23

Notably, X _i is different from X in the solution of single-layer IHCP. X _i is calculated by using (I − i + 1) layers below the ith layer. The estimated interface temperature ${\widehat{\mathbf{T}}}_i$ yields

$${\stackrel{⌢}{\mathrm{T}}}_i={\mathbf{X}}_i{\widehat{\mathrm{q}}}_i$$ E24

The interface heat fluxes ( ${\widehat{\mathbf{q}}}_i$ ) and temperature ( ${\widehat{\mathbf{T}}}_i$ ) are used as the known boundary condition to solve the solution for the (i−1)th layer, where the heat fluxes are unknown at y = H ₁ + H ₂ + ··· + H _i−2. The interface heat fluxes ${\widehat{\mathbf{q}}}_{i-1}$ between the (i−1)th and (i−2)th layers can be expressed as

$${\widehat{\mathrm{q}}}_{i-1}={\left[{{\mathbf{X}}_{i,\mathrm{d}}}^{\text{'}}{\mathbf{X}}_{i,\mathrm{d}}+{\alpha }_t{{\mathbf{H}}_{\mathrm{t}}}^{\text{'}}{\mathbf{H}}_{\mathbf{t}}+{\alpha }_s{{\mathbf{H}}_{\mathbf{s}}}^{\text{'}}{\mathbf{H}}_{\mathbf{s}}\right]}^{-1}{{\mathbf{X}}_{\mathrm{d}}}^{\text{'}}{\stackrel{⌢}{\mathrm{T}}}_i-{\mathbf{X}}_{i,\mathrm{u}}{\widehat{\mathrm{q}}}_i={\mathbf{F}}_{i,\mathrm{d}}{\stackrel{⌢}{\mathrm{T}}}_i-{\mathbf{X}}_{i,\mathrm{u}}{\widehat{\mathrm{q}}}_i$$ E25

where X _i,d and X _i,u have a similar form as the sensitivity matrix and X is for single-layer IHCP. However, excessive temperature ϕ in X is different. For X _i,d, ϕ yields

$${\varphi }_{\mathit{jk}}\left(n\right)=T\left(x_{{\mathrm{TC}}_j},H_1+H_2+\cdots +H_{i-1},{\mathit{nt}}_{\mathrm{d}}\right)\left|{}_{q_k=1,q_1=q_2=\cdots \cdots =q_{\mathrm{K}}=0}\right.-T_0$$ E26

For X _i,u, ϕ can be described as

$${\varphi }_{\mathit{jk}}\left(n\right)=T\left(x_{{\mathrm{TC}}_j},H_1+H_2+\cdots +H_{i-1},{\mathit{nt}}_{\mathrm{d}}\right)\left|{}_{q_k=1,q_1=q_2=\cdots \cdots =q_{\mathrm{K}}=0}\right.-T_0$$ E27

By substituting Eq. (23) into Eq. (25), the estimated interface heat flux is derived as

$${\widehat{\mathrm{q}}}_{i-1}={\mathbf{F}}_{i,\mathrm{d}}\left({\mathbf{X}}_i-{\mathbf{X}}_{i,\mathrm{u}}\right){\mathbf{F}}_i\mathbf{Y}={\mathbf{F}}_{i-1}\mathbf{Y}$$ E28

Interface heat flux ${\widehat{\mathbf{q}}}_{i-2}$ , ${\widehat{\mathbf{q}}}_{i-3}$ ,··· ${\widehat{\mathbf{q}}}_2$ can be calculated by repeating Eqs. (25)-(28). Thereafter, surface heat flux matrix $\widehat{\mathbf{q}}$ can be estimated. The solutions for 2D single-layer and multilayer IHCPs are then obtained.

Given that most filter coefficients can be disregarded except for those of the (m_p + m_f + 1) time step, the solutions for IHCPs (Eqs. (22) and (28)) can be simplified into a general filter solution, thereby allowing the real-time heat flux monitoring and small computational load [8, 19, 38]:

$$q_k\left(n\right)=\sum _{m=1}^{\left({\mathrm{m}}_{\mathrm{p}}+{\mathrm{m}}_{\mathrm{f}}+1\right)\times \mathrm{J}}\left(f_{k,m}{Y}_{\left(n-{\mathrm{m}}_{\mathrm{p}}-1\right)\times \mathrm{J}+m}\right)$$ E29

where f _k,m denotes the filter coefficient in the mth column from one row of F associated with the kth unknown heat flux. Notably, α_t and α_s significantly affect the accuracy of the estimated heat flux, which should be determined for a specific IHCP [9, 18, 19]. These parameters can be optimized by the optimal comparison criterion developed for solving 2D single-layer and multilayer IHCPs [36]. The sum of the squares of the errors, R_q, between the estimated and real heat fluxes was used to examine the accuracy of the calculation:

$${\mathrm{R}}_{\mathrm{q}}\mathbf{=}{\left(\widehat{\mathrm{q}}-\mathbf{q}\right)}^{\text{'}}\left(\widehat{\mathrm{q}}-\mathbf{q}\right)\mathbf{=}{\left(\mathbf{FY}-\mathbf{q}\right)}^{\text{'}}\left(\mathbf{FY}-\mathbf{q}\right)$$ E30

A more efficient means is to minimize the expected value of R_q [19, 38] as follows:

$$E\left({\mathrm{R}}_{\mathrm{q}}\right)\mathbf{=}{\left[\left(\mathrm{FX}-\mathrm{I}\right)\mathrm{q}\right]}^{\text{'}}\left[\left(\mathrm{FX}-\mathrm{I}\right)\mathrm{q}\right]+{\sigma }_Y^2\mathit{tr}\left({\mathbf{F}}^{\text{'}}\mathbf{F}\right)$$ E31

where E is the expected value, tr denotes the sum of the diagonal elements, and σ_Y (0.01°C [38]) is the uniform measurement errors. E(R_q) contains two components: E _q,bias :

$$E_{q,\mathit{\text{bias}}}={\left[\left(\mathrm{FX}-\mathrm{I}\right)\mathrm{q}\right]}^{\text{'}}\left[\left(\mathrm{FX}-\mathrm{I}\right)\mathrm{q}\right]$$ E32

and E _q,rand in the filter form

$$E_{q,\mathit{rand}}={\sigma }_Y^2\mathit{tr}\left({\mathbf{F}}^{\mathrm{T}}\mathbf{F}\right)={\sigma }_Y^2\left(\mathrm{N}-\left({\mathrm{m}}_{\mathrm{p}}+{\mathrm{m}}_{\mathrm{f}}+1\right)\right)\sum _{k=1}^{\mathrm{K}}\sum _{j=1}^{\left({\mathrm{m}}_{\mathrm{p}}+{\mathrm{m}}_{\mathrm{f}}+1\right)}{f}_{k,j‐1\ast \mathrm{J}+k}^2$$ E33

To examine the accuracy of the estimated heat fluxes, the mean relative error (MRE) between estimated and hypothetical heat fluxes was employed as follows:

$$q_{\mathit{MRE}}=\frac{1}{\mathrm{N}-{\mathrm{m}}_{\mathrm{f}}}\sum _{n=1}^{\mathrm{N}-{\mathrm{m}}_{\mathrm{f}}}\left|\frac{\widehat{q}\left(n\right)-q\left(n\right)}{q\left(n\right)}\right|$$ E34

3.2 Validation of 2D filter solution

The 2D three-layer geometry (Figure 6) with six FTC sensors was used to validate the accuracy of surface heat flux estimated by 2D filter solution. For FTC measurement, the multiple surface heat fluxes were assumed to be uniform in a specific upper lateral surface area (e.g., 0 ≤ x ≤ x ₁) and different from each other. Temperature was measured from the spray center (x = 0 mm, TC₁) to the periphery (x = 10 mm, TC₆), and the lateral distance between two FTCs was ∆x = 2 mm. The geometry width with TFTC measurement was W = 11 mm, which is half that of the spray diameter because the computational domain is symmetric.

Figure 6.

Geometry model of three-layer IHCP with six FTCs [36].

The optimal comparison criterion developed for 2D IHCPs was employed in this study to optimize the Tikhonov regularization parameters (α_t and α_s ). As shown in Figure 7, logarithmic random component E _q,rand increased as the regularization parameters α_t and α_s decreased. E _q,bias and E(R_q) reached their minimum value at α_t  = α_s  = 10⁻⁹. Eventually, the optimum regularization parameters were determined to be α_t  = α_s  = 10⁻⁹.

Figure 7.

Geometry model of three-layer IHCP with six FTCs [36].

As shown in Figure 8, the estimated heat fluxes agreed well with the hypothetical ones. However, small deviations were observed at the descent stage after the maximum heat flux for q ₁ to q ₆. The q_MRE of q ₁–q ₆ for TFTC measurement was 2.63%, 2.66%, 2.76%, 2.78%, 3.00%, and 5.18%, respectively. Adding the random noise at the measured point, the maximum q_MRE increased to 3.71%, which indicated that the accuracy and stability of the filter solution are satisfactory.

Figure 8.

Results of the comparison between estimated heat fluxes (scatters) and hypothetical ones (solid lines) for three-layer IHCP with FTCs [36].

3.3 Comparison between 2D filter solution and 1D improved Duhamel’s theorem

To investigate the importance considering lateral heat transfer, the estimated surface heat flux (Figure 9(a)) and simulated temperature (Figure 9(b)) using the estimated surface heat flux as boundary, calculated by 2D filter solution and 1D improved Duhamel’s theorem, were compared taking the example of 2D three-layer geometry with FTC measurement. For FTC measurement containing aluminum film with high heat conductivity coefficient (λ = 236 W∙m⁻¹∙K⁻¹), the maximum heat flux calculated by the 1D method is underestimated by 60% than that calculated by 2D filter solution, indicating that the lateral heat transfer cannot be disregarded, especially when the heat conductivity coefficient of the material is large. As shown in Figure 9(b), the simulated temperature using estimated surface heat flux as boundary calculated by 2D filter solution agreed well with measured temperature, indirectly indicating the accuracy of 2D filter solution. However, the large simulated temperature deviation was observed using 1D improved Duhamel’s theorem, owing to the inaccurate estimated surface heat flux disregarding lateral heat transfer.

Figure 9.

Comparative results at spray center calculated by 2D filter solution and 1D improved Duhamel’s theorem [36]. (a) surface heat flux and (b) temperature at sensor location.

The dynamic internal temperature measured by six FTCs with ∆t = 50 ms and L = 30 mm is depicted in Figure 10(a). The temperature histories were similar, but differences existed in the specific values. The minimum temperatures (T _min) at r = 0, 2, 4, 6, 8, and 10 mm were −43.42, −36.13, −29.55, −28.22, −27.32, and −22.45°C, respectively. Additionally, the measured temperature was lower at the spray center (r = 0 mm) than periphery. Figure 10(b) presents the estimated heat fluxes calculated by the filter solution for 2D three-layer IHCP with FTC measurement. The estimated heat flux profiles at different locations in this figure were also similar. However, a large difference was observed in heat fluxes at different lateral locations. The best cooling capacity was found at the spray center (r = 0–2 mm).

Figure 10.

Variations of surface temperature and heat flux with FTC measurement [36]. (a) measured surface temperature and (b) estimated surface heat flux.

4. Clinical potential of spray cooling by low GWP R1234yf

The extremely high global warming potential (GWP = 1430) of commercially used cryogen R134a with boiling point of −26.1°C will cause severe environmental hazards, and the substitution of R134a in clinical application is urgent. R1234yf with boiling point of −29.5°C may be a potential candidate for environment protection due to its low GWP (<1). In this chapter, the clinical potential of R1234yf substitution for R134a was investigated.

Using the maximum surface heat flux correlation obtained by experimental spray characteristics (droplet temperature, velocity, and diameter) and surface heat transfer performance (surface heat flux calculated by 2D filter solution), Figure 11(a) shows the variations of effective surface heat flux with different R134a and R1234yf as a function of spray distance. The effective surface heat flux (q_e ) was obtained by multiplying the maximum surface heat flux, and cooling concentration within the radius of 2 mm is the interested area [39]. The effective surface heat flux by using R134a and R1234yf increased firstly due to the droplet temperature reduction as spray distance increased. Then, q_e reaches their maximum value and finally decreases slowly. The maximum q_e of R134a and R1234yf is 262.1 and 225.8 kW/m² at the optimal spray distances of 25.6 and 25.1 mm (see Figure 11(a)), respectively. The substitution of R1234yf for R134a can produce remarkable reduction of global warming potential to <1. However, the cooling capacity should be enhanced for the clinical application in laser treatment, owing to the 13.8% reduction in effective heat flux (from 262.1 to 225.8 kW/m²). Therefore, the enhancement of cooling capacity is necessary for the implementation of R1234yf in clinical laser treatment of PWS.

Figure 11.

Variations of effective surface heat flux as a function of spray distance [39]. (a) effective surface heat flux with R134a and R1234yf and (b) effective surface heat flux before and after cooling enhancement with R1234yf.

According to our experience, two simple ways are available to enhance the cooling capacity, i.e., changing the nozzle diameter and decreasing the back pressure by decreasing the boiling point of cryogens [40, 41]. As shown in Figure 11(b), the enhancement of effective surface heat flux at different spray distances was remarkable. After reducing the nozzle diameter and decreasing the back pressure, q_e increases rapidly due to violent evaporation, and the peak value of the effective heat flux (268.3 kW/m²) is increased by 18.8%. This result is comparable with that of R134a under 1 atm, which proves the potential of R1234yf in the clinic CSC for the laser treatment of PWS.

5. Conclusions

Several algorithms including the SFS, TF, and Duhamel’s theorem methods were analyzed and compared in predicting time-varying surface heat flux during CSC. Duhamel’s theorem was improved to get the accurate results through the transformation of internal temperature into surface temperature, when the indirect surface temperature measurement (FTC) method is used. A hypothetical triangular pulse heat flux was employed to analyze the accuracy and sensitivity to noise of the algorithms under TFTC and FTC measurements. The estimated result of Duhamel’s theorem and SFS method widely deviated from the given heat flux under the three-layer FTC measurement method, whereas the transfer function and improved Duhamel’s theorem all provided the exact estimated heat flux.

The 2D filter solution was proposed to solve a general 2D multilayer IHCP for the estimation of surface heat flux. An optimal comparison criterion was employed to optimize the key parameters, namely, α_t and α_s . Six hypothetical triangular heat fluxes and random temperature errors of 1°C were employed to analyze the accuracy and sensitivity of the filter solution for 2D three-layer IHCPs with FTC measurement. The q_MRE values for FTC measurement with and without the random temperature errors were all within the acceptable range, which validates the good accuracy and stability of the filter solutions. The maximum heat flux calculated by the 1D method was underestimated by 60% than that calculated by 2D filter solution considering lateral heat transfer. The 2D filter solution was more accurate than the 1D method. Moreover, lateral heat transfer should not be ignored, especially when the heat conductivity coefficient of the material is large.

The surface heat transfer characteristics of spray cooling with R134a and R1234yf were investigated based on 2D filter solution. The maximum effective heat fluxes q _e,max were 262.1 and 225.8 kW/m² for R134a and R1234yf at different spray distances of 25.6 and 25.1 mm. Through the cooling enhancement of reducing the nozzle diameter and decreasing the back pressure, q _e,max of R1234yf was increased by 18.8% (D = 0.4 mm and P_b  = 0.04 MPa). The enhanced q _e,max is a bit higher than that of R134a in normal condition, which provides a theoretical basis for potential application of low GWP R1234yf in clinics.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (51727811) and Fundamental Research Funds for the Central Universities.

Conflict of interest

The authors declare no conflict of interest.