Approximating the solution of the discharging process in a domestic hot water storage tank

ABSTRACT In this article a new approximation procedure for the temperature during the discharging process in a domestic hot water storage tank is developed. The main idea of this procedure is the mathematical description of the dynamic behaviour of a boundary layer that appears inside the storage tank. This leads to an approximate solution of a simple structure but a high accuracy throughout the whole discharging process for the temperature distribution inside the storage tank. For validation the approximate solution is compared with the exact solution that is constructed with help of generalized Fourier coefficients and a numerical solution that is obtained by a semi-discretization method. Adapting the method leads also to an approximation procedure for charging processes.


Introduction
In this article a domestic hot water storage tank is considered. Such storage tanks are typically used in one family houses for supplying people with hot water. The advantage of such storages is that they can also be used for buffing thermal energy if a solar plant or a micro combined heat and power plant is installed. So there is a great interest to simulate the dynamic thermal behaviour of such a storage tank to test control strategies for such power plants for example.
There are many tools for the numerical simulation of the thermal behaviour of such a storage in three space dimensions. Such a simulation can give information of the distribution of the temperature inside the storage and can be used in the process of constructing the geometry of the storage. Because of the high numerical costs of such a simulation these three dimensional simulations can be used only for short time intervals.
For long time simulations one dimensional models in space are more appropriate. In such a model it is assumed that the temperature inside the storage depends only on one space variable. This is a simplifying assumption so that certain phenomena cannot be described by such a model as for example colder water at the wall and warmer water in the middle of the storage tank. There are various approaches for the numerical treatment of such a one dimensional model. Kreuzinger, Bitzer and Marquardt [1] consider a stratified storage tank and use a method of lines approach [2]. For the spatial discretization a finite volume up-and downwind scheme depending on the flow direction as well as a central difference scheme is used. Because of the occurrence of oscillations for sharp spatial transitions, furthermore high resolution slope limiter schemes are considered to avoid these oscillations. Numerical tests with various limiter functions yield the best results using a Superbee limiter with a certain limiter function [3]. The resulting ODE system is treated with the ode23tb solver, a solver for stiff ODEs which uses the trapezoidal rule in combination with a backward difference scheme.
Mojtabi and Deville [4] use a finite element approach for the spatial discretization of a linear advection-diffusion equation with homogeneous Dirichlet boundary conditions and a sine profile for the initial condition. The resulting ODE system is separated in a viscous part and an advection term. For the viscous part an implicit Crank-Nicolson time scheme (combination of forward and backward differences) is used and the advection term is treated with a second order Adams-Bashforth scheme.
Nash, Badithela and Jain [5] consider a sensible thermal energy storage tank with an immersed coil heat exchanger. The tank is discretized vertically and in every finite volume the conservation of energy is used to derive an ODE system (finite volume method). The ODE system is treated by a finite difference scheme with a temperature inversion correction method [6].
Walter, Strohmayer and Hameter [7] also use a finite volume method for the numerical treatment of a thermocline energy storage device.
An alternative to the numerical treatment is the usage of analytical solutions if such solutions are available and suitable for use. The advantages and disadvantages of theoretical calculations are already discussed in [8]. If an analytical solution is available it can be used for validation or comparison of numerical results. It can also be used for the validation of analytical approximations of the exact solution if the exact solution is not suitable for practical use. The usage of simple analytic approximations has the advantage of low numerical cost that is important especially for long time simulations.
There are approaches for approximating the temperature profile in a storage tank with a smooth cumulative distribution function of a probability distribution as for example the logistic distribution function [9][10][11] or the Fermi-Dirac distribution function [12]. There are also approaches with sets of splines or polynomials [13] and with sets of Fourier polynomials [14].
The storage tank that is considered here will be described in space by an ideal circle cylinder with height L and radius R. Figure 1 provides a sketch of the storage tank.
The temperature distribution inside the storage is assumed to be only dependent on the spatial variable x, so it is a one dimensional model in space. Describing the thermal behaviour of the storage there are essentially three operational modes that have to be considered: discharging, charging and stand-by.
Discharging takes place if hot water is requested by tapping. During a discharging process cold water is filled in at the bottom of the storage (position 1) and the requested hot water is taken from the top of the storage (position 2). In a charging process hot water with a desired set temperature is filled in on top (position 3) and cold water is taken from the bottom of the storage (position 4). The hot water for the charging process can be provided from a primary heat source (a boiler for example) or a secondary heat source that can be a solar plant or a micro combined heat and power plant.
The scientific contribution of this article is the development of a new approximation procedure for the discharging process that provides an approximate solution for the temperature profile in the storage tank in terms of simple functions. This new procedure uses a mathematical description for the growing process of the thickness of the boundary layer that appears during the discharging process. So, a high accuracy of the approximate solution can be achieved every time in the process, although the approximate solution has a quite simple structure. This combination of simplicity and accuracy can be beneficial in long time simulations as for example for testing new control strategies.
The article is structured as follows. In Section 2 an initial and boundary value problem for the temperature inside the storage during the discharging process is presented. For this problem uniqueness and existence results are provided and the unique solution is constructed with help of generalized Fourier coefficients.
Section 3 first provides a detailed discussion of the numerical problems that appear by evaluating the exact solution for larger values of the velocity in the storage. An explicit condition for the number of available digits is derived to avoid instability of the exact solution. These considerations give a motivation for constructing an approximate solution. In a first step the non-diffusion case leads to a non-continuous solution that can be described with help of the Heaviside function. Because of the temperature difference of the water in the storage tank and the incoming cold water, diffusion effects lead to a smoothening of the jump in the solution of the non-diffusion case and hence to the appearance of a boundary layer. The thickness of this boundary layer is growing in time up to a stationary thickness. The main idea of the approximation procedure is the mathematical description of the growing process of the boundary layer. For this the Heaviside function in the non-continuous solution is approximated by a suitable function. This approximation of the Heaviside function uses a mathematical expression that describes the growing process of the thickness of the boundary layer and a suitable chosen sigmoid profile inside the boundary layer. For describing the growing process, three parameters are used that can be interpreted as an initial thickness, a stationary thickness and a time constant for the process. With help of the approximate Heaviside function and the initial profile of the temperature now the approximate solution for the temperature inside the storage tank during the discharging process can be constructed. For the determination of the three parameters two alternative methods are presented. The first method uses a minimization of a cost function and the second method uses a matching of the time derivative of the approximate solution to the time derivative of the exact solution in suitable chosen points.
In Section 4 numerical examples are used for a validation of the approximation procedure. For this on the one hand the approximate solution is compared with the exact solution and on the other hand with a numerical solution that is obtained by a semi-discretization method. Although the main focus of this article lies on the discharging process, in Section 5 the operational modes charging and stand-by are discussed shortly. The results obtained for the discharging process are used to adapt the method to charging processes. Furthermore a simple description for stand-by periods is given. Finally, Section 6 gives a short conclusion.

