Local heat flux estimation inside tubes through conjugate gradient method with adjoint operator: application to the pulsating heat pipes case

Abstract

Purpose

The purpose of this paper is to apply the conjugate gradient (CG) method, together with the adjoint operator (AO) to the pulsating heat pipe problem, including some quite interesting experimental results. The CG method, together with the AO, was able to estimate the unknown functions more efficiently than the other techniques presented in this paper. The estimation of local heat transfer coefficients, rather than the global ones, in pulsating heat pipes is a relatively new subject and presenting a robust, efficient and self-regularized inverse tool to estimate it, supported also by some experimental results, is the main purpose of this paper. To also increase the visibility and the general use of the paper to the heat transfer community, the authors include, as supplemental material, all numerical and experimental data used in this paper.

Design/methodology/approach

The approach was established on the solution of the inverse heat conduction problem in the wall by using as starting data the temperature measurements on the outer surface. The procedure is based on the CG method with AO. The here proposed approach was first verified adopting synthetic data and then it was validated with real cases regarding pulsating heat pipes.

Findings

An original fast methodology to estimate local convective heat flux is proposed. The procedure has been validated both numerically and experimentally. The procedure has been compared to other classical methods presenting some peculiar benefits.

Practical implications

The approach is suitable for pulsating heat pipes performance evaluation because these devices present a local heat flux distribution characterized by an important variation both in time and in space as a result of the complex flow patterns that are generated in this type of devices.

Originality/value

The procedure here proposed shows these benefits: it affords a general model of the heat conduction problem that is effortlessly customized for the particular case, it can be applied also to large datasets and it presents reduced computational expense.

Citation

Colaço, M., Bozzoli, F., Cattani, L. and Pagliarini, L. (2023), "Local heat flux estimation inside tubes through conjugate gradient method with adjoint operator: application to the pulsating heat pipes case", International Journal of Numerical Methods for Heat & Fluid Flow, Vol. 33 No. 5, pp. 1754-1774. https://doi.org/10.1108/HFF-09-2022-0547

Nomenclature

SymbolQuantitySI unit

c_p

Specific heat at constant pressure;       J/kg · K

d

Direction of descent;

g

Source term;K/m²

h

Heat transfer coefficient;       W/m² · K

k

Thermal conductivity        W/m · K

M

Number of measurements on space;

q

Convective heat flux per unit surface;  W/m²

r

Radial coordinate;          m

S

Objective function;

t

Time;              s;

T

Temperature;            K, °C

z

Axial coordinate;          m

α

Thermal diffusivity;        m²/s

β

Search-step size;

ϵ

Tolerance for stopping criteria;

λ

Lagrange multiplier;

ρ

Density;              kg/m³

σ

Noise level; and              K, °C

γ

Coefficient of conjugation.

1. Introduction

A tricky issue in several engineering cases comes from the need of capturing what occurs within a heat transfer device by measuring a physical quantity on the external surfaces.

It is a very powerful task when it is not conceivable to perform invasive measurements to examine the inner parts of the device, for technical or security causes, as, for example, in nuclear plants or in any inner flow case when the inner walls could not be directly investigated. A simple answer to this challenging task is provided by devising an inverse heat conduction problem (IHCP) in the boundary walls of the apparatus, by measuring, conceivably adopting contactless measurements, the temperature on the outer wall. Between the numerous imaginable implementations of this procedure, the one that aims at assessing the heat flux or the convective heat transfer coefficient distribution on the inner surface of a tube, adopting as initial data the outer surface temperature, has demonstrated to be worthwhile between experimental heat exchange studies.

Some profitable adoptions of this technique have been proposed for the estimation of local heat transfer between the internal fluid and the walls in pulsating heat pipes (PHPs) (Cattani et al., 2019; Pagliarini et al., 2021a). In these devices, indeed, the convective heat flux significantly changes both spatially and temporally, as a consequence of the flow patterns that develop because of the complex thermofluidic interaction between the walls and the two-phase fluid inside the device (Pagliarini et al., 2021b). PHPs are reasonably new two-phase passive heat transfer devices that are obtaining resounding interest because of their reduced dimensions, versatility, capability to perform in absence of external energy inputs and the adaptability to run in absence of gravity. Even though the substantial benefits, the PHPs’ working principles are not yet completely clarified. To analyze and profoundly comprehend their governing phenomena, several investigations have been performed in the past years (Mangini et al., 2018; Zhang and Faghri, 2008; Han et al., 2016). Among them, some recent studies focused on local thermal behavior of these devices because its understanding is crucial for the PHP design (Cattani et al., 2019; Pagliarini et al., 2021a; Mameli et al., 2014; Jo et al., 2019). Data concerning the local wall to fluid heat flux can give an improved awareness of the heat transfer increase mechanisms and the causal correlation among it and the produced oscillating flow.

