Inverse geometry problem of estimating the phase front motion of ice in a thermal storage system

Abstract

The moving phase change interface in the latent heat energy storage system is estimated by applying a two-dimensional inverse geometry problem. The energy storage system of the vertical tube type is considered for the present inverse geometry problem. To solve the phase change problem, the boundary element method is adopted. The moving phase change interface is estimated by using the conjugate gradient method. Estimation of the phase front motion is verified by conducting the inverse analysis for an assumed phase front motion. An inverse analysis for the desired temperature distributions is executed to investigate the possibility of desired front motion control or ice monitoring. The effects of the noise levels and the thermocouple spacing on the inverse solutions are also examined.

Keywords:

• Inverse geometry problem
• Phase front motion
• Thermocouple spacing

Introduction

By using inverse analysis, engineers can obtain the thermal quantities that are difficult to measure or analyze. When the shape of the problem is unknown, the inverse geometry problem can be applied to the shape estimation.

Thermal energy storage is becoming popular as an energy saving method to solve the imbalance between day and night energy demand patterns. In particular, the latent heat energy storage system has advantages of low volume/energy ratio together with its small temperature deviation as compared with the sensible heat storage system.

In the previous work on the thermal storage system by Ismail et al. [1–4], a direct phase change problem on the solidification of phase change material was studied numerically.

Zabaras et al. [5] estimated the thermal boundary conditions for desired phase front motion. Huang and Chen [6] carried out the inverse geometry problem successfully by using the boundary element method and conjugate gradient method. Huang and Hsiung [7] also conducted the optimal design problem of the cooling passage inside a turbine blade by using an inverse geometry problem.

The objective of this work is to extend the inverse geometry problem to the phase change problem in estimating the ice growth process in a vertical tube type energy storage system.

Direct Problem

Two-dimensional Phase Change Problem

In the thermal storage system shown in Fig. 1, the liquid water is contained inside the external shell, while the working fluid flows into a set of tubes inside the shell. When the temperature of the working fluid is lowered below the phase change temperature, the solid–liquid interface is growing on the outside surface of the inner tube.

FIGURE 1 Thermal storage unit and tube arrangement inside a storage tank.

The mathematical modeling is conducted using only one inner tube of the thermal storage system as shown in Fig. 2. The thermal resistance of the inner tube is assumed to be negligible. To transform the formulation to the dimensionless form, the following dimensionless quantities are defined.

FIGURE 2 Solid and liquid domain with boundary shape.

where the ‘’ denotes dimensional quantities, means the latent heat, x, y are dimensionless geometrical coordinates, and α is dimensionless thermal diffusivity. are reference length and initial temperature, respectively. Bi is the Biot number that means heat transfer between tube surface and working fluid, and Ste is the Stefan number that means the effect of latent heat. In the formulation, it is assumed that the phase change is dominated by heat conduction. Considering the above assumptions, the dimensionless governing equation for a two-dimensional phase change problem is given as: (--1-) (--1-) (--1-) (--1-) (--1-) (--1-) (--1-) where n is normal outward direction to the boundary, v_n is the n-direction velocity of moving interface, and T_m is the phase change temperature.

Boundary Element Method

Since the BEM [8] can readily accommodate the change of unknown boundary, the BEM is adopted for the estimation of interface shape.

After multiplying the conduction equation by the Green function, integrating the result with respect to all space and time domain, and applying Green’s second identity, the following equation is obtained. (--2) (--3) (--4) where K is the thermal diffusivity, t_J denotes the final time step, and N is a point to be analyzed. T * is Green function, and q * is the derivative of T * with respect to the normal direction; r is the length to a point on the boundary from N, and τ = t_J − t; d is the partial derivative of r with respect to normal direction.

Since the initial temperature is zero, the last integral term of the right-hand side of Eq. (2(--2) ) is eliminated. Reformulating the result, the boundary element equation is obtained. By solving the boundary element equation, the unknown temperatures or heat fluxes can be obtained.

Inverse Problem

The present inverse estimation of the phase change interface is to be performed in such a way that the following functional is minimized. (--5) Here T_S,m(t) is the calculated temperatures by the estimated interface, and Y_m(t) is the measured temperatures.

Conjugate Gradient Method

The CGM iterative process is used for the estimation by minimizing the functional. The phase change interface is given as (--6-) (--6-) (--6-) where ϕ is the angle between the horizontal and outward normal direction of the unknown boundary. is the direction of descent (i.e., search direction) given by (--7) where γⁿ is the conjugate coefficient defined as (--8)

Sensitivity Problem and Search Step Size

The sensitivity problem is obtained by replacing Γ₂ by Γ₂ +ΔΓ₂ and T by T+ΔT in the direct problem, then subtracting the direct problem from the resultant expression, and neglecting the second-order terms.

(--9-) (--9-) (--9-) (--9-) (--9-)

Only one of the two boundary conditions is applied to unknown interface, since the measured temperatures are used additionally. If Eq. (9d(--9-) ) is used for the boundary condition, the iteration cannot be performed because of the fluctuation resulting from the condition including a ‘v_n’ term coupled with an unknown interface that is changed with iterations. Therefore, Eq. (9e(--9-) ) is used for the interface condition in this study.

