Predictor-corrector algorithm for solving quasi-separated-flow and transient distributed-parameter model for heat exchangers

Abstract

The successive sub〈stitution (SS) method is a suitable approach to solving the transient distributed-parameter model for heat exchangers. However, this method must be enhanced because its convergence heavily depends on the iterative initial pressure. When the iterative initial pressure is improperly assigned, the calculated flow rates become negative values, causing the state parameters to exhibit negative values as well. Therefore, a predictor-corrector algorithm (PCA) is proposed to improve the convergence of the SS method. A predictor is developed to determine an appropriate iterative initial pressure. Total fluid mass is adopted as the convergence criterion of pressure iteration instead of outlet flow rate as is done in the SS method. Convergence analysis and case study of the PCA and SS method are conducted, which show that the PCA has better convergence than the SS method under the same working conditions.

Keywords

algorithm / convergence / heat exchanger / modeling / transient

Introduction

Almost every heat exchanger operates under variable working conditions; thus, an absolute steady state does not exist, raising the necessity of investigating the transient characteristics of heat exchangers [1]. Because experimental tests are time consuming and expensive, a faster and less costly alternative that uses mathematical models to simulate the transient characteristics of heat exchangers is needed. According to parameter-lumped characteristics, three kinds of mathematical models, namely, the lumped parameter model [2], the zone model [3-5], and the distributed parameter model [6-16], have been successively developed. Of these, the distributed parameter model has been drawing increased attention in recent years because it can provide the most complete information and best insight into the dynamic behavior of heat exchangers. Table 1 summarizes the transient distributed-parameter models for heat exchangers available in literature. There are three types of distributed-parameter models according to the difference in simplifying two-phase fluid: the homogeneous-flow model [6-8], the separated-flow model [9], and the quasi-separated-flow model [10-16]. The homogeneous-flow model assumes thermodynamic equilibrium and nil slip between liquid and vapor phases. The algorithm of this model can easily be devised; however, the calculated mass of the two-phase fluid is different from the actual mass, even causing errors in multiple orders of magnitude [17]. The separated-flow model has fewer assumptions and can reflect the true character of the working fluid. Although the simulation results may be more precise; however, considerably more computing time is required; poor stability exists in solving the governing equations of each phase fluid [18]. For the quasi-separated-flow model, an appropriate void fraction model is adopted to reflect the velocity difference between two fluids. Both the accuracy and calculation speed of this model are between those of the homogeneous-flow and separated-flow models. Therefore, the quasi-separated-flow and distributed-parameter models are recommended for transient simulation of heat exchangers [10-16].

Table 1 shows that two types of algorithms have been reported for solving the quasi-separated-flow and transient distributed-parameter model. The first is the simultaneous equation solving (SES) method [10-15], in which all governing equations of one control volume are combined into one equation set and solved simultaneously using the Euler or Newton-Raphson method, etc. The calculations are performed from the inlet to the outlet of the heat exchanger, one control volume after another. During the equation numerically solving process, there is no physical meaning; thus, the calculation stability of this method is not easy to ensure [18]; it is difficult for users to detect the divergence cause that may have happened in the calculation. The second algorithm, the successive substitution (SS) method [13-15], uses a number of balance conditions, such as mass or energy balance, as convergence criteria. A set of initial values are assumed, calculation begins from the innermost cycle, and other parameters are generated. If the convergence criterion is not satisfied, the previously assumed initial values are updated, and then, the iteration is repeated. This method has obvious physical meanings; thus, debugging and ensuring calculation stability are easy tasks. Theoretically, the SS method is more effective and straightforward than the SES method. However, in practice, the problem on convergence persists when using the former. The convergence of the SS method heavily depends on the iterative initial pressure, especially when the outlet fluid flow rate is set as the convergence criterion of the pressure iteration. When the iterative initial pressure is improperly assigned, this may generate negative values for the calculated flow rates. As a consequence, the state parameters likewise become negative values, making the method unsuitable. Therefore, the SS method requires further improvement.

In this letter, we propose a predictor-corrector algorithm (PCA) to improve the convergence of the SS method, as well as present the convergence analysis of the PCA case studies.

Heat exchanger model and problem analysis of the SS method

Heat exchanger model and the SS method

The heat exchanger model, which can be simplified as one equivalent tube model (Fig. 1), is not limited to a specific type of heat exchanger. Generally, the first fluid (two-phase) is assumed to flow inside the tube, and the second fluid (single-phase) flows outside the tube. The flow configurations of the two fluids may arbitrarily be one of counter, parallel, or cross flow.