Concerning IHCPs, many solving techniques have been proposed (Beck and Arnold, 1977) to assess the local heat flux in ducts. Between these methods, the function specification methods (Beck and Arnold, 1977; Beck et al., 1985), iterative methods (Alifanov, 1974; Alifanov, 1994; Orlande et al., 2011), filtering techniques (Murio, 1993; Bozzoli et al., 2013; Massard et al., 2012; Bozzoli and Rainieri, 2011) and regularization method (Kim and Lee, 2002; Tikhonov and Arsenin, 1977; Hansen, 1990) are present.

In the present work it is proposed and validated, numerically and experimentally, a promising procedure for the assessing of the local heat flux at the inner wall of pipes by using as input data the temperature measured on the outer wall. The adopted method is based on the conjugate gradient (CG) method with adjoint operator (AO). There are two main advantages of using the CG with AO: the first one relies on the fact that the method works on the function space. Therefore, the entire unknown function is estimated, rather than a series of parameters. The second one is because the adjoint problem gives an approximation of the gradient of the functional being minimized, therefore avoiding the necessity of calculating very large sensitivity matrices. Besides that, using the Morozov’s discrepancy principle, the method itself is self-regularized. The implementation of this method is encouraging because it allows to avoid invasive measurements, the computational time and cost are significantly small and that is not restricted, at least for the practical applications shown here, by the amount of data used and by the thermal properties of the material.

2. Direct problem

In this paper, for the development of the estimation approaches, we consider two different models for the physical problem: a thin-wall model and a complete model (Cattani et al., 2019; Pagliarini et al., 2021a). Considering a long duct, the thin-wall model assumes that axial diffusion is more important than in the radial direction, lumping the temperature on the r direction. This is valid when the heat conduction within the solid object is stronger than the heat convection to fluid flowing on the object walls.

Both models consider a hollow cylinder (Figure 1), where the internal boundary condition (r = r_i) consists of an applied heat flux because of the internal convection to the inner fluid (a quantity whose estimation is the objective of this paper), and the external boundary condition (r = r_e) can be represented as a combined convection and radiation boundary condition to the environment. To check the validity of the thin-wall model, we will present a comparison against both models for the materials considered in this paper, to verify the discrepancy between them.

Two different sets of boundary conditions on the end lengths of the hollow cylinder (z = 0 and z = L) are considered: adiabatic on both ends and prescribed temperature on both ends. Both models also consider that the heat conduction in the circumferential direction is small compared with the ones on the radial and longitudinal directions.

The complete model, considering adiabatic boundary conditions on z = 0 and z = L, can then be written as:

(1.a) 1α∂T∂t=1r∂∂r(r∂T∂r)+∂2T∂z2 in ri<r<re; 0<z<L;t>0

(1.b) −k∂T∂r=q(z,t) on r=ri; 0<z<L;t>0

(1.c) −k∂T∂r=h(T−T∞) on r=re; 0<z<L;t>0

(1.d) ∂T∂z=0 on ri<r<re;z=0;t>0

(1.e) ∂T∂z=0 on ri<r<re;z=L;t>0

(1.f) T=T0(r,z,0) in ri<r<re; 0<z<L;t=0

where h is a combined heat transfer coefficient, because of convection and radiation to the environment, k is the thermal conductivity of the hollow cylinder and T_∞ is the bulk temperature of the external fluid.

For the thin-wall model, we can integrate equation (1.a) on the radial direction considering that T˜def¯¯T(z,t), use the boundary conditions equations (1.b) and (1.c), therefore lumping the temperature on r (Pagliarini et al., 2021a) and obtain:

(2.a) 1α∂T˜∂t=∂2T˜∂z2+CT˜+g(z,t) in 0<z<L;t>0

(2.b) ∂T˜∂z=0 on z=0;t>0

(2.c) ∂T˜∂z=0 on z=L;t>0

(2.d) T˜=T˜0(z,0) in 0<z<L;t=0

where the constant C is given as:

(3) C=−2rehk(re2−ri2)

and the heat flux on r = r_i now appears as a source term of equation (2.a), defined as:

(4) g(z,t)=2k(re2−ri2)[rehT∞−riq(z,t)]

These two models will be used in the next section to estimate the unknown quantity q(z, t) on the inner radius r = r_i. For the inverse problem point of view, however, it is interesting to notice that the full model involves a boundary condition estimation, whereas the thin-wall model involves a source term estimation.

3. Inverse problem solution

In the current work, we used the CG method with AO. This is a well-established method (Beck and Arnold, 1977; Beck et al., 1985; Alifanov, 1974; Tikhonov and Arsenin, 1977; Hansen, 1990) where the gradient of the objective function and the search-step size are given as the solutions of two auxiliary problems, like the original physical model.

Supposing that some measurements Y(z_m, t) and some estimated temperatures T(z_m, t;q), obtained by a model of the physical problem, are available at r = r_e, we can define the following objective function, which must be minimized to make the physical model match the experimental data:

(5) S[q(z,t)]=12∫t=0tf∑m=1M[T(zm,t;q)−Y(zm,t)]2

where t_f is the total time of measurements, M is the number of measurements on space and the subscript m refers to the position of the measurement. On this expression, it is emphasized that the temperature T(z, t), taken on the external radius, which comes from the physical model, is a function of the unknown heat flux q(z, t).