The functional can be rewritten as follows:

(--10)

The following expression on the search step size is obtained by differentiating Eq. (10(--10) ) with respect to β_n and letting by zero. (--11) The sensitivity function can be obtained from the solutions of Eq. (9(--9-) ) by letting .

Adjoint Problem and Gradient Equation

To obtain the adjoint problem, Eq. (1(--1-) ) is multiplied by the , then integrated over all space and time domain, and added the functional. In the resulting equation, the ΔJ is obtained by perturbing Γ₂ by ΔΓ ₂ and T by ΔT, then subtracting the original equation, and neglecting the second-order terms. This gives (--12) where is the Dirac delta function. In Eq. (12(--12) ), the domain integral term is reformulated based on Green’s second identity, the boundary and the initial conditions of the sensitivity problem are utilized, and the integral including ΔT becomes zero. Then the following adjoint problem is obtained (--13-) (--13-) (--13-) (--13-) After the introduction of the adjoint problem, the following integral term is left. (--14) By definition [6], the functional increment can be represented as

(--15)

By comparing Eq. (14(--14) ) with (15(--15) ), the gradient of functional is obtained as

(--16)

In Eq. (16(--16) ), the estimation cannot be performed since the initial and final values are always equal to zero. Therefore the following assumptions are introduced. (--17) (--18) where Δt is the numerical time increment.

Stopping Criterion

If the problem contains no measurement errors, the traditional check condition is specified as (--19) where ε is a small value. However, the measured temperature data may contain measurement errors. Therefore the stopping criteria ε is obtained by using the discrepancy principle [9]. (--20) where σ is the standard deviation of the measurements, which is assumed to be a constant.

Computational Procedure

When is available, the calculation for the iteration n + 1 is summarized as below.

Table

Step 1	Solving the direct problem given by Eq. (1), and compute the functional of Eq. (5)
Step 2	Examine the stopping criteria given by the Eq. (19) and (20).
Step 3	Solve the adjoint problem of Eq. (13).
Step 4	Compute by Eq. (16).
Step 5	Compute and from Eq. (7) and (8), respectively.
Step 6	Solve the sensitivity problem of Eq. (9).
Step 7	Compute the search step size of Eq. (11).
Step 8	Estimate the new from Eq. (6), and return to Step 1.

Results

We consider one specific inverse problem whose phase change interface is assumed to be the known function to verify the methodology considered in the present work. The specified moving velocity of the interface and the relational equation between the moving velocity and the phase change interface are used for assuming the known function of the interface.

The temperature readings for this inverse problem are produced numerically by adding numerical noises to the exact solution of the direct problem. The numerically produced measurement data can be expressed as follows. (--21) where Y_exact is the exact solution of the direct problem, σ is the standard deviation of the measurements, and ω is a random variable generated by IMSL subroutine DRNNOR; ω has a value within −2.576 and 2.576 for the 99% confidence bounds.

The reference length representing the radius of inner tube is given as unity. The thermal properties of ice in the solid and liquid regions are given as:

In order to specify an arbitrary phase change configuration, the moving velocity of a phase change interface is assumed to be (--22)

The explicit relationships between the moving velocity and the phase change interface are given as

(--23-) (--23-)

From the interface velocity given by Eq. (22(--22) ) and the explicit relationships given by Eq. (23(--23-) ), the phase change interface for the validation is specified. The interface configuration specified by Eqs. (22(--22) ) and (23(--23-) ) is shown in Fig. 3.

FIGURE 3 Exact interface shape of ice.

The initial guessed interface configuration for the iteration is assumed to be the surface shape of inner tube, .

The difference between the exact and the estimated interface is expressed by the following shape error. (--24) where I and J represent the number of elements discretized by space and time, respectively.

The relative quantity of measurement errors is calculated by the following measurement error. (--25) Here Y_avg is the average of measured data.

Validation Case Using 9-Thermocouples

The inverse analysis using 9-thermocouples, denoted by the dots in Fig. 2, is performed. Without measurement errors (σ = 0), the given stopping criteria ε = 0.002 is satisfied after the estimation with 33-iterations. The estimated result of the phase change interface is shown in Fig. 4. In this case, the functional is converged to 0.0018 from 1900, the shape error to 0.02% from 11.0%. From Figs. 3 and 4, it can be seen that the phase change interface of ice is obtained accurately in the case of no measurement errors.

FIGURE 4 Estimated interface shape of ice using σ = 0.0 and M = 9.

Since measurement errors are always introduced in any real measurement, it is required to observe the influence of measurement errors on the inverse estimation. Figure 5 shows the estimated interface after 13-iterations in the case of σ = 0.05. In this case, the measurement error is calculated as 1.40%. The functional is converged to 0.37 from 1900, and the shape error to 0.27% from 11.0%. In Fig. 6, the estimated interface converged after 8-iterations is shown in the case of σ = 0.2. Here the measurement error is calculated as 5.70%. The functional is converged to 2.2 from 1900, and the shape error on 0.75% from 11.0%. The above results indicate that the phase change interface estimation is accurately performed even though measured temperatures contain measurement errors.