The heat exchanger can be divided into n control volumes along the tube. Each control volume includes the first fluid, the tube wall, and the second fluid. The following recognized assumptions [10,11] are listed: ① the first fluid is one-dimensional flow; ② the pressure drop in the two-phase fluid can be neglected; ③ the phases in the two-phase fluid are in thermodynamic equilibrium; and ④ the vapor and liquid have an average cross-sectional velocity, and the void fraction correlation is used to describe the ratio of vapor cross-sectional area to the total cross-sectional area. Based on these assumptions, the governing equations of each control volume are simplified as follows:

The mass equation of the first fluid is given as

(1){\displaystyle \frac{\mathrm{d}{M}_{\mathrm{i}\mathrm{s},i}}{\mathrm{d}t}}=m_{\mathrm{i}\mathrm{s},i}-m_{\mathrm{i}\mathrm{s},i+\mathrm{1}}.

The energy equation of the first fluid is given as

(2){\displaystyle \frac{d{(Mh)}_{is,i}}{dt}}=m_{is,i}{h}_{is,i}-m_{is,i+\mathrm{1}}{h}_{is,i+1}+Q_{is,i}.

The energy equation of the second fluid is given as

(3){\displaystyle \frac{d{(M{c}_{p}T)}_{os,i}}{dt}}=m_{\mathrm{os}}{c}_{p,os}{T}_{os,i+1}-m_{\mathrm{os}}{c}_{p,os}{T}_{os,i}-Q_{os,i}.

The energy equation of the tube wall is given as

(4){\displaystyle \frac{d{(M{c}_{p}T)}_{wall,i}}{dt}}=Q_{os,i}-Q_{is,i}.

For the heat exchanger model, the known values of the first fluid are inlet mass flow rate m_is,in and inlet enthalpy h_is,in. The known values of the second fluid are mass flow rate m_os and inlet temperature T_os,in. There are four unknown values in the governing equations: the outlet mass flow rate m_is,i+1, the outlet enthalpy h_is,i+1 of the first fluid, the temperature T_wall,i of the tube wall, and the temperature T_os,i of the second fluid. Other unknown variables can de deduced from the aforementioned four unknown variables. The number of unknown values is equal to the number of governing equations; thus, the set of governing equations can be solved.

The SS method is used to solve the abovementioned heat exchanger model. In the numerical formulation, a fully implicit scheme is chosen to ensure computational stability. All first-order spatial derivatives and time derivatives in the governing equations are discretized using the backward difference scheme. Thus, the discretization of Eqs. (5)-(8) is obtained as follows:

Discretizing Eq. (1), the outlet mass flow rate of the first fluid is obtained as

(5);m_{is,i+\mathrm{1}}=m_{is,i}-(M_{is,i}-M_{is,i}^0)/\mathrm{\Delta }t.

Discretizing Eq. (2), the outlet enthalpy of the first fluid is obtained as

(6);h_{is,i+1}={\displaystyle \frac{M_{is,i}^{\mathrm{0}}{h}_{is,i+\mathrm{1}}^{\mathrm{0}}+m_{is,i}{h}_{is,i}\mathrm{\Delta }t+Q_{is,i}\mathrm{\Delta }t}{M_{is,i}+m_{is,i+1}\mathrm{\Delta }t}}.

Discretizing Eq. (3), the outlet temperature of the second fluid is obtained as

(7);T_{os,i}={\displaystyle \frac{M_{os,i}^{\mathrm{0}}{c}_{p,os}{T}_{os,i}^{\mathrm{0}}+m_{os,i+\mathrm{1}}{c}_{p,os}{T}_{os,i+\mathrm{1}}\mathrm{\Delta }t+f_{\mathrm{os}}{A}_{\mathrm{os}}{T}_{wall,i}\mathrm{\Delta }t}{M_{os,i}{c}_{p,os}+m_{os,i}{c}_{p,os}\mathrm{\Delta }t+f_{\mathrm{os}}{A}_{\mathrm{os}}\mathrm{\Delta }t}}.

Discretizing Eq. (4), the tube wall temperature is obtained as

(8).T_{wall,i}={\displaystyle \frac{M_{wall,i}{c}_{p,wall}{T}_{wall,i}^{\mathrm{0}}+f_{os,i}{A}_{os,i}{T}_{os,i}\mathrm{\Delta }t+f_{is,i}{A}_{is,i}{T}_{is,i}\mathrm{\Delta }t}{M_{wall,i}{c}_{p,wall}+f_{os,i}{A}_{os,i}\mathrm{\Delta }t+f_{is,i}{A}_{is,i}\mathrm{\Delta }t}}.