The method starts by proposing the following iterative procedure for the estimation of the unknown heat flux q(z, t) (Beck and Arnold, 1977; Beck et al., 1985; Alifanov, 1974, Tikhonov and Arsenin, 1977; Hansen, 1990):

(6) qk+1(z,t)=qk(z,t)+βkdk(z,t)

where β^k is the search-step size and d^k(z, t) is the direction of descent. There are several versions of the CG method, depending on the choice of d^k(z, t) (Hansen, 1990). In this paper, we choose the Fletcher–Reeves (Özisik and Orlande, 2021) version of it, where d^k(z, t) is given by:

(7) dk+1(z,t)=−∇S[qk(z,t)]+γkdk−1(z,t)

and the coefficient of conjugation γ is given as:

(8) γk=∫t=0tf∫z=0L‖∇S[qk(z,t)]‖2dzdt∫t=0tf∫z=0L‖∇S[qk−1(z,t)]‖2dzdt for k=1,2,

(9) γk=0 for k=0

In the CG method with AO the search-step size is achieved by the solution of an auxiliary sensitivity problem and the gradient of the objective function is attained from the solution of another auxiliary problem, called adjoint problem. These two auxiliary problems depend on the physical model used to represent T(z_m, t;q) in equation (5). Because in this paper we consider two different models: a thin-wall and a complete model, two different sets of auxiliary problems will be obtained. Those are described next.

3.1 Sensitivity problem

The sensitivity problem is used to define the search-step size β^k that appears in equation (6). If the heat flux q(z, t) is perturbed by a quantity εΔq(z, t), we can suppose that the temperature of the complete model T(r, z, t) will also be perturbed by some quantity ξΔT(r, z, t).

Writing the direct problem as an operator F[.], we can obtain the sensitivity problem applying the following limiting process:

(10) lim⁡ε→0F[q(z,t)+εΔq(z,t)]−F[q(z,t)]ε=0

Therefore, for the complete model, we obtain:

(11.a) 1α∂ΔT∂t=1r∂∂r(r∂ΔT∂r)+∂2ΔT∂z2 in ri<r<re; 0<z<L;t>0

(11.b) k∂ΔT∂r=Δq(z,t) on r=ri; 0<z<L;t>0

(11.c) −k∂ΔT∂r=h(ΔT) on r=re; 0<z<L;t>0

(11.d) ∂ΔT∂z=0 on ri<r<re;z=0;t>0

(11.e) ∂ΔT∂z=0 on ri<r<re;z=L;t>0

(11.f) ΔT=0 in ri<r<re; 0<z<L;t=0 

whereas for the thin-wall model, we have:

(12.a) 1α∂ΔT˜∂t=∂2ΔT˜∂z2+CΔT˜+Δg(z,t) in 0<z<L;t>0

(12.b) ∂ΔT˜∂z=0 on z=0;t>0

(12.c) ∂ΔT˜∂z=0  on z=L;t>0

(12.d) ΔT˜=0  in 0<z<L;t=0

The search-step size β^k can be obtained by minimizing the objective function provided by equation (5) with regard to it, at some iteration k + 1:

(13) min⁡βkS[qk+1(z,t)]=min⁡βk12∫t=0tf∑m=1M[T(zm,t;qk+1)−Y(zm,t)]2

where T(z_m, t;q^k+1) is the temperature obtained either from the thin-wall or from the complete model on r = r_e for the measurement locations z_m. However, q^k+1 is given by equation (6) where β^k appears explicitly. Therefore, introducing equation (6) into equation (13) and applying a first-order Taylor expansion to linearize T(z_m, t;q^k+1) and make β^k appears explicitly in equation (13), we obtain, after making ∂S[qk+1(z,t)]∂βk=0:

(14) βk=−∫t=0tf∑m=1M[T(zm,t;qk)−Y(zm,t)]ΔT(zm,t;qk)∫t=0tf∑m=1M[ΔT(zm,t;qk)]2

where, if the complete model is being considered, ΔT is the solution of the sensitivity problem provided by equation (11) considering Δq(z, t) = d^k(z, t), the direction of descent appearing in equation (7), and if the thin-wall model is considered, ΔT is the solution of equation (12) considering Δg(z, t) = d^k(z, t).

3.2 Adjoint problem

The adjoint problem is adopted to determine the gradient of the objective function ∇S[q^k(z, t)] that appears in equations (7) and (8). To do that, we can rewrite the original objective function, given by equation (5) as an extended Lagrangian, where in this case the minimization of equation (5) cannot produce any solution, but instead must result in an estimated heat flux that satisfies the solution of the direct problem. Therefore, the minimization problem is transformed into a constrained one. The Lagrangian for the complete and thin-wall model can be written, respectively, as:

(15) J[q(z,t)]=12∫t=0tf∑m=1M[T(zm,t;q)−Y(zm,t)]2+∫t=0tf∫z=0L∫r=rire[1α∂T∂t−1r∂∂r(r∂T∂r)−∂2T∂z2]λrdrdzdt

