Department of
Physics and Astronomy Uppsala University Box 516
SE-75120 Uppsala Sweden
Department of Physics and Astronomy Uppsala University
FREIA Report 2020/01 July 9, 2020
Benchmarking a Cryogenic Code
for the FREIA Helium Liquefier
Elias Waagaard
Supervisor: Volker Ziemann
Subject reader: Roger Ruber
Bachelor Thesis, 15 credits
Uppsala University
Abstract
The thermodynamics inside the helium liquifier in the FREIA laboratory still contains many unknowns. The purpose of this project is to develop a theoretical model and im-plement it in MATLAB, with the help of the CoolProp library. This theoretical model of the FREIA liquefaction cycle aims at finding the unknown parameters not specified in the manual of the manufacturer, starting from the principle of enthalpy conservation. Inspiration was taken from the classical liquefaction cycles of Linde-Hampson, Claude and Collins. We developed a linear mathematical model for cycle components such as turboexpanders and heat exchangers, and a non-linear model for the liquefaction in the phase separator. Liquefaction yields of 10% and 6% were obtained in our model simula-tions, with and without liquid nitrogen pre-cooling respectively - similar to those in the FREIA liquefier within one percentage point. The sensors placed in FREIA showed simi-lar pressure and temperature values, even though not every point could be verified due to the lack of sensors. We observed an increase of more than 50% in yield after adjustments of the heat exchanger design in the model, especially the first one. This constitutes a guideline for possible future improvements of the liquefier.
Sammanfattning
Contents
1 Introduction 1 2 Fundamental Thermodynamics 2 3 Cycle Components 6 3.1 Compressor . . . 6 3.2 Heat Exchanger . . . 63.2.1 Counterflow Two-fluid Heat Exchanger . . . 7
3.2.2 Counterflow Three-fluid Heat Exchanger . . . 11
3.3 Turboexpander . . . 15
3.4 Joule-Thomson Valve . . . 17
4 Classical Liquefaction Cycles 18 4.1 Linde-Hampson Cycle . . . 18
4.2 Claude Cycle . . . 20
4.3 Collins Cycle . . . 23
5 The FREIA Helium Liquefier 25 5.1 Simulations and Results of the FREIA Model . . . 28
6 Discussion 33 7 Future Prospects 35 8 Conclusion 35 9 Popular Science Summary 36 10 Acknowledgements 37 References 38 11 Appendices 39 11.1 FREIA LN2 Enthalpy Model Function . . . 39
1
Introduction
Cryogenics, the study of material behaviour at low temperatures, has a long tradition in Upp-sala. An important step to better understand the concept of temperature was taken there in 1742, when Anders Celsius launched his centigrade thermometer scale. This scale was based on phase transitions: the freezing and boiling point of water at a particular atmospheric pres-sure. This standardized temperature measurement around the world until this day, although he initially proposed that water should boil at 0◦C and freeze at 100◦C. This scale was reversed after his death and then accepted worldwide [1]. Today, the FREIA Laboratory at the Uppsala University carries out experiments for accelerator physics and instrumentation at low temper-atures with liquid nitrogen and liquid helium as coolants. However, lighter gases such as these are not trivial to liquefy. Carl von Linde and William Hampson independently succeeded with the liquefaction of air in 1895 using the cycle that was later named after them [2]. It would take until 1908 for helium to be liquefied for the first time, achieved by Heike Kamerlingh Onnes with a pre-cooled Linde-Hampson cycle [3].
Nowadays, cryogenics play a fundamental role to supply coolants to a range of scientific experi-ments, medical applications and industrial machinery. For instance, liquid helium is often used for the cooling of superconducting magnets in Functional Magnetic Resonance Imaging (fMRI) machines [4]. The need for low-temperature equipment in industry is also growing rapidly. The industrial gas market size is estimated to grow from USD 80 Billion in 2015 to reach more than USD 100 Billion by 2020, dominated by large international corporations such as the Linde Group and Air Liquide [5]. Only the system of superconducting magnets at the Large Hadron Collider at CERN use a total of 120 metric tonnes of helium to be cooled down to the working temperature of 1.8 K [6].
Due to the abundant applications and few deposits on Earth, the supply of liquid helium is limited and causes very fluctuating market prices. As pointed out by Sophia Hayes in Physics Today, the market price of liquid helium was USD 5 per litre in 2010, but this had almost quadrupled by 2016 [7]. This volatility of the global helium supply stresses even further the need for robustly designed and leakproof helium liquefaction systems. Moreover, simulations help to optimize and guide towards improving the efficiency where the information provided by the manufacturer is not sufficient.
In the FREIA Laboratory, the helium liquefier Linde L140 built by the Linde Group delivers over 140 l/h of liquid helium into a 2000 l storage dewar, and includes liquid nitrogen pre-cooling and a helium recovery system [8]. However, the process of helium liquefaction still contains many unknown parameters at several points in the cycle, where sensors are missing. These unknowns include gas flows, temperatures and pressures. The purpose of this project is to construct a theoretical model of the liquefier, starting from drawings of the vendor. This theoretical model is based on balancing the flow of enthalpies through the system at thermodynamic equilibrium, and is implemented in a model simulation in MATLAB. The system parameters from the FREIA control system, such as temperatures and pressures, are compared with output parameter values of the simulation.
In this report, we first describe some essential terminology in thermodynamics in Section 2, followed by mathematical models for the individual cycle components in Section 3. These com-ponents will then be put together in liquefaction cycle models that are simulated in MATLAB. We will gradually increase the complexity of the models, starting from the historically impor-tant liquefaction cycles of Linde-Hampson, Claude and Collins in Section 4, to finally simulate the FREIA liquefaction cycle in Section 5. The MATLAB code to simulate the FREIA lique-faction cycle can be found in Appendices 11.1 and 11.2. A popular science summary can be found in Section 9.
2
Fundamental Thermodynamics
To start with, several recurrent thermodynamic concepts need to be defined. We will often come back to the first law of thermodynamics, which states the conservation of the internal energy U , or the energy associated to the random motion of molecules in a system. According to the first law, an infinitesimal change dU in the internal energy is
dU = δQ + dW, (1)
where δQ is the infinitesimal heat of the process and dW is the infinitesimal work. In other words, the sum of all incoming energies in a system equals the outgoing energies or what is accumulated [2].
To describe a system, we use principal properties such as pressure P , volume V , molar quantity n, and temperature T . P and T do not depend directly on the size or the amount of material in the system, and are called intensive quantities. P is the force exerted per unit area of the system, and T is a quantitative expression of the heat in the system. For such a system, we express the infinitesimal work as dW = −P dV . In this report, we shall always give P in terms of absolute pressure, in unit [bar(a)]. Conversely, V and n are extensive quantities and depend on the size of the system, or how much is contained in it. The thermodynamic relation between these quantities is called an equation of state. In classical thermodynamics, the equation of state for an ideal gas is
P V = nRT, (2)
where n is the amount of substance of the gas (in moles) and R = 8.314 [J/(mol K)] is the ideal gas constant. An ideal gas is a theoretical gas for which Equation (2) holds for all pressures P and temperatures T [10]. Gases such as helium, nitrogen, oxygen but also some heavier gases such as CO2 behave like ideal gases to the first order, and approach a perfect ideal gas at higher
EK = 1 2mv 2 rms = 3 2kBT, (3)
where vrms =phv2i is the root mean square velocity of the gas particles and m its individual
mass [10]. Consequently, EK,T ot ≡ U = 32N kBT if summing up the energies of all the N
particles. For ideal gases, we can also easily relate heat and temperature. We define the heat capacity CV of a gas at constant volume (isochoric) as the heat differential dQ needed to
increase the system temperature with dT
CV = ∂Q ∂T V , (4)
where the unit is [J/K]. From Equation (1), (3) and the fact that W = −P dV = 0 for isochoric processes, we have in particular for ideal monatomic gases
CV = ∂ ∂T(U − W )V = ∂ ∂T 3 2N kBT = 3 2N kB = 3 2nR. (5)
The specific heat capacity cV [J/(mol K)] is then
cV ≡
CV
n =
3
2R (6)
for an ideal monatomic gas. If we instead are interested in the specific heat capacity at constant pressure cP, we resort to Mayer’s relation [10]
cP − cV = R, (7)
which results in cP = 52R [J/(mol K)]. These values of the specific heat capacities shall be
frequently used when we have gases that behave like ideal gases, such as helium at higher T . In the thermodynamic cycles that are to come, cP is more frequently used as gas may change
volume when liquefied, but the pressure stays the same. Specific heat capacities cP of helium
at various pressures, including those that will later be observed in the FREIA liquefier, are presented in Figure (1).
Let us elaborate further on the concept of heat. If we have a totally isolated system, we can express the internal energy and the pressure-volume P V term in terms of the so-called enthalpy H, defined as
H = U + P V. (8)
Figure 1: Specific heat capacities at constant pressure cP of helium, generated with CoolProp.
cP = 52R is also marked out as a black horizontal line.
L = Hvap− Hliq, (9)
where Hvap is the enthalpy of the vapour and Hliq is the enthalpy of the liquid. We obtain the
specific latent heat l by simply dividing both sides in Equation (9) by the concerned mass m. Another essential concept is the entropy S, also an extensive quantity. The formal definition of an infinitesimal entropy increment dS of a process on a closed system (i.e. where heat but not matter can enter) is
dS ≥ δQ
T , (10)
For this study, we are particularly interested in the thermodynamics regarding liquefaction. We need to consider that there are situations in which vapour and liquid phases can coexist. The critical point with temperature TC and pressure PC is the point in a T S diagram from
which coexistence of phases is possible. Below this point, the so-called saturation dome starts, a bell-shaped region in the T S diagram where vapour and liquid phases co-exist. The percentage of vapour and liquid respectively is determined by the vapour quality, where a fully saturated vapour has a vapour quality of 100% and fully saturated liquid has a vapour quality of 0%. As we shall see later in this study, these phenomena are crucial for liquefaction. An illustration of the saturation dome can be seen in Figure (2) to the right. The dashed red line to the right of the critical point is called the saturated vapour line (100% vapour quality), and the red dashed line to the left of the critical point is the saturated liquid line (0% vapour quality). The zone where phase coexistence occurs is underneath the critical point inside the bell-shaped curve, between these saturation lines.
Step 3 describes the isenthalpic process that leads to liquefaction, if temperatures are low enough. What is liquefied of the gas then goes towards the left of the saturation dome (as saturated liquid), and what remains gaseous progresses to the right as saturated vapour. This gas is then reheated along step 4 isobarically, before it reaches the initial point of the cycle. If the temperature is too high at the start of step 3, one can include isentropic expansion as a part of a multi-stage cooling process of the gas to achieve liquefaction. This process is illustrated in step 5 in Figure (2), and will be further explained along with the Claude cycle in Section 4.2. One of the objectives of the project is to extract such a T S diagram for the FREIA cycle.
Figure 2: Sketch of Carnot cycle displayed in a T S diagram to the left, and a typical Claude liquefaction cycle displayed to the right. The liquefaction cycle is represented in the blue solid line, and the saturation dome in dashed red.
3
Cycle Components
The various thermodynamic components that build up a cycle are described in the following section, including code snippets in MATLAB with examples of how they practically can be implemented. Some involve more theoretically complex physical phenomena, whereas others are simpler.
3.1
Compressor
A compressor (CMP) is a device that increases the pressure P of a gas by reducing its volume V , by compression [2]. In the thermodynamic cycles that will be presented, the compressor is often the first step. Ideally, this process is isothermal, such that the gas keeps its original temperature while P is increased. This is achieved by connecting the compressor system to a heat bath with large heat capacity. In this study, we often use the high-pressure outlet of the compressor as the start of the cycle.
3.2
Heat Exchanger
Heat exchangers are devices that transfer enthalpy between two or more flows of fluids or gases, without transferring the gas or fluid itself. These hot and cold flows can be parallel, anti-parallel, or arranged in a more intricate way. In cryogenics, counterflow heat exchangers are often preferred to parallel flows as larger temperature drops can be achieved. This is due to the fact that the outlet temperature T4 of the cold fluid can be much higher than the outlet
temperature T2 of the hot fluid in a counterflow system, which is not possible in a parallel flow
[2]. In Figure (3), a simple one-dimensional counterflow heat exchanger is presented.
Figure 3: Illustration of a simple counterflow heat exchanger. The hot flow from 1 to 2 is marked in red, and the cold flow from 3 to 4 is marked in blue.
The principle of enthalpy transfer is always the same regardless of the type. If enthalpies at the inlet points H1, H3 and the enthalpy transfer differential dH are known, we simply add to
the cold flow and subtract from the hot flow to solve for H2 and H4, as presented in the code
snippet below. However, we are often interested in the total enthalpy transferred per second across the whole heat exchanger, which we denote ∆ ˙H.
3.2.1 Counterflow Two-fluid Heat Exchanger
If only the inlet temperatures T1and T3 for a heat exchanger are known, how do we calculate the
final temperatures T2 and T4, and the total enthalpy ∆ ˙H transferred? We start with a simple
one-dimensional counterflow heat exchanger (CX) with two flows, assuming an incoming hot flow with total heat capacity per second Ca= ˙maca and a return flow with total heat capacity
per second Cb = m˙bcb, where ci are the specific heat capacities at constant pressure of the
fluid and ˙mi are the mass flows (i = a, b). We have temperatures T1, T2, T3 and T4 as seen
in Figure (4). Let us also assume that the length of the heat exchanger is L, and that C0 is the heat exchanger design parameter. C0 constitutes the heat transfer coefficient times the contact area, often also denoted U A, indicating how easily heat is transferred between the hot and the respective cold streams. In addition, we do not know the exact nature of the flow inside the heat exchanger. The so-called Reynolds number is a parameter that characterizes the flow patterns in fluids. This behaviour affects how much the fluid comes in contact with the inner walls, and consequently how much enthalpy the fluid absorbs or emits. At low Reynolds numbers, the flow is laminar (sheet-like), whereas high Reynolds numbers indicate a turbulent flow [2]. Turbulent flows increase chances that all the fluid will come in contact with the inner walls of the heat exchanger, meaning a higher transfer of enthalpy. In our case however, this effect is included in the parameter C0. We once again look into an idealized one-dimensional case, but it can easily be applied to three dimensions if U A is used instead of C0.
Figure 4: Illustration of simple one-dimensional counterflow heat exchanger with incoming flow a and return flow b, with inlet temperatures T1 and T3, outlet temperatures T2 and T4, and
heat capacities per second Ca and Cb.
For the flow in a short segment of length dx, we set up the system of equations for the differential heat transfer rate
d ˙Ha= −CadTa (from hot fluid) (11)
d ˙Hb = −CbdTb (to cold fluid) (12)
d ˙H = C0(Ta− Tb)dx (between fluids), (13)
where d ˙Ha, d ˙Hb and d ˙H are the instantaneous enthalpy transfers from the incoming flow, the
is absorbed in the cold flow, which gives a plus sign. However, a minus sign arises as we move in the direction of −dx, as can be seen in Figure (4). The heat emitted from one side corresponds to the absorbed heat on the other side as energy is conserved, so we have d ˙Ha = d ˙H = d ˙Hb.
Hence, we can manipulate Equation (11), (12) and (13) and obtain dTa dx = − C0 Ca (Ta− Tb) and dTb dx = − C0 Cb (Ta− Tb), (14)
or expressed in matrix form
d dx " Ta Tb # ="− C0 Ca C0 Ca −C0 Cb C0 Cb # " Ta Tb # . (15)
The eigenvalues α to the matrix in Equation (15) are obtained from solving the characteristic equation α + C 0 Ca α − C 0 Cb + C 02 CaCb = 0. (16)
The solutions to Equation (16) are α1 = 0 and α2 = −C0
1 Ca− 1 Cb
, with respective eigenvectors ~ v1 = 11 and ~v2 = C0 Ca C0 Cb
. As Equation (15) is a linear system of first-order homogeneous differential equations, we know that its solution can be expressed as
" Ta Tb # = A1~v1eα1x+ A2~v2eα2x = "−C0 Ca C0 Cae −αx −C0 Cb C0 Cbe −αx # " A1 A2 # , (17)
where A1 and A2 are integration constants and α = C0 C1a −C1b. These are fixed by imposing
the boundary conditions Ta(x = 0) = T1 and Tb(x = L) = T3, which means
" T1 T3 # ="− C0 Ca C0 Ca −C0 Cb C0 Cbe −αL # " A1 A2 # ⇒ " A1 A2 # ="− C0 Ca C0 Ca −C0 Cb C0 Cbe −αL #−1" T1 T3 # . (18)
Solving for A1 and A2, we find the integration constants to be
A1 = C0 Cbe −αLT 1− C 0 CaT3 C0 Cbe −αL− C0 Ca and A2 = T3 − T1 C0 Cbe −αL− C0 Ca . (19)
T2 = T1− 1 − e−αL 1 −Ca Cbe −αL(T1− T3) ≡ T1− η(T1− T3) T4 = T3+ Ca Cb 1 − e−αL 1 − Ca Cbe −αL(T1− T3) ≡ T3+ Ca Cb η(T1− T3) (20) where η = 1−e−αL 1−Ca Cbe
−αL is defined as the efficiency of a counterflow heat exchanger. Therefore,
we see that the heat exchanger temperatures constitute a system of linear equations that can easily be solved.
If, on the other hand, we are interested in the total enthalpy transfer across a counterflow heat exchanger as a function of the inlet temperatures T1 and T3 as in Figure (4), we resort once
again to Equation (13) and our solutions from Equation (17)
d ˙H = C0(Ta(x) − Tb(x))dx = C0αA2e−αxdx. (21)
Integrating over the whole heat exchanger gives
∆ ˙H = C0A2 Z L 0 αe−αxdx = CaCb(1 − e −αL ) Cae−αL− Cb (T3− T1) ≡ CH(T3− T1) (22)
if Ca 6= Cb. This constant CH = CaCb(1−e
−α)
Cae−α−Cb is called the heat conduction constant of a heat
exchanger.
However, if Ca = Cb, the reasoning thus far breaks down as the matrix in Equation (15)
is singular with degenerate eigenvalues, so it is not invertible. For this case, we manipulate Equation (11), (12) and (13) yet again but in a slightly different order
dTa− dTb = d(Ta− Tb) = −d ˙H 1 Ca − 1 Cb = −C0(Ta− Tb) 1 Ca − 1 Cb dx = 0 (23)
which implies that
d(Ta− Tb)
dx = 0, ⇒ Ta− Tb = D (24)
where D is an integration constant. Now we once again employ Equation (11) and (12), meaning
− dTaCa= C0(Ta− Tb)dx dTa dx = C0 Ca (Ta− Tb) = − C0 Ca D ⇒Ta= E − C0 Ca Dx, (25)
Tb = Ta− D ⇒ Tb = E − D − C0 Ca Dx = E − D 1 + C 0 Ca x . (26)
With the boundary conditions Ta(x = 0) = T1 and Tb(x = L) = T3, we immediately observe
that E = T1 and consequently deduce that
T3 = T1− D 1 + C 0 Ca L so D = T1− T3 1 + CC0 aL . (27)
This finally leads to
Ta(x) = T1+ C0 Ca(T1− T3) 1 + CC0 aL x, Tb(x) = T1 + 1 + C 0 Ca (T1− T3) 1 + CC0 aL x. (28)
If we now solve for the total enthalpy flow
d ˙H = C0(Ta(x) − Tb(x))dx = C0Ddx ∆ ˙H = Z L 0 C0Ddx = C0LT1− T3 1 + CC0 aL ≡ CH(T1− T3) ⇒ 1 CH = 1 Ca + 1 C0L (29)
if Ca= Cb. No matter if the heat capacities are equal or not, the total transferred enthalpy per
second ∆ ˙H over the heat exchanger can be found.
The code snippet below shows how CH can be calculated in the various cases, to finally find
∆ ˙H if Ca, Cb, C0 and L are known. The if-statement in line 6 makes sure to avoid numerical
singularities if the exponent is too small.
3.2.2 Counterflow Three-fluid Heat Exchanger
When three flows of fluid or gas are involved, the situation is more complex than for only two flows. In a general three-fluid counterflow heat exchanger, we assume a hot incoming flow with heat capacity per second Ch between two parallel cold return flows with respective heat
capacities per second C1 and C2, as can be seen in Figure (5). C10 and C 0
2 represent the heat
exchanger design parameters on the respective sides, or the respective heat transfer coefficient times the contact area. Also the input temperatures Th,i, Tc1,i and Tc2,i are known, and will
constitute our boundary conditions.
Figure 5: Three-fluid heat exchanger with one hot flow with temperature Th in the centre,
surrounded on each side by two parallel cold flows with temperatures Tc1and Tc2, with direction
opposite to the hot middle flow.
Just like in the case with a two-fluid counterflow heat exchanger, we set up the equations for the heat balance for a small segment dx
C1 dTc1 dx = −C 0 1(Th− Tc1) (30) Ch dTh dx = −C 0 1(Th− Tc1) − C20(Th− Tc2) (31) C2 dTc2 dx = −C 0 2(Th− Tc2). (32)
d dx Tc1 Th Tc2 = C10 C1 − C10 C1 0 C10 Ch − (C10+C02) Ch C20 Ch 0 −C20 C2 C20 C2 Tc1 Th Tc2 , (33)
which, just like in the two-fluid case, leads us to deduce the ansatz Tc1 Th Tc2 = X i Ai~vieαix, (34)
where αi represent the eigenvalues of the matrix M in Equation (33), ~vi the corresponding
eigenvectors and Aiintegration constants (i = 1, 2, 3). The characteristic equation |M −αI| = 0
gives α −α2− C 0 1 C1 +C 0 2 C2 +(C 0 1+ C 0 2) Ch α−C 0 1C 0 2 C1C2 +(C 0 1+ C 0 2) Ch C0 2 C2 +C 0 1 C1 − C202 C2Ch + C 02 1 C1Ch = 0, (35) where I is the 3×3 identity matrix. We immediately see that α1 = 0 forms one solution. α2
and α3 are found by solving the second-order polynomial equation
− α2− C 0 1 C1 +C 0 2 C2 +(C 0 1+ C 0 2) Ch α − C 0 1C 0 2 C1C2 +(C 0 1+ C 0 2) Ch C0 2 C2 + C 0 1 C1 − C202 C2Ch + C 02 1 C1Ch = 0, (36) that is a quadratic equation and can be solved easily, though the expressions become rather long. To each of these eigenvalues α1, α2 and α3, there is a correlated eigenvector ~v1, ~v2 and ~v3.
As α1 = 0, we conclude that ~v1 = v11 v21 v31 =11 1 in order to fulfil (M − α1) ~v1 = ~0. ~v2 = v12 v22 v32 and ~v3 = v13 v23 v33
can easily be found with numerical methods, such as the eigenvalue and eigenvector finder function in MATLAB. Consequently, the general solution to the temperature profile described by Equation (33) is
Tc1 Th Tc2 = A1v~1eα1x+ A2v~2eα2x+ A3v~3eα3x = v11eα1x v12eα2x v13eα3x v21eα1x v22eα2x v23eα3x v31eα1x v32eα2x v33eα3x A1 A2 A3 (37)
where A1, A2 and A3 are integration constants determined the boundary conditions Tc1(x =
L) = Tc1,i, Th(x = 0) = Th,i and Tc2(x = L) = Tc2,i. Imposing these boundary conditions we
Tc1,i Th,i Tc2,i = v11eα1L v12eα2L v13eα3L v21 v22 v23 v31eα1L v32eα2L v33eα3L A1 A2 A3 ⇒ A1 A2 A3 = v11eα1L v12eα2L v13eα3L v21 v22 v23 v31eα1L v32eα2L v33eα3L −1 Tc1,i Th,i Tc2,i , (38)
which lets us determine the integration constants. At the end of the day, just like in the case of a two-fluid counterflow heat exchanger, we find that the relation between the temperatures is described by a system of linear equations that can easily be solved, at least numerically. If we are interested in the actual enthalpy transfer between the flows, we recall Equation (30) and apply it on the hot stream and cold return stream 1
d ˙H1 = C10(Th(x) − Tc1(x))dx. (39)
If we substitute the temperatures from Equation (37), we find that
d ˙H1 = C10
A1(v21− v11)eα1x+ A2(v22− v12)eα2x+ A3(v23− v13)eα3x
dx. (40)
Integrating over the whole heat exchanger gives
∆ ˙H1 = C10 Z L 0 A1(v21− v11)eα1x+ A2(v22− v12)eα2x+ A3(v23− v13)eα3x dx = C10 A1(v21− v11)(e α1L− 1) α1 + A2(v22− v12)(e α2L− 1) α2 +A3(v23− v13)(e α3L− 1) α3 , (41) which is the total enthalpy transfer between the hot stream and cold stream 1. Analogously we know that d ˙H2 = C20(Th(x) − Tc2(x))dx, so one obtains a similar expression for the total
transfer between the hot stream and cold stream 2
∆ ˙H2 = C20 A1(v21− v31)(eα1L− 1) α1 + A2(v22− v32)(e α2L− 1) α2 +A3(v23− v33)(e α3L− 1) α3 . (42) For this type of system, one of the eigenvalues αi is zero-valued, with corresponding eigenvector
~ vi = 1 1 1
. This means that any vector element difference vi,j − vi,k = 0 before any integration
takes place. Hence, the term relating to the zero-valued eigenvalue will vanish from Equation (41) and (42). In contrast to the case of a two-fluid counterflow heat exchanger, the matrix in Equation (33) is not singular even in the case when Ch = C1 = C2. In the latter case of equal
heat capacities, it can be found from Equation (35) that α1 = 0, α2 = −α3, so the system is
Equation (41) and (42) constitutes the general solution to the enthalpy transfer in a three-way heat exchanger. Any difference in order or direction of the cold streams and hot stream will simply add or remove minus signs in the heat balance in Equation (33), whether we go along or against the dx element. One of the three-fluid heat exchangers in the FREIA liquefier cycle resembles to Figure (5), but there is also one three-fluid heat exchanger where the cold flow 2 is parallel with the central hot flow, but is still anti-parallel to the cold flow 1. Consequently, due to the fact that we move along dx, the minus sign will vanish in Equation (32) and the CX temperatures are instead described by
d dx Tc1 Th Tc2 = C0 1 C1 − C0 1 C1 0 C0 1 Ch − (C0 1+C20) Ch C0 2 Ch 0 C02 C2 − C0 2 C2 Tc1 Th Tc2 . (43)
We must not forget to adjust the boundary condition such that Tc2(x = 0) = Tc2,i, and the
conditions for Tc1 and Th remain unchanged. However, the final solution will still be of the
form expressed in (37). The two-fluid and three-fluid cases are compared in Figure (6) with numerical test values and length L = 1 m. The left-most and the central plot represent the systems as illustrated in Figure (4) and (5) respectively. The right-most plot corresponds to Figure (5), with the difference that the central hot flow is only anti-parallel to cold flow 1 and parallel to cold flow 2. All of these occur in the FREIA liquefaction cycle.
Figure 6: Example plot of temperatures from two-fluid and three-fluid counterflow heat ex-changers of length L = 1 m, with corresponding layout above.
For the three-fluid heat exchangers, we considered C10 = 300, C20 = 300, C1 = 200, C2 = 500 and
Ch = 500 as test values (in [J/K]). The incoming cold temperatures are 200 K and the incoming
are kept but the second cold flow is simply removed, so C0 = C10, Ca = Ch and Cb = C1. The
hot and the cold input temperatures are also kept. A much larger temperature decrease for the hot flow is achieved with two cold flows with respect to a single cold flow. Maximum cooling occurs in the central plot, when the hot flow is opposing both cold flows.
3.3
Turboexpander
Turboexpanders (TXP), or expansion engines, are devices allowing the gas to perform work on a spinning rotor and expand. This expansion results in large temperature reductions, so turboexpanders act as very efficient means of refrigeration. The process occurs ideally isentrop-ically, so ∆S = 0. In real processes, there will however be a small change in entropy. For the mathematical background on turboexpanders for ideal gases, we consult Flynn [2]. A schematic is presented in Figure (7). The gas first enters at lower speed at the inlet (point 1) before it is accelerated at the nozzles (point 2), then working its way through the rotor into the exhaust (point 3). The gas typically enters via the inlet and is directed tangentially by the turbine nozzles, accelerated by the peripheral ends of the rotor blades. The gas gradually has to work its way up to the exhaust, against the centrifugal force it experiences inside the rotor.
Figure 7: Schematic of a turboexpander.
For a loss-free turbine operating at sonic speed vc, the kinetic energy EK per mole of gas leaving
the nozzle is EK = 1 2M v 2 c, (44)
where M is the molar weight. By energy conservation, the work to expel one mole from the nozzle is the difference in the specific enthalpy: h1− h2 = ∆h = M v2c/2 [2]. We also recall that
∆h = cP∆T = cP(T1− T2), (45)
where cP is the heat capacity at constant pressure. The next step in the turboexpander is for
the gas to overcome the centripetal potential energy barrier between point 2 and 3. We know that the centripetal force F is expressed as F = M v2
c/r. Employing the relation vc = ωcr, we
W = Z r=R r=0 F dr = Z r=R r=0 M ω2crdr = 1 2M ω 2 cR2 = 1 2M v 2 c, (46)
assuming constant angular velocity of the turbine. The work done by the gas has to originate from somewhere for the gas to work against the centrifugal force. Thus, the deposited energy in the shaft of the rotor is once again the difference in the specific enthalpy: h2− h3 = M v2c/2.
This is true if the exhaust velocity is very small. Thus, the energy absorbed from one mole of gas at sonic velocities is M v2c/2. With the same argument as in Equation (45) between point 2 and 3, we obtain in total for the whole turboexpander process between point 1 and 3
cP(T1 − T3) = cP(T1− T2) + cP(T2− T3) = 1 2M v 2 c + 1 2M v 2 c = M v 2 c. (47)
For an ideal gas, we can take this reasoning one step further. We use the well-known thermo-dynamic relations for isentropic processes
P Vγ = constant ⇒ P1 P2 = V2 V1 γ , (48)
and for the speed of sound in an ideal gas
v2c = γP V
M , (49)
where γ = cP/cV [2]. From the equation of state for one mole of an ideal gas in Equation (2),
we obtain between point 1 and 2
V1
V2
= P2T1 P1T2
. (50)
From Equation (50) in (48), we get
P1 P2 = T1 T2 γ/(γ−1) , (51)
which can also be applied between 1 and 3. Plugging in the ideal gas law Equation (2), together with (47) and (49), we see that
T1 T2 − 1 = 1 − T3 T2 = γR 2cP ⇒T1 T3 = 1 + γR 2cP 1 − 2cγR P . (52)
For a monoatomic ideal gas such as helium, cP = 52R and γ = 5/3, so from Equation (52) and
T1
T3
= 2 and P1
P3
= 5.64, (53)
where 1 and 3 are the initial and final stages of the turboexpander. In most cases, it is essential that the gas does not undergo too much cooling and liquefy, which can damage the turbine. In addition, supersonic speeds can cause other unwanted effects that reduce the efficiency of the turboexpander. Here below is a code snippet, exemplifying how the CoolProp library can be used to describe a turboexpander for helium and calculate its output work ˙We, if the input
temperature T1, pressure P1 and mass flow ˙m are known. The PropsSI function provided by
CoolProp takes an input pair of properties and returns a third property from the equation of state. In this case, the input pair is P and T and the returned property is the specific enthalpy h, although it is capitalized in the code.
1 % T u r b o e x p a n d e r for h e l i u m 2 h1 = py . C o o l P r o p . C o o l P r o p . P r o p s S I ('H','P', P1 ,'T', T1 ,'Helium ') ; 3 T3 = T1 /2; 4 P3 = P1 / 5 . 6 4 ; 5 h3 = py . C o o l P r o p . C o o l P r o p . P r o p s S I ('H','P', P3 ,'T', T3 ,'Helium ') ; 6 d e l t a h = h1 - h3 ; 7 We = m * d e l t a h ;
3.4
Joule-Thomson Valve
An important phenomenon for liquefaction cycles is the Joule-Thomson (or throttling) process. It describes the temperature change of a real gas when it is forced through a porous plug or valve. This process is isenthalpic, meaning that ∆H = 0 and no heat is exchanged with its environment. For ideal gases, the Joule-Thomson (JT) expansion does not change P V , meaning that U is unchanged and the temperature is constant. For real gases however, the P V factor does change. The rate of change of temperature is described by the Joule-Thomson coefficient µJ T, defined as µJ T = ∂T ∂P H=const . (54)
In a gas expansion, ∂P is negative by definition. Hence, if µJ T < 0, the change in temperature is
positive upon Joule-Thomson expansion, and conversely µJ T > 0 means that the gas experiences
a temperature drop. The cooling produced by this effect is often used for the liquefaction of gases. The temperature at which µJ T changes sign is called the inversion temperature of a gas.
Several examples of µJ T at atmospheric pressure can be seen in Figure (8).
At room temperature, only hydrogen, helium and neon increase their temperature upon expan-sion, and require lower temperatures for throttling to act as cooling. Helium has the lowest inversion temperature of all gases, around 40 K depending on the pressure [2]. Therefore, in order to achieve liquefaction of helium, temperatures below this limit are needed.
If the initial temperature T1, the pressure P1 before and the pressure after P2 are known, we
can easily find the final temperature T2 as we know the isenthalpic character of the throttling
Figure 8: Joule-Thomson coefficient µJ T values for some gases at atmospheric pressure. The
inversion temperature for each gas at this pressure is where the x-axis is intersected. Image courtesy: Hankwang, Wikimedia Commons.
the PropsSI function provided by CoolProp, with P1 and T1 as the input pair. As we know
that enthalpy and mass flow are conserved, we simply use PropsSI again with h2 = h1 and P2
as the input pair.
1 % JT - v a l v e - i s e n t h a l p i c p r o c e s s
2 h1 = py . C o o l P r o p . C o o l P r o p . P r o p s S I ('H','P', P1 ,'T', T1 ,'Helium ') ; 3 h2 = h1 ;
4 T2 = py . C o o l P r o p . C o o l P r o p . P r o p s S I ('T','P', P2 ,'H', h2 ,'Helium ') ;
4
Classical Liquefaction Cycles
In this section, some of the classical liquefaction cycles and the theory behind are presented by gradual build-up. They also include some results from simulations of the theoretical model of Claude and Collins, with example input values of e.g. P and T similar to those in the FREIA Liquefier.
4.1
Linde-Hampson Cycle
pressures. However, the inlet temperature of the JT valve needs to be even lower for liquefaction to take place. In the steady state, a portion of the gas in 3 is then liquefied and collected from the phase separator. The phase separator is the device in the liquefaction cycle preceded by the Joule-Thomson expansion valve, typically operating at lower pressures. If the inlet temperature in the throttling process is low enough, the temperature reduction will be enough to liquefy part of the gas.
Figure 9: Illustration of the Linde-Hampson cycle with initial mass flow ˙m and output fluid mass flow ˙mf after point 3.
What amount of liquid comes out of the process? We define the yield y as the outgoing fluid mass flow ˙mf divided by the total mass flow ˙m. The return gas in 4 is then heated up
again by the counterflow heat exchanger and enters the compressor at 5. In order to derive y analytically, we set up a control volume confined to the dashed cyan box. Applying the first law of thermodynamics, the sum of the ingoing enthalpies is equal to the outgoing enthalpies of the volume
˙
mh1 = ˙mfhf + ( ˙m − ˙mf)h5, (55)
where hi is specific enthalpy and Hi = ˙mihi is the total enthalpy, for i = 1, 5, f (f stands for
“fluid”). Solving for y we get
˙ m(h1 − h5) = ˙mf(hf − h5) ⇒ y ≡ ˙ mf ˙ m = h1− h5 hf − h5 (56)
Equation (56) tells us that the proportion of the total mass flow ˙m is simply a function of the specific enthalpies at various points. In addition, this crucial equation states exactly how much of the incoming mass flow in the phase separator that is locally turned into liquid. Below the critical temperature, we are within the saturation dome with mixed phases. The main purpose of the phase separator is to discharge all remaining gas back into the liquefaction cycle and deliver the liquid separately into its cryogenic storage dewar, a well-isolated vessel for the liquid. The following code snippet illustrates how we use y to simulate what happens inside the phase separator. If T < Tcritical, parts of the gas will be turned into liquid. We use the
gas hgas. The yield y is found from Equation (56), and the total gas enthalpy Hgas is simply
Hin− Hf = Hin− y ˙mhf, as we know that the outgoing liquid mass flow ˙mf = y ˙m.
1 % P h a s e s e p a r a t o r 2 Hin = m * hin ; 3 if Tin < 5 . 1 9 5 3 % c r i t i c a l T of h e l i u m 4 hf = py . C o o l P r o p . C o o l P r o p . P r o p s S I ('H','T', Tin ,'Q',0 ,'Helium ') ; % Q = v a p o u r q u a l i t y 5 h g a s = py . C o o l P r o p . C o o l P r o p . P r o p s S I ('H','T', Tin ,'Q',1 ,'Helium ') ; 6 y =( hgas - hin ) /( hgas - hf ) ;
7 y = min (1 , max (0 , y ) ) ; % to e n s u r e 0 < y < 1 8 H g a s = Hin - y * m * hf ; % t o t a l e n t h a l p y of o u t g o i n g gas 9 e l s e 10 H g a s = Hin ; 11 end
4.2
Claude Cycle
In 1902, Georges Claude further improved the Linde-Hampson cycle by adding a piston, even though turboexpanders are often used nowadays [2]. The Claude cycle combines isentropic and isenthalpic expansions to achieve liquefaction. As the gas performs work in the rotor of the turboexpander isentropically, greater pre-cooling is achieved and also a higher yield y. The Claude cycle is displayed in Figure (10), with control volume marked in cyan. The Claude cycle relies on the same principles as the Linde-Hampson cycle, with the important addition of the isentropic expansion in the turboexpander (TXP). Some of the gas between point 2 and 3 from the main path is led through the turboexpander, and then returned in the mixer between point 8 and 9, where they merge again. This deviation of gas from the main flow is also illustrated in step 5 to the right in Figure (2). We also have three counterflow heat exchangers instead of one, as in the Linde-Hampson cycle.
Figure 10: Illustration of Claude cycle with initial mass flow ˙m entering at point 1, mass flow ˙
The output work of the turboexpander ˙We is the result of the isentropic expansion of the gas.
We define x as the fraction x = m˙e
˙
m , where ˙me is the mass flow through the expander and ˙m is
the total incoming mass flow at 1. The output fluid mass flow from the cycle after point 6 in Figure (10) is denoted ˙mf. Applying once again the first law of thermodynamics in the control
volume, we have that the ingoing enthalpy must sum up to the total outgoing enthalpy
˙
mh1 = ˙We+ ( ˙m − ˙mf)h11+ ˙mfhf. (57)
We take into account that the expander work output is ˙We = ˙meh2− ˙mehe = ˙me∆h, where he
is the specific enthalpy of the expanded gas that has passed through the turboexpander. We therefore have, solving for the yield y
˙ mh1 = ˙me∆h + ( ˙m − ˙mf)h11+ ˙mfhf ˙ mh1+ ˙me∆h − ˙mh11= ˙mf(hf − h11) ⇒y ≡ m˙ f ˙ m = h11− h1 h11− hf + x ∆h h11− hf (58)
as long as x + y < 1. If we compare Equation (58) to Equation (56) from the Linde-Hampson cycle, we observe that there is an extra term depending on the proportion x of the mass flow that enters the turboexpander, and a higher y can be achieved with respect to the Linde-Hampson cycle.
However, we need to be aware of what now happens in the phase separator. As the liquefaction occurring in the phase separator is purely local, the turboexpander flow x ˙m also reduces the mass flow into the JT valve and the phase separator to ˙m0 = (1 − x) ˙m. We therefore introduce the concepts local and global yield. By local, we mean at the point precisely after the JT valve, where liquefaction takes place on a molecular level. On the other hand, global means over the whole cycle, including the separate turboexpander flows. The global yield y2 for the Claude
cycle is the yield we saw from Equation (58). For the local yield however, we keep our focus on point 6 and 7 in Figure (10). By looking at the conservation of enthalpy and mass locally around the phase separator, we see that what comes in must also come out:
˙ m0h5 = ˙mfhf + ( ˙m0− ˙mf)h7 ⇒y1 ≡ ˙ mf ˙ m0 = h5− h7 hf − h7 (59)
Consequently, the local yield y1 is simply obtained from Equation (59), at first sight similar to
the case of the Linde-Hampson cycle in Equation (56). The difference is that Equation (59) depends on the specific enthalpies just around the phase separator and not at the inlet/outlet points of the whole cycle, as is the case in Equation (56). The amount liquefied in the phase separator is thus the local yield times the local mass flow: y1(1 − x) ˙m. No matter how
com-plicated the liquefaction cycle may be, y1 will always be of the simple form in Equation (59).
˙
m and ˙m0. As we later shall see for numerical iterative simulations in MATLAB, first solving for y1 locally is crucial to find the global yield y2 for the whole cycle.
Results from the simulations of the theoretical model can be seen in Figure (11a) and (11b), where the points on the x-axis correspond to the points in Figure (10). In Figure (11a), the temperature and gas flow profile of the Claude cycle is displayed. The top diagram shows the temperature T at that point, and the bottom diagram shows the gas flow Q at the same point. The gas flow Q has the same unit [kg/s] as mass flow ˙m and has identical meaning -not to be confused with vapour quality that is also de-noted Q! All the red boxes correspond to the heat exchangers, where T (the black line) gradually decreases or increases. The steep temperature drop occurs in the Joule-Thomson valve (dark-green box). The kink in the gas flow Q between point 6 and 7 in the cyan-coloured phase separator corresponds to the gas that has been liquefied. Step 12-13 correspond to the before and after the turboexpander (light-green box), respectively connected (marked by dashed lines) to the division of the flow between point 2 and 3, and the mixer between point 8 and 9 (small dark-blue box). The mixer is simply where the main gas flow and the turboexpander flow merge.
In Figure (11b), the local yield y1, the normalized local yield (1 − x)y1 and the global yield
y2 are shown for each iteration. The x-axis describes the iterative process towards equilibrium
conditions, and not the temporal evolution of the system. In particular, the initial spike is purely numerical. The normalized local yield and global yield actually overlap entirely to form one line, which they should if global and local calculations agree. Without any additional pre-cooling from ambient temperature, the global yield y2 is about 2.3% with a fraction x = 0.65
going into the turboexpander flow.
(a) Temperature and gas flow profile of the Claude cycle. (b) Global and local yields of the Claude cycle, where the close-up in the red box shows how the global and normalized local yield overlap.
4.3
Collins Cycle
The Collins cycle is an extension of the Claude cycle, but with two turboexpanders instead of one and two additional heat exchangers, as seen in Figure (12), with the control volume marked in cyan. The extra turbine and heat exchangers lead to further pre-cooling before liquefaction, meaning higher yields. It was invented by Samuel Collins in 1946, which revolutionized the commercially available production of liquid helium [14].
Figure 12: Illustration of the Collins cycle, with mass flow ˙m entering at point 1 and mass flows ˙
me1 and ˙me2 through the respective turboexpanders. The output fluid mass flow ˙mf occurs
after point 9.
The output work of the turboexpanders TXP1 and TXP2 are denoted ˙We1and ˙We2, respectively
with mass flows ˙me1 and ˙me2. Turboexpander flow 1 starts after point 2 and is reunited with
the main path again after point 14, whereas turboexpander flow 2 starts after point 5 and merges with the main path again after point 11. We also know that
˙
We1 = ˙me1∆he1
˙
We2 = ˙me2∆he2,
(60)
where ∆he1 and ∆he2 are the drops in specific enthalpy across the respective turboexpanders.
˙
mh1 = ˙me1∆he1+ ˙me2∆he2+ ( ˙m − ˙mf)h17+ ˙mfhf. (61)
Dividing by ˙m and rearranging, we get
˙ mf ˙ m (h17− hf) = h17− h1 + ˙ me1 ˙ m ∆he1+ ˙ me2 ˙ m ∆he2 ⇒y ≡ m˙f ˙ m = h17− h1 h17− hf + x1 ∆he1 h17− hf + x2 ∆he2 h17− hf , (62)
where x1 = m˙m˙e1 and x2 = m˙m˙e2. As can be seen in the final expression of Equation (62), we have
two additional turboexpander-related terms that contribute to increased global yield.
Results from the simulations of the theoretical model can be seen in Figure (13a) and (13b). In Figure (13a), the temperature and gas flow profile of the Collins cycle is displayed, where the points on the x-axis correspond to the points in Figure (12). The top diagram shows the temperature T at that point, and the bottom diagram shows the gas flow Q at the same point. All the red boxes correspond to the heat exchangers, where T gradually decreases or increases, just like for Claude but with two more heat exchangers. The steep temperature drop occurs in the Joule-Thomson valve. The kink in the gas flow between point 9 and 10 in the cyan-coloured phase separator corresponds to the gas that has been liquefied. Some of the gas from point 2.5 goes through TXP1 and is then returned to point 14.5. This separate mass flow is denoted step 18-19 in the graph. Similarly some of the gas from point 5.5 goes through TXP2, and is returned at point 11.5. This step is marked 20-21.
(a) Temperature and gas flow profile of the Collins cycle. (b) Global and local yields of the Collins cycle. Figure 13: Temperatures, gas flows and yields of the Collins cycle.
In Figure (13b), the local yield y1, the normalized local yield (1 − x1)(1 − x2)y1 and the
equilibrium conditions, and not the temporal evolution of the system. In particular, the initial spike is purely numerical. Also here, the normalized local yield and global yield actually overlap entirely to form one line, which they should if global and local calculations agree. Without any additional pre-cooling from ambient temperature, the maximum global yield y2 is about 5.2%
with fractions x1 = 0.45 and x2 = 0.35 of the gas flow going into the respective turboexpander
flows.
5
The FREIA Helium Liquefier
The liquefier Linde L140 in the cryogenic facility of the FREIA Laboratory has been installed by the Linde Group, and is capable of producing 140 l/h of liquid helium [8]. The liquefied helium is then collected in a 2000 l storage dewar. It is an important part of the research infrastructure to ensure a steady supply of coolants, needed for instance at the experiments in the low-temperature cryostats, but also at other institutions of the Uppsala University. A technical schematic of the system when not in operation is shown in Figure (14). The cycle starts with the compressor at the very top right, and then goes clock-wise in the schematic. High-pressure gas from the compressor is led through the first heat exchanger E3110, where the red loop corresponds to the liquid nitrogen pre-cooling system. The gas then goes through more heat exchangers and turboexpanders (central right), before reaching the JT-valve and then the phase separator (bottom left). The gas either goes in the return flow back to the compressor, or through the purification system located in the whole top left. Temperature sensors are the small boxes in red containing numbers, whereas the pressure sensors are the boxes coloured in blue.
In its structure, the FREIA liquefier is a hybrid of the Claude and Collins liquefaction cycles. A simplified thermodynamic cycle drawing is shown in Figure (15). Just like the Collins cycle, it has two turboexpanders for additional pre-cooling of the gas. However, the crucial distinction is that a single mass flow goes through both turboexpanders. As in the Claude cycle, there is only one separate mass flow going through the turboexpanders. In the case of the FREIA cycle, this intermediate step between the turboexpanders is connected to a three-fluid heat exchanger CX3. It is worth to point out that the turboexpander flow is parallel to the incoming hot flow on the high-pressure side, but anti-parallel to the outgoing cold flow on the low-pressure side of CX3. We denote Mixer 1 as the point between 2 and 3 that separate the mass flows, and Mixer 2 that reunite the flows between point 11 and 12.
Another important feature of the FREIA liquifier is the possibility of liquid nitrogen (LN2)
pre-cooling at around 77 K, connected to the three-fluid heat exchanger CX1. The mass flow of the liquid nitrogen over CX1 is ˙mN. As the nitrogen used as a coolant is heated up from
the hot flow, there is a specific enthalpy gain ∆hN. Differently to CX3, here both the LN2 flow
and the cold return flow are anti-parallel to the incoming hot flow. In addition, some part of the mass flow xpur =
˙ mpur
˙
m between point 6 and 7 may be discharged to cool the purification
system. This helium gas is not purified itself, but cool down impure helium gas that has been extracted from the recuperative helium system.
To derive the yield of the FREIA cycle, we once again look into the conservation of enthalpy in a control volume marked by dashed cyan lines in Figure (15). Just like for the Collins cycle, we assume a specific enthalpy drop ∆h1 and ∆h2 over TXP1 and TXP2 respectively, through
Figure 15: Simplified drawing of the FREIA liquefaction cycle. Mass flow ˙m enters at point 1, and mass flow ˙me goes through the turboexpanders. The output fluid mass flow ˙mf occurs
after point 9.
With respect to previous formulae for yield, Equation (63) contains a negative contribution from the purification process, and a positively contributing term from the LN2 pre-cooling.
However, we will neglect any purification processes for now as it decreases performance of the liquefaction and only needs to be done for small time intervals on a regular basis.
In order to simulate the LN2pre-cooling, we implement the theoretical model for the three-fluid
heat exchanger as described in Section 3.2.2. We model the flow of liquid nitrogen as a heat bath, absorbing as much enthalpy as possible from the hot flow coming from the compressor. A heat bath, or thermal reservoir, is a part of a thermodynamic system whose heat capacity C is large enough such that its temperature remains almost constant in contact with any other part of the system. The specific heat capacity hN of LN2 is known and can easily be found
with CoolProp, so to achieve a heat bath we assume that the mass flow of LN2 is 100 times
greater than the incoming hot flow of helium. How much enthalpy that is transferred across CX1 depends on the temperatures of the incoming flows, but also on the design of the heat exchanger. We can adjust the heat exchanger design parameters C10 and C20, that determine the coupling in the heat equations between the hot flow and the cold return flow, and between the hot flow and the LN2 flow respectively.
manual of the manufacturer. We know for instance that the compressor in FREIA is designed for an initial mass flow of around 41 g/s operating at 14 bar, but the mass flow is slightly lower during normal operation at 12 bar. We also know that the low-pressure side of the cycle has a pressure of 1.12 bar. There are several pressure and temperature sensors in the liquefier that can give us guidance, but far from everywhere. The global yield of the actual liquefier while running with LN2 pre-cooling is around 10%, whereas it drops to around 5-6% without.
Thus, there are several unknown parameters that we can adjust. First of all, we have to entirely guess the design of the heat exchangers by adjusting the design parameters C10 and C20 such that there is enough pre-cooling for liquefaction to occur. There are no mass flow sensors in the liquefier, so we do not know the turboexpander mass flow proportion x. The turbine expansion process was also assumed to be isentropic, which we do not know for the real liquefier. To imitate the real liquefier as well as possible, these degrees of freedom can be adjusted accordingly. Nevertheless, the theoretical model can also be used to simulate the hypothetical maximum efficiency of the liquefier.
5.1
Simulations and Results of the FREIA Model
Similar to the procedure of the other liquefaction cycles, we use the theoretical model based upon the assumption of enthalpy conservation in the cycle. A MATLAB function takes the input values of inlet pressure P1, gas flow Q1 and temperature T1. Q1 has the same unit
[kg/s] as mass flow ˙m and has identical meaning - not to be confused with vapour quality that is also denoted Q! This function then returns the temperatures and other parameter values at each point, which are iteratively used as input in a tester script for the next execution of the function. This process is then repeated until we have found numerical equilibrium, which we assume to correspond to the steady-state operation of the liquefaction cycle with no major changes in enthalpy flows. The temperature and gas flow profiles of the FREIA model simulations without and with LN2 pre-cooling are presented in Figure (16a) and (16b),
where the points on the x-axis correspond to the points in Figure (15), except step 17-20 that describe the turboexpander flow. The top diagram shows the temperature T at that point, and the bottom diagram shows the gas flow Q at the same point. Just like in the temperature profile of Claude and Collins, there is a gradual decrease in T through the heat exchangers, a sudden drop in T in the JT-valve and then gradual reheating again in the heat exchangers. Some of the gas is led from point 2.5 to the unified turboexpander flow, participating in the three-fluid heat exchanger CX3 between TXP1 and TXP2, then returning into the main flow at point 11.5. Once again, step 17-20 do not exist in Figure (15) but we include them to describe the temperature and gas profile of the turboexpander flow. A larger kink in the gas flow between point 9 and 10 can be seen for the model with LN2 pre-cooling, meaning that the yield is higher.
The yields of the FREIA model simulations without and with LN2 pre-cooling can be seen in
Figure (17a) and (17b) for each iteration. The local yield y1, the normalized local yield (1−x)y1
and the global yield y2 are shown for each iteration. The x-axis describes the iterative process
(a) Temperature and gas flow profile of the FREIA cycle without LN2
pre-cooling.
(b) Temperature and gas flow profile of the FREIA cycle with LN2
pre-cooling.
Figure 16: Temperature and gas flow profiles of the FREIA liquefier model simulations.
for these results, can be found in Appendices 11.1 and 11.2.
(a) Yields of the FREIA cycle without LN2 pre-cooling. (b) Yields of the FREIA cycle with LN2 pre-cooling.
Figure 17: Yields of the FREIA liquefier model simulations.
around order of 10−12 J/s and ∆Q stays flat at 0 kg/s. Therefore, we can conclude that both the enthalpy and gas flow are conserved over the cycle. This difference acts as a sanity check to ensure that no enthalpy or mass flow “disappears” in the cycle.
Figure 18: In/out difference in gas flow ∆Q and enthalpy ∆H in the FREIA model. A temperature-entropy (T S) diagram of the results with LN2 pre-cooling in Figure (16b) and
main path in Mixer 2. Point 9 to 16 represent the return flow of the gas, which is reheated isobarically. Compare the T S diagram of FREIA in Figure (19) to the liquefaction cycle sketch to the right in Figure (2).
Figure 19: T S diagram of the FREIA cycle with LN2 pre-cooling. The main path of the gas in
the FREIA cycle is marked out, with all points in Figure (15).
We also compared our simulation results to the actual sensor values inside the FREIA liq-uefaction cycle. The sensor values of the turboexpander path during liqliq-uefaction with LN2
pre-cooling can be seen in Figure (20a) and (20b). The average inlet temperature to TXP1 under operation was 40 K, and the average outlet temperature after TXP2 was 9.3 K, whereas our simulation results in Figure (16b) gave 49 K and 7.4 K respectively. However, we did observe a big discrepancy in the pressure drop over TXP1 in FREIA. The average value of P went from 12 bar to 5.8 bar over TXP1, a reduction by slightly more than a factor 2. For the isentropic expansion, we had assumed a pressure reduction of a factor 5.64. On the other hand, the average value of P decreases by a factor 5.18 over TXP2 if we consider the P = 1.12 bar of the return flow. This decrease factor is closer to our estimate. The sensor values, especially in pressure, wiggled in a periodic way. For instance, there is a clear drop at around 24 minutes into the cycle for the pressure sensors PI3130 and PI3150. This is due to the fact that a small portion of the gas is regularly used to cool the purification circuit, and the large drop in pressure occurs when this purification has finished. For the comparison with our simulation results, we took the average sensor values for TXP1 and TXP2 above.
In our calculations for the turboexpanders in FREIA, the total difference in entropy over both turboexpanders is ∆S = 2900 [J/(kg K)]. If we assume isentropic expansions in both turboexpanders, the enthalpy transfer ∆ ˙H from the turbine flow to the hot flow in CX3 needs to be 1660 J/s, but in our simulation this enthalpy transfer was 1300 J/s at steady state. This discrepancy means that TXP1 probably adds some entropy to the process, so it is very likely that TXP2 does so too. We can thus state that the real turboexpanders are not entirely isentropic, but almost.
(a) Sensor data of pressures before TXP1 (PI3130) and before TXP2 (PI3150).
(b) Sensor data of temperatures before TXP1 (TI3130) and after TXP2 (TI3155).
Figure 20: Example data from sensors in the FREIA liquefier during operation, here from the turboexpander section.
We found a maximum yield of 17.6% by setting the fraction x = 0.62 into the turboexpander flow and add a strong coupling between the hot incoming flow and the cold return flow in CX1 (i.e. increase design parameter C10), and slightly weakening the coupling between the hot flow and the LN2 pre-cooling (i.e. decrease C20). We also increased the design parameter C0 for all
the other heat exchangers. We found that the maximum global yield y2 = 17.6%. This was
done by varying the degrees of freedom until the maximum value of y at equilibrium was found. The exact parameter values can be found in Appendix 11.1 in the commented code, right next to the three-fluid heat exchangers.
The final results also indicated that the helium did not always behave like an ideal gas. De-pending on the pressure, the heat capacity at constant pressure cP is no longer constant below
roughly 20 K, as we saw in Figure (1). In our theoretical model for FREIA, this affects in particular the last two heat exchangers: CX4 and CX5. Surprisingly enough, the ideal gas approximation also holds for the liquefaction temperature at point 10 in Figure (15). However, the only point where cP actually is relatively far off from the ideal gas approximation is on the
high-pressure side at point 7. This deviation is about 33% from cP = 52R. On the other hand,
this effect can be compensated for in the design of the heat exchanger, as the heat exchanger design parameter C0 and length L are unknown and can be chosen freely in the model.
6
Discussion
In this study, a theoretical model of the FREIA helium liquefaction cycle has been constructed on the same foundations as the classical liquefaction cycles of Linde-Hampson, Claude and Collins. Despite the simplicity of enthalpy conservation, the interplay of the numerous thermo-dynamic components in the FREIA model is not always trivial. In the theoretical model, we had to ensure that not only enthalpies must be conserved, but especially that the processes were reasonable and physical: heat should not flow from cold to hot only by itself, the gas should not be liquefied inside the turboexpanders, the model must start with room temperature, and so on. The list of physical demands can be made long. The simulations of the theoretical model must converge to steady-state values at the equilibrium of the liquefaction process.
Surprisingly enough, it turns out that a major part of the interactions, especially in the heat ex-changers, are linear and can be solved relatively easily. As was shown in Section 3.2 and 3.3, the relation between the inlet and outlet temperatures of all heat exchangers and turboexpanders are purely linear. Concerning the theoretical background of the heat equations, we consulted Flynn [2] for the two-fluid heat exchangers and Sekulic et al. [12] for the three-fluid heat ex-changers. Flynn was also consulted for the turboexpanders. Our approach of solving systems of linear first-order homogeneous differential equations differed from the Effectiveness - Num-ber of Transfer Units (NTU) method used by Flynn and Sekulic, but the results of both heat exchanger types are identical. Our solution with eigenvectors and eigenvalues with boundary conditions turned out to be more mathematically digestible, with a consistent add-on strategy of the matrix describing the heat equations, no matter the number of flows. In addition, we can easily generalize our heat exchanger model to any three-dimensional shape by replacing our heat exchanger design parameter C0 with the heat transfer coefficient times the contact area U A, as mentioned in Section 3.2.
The three-fluid heat exchangers may result in long analytical expressions, but the heat equations governing the interactions are nonetheless linear with a unique solution. In all cycles, the only non-linear interaction is the liquefaction process in the phase separator. These are governed by the equations for the yield y, that in turn depends on the specific enthalpies hiat various points.
At this stage, CoolProp [9] is an invaluable tool for easily switching between thermodynamic quantities even when they are not linearly related. The thermodynamic cycles as entities are thus non-linear, but this non-linearity is very weak as most components are in fact linear with a unique solution. Therefore, the convergence to a steady-state liquefaction should be unique with no other stable points. Nevertheless, this reasoning of uniqueness is important to keep in mind, especially since our theoretical model is solved numerically.
values, such as the heat exchanger design parameter C0, we can swiftly change the output yield as desired. From the configuration presented in Figure (16b), (16a), (17b) and (17a), our simulations result in slightly more than 10% with LN2 pre-cooling, whereas it drops to 6.1%
by simply removing the coupling to the LN2 flow. We used 12 bar for the high-pressure side
in the simulations, just like the real liquefier does. However, the cycle inlet gas flow Q1 of
41 g/s of gaseous helium from the compressor manufacturer’s guide is specified for 14 bar, a slightly higher pressure than at which the real liquefier in FREIA operates. Thus, the real gas flow is slightly lower than 41 g/s. As no gas flow sensors are installed in the liquefier, we do not know this quantity. However, by operating the liquefier compressor at different pressures and investigate the liquid output, one could in theory interpolate what initial mass flow ˙m corresponds to 12 bar. In addition, the fact that the real turboexpanders are not entirely isentropic, as we assumed in the model, may also contribute to discrepancies. Nonetheless, the yields of the simulations could be adjusted from design parameters and agree with the real liquefier, within less than one percentage point.
Not only did the simulations agree with real liquefaction yields, but also demonstrated an upper limit y = 17.6% without making the processes non-physical or technically infeasible. We observed that the major improvement in yield could be achieved by setting the turboexpander flow fraction x of the initial mass flow to 0.62. This also includes adding a strong coupling between the hot incoming flow and the cold return flow in CX1, and slightly weaken the coupling between the hot flow and the LN2 pre-cooling. We also increased C0 for all heat exchangers.
Even though the real liquefier does not reach 17% yield, this may act as a guidance for what to focus on in improving the performance of the machine. Even though heat exchangers are expensive to replace, it may also indicate where possible future hardware upgrades of the liquefier may have the largest impact.
Throughout this study, several approximations of real physical phenomena had to be made. In our theoretical model, we assumed isentropic expansion of monatomic gas in the turboex-panders. In reality, the turbines in use inside the FREIA liquefier are not entirely isentropic. Due to the lack of sensors at some points in the liquefier, we cannot verify this. As can be seen in the schematic in Figure (14), temperature sensors are not placed anywhere between CX2 and CX4, nor between the turboexpanders. This means that even though we know the difference in entropy between the inlet to TXP1 and the outlet of TXP2, we do not know how much enthalpy is transferred in CX3. For the turboexpanders in FREIA to be isentropic, more enthalpy had to go through the heat exchanger CX3 than allowed for in our simulations, which shows that they are not fully isentropic. In order to investigate this further, we would need more detailed manufacturing information and design.
Another approximation was to assume that helium behaves like an ideal gas, with specific heat capacity at constant pressure cP = 52R. As we saw in Figure (1), this approximation holds
for higher temperatures and precisely at liquefaction temperature inside the FREIA liquefier. The sole large discrepancy from the ideal gas approximation occurred at point 7 in Figure (15), but this effect is small and can be considered in the heat exchanger design, as the heat exchanger design parameter C0 is unknown. The important principle of enthalpy conservation remains. For instance, a small increase in cP for the return flow simply means that C0 decreases
slightly. As we can adjust C0 freely in the model, this effect is probably due to numerics. The main flaw of the model is probably the fact that we assume cP to be constant inside the heat
exchangers, which we know is not the case for temperature intervals spanned by CX5. However, this correction might belong to the higher order.
in the introduction, the market size of cryogenics is rapidly growing as more applications are discovered. Medical equipment, organ transplantation, food storage, power transmission, super-conductors, electronics and many interdisciplinary scientific experiments all require cryogenic temperatures. Liquid nitrogen at 77 K suffices in many cases, but many experiments that require lower temperatures must use liquid helium. To better understand the liquefaction pro-cess, not only of helium but also of other gases, is thus crucial for an efficient production. In addition, an optimized liquefaction cycle allows for electricity savings and lower maintenance costs. For instance, a slight change to the heat exchanger design may allow the compressor to work at lower power, or to reduce the amount of liquid nitrogen needed for cooling, while keeping high liquefaction yields. Theoretical models and simulations allow laboratories, big and small, to better understand and take the command of the thermodynamic processes in their own facilities.
Hopefully, these types of simulations may also lead to improved methods to increase the yield and maintain higher liquefaction rates. As we showed with the inclusion of LN2 pre-cooling,
a great increase in yield was obtained as the helium gas could be cooled even further before reaching the Joule-Thomson valve, where the temperature-dependent liquefaction process oc-curs. Adjusting the design heat exchanger CX1 that contained the pre-cooling suddenly led to yet another increase in yield. Before changing the hardware of the liquefier itself, it is faster and less expensive to optimize the design in a proper theoretical model. Benchmarking a first version of cryogenic code for the FREIA liquefier is hopefully a step in the right direction.
7
Future Prospects
There are several aspects of the model that can be improved or extended. As we saw for the heat capacity at constant pressure cP of helium, it is not constant for lower temperatures.
This affected mostly the high-pressure inlet at the last heat exchanger, and could probably be considered in the heat exchanger design. In our theoretical approximation for the heat exchangers, we assumed cP as constant and that of an ideal gas. Even though this would
be more accurate for the last heat exchanger, a theoretical model with non-constant cP would
probably only provide a higher-order correction and be very complicated. It would be interesting to investigate how the Linde Group deals with heat exhangers for helium spanning over this domain of irregular cP and if the design is drastically different.
Another outlook for our theoretical model is to install further sensors in the FREIA liquefier, especially for mass flows and temperatures between the turboexpanders, still a terra incognita. For the maximum yield of liquefaction, we saw that a slight change in the heat exchanger design parameter C0 led to an increased yield of several percentage points. This was especially true for the third heat exchanger in FREIA connected to the turboexpander flow. If possible, after theoretical optimization with our model, one could try to adapt this design to see if the actual liquefaction increases.