Figure 2 shows the one-time step flow chart of the SS method. The main idea is to use the outlet flow rate from the last control volume to check whether the assumed pressure at the beginning of the computational step is correct. However, there is a problem on convergence for the SS method, as explained in Section 2.2 below.

Problem analysis of SS method

Whether the SS method can reach convergence successfully depends on two conditions: ① the monotone function between the pressure and the outlet flow rate and ② the appropriate iterative initial pressure. Condition ① is satisfied for the heat exchanger model. This monotone function ensures that there must be one pressure that can enable the pressure iteration to reach convergence. Therefore, an appropriate iterative initial pressure is the key to the convergence of the SS method. In practice, the convergence of the SS method heavily depends on the iterative initial pressure.

The convergence analysis of the SS method was performed in two steps. First, the method of reduction to absurdity is used for the heat exchanger with one control volume. Second, the analysis result is extended to the heat exchanger with multiple control volumes.

For the heat exchanger with one control volume, the boundary conditions of the first fluid are inlet enthalpy h_is,in, inlet mass flow rate m_is,in, and outlet mass flow rate m_is,out, which are all positive values.

The other way of calculating outlet mass flow rate m_{is ,out} is obtained with Eq. (9), which is compared with the outlet mass flow rate of boundary condition, which is a convergence criterion (Fig. 2).

(9)m_{is,out}=m_{is,in}-(M_{\mathrm{is}}-M_{\mathrm{is}}^{\mathrm{0}})/\mathrm{\Delta }t.

As the inlet mass flow rate m_is,in, the total mass of the previous time step M_{\mathrm{i}\mathrm{s}}^{\mathrm{0}}, and the time step Δt are determined; the outlet mass flow rate m_is,out is related only to the total mass M_is of the current time. Here, the method of reduction to absurdity is adopted. Let m_is,out of Eq. (9) be less than zero; then Eq. (10) is obtained.

(10)M_{\mathrm{is}}>M_{\mathrm{is}}^{\mathrm{0}}+m_{is,in}\mathrm{\Delta }t\mathrm{.}

As (m_is,inΔt) and M_{\mathrm{is}}^{\mathrm{0}} have the same dimensions, the right side of Eq. (10) is combined into unity, namely, M_{\mathrm{Lim}}=M_{\mathrm{is}}^{\mathrm{0}}+m_{is,in}\mathrm{\Delta }t.

Equation (10) is arranged, and Eq. (11) is obtained.

(11)M_{\mathrm{is}}>M_{\mathrm{Lim}}\mathrm{.}

The fluid mass is the function of the pressure, increasing with the rise in pressure when its enthalpy is constant. In other words, M_is = f (p_is) is a monotone function. Thus, Eq. (11) can be changed into the representation form of the corresponding pressure:

(12)p_{\mathrm{is}}>p_{\mathrm{Lim}}\mathrm{.}

Thus, we conclude that when the assigned iterative initial pressure p_is is higher than the limit pressure p_Lim, the calculated outlet mass flow rate m_is,out will be negative, which is in direct contradiction with the positive value m_is,out of the boundary condition.

For the heat exchanger with multiple control volumes, the analysis result above is extended. When Eq. (12) is met, the mass flow rate in the middle of the heat exchanger may be a negative value. Thus, the fluid flows backward relative to the flow direction of inlet and outlet fluids from the physical meaning, which cannot happen practically. In addition, because of the negative mass flow rate, the negative enthalpy is obtained from Eq. (6) in the calculation of the succeeding control volumes, making the state parameters lose their physical meaning. Thus, the SS method becomes invalid.

In other words, the convergence of the SS method heavily depends on the iterative initial value. When the iterative initial value of pressure is improperly assigned, this can lead to a negative value of the mass flow rate, which causes the problem of nonconvergence. To improve convergence, the SS method should also be improved.

New method: PCA

The newly proposed method is named PCA, which includes two parts: the predictor and corrector. An appropriate iterative initial pressure is obtained in the predictor, in which the total mass is set as the convergence criterion. The final mass flow rates of each control volume and other state parameters are calculated in the corrector. The flow chart of the PCA is illustrated in Fig. 3. The detailed steps of the PCA are listed below:

1) Predictor

In this part, the main objective is to determine an appropriate iterative initial pressure for the corrector.