(16) J˜[q(z,t)]=12∫t=0tf∑m=1M[T˜(zm,t;q)−Y(zm,t)]2+∫t=0tf∫z=0L[ 1α∂T˜∂t−∂2T˜∂z2−CT˜−g(z,t)]λ˜dzdt

where λ is a Lagrange multiplier. The integral in time can be transformed into an integral in both space and time by considering a Dirac delta function, such as, considering the complete model as an example:

(17) 12∫t=0tf∑m=1M[T(zm,t;q)−Y(zm,t)]2=12∫t=0tf∫z=0L∫r=rire∑m=1M[T(zm,t;q)−Y(zm,t)]2δ(z−zm)δ(r−re)rdrdzdt

It is possible now to specify the directional derivative of J[q(z, t)] in the direction of the perturbation in the heat flux as (the same procedure is valid for J˜[q(z,t)]):

(18) ΔJq[q(z,t)]=lim⁡ε→0J[q(z,t)+εΔq(z,t)]−J[q(z,t)]ε

in such a way that the minimization of equation (16) is obtained when the limiting process of equation (18) is equal to zero. After a lengthy procedure, where the boundary conditions of both the direct and the sensitivity problems are used, we can find the subsequent adjoint problem for the complete model:

1α∂λ∂τ=1r∂∂r(r∂λ∂r)+∂2λ∂ζ2+∑m=1M[Y(zm,t)−T(zm,t;q)]δ(z−zm)δ(r−re)

(19.a) in ri<r<re; 0<ζ<L;t>0

(19.b) ∂λ∂r=0  on r=ri; 0<ζ<L;τ>0

(19.c) −k∂λ∂r=hλ  on r=re; 0<ζ<L;τ>0

(19.d) ∂λ∂ζ=0  on ri<r<re;ζ=0;τ>0

(19.e) ∂λ∂ζ=0  on ri<r<re;ζ=L;τ>0

(19.f) λ=0  in ri<r<re; 0<ζ<L;τ=0

and also, for the thin-wall model:

1α∂λ˜∂τ=∂2λ˜∂ς2+Cλ˜+∑m=1M[Y(zm,t)−T(zm,t;q)]δ(z−zm)

(20.a) in 0<ζ<L;τ>0

(20.b) ∂λ˜∂ζ=0  on ζ=0;τ>0

(20.c) ∂λ˜∂ζ=0 on ζ=L;τ>0

(20.d) λ˜=0 in 0<ζ<L;τ=0

where:

(21.a) τ=tf−t

(21.b) ζ=L−z

During the above procedure, the following term remain for the complete model:

(22.a) ΔJq[q(z,t)]=∫t=0tf∫z=0LΔqkλr|r=ridzdr

By assuming that q(z, t) belongs to the Hilbert space of square integrable functions in the domain (0,t_f)x(0,L), it is possible to write:

(22.b) ΔJq[q(z,t)]=∫t=0tf∫z=0L∇S[q(z,t)]Δq(z,t)dzdr

Therefore, comparing equations (22.a) and (22.b), we obtain, for the complete model:

(22.c) ∇S[q(z,t)]=λrk|r=ri

Following a similar procedure, we can obtain, for the thin-wall model:

(23) ∇S[q(z,t)]=−λ˜(z,t)

3.3 Boundary condition of prescribed temperature on z = 0 and z = L

The previous procedure can be repeated for the case where, instead of adiabatic boundary conditions on z = 0 and z = L, a prescribed temperature is specified. This can be the case where the axial conduction on the end walls of the pipe is important and, therefore, some measured temperatures should be used instead of considering that the end walls are in a region away from the major heat transfer section. Therefore, just for the completeness of the theory, the following equations present the modifications necessary:

1. Equations (1.d) and (e) become:

(24.a) T=Ta(t)  on ri<r<re;z=0;t>0

(24.b) T=Tb(t)  on ri<r<re;z=L;t>0

2. Equations (2.b) and (c) become:

(25.a) T˜=Ta(t)  on z=0;t>0

(25.b) T˜=Tb(t)  on z=L;t>0

3. Equations (11.d) and (e) become:

(26.a) ΔT=0  on ri<r<re;z=0;t>0

(26.b) ΔT=0  on ri<r<re;z=L;t>0

4. Equations (12.b) and (c) become:

(27.a) ΔT˜=0  on z=0;t>0

(27.b) ΔT˜=0  on z=L;t>0

5. Equations (19.d) and (e) become:

(28.a) λ=0  on ri<r<re;ζ=0;τ>0

(28.b) λ=0  on ri<r<re;ζ=L;τ>0

6. Equations (20.b) and (c) become:

(29.a) λ˜=0  on ζ=0;τ>0

(29.b) λ˜=0  on ζ=L;τ>0

3.4 Stopping criteria

We stop the iterative process of the CG method once the functional provided by equation (5) is sufficiently small:

(30) S[q(z,t)]<ϵ

If the measurements are considered without any experimental error, it is possible to define ∈ as a reasonably small number. Nevertheless, the real measurements are not noisy-free, which will produce an unstable inverse problem solution as the assessed temperatures become closer to those experimentally obtained. This complication could be lessened by applying the Morozov’s discrepancy principle (Alifanov, 1974; Alifanov, 1994; Orlande et al., 2011) to end the iterative procedure and to provide the CG method with the required regularization for a robust solution.

