Heat Transfer in Evaporator of Thermal Sink in Presence of Subcooled Boiling Section

One of the most important tasks of spacecraft design is providing their thermal regime. Currently, there is a tendency to use in spacecraft with heat generation over 5 kW thermal control systems (TCS) based on mechanically pumped two-phase loops (MPTL) [1-3].

MPTL have several advantages over single-phase thermal control systems. In particular, MPTL can transfer much more heat per unit of mass flow rate than with single-phase heat transfer systems. Energy consumption of pump for working fluid pumping is low. Heat transfer by means of boiling allows maintaining temperature of objects at almost entire length of the loop close to saturation temperature. Heat transfer processes where aggregate state of substance changes take place (boiling, condensing) are much more intense than those where convective heat transfer in a single-phase liquid occurs. These features of MPTL result in a significantly lower mass of thermal control system compared to single-phase system. Therefore, the task of creating TCS on the basis of MPTL is actively solved all over the world. For example, this task is included in the program of NASA prospective works [4] for 2015-2035.

The main element of MPTL heat acquisition subsystem is a thermal sink - a contact heat exchanger, on surface of which the cooled components are fixed. The heat from components is drawn to the thermal sink and then to the thermal fluid, which is pumped through the channel evaporator. There is a section of "subcooled boiling" in case of subcooled liquid at evaporator inlet, where heat is transferred by surface boiling of fluid. At present, however, the calculation of thermal sinks for both space and ground applications has not been sufficiently worked out. Particularly this is the case of thermal sink, where subcooled single-phase liquid at inlet and two-phase flow at outlet. In this case, depending on heat flux density and parameters of working fluid, there may be a section of "subcooled boiling" - a section where the heat transfer takes place due to surface boiling in subcooled liquid.

The research aims to develop a model for calculating heat transfer and average evaporator wall temperature in presence of subcooled boiling section in a wide range of parameters specific for thermal control systems of space vehicles. In particular:

- ammonia saturation temperature - from sub-zero temperatures to +75℃; 

- liquid ammonia subcooling at the evaporator inlet - from 0 to 30 K;

- mass velocity - from 27 to 200 kg/(s*m²);

- average heat flux density on the evaporator wall - from 0 to 18 W/cm².

Problems of heat transfer at forced convection [5, 6] and saturated boiling in channels [7-9] are well investigated, numerous correlations and methods are developed for calculation of local heat transfer coefficients and evaporator wall temperature. However, calculating the heat transfer during surface boiling in subcooled liquid still have certain difficulties at present. Performed experimental research is mainly related to boiling of water. Data on ammonia boiling are limited [10]. To determine the heat transfer intensity in the section of subcooled boiling it is necessary to perform additional experimental studies in the range of parameters specific for spacecraft thermal control system.

Figure 1 shows the modern outlooks of heat transfer mechanism of liquid in uniformly heated channel in a presence of subcooled boiling [11].

Heat transfer in zone A is determined by forced convection of liquid.  The convective heat transfer coefficient of liquid is almost constant (without an account of properties change with temperature), wall temperature increases linearly and parallel to liquid temperature. The first bubbles appear on the wall at ONB point (beginning of zone B), which is identified as beginning of surface boiling. There is usually a temperature jump associated with a sharp increase in heat transfer intensity in this cross-section. The contribution of boiling to heat transfer is higher with flow downstream following, while the single-phase convective heat transfer contribution decreases. It's a partial subcooled boiling interval. At the FDB point the convective contribution to the total heat transfer becomes insignificant and a fully developed subcooled boiling is established (zone C₁). Following on the flow, the average temperature of liquid reaches saturation temperature and developed boiling take place in already saturated liquid (zone С₂). The wall temperature remains almost constant in zones С₁ and С₂.

1.png

Figure 1. Distribution of wall and liquid temperature, mass vapor quality along the length of a uniformly heated channel

