A novel 1D/2D model for simulating conjugate heat transfer applied to flow boiling in tubes with external fins

This study presents a novel, simplified model for the time-efficient simulation of transient conjugate heat transfer in round tubes. The flow domain and the tube wall are modeled in 1D and 2D, respectively and empirical correlations are used to model the flow domain in 1D. The model is particularly useful when dealing with complex physics, such as flow boiling, which is the main focus of this study. The tube wall is assumed to have external fins. The flow is vertical upwards. Note that straightforward computational fluid dynamics (CFD) analysis of conjugate heat transfer in a system of tubes, leads to 3D modeling of fluid and solid domains. Because correlation is used and dimensionality reduced, the model is numerically more stable and computationally more time-efficient compared to the CFD approach. The benefit of the proposed approach is that it can be applied to large systems of tubes as encountered in many practical applications. The modeled equations are discretized in space using the finite volume method, with central differencing for the heat conduction equation in the solid domain, and upwind differencing of the convective term of the enthalpy transport equation in the flow domain. An explicit time discretization with forward differencing was applied to the enthalpy transport equation in the fluid domain. The conduction equation in the solid domain was time discretized using the Crank–Nicholson scheme. The model is applied in different boundary conditions and the predicted boiling patterns and temperature fields are discussed.

Abstract This study presents a novel, simplified model for the time-efficient simulation of transient conjugate heat transfer in round tubes. The flow domain and the tube wall are modeled in 1D and 2D, respectively and empirical correlations are used to model the flow domain in 1D. The model is particularly useful when dealing with complex physics, such as flow boiling, which is the main focus of this study. The tube wall is assumed to have external fins. The flow is vertical upwards. Note that straightforward computational fluid dynamics (CFD) analysis of conjugate heat transfer in a system of tubes, leads to 3D modeling of fluid and solid domains. Because correlation is used and dimensionality reduced, the model is numerically more stable and computationally more time-efficient compared to the CFD approach. The benefit of the proposed approach is that it can be applied to large systems of tubes as encountered in many practical applications. The modeled equations are discretized in space using the finite volume method, with central differencing for the heat conduction equation in the solid domain, and upwind differencing of the convective term of the enthalpy transport equation in the flow domain. An explicit time discretization with forward differencing was applied to the enthalpy transport equation in the fluid domain. The conduction equation in the solid domain was time discretized using the Crank-Nicholson scheme. The model is applied in different boundary conditions and the predicted boiling patterns and temperature fields are discussed.

Introduction
Nowadays, in the era of miniaturization, designers tend to minimize the size while maximizing the thermal efficiency and performance of heat exchange devices. This results in a considerable reduction in unit and maintenance costs. Consequently, the flow boiling process has attracted considerable scientific interest, because it allows thermal engineers to achieve very high heat transfer coefficients, compared to single phase convection [1][2][3][4][5][6]. The heat transfer augmentation associated with flow boiling, leads to the design of high performance heat exchange apparatus with a considerably smaller area than a single phase flow operation. However, flow boiling is a very complex process (including the nonlinear dynamics of bubble formation, unstable flow structures and different heat transfer regimes), which depends on many factors such as: mass flow rate, heat flux, pressure, surface tension and wall roughness. Therefore experimental investigation and numerical simulation were performed to explore the fluid flow and heat transfer phenomena. Recent developments in measurement and visualization techniques enable research at mini and micro scales as well as investigating complex flow structures. Recent scientific interest has focused particularly on two-phase flows in mini-and micro-channels [7][8][9][10][11] as well as on the boiling of nano-fluids to enhance heat transfer capability [12][13][14][15][16]. Experimental investigations have also been carried out on two-phase flows in the components of conventional heat exchange devices such as tubes, channels and walls. These investigations determined the critical heat flux [17,18], heat transfer coefficient and pressure drops in the horizontal and vertical tubes [18][19][20][21][22][23][24][25][26][27][28][29]. The improved correlations can be applied in the design of many devices, such as heat exchangers, evaporators and heat pipes.
Recent increases in computational power enables advanced numerical analysis, e.g. the three-dimensional transient simulations of heat and mass transport processes associated with two-phase flows. Because of its large heat transfer efficiency, there has been particular scientific interest devoted to modeling the nucleate boiling process [30][31][32][33][34][35], including the analyses of flow structures [36] and the nonlinear dynamics of bubble growth [37][38][39].
Nevertheless, the computational fluid dynamics (CFD) models have so far typically only been applied to simple geometries such as isolated tubes or cavities [38][39][40][41], and not to complex systems. This is because of the extreme physical complexity of the boiling phenomenon, and the associated numerical difficulties in stability and convergence. In realistic heat exchanger configurations comprising a large number of tubes exhibiting different stages of boiling, a coupled CFD simulation of the system as a whole is too difficult because of convergence and stability problems, along with the huge computational overhead. The problem is further complicated by the fact that in many heat exchanger applications, a coupled treatment of the heat conduction in the tube wall, i.e. the conjugate heat transfer modeling, is used. Obviously, such calculations would also demand very long computational times. To circumvent these problems, and enable a time-efficient, but still realistic and useful simulation of the conjugate heat transfer in the tubes of heat exchangers with flow boiling, a novel simplified numerical model is proposed in this study. A significant feature of the model is one-dimensional treatment of the flow side, coupled with two-dimensional treatment of the conduction in the tube wall. On the flow side, 1D modeling of the boiling heat transfer is achieved, using empirical correlation. Because correlation and the reduction in dimensionality are used, the model is numerically much more stable than a CFD simulation, and computationally much more time-efficient. With these features, the model itself tends to be a convenient alternative for the modeling of large systems of tubes. Note that, the model is general in that boiling heat transfer and single phase convection can also be treated by the methodology. However, in this study, the emphasis is placed on boiling heat transfer.
In this study, the computational procedure is demonstrated using the example of water boiling in a vertical tube, with a uniformly heated outer surface. The external fins along the tube wall increase the heat transfer capability.
As already mentioned above, the proposed computational procedure is based on the empirical correlations both for the two-phase pressure drop (i.e. the Friedel model [42][43][44] for frictional pressure drop) and the heat transfer coefficient (Steiner-Taborek model [43,45]). These correlations, validated in the past over an extensive data range, coupled with a simple 1D homogeneous mixture model, allow the determination of flow parameters such as pressure, enthalpy, steam quality, temperature and mass velocity.
The fluid temperature distribution is determined by solving the one-dimensional transient energy equation coupled with the pressure drop formula. The two-dimensional axisymmetric heat conduction equation is solved to obtain the transient temperature distribution in the tube wall. The system of governing equations is solved using the Finite Volume Method [46][47][48][49]. The heat transfer between the solid and fluid domains is conjugated at the fluid-solid interface. At this location, the heat flux is determined as a product of the heat transfer coefficient from the fluid to a solid wall, and the local temperature difference between the wall and bulk fluid temperatures. Because the mass, momentum and energy transport equations for two-phase flow are replaced by the one-dimensional model and the empirical correlation; the computational time is much shorter than in the case of complex CFD simulations. Moreover, the results obtained can then be used as the starting values for more advanced multidimensional modeling.
2 Empirical correlation used for the heat transfer coefficient and pressure drop for flow boiling in vertical tubes Empirical correlation for the heat transfer coefficient and pressure drop are used to determine the heat flux at the fluid solid interface and the pressure level for the flow boiling in vertical round tubes ( Fig. 1).
At the flow inlet, single phase (water) flow is assumed. If the fluid flowing through the tube is heated with a uniform heat flux at a certain distance from the tube inlet, it is possible for the water temperature to reach the saturation temperature. Starting from this location the flow boiling process begins and progresses through the length of the tube. Hence, the proper correlation of the two-phase flow heat transfer coefficient and pressure drop must be applied.

Heat transfer coefficient calculation for the flow boiling process
Conjugate heat transfer occurs at the fluid-solid interface. This study considers the transient heat transfer analysis, therefore all process parameters vary with time s. The heat transfer rate for the two-phase flow depends strongly on the following parameters: the tube inner diameter d i , mass flow rate m p , bulk temperature of fluid T ? (z, s), steam quality X, surface tension r and wall temperature T(r i , z, s). The heat flux at the inner surface of the tube wall q i is given by the following equation: where h is the fluid side heat transfer coefficient (water, wet steam or dry steam), which can be obtained from the following formulae [43,45]: The symbol h sp,l , which appears in Eq. (2) refers to the heat transfer coefficient for single-phase flow, calculated according to the Gnielinski formula [43]. The Reynolds and Prandtl numbers of liquid are denoted as Re l and Pr l respectively, and k l is thermal conductivity of the liquid. The Darcy-Waisbach friction coefficient f D,l for the liquid is defined as: where the R p is the height of the roughness. The Steiner-Taborek model covers the range of R p from 0.1 to 18 lm [43]. In the numerical computations presented in this study, it is assumed that R p = 18 lm. The heat transfer coefficient for the two-phase flow h tp is obtained using the asymptotic Steiner-Taborek model [the second expression of Eq. (2)] where: h l is the local liquid-phase forced convection coefficient based on the total flow for the liquid obtained using the Gnielinski correlation, h nb is the local nucleate pool boiling coefficient at the reference heat flux q 0 and at the reduced pressure p r = 0.1, which is the ratio of the fluid pressure and the critical pressure (p r = p/p crit ), F nb is the nucleate boiling correction factor, F tp is the twophase flow correction factor that accounts for the enhancement of liquid convection resulting from the higher velocity of the two-phase flow as compared with the singlephase flow of liquid in the channel. The asymptotic Steiner-Taborek model assumes [45], that convective and nucleate boiling occurs if q i [ q onb . The critical heat flux for onset of nucleate boiling q onb is defined as: where r denotes the surface tension, T sat is the saturation temperature and r b is the critical nucleation radius, which is assumed to be 0.3 9 10 -6 m. Specific enthalpy and density of vapor are denoted as i g and q g , respectively and i l is the specific enthalpy of saturated liquid. If q i \ q onb only convective boiling occurs. More detailed description of the asymptotic Steiner-Taborek model can be found in [43,45]. For flow boiling processes occurring in a vertical tube, the model presented in this study assumes, that the critical steam quality X cr is equal to 0.5. It indicates the value of X above which the vapor phase starts to dominate in the heat transfer process. For X [ X cr it is assumed, that the heat transfer coefficient for the two-phase flow h(z) is equal to the heat transfer coefficient h sp,g of the vapor phase. The Darcy-Waisbach friction coefficient for vapor f D,g can be obtained from Eq. (3) by replacing Re l with Re g .
The properties for the vapor and liquid phases are determined for each control volume of fluid as a function of pressure and specific enthalpy. The procedure for determining these two variables is presented in the following subsections.

Pressure drop for two-phase flow in vertical tubes
The separated flow model [43] is used to determine the two-phase pressure drop. Three contributions of pressure drop Dp are considered: the hydrostatic pressure drop Dp hstatic , the momentum pressure drop Dp mom and the frictional pressure drop Dp frict : The hydrostatic pressure drop can be calculated as: where Dz is the vertical distance between the inlet and outlet of the control volume (for equally spaced grid points Dz = L/(N -1)), L is the tube length, g is gravitational acceleration and q tp is the estimated density of the twophase mixture, which can be obtained using the following formula: The separated flow model assumes that the two phases are artificially separated into two streams, each flowing inside an individual tube. The fraction of the channel cross-sectional area that is occupied by the gas phase is the so-called void fraction e. According to the Rouhani and Axelsson model [43,50], for e [ 0.1, the void fraction can be obtained from the following expression: The momentum pressure drop Dp mom can be calculated as follows [43]: where out and in refer to the outlet and inlet of the control volume respectively, and G = m p /A is the mass velocity defined as the ratio of the mass flow rate and cross-sectional area of the fluid flow. The frictional pressure drop Dp frict is determined under the assumption that within the zone occupied by a phase, the velocities of each phase are constant across the crosssection. As mentioned previously, the Friedel model, which is applicable when the ratio of dynamic viscosities of the liquid and vapor (l l/ l g ) is less than 1,000 [43], is used to predict the frictional pressure drop for two-phase flow. This condition for the dynamic viscosity ratio of the liquid and vapor phases is satisfied for the computational cases in this study.
The Friedel model correlates the two-phase frictional pressure drop Dp frict with the pressure drop within the liquid phase Dp l as follows: using the multiplier U Fr 2 defined as: where the dimensionless parameters E, F, H and the Froude Fr H and the Weber We L numbers are expressed as [42]: The homogenous density of the two-phase mixture q H is calculated as: The frictional pressure drop for the liquid phase is obtained from: If X = 0, then Dp frict = Dp l .

Governing equations
This study presents the transient analysis of the heat transfer processes occurring in a vertical tube with external fins. External surfaces are often used in heat exchangers to increase the heat transfer coefficient from the gas side [51][52][53][54][55][56]. It is assumed that the flow boiling process occurs inside the tube, which has an inner diameter of d i = 2r i and wall thickness of t = r or i . The part of the computational domain used in the heat transfer analyses is presented in Fig. 2. The fin, w wide and H f = r fr o high, is fixed to the outer surface of the tube wall. The spacing between two consecutive fins is denoted as s. The axisymmetric heat transfer model is applied to the solid domain (tube wall). The solid domain consists of Finite Volumes Dz long and Dr wide. A constant heat flux q is applied to the outer surface of the tube wall (C 1 ). The fluid domain is subdivided into N Control Volumes, located along the length of the tube. The spacing between two consecutive cell centers is denoted Dz. The onedimensional heat transfer model is applied to the fluid domain. It is assumed, that the single phase (water) flow occurs at the tube inlet, the mass flow rate of water is denoted m p , and the inlet temperature and pressure are denoted T in and p in , respectively. The heat transfer between the two domains occurs at the Fluid-Solid interface (C 2 ). The computations are carried out considering the temperature dependent fluid properties (density, thermal conductivity, viscosity and surface tension) but assuming constant thermal properties for the tube wall material (density, specific capacity and thermal conductivity). It is assumed that the initial temperatures of the fluid and tube wall are equal.

Energy equation for the fluid domain
The transport equation for the static fluid enthalpy i can be given as [57]: along with an assumed initial temperature and pressure field: The pressure and enthalpy at the inlet are given by: The initial distribution of the heat transfer coefficient is determined using Eq. (2) for a given temperature and pressure: T 1 ðz; s ¼ 0Þ and p(z, s = 0), respectively.

Heat conduction equation for the solid domain
The axisymmetric heat conduction equation in a cylindrical coordinate frame is given by [49]: The heat flux in the j direction is a function of temperature gradient and can be modeled according to Fourier's law: q j ¼ Àk oT ox j . In the cylindrical coordinate frame x j (j = 1 or 2) corresponds to the r and z coordinates respectively. The boundary conditions are assumed as: where n j denotes the outward normal direction cosines. The initial temperature field is defined as: 4 Computational procedure The initial temperature distribution in the fluid and solid domains is assumed according to Eqs. (16) and (23), respectively. The initial pressure field is given by Eq. (17). Equation (20) is discretized spatially for the solid domain using the Finite Volume Method, which results in a system of ordinary differential equations (ODE): A detailed description of the discretization of the axisymmetric heat condition equation for the solid domain using the Finite Volume Method can be found in [49]. Equation (24) is discretized in time and solved at consecutive time instances. Using the common ''Finite Element'' nomenclature, in Eq. (24), M denotes a diagonal global capacitance matrix, dT ds is the vector of the temporal derivatives of the nodal temperatures, K is a global thermal stiffness matrix, which is the sum of the conductivity matrix K cond and the diagonal matrix of convective loads K conv : T denotes the nodal temperatures vector, F refers to the thermal load or forcing vector, which is the sum of the heat flux load vector F q and convective load vector F conv : To solve the ODE system (24) for nodal temperatures T sþDs at time instance s ? Ds, the Crank-Nicolson time stepping scheme is applied: Assuming that Ds approaches zero, the following simplification can be made: Applying the central differencing scheme for dT ds sþDs=2 yields: dT ds The nodal temperatures vector T sþDs=2 is approximated as: Substituting Eqs. (28)-(31) into Eq. (27), and rearranging, the following formula for T sþDs is obtained: where: If the temperature distribution in the solid domain T sþDs is calculated, the temperature field in the fluid domain is solved using Eq. (15). The discretization of Eq. (15) is discussed in the following section. The conductivity term can be ignored in Eq. (15), resulting in: At time instance s ? Ds the enthalpy at the z location of the jth node i j s?Ds can be determined by employing the forward time stepping scheme for temporal derivatives and backward differencing scheme for spatial derivatives: i sþDs The value of the pressure at the z location of the jth node is calculated as follows: The fluid temperature T j,?
s and density q j s are determined using steam tables XSteam v2.6 [58] for the previously determined p j s and i j s : Heat Mass Transfer (2015) 51:553-566 559 The steam quality can be obtained using the following formula: The heat flux q j s at the location of the fluid-solid interface is calculated as: where T j s denotes the tube wall temperature at the corresponding z location of the jth node, determined using Eq. (32). This temperature is calculated at the location of the fluid-solid interface (C 2 , Fig. 2).
The explicit method for solving the energy equation for the fluid domain, Eq. (35), is valid only for a CFL number Dz ) less than 1. Therefore, for all control volumes, the following stability criterion for time step Ds must be satisfied [57]: where v j is the velocity calculated at the corresponding z coordinate of the cell node: 5 Results and discussion The computational procedure, described in Sect. 4, was implemented in numerical code written in MATLAB r2012b [59]. The transient and steady-state computations were performed using the code to obtain the 1D transient temperature field of the flow, and 2D transient temperature field in the wall material. The results are presented in the following form: • Transient temperature variations (fluid and solid temperatures) at specified distances measured from the inlet plane, • Steady state distributions of temperature, pressure, steam quality and heat transfer coefficient within the tube wall.
The influence of the mass flow rate and heat flux at the outer surface of tube wall on the transient temperature distribution in the fluid and solid wall was studied first. The computations were performed for the settings defined in Table 1.   In Fig. 3a, for the low mass flow rate m p = 0.03 kg/s, it is observed that flow boiling occurs at z = 0.2 m. The corresponding time for the start of the boiling process is s b = 37 s. This results in a rapid decrease in the wall temperature (see Fig. 3b) for: the inner and outer surfaces of the tube wall, and the fin tip. The largest temperature drop is observed for the inner tube surface (from 217 to 178°C), the smallest temperature drop occurred for the fin tip (from 262 to 238°C).  15) and (39) we can deduce that the large mass flow rates lower the rate of enthalpy transport along the flow direction. This conclusion seems to be obvious because if the flow is faster, its temperature increase along the flow direction must be lower, because the tube wall is cooled down intensively. On the other hand, it should be noted that the single phase flow heat transfer coefficient [see Gnielinski correlation-Eq. (2)] increases with an increase in the mass flow rate. In turn, the wall-to-liquid temperature difference (T wall -T ? ) increases faster with time for low mass flow rates (i.e. m p = 0.03 kg/s) than for the higher flow rates (i.e. m p = 0.3 kg/s)- Fig. 3b-d. All these factors directly influence the heat flux value at the fluid-solid interface, affecting the shape of the transient temperature variations for the fluid (Fig. 3a). The rapid increase in the heat transfer coefficient, associated with the flow transition from single-phase to two-phase types, causes the wall temperature to drop (especially at the inner tube surface). This can be observed for z = 1 m (Fig. 3c) and 1.8 m (Fig. 3d) (Fig. 3c, d).
Besides the mass flow rate, the other parameter which influences the boiling process in the analyzed fin and tube system is the heat flux at the outer surface of the tube wall. In Fig. 3e we observe that for the constant mass flow rate m p = 0.3 kg/s, the flow boiling starts faster, if the heat flux q is large. For the flow cross section at z = 0.2 m, the single phase flow occurs for all the analyzed values of q: q = 7,500 W/m 2 , q = 15,000 W/m 2 and q = 25,000 W/ m 2 . The fluid temperature at these locations is lower than the temperature of the saturated liquid for the whole period of time. The wall temperature variations corresponding to the location z = 0.2 m are depicted in Fig. 3f. Because the phase change process does not occur at this location, the wall temperature increases with time, and reaches steady state after 130 s. Moving in the flow direction, for z = 1 m and q = 25,000 W/m 2 , we observe that flow boiling occurs earlier (s b = 22 s) compared to the case q = 15,000 W/m 2 with s b = 54 s (Fig. 3e). The rapid change in the heat transfer mechanism by boiling causes a substantial drop in temperature in time (DT), from its value at the onset of boiling to its value at steady state, as discussed above. These temperature drops corresponding to the rapid change in the heat transfer mechanism are for q = 15,000 W/m 2 DT(r i , z = 1 m) = 6°C, DT(r o , z = 1 m) = 5°C, DT(r f , z = 1 m) = 3°C (Fig. 3g).
If q = 25,000 W/m 2 then the fin-tip temperature does not decrease during the analyzed time period. At the flow cross-section, where z = 1.8 m, fluid starts to boil when s b = 12.5 s for q = 25,000 W/m 2 , when s b = 19.5 s for q = 15,000 W/m 2 , and when s b = 62 s for q = 7,500 W/ m 2 . For z = 1.8 m, the temperature decrease of the fin tip and the tube outer surface is not observed if s [ s b when q = 7500 W/m 2 and q = 25,000 W/m 2 (Fig. 3h). The slight decrease in the temperature of the inner surface of the tube wall DT(r i , z = 1.8 m) = 1.5°C is observed for all analyzed values of q if s [ s b . Figure 4 presents the steady state distributions of steam quality, heat transfer coefficient and pressure along the tube length for different values of heat flux q and mass flow rate m p . For a constant heat fluxq = 20,000 W/ m 2 and variable mass flow rates: m p = 0.03 kg/s, m p = 0.1 kg/s, and m p = 0.3 kg/s, it can be observed that the z location of the cross section where flow boiling begins increases with mass flow rate (Fig. 4a). For m p = 0.03, 0.1 and 0.3 kg/s, the boiling onset positions are at z = 0.1, 0.25 and 0.68 m, respectively (see Fig. 4a). The phase change process results in a significant increase in the heat transfer coefficient (typical for the nucleate boiling process): from 800 W/(m 2 K) to 19,000 W/(m 2 K) for m p = 0.03 kg/s, from 1,900 W/(m 2 K) to 15,200 W/ (m 2 K), for m p = 0.1 kg/s, and from 4,800 W/(m 2 K) to 13,600 W/(m 2 K) for m p = 0.3 kg/s (Fig. 4b). Therefore, a large spatial temperature difference (DT s ) between the  Fig. 4d). The nucleate boiling occurs for q = 25,000 and 15,000 W/m 2 . The heat transfer coefficient increases rapidly from 4,800 W/(m 2 K) to 15,850 W/(m 2 K) if q = 25,000 W/m 2 and to 11,150 W/(m 2 K) if q = 15,000 W/m 2 . The mass velocity of steam increases with the flow direction, which results in a further increase in the heat transfer coefficient. The nucleate boiling does not occur in the case of q = 7,500 W/m 2 . The slight increase in the heat transfer coefficient from 4,800 to 6,500 W/(m 2 K), which is typical for convective boiling, is observed for z = 1.65 m to z = 2 m.
The pressure variations along the flow direction (Fig. 4c, f) slightly affect the local saturation temperature (referred to the corresponding z coordinate of the flow cross-section), which can be observed in Fig. 5. The hydrostatic pressure drop is the largest portion of the

Summary
This study presents a novel, simplified mathematical model for the simulation of conjugate flow boiling heat transfer in vertical round tubes with external fins. The model uses the experimental correlation developed by Steiner and Taborek, which allows the heat transfer coefficient for flow boiling to be determined. For two-phase flow, the pressure level was determined using empirical correlation (i.e. the Friedel model for the frictional pressure drop, the Rouhani and Axelson model for the momentum pressure drop, and the separated flow model for the hydrostatic pressure drop). The one-dimensional heat transfer model was applied to the fluid domain, while the axisymmetric heat conduction model was used for the tube wall. The Finite Volume Method was employed to discretize the heat transfer equations for fluid and solid domains. Transient heat transfer analyses were performed. The combination of the empirical correlations and simplified (1D/2D) mathematical models for the heat transfer processes occurring in the fluid domain and the tube wall enable efficient simulations of the conjugate boiling processes. The developed model is simplified but very timeefficient compared to the CFD approach. Thus it can conveniently be used in the analysis of large and complex tube systems that are frequently encountered in practical engineering applications. If required, for a more detailed analysis of the boiling process in certain tubes or parts of the system, a CFD analysis can be performed for an isolated region, using the values delivered by the present model as initial and boundary conditions. Future work will include the experimental validation of the present methodology as well as comparisons with CFD simulations for isolated tube parts.
Open Access This article is distributed under the terms of the Creative Commons Attribution License which permits any use, distribution, and reproduction in any medium, provided the original author(s) and the source are credited.