Discharging process
In this Section an initial and boundary value problem for the temperature inside the storage tank during a discharging process is presented and the unique solution of the problem is constructed.
Starting at time t 0 water with the constant temperature # S;in is flowing into the storage at position x ¼ 0. Considering heat conduction and convection, the energy balance leads to a convection diffusion equation for the temperature # S in the storage. Given a continuous initial temperature profile # 0 : ½0; L� ! R and a time T > t 0 , the discharging process can be described by the initial and boundary value problem [1]: Table 1 collects the used nomenclature in (1) -(4). For constructing an approximate solution of the problem (1) -(4), it has to be shown that there exists a unique solution that can be approximated. The following Theorem 1 provides the uniqueness of a solution in the class of classical differentiable functions. We set: Furthermore let C 2 ðU T Þ be the class of two times continuously differentiable functions on U T and CðU T Þ the continuous functions on U T . Here U T denotes the closure of the set U T . Hence the class C 2 ðU T Þ \ CðU T Þ consists of all functions that are continuous on U T and two times continuously differentiable in U T .

Theorem 1. The initial and boundary value problem (1) -(4) has at most one solution in the class C
Proof. Let # S and # S be two solutions in the class C 2 ðU T Þ \ CðU T Þ. Setting: we get: D#ðx; tÞ ¼ D# S ðx; tÞ À D# S ðx; tÞ ¼ 0 and hence: The boundary conditions lead to: Furthermore the initial condition yields for all x 2 U: Multiplying (8) with #ðx; tÞ and integrating over U T leads to: #ðx; tÞ @ @t #ðx; tÞ dx dt |ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl {zffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl } #ðx; tÞ @ @x #ðx; tÞ dx dt |ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl {zffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl } #ðx; tÞ @ 2 @x 2 #ðx; tÞdx dt |ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl {zffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl } ¼: I 3 (11) and hence: Evaluating the integrals I 1 ; I 2 ; I 3 we will get the boundary condition #ðL; �Þ ¼ 0. Under the regularity assumptions we can conclude: Here k �k L 2 ðUÞ denotes the L 2 -norm [15]. Furthermore we get: and: Equation (12) now yields: Since all the terms on the right hand side are nonnegative and # is continuous, we can conclude: The continuous function # takes its maximum and its minimum value on the compact set U T . The operator D is a parabolic operator with vanishing absolute term. From the strong maximum principle and the strong minimum principle for parabolic operators [15] it can be concluded that # takes its maximum and its minimum value only on the boundary and: For every k 2 N then λ k is an eigenvalue for the Sturm Liouville eigenvalue problem: and by an eigenfunction φ k to the eigenvalue λ k is given. Furthermore every continuous function g : ½0; L� ! R can be represented in the form: with the generalized Fourier coefficients: Proof. Please see for example [16]. □ Theorem 3. Let ðλ k Þ be the sequence of eigenvalues from Lemma 2. Furthermore let # 0 be a continuous initial profile. Then the function # S defined by solves the initial and boundary value problem (1) - (4).
Proof. Since # 0 is continuous, the function g defined by is continuous and hence Lemma 2 yields the representation: with c k as in (27). Now we consider the function u defined by: Because of λ k ! 1 for k ! 1 with OðkÞ the factor e À a W λ k ðtÀ t 0 Þ is rapidly decaying for every t > t 0 . So the term under the sum is a C 1 function [17] Since u is a solution to (31), # S is a solution to (1). The initial condition (32) and the boundary conditions (33) and (34) lead to: So # S defined by (26) solves the initial and boundary value problem (1) -(4). □