FIGURE 5 Estimated interface shape of ice using σ = 0.05 and M = 9.

FIGURE 6 Estimated interface shape of ice using σ = 0.2 and M = 9.

Validation Case Using 5-Thermocouples

In most real systems, there are some difficulties in installing enough sensors. Therefore, it is necessary to verify the effect of having a small number of sensors. To examine the influence of the number of sensors, the inverse analysis is performed using 5-thermocouples.

In the case with 5-thermocouples, only five-spatial information can be obtained. But the five-spatial information is not enough to express the interface for the inverse calculation. Therefore, the phase change interface is expressed more closely by interpolation between each thermocouple.

Figure 7 shows the estimated result of phase change interface in the case with no measurement error. It satisfies the given stopping criteria ε = 0.002 after 33-iterations. The functional is converged to 0.0018 from 1900, the shape error to 0.14% from 11.0%. Using 5-thermocouples without measurement errors, the phase change interface of ice is estimated accurately.

FIGURE 7 Estimated interface shape of ice using σ = 0.0 and M = 5.

In Fig. 8, the estimated interface after 18-iterations is shown in the case of σ = 0.05. The measurement error is 1.40%, the functional is converged to 0.19 from 1100, and the shape error to 0.35% from 11.0%. As expected, the estimated interface with measurement errors is worse than that with no measurement error. But the estimated interface shows a good agreement with the exact interface.

FIGURE 8 Estimated interface shape of ice using σ = 0.05 and M = 5.

In Fig. 9, the estimated interface after 12-iterations is shown in the case of σ = 0.2. Here the measurement error is 5.80%, the functional is converged to 3.4 from 1100, and the shape error to 1.1% from 11.0%. Though there are a small number of thermocouple and a large measurement error, the interface can also be estimated with good agreement.

FIGURE 9 Estimated interface shape of ice using σ = 0.2 and M = 5.

Results on the Desired Temperature Distributions

To illustrate the ability of the present inverse estimation in operating the practical thermal storage system and monitoring the generated ice, we consider arbitrary temperature distributions on the thermocouple location. The estimations for the arbitrary temperature distributions are performed.

Figure 10(a) shows the estimated interface configuration when imposing relatively higher temperature distribution on the first thermocouple and lower temperature distributions on the fifth and sixth. It can be seen that a thinner ice layer is generated around the first thermocouple, and a thicker ice layer around the fifth and sixth. Figure 10(b) is the estimated result for the similar temperature distribution as in Fig. 10(a), but temperature increases more with time compared to that of Fig. 10(a). Figure 10(b) shows that the ice growth is more active than that of Fig. 10(a) because of increasing temperature differences.

FIGURE 10 Estimated ice shape for higher desired temperatures at 1 and lower at 5, 6: (a) lower temperature increment; (b) higher temperature increment.

Figure 11(a) illustrates the estimated interface configuration when imposing relatively lower temperature on the ninth point and increasing the temperature clockwise from that point. It can be seen that the thick ice layer is generated around the ninth thermocouple, and the thickness of ice layer decreases clockwise from that location. Figure 11(b) shows the estimated result with the similar temperature distributions as in Fig. 10(a), but temperature increases more with time compared to that of Fig. 11(a). It is observed that the thickness of ice layer in Fig. 11(b) is thicker than that of Fig. 11(a).

FIGURE 11 Estimated ice shape for higher desired temperatures along 1, 2, 3, medium along 4, 5, 6 and lower along 7, 8, 9: (a) lower temperature increment; (b) higher temperature increment.

Conclusion

In this work, the phase change interface is estimated using the inverse geometry technique by applying the conjugate gradient method (CGM) along with the boundary element method (BEM). Several cases for different numbers of sensors and measurement errors are successfully performed. Physical validity for an arbitrary temperature distribution imposed on the thermocouples is also examined. The present results show that the phase change interface can be estimated with a good accuracy by using the method considered in this work. The present work may also be used for controlling and monitoring the thermal storage system.

Nomenclature

Table

	= constant-pressure specific heat
J	= functional defined by Eq. (5)
J′	= gradient of functional
k	= thermal conductivity
La	= latent heat
M	= total number of thermocouples
T	= dimensionless temperature
	= sensitivity function
Ste	= Stefan number
Y	= dimensionless measured temperature

Greek Symbols

Table

β	= search step size
γ	= conjugate coefficient
Ω	= computational domain
Γ	= boundary of computational domain
λ	= Lagrange multiplier of adjoint problem
ω	= random variable in measured data
σ	= standard deviation of the measurement errors

Subscripts

Table

L	= liquid region
m	= location of thermocouple
S	= solid region
∞	= working fluid

Superscripts

Table

⁁	= estimated variables
n	= number of iteration
*	= Green function

Acknowledgement

This work was supported by the research fund of Hanyang University (HY-2001).