Step 1: The inlet and outlet mass flow rates of the previous time step are adopted as the corresponding values of the current time step for each control volume.

Step 2: The pressure p_is and the tube wall temperature T_wall,i of each control volume are assumed.

Step 3: The outlet enthalpy of each control volume is calculated according to Eq. (6). Other state parameters, such as density, etc., can be obtained according to the pressure and the enthalpy.

Step 4: The mass of each control volume M_is,i is calculated and accumulated into total mass, that is, M_is(1).

Step 5: The temperature T_os,i of the second fluid and the temperature T_wall,i of the tube wall are calculated in sequence, according to Eqs. (7) and (8).

Step 6: The tube wall temperatures obtained from Steps 2 and 5 are compared. If they are not equal, the tube wall temperature is adjusted, and we return to Step 3 until the convergence of the tube wall temperature is reached.

Step 7: The total mass M_is(2) of the first fluid is calculated according to Eq. (13).

Equation (1) is integrated until the total mass of another form is obtained.

(13)M_{\mathrm{is}(\mathrm{2})}=M_{is,total}^{\mathrm{0}}+(m_{is,in}-m_{is,out})\mathrm{\Delta }t

where m_is,in, m_is,out, and M_{is,total}^{\mathrm{0}} are known values in Eq. (13). When the time step is given, the total mass is fixed and unique.

Step 8: The total masses of Steps 4 and 7 are compared. If they are not equal, the pressure p_is is adjusted, and we return to Step 2 until the convergence of the total mass is reached.

The determined pressure in the predictor is an appropriate iterative initial value only for the corrector.

2) Corrector

In this part, the mass flow rates and other state parameters of each control volume are calculated correctly. The flow chart of the corrector is the same as that in Fig. 2. Then, the one-time step PCA is completed.

Convergence analysis of PCA

The three factors ensuring the better convergence of the PCA compared with that of the SS method are listed: ① the total mass of the first fluid is unique and constant in each time step; ② the total mass is independent of the mass flow rate distribution inside the heat exchanger and is only related to the inlet mass flow rate, outlet mass flow rate, and the total mass of the previous time step, which are all known values; ③ the total mass and pressure are monotone functions [19], which implies that only one mass corresponds to one pressure.

In the PCA predictor, Eq. (13) is used to calculate the total mass. In the SS method, when the negative outlet mass flow rate is obtained because of a higher iterative initial pressure, the total mass is satisfied with Eq. (10).

Compare the right side of Eqs. (10) and (13), and the following result can be obtained:

(14)

Equation (14) shows that the total mass of the PCA must be less than that of the SS method. As M_is = f (p_is) is a monotone function, the iterative initial pressure determined in the predictor must be satisfied with Eq. (15):

(15)

Equation (15) guarantees that the calculated outlet mass flow rate of any control volume will not be a negative value. Thus, the PCA can reach convergence successfully.

Application of PCA into the heat exchanger simulation

A specified double-tube heat exchanger is chosen as the simulation object, whose detailed structure parameters are referred to in Ref. [15]. Refrigerant 134a flows inside the inner tube and water flow outside. The dynamic condition is obtained by introducing a step decrease in the refrigerant flow rate. The boundary conditions are shown in Table 2, showing that the auxiliary equations are needed in the simulation. The Dittus-Boelter correlation is adopted for the heat transfer coefficient in a single phase. The heat transfer coefficient in two phases adopts the correlation of Ref. [20]. The calculation method of the refrigerant property is referred to in Ref. [21]. The total heat exchanger is divided into 200 control volumes. The time step is set to 0.5 s.

Transient behavior was simulated in the SS method and the PCA, respectively. The iterative initial pressure was assigned as 800 kPa, which was higher than the limited pressure p_lim and satisfied Eq. (12). For the SS method, the obtained mass flow rates decrease along the tube to negative values because of the higher iterative pressure in the first time step. With the negative mass flow rate, the calculated enthalpy decreases to negative values in the succeeding calculation of control volumes. The SS method thus becomes invalid. The mass flow rate and enthalpy distribution are shown in Fig. 4(a) according to the SS method, which is coincident with the analysis of Section 2.2.

However, the PCA can work normally under the same working conditions and iterative initial pressure. The mass flow rate and enthalpy distribution are shown in Fig. 4(b) after reaching a steady state. Although the assigned iterative initial pressure is higher, it can be adjusted and finally converge to an appropriate value according to the monotone function between the total mass and the pressure in the predictor, which ensures the successful convergence of the PCA.

Conclusion

