A robust and fast algorithm for three-dimensional transient inverse heat conduction problems

International Journal of Heat and Mass Transfer 55 (2012) 7865–7872

Contents lists available at SciVerse ScienceDirect

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

A robust and fast algorithm for three-dimensional transient inverse heat conduction problems Shuai Lu c, Yi Heng a,b,⇑, Adel Mhamdi b a

AixCAPE e.V., 52062 Aachen, Germany AVT-Process Systems Engieering, RWTH Aachen University, 52056 Aachen, Germany c School of Mathematical Science, Fudan University, 200433 Shanghai, People’s Republic of China b

a r t i c l e

i n f o

Article history: Received 9 February 2012 Received in revised form 27 July 2012 Accepted 7 August 2012 Available online 5 September 2012 Keywords: 3D inverse heat conduction problems Conjugate gradient method Tikhonov regularization Model function approach

a b s t r a c t Despite numerous studies of inverse heat conduction problems (IHCP) over the last several decades, their solutions still suffer from the mathematical difficulties and the bottleneck of currently available numerical methods for large-scale problems. In this paper, we present a robust and efficient algorithm for the solution of a specific type of three-dimensional (3D) IHCP commonly involved in various engineering applications. The solution method incorporates the Tikhonov regularization for tackling the severe ill-posedness and the conjugate gradient (CG) method for solving the resulting minimization problems. A model function approach is used to significantly reduce the effort needed to find the optimal Tikhonov regularization parameter. The proposed solution method requires no a priori knowledge of the measurement noise and is much more computationally efficient than the traditional Tikhonov regularizationbased inversion approaches. Thus, it can be used for the efficient solution of large-scale practical problems. Two simulation case studies of practical significance are presented to validate and assess the performance of the proposed method. Finally, the solution method is successfully applied to the reconstruction of instantaneous heat fluxes from experimentally measured temperature data. Ó 2012 Elsevier Ltd. All rights reserved.

1. Introduction Over the last several decades, the study of inverse heat transfer problems (IHTP) has gained great research interests owing to the computational challenges and practical significance. Different types of IHTP such as boundary IHTP, retrospective IHTP, coefficient IHTP and geometric IHTP, etc., have been extensively studied. In this paper, we study a specific type of IHTP that belongs to the class of inverse heat conduction problems (IHCP), namely the reconstruction of unknown heat flux on an inaccessible boundary from temperature observations on an accessible boundary. This kind of problem has various applications in different fields of engineering. Since the pioneer papers [1–3], numerous numerical and application-oriented IHCP papers have been published. They are mostly based on solution strategies such as the Tikhonov regularization [4], space marching [5], function specification [6], iterative regularization methods [7–9] or Bayesian inference approach based on the maximum likelihood estimator [10]. Other type of approaches such as a stochastic algorithm known as the quantum-behaved particle swarm optimization [11] or a decentralized fuzzy inference method [12] can also be applied to solve smallscale IHCP. ⇑ Corresponding author. E-mail address: [email protected] (Y. Heng). 0017-9310/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ijheatmasstransfer.2012.08.018

The numerical solution of IHCP is usually nontrivial due to both the severe ill-posedness and the high computational cost when a three-dimensional (3D) computational domain is employed. Most previous IHCP papers are restricted to one or two space dimensions. However, the reduction in geometric complexity usually sacrifices solution accuracy. A few recent works on 3D IHCP are briefly reviewed below. A procedure involving dynamic state observers was reported to solve multi-dimensional IHCP [13]. In [14] specialized methods for the solution of anisotropic 3D IHCP were studied, which have been applied to investigate the influence of the heat transfer in falling films by the wave characteristics. An iterative regularization strategy based on the conjugate gradient method for the normal equation for the efficient solution of a 3D IHCP in pool boiling has been developed in [15]. An analytical solution based on Green’s function was presented in [16] to solve a 3D IHCP. Among currently available approaches, Tikhonov regularization and iterative regularization are commonly used to solve IHCP and iterative regularization is considered very efficient [9]. In [17] an iterative regularization approach based on gradient methods to reconstruct space–time heat-flux distribution has been presented. A conjugate gradient (CG) type method was developed in [18] for the iterative solution of a transient boundary IHCP. Such an approach has been successfully applied to various practical problems, e.g. the estimation of transient heat flux at the inaccessible film side in a falling film experiment [19], the reconstruction

7866

S. Lu et al. / International Journal of Heat and Mass Transfer 55 (2012) 7865–7872

Nomenclature Y z

Latin symbols a temporal function b spatial function cp specific heat capacity C constant parameter dk data error in iteration k f frequency of laser heat flux F function operator of forward heat conduction J objective functional k iteration number m model function n outer normal on the object boundaries pk conjugate search direction in iteration k q unknown heat flux on the boundary C3 q0 initial guess for the unknown heat flux q1 heat flux on the boundary C1 q2 heat flux on the boundary C2 qa minimizer of the Tikhonov cost functional ^ q heat-flux estimate r radius of Gaussian laser beam Sk solution of sensitivity problem in iteration k t time variable t0 backwards time tf final time of observation T temperature variable; constant parameter x horizontal coordinate on the surface xc constant parameter x 2D or 3D spatial coordinate X heat-flux solution space y vertical coordinate on the surface yc constant parameter