Approximate solution
From Section 2 we know that the initial and boundary value problem (1) -(4) has a unique solution in the class C 2 ðU T Þ \ CðU T Þ and with (26) a representation of the solution is given. Representing (26) in the form: it can be seen that for and large values of the velocity v D the term on the left hand side of (40) is getting very small and hence the limit value of the infinite sum on the right hand side. Let the largest term in the infinite sum be of order 10 À N 1 and the machine accuracy be of order 10 À N 2 . Then the smallest value that can be computed numerically is of order 10 À ðN 1 þN 2 Þ . If now the exponential term on the left hand side in (40) is of order 10 À ðN 1 þN 2 þN 3 Þ , the computed value for # S ðx; tÞ À # S;in is of order 10 N 3 . If for example N 3 ¼ 10 the computed value for # S ðx; tÞ À # S;in is of order 10 10 while the exact value is of order 10 1 . This means, if the machine accuracy is not enough for the considered values of the velocity, the just described process leads to an instable behaviour.
To get a condition, for N 2 N and x ¼ L we consider: The number of available digits then must be at least N plus the required number of exact digits. Here it has to be taken into account that the computation of the (finite) sum requires a certain number of exact digits. For 16 available digits (double precision) and 4 required exact digits, condition (42) leads to: Typical discharging velocities are in the region of 1m=h. With the value for a W from Table 1 and a storage height of L ¼ 1m we get: Condition (42) leads to N > 434 and the number of required exact digits has to be added for computing the number of available digits. Mojtabi and Deville [4] investigate the same problem by treating an advection diffusion equation and take 500 digits for a realistic Reynolds number of 2000.
The previous considerations show that the representation (26) is not suitable for the numerical computation with large velocities and hence for an application in a simulation. Therefore we will derive an approximation of the exact solution which is more appropriate for numerical simulations. Nevertheless, representation (26) will serve for the validation of the approximation.
During discharging processes with large velocities the convective power is dominant so that heat conduction does not play an important role concerning the energy balance. Neglecting heat conduction in (1), that means setting a W ¼ 0, Equation (1) turns to a linear transport equation [7] with the solution: Then the initial profile is transported through the storage with the velocity v D . If # 0 ð0Þ�# S;in the jump at x ¼ v D ðt À t 0 Þ is also transported through the storage. Since the solution of the transport equation does not consider heat conduction this solution represents for # 0 ð0Þ�# S;in a non-continuous solution with a jump at x ¼ v D ðt À t 0 Þ. The solution (45) can be written in the form: with the Heaviside function H. Caused by heat conduction a diffusion process leads to a smoothening of the jump and so there appears a boundary layer that is growing up in time starting at the initial time t 0 to a (quasi) stationary boundary layer. To describe the growing process of the thickness of the boundary layer we consider ε 0 ; ε 1 ; τ > 0 with ε 0 < ε 1 and define: Here the value ε 0 can be interpreted as the thickness of the boundary layer at the starting time t 0 , the value ε 1 describes the stationary thickness and τ is a time constant. The value εðtÞ describes the thickness of the boundary layer at time t � t 0 . For numerical reasons, the value ε 0 cannot be chosen arbitrary small. With help of εðtÞ the boundary layer can be described by the following approximation of the Heaviside function: Figure 2 shows the approximated Heaviside function for the initial case, the stationary case and an intermediate value.
Since for t > t 0 the value H app ðx; tÞ�0 in ½À εðtÞ; 0� the initial profile has to be prolonged. Therefore we set: Now the approximation # S;app can be defined by: For evaluating εðtÞ and hence H app ð�; tÞ the parameters ε 0 ; ε 1 ; τ must be determined. To do so we consider for m; n 2 N the set: The values # S;fit ðx i ; t j Þ can be obtained by measuring or from the exact solution (26). Now we define the cost function F by: Minimizing the cost function F leads to an optimal set of the parameters ε 0 ; ε 1 ; τ. Such a nonlinear optimization problem usually can be solved for example with the method of Gauß-Newton [18]. In this case, an application of the Gauß-Newton method would have the problem that H app is not continuously differentiable with respect to the parameters and so a requirement for the application of the Gauß-Newton method is not given. Therefore, the following procedure can be applied to minimize the cost function F. Given a set of possible values of the parameters and K 0 ; and set: This means, the chosen values of the parameters are those that minimize the cost function in a given discrete set of possible values.
As an alternative to the minimization of the cost function F, the following approach leads to explicit formulas for the determination of the parameters. Since the time derivative of # S is an indicator for the thickness of the boundary layer, we will match the time derivatives of # S and # S;app in the two points: These two points are chosen to get a value near the bottom and a value on top of the storage tank. Because of the matching leads to: With @ @t and (47) for i ¼ 1; 2 we get: Here it is assumed that the one side derivative # 0 0 ð0Þ exists. The system (61) can be interpreted as a linear system for the parameters ε 0 and ε 1 with the solution: In (62) the time constant τ is still unknown. The duration of the discharging process can be estimated by: Numerical experiments with various velocities indicated that an appropriate choice for the time constant τ seems to be given by T D =2. Hence we set τ ¼ T D =2. With (62) the parameters τ; ε 0 ; ε 1 now can be computed by: While the just developed method requires the knowledge of the exact solution, the minimization of the cost function F can be made with measured data. Since the computation of the exact solution takes a long time, especially for large values of the velocity, the computation of the cost function F will produce high numerical costs if the exact solution is used for the data # S;fit ðx i ; t j Þ. So the minimization of the cost function should be applied if measured data are used and the method of matching the derivatives is more appropriate if the exact solution is used. In the following Section the method of matching the derivatives is used for the determination of the parameters. For a better overview, in Table 2 the parameters and variables concerned with the approximate solution # S;app are collected.

