Mathematical modeling of unsteady response of plate and fin heat exchanger to sudden change in liquid flow rate

A mathematical simulation of plate fin and tube heat exchanger will be presented. The transient operation of the heat exchanger was carried out using general numerical model developed in [1]. Reynolds number of the water flowing inside the tubes varied in the range from 4000 to 12000. A detailed analysis of transient response was modeled for growth in the water volume flow in time. At first, heat transfer correlations for air and water were determined based on the experimental data. Unknown parameters appearing in the relationships for the Nusselt numbers on the air and water-sides were estimated using the least squares method. The power-type form of the relationship for the air-side Nusselt number was used while two correlations of different form were chosen for the water-side Nusselt number. The first is the exponential correlation, and the second form is similar to the relationships of Petukhov-Kirillov and Gnielinski. These correlations were used in the mathematical model of the heat exchanger for the simulation of its transient operation. The results of the numerical simulations of a heat exchanger using experimentally determined air and water-side heat transfer formulas for calculation of heat transfer coefficient were compared with the experimental data. Very good agreement of computation results (i.e. air and water temperature at the outlet of the heat exchanger) with the experimental data was obtained.


Introduction
Publications on the dynamics on the plate fin and tube heat exchangers (PFTHE) are not numerous.Mathematical models of heat exchangers to simulate the non-steady state operation are needed to analyze heat exchangers startup and shutdown, and in PID control systems as well as in control systems using mathematical models of the controlled heat exchanger.
However, most of the works refers to a steady state of these exchangers [2].
The transient response of PFTHE was modeled in [1,[3][4][5][6].Usually, fins are modeled as elements with lumped thermal capacity while the equivalent heat transfer coefficient on the air side is calculated as in the steady-state [1,[3][4][5][6].The steady-state temperature distribution in fins is used to calculate the equivalent heat transfer.This simplification can lead to some errors, especially for higher fins.
The system of the governing differential equations is solved by the Laplace transform or the finite difference method.The Laplace transform method is widely used by Roetzel and Xuan [3] to model the transient operation of various heat exchangers.Taler [5] modeled the transient response of the PFTHE for a stepwise increase in liquid or gas temperature using the Laplace transform method.The transient temperature of the fluids was compared with temperature obtained by the finite volume method.
Korzeń and Taler [1] developed a new mathematical model of the PFTHE to simulate its transient operation.In contrast to the existing methods for the modeling transient response of heat exchangers with extended surfaces, in which the weighted steady-state heat transfer coefficient on the finned tube side is used, the transient temperature distribution will be calculated in each fin.This allows for computing more exactly the heat flow rate from the fins to the flowing gas.
Usually, modeling PFTHEs assumed that the fluid flow inside the tubes is turbulent.However, many low-duty heat exchangers operate in the transitional region when the Reynolds number varies in the range of 2 300 to about 12 000.In this paper, numerical modeling of transient work of a PFTHE for low Reynolds numbers on the liquid side will be carried out.

Mathematical formulation of the problem
The system of partial differential equations describing the space and time changes of liquid T1, tube wall Tw, and air T2 temperatures in one-row plate-fin and tube heat exchanger are [1]: • liquid • air where 2 T denotes the mean air temperature over the row thickness, defined as The symbols ch x x L   and 2 y y p   in equations (1-3) stand for dimensionless coordinates.The numbers of heat transfer units N1 and N2 are given by where: The time constants 1 The symbols in equations (1-6) denote: The subscript w refers to the wall, f to the fin, and m to the mean value.
The weighted heat transfer coefficient ho is defined by The initial temperatures of both fluids are equal and amount to T0.The initial conditions are The boundary conditions have the following form where f1(t) and f2(t) are functions of time describing the variation of the liquid and air temperature at the inlets to the exchanger.The fin efficiency appearing in the relationship (7) was calculated using the Finite Volume Method based on the Finite Element Method.The division of the fin into finite volumes with the number of nodes is shown in Fig. 2.
The initial-boundary value problem defined by equations (1-14) applies to heat exchangers made of bare or finned tubes.For bare tubes mf is equal to zero because there are no fins.
The transient temperature of the fluids and tube wall in the one-row heat exchanger (Fig. 1) were determined by the explicit finite difference method.To calculate timedependent efficiency f  of the rectangular fin attached to an oval tube the finite volume - finite element method (FVM-FEM) was used.Based on the mathematical model of the onerow heat exchanger, a mathematical model of the whole heat exchanger with a complex flow system was built.
The automotive radiator for the spark-ignition combustion engine with a capacity of 1580 cm 3 is a double-row, two-pass plate-finned heat exchanger.The radiator consists of aluminium tubes of oval cross-section.The cooling liquid flows in parallel through both tube rows.In the developed mathematical model of the exchanger, the heat transfer coefficients on the air and water side are calculated from experimental correlations.The correlations for the heat transfer coefficient were derived assuming a fully developed turbulent flow in straight tubes.One of the most common correlations is the Dittus-Boelter formula [7][8][9]: where n = 0.4 when the fluid is heated and n = 0.3 when it is cooled.In recent years, the Gnielinski correlation is becoming more and more popular [10][11].
where the friction loss coefficient  is defined as: In Formula ( 16), symbols dr and Lch denote the tube hydraulic diameter and length, respectively.The reason for the growing popularity of the Gnielinski formula ( 16) is its greater accuracy in the transition area of 2300 < Rew < 10000.Despite the fact that the Gnielinski proposal gives smaller values of the heat transfer coefficient in the transition area compared to the Dittus-Boelter formula, these values are still inflated [12][13].If Rew < 2300, the flow is laminar; the Nusselt number is Nuw = 4.364 and does not depend on the Reynolds number [7][8].It can be checked easily that both the Dittus-Boelter and the Gnielinski formulae give values of Nuc much higher than 4.364.For this reason, the correlations for the Nusselt numbers on the water and air side were determined based on experimental testing.
The following form of the correlation for the Nusselt number at the air side was taken: Considering that the range of changes in the Reynolds number at the air side is not too big, it is possible to approximate experimental data well using the exponential equation written above (18).
The correlation for the Nusselt number at the water side was taken in a form similar to the formulae developed by Petukhov and Kirillov [14] and Gnielinski [10].
The values of coefficients 1 2 , , ..., m x x x will be selected so that the sum of the squares of differences in temperature [%] in the following intervals [15][16]: ,..., ,..., Where n -number of measuring points and m -number of sought parameters.
Quantities cii are diagonal elements cii of the covariance matrix.
In this case, the number of measuring points is n = 47, and the number of unknown parameters is m = 4.

Experimental correlations
The thermal and flow measurements of the exchanger were performed for different air and water flow velocities.The measurement results are presented in [17].
The Reynolds number values at the air side: .The air physical properties were calculated at the air average temperature.Similarly, the thermophysical properties of water were calculated at the average temperature.The average temperature is understood as the arithmetic mean of the average temperature at the exchanger inlet and outlet.The water velocity in the upper pass wwg is smaller than the water velocity in the lower pass wwd because in the upper pass water flows in parallel through 20 tubes whereas in the lower one -through 18.
Based on the data presented in In the mathematical model of the radiator operation under unsteady-state conditions the relationships (24) and (25) were used.

Increase in water mass flow in time
The comparison of the results obtained from measurements and calculations is presented in Fig. 3.The increase in the water volume flow causes a rise in the water and air temperatures at the exchanger outlet (Fig. 3).The increment in the air temperature is bigger than the increment in the water temperature.This results from the smaller flux of the air heat capacity a pa m c  compared to the heat capacity flux of water w w m c  .If the heat flux given up by water is the same as the one absorbed by air, the increment in the air temperature is higher.Fig. 4 presents changes in the water and the tube wall temperatures on the exchanger length in time t = 119 s (Fig. 3).Analysing the results presented in Fig. 4, it can be seen that water in the first row tubes cools down more than in the tubes placed in the second row.This concerns both the upper and lower pass.It can be noticed that the tube wall temperature is lower than the water temperature due to the fact that heat is absorbed by cold air.The curves illustrating time-dependent changes in the temperature of selected fins are presented in Fig. 5a -d.An analysis of the results indicates that the highest temperature is at the fin base, and the lowest in node 3 (Fig. 2), which is the most distant from the tube axis.It can be seen that the temperatures of the fins in the first tube row (Fig. 5a and 5c) are lower than in the second row (Fig. 5b and 5d).Cold inflowing air cools the fins and tubes located in the first row better.The air temperature after the first tube row is much higher compared to the inlet air temperature.Due to a reduction in the difference between the temperatures of the flowing air and the tubes located in the second row, the temperature of the water and fins in the second tube row is higher than in the first row.The calculations and measurements carried out in the paper confirm that the developed numerical model of the plate fin and tube heat exchanger has a high accuracy.This model is particularly useful to simulate exchanger, where the temperature of the fins can be very high.Thanks to the experimental determination of heat transfer correlations on the air and water sides a very good agreement of calculations and measurements was reached.The developed model can be used to simulate the heating and cooling heat exchanger and in automatic control

Fig. 1 .
Fig. 1.Scheme of the one-row heat exchanger (a) and control volume (b).

Fig. 2 .
Fig. 2. Division of the fin model into finite volumes with node numbers.

Fig. 3 .
Fig. 3.The exchanger unsteady-state operation caused by a sudden reduction in the water volume flow; 1w V   -water volume flow at the exchanger inlet, l/s 2am T  -average air temperature before the exchanger, ˚C 3am T  -average air temperature after the exchanger, ˚C 4w T  -water temperature before the exchanger, ˚C 5w T  -water temperature after the exchanger, ˚C 6 -0 w -air velocity before the exchanger, m/s 7w T  -calculated water temperature after the exchanger, ˚C 8am T  -calculated average air temperature after the exchanger, ˚C.

Fig. 4 .
Fig. 4. Distribution of the water and the tube wall temperatures on the exchanger length for time t = 119 s in the steady state.

Fig. 5 .
Fig. 5. Changes in the fin base temperature Tb and the fin temperature in 12 nodes, that are marked in Fig. 2, at a distance of 247 mm (9.5 cell lengths) from each pass inlet: (a) -first tube row in the first (upper) pass, (b) -second tube row in the first (upper) pass, (c) -first tube row in the second (lower) pass, (d) -second tube row in the second (lower) pass.