Greek symbols Tikhonov regularization parameter optimal Tikhonov regularization parameter b surface absorptivity l constant parameter lk search step length in iteration k ck conjugate coefficient in iteration k C1 surface boundary C2 surface boundary C3 surface boundary d upper bound of measurement error  error threshold H temperature variable H0 initial temperature distribution Hm noise-free temperature measurements Hdm corrupted temperature measurements k thermal conductivity q material density wk solution of adjoint problem in iteration k Wl, WC,T,l objective functionals x zero mean normal distribution with variance one; contraction factor X 3D computational domain

a a⁄

Mathematical notations L2 square integrable functions r gradient operator kk function norm

of the instantaneous local boiling heat flux in a pool boiling experiment [20] and the recovery of high intensity periodic laser heat flux on the front surface of a target object [21]. To simplify the following discussion and the motivation for this work, we introduce briefly the considered transient heat conduction equation representing the forward problem:

@H qcp ¼ r  ðkrHÞ; in X  ð0; t f Þ; @t Hð; 0Þ ¼ H0 ðÞ; on X; @H ¼ q1 ; on C1  ð0; t f Þ; k @n @H k ¼ q2 ; on C2  ð0; t f Þ; @n @H k ¼ q; on C3  ð0; tf Þ; @n

data space coordinate directing from one surface to another

on the boundary C1. The error bound d of the measured temperatures is usually unknown or at least cannot be estimated accurately in practice. Let X, Y denote appropriate function spaces for the heat flux q and the temperature data Hdm , respectively. The 3D IHCP in the operator form then reads as

Find q 2 X such that FðqÞ ¼ Hdm ; ð1Þ ð2Þ ð3Þ ð4Þ ð5Þ

where H denotes the temperature defined on an observation time interval [0, tf] and on a 3D spatial domain X. q1, q2 and q denote the time- and space-dependent heat fluxes through different conductor boundaries C1, C2 and C3, respectively. q, cp and k denote the density, the heat capacity and the thermal conductivity, respectively. It is well-known that Eqs. (1)–(5) represents a well-posed parabolic partial differential equation (PDE) in the sense of Hadamard [22]. Issues of existence, uniqueness and stability of the solution of parabolic PDE have been well addressed in literature, e.g. in [23,24]. The 3D IHCP studied in this paper arises from various industrial applications. It corresponds to the reconstruction of the unknown heat flux q on the boundary C3 using temperatures Hdm observed

ð6Þ

where the function operator F : X ? Y is implicitly given by Eqs. (1)– (5). For the solution of the IHCP (6), either the CG-based iterative regularization method or the Tikhonov’s regularization method can be applied. Below, we first point out their major drawbacks and then present our efficient solution strategy to cope with these issues. The CG-based iterative regularization approach aims to solve a minimization problem in which the least squares error between the temperature measurements and model predictions

 2 JðqÞ :¼ FðqÞ  Hdm 

ð7Þ

is considered as the objective functional. k  k denotes for instance the norm in the set of L2 functions. In order to cope with the illposedness, an additional stopping rule, e.g. the discrepancy principle [25], has to be used to terminate the iterative process properly. Namely, the minimization procedure stops if the temperature prediction error matches roughly the noise level, which is assumed to be a priori known. Other stopping rules such as the L-curve criterion [26] can be applied if the measurement noise is unknown. However, it requires to perform many more iterations (cf. [20],

7867

S. Lu et al. / International Journal of Heat and Mass Transfer 55 (2012) 7865–7872

more than 1000 iterations were applied) in order to construct the entire L-curve for further analysis. Further, there is only little theoretical justification of the technique though it works well for a variety of problems. In fact, the CG-based iterative regularization approach using L-curve rule as a stopping criterion might be mathematically not stable. The Tikhonov regularization [4] targets at minimizing an objective functional which involves not only the term of data fitting error but also a term of solution norm, i.e.

 2 JðqÞ :¼ FðqÞ  Hdm  þ akqk2 ;

a > 0:

ð8Þ