Validation
In this Section the approximate solution # S;app will be validated. For this the approximate solution will be compared on the one hand with the exact solution (26) and on the other hand with a numerical solution. For the numerical solution we use a semi-discretization (Method of Lines [2]), that means only the space variable x will be discretized. The interval ½0; L� will be partitioned with the equidistant step size Δx, that means: The first derivative in space will be approximated by a backward difference and the second space derivative will be approximated by a central difference, hence: In the following Examples we consider a storage tank with the data: All computations are made with the software Mathematica Version 9. For the numerical computation of the eigenvalues λ k the FindRoot routine is used and the computation of the generalized Fourier coefficients c k is done with the Integrate routine. The initial value problem (70) -(73) is treated with the NDSolve routine.

Example 1.
We already know from Section 3 that for 16 available digits and 4 reserved digits for accuracy, condition (42) leads to velocities v D < 55a W =L that can be treated. Hence, in this first example we consider the velocity: This velocity corresponds to the volume flow _ V ¼ π R 2 v D ¼ 8l=h and hence to a very low tapping flow rate. With (64), (65) and (66) we get the following values for the parameters of the approximate solution: For the computation of the exact solution 200 terms in the sum are considered and for the numerical solution n ¼ 50 is taken. In Figure 3 the comparison at the points x ¼ L=10, x ¼ L=2 and x ¼ L can be seen, where # S;app and # S;num are dashed indicated. Figure 4 shows the errors in the three considered points.