In contrast to free convection (Pool boiling), forced convection affects the heat transfer in addition to the boiling heat transfer when liquid flow in the channel [12, 13]. The combined effect of these two mechanisms of heat transfer takes place in all regions of the heated channel: subcooled boiling and boiling in saturated liquid. The contribution of forced convection to the total heat transfer is mainly determined by the flow rate: it could be close to zero or almost completely determines the heat transfer intensity. In the subcooled boiling section, these two mechanisms are always compatible: at the beginning of the section, the heat transfer is fully determined by convection, at the end of the section, the boiling contribution becomes the determining factor [14, 15].

Numerous experiments show that in channels at developed subcooled boiling or boiling in saturated liquid up to the volume vapor quality of the flow ~ 0.7 the contribution of the convective component can usually be neglected [5]. The heat transfer coefficient in this section does not depend either on subcooling or on the mass flow rate of the working fluid. It can be calculated as a developed saturated boiling (Pool boiling) [16, 17].

Various models have been proposed to calculate the heat transfer for the forced flow where the boiling take place. The paper [18] provides a review of the most common models. In all models, the flow region is divided into characteristic heat transfer zones: single-phase convection, subcooled boiling and developed boiling. The boundaries of each area are set differently in different models. Calculation of heat transfer in the zone of subcooled boiling is based on the superposition principle. For example, in Labuntsov model [5], the whole range of possible values of heat transfer coefficients is divided into three zones distinguished by different contribution of heat transfer convective component. In this case dimensionless coordinates are used h_q/h_L.

The author identifies the following zones:

In zone 1, h_q/h_L <0.5, h_TP = h_L – boiling does not affect heat transfer. The heat transfer intensity is fully determined by the forced convection of the liquid.

In zone 2, h_q/h_L> 2, h_TP = h_q – the heat transfer intensity is completely determined by developed pool boiling; forced convection does not affect heat transfer. 

In zone 3, 0.5 ≤ h_q/h_L ≤ 2 – transient boiling zone with the mutual influence of convective heat transfer and boiling heat transfer; the expression for the heat transfer coefficient is described by the interpolation formula:

$h_{T P}=h_{L} \cdot\left(4 \cdot h_{L}+h_{\varphi}\right) /\left(5 \cdot h_{L}+h_{q}\right)$     (1)

Also note that most models [8, 19] consider local heat transfer coefficients and wall temperatures along the length of a uniformly heated channel. However, in engineering tasks of thermal state calculation of thermal sink, especially thermal sink with high longitudinal thermal conductivity, it is of interest not to describe in detail the wall temperature changes along the length of the heated channel, but to determine the average or maximum wall temperature of the heat transfer device. For this purpose, it is rational to use the concept of average heat transfer coefficient along the perimeter and length of evaporator and the simple ratios for calculation of average wall temperature T_w.

The paper proposes a model for calculating the average wall temperature of a thermal sink with large longitudinal thermal conductivity. The model makes it possible to calculate the intensity of heat transfer in each zone. At the same time, the boundaries of each zone are justified experimentally for the parameters, typical for TCS of spacecrafts. Ammonia is used as a working fluid. Adequacy of the model is confirmed by experimental results.

Ammonia with a given pressure, mass flow rate and temperature were supplied to the thermal sink inlet in experiment.

 The experiments are performed in the following range of parameters:

• saturation temperature T_sat = 35 … 75℃;

• inlet fluid subcooling ∆T_sub = 0…30℃;

• mass flow rate m = 0.2…7.5 g/s;

• mass velocity: G = 5…200 kg/(s*m²);

• average heat flux density q = 0…18 W/cm².

The experiments are performed in three thermal sink orientation: two horizontal (heater installed at the top or bottom) and one vertical (upstream).

The parameters were measured in the experiments:

• inlet pressure of the thermal sink р_HCA;

• mass flow rate m;

• heater power Q;

• inlet and outlet thermal sink temperatures T_in, T_ex;

• temperatures on thermal sink surface (9 sensors) T_w1...T_w9.

Instrumental measurement errors of parameters:

- absolute pressure ± 0.2 %;

- mass flow rate of liquid ± 0.2%;

- temperature: ± 0.15°C (0℃), ± 0.35℃ (100℃).

To avoid heat loss, the thermal sink was insulated. Heat loss was not more than 3% and was not taken into account under the analysis of results.

The thermal state of the thermal sink was analyzed using SolidWorks Flow Simulation tool for different heater powers. The analysis showed that in the fixed cross-section of the thermal sink the average surface temperature D is almost equal to the average measured temperature of three sensors installed in points 1, 2, 3 on the profile surface (Figure 2). At the same time, the analysis showed high longitudinal heat conductivity of the structure: the average temperature difference in three cross-sections did not exceed 1K. Therefore, the average temperature of the heat transfer surface along the perimeter and length of the T_w evaporator was determined by averaging the readings of 9 sensors on the profile of the thermal sink. With that, the methodological error is about 1K.

Based on measured results next parameters were calculated: 

• average heat flux density over surface "D":

$q=Q / F_{\text {ev }}, \mathrm{W} / \mathrm{cm}^{2}$;     (2)

• saturation temperature T_sat as a function of pressure;

• liquid subcooling at the evaporator inlet ∆Т_sub:

$\Delta T_{\text {suk }}=T_{\text {sat }}-T_{\text {in }}, \mathrm{K}$;    (3)

• vapor quality х_ex at the evaporator outlet:

$x_{e x}=\left(Q / m-C_{p} \cdot \Delta T_{s u b}\right) / r$,     (4)

where, Cp – is heat capacity at constant pressure of the liquid working fluid, J/kg/K, r – latent heat of vaporization, J/kg.

If the balance vapor quality х_ex was negative, the balance temperature of liquid at the outlet channel was also calculated:

$T_{e x . b}=T_{i n}+Q /\left(C_{p} \cdot m\right)$    (5)

The heat transfer equation was used to find an effective average heat transfer coefficient $\mathrm{h}_{\mathrm{ef}}$. In this case, the effective average temperature of working fluid was specially indicated $\mathrm{T}_{\mathrm{ef}}$:

$h_{e f}=q /\left(T_{w}-T_{e f}\right)$    (6)

The heat transfer model in evaporator of thermal sink for specific heat flux density q varies from zero to the maximum value, which correspond to the developed saturated boiling mode, is given below.

The model is illustrated in Figure 7. The full range of heat transfer in the thermal sink evaporator T_w=f(q) s divided into three zones: А, В, С. There are calculated boiling curve and heat transfer zones are applied on the experimental results on figure. The calculated balance vapor quality x_ex is also given. As can be seen, there is a liquid at thermal sink outlet up to q_x0 ≈ 7.1 W/m². The section of the boiling curve from q_a to q_x0 can be conventionally called as section of "subcooled boiling".

7.png

Figure 7. Heat transfer zones

The model includes determination of boundaries of zones A, B, C, effective average temperature value of working fluid T_ef and sequential computation of heat transfer in each zone (calculation of T_w and h_ef).

As input data for model are specified:

• T_sat – saturation temperature in thermal sink evaporator, ℃;

• ∆T_sub – liquid subcooling at the evaporator inlet, К;

• m – mass flow rate, kg/s;

• F_ev – evaporator heat transfer surface, m²;

• d,L− diameter and length of the evaporator, m;

• Q – thermal sink heat load, W.

Calculated parameters:

• T_w – average heat transfer surface temperature, ℃;

• T_ef – effective average temperature of the working fluid, ℃;

• h_ef – calculated heat transfer coefficient, W/(m²·K);

• х_ex – mass vapor quality at the evaporator outlet.

• average heat flux density q (Eq. (2)), W/m²;

• evaporator inlet temperature T_in (Eq. (3)), ℃;

• evaporator outlet temperature T_ex (Eq. (5)), ℃.

Hereinafter, the index "calc" was omitted.

Note: calculation from the given input data allows to determine only T_w and х_ex. To find the heat transfer coefficient h_ef according to the Eq. (6), you first need to determine the effective temperature of the coolant T_ef. In different zones of the boiling curve, it is determined in different ways.

Since it is not possible to determine the heat transfer zone at once at a given value Q, a sequential calculation of boundaries and heat transfer in all three zones of the boiling curve is implemented.

Zone A - convective single-phase heat transfer, linear dependence T_w = f(q)

 The calculation is performed according to traditional methods for heat transfer of single-phase liquid in channels [20, 21]. If some simplifying assumptions are accepted, then can be written down:

• effective average temperature of the working fluid:

$T_{\text {ef. } A}=\left(T_{\text {in }}+T_{\text {ex }}\right) / 2=T_{\text {in }}+\left(q \cdot F_{\text {ev }}\right) /\left(2 \cdot C_{p} \cdot m\right)$     (7)

• thermophysical properties of a liquid are considered as constant and correspond to Т_ef.А;

• heat transfer coefficient: $\mathrm{h}_{\mathrm{L}, \mathrm{A}}=\mathrm{Nu} \cdot \lambda_{\mathrm{L}} / \mathrm{d}$, where Nu is determined according to the flow pattern by known correlations for single-phase flow [20, 21];

the average wall temperature in zone "A" can be determined from the equation:

$T_{w \cdot A}=T_{i n}+q \cdot\left(F_{e v} /\left(2 \cdot C_{p} \cdot m\right)+1 / h_{L}\right)$     (8)

Essentially, this equation is a straight line Т_w = f(q), as the heat transfer coefficient in this zone changes slightly, what is confirmed by the experiment.

The boundary of zone "A" is defined from the condition: T_w.a = T_sat.

At the border (at point "a"):

$T_{w \cdot a}=T_{\text {sat }}$.    (9)

$q_{a}=\left(T_{w . a}-T_{i n}\right) /\left(F_{e x} /\left(2 \cdot C_{p} \cdot m\right)+1 / h_{L, A}\right)$     (10)

Zone C - developed boiling zone (С1 - in subcooled liquid; С2 – in saturated liquid)

The effective heat transfer coefficient in zone C can be calculated by correlations for developed boiling in channels with correction for flow velocity (mass flow rate) [9]. The influence of flow velocity on heat transfer in this zone up to the volume vapor content at the outlet ~ 0.8 cannot be taken into account and it is acceptable to use heat transfer coefficients for pool boiling. 

The heat transfer coefficient for pool boiling is proposed to be calculated by Kupriyanova formula [10].

$h_{e f . C}=h_{q \cdot C}=2.2 \cdot p^{0.21} \cdot q^{0.7}$     (11)

In spite of the fact that this correlation was obtained for boiling on the surface of a horizontal tube in the temperature range from -40…+20℃, and our experiments were carried out in the temperature range of +35...+75℃, comparison with the experiment showed that the formula quite well describes all performed series of experiments. The difference between experimental and calculated T_w, values did not exceed 2 K.

The effective temperature in zone "C" is:

$T_{e f . C}=T_{\text {sat }}$     (12)

The average wall temperature in zone C at any heat flux density is determined from the heat transfer equation (Eq. (6)):

$T_{w . C}=T_{\text {ef.C }}+q / h_{\text {ef.C }}$    (13)

The top border of zone C q_max is limited by volume vapor quality ~ 0.7...0.8.  Nevertheless, as experimental results show, satisfactory matching of the theory and experiment was observed up to х_ex ≈ 0.7…0.8.

The lower boundary of the zone q_b is determined by the algorithm below.

Zone B - transition zone of subcooled boiling in part of evaporator channel

The boundaries of zone "B" are indicated by points "a" and "b" on the boiling curve. The parameters in this zone change smoothly without jumps and curve breaks from point "a" to point "b". Therefore, to describe the function T_w =f(q) in the zone "B", it is acceptable to use approximation in the form of a Bezier spline function [22] without derivative discontinuity in the points "a" and "b".

Note: Some outliers of experimental values T_w in zone B (Fig. 3 and 7) are explained by different number of active centers of steam formation on the heat-exchange surface at increase and decrease of specific heat flux. This phenomenon is known and called "heat flux hysteresis" [10]. In experiments with aluminum thermal sink, the hysteresis was small. However, in stainless steel smooth channels for ammonia heat transfer, the hysteresis can be more than 10K.Therefore, our proposed methodology for calculating the parameters in the transition zone B refers to the process of reducing the heat load from the developed boiling zone, when most of the vaporization centers are activated.

 The tangents in the points "a" and "b" should be drawn up to theirs intersect at point "0" in order to trace the Bezier curve. By the three points "a", "b" and "0" the Bezier function (Figure 7, line T_w.b) could be found. Thus, to find the function T_w =f(q) in the zone B, the coordinates of the points "a", "b" and "0" should be available.

The coordinates of point "a" are determined according to Eq. (9) and (10).

To determine the width of the transition zone "B" and the coordinates of point "b" we use Labuntsov approach - use dimensionless coordinates h_q/h_L [5].

Determine the position of points "a", "0" and "b" relative to each other in coordinates h_q/h_L.

The difference in the "a" and "0" points position in the non-dimensional coordinates is equal:

$\Delta\left(h_{q} / h L\right)_{a-0}=\left(h_{q} / h_{L}\right)_{0}-\left(h_{q} / h_{L}\right)_{a}$    (14)

We will estimate the position of point "0" by the results of experimental studies. It is determined that $\Delta\left(\mathrm{h}_{\mathrm{q}} / \mathrm{h}_{\mathrm{L}}\right)_{\mathrm{a}-0}$ changes depending on parameters of the working fluid at the thermal sink inlet in very insignificant limits, from 0.65 to 0.70.

Since the curve T_w.q=f(q) in the area of point "b" is slightly sloping, its location can be assigned approximately by taking the offset of points "a" and "b" in dimensionless coordinates:

$\Delta\left(h_{q} / h_{L}\right)_{a-b} \approx 2 \Delta\left(h_{q} / h_{L}\right)_{a-0} \approx 1.5$      (15)

$\left(h_{q} / h_{L}\right)_{b}=\left(h_{q} / h_{L}\right)_{a}+1.5$      (16)

The further algorithm for calculating the coordinates of point "b" is as follows:

- determine the heat transfer coefficient in point "a" h_Lа by known correlations for single-phase flow;

- knowing the heat flux density q_a at point "a" (Eq. (10)), the heat transfer coefficient at boiling h_qa should be calculated by the Eq. (11).

- find a dimensionless relation $\left(\mathrm{h}_{\mathrm{q}} / \mathrm{h}_{\mathrm{L}}\right)_{\mathrm{a}}$ at point "a";

- calculate a dimensionless relation $\left(\mathrm{h}_{\mathrm{q}} / \mathrm{h}_{\mathrm{L}}\right)_{\mathrm{b}}$ at the point "b" by the Eq. (16);

- determine the heat transfer coefficient of boiling at the point "b" considering that $\mathrm{h}_{\mathrm{L.b}} \approx \mathrm{h}_{\mathrm{L.a}}$:

$h_{q . b}=h_{L . b} \cdot\left(h_{q} / h_{L}\right)_{b}$     (17)

- calculate the heat flux density at point "b" by transforming the Eq. (11):

$q_{h}=\left(h_{q, b} / 2 \cdot 2 \cdot p^{0.21}\right)^{\frac{1}{0.7}}$   (18)

- wall temperature at this point:

$T_{w \cdot b}=T_{s a t}+q_{h} / h_{q . b}$    (19)

Find the derivatives at points "a" and "b": $\alpha=\left(d T_{w . A} / d q\right)_{a}=F_{e v} /\left(2 \cdot C_{p} \cdot m\right)+1 / h_{L}$ (see. Eq. (8));

$\beta=\left(d T_{w . C} / d q\right)_{b}=0.3 / h_{q . b}$ (see. Eq. (13)).

The coordinates of the point "0" should be found as intersection of tangents:

- heat flux density at point "0":

$q_{0}=q_{a}+\left[\left(T_{w . b}-T_{s a t}\right)-\left(q_{b}-q_{a}\right) \cdot \beta\right] /(\alpha-\beta)$

- wall temperature at point "0":

$T_{w .0}=T_{s a t}+\left(\left[\left(T_{w . b}-T_{s a t}\right)-\left(q_{b}-q_{a}\right) \cdot \beta\right] /(\alpha-\beta)\right) \cdot \alpha$

Having coordinates of three points "a", "b" and "0" the Bezier quadratic curve could be drawn [22], by which we determine the temperature of the wall in any point of zone B.

The heat transfer coefficient in zone "B" does not make sense to find, because it is not defined how the T_ef value in this zone changes.