Although Tikhonov’s regularization method is mathematically rigorous, its performance essentially depends on the choice of the regularization parameter a. In order to find the optimal regularization parameter, it usually requires solving many minimization problems for a large set of a, which leads to a significant increase of computational effort. In this paper, we present a robust and efficient solution method for the 3D IHCP (6). It combines the advantages of the iterative-type and Tikhonov-type of regularization methods. The solution method aims to minimize the objective functional (8) and uses the model function approach [27] to reduce the number of a tests. Further, an efficient CG method is proposed to solve each resulting minimization problem for a fixed a value. In such a way, the total computational effort can be significantly reduced. The proposed method is much faster than traditional Tikhonov regularization-based solution approaches. Compared to the previous CG-based iterative regularization approach, the method proposed in this paper is more robust and requires no a priori knowledge of measurement noise, which is useful for a wide range of practical problems. The paper is organized as follows. In Section 2, the enhanced solution procedure via Tikhonov regularization is briefly reviewed. Details of the efficient CG iteration used to minimize the objective functional (8) for a fixed a are given in Section 3. In Section 4, we validate the performance of the proposed solution method with two simulation examples of practical significance. The reconstruction of local heat fluxes on the inaccessible boiling surface with experimentally measured temperature data is given in Section 5. Conclusions are given in the final section. 2. Enhanced Tikhonov regularization for the continuous IHCP The solution of the 3D IHCP (6) by means of the Tikhonov regularization is obtained by minimizing the objective functional (8). From the definition, the Tikhonov cost functional (8) considers  2 the trade-off between the data fitting error FðqÞ  Hdm  and 2 the variation of the heat-flux estimate kqk . The a near the penalty term kqk2 is the Tikhonov regularization parameter used to control the reconstruction quality in case that temperature measurements are corrupted with noise. It can be shown that, for an appropriately fixed a, the minimizer of the Tikhonov cost functional (8) is stable with respect to perturbations of the data Hdm . Since the error bound d is generally unknown in practice, the traditional approach to choose such an a performs a large number of tests aiming at finding the optimal regularization parameter among them. Such a procedure is computationally very expensive. In this paper, we consider the use of an efficient method proposed in [27] to significantly speed up the solution procedure via Tikhonov regularization. Its basic idea is to update a by minimizing an alternative objective functional [28], namely the modified L-curve method,



2

Wl ðaÞ ¼ Fðqa Þ  Hdm  kqa k2l ;

l > 0:

ð9Þ

In [28] the author has proven that if the regularization parameter a⁄ chosen by the standard L-curve method, which has the   2 maximum curvature along the whole L-curve log Fðqa Þ  Hdm  ; log kqa k2 Þ, has a slope  l1 , then the minimizer of the functional (9) coincides with the a⁄. Nevertheless, since there is no explicit way computing the minimizer in [28], one still needs to consider a large a set and plot the functional Wl(a) for finding the minimum point. The fast algorithm proposed in [27] aims at replacing the data fitting error and the solution norm with appropriate model functions which essentially reduces the computational effort for finding the minimizer of (9). For completeness, we briefly review the idea of the model function approach. By taking q :¼ q(a) and considering the Euler equation of (8), we reformulate the original Tikhonov cost functional (8) in the following form

 2 JðaÞ :¼ FðqðaÞÞ  Hdm  þ akqðaÞk2  2 ¼ Hdm   kFðqðaÞÞk2  akqðaÞk2 :

ð10Þ

Note that the first derivative of J(a) with respect to a is J0 (a) = kq(a)k2, (i.e. [29]). We approximate the term kF(q(a))k2 in (10) locally by Tkq(a)k2 where the constant T is iteratively determined in the algorithm. The Eq. (10) then yields

 2 JðaÞ þ ða þ TÞJ 0 ðaÞ  Hdm  : The model function m(a) which mimics the performance of J(a) should satisfy the differential equation

 2 mðaÞ þ ða þ TÞm0 ðaÞ ¼ Hdm  ; that is m(a) can be represented in the form of a simple family of functions

 2 mðaÞ ¼ Hdm  þ

C

aþT

ð11Þ

;

where the constants C and T should be determined locally. We use the first derivative of the model function m(a) to approximate the first derivative of J(a), i.e J0 (a)  m0 (a). It is easy to observe that the original function Wl(a) can be (locally) approximated by another simple function WC,T,l(a) defined as

WC;T;l ðaÞ :¼ ðmðaÞ  am0 ðaÞÞðm0 ðaÞÞl ¼

 d 2 H  þ m

C aC þ a þ T ða þ TÞ2

!

C

!l ð12Þ

ða þ TÞ2

by introducing the model function m(a) whose minimizer can be explicitly calculated. A rigorous mathematical discussion of this approach can be found in [27] by the same authors. Here, we only give the model function Algorithm 1 in the modified L-curve method. There the constants C and T are updated by

Ck ¼ 

ðkFðqðak ÞÞk2 þ ak kqðak Þk2 Þ2 2

kqðak Þk

;

Tk ¼

kFðqðak ÞÞk2 kqðak Þk2

:

ð13Þ

We mention that the Step 5 in Algorithm 1 might be restarted with an iterative ak+1 as

akþ1 ¼ x T k þ

!  d 2 2 ð2 þ 2lÞC k T k 2lHm  ðT k þ aÞ ;  2C k þ 4lC k 2C k þ 4lC k

with a contraction factor x 2 (0, 1). A theoretical justification can be found, for instance, in Remark 3.8 [27].

7868

S. Lu et al. / International Journal of Heat and Mass Transfer 55 (2012) 7865–7872

Note that if a is set to zero, the aforementioned CG iteration turns out to be a particular case (cf. [20]). Thus, the results presented in this section generalize the CG-based solution procedure such that it can also tackle the Tikhonov cost functional (8).

Algorithm 1: Model function approach in modified L-curve [27] Input:  > 0; Hdm ; l > 0. 1: Choose initial guess a1 > a⁄, and set k = 1. 2: Do 3: Solve the minimization problem (8) to find q = q(ak) 4: Update Ck and Tk according to Eq. (13) and construct the corresponding model function

 2 mk ðaÞ ¼ Hdm  þ

Ck

a þ Tk

4. Simulation case studies In this section, we apply the proposed solution approach to two case studies of practical importance. One is the estimation of the laser heat flux on the front surface of a target object. It aims at accurately assessing the resulting thermo-mechanical response after the tests of high-energy laser beams. Note that the sensors cannot be used to directly measure the unknown quantity since they can be easily destroyed or interfered with the laser beam. The other example is the reconstruction of the instantaneous local heat flux on the boiling surface. The estimation results are considered as a crucial prerequisite in boiling heat transfer modeling. Similarly, the direct placement of temperature sensors on the boiling surface is also infeasible due to technical reasons. Thus in both cases the use of efficient inversion approaches for the solution of unmeasurable heat fluxes is indispensable.

:

Insert mk(a) into Eq. (12), Update ak+1 as the minimizer of WC k ;T k ;l ðaÞ, set k :¼ k + 1.     6: While akþ1akak  6 .

5:

7: If the stopping rule is fulfilled return q, ak. 3. Conjugate gradient iteration Step 3 in Algorithm 1 is implemented by using the efficient CG ^kþ1 of q method. The CG iteration calculates a heat-flux estimate q on the boundary C3 according to

^kþ1 ðx; tÞ ¼ q ^k ðx; tÞ  lk pk ðx; tÞ; q

for k ¼ 0; 1; 2; . . . ;

4.1. Case study 1: The estimation of the front-surface laser heat flux The first simulation example for the estimation of front-surface laser heat flux has been recently reported [21]. The target object is assumed to be made of Aluminum Alloy 2024-T6 and its size is 68.4  68.4  4.8 mm3. The simulated thermal properties of the material can be found in [30]. The simulated heat flux q(x, y, t) imposed on the front surface C3 is of Gaussian type in space and explicitly defined by

ð14Þ

0

^0 ¼ 0; p0 ¼ rJ jC . lk is the search step length in each iterawith q 3 tion k. The conjugate search direction pk(x, t) is updated by

pk ðx; tÞ ¼ rJ k ðx; tÞjC3 þ ck pk1 ðx; tÞ;

ð15Þ

where the gradient is determined by

rJk ðx; tÞjC3 ¼ wk ðx; tÞjC3 þ aq^k :

q ¼ qmax  b  expf½ðx  xc Þ2 þ ðy  yc Þ2 =r 2 g  f1 þ sin½2pf

ð16Þ

 ðt  0:1Þg

The conjugate coefficient ck in (15) is determined from

R tf R

c0 ¼ 0; ck ¼ R tf0 R 0

C3 ½ r J

C3 ½rJ

k 2

 dx dt

k1 2

 dx dt

ð17Þ

:

In each iteration, lk and wk are determined by solving the sensitivity problem and the adjoint problem. The sensitivity problem takes the form

8 k > qcp @S@t ¼ r  ðkrSk Þ; > > > > < Sk ð; 0Þ ¼ 0; k > ¼ 0; > k @S > @n > > : @Sk k @n ¼ pk ;

in X  ð0; t f Þ; on X;

ð18Þ

on C1 [ C2  ð0; t f Þ; on C3  ð0; tf Þ

and the search step length is determined by

h

R tf R k

l ¼

0

C1

i

Hðx; t; q^k Þ  Hdm ðx; tÞ Sk ðx; tÞdx dt þ a h

R tf R 0

C1

R tf R 0

C3

^k pk dx dt q ;

i2 Rt R Sk ðx; tÞ dx dt þ a 0f C3 ½pk 2 dx dt

ð19Þ ^k Þ denotes the temperature solution of forward probwhere Hðx; t; q ^k . The resulting adjoint problem corresponds to lem (1)–(5) for q ¼ q

8 k > qcp @w@t ¼ r  ðkrwk Þ; > > > > < wk ð; t f Þ ¼ 0; k > > k @w ¼ 2ðHdm  HÞ; > @n > > : @wk k @n ¼ 0;

in X  ð0; t f Þ; on X; on C1  ð0; t f Þ;

ð20Þ

on ðC2 [ C3 Þ  ð0; tf Þ:

The final-time value problem (20) can be transformed to a standard initial value problem by introducing a new time variable t0 = tf  t.

ð21Þ

on the time interval [0, 0.7] s. qmax is the maximum heat flux at the laser Gaussian beam center, b denotes the surface absorptivity, xc, yc are the center coordinates on the front surface, r stands for the radius of Gaussian laser beam, and f is the frequency of laser heat flux. These parameters are chosen as qmax = 8 W mm2, b = 0.018, xc = yc = 34.2 mm, r = 13.4 mm, and f = 5 Hz according to [21]. For the numerical simulation, the noise-free temperature measurement Hm on C1  [0, 0.7] is obtained by solving the well-posed forward problem (1)–(5) with H0 = 300 K and q1 = q2 = 0 (the boundaries C1 and C2 are assumed adiabatic [21]). The fine spatial mesh used for computation is generated using piece-wise linear finite elements on a tetrahedral grid and the grid sizes in x-, y-directions are chosen as 0.45 mm and that in z-direction is chosen as 0.4 mm. A one-step Crank–Nicolson scheme with a step size 0.01 s is applied for time discretization to gain a second-order accuracy. Such a fine discretization model, a total number of 21,606,507 spatial and temporal points, is sufficient to yield accurate temperature values with negligible error. The simulated noisy measurements are constructed by adding an artificial error d to Hm, namely, Hdm ¼ Hm þ dx, where x is generated from a zero mean normal distribution with variance one. In this case study d = 0.1 is applied. It has to be mentioned that the information of the measurement noise is not required for the solution procedure addressing the inverse problem. To avoid inverse crimes [31], the IHCP is solved on a coarser spatial mesh with the grid sizes in x-, y-directions chosen as 0.9 mm and that in z-direction chosen as 0.8 mm. For running Algorithm 1,  = 0.01, l = 1 and a1 = 0.1 are applied. As shown in Fig. 1, the Tikhonov regularization parameter converges quickly and the optimal Tikhonov regularization parameter is obtained after only five iterations, i.e. a⁄ = 0.048. For lack of space, only

7869

S. Lu et al. / International Journal of Heat and Mass Transfer 55 (2012) 7865–7872

Table 1 Case study 1: A quantitative study of the heat-flux reconstruction for the noisy ^ measurement data Hdm . While q denotes the simulated heat flux (in W mm2), q represents the estimated heat flux (in W mm2). Time

Heat flux

Max value

Min value

Mean value

t=0s

q ^ q q ^ q q ^ q q ^ q q ^ q q ^ q q ^ q q ^ q q ^ q q ^ q q ^ q q ^ q q ^ q q ^ q q ^ q

0 0.0138 0 0.0236 0.1440 0.1620 0.2880 0.2504 0.1440 0.1362 0 0.0431 0.1440 0.1442 0.2880 0.2491 0.1440 0.1440 0 0.0399 0.1440 0.1434 0.2880 0.2520 0.1440 0.1577 0 0.0255 0 0

0 0.0245 0 0.0089 0 0.0085 0 0.0070 0 0.0072 0 0.0077 0 0.0082 0 0.0138 0 0.0080 0 0.0060 0 0.0077 0 0.0091 0 0.0074 0 0.0042 0 0

0 0.0024 0 0.0022 0.0169 0.0188 0.0338 0.0296 0.0169 0.0161 0 0.0041 0.0169 0.0177 0.0338 0.0296 0.0169 0.0163 0 0.0043 0.0169 0.0169 0.0338 0.0291 0.0169 0.0181 0 0.0029 0 0

t = 0.05 s t = 0.1 s t = 0.15 s t = 0.2 s t = 0.25 s Fig. 1. Case study 1: Performance of Algorithm 1 with a1 = 0.1, simulated noisy data Hdm .

 = 0.01

and

t = 0.3 s t = 0.35 s t = 0.4 s t = 0.45 s t = 0.5 s t = 0.55 s t = 0.6 s t = 0.65 s t = 0.7 s

4.2. Case study 2: The estimation of instantaneous local boiling heat flux The second example aims to simulate the dynamics of a ringshaped boiling heat flux occurring in a single-bubble nucleate boiling experiment [32]. According to the experimental settings, the considered stainless steel heating foil is modeled as a 3D spatial domain 1  1  0.05 mm3. The time interval of observation is [0, 0.05] s, which corresponds to a single-bubble cycle. The simulated heat flux (in W mm2) is

qðx; y; tÞ ¼ aðtÞ  bðx; yÞ;

ðx; y; tÞ 2 C3  ½0; 0:05;

ð22Þ

where the temporally varying term is

 aðtÞ ¼

0:1 sinð100pðt  0:02ÞÞ þ 0:1; t 2 ½0:015; 0:035; 0; otherwise

and the spatially varying term is

 bðx; yÞ ¼ Fig. 2. Case study 1: A comparison of simulated heat flux (left column) and estimated heat flux (right column) at representative time instants.

the estimated laser heat fluxes at representative time instants are shown in Fig. 2. As we can see, the estimated heat flux is close to the simulated heat flux and captures the periodic dynamics accurately. Estimation results at other time instants show similar quality (cf. Table 1). In general, for solving such an ill-posed inverse problem, the reconstruction quality of the unknown quantity depends on the level of measurement noise (cf. [32] for a detailed analysis). The small difference between the extreme values of simulated heat flux and those of estimated heat flux (cf. Table 1) is therefore a result of the ill-posedness of the IHCP and the corrupted temperature data available.

0:5sinð5pð0:3  f ðx; yÞÞÞ þ 0:5; f ðx; yÞ 6 0:4; 0;

otherwise;

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi with f ðx; yÞ ¼ ðx  0:5Þ2 þ ðy  0:5Þ2 . Simulated material properties, initial and boundary conditions are chosen the same as those used by [32]. The noise-free temperature measurement is generated by solving the forward problem on a sufficiently fine mesh with grid sizes in x-, y- and z- directions all chosen as 0.005 mm. The inverse problem is solved on a coarser spatial mesh with the grid sizes in x-, y- and z- directions all chosen as 0.01 mm. The time step 0.001 s is applied for the time discretization. Parameters  = 0.01, l = 1 and a1 = 103 are applied to run Algorithm 1. The optimal Tikhonov regularization parameter is obtained after four iterations (a⁄ = 3.99  104; cf. Fig. 3). The computational efficiency of the method proposed in this paper is comparable to the CG-based iterative regularization method that considered only

7870

S. Lu et al. / International Journal of Heat and Mass Transfer 55 (2012) 7865–7872 Table 2 ^A (in W mm2) Case study 2: A quantitative study of the heat-flux reconstruction q ^B (in W mm2) obtained in the present study and the heat-flux reconstruction q obtained by Heng et al. [32]. Time

Heat flux

Max value

Min value

Mean value

t=0s

q ^A q ^B q q ^A q ^B q q ^A q ^B q q ^A q ^B q q ^A q ^B q q ^A q ^B q q ^A q ^B q q ^A q ^B q q ^A q ^B q q ^A q ^B q q ^A q ^B q

0 0.0089 0.0131 0 0.0065 0.0094 0 0.0043 0.0060 0 0.0046 0.0061 0.1000 0.0988 0.1033 0.2000 0.1953 0.2055 0.1000 0.0998 0.1049 0 0.0071 0.0093 0 0.0046 0.0060 0 0.0050 0.0068 0 0 0

0 0.0077 0.0107 0 0.0054 0.0075 0 0.0056 0.0076 0 0.0055 0.0077 0 0.0072 0.0111 0 0.0140 0.0166 0 0.0091 0.0092 0 0.0054 0.0070 0 0.0054 0.0061 0 0.0051 0.0072 0 0 0

0 3.1576  105 1.0439  104 0 1.2850  104 2.2184  104 0 2.6086  104 2.4568  104 0 1.5918  104 8.4000  106 0.0246 0.0248 0.0248 0.0493 0.0491 0.0492 0.0246 0.0248 0.0249 0 1.7190  104 1.8532  104 0 1.5644  104 2.7386  104 0 2.4816  105 1.9194  104 0 0 0

t = 0.005 s

t = 0.01 s

t = 0.015 s

Fig. 3. Case study 2: Performance of Algorithm 1 with a1 = 0.001, simulated noisy data Hdm .

 = 0.01

and

t = 0.02 s

t = 0.025 s

data error in least squares and applied the L-curve stopping criterion [32]. While the oscillation of heat-flux estimates obtained in the present study is slightly smaller, the obtained estimation results by using the two different methods are similar (cf. Table 2 and Fig. 4). However, for tackling the ill-posedness, the solution method proposed in this paper is considered more robust since the L-curve criterion used in previous approach has weak theoretical justification.

5. Heat flux estimation with experimental temperature data In the nucleate boiling, temperature superheat is large enough to produce small bubbles along the heating surface. As temperature superheat further increases, it is sufficient to sustain the creation of the bubbles such that they depart and rise through the liquid regardless of the condensation rate until the critical heat flux is reached. Many literature emphasizes the importance of nucleate boiling heat transfer in a wide range of industrial processes. However, its underlying physical transport mechanisms is still not fully understood yet. Among all kinds of studies, the accurate estimation of local boiling heat fluxes are fundamental but very important. In this section we focus on real problems arising from singlebubble nucleate boiling experiments operated by our partners [33]. There, thin stainless steel heating foils (50 and 25 lm) were used to conduct two different nucleate boiling experiments. The refrigerant HFC 7100 was chosen as the boiling fluid. The explicit values of thermal material properties according to [32] are applied for the computation in the present study. For the first experiment, a 3D spatial domain X :¼ 2.416  2.288  0.05 mm3 is used to represent the heating foil. The chosen space and time discretization follows the spatial and temporal resolution of measuring system (152  144 pixels with a pixel size of 16  16 lm in the x-y-planes, and a frame rate of 987 Hz in time). The supplied heat flux on the heated boundary C1 is set to 5356 W m2. For the second experiment, a 3D spatial domain X :¼ 5.104  4.08  0.025 mm3 is considered, and the space and time discretization also follows the spatial and temporal resolution of measuring system (320  256 pixels with a pixel size of 16  16 lm in the xy-planes, and a frame rate of 987 Hz in time). Differently, a relatively larger input heat flux, q1 = 12,900 W m2, is applied in this case. In both cases, the temperature fields measured at the first time frame are used as the initial temperature distributions by assuming equal values across the foil thickness due to lack of information. For the experiment using the heating foil with a thickness of 50 lm, the optimal Tikhonov regularization parameter is obtained

t = 0.03 s

t = 0.035 s

t = 0.04 s

t = 0.045 s

t = 0.05 s

after ten iterations, namely a⁄ = 0.00437, and the estimated boiling heat flux (in W m2) is shown in Fig. 5. From the estimation results, it is apparent that the transient boiling heat flux had a significant change during a single bubble cycle. At time frames 22, 23, 24, the vapor bubble was formed and grew, and resulted in an increasingly larger Gaussian-type high heat-flux distribution. At subsequent time frames the bubble continued to expand and rise, and finally detached from the contact surface completely. During the process of detaching, the estimated ring-shaped high heat-flux region had a peak value that was nearly 30 times larger than the supplied input heat flux. These heat-flux estimates, which could not be measured during the experiment directly, confirm the microlayer theory [34]. The theory states that most of the heat during boiling is transferred in the micro-region of the three-phase contact line by evaporation. For the experiment using the heating foil with a thickness of 25 lm, the optimal Tikhonov regularization parameter is obtained after nine iterations, namely a⁄ = 7.046  104. Fig. 6 shows the estimated heat fluxes at a few representative time frames. While the ring-shaped contour is much clearer than that for the first experiment, a similar high heat-flux pattern is observed. Furthermore, it is worth to note that, by far this is the first reconstructed solution for the 25 lm measurement data by means of our 3D IHCP solution technique, which is independent of the chosen boiling fluid and applicable to heated bodies of arbitrary geometry.

6. Conclusions In recent years, the authors at RWTH Aachen University have published research work to address the issue of efficient solution

S. Lu et al. / International Journal of Heat and Mass Transfer 55 (2012) 7865–7872

7871

Fig. 4. Case study 2: A comparison of simulated heat flux (left column), estimated heat flux in the present study (middle column) and estimated heat flux obtained in [32] (right column).

Fig. 5. Real experiment 1: The estimated boiling heat fluxes (in W m2) at time frames 21–36 (from top to bottom, left to right).

of 3D IHCP arising in pool boiling by means of advanced computational techniques involving the CG-based iterative regularization, the multi-level computational strategy and the adaptive mesh refinement via a posteriori error estimation [15,35]. Those ap-

proaches depend on the choice of measurement error bound d and focus on improving the computational efficiency of the solution procedure by selecting smaller discretization models in an automated manner.

7872

S. Lu et al. / International Journal of Heat and Mass Transfer 55 (2012) 7865–7872

Fig. 6. Real experiment 2: The estimated boiling heat fluxes (in W m2) at representative time frames 14, 15, 16.

In this paper, we have presented a robust and efficient approach for the solution of large-scale 3D IHCP on the continuous domain. It uses enhanced Tikhonov regularization via an efficient CG iteration to ensure both robustness and computational efficiency. The solution method is independent of the choice of discretization models and robust in practical applications since no information of measurement error is needed for the solution procedure addressing the inverse problem. The proposed solution method has been validated by two interesting practical examples and also successfully applied to the estimation with experimental temperature data obtained from single-bubble nucleate boiling experiments. The direct application of the method to a wide range of similar 3D IHCP is straightforward. Acknowledgments Dr. S. Lu is supported by the National Natural Science Foundation of China (Key project Nos.91130004 and 11101093) and Shanghai Science and Technology Commission (Nos.11ZR1402800 and 11PJ1400800). He is currently a research fellow of the Alexander von Humboldt Foundation. Dr. Y. Heng is very grateful to Prof. W. Marquardt for providing research facilities for his scientific work. The authors particularly thank Prof. P. Stephan and other colleagues for providing the data of their boiling experiments. References [1] J.V. Beck, Calculation of surface heat flux from an internal temperature history, ASME Paper, 1962, 62–HT–46. [2] N.V. Shumakov, A method for the experimental study of the process of heating a solid body, Sov. Phys. Tech. Phys. 2 (1957) 771–781. [3] G. Stolz, Numerical solution to an inverse problem of heat conduction for simple shapes, J. Heat Transfer 82 (1960) 20–26. [4] A.N. Tikhonov, V.Y. Arsenin, Solutions of Ill-Posed Problems, V.H. Winston and Sons, Washington, 1977. [5] M. Raynaud, J. Bransier, A new finite-difference method for the nonlinear inverse heat conduction problem, Numer. Heat Transfer (Part A) 9 (1986) 27– 42. [6] J.V. Beck, B. Blackwell, C.S. Clair, Inverse Heat Conduction. Ill-Posed Problems, Wiley, New York, 1985. [7] O.M. Alifanov, Inverse Heat Transfer Problems, Springer, Berlin, 1994. [8] Y. Jarny, M.N. Ozisik, J.P. Bardon, A general optimization method using adjoint equation for solving multidimensional inverse heat conduction, Int. J. Heat Mass Transfer 34 (1991) 2911–2919. [9] B. Kaltenbacher, A. Neubauer, O. Scherzer, Iterative Regularization Methods for Nonlinear Ill-Posed Problems, de Gruyter, Berlin, 2008. [10] J.B. Wang, N. Zabaras, A bayesian inference approach to the inverse heat conduction problem, Int. J. Heat Mass Transfer 47 (2004) 3927–3941. [11] N. Tian, J. Sun, W.B. Xu, C.-H. Lai, Estimation of unknown heat source function in inverse heat conduction problems using quantum-behaved particle swarm optimization, Int. J. Heat Mass Transfer 54 (2011) 4110–4116. [12] G.J. Wang, L. Zhu, H. Chen, A decentralized fuzzy inference method for solving the two-dimensional steady inverse heat conduction problem of estimating boundary condition, Int. J. Heat Mass Transfer 54 (2011) 2782–2788.

[13] P.F.B. Sousa, S.R. Carvalho, G. Guimaraes, Dynamic observers based on green’s functions applied to 3D inverse thermal models, Inverse Prob. Sci. Eng. 16 (2008) 743–761. [14] M. Soemers, Numerical methods for the solution of a three-dimensional anisotropic inverse heat conduction problem, Ph.D. thesis, RWTH Aachen University, 2008. [15] H. Egger, Y. Heng, W. Marquardt, A. Mhamdi, Efficient solution of a threedimensional inverse heat conduction problem in pool boiling, Inverse Prob. 25 (2009) 095006. 19pp. [16] A.P. Fernandes, P.F.B. Sousa, V.L. Borges, G. Guimaraes, Use of 3D-transient analytical solution based on green’s function to reduce computational time in inverse heat conduction problems, Appl. Math. Model. 34 (2010) 4040–4049. [17] O.M. Alifanov, A.V. Nenarkomov, Three-dimensional boundary inverse heat conduction problem for regular coordinate systems, Inverse Prob. Eng. 7 (1999) 335–362. [18] C.-H. Huang, S.-P. Wang, A three-dimensional inverse heat conduction problem in estimating surface heat flux by conjugate gradient method, Int. J. Heat Mass Transfer 42 (1999) 3387–3403. [19] S. Groß, M. Soemers, A. Mhamdi, F.A. Sibai, A. Reusken, W. Marquardt, U. Renz, Identification of boundary heat fluxes in a falling film experiment using high resolution temperature measurements, Int. J. Heat Mass Transfer 48 (2005) 5549–5562. [20] Y. Heng, A. Mhamdi, S. Groß, A. Reusken, M. Buchholz, H. Auracher, W. Marquardt, Reconstruction of local heat fluxes in pool boiling experiments along the entire boiling curve from high resolution transient temperature measurements, Int. J. Heat Mass Transfer 51 (2008) 5072–5087. [21] J.H. Zhou, Y.W. Zhang, J.K. Chen, Z.C. Feng, Inverse estimation of surface heating condition in a three-dimensional object using conjugate gradient method, Int. J. Heat Mass Transfer 53 (2010) 2643–2654. [22] J. Hadamard, Lectures on Cauchy’s Problem in Linear Partial Differential Equations, Yale University Press, New Haven, 1923. [23] L.C. Evans, Partial differential equations, Graduate Studies in Mathematics, vol. 19, American Mathematical Society, Providence, RI, 1998. [24] O.A. Ladyzenskaja, V.A. Solonnikov, N.N. Ural’ceva, Linear and quasilinear equations of parabolic type, Translations of Mathematical Monographs, vol. 23, American Mathematical Society, 1968. [25] H.W. Engl, M. Hanke, A. Neubauer, Regularization of Inverse Problems, Kluwer Academic Publishers, Dordrecht, 1996. [26] P.C. Hansen, Rank-Deficient and Discrete Ill-posed Problems: Numerical Aspects of Linear Inversion, SIAM, Philadelphia, 1998. [27] Y. Heng, S. Lu, A. Mhamdi, S.V. Pereverzev, Model functions in the modified Lcurve method – case study: the heat flux reconstruction in pool boiling, Inverse Prob. 26 (2010) 055006. 13pp. [28] T. Regin´ska, A regularization parameter in discrete ill-posed problems, SIAM J. Sci. Comput. 17 (1996) 740–749. [29] K. Kunisch, J. Zou, Iterative choices of regularization parameters in linear inverse problems, Inverse Prob. 14 (1998) 1247–1264. [30] F.P. Incropera, D.P. Dewitt, T.L. Bergman, A.S. Lavine, Fundamentals of Heat and Mass Transfer, sixth ed., Wiley, New York, 2007. [31] J.P. Kaipio, E. Somersalo, Statistical and Computational Inverse Problems, Springer Verlag, New York, 2004. [32] Y. Heng, A. Mhamdi, E. Wagner, P. Stephan, W. Marquardt, Estimation of local nucleate boiling heat flux using a three-dimensional transient heat conduction model, Inverse Prob. Sci. Eng. 18 (2010) 279–294. [33] E. Wagner, A. Sprenger, P. Stephan, O. Koeppen, F. Ziegler, H. Auracher, Nucleate boiling at single artificial cavities: bubble dynamics and local temperature measurements, in: Proc. sixth Int. Conf. Multiphase Flow, 2007. [34] P. Stephan, J. Hammer, A new model for nucleate boiling heat transfer, Wärmeund Stoffübertragung 30 (1994) 119–125. [35] Y. Heng, A. Mhamdi, W. Marquardt, Efficient reconstruction of local heat fluxes in pool boiling experiments by goal-oriented adaptive mesh refinement, Heat Mass Transfer 46 (2010) 1121–1135.