Example 2.
In this second example we take the velocity: This velocity corresponds to the volume flow _ V ¼ 29l=h and hence to a still very low flow rate.
Condition (42) here yields N > 43. We choose a precision of 60 digits and take 500 terms in the sum for the computation of the exact solution.
With (64), (65) and (66) we get the following values for the parameters of the approximate solution: For the numerical solution n ¼ 300 is taken. Figures 5 and Figures 6 show the corresponding values as in Example 1.

Example 3.
In this third example we take the velocity: and hence _ V ¼ 290l=h. This volume flow is a value in the middle region of typical tapping flow rates.  Condition (42) here yields N > 434. We choose a precision of 600 digits and for the computation of the exact solution 5000 terms in the sum are considered.
With (64), (65) and (66) we get the following values for the parameters of the approximate solution: The Examples show that the approximate solution can fit the form of the curves of the exact solution throughout the whole discharging process for all the values of the velocity that are considered here, although it has a quite simple structure. Especially, it can be seen that the approximate solution improves the accuracy for the larger values of the velocity, and hence for realistic tapping flow rates. The number n for the numerical solution was chosen such that the errors are in the same region. Larger values of n lead to better results of the numerical solution.
The results of this Section show that the approximate solution is an alternative to numerical solutions for simulating a discharging process. In contrast to numerical procedures for the approximate solution no stability analysis must be made, the  approximate solution is condition-less stable and independent from any step size. While for implicit numerical procedures in every time step a system of equations has to be solved, the approximate solution can be evaluated instantly. Because of its simple structure the approximate solution has very low numerical costs and yields a remarkable accuracy.

Other operational modes
A discharging process can be treated with the results that have been developed in the previous Sections. In an analogous way a charging process leads to the initial and boundary value problem:    with the charging velocity v C and the charging temperature # Load . The complete analysis can be adapted and leads to the unique solution # S defined by: Analogously, here we get the approximation: with ωðx; tÞ ¼ L À x À v C ðt À t 0 Þ and: Neglecting heat conduction, a stand-by period can be described by the initial value problem: Here # A is the constant ambient temperature in the room where the storage is placed and b S describes the thermal losses through the wall of the storage. The unique solution of this initial value problem is given by: Solution (92) matches exactly with the benchmark solution given in [19], this means that this solution would pass the benchmark test given in [19] for every time step. So for the three operational modes that are mentioned in Section 1 analytical approximations can be used for describing the process.
It should be noted that in reality also situations can appear where one of the heat sources (boiler or plant) is charging the storage and at the same time a request of thermal energy is demanded by tapping hot water. Also it can appear that the thermal energy of the storage is not sufficient for a large demand of hot water and so a heat source (the boiler for example) must supply the requested thermal energy. For modelling such situations the control strategy of the underlying heating system must be known. So these situations are not considered here.

Conclusion
The approximate solution constructed in this article can be used for applying in a simulation of the dynamic thermal behaviour of a storage tank as an alternative to numerical procedures. Because of the simple structure of the approximation a simple implementation, comparably low numerical costs and condition-less stability will be advantages.
The combination of low numerical costs and although a high accuracy of the approximate solution is a property that could be very helpful for testing new control strategies. Because of the low numerical costs, annual simulations can be made in an adequate time and so various control strategies could be compared by evaluating the energy savings over the whole year, for example.
For future research it will be interesting to apply the approximation method to hot water storage tanks with an immersed coil [5]. Such storage tanks are also often used in domestic installations and in combination with a solar plant or a micro combined heat and power plant. In such a storage tank usually the immersed coil is positioned in the lower region of the storage and works as an internal heat exchanger. So the storage can be divided in two parts, an upper part and a lower part. The upper part can be treated as a storage tank that is investigated in this article. The lower part can be considered as a heat exchanger with a large thermal mass.
In contrast to conventional heat exchangers with low thermal masses, here the thermal mass cannot be neglected and hence a PDE system for the temperature in the coil and the temperature in the storage has to be solved. Since in real storage tanks natural convection leads to a natural stratification, a mathematical model has to avoid the phenomenon of temperature inversion. Numerical approaches for such storage tanks usually apply temperature inversion correction methods [5,6].
For applying an analytical approximation, the storage tank could be treated with the results of this article for the discharging process and the stand-by mode. Since the charging process is realized by the immersed coil, the method from this article cannot be taken. Here the PDE system in the lower part has to be approximated and the interaction between the lower and the upper part during a charging process has to be described adequately. The analytical solution of the initial and boundary value problem for the PDE system seems to be a quite difficult problem. So the investigation of a storage tank with an immersed coil will be a very interesting challenge for future research.

Disclosure statement
No potential conflict of interest was reported by the author.