Author: Fernando Zigunov

I'm a brazilian Mechanical Engineer and PhD, currently a Postdoctoral Scholar at Florida State University, working in various projects in aerodynamics, aeroacoustics and flow control. In the past I've worked as a refrigeration systems designer and later as an R&D specialist at a refrigeration contracting company, researching for new products to push the industrial refrigeration market technologies forward. Currently I'm working with experimental aerodynamics at FSU, using state-of-the-art equipment and techniques to observe the most intricate flow physics!

Resonant Coupling Efficiency for Wireless Power Transfer

Posted on February 24, 2017October 22, 2020 by Fernando Zigunov in Hobby Electronics

So I’ve been working on a hobby project that involves transferring power to a circuit wirelessly. The distance is quite short (<10mm), as I just need to have one part of the device spinning and the other stationary.

I decided to stick with LC resonant coupling circuits as a solution to transfer this power, as I can admit for some power loss and having low audible noise is a requirement. But while I was studying this circuit, I wanted to know what would be the best conversion efficiency attainable, and also I was interested in understanding why in the first place a resonant circuit is needed (why not just use a regular circuit like a transformer?).

To answer these questions, I was thinking of building the circuit model in a simulator (as in the picture shown below), and simply playing with the parameters should give me a good grasp on what was happening.

While it did help me understand some concepts, it still was not very clear what was the effect of the coupling coefficient (k) on the power dissipated in the load and the overall power transfer efficiency of the system. So I decided to step up a little bit and use Matlab. The simple code I wrote is in this link and you can also play with the parameters for your circuit.

Basically, I used the control systems toolbox to write the transfer function of the circuit. The transfer function I’m interested on, eta, is the power dissipated in the load ZL divided by the power generated in the source Vin. The remaining power should be dissipated in the source itself, through its own impedance ZS.

The plot of frequency response for many coupling factors (as shown below) is rather interesting: It shows there is indeed a dependence of peak power transfer efficiency as the coupling factor increases. But it shows also that you don’t need a very large coupling factor to get interesting efficiencies (for k=0.4, eta peaks at 77%, for example). This means we don’t need astonishingly good coupling (as in a transformer, for example, where k>0.98) to transfer power. I know, this conclusion has already been made by smarter researchers than me long ago, but it’s nice to arrive at it by myself anyways.

Frequency domain response of the power transfer efficiency for ZL=20Ω, ZS=1Ω, C1=C2=100nF, L1=5.5µH, L2-6.0µH.

The following picture shows the same information, but complies the peak efficiencies against the coupling factor k. It shows that below a certain coupling factor (~0.4) the power transfer efficiency plummets. Unfortunately, this means that for large distances the resonant circuit is not the solution (as everyone already knows!), since the coupling factors are absurdly small for these distances.

Even so, for my simple peasant application of transferring power through a 10mm air gap that’s more than enough enlightenment in the topic.

eff2 — Peak efficiencies for the previous circuit parameters plotted versus k.

Just for the record, I measured the inductances of my homemade pair of coils (5.5µH and 6.0µH, as stated in the figure) and the coupling factor against distance (no capacitors, just the coils, using the open circuit voltage method, f=50 to 800 kHz) and the results are as shown below. According to the charts I generated, then, I should be fine, right? Well, it will depend whether I can make a MOSFET driver that will be able to attain a nice efficiency factor to multiply these maximum eta values I obtained!

k-vs-dist — Homemade coils and the coupling factor measurement results. The distance d is the gap between coils (I used a known thickness wood spacer)

A freezing tunnel (VRT) thermal model

Posted on February 9, 2017October 22, 2020 by Fernando Zigunov in Freezing Tunnel

As many refrigeration engineers around the world, I’ve always been curious to understand better how this machinery behaves transiently. For the ones who do not know the Variable Retention Tunnel, take a brief stop at Google before reading on, so you can appreciate better what I’m talking about!

In practice, we see that when the production rate increases, the heat load rises and the compressors might struggle to keep the evaporation temperature low. Therefore, the evaporator temperature rises, which makes the air temperature inside the freezing tunnel rise as well. The system is constantly looking for some equilibrium point, where the evaporator temperature fits for the current product mass flow rate.

But as each product cools down and eventually freezes, the heat transfer at its surface slows down because the surface temperature drops. This means some of the product inside the tunnel is not contributing to the heat load as much as the fresh, hot product that just entered the device. And all of this is happening simultaneously, which makes the appreciation and understanding of the process remarkably difficult.

So, recently while working on a project I had the idea: Why don’t I simulate the system so I could grasp better what is actually happening? I began, then, an endeavor towards this objective, and the concepts I could get in the path were quite enlightening.