This letter provides a new algorithm, the PCA, for solving quasi-separated-flow andtransient distributed-parameter model for heat exchangers. The PCA includes two parts: the predictor and corrector. The new method can overcome the divergence problem of the SS method. In the predictor, an appropriate iterative initial pressure is obtained according to the monotone function of the total mass and pressure; in the corrector, the mass flow rates and state parameters are calculated correctly. The case study shows that the PCA improves algorithm convergence and enhances the applicability of the algorithm. The PCA is not developed for one fixed heat exchanger type; thus, the new method can be extended to any type of heat exchanger.

Notation

References

[1] Wedekind G L, Stoecker W F. Transient response of the mixture-vapor transition point in horizontal evaporating flow, ASHRAE Transactions, 1966, 72(II): IV.2.1–IV.2.15

[2] Chi J, Didion D. A simulation model of the transient performance of a heat pump, International Journal of Refrigeration, 1982, 5(3): 176–184

[3] Fu L, Ding G L, Zhang C L. Dynamic simulation of air-to-water dual-mode heat pump with screw compressor. Applied Thermal Engineering, 2003, 23(13): 1629–1645

[4] Pettit N B O L, Willatzen M, Plong-Sørensen L. A general dynamic simulation model for evaporators and condensers in refrigeration part II: simulation and control of an evaporator. International Journal of Refrigeration, 1998, 21(5): 404–414

[5] Zhang W J, Zhang C L. A generalized moving-boundary model for transient simulation of dry-expansion evaporators under larger disturbances. International Journal of Refrigeration, 2006, 29(7): 1119–1127

[6] Nyers J, Stoyan G. A dynamic model adequate for controlling the evaporator of a heat pump. International Journal of Refrigeration, 1994, 17(2): 101–108

[7] Jia X, Tso C P, Chia P K, Jolly P. A distributed model for prediction of the transient response of an evaporator. International Journal of Refrigeration, 1995, 18(5): 336–342

[8] Aprea C, Renno C. An air cooled tune-fin evaporator model for an expansion valve control law. Mathematical and Computer Modeling, 1999, 30(7–8): 135–146

[9] Wang C C. A numerical method for thermally non-equilibrium condensing flow in a double-pipe condenser. Applied Thermal Engineering, 1997, 17(7): 647–660

[10] Wang H, Touber S. Distributed and non-steady-state modeling of an air cooler. International Journal of Refrigeration, 1991, 14(2): 98–111

[11] Jia X, Tso C P, Jolly P, Wong Y W. Distributed steady and dynamic modeling of dry-expansion evaporators. International Journal of Refrigeration, 1999, 22(2): 126–136

[12] Tso C P, Cheng Y C, Lai A C K. Dynamic behavior of a direct expansion evaporator under frosting condition Part I. Distributed model. International Journal of Refrigeration, 2006, 29(4): 611–623

[13] Koury R N N, Machado L, Ismail K A R. Numerical simulation of a variable speed refrigeration system. International Journal of Refrigeration, 2001, 24(2): 192–200

[14] Mithraratne P, Wijeysundera N E, Bong T Y. Dynamic simulation of a thermostatically controlled counter-flow evaporator. International Journal of Refrigeration, 2000, 23(3): 174–189

[15] Mithraratne P, Wijeysundera N E. An experimental and numerical study of the dynamic behavior of a counter-flow evaporator. International Journal of Refrigeration, 2001, 24(6): 554–565

[16] Jiang H B, Radermacher R. A distributed model of a space heat pump under transient conditions. International Journal of Energy Research, 2003, 27(2): 145–160

[17] Rice C K. The effect of void fraction correlation and heat flux assumption on refrigerant charge inventory predictors. ASHRAE Transactions, 1987, 93(part 1): 341–367

[18] Ding G L. Recent developments in simulation techniques for vapour-compression refrigeration systems. International Journal of Refrigeration, 2007, 30(7): 1119–1133

[19] Ding G L. Dynamic simulation and optimization of small refrigeration equipments. Dissertation for the Doctoral Degree. Shanghai: Shanghai Jiao Tong University, 1993, 45–67

[20] Gungor K E, Winterton, R H S. A general correlation for flow boiling in tubes and annuli. International Journal of Heat Mass Transfer, 1986, 29(3): 351–358

[21] Ding G L, Wu Z G, Liu J, Inagaki T, Wang K J, Fukaya M. An implicit curve-fitting method for fast calculation of thermal properties of pure and mixed refrigerants. International Journal of Refrigeration, 2005, 28(6): 921–932

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg