The freezing heat exchanger: Part I – Introduction

Abstract

Here I’m going to show the methods I used to tackle the problem of sizing a freezing heat exchanger. I think this method can be used for many heat exchangers in the industry, and can provide insight on whether the device clogs up, what is the operating point and how the solid layer evolves through time.

1. Introduction

Heat exchangers are ubiquitous. They are present in computers, fridges and all sorts of machines. More times than not, the heat exchanger is not even intended to be there!

Each HX has its own design difficulties and particularities, but the one I’m going to talk about today is a quite unusual one. Sometimes, an industrial process needs to cool down a liquid as close as possible to its freezing point. Or maybe, it just so happens that a given piece of pipeline is carrying a molten metal that should not freeze or the process will stop. Whichever the application, a heat exchanger that works below solidification temperature will always run the risk of clogging up, which normally results in a catastrophic failure. In the case of freezing water HXs, the reduction in density will most certainly break the tube due to astronomic pressures built up inside of it. So we, engineers, don’t want that.

So which tools do we need to predict what’s going to happen? Well, if you’re reading this you should probably already be familiar with standard heat exchanger sizing procedures, such as the LMTD and the ε-NTU methods. If not, I strongly suggest you to review that, so you can appreciate the content I’m going to talk about!

2. State of the art

These two methods aforementioned look at the heat exchanger under a “global” point of view, where the external and internal convection coefficients are constant and the fouling thermal resistance is uniform throughout the HX surfaces. Although simplistic, these assumptions work pretty well for the sizing of most HX’s, and everybody in the field knows the errors are in the order of 25%, so oversizing is a common practice.

But when it comes to freezing exchangers, these assumptions are not quite true, to the point where it might be relevant to consider the heat exchanger in more detail. Why? As the liquid cools down further in the HX, more and more solid forms, thickening the “solid fouling” layer. This impairs the local heat transfer coefficient, but if one assumed a constant thickness to consider as a fouling factor, the HX might be overly oversized. Actually, oversizing a freezing heat exchanger can be even more dangerous, as the risks of clogging increase with the HX surface area.

ice-buildup — Ice buildup in the external surface of a water chiller

Apart from the freezing fouling factor, there is the secondary effect of variable flow speed. Let’s analyze the case of a tubular HX where the freezing liquid runs inside the tube (remember the number one rule of HX design: the fouling fluid must be in the tube!). As the frozen layer thickens, the inner diameter available for the liquid to flow is gradually reduced. This has two main effects: The local speed of the fluid increases, but the pressure drop of the overall HX will increase, too. Depending on the pumping system design, the fluid speed might actually decrease. In any case, the internal convection coefficient, heavily dependent on the flow speed, will be higher at the regions where the tube is most frozen.

Dr. Lee, from Iowa State University, did a great job in reviewing the problem of freezing in pipe flows (check out his PhD thesis here). In summary, the problem of freezing is a combination of the Stefan-Neumann problem and the Graetz problem. The former analyzes how the boundary evolves in a freezing fluid that does not flow, while the latter equates how to come up with a description of the thermal boundary layer. The equations get pretty complex and one can easily “clog up” their brain. No pun intended.

3. Complicating our lives

So the equations I’m talking about are the ones shown below. They are an excerpt from Dr. Lee’s work, and assume 2D internal flow, in radial coordinates. It also assumes the liquid never supercools, and therefore its temperature is never below freezing. The boundary solid-liquid is expressed in the last equation, named “interface”.

As it is rather unpleasant when the variables that comprise the equations are not defined, I’m going to define them below for the sake of completeness:

$r$ is the radial coordinate and $z$ is the axial coordinate. The location $R$ is the radial coordinate of the solid-liquid boundary. $\dot R=\frac{\partial R}{\partial t}$ is the time derivative of the boundary coordinate, while $R^{\prime}=\frac{\partial R}{\partial z}$ is its axial derivative.

$T_s$ and $T_l$ are, respectively, the local solid and liquid temperatures, $k_s$ and $k_l$ are their thermal conductivities, $\alpha _s$ and $\alpha _l$ are their thermal diffusivities and $h_{sl}$ is the latent heat of solidification.

Finally, $t$ is the time coordinate and $u$ and $w$ are the fluid velocities in the $z$ and $r$ directions, respectively. The local pressure $P$ , liquid density $\rho$ and liquid kinematic viscosity $\nu$ are needed to solve the flow momentum equations.

So now that we’ve defined all the variables, I can talk about what we can do to make the problem simpler. Of course, one can completely solve the six equations simultaneously in a CFD package, but apart from expensive, time-consuming and computationally difficult, these solutions will not give much insight about the problem. So here I’ll present an alternative, much simpler solution, that might not be “mathematically accurate” but will probably be sufficient for an engineering application.

4. Simplifying our lives

As always, be aware of the simplifications and assumptions and evaluate whether they are true in your application. So, let’s first state them:

The equations presented will only work for a tubular HX where the internal fluid is the one which freezes.
The thermal boundary layer will fully develop in a small fraction of the HX tube length and will change slowly as the fluid cools down.
The ice thickness gradient in the axial direction is very small.
The fluid will not supercool. This means the wall temperature will never drop below freezing temperature. Either the wall is frozen ( $T_{l, @r=R}=T_f$ ) or above freezing temperature ( $T_l>T_f$ ).

These very reasonable assumptions for the engineering application of a shell-and-tube heat exchanger will basically eliminate the need to solve 5 out of 6 of the differential equations shown before. How?

If the third assumption holds, we can simplify the interface equation by dropping the $R^{\prime}$ term. This will make it look like the following:

$\rho_s h_{sl} R^{\prime}=k_s \frac{\partial T_s}{\partial r}-k_l \frac{\partial T_l}{\partial r}$

Let’s now analyze what this equation means. It is equivalent to an energy balance in terms of heat flux. We could replace it by the equation shown below, where each term corresponds to a heat flux. Basically, what it’s saying is that there must be an imbalance between the overall heat transferred through all the thermal resistances ( $q_s''$ ) and the heat absorbed by the liquid ( $q_l''$ ). And this imbalance, as energy is conserved, must go somewhere. It will go to freeze some of the liquid ( $q_{sl}''$ ).

$q_{sl}''=q_{s}''-q_{l}''$

qsl — Simple diagram showing the heat flux conservation

What is beautiful about this equation the way it is stated is that $q_{s}''$ is only dependent on the five typical thermal resistances of a heat exchanger, while $q_{l}''$ is dependent only on the convection of the fluid size, if the tube wall temperature is smaller than zero. In mathematical words:

$q_{s}''=(T_l-T_{freeze})(\frac{1}{h_{ext}}+\frac{t_{frozen}}{k_{s}}+\frac{t_{wall}}{k_{wall}})^{-1}$

and

$q_{l}''=\frac{R_{ext,tube} h_{int} }{R}(T_l-T_{freeze})$

(These equations above only hold if the wall is already frozen. Otherwise, the correct equation to be used is the regular heat exchanger equation.)

If the second assumption holds, we can use the correlation of our choice to solve for finding $h_{int}$ , such as the Sieder and Tate correlation:

$Nu_D=0.027 Re_D^{0.8} Pr^{1/3} (\frac{\mu}{\mu_s})^{0.14}$

This equation basically will eliminate the need to solve for the fluid pressure field (the first 3 equations) and for the temperature field (the 4th and 5th equations). Of course, both the use of the correlation and the assumptions made will make the results obtained by this method differ from reality. As I evolve in this work, I’ll be able to compare the results with experimental ones.

Continue to Part II >>