The discrepancy principle states that IHCP solution is adequately precise if the residuals among estimated and experimental temperatures are of the order of the standard deviation σ of the measurements. The tolerance ∈ is then obtained as:

(31) ϵ=12Mσ2tf

3.5 Iterative procedure

The subsequent iterative method summarizes the steps needed to apply the CG method with AO, used for the complete or the thin-wall model:

STEP 1: Using an estimate of q(z, t), solve the direct problem provided by equation (1) or equation (2) to achieve the assessed temperatures T(r, z, t) or T˜(z,t) on r = r_e.

STEP 2: Test the stopping criteria given by equation (30). Go on if not verified.

STEP 3: Solve the adjoint problem given by equation (11) or equation (12) to achieve the Lagrange multiplier λ(r, z, t) or λ˜(z,t).

STEP 4: Calculate the gradient of the function ∇S[q(z, t)] from equation (22.c) or equation (23).

STEP 5: Determine the direction of descent d^k+1(z, t) with equation (7).

STEP 6: Solve the sensitivity problem given by equation (11) or equation (12) and obtain the search-step size β^k from equation (14).

STEP 7: Calculate the new estimate q^k+1(z, t) using equation (6) and return to STEP 1.

If the objective function given by equation (5) is convex, the initial guess used on STEP 1 has little influence on the estimate, except by the fact that the method can reach the minimum faster or slower depending on how close to the optimum value the initial guess was. In all cases presented in this paper, the initial guess was adopted as zero and the method converged satisfactorily to the minimum of equation (5), as it can be verified by the good results obtained.

4. Results and discussion

4.1 Introduction

In this paper, we applied the above procedure to the estimate of local heat fluxes inside circular pipes undergoing convection. As applicative test, the case of PHPs has been adopted because in these devices the local heat fluxes strongly vary and change with high frequencies in time. They therefore represent a challenging test to evaluate the estimation performance of the proposed procedure. In particular, two cases regarding PHP devices with different performances and pipe materials were considered. In the first case, the walls were made by aluminum, which is a highly thermal conductive material and, therefore, the thin-wall model can be correctly applied. On the other hand, the second case concerned the study of a sapphire glass tube, where large temperatures gradients are present across the wall and the complete model would be expected to give better results. These two cases were studied by Pagliarini et al. (2021a) and Cattani et al. (2019), respectively. Both studies aimed at estimating the local heat fluxes exchanged between the working fluid and the PHPs walls. In the former work, a Gaussian filter was adopted during the inverse problem analysis and adiabatic boundary conditions were assumed on z = 0 and z = L, whereas a Tikhonov zero order regularization scheme was adopted in the later together with prescribed temperature boundary conditions on z = 0 and z = L.

4.2 Synthetic results

The results obtained by the CG method will be compared with the ones obtained by the techniques presented in Cattani et al. (2019) and Pagliarini et al. (2021a), adopting both synthetic and experimental data. Synthetic data were obtained by solving the complete direct problem defined by equation (1) through the implementation of a numerical model within COMSOL Multiphysics^® environment with an established distribution of convective heat flux (q_exact) at the inner wall. The obtained temperature distribution, spoiled by a random noise, was subsequently used as the starting data of the IHCP. In line with the awaited outcomes in the real conditions, q was assumed to be characterized by a significant variance along time and the axial coordinate:

(32) qexact(z,t)=A cos⁡(pπt)[1+(zΔz+1) L8]

where the constants A and p varied according to the case analyzed.

Gaussian random noises were introduced in the synthetic temperatures obtained from the numerical models, according to the following equation:

(33) Ymeasured=Ymodel+ϵσ

where σ is the standard deviation of the synthetic measurements and ∈ is a Gaussian variable with zero mean and unity standard deviation. The general physical properties of the two test cases are listed in Table 1, where it is also included the frequency of the measurements considered and the space between measurements. In this table, M is the number of measurements in space and N the number of measurements in time.

This table also presents the initial temperature, which was taken equal to T_∞, and the discrepancy between the two models (thin wall and full model), measured as the maximum and average difference for the measured temperature on r = r_e, from t = 0 up to t = t_f, considering A = 5,000 and p = 1.2 in equation (32). As it can be verified, the thin-wall model is suitable for both materials considered, which will also be shown next for the inverse analysis.

The differences between the estimated and the exact heat fluxes were calculated considering the following expression for the round mean squared (RMS) error:

(34) Eq=‖qestimated−qexact‖2AM˜⋅N˜

where M˜ and  N˜ are the number of measurements in space and time used. To avoid possible estimation problems close to the initial and final times, as well as from the upper and lower boundaries of the PHP, we discarded some parts of the spatial and temporal domains for the calculation of the RMS error. Therefore, M˜ and  N˜ represent the number of measurements after this cutting procedure was used.

The space window consisted of cutting 1 mm from each side of the domain. However, for the time window, to use a more physically appealing procedure, we chose to eliminate, at the beginning and at the end of the observation window, the time interval needed for an increase of ΔT = 0.1K, according to the equivalent lumped-capacitance model for a hypothetical heat flux of q = 10³W/m²:

(35) tcut≥ρLCpΔTq

Table 2 shows the time and space intervals used for the calculation of the RMS error for the two setups used in this paper.

Table 3 shows a comparison for the RMS error obtained using a Gaussian filter (Pagliarini et al., 2021a) and the present methodology, for a tube made of aluminum (Pagliarini et al., 2021a), for both the complete and the thin-wall model, considering adiabatic boundary conditions on z = 0 and z = L.

Because the estimates are based on synthetic data with random noises, these results are averaged values over 100 runs for the CG method with the thin-wall model and 10 runs for the CG method with the complete model. Results obtained by the filtering technique used only the thin-wall model and were averaged over 100 runs.

As a general trend, the current methodology produces results with much lower errors, compared with the filtering technique presented previously (Pagliarini et al., 2021a). Also, although the thermal conductivity of the aluminum is very high, the CG method based on the complete model still represents a little better the physical problem than the thin-wall one, therefore being able to better recover the unknown function. Nevertheless, the difference is not so large to justify the use of a complete model, which is computationally more expensive.

Figure 2 shows the exact and estimated function using the thin-wall model for a tube made of aluminum, with A = 4,000 W/m² and for four different values of p (i.e. 0.6, 0.8, 1.0, 1.2), considering σ = 0.1 K. These values have been chosen representative of the cases experimentally studied in Pagliarini et al.(2021a). The residuals of the temperatures are also shown, where one can see that they are uncorrelated and of the same order of magnitude of the standard deviation.

As a second synthetic set of data, we compared the presented methodology with the one presented in Cattani et al. (2019) for a sapphire tube (i.e. Tikhonov zero order regularization). Table 4 shows a comparison for the RMS error [equation (34)] obtained using a Tikhonov regularization (Cattani et al., 2019) and the present methodology, considering prescribed boundary conditions on z = 0 and z = L. Because the estimates are based on synthetic data with random noises, these results are averaged values over 100 runs for the CG method and 100 runs for the Tikhonov regularization. The results highlight that the here proposed methodology performs almost equivalently compared with the Tikhonov regularization technique adopted in Cattani et al. (2019).

In Figure 3, the exact and estimated heat fluxes using the thin-wall model for a tube made of sapphire, with p = 1.2 and for four different values of A (i.e. 2,000, 3,000, 4,000 and 5,000 W/m²), are shown, considering adiabatic walls on z = 0 and z = L. Also, in this case the values have been selected as representative of the real physical behavior in Cattani et al. (2019). The reported graphs underline the good performance of the current method in restoring the heat flux signal. Also, in this case, the random distribution of the temperature residuals confirms the appropriateness of the adopted solution technique.

Once the method was compared for synthetic data, we proceeded to some experimental analyses.

4.3 Experimental results

This section presents the application of the current methodology for a real case taken from Cattani et al. (2019) and Pagliarini et al. (2021a).

The test section experimentally investigated in Cattani et al. (2019) was represented by a single loop PHP comprised of an evaporator and a condenser section fabricated in copper (2 mm internal diameter and 1 mm of wall thickness), attached to two sapphire pipe (110 mm long) with the same diameters. The device was partly loaded (filling ratio = 60%) with pure ethanol. During the experiments, the input heat power at the evaporator was changed to obtain, in the sapphire portion, the conventional flow patterns normally detected in PHPs, that is, oscillating, slug/plug and semiannular flow.

Figure 4 presents a comparison between the present results and the one described in Cattani et al. (2019) for a tube made of sapphire for the three different flow conditions experimentally observed.

The alternating feature of the evaporative phenomena from the evaporator (red stripes in Figure 4 – oscillating working condition) and flow reversals from the condenser (blue lines in Figure 4 – oscillating working condition) are responsible for an incessant heat transfer among the fluid and the pipe surface and for the oscillating motions. The annular flow presents a continually rising vapor column from the evaporator to the condenser enveloped by a liquid film attached on the inner wall, whereas the start-up corresponds to the working condition that follows the activation of the device.

In Cattani et al. (2019), the selection of the proper regularization factor in the Tikhonov procedure was performed by implementing the L-curve technique defined by Hansen and O’Leary(1993).

The residual values between measured and estimated temperatures found adopting Tikhonov + L-curve method were used to end the iterative process in Cattani et al. (2019), adopting the discrepancy principle. For the CG with AO, there is no need to use the Tikhonov regularization and the iterative process was stopped using the Morozov’s discrepancy principle.

The color scale was adjusted to match the maximum and minimum value of each figure. As it can be seen, both results are qualitatively close to each other, although some small discrepancies on the numerical values can be found. Nevertheless, because this is a real test case, an exact value of the heat flux is not established. One advantage of the current approach is that it can handle large amounts of data at once, whereas the one used in Cattani et al. (2019) had to perform the estimate dividing the total time into small fractions of 2 s each because of the significantly high computational cost to build the sensitivity matrix.

Finally, Figure 5 presents a comparison among the present outcomes and the one presented in Pagliarini et al. (2021a) for a tube made of aluminum, where the same physical behavior was observed.

In Pagliarini et al. (2021a), a PHP consisting in an aluminum tube (inner/outer diameter = 3/5 mm) closed in a 14 turns loop was tested in microgravity conditions. The analysis was able to detect two flow regimes, that is, intermittent flow and full activation. In Figure 5, it is reported a case of intermittent flow where sporadic fluid motion is observed in only parts of the 14 turns. The reported heat flux refers to one of the 14 branches and it was computed using the Gaussian filtering technique. Under this approach the cutoff frequency was selected adopting discrepancy principle considering a standard deviation σ = 0.06 K that was experimentally measured in Pagliarini et al. (2021a).

Despite discrepancies at the edges of the two considered maps, probably because of the highly smoothing effect of the filtering method, both approaches provide qualitative similar results for the estimated heat flux.

The main contribution of this paper, besides comparing different methodologies to estimate heat fluxes on PHPs is to present a fast methodology that is not restricted, at least for the practical applications shown here, by the amount of data used and by the thermal properties of the wall.

Also, supplementary material containing all data used in the figures is provided and readily employable by other authors for the evaluation of the effectiveness of their methodologies. The data are sorted as follows:

• Fig2-p0.6-estimated.dat: estimated values appearing on Figure 2, for A = 4,000 and p = 0.6.

• Fig2-p0.8-estimated.dat: estimated values appearing on Figure 2, for A = 4,000 and p = 0.8.

• Fig2-p1.0-estimated.dat: estimated values appearing on Figure 2, for A = 4,000 and p = 1.0.

• Fig2-p1.2-estimated.dat: estimated values appearing on Figure 2, for A = 4,000 and p = 1.2.

• Fig2-p0.6-exact.dat: exact values appearing on Figure 2, for A = 4,000 and p = 0.6.

• Fig2-p0.8-exact.dat: exact values appearing on Figure 2, for A = 4,000 and p = 0.8.

• Fig2-p1.0-exact.dat: exact values appearing on Figure 2, for A = 4,000 and p = 1.0.

• Fig2-p1.2-exact.dat: exact values appearing on Figure 2, for A = 4,000 and p = 1.2.

• Fig3-A2000-estimated-adiabatic-thin-wall-model.dat: estimated values appearing on Figure 2, for A = 2,000 and p = 1.2.

• Fig3-A3000-estimated-adiabatic-thin-wall-model.dat: estimated values appearing on Figure 2, for A = 3,000 and p = 1.2.

• Fig3-A4000-estimated-adiabatic-thin-wall-model.dat: estimated values appearing on Figure 2, for A = 4,000 and p = 1.2.

• Fig3-A5000-estimated-adiabatic-thin-wall-model.dat: estimated values appearing on Figure 2, for A = 5,000 and p = 1.2.

• Fig3-A2000-exact.dat: exact values appearing on Figure 2, for A = 2,000 and p = 1.2.

• Fig3-A3000-exact.dat: exact values appearing on Figure 2, for A = 3,000 and p = 1.2.

• Fig3-A4000-exact.dat: exact values appearing on Figure 2, for A = 4,000 and p = 1.2.

• Fig3-A5000-exact.dat: exact values appearing on Figure 2, for A = 5,000 and p = 1.2.

• Fig4-Annular-estimated.dat: estimated values appearing on Figure 4 for the CG method, considering the annular working condition.

• Fig4-Oscillating-estimated.dat: estimated values appearing on Figure 4 for the CG method, considering the oscillating working condition.

• Fig4-Startup-estimated.dat: estimated values appearing on Figure 4 for the CG method, considering the startup working condition.

• Fig4-Annular-meas.dat: experimentally measured temperatures and heat fluxes by estimated by Cattani et al. (2019), considering the annular working condition.

• Fig4-Oscillating-meas.dat: experimentally measured temperatures and heat fluxes by estimated by Cattani et al. (2019), considering the oscillating working condition.

• Fig4-Startup-meas.dat: experimentally measured temperatures and heat fluxes by estimated by Cattani et al. (2019), considering the startup working condition.

• Fig5-estimated.dat: estimated values appearing on Fig. t for the CG method.

• Fig5-meas.dat: experimentally measured temperatures and heat fluxes by estimated by Pagliarini et al. (2021a).

All data files can be read according to the following pseudo-code, where dimensions are seconds for time, meters for distances, Kelvin for temperatures and Watts per squared meters for heat fluxes:

read M, N (M is the number of points in space, N is the number of points in time)do j = 1, N   do i = 1, M      read time(j), space(i), temperature(i,j), heat flux(i,j)   end doend do

6. Conclusions

The aim of the current research was to propose and validate an estimation method to assess the local convective heat flux on the inner wall of a duct, under a convection problem in transient conditions. The approach was established on the solution of the IHCP in the wall by using as starting data the temperature measurements on the outer surface. The procedure is based on the CG method with AO. The here proposed approach was first verified adopting synthetic data and then it was validated with real cases regarding PHPs taken from Cattani et al.(2019) and Pagliarini et al.(2021a). These study cases are especially appropriate for assessing the proposed solution method because they present a local heat flux distribution characterized by an important variation both in time and in space as a result of the complex flow patterns that are generated in this type of devices. The experimental analysis allowed the comparison with other two recognized and robust approaches: Tikhonov regularization method and Gaussian filter. The comparative highlighted an estimating performance improved respect to Gaussian filter and comparable with Tikhonov.

Eventually, it must be stressed that even if the here proposed method and the Tikhonov regularization method present almost equivalent performance, the procedure here proposed shows some benefits as it affords a general model of the heat conduction problem that is effortlessly customized for the particular case, it is not restricted by the amount of data used and it presents reduced computational expense.

Figures

Figure 1.

Sketch of the test section

Figure 2.

Exact and estimated heat fluxes considering the thin-wall model with A = 4,000 W/m² for a tube made of aluminum

Figure 3.

Exact and estimated heat fluxes considering the thin-wall model for a tube made of sapphire, with adiabatic boundary conditions on z = 0 and z = L

Figure 4.

Comparison between the current results and the ones presented in Cattani et al. (2019) for different working conditions of a pulsating heat pipe made of sapphire

Figure 5.

Comparison between the current results and the ones presented in Pagliarini et al. (2021b) for a pulsating heat pipe made of aluminum

Table 1.

Physical properties and general constants for the direct problem

Variable	Aluminum tube (Pagliarini et al., 2021a)	Sapphire tube (Cattani et al., 2019)
L [m]	0.1	0.04
tf [s]	20	2
ri [m]	0.0015	0.001
re [m]	0.0025	0.002
k [W/m · K]	201	27
h [W/m2 · K]	10	5
T∞ [K]	295.15	294.18
ρ [kg/m3]	2,700	3,970
cp [J/kg · K]	900	419
Δzmeas [m]	0.001	0.001
Δtmeas [s]	0.02	0.02
M	101	41
N	1,001	101
Maximum discrepancy between the thin-wall and the full model for T on r = re from t = 0 to tf [K]	0.025	0.034
Average discrepancy between the thin-wall and the full model for T on r = re from t = 0 to tf [K]	0.008	0.017

Table 2.

Windows considered for the calculation of RMS error

Cutting window	Aluminum tube	Sapphire tube
Space window	[0.001, 0.099] m	[0.001, 0.039] m
Time window	[0.26, 19.74] s	[0.18, 1.82] s

Table 3.

E_q [%] for the results obtained with the CG method with adjoint operator and the filtering technique for a tube made of aluminum

A	p	Method	σ = 0.01 K	σ = 0.06 K	σ = 0.1 K
2,000	1.2	CG method (thin-wall model)	5.15	11.25	17.34
				CG method (complete model)	3.53	10.37	16.65
				Filtering (Pagliarini et al., 2021a)	14.35	20.32	24.28
		1.0	CG method (thin-wall model)	4.56	9.10	13.39
				CG method (complete model)	3.33	8.33	12.74
				Filtering (Pagliarini et al., 2021a)	12.10	17.02	20.22
		0.8	CG method (thin-wall model)	4.17	7.52	10.39
				CG method (complete model)	3.20	6.80	9.64
				Filtering (Pagliarini et al., 2021a)	9.88	13.87	16.44
		0.6	CG method (thin-wall model)	3.88	6.10	8.17
				CG method (complete model)	3.06	5.74	7.66
				Filtering (Pagliarini et al., 2021a)	7.63	10.74	12.72
4,000	1.2	CG method (thin-wall model)	9.44	14.76	19.70
				CG method (complete model)	5.54	12.29	17.63
				Filtering (Pagliarini et al., 2021a)	13.44	16.97	19.28
		1.0	CG method (thin-wall model)	8.66	12.71	16.35
				CG method (complete model)	5.25	10.49	14.73
				Filtering (Pagliarini et al., 2021a)	11.35	14.30	16.23
		0.8	CG method (thin-wall model)	7.70	10.83	13.51
				CG method (complete model)	5.07	9.17	12.20
				Filtering (Pagliarini et al., 2021a)	9.24	11.65	13.18
		0.6	CG method (thin-wall model)	6.75	9.24	11.08
				CG method (complete model)	5.03	8.77	10.68
				Filtering (Pagliarini et al., 2021a)	7.11	9.04	10.29

Table 4.

E_q [%] for the results obtained with the CG method with adjoint operator and the Tikhonov regularization for a tube made of sapphire, considering prescribed temperature boundary conditions

A	p	Method	σ = 0.1 K
			Prescribed temperature at z = 0 and z = L
2,000	1.2	CG method	17.77
		Tikhonov (Cattani et al., 2019)	19.40
3,000		CG method	14.99
		Tikhonov (Cattani et al., 2019)	15.30
4,000		CG method	13.90
		Tikhonov (Cattani et al., 2019)	13.60
5,000		CG method	13.16
		Tikhonov (Cattani et al., 2019)	12.80

Supplementary material

The supplementary material for this article can be found at: http://colaco.org/share/IJNMHFF2022/supplemental.zip