1. The product comes in waves

Probably you are familiar with the product freezing curve shown below. It tells us the temperature evolution in time for some unfrozen product that is subject to a freezing process. The surface freezes first, then the low temperatures slowly (depending on the product’s thickness) reach the center of the product.

prodfreezing — Typical food product freezing curve

Another thing this curve also tells us, although indirectly, is the power the product releases onto the environment air. As $q'=h*A*(T_s-T_\infty)$ , and assuming the convection coefficient and the product surface area don’t change, we can assure that the power released $(q')$ by the product is proportional to the red curve in the figure.

The implications of this are quite interesting: As the curve shows, the temperature difference (and therefore, the power) drops by half as much in the first hour. So the product contributes the most to the heat load in the first hour, and then it “locks” the inner heat inside the product, making the heat transfer process much slower. As product is constantly coming in the tunnel during the loading hours, the effect looks like a wave of products, as shown in the video below:

The tunnel in this video has a 12h minimum retention time for any given product. This means during the working hours it is almost always full of “hot” product. When product stops coming, then it gradually stops “heating up” the tunnel and global heat load drops. But how does this translate in the total heat load transient profiles?

2. The daily oscillating heat load profile

Unfortunately, in practice we cannot measure the total heat load except with quite expensive equipment (and even so it would be quite inaccurate). But after simulating a whole tunnel, one can just integrate the data and get the total instantaneous heat load at each instant simulated.

Of course, the heat load profile will depend on the production schedule and on how long and when the compressors are turned off to save energy. In some countries, there is a strong variation on the energy fees depending on the hour of the day and on the season of the year. These variations are placed so the energetic system is not overloaded during the hours most people are home using electricity (normally early in the night), and companies are stimulated to turn off their machines when the fees are higher.

These factors make the heat load profile very rich and detailed. Below there is a picture of the same tunnel’s simulated profile during five days. The profile looks very repetitive, with the same peaks and troughs daily. During region A in the graph, the tunnel is empty and heat load builds up as product enters the tunnel. The heat load drops to zero in region E because the fans and the machine room are turned off between 18:00 and 21:00. Since the production is still ongoing during these times, the heat load peaks to B briefly during the machine room retake. This far exceeds the capacity of the compressors at the setpoint, and the machine room will struggle during this time. As product stops being fed into the system (at 0:00), the heat load stabilizes and, during the off-production hours (from 0:00 to 8:00, it drops until point D in the graph. This cycle runs around every day, and during the weekend break the power drops to zero (remaining only the conduction through insulation and the fan power). The weekly cycle begins again.

Simulated heat load profile of a typical VRT installation in Brazil.

3. The machine room struggles… But sometimes that’s okay

I thought that instead of assuming the machine room was “sized” to properly maintain the evaporator pressure at all times, I’d rather have it to “struggle” maintaining it at the setpoint during these high heat load periods of the day. “How would this affect the product?”, I asked myself. So let’s examine the chart below:

The tunnel starts the week, again, with virtually no heat load as the product’s temperature is already homogenized with the air temperature for staying so long during the weekend (region A). As product slowly loads the tunnel, the air temperatures (region B) rise but the machine room is still partially loaded and still keeping the evaporator temperature. In region C, the heat load of the product exceeded the machine room heat load, increasing the pressure of the evaporating fluid (and therefore, $T_{evap}$ ) even though the compressors are working at full load.

Simulated Evaporator and Air Temperatures for the same VRT installation.

In region D, the compressors are turned off for 3 hours, and the temperature readings are virtually meaningless. After turning on the compressors (region E), the huge heat load hits the machine room hard and makes the evaporator’s pressure rise all the way to a saturation temperature of -26ºC. Then as the heat load slowly settles down to a smaller value, the compressors begin to catch up.

This is a beautiful picture of the freezing tunnel and one that corresponds quite closely to what I saw in my experience. But if the compressors cannot keep the evaporator setpoint temperature, isn’t it a compressor sizing problem?

Well… It’s complicated. On one hand, the compressors normally are sized for a heat load that meets some average or “guessed” instantaneous criterion. Let’s say, for example, that one sizes the machine room capacity to meet the instantaneous heat load of an hour. The data used for this simple calculation is:

$\dot{m}=24000 kg/h$

$\Delta H=275 kJ/kg$

$ExtraPower (OtherDevices)=257 kW$

$Safety Factor=f_s=1.15$

The total power required for these compressors would, then, be:

$P=f_s * (\dot{m} * \Delta H + Extra Power) = 2395 kW$

Which corresponds to the magenta dotted line in the heat load chart. As the chart shows, the magenta line is below the peaks each and every day. But does this mean the compressors are undersized? This would mean they should actually be sized to a 9MW capacity, that occurs for a short half-an-hour. Is it worth it?

Well, it turns out, in practice, it’s not. The whole system has quite a lot of inertia. This tunnel stores 300 tons of product, which slows down the reaction time. Therefore, some transients can be accepted. But how far can we go?

Let’s say that another designer “wants” to size the same system for a daily average heat load. The system design heat load now will be reduced to $2395 kW*(21h/24h)=2095 kW$ . Let’s say we size our compressors for this heat load. This would mean we now are 300 kW (12.5%) shorter in compressor juice than we were before. So, obviously, we’re going to stay for longer in the “undersized” region, like the chart below displays. But what effect does this have in the product temperatures?

te2 — New evaporator temperature chart after reducing the compressor capacity by 12.5%

Well, let’s take a look and see what the product does. As different products stays different times in the tunnel (I’ll talk about this), then some products will exit at a lower temperature than others. This means there is some “statistical distribution” of the products, which is displayed in the next chart.

The chart actually shows the integral of statistical distribution, shows as a percentile. Let’s take the point (-30.2 ºC, 40%) in the upper chart as an example of how we read this chart. This point says that if I take a sample of the lowest 40% of the products outlet temperatures, the largest temperature would be -30.2ºC. In other words, 40% of the product is at -30.2ºC or lower.

The two charts above show that when we reduce the machine room compressor’s capacity, the average outlet temperature at the center of the products will be larger. But we need to compare it with a threshold, which is shown as a black dashed line in the chart (-18ºC). In the upper case (larger compressor), 99.6% of the product meets or exceeds the -18ºC threshold line. In the lower scenario (smaller compressor), a lower portion (97.1%) of the product meets the threshold line. To define whether these numbers are satisfactory or not is a decision the plant operator must make.

4. But I, the designer, shall dictate the capacity of the system!

The takeaway from the previous section is that the heat load of a given system (and this applies to other thermal systems as well) is not something we “decide”, but a product of the operational variables. Of course, we can control most of the operational variables, but we shouldn’t assume that’s always the case.

Another designer might say, for example, he “wants” to size the compressor to the weekly averaged heat load. He argues that the weekend time is idle, compressors are shut off unnecessarily during this time and he ought make the initial investment lower by putting more machines that will operate the whole week, therefore saving millions of dollars! Is that possible? Let’s see…

So the figure below shows how the temperatures behave for this compressor sizing. We can see that the oscillations are wilder, the peaks are roughly 5ºC above the previous simulations, and the machine room saturation temperature is constantly oscillating according to the tunnel heat load. But still at the weekend, the heat load drops and the compressors need to be shut off.

te3 — Temperature profile for a 1496 kW compressor (weekly sized).

So no improvement happened here, the compressors are still struggling during the week and being shut off during the weekend. Why? Well, that’s because we, as I previously mentioned, cannot dictate the heat load. The system does that. We just size for it.

So you’re probably curious to see how did the product temperatures perform under these “undersized” conditions. As shown in the figure below, we have now 75.4% of the product meeting the -18ºC center temperature criterion. It’s not remarkably bad, actually – but some clients will not be remarkably happy with some -10ºC measurements in their quality control spreadsheets. That’s the reason why we oversize!

temp2 — Effect of a 37.5% reduction of compressor capacity in the product temperature statistical distribution.

5. Conclusion

This is sort of a results article. I believe I’ll write some other articles on the mathematical devices I used to put the simulator together. I hope it was enlightening for you, the reader, as much as it was for me. If you think that kind of simulation would be useful for you, feel free to contact me.

Thanks for reading so far!

Mass file renaming and relinking in Inventor Part I: The results

Posted on November 16, 2016October 22, 2020 by Fernando Zigunov in Autodesk Inventor Tips

*Edit: New, translated, version posted on GitHub with install instructions. Check it out at https://github.com/3dfernando/Ferramentas-de-Montagem

1. Introduction

If you’re a mechanical designer for a living, chances are that you already had to copy an existing design. Even higher are the chances that you already had to rename all the files in a given assembly to match an ERP needs or a coding system needs.

Well, there are many ways of doing this, but I found that all the built-in ways in Autodesk Inventor have some sort of caveat:

Manual renaming in explorer: This is a pain in the ***. Apart from renaming the files one by one, you have to keep track of the old-new name correlations so when you open the assemblies in Inventor, you can redo the links manually.
Excel VBA to rename files: This is a little bit less cumbersome. You relieve the pain of renaming the files one by one, but the relinking still needs to be done manually in Inventor. At least you can print the correlation table from Excel.
Using the Design Assistant: This tool is nice, but there is a missing functionality that keeps it from being something usable for large assemblies: One cannot paste a table in the treeview. Imagine you as a designer already having the old-new naming correlations in an Excel table, but unable to paste them in Design Assistant. You still have to do the renaming manually. Sometimes it hurts productivity even more as the files are organized according to the tree, so finding them is another issue.
Using the CopyDesign add-in sample from Inventor SDK: This is an example piece of software that comes with the Inventor API SDK. It does what it’s supposed to do: Copies a design and renames the files. The problem is, you need to have the new names beforehand. And new files (which are probably going to happen in a copied design) will have to be named correctly during design, which is not always possible under an ERP system.

2. My solution

So, now that I showed the available options and how they did not meet my needs, let me state what I want to do and what I came up with using the Inventor API.

I want to make a file list using Excel VBA (which will list the file names for the current assembly) and, using some made-up rule, assign new names to the existing files. Here I think Excel will always be way more flexible than any piece of software I would ever write, so I decided to leave the rule creation to the freedom of the user. A table like the one shown below will be the output from Excel as a CSV file and read by this program:

existing-new — Input CSV file format for the renaming system

I open the main assembly file and import the correlation table on my API-brew software. It figures out what files need to be renamed if they match the name in the “from” column of the correlation table. It finds out whether there are any equally-named *.idw files in the path. It renames all the files and then redoes the links in all levels of the assembly. Finally, it saves the assembly.

That’s the procedure in a nutshell. The heavy work is done by the software. So how does it look in Inventor? Follow the step-by-step screenshots shown below.

2.1. Main assembly loaded in my API interface

***Edit: I changed the button name to “Batch Renaming“.

Screen of the “Batch Renaming” add-in showing the list of components referenced by the assembly

2.2. Loading the CSV file

CSV file being selected in the file browser

2.3. Correspondences loaded in the file being displayed on-screen

In this case, not all files were selected to be renamed as some parts belong to the library.

Files to be renamed are shown in the table in bold

2.4. Renamed files in Explorer

fifthscreen — Files after the renaming in the explorer view…

sixthscreen — … and in the Inventor treeview!

3. Conclusion

So the use of this tool saved us quite a lot of time on the renaming of files. For larger assemblies, of >700 parts, a reliable manual renaming took up to 3-4 hours and could be reduced to something like 10 minutes. If you think this tool might be useful for you too, please let me know in the comments!

I published it on GitHub in this address: https://github.com/3dfernando/Ferramentas-de-Montagem. Please check it out and let me know if it was useful for your job!

If you guys understand Portuguese, you might find this video I made in my language for my comrades useful. I’ll soon produce an English version with the updated code so more audience can understand it. **EDIT: there it is =)

**Edit 2: Watch the installation process if you’re having trouble:

The freezing heat exchanger: Part III – Some insight

Posted on October 24, 2016November 18, 2016 by Fernando Zigunov in Heat Exchanger Design

So in Part I and Part II we spoke about how to solve numerically the freezing heat exchanger. After some tinkering with the results I got, I realized something very simple: The maximum heat transfer rate in a freezing heat exchanger is not dependent on:

The ice thickness;
The non-freezing fluid heat transfer coefficient;
The material/thickness of the tube wall.

Which is quite unintuitive! Basically, the maximum heat transfer rate of a freezing HX happens when it is completely frozen. And then that happens, we know what’s the temperature of the liquid boundary: $T_{wall}=T_{freeze}$ (Of course, assuming no supercooling).

This allows us to analyze a freezing heat exchanger based on the maximum heat transfer rate possible, in a way similar to the ε-NTU method. Let’s call it the ice-efficiency method, and name an ice-efficiency factor $\varepsilon_{ice}$ , defined as:

$\varepsilon_{ice}=\frac{q_{HX}}{q_{fully-frozen}}$

Where $q_{HX}$ is the effective heat transfer rate of the heat exchanger, while $q_{fully-frozen}$ is the heat transfer rate of the same heat exchanger, should it be fully frozen. Of course, $q_{HX}<q_{fully-frozen}$ .

Actually, a similar approach has already been suggested as a “dimensionless heat transfer rate variable” $q^*$ by Zerkle (1964) in his PhD thesis, so it’s not something really that innovative. The difference here is that I’m using the maximum heat exchangeable by the HX the size it is, instead of the maximum heat exchangeable by an “infinitely long” HX.

$q_{fully-frozen}$ can be calculated using the ε-NTU method, assuming $T_{wall}=T_{freeze}$ . First, let’s define $q_{max}=q_{HX}/\varepsilon$ :

$q_{max}=C_{min}(T_{h,i}-T_{freeze})$

And the number of transfer units:

$NTU=\frac{h_i A}{C_{min}}$

An important thing to notice in the NTU equation is the fact that U, the global heat exchanger coefficient, has been replaced by $h_i$ , the internal convection coefficient. This is because given a wall temperature ( $T_{freeze}$ ), the heat transfer rate is now independent on the other thermal resistances. What this means is that, once fully frozen (which is the best case), the maximum heat transfer possible is only dependent on the liquid convection intensity.

Two takeaways come from this interpretation: The first is that (mild) solidification is desirable in a freezing heat exchanger. The second is that this sets up the operational limit. Of course, different from the ε-NTU method, the way to find $\varepsilon_{ice}$ is not so clear. Maybe we can come up with a correlation for that?

The freezing heat exchanger: Part II – Solution

Posted on October 24, 2016November 18, 2016 by Fernando Zigunov in Heat Exchanger Design

1. Introduction

So in Part I (link here) we devised a method to possibly solve the freezing heat exchanger problem. Here we’re going to look into the numerical method itself and what are the results we can get with this tool.

2. Solving it numerically

We began by stating the fundamental equations of this 2D flow, which need to be integrated in three dimensions (two space dimensions and one time dimension). By eliminating the need to solve for the flow field, the integrator will only need to solve for one space and one time dimension, greatly simplifying both the integration method needed and the number of computations to be made.

The solution method I devised is explained graphically in the diagram below:

solution-diagram — Integrator diagram to solve for the freezing HX

This is the basic skeleton of the integrator. We can add a lot of stuff there to adapt the solution to our needs. For example, in my case I added the dynamic behavior of the shell side, so I could implement a controller for a motorized flow valve. This allowed me to understand the control behavior and how the HX responded to different conditions.

I implemented this method in VBA for MS Excel, and although I know this isn’t the best tool for this job it’s the one we can afford here. It doesn’t actually matter which programming language you use, as the calculations aren’t that intensive.

3. The results

Below I’m showing some sample results of what we can get with this simple methodology. The graphs are shown for a heat exchanger comprised of stainless steel piping, $d_e=38.1mm$ and $d_i=35.1mm$ , $L=5.0m$ and 16 tubes in total. It’s arranged in 4 circuits of 4 passes each, and the shell-side fluid is ammonia evaporating at -5ºC. The tube side fluid is water, which is pumped at a rate of 70 m³/h, entering at +2ºC.

From the graphs below, we can see that the ice builds up to a maximum thickness of 0.66mm at the outlet, which stabilizes after ~200 s. This means that in this condition, the heat exchanger will not clog with ice.

icethickness — Graphs generated by solving the equations presented

So what happens if we reduce the water flow rate to half as much? The results for this condition are presented below. Look at how the ice layer builds up mich thicker now, going towards a thickness of 2.15mm and growing. Also, note in the instantaneous power graph (first graph) how the shell side power (power delivered to the fluid in the shell) is now much different to the tube side power (power delivered to the fluid in the tube). This is exactly due to the ice buildup, which requires heat to be released.

Even in this condition, it is not very likely that the HX will shut off with an ice plug. We can see that the process conditions must be very harsh in order for such a failure to occur.

So let’s see what happens if, due to a hypothetic process heat load drop, the temperature of the water drops to 1.0ºC instead. The simulation time now was increased to 1000s (~17min) to show the long-term effects. The ice now builds up to 4.5mm thick in the outlet. The entire heat exchanger is now frozen to some thickness. Even so, it seems to be in a stable condition, and if the pump can still move those 35m³/h through the HX, it’ll probably not clog up.

Well, we cannot assure anything about the integrity of the pipe, as the stresses in the pipe wall increase as ice builds up. For this, we need to perform a different calculation.

As a last graph I wanted to show, take a look at what happens if we simulate this HX, with 70m³/h of water coming in at +2ºC, if we begin with a constant ice thickness of 5mm:

The ice thickness, as one might have expected, decrease to the operational condition. This means that we can simulate both HX freezing and thawing with this method, possibly opening up a way to simulate the behavior of ice banks.

I hope these plots displayed the capabilities of the method and inspired you to make your own version. If you have any doubts, feel free to send me an e-mail or a comment!

Continue to Part III >>

The freezing heat exchanger: Part I – Introduction

Posted on October 21, 2016November 18, 2016 by Fernando Zigunov in Heat Exchanger Design

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 >>