# Fluid mechanical response of a pulse tube cryocooler: modelling and experimental validation

Kévin Lauzier, Federico Bribiesca Argomedo, Jean-Yves Gauthier, S Sesmat, X Lin-Shi, D Carvalho Lopes

► To cite this version:

Kévin Lauzier, Federico Bribiesca Argomedo, Jean-Yves Gauthier, S Sesmat, X Lin-Shi, et al., Fluid mechanical response of a pulse tube cryocooler: modelling and experimental validation. IOP Conference Series: Materials Science and Engineering, 2020, Advances in Cryogenic Engineering: Proceedings of the Cryogenic Engineering Conference (CEC) 2019, 21-25 July 2019, Hartford, Connecticut, USA, 755, pp.012012. 10.1088/1757-899X/755/1/012012. [hal-02934490](https://hal.science/hal-02934490)

HAL Id: hal-02934490

https://hal.science/hal-02934490

Submitted on 9 Sep 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
PAPER • OPEN ACCESS

# Fluid mechanical response of a pulse tube cryocooler: modelling and experimental validation

To cite this article: K Lauzier et al 2020 IOP Conf. Ser.: Mater. Sci. Eng. 755 012012

View the article online for updates and enhancements.
# Fluid mechanical response of a pulse tube cryocooler: modelling and experimental validation

K Lauzier¹,², F Bribiesca-Argomedo², J-Y Gauthier², S Sesmat², X Lin-Shi² and D Carvalho Lopes¹

¹Space business unit, Air Liquide Advanced Technologies, Sassenage, 38360, France

²University of Lyon, INSA Lyon, CNRS, Ampère, Villeurbanne, F-69621, France

E-mail: kevin.lauzier@airliquide.com

**Abstract.** Earth observation satellites require cryocoolers to cool down their infrared imagers at very low temperatures. Besides having very good thermodynamic performances, satellite cryocoolers are expected to generate as little vibrations as possible. In order to better understand vibration causes between 50 and 500 Hz, a precise model of the whole cryocooler is necessary. In the literature, two main modelling approaches for pulse tube cryocoolers exist: compressor-oriented models reduce the thermodynamic system to a linear mass-spring-damper system acting on compressor's pistons; and thermodynamically oriented models aimed at understanding and predicting thermodynamic performances.

In this paper, the mechanical behavior of the thermodynamic system is modelled. Assumptions concerning gas properties and thermodynamic behavior are made based on SAGE software simulations. A simplified Redlich Kwong equation is used. When necessary, polytropic coefficients were identified using simulation data, as were the time-averaged temperatures. The thermodynamic system is split into volumes, pipes and regenerators: conservation laws in volumes are integrated, dynamic mass and momentum conservation equations in pipes are solved using the method of characteristics and equations in regenerators are solved using a finite difference method. Three friction laws are used: one for straight pipes (Moody chart), another for wound pipes (White's correlation) and a last one for porous media (Modified Ergun equation). Porosity is measured by weighing. The model is built to respect integral causality and propagation phenomena.

The developed model is validated using experimental data. The simulations highlight the non-linear mechanical behavior of the thermodynamic system: a sinusoidal motion of the pistons induces a non-sinusoidal pressure in the compression chamber. Among other, first and second harmonics amplitudes are about 3% of fundamental pressure amplitude.

This model can now be integrated into a global cryocooler model to predict compressor's vibrations, power consumption or electrical harmonics. It could also be extended to predict vibrations from the thermodynamic system.

## 1. Introduction

For Earth observation purposes, satellites carry infrared cameras. In order to minimize the picture's noise, focal planes must be cooled down at very low temperatures: typically between 150 K and 50 K. Pulse tube cryocoolers are among the most commonly used cryocoolers for such missions. Some of the advantages of pulse tube technology are its reliability and low vibration level. For each harmonic of the 50 Hz driving frequency, up to 500 Hz, vibrations have to meet stringent requirements in terms of amplitude. Since the only moving parts are the compressor's pistons, they are at the origin of most of

<figure><img src="image_1.png" /></figure>
the vibrations in the system. Those vibrations are mostly due to compressor asymmetries and are affected by non-linearities and the frequency response of the whole system. In order to understand and predict vibrations, each phenomenon acting on pistons must be modelled properly in the frequency range 50-500Hz. Thus, particular care must be taken to model the interaction between the compressor and the thermodynamic system.

In the literature, there are two main modelling approaches for pulse tube cryocooler. First, purely compressor-oriented models reduce the thermodynamic system to a linear mass-spring-damper system [1] [2]. This leads to a linear first order frequency behavior, which is reductive and does not adequately represent the harmonic content of the measured forces. Second, more thermodynamically oriented models aimed at understanding and predicting thermodynamic performances: efficiency, thermal behavior and available cooling power. Since pulse tube technology has been invented in the 1960’s [3], thermodynamic improvements were mostly performed through experiments until the 1990’s, when the first numerical models started to appear for this technology [4] [5]. Nowadays, a lot of different models are available and SAGE software, developed by Gedeon [6], comes out of those works. More recently, system-oriented models have been developed to study the performances of the thermodynamic system coupled with a compressor but they do not take harmonics into account and assume a linear behavior [7].

In this paper, an intermediate model based on conservation laws is developed. It is not aimed at predicting thermodynamic behavior or efficiency of the system, but rather reproduce the mechanical interaction between the compressor and the thermodynamic system at fundamental frequency and harmonics. The achieved model is validated using experimental data.

The paper is organized as follows: the whole system will first be described in section 2, the model establishment will next be explained in section 3, the results of experimental validations will be presented in section 4 and finally some conclusion and perspectives will be given.

## 2. System description

A cryocooler unit (Figure 1) is composed of a cold finger assembly (CFA), a compressor part assembly (CPA) and a control electronics with sensors.

<figure><img src="image_2.png" /><figcaption>Figure 1: Picture of the cryocooler unit. The thermodynamic system is inside the CFA and CPA as shown in Figure 2</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>Figure 2: Limits of the studied system. Its components are numbered as follow: 1:Compression chamber; 2:Split pipe; 3:Regenerator warm heat exchanger; 4: Warm regenerator; 5:Cold regenerator; 6:Cold heat exchanger; 7:Cold flow straightener; 8:Pulse tube; 9:Warm flow straightener; 10: Pulse tube warm heat exchanger; 11:Diffuser; 12:Inertance pipe; 13:Wound inertance; 14:Buffer.</figcaption></figure>

In the thermodynamic system, the compressor makes the fluid (helium) oscillate around approximately 30 bar. Lopes summarizes in [8] the history and the aim of each component of the thermodynamic system: thanks to the inertance and the buffer, there is a phase shift between the flow and the pressure in the pulse tube. This way, the gas relaxes adiabatically inside the pulse tube, so its temperature varies in the same way as the pressure, and evolves almost isothermally inside the regenerator as the latter has a very large specific heat exchange area and a strong heat capacity. The aim of the regenerator is also to store cold power in order to cool down the gas before entering in the pulse tube.

This present paper focuses on the fluid mechanical behavior of the thermodynamic system that is inside the CPA and the CFA as defined in Figure 2.

Nomenclature

E Pressure relative error (%)

L Component length (m)

M Molar mass (kg·mol⁻¹)

P Pressure (Pa)

R₁ First Riemann variable

R₂ Second Riemann variable

Re Reynold's number

S Pipe section (m²)

T Temperature (K)

U Voltage signal of the sensor (V)

V Volume (m³)

b Redlich Kwong molecular size const.

c Wave propagation speed (m·s⁻¹)

dₕ Hydraulic diameter (m)

e Specific energy (J·m⁻³)

k Polytropic coefficient

qₗ Volume flow rate (m³·s⁻¹)

r Specific gas constant (J·kg⁻¹·K⁻¹)

t Time variable (s)

u Uncertainty (same as the variable)

v Flow velocity (m·s⁻¹)

x Position variable (m)

γ Heat capacity ratio

ρ Density (kg·m⁻³)

ψ Darcy friction factor

Subscripts

0 Time average

A Component closer from compressor

B Component farther from compressor

c Combined

P Relative to pressure

t Partial derivative with respect to time

tot Total

U Relative to voltage

x Partial derivative with respect to position

x0 Inlet position

xL Outlet position
## 3. Modelling

### 3.1. Methodology and assumptions

To model the thermodynamic system of Figure 2, each component is identified either to a volume, or a line or a regenerator. A section change model is sometimes necessary to interconnect models. Models are then connected to comply with integral causality or propagation phenomenon. Figure 3 presents how the submodels are interconnected in the global thermodynamic model. In volumes, pressure is a space constant and section is wide enough to neglect friction for typical mass flows. In lines, section is small enough to consider friction and eventually propagation if it is long enough. The regenerators are considered as porous lines with a temperature gradient.

<figure><img src="image_4.png" /><figcaption>Figure 3: Submodels interconnection and causality variables in the global thermodynamic model.</figcaption></figure>

Thermodynamic hypotheses for each submodel are made based on SAGE software simulations. Time-averaged temperatures as well as equivalent polytropic coefficients are extracted.

### 3.2. Generic models

#### 3.2.1. State relation for the gas.
As the gas temperature varies between 293K and 50K for a pressure of about 30 bar, perfect gas state relation is not adequate. Starting from Redlich Kwong equation [6] [9] with adapted coefficients to helium [6] [10], it is possible to neglect intermolecular forces for the considered temperature range (the same result can be achieved using the Van Der Waals equation) in order to obtain:

$$ \frac{P}{\rho r T} = 1 + \frac{\rho \frac{b}{M}}{1 - \rho \frac{b}{M}} \qquad (1) $$

This simplified Redlich Kwong equation as well as initial Redlich Kwong equation and perfect gas state relation have been compared to values from SAGE tabulated gas properties. The pressure relative errors are presented in Figure 4 and enable to validate the relevance of the simplified Redlich Kwong equation in the considered operating conditions.
<figure><img src="image_5.png" /><figcaption>Figure 4: Pressure relative errors ($E = 100 \cdot (P_{equation} - P_{SAGE_{tabulated}})/P_{SAGE_{tabulated}}$). The thermodynamic system has been simulated with SAGE. The pressure $P_{SAGE_{tabulated}}$ at any point and any given time has been compared to state equations results for each given state ($\rho_{SAGE}$, $T_{SAGE}$).</figcaption></figure>

### 3.2.2. Conservation laws.
Equations (2) show the mono-dimensional adiabatic Euler equations for viscid fluid in a pipe of time variable section [6] [11]. It is considered that only pressure and friction forces are acting on the fluid.

$$ 
\left\{ 
\begin{aligned}
\partial_t(\rho S) &= -\partial_x(\rho v S) && (a) \\
\partial_t(\rho v S) &= -\partial_x(\rho v^2 S) - \partial_x(PS) - \frac{\psi S \rho v |v|}{2d_h} && (b) \\
\partial_t(eS) &= -\partial_x(veS) - \partial_x(PSv) - P \partial_t S && (c)
\end{aligned}
\right. 
\qquad (2) $$

From an energy point of view, it is possible to consider the gas as perfect. Sage tabulated gas properties reveal this hypothesis leads to an internal energy maximal error of 2 %. The total specific energy e can thus be expressed as a function of pressure, density and flow velocity [12]:

$$ e = \frac{1}{2} \rho v^2 + \frac{P}{\gamma - 1} \qquad (3) $$

Equations (2) can then be expressed using $\rho$, $v$ and $P$ as follows:

$$ 
\left\{ 
\begin{aligned}
\partial_t(\rho S) + \partial_x(\rho v S) &= 0 && (a) \\
\partial_t v + \frac{1}{2}\partial_x v^2 + \frac{1}{\rho}\partial_x P + \frac{\psi v |v|}{2d_h} &= 0 && (b) \\
\partial_t \ln(PS^\gamma) + v\partial_x \ln(P) + \gamma\partial_x v + \frac{(1-\gamma)\psi\rho|v|^3}{P} &= 0 && (c)
\end{aligned}
\right. 
\qquad (4) $$

### 3.2.3. Polytropic behavior.
For some components, it is possible to make additional assumption concerning the thermodynamic behavior. This means replacing equation (4.c) with an algebraic relation between $P$ and $\rho$ called the polytropic equation:

$$ P = \frac{P_0}{\rho_0^k} \rho^k \qquad (5) $$

When an adiabatic process is assumed, the polytropic coefficient $k$ is set to the gas heat capacity ratio $\gamma$. As for isothermal process, $k$ is set to one.

### 3.2.4. Darcy friction factor $\psi$.
Three different dependences of Darcy friction factor on Reynolds number are considered: one for smooth straight pipes, another one for smooth wound pipes and a last one for porous pipes. In the case of a smooth straight pipe, the Moody diagram is piecewise approximated [8]. As for smooth wound pipes, White correlation is used [13] [14]. Finally, Gedeon's
works on porous media friction factors for woven screen matrices are used in porous pipes [15]. These friction factors are plotted in Figure 5 as a function of Reynolds number.

<figure><img src="image_6.png" /><figcaption>Figure 5: Darcy friction factors used in this paper</figcaption></figure>

### 3.3. Variable adiabatic volume model

Considering pressure gradient and friction to be negligible in the energy conservation equation (4.c), pressure as a function of volume as well as incoming and outgoing flow rate is obtained in equation (6) whatever the temperature gradient. In addition, the change in density is supposed to be negligible when entering or outgoing the volume.

$$P(t) = P(0) \left(\frac{V(0)}{V(t)}\right)^{\gamma} \exp\left(\gamma \int_{0}^{t} \frac{q_{v_{x0}} - q_{v_{xL}}}{V(t)} dt\right) \quad (6)$$

Such a model complies with integral causality and represents, as expected, a capacitive storage.

### 3.4. Line model

To model lines of constant section, the method of characteristics [15] [12] is used. The matrix of the two first conservation equations (4.a) and (4.b) is diagonalized and polytropic equation (5) is considered to calculate Riemann's variables (7). The advantage of this method is to use propagation phenomenon to maximize the simulation step size and minimize the simulation time.

$$\left\{ \begin{array}{l} \left(\begin{pmatrix} R_1 \\ R_2 \end{pmatrix}_t + \begin{pmatrix} v+c & 0 \\ 0 & v-c \end{pmatrix} \begin{pmatrix} R_1 \\ R_2 \end{pmatrix}_x + \frac{\psi v |v|}{2d_h} \begin{pmatrix} 1 \\ 1 \end{pmatrix} = 0 \quad (a) \\ \\ \left(\begin{pmatrix} R_1 \\ R_2 \end{pmatrix}_t = \begin{pmatrix} v + \frac{2}{k-1}c \\ v - \frac{2}{k-1}c \end{pmatrix} \quad (b) \quad \text{with } c = \sqrt{\frac{k}{\rho}} \quad (7) \end{array} \right. \\ \text{It can be remarked there are two waves traveling in opposite direction with quasi-equal velocity: } c - v \text{ and } c + v \text{ with } v \ll c \approx \sqrt{rT} \approx 800 \text{ m/s.}$$

To solve equation (7), directional derivative is used, considering waves propagating at constant velocity c. The equations are discretized to be solved before recovering physical quantities ρ and v from Riemann's variables by using the relation:

$$Y = \begin{pmatrix} \rho \\ v \end{pmatrix} = \left( \left[ \left( \frac{k-1}{4} (R_1 - R_2) \right)^2 \frac{\rho_0^k}{k P_0} \right]^{\frac{1}{k-1}} \frac{R_1 + R_2}{2} \right) \quad (8)$$

To respect propagation phenomenon, causality of such lines will be on Riemann's variables whenever possible as shown in Figure 3. When there is a section change from a component to another, it is assumed the transition occurs in an incompressible way so that there is no volumetric mass change
and volumetric flow is conserved. With those assumptions, a section change from $S_A$ to $S_B$ corresponds to the following change in Riemann's variables.

$$ egin{cases} R_{2A} = 2R_{2B} \frac{S_B}{S_A + S_B} + R_{1A} \frac{S_B - S_A}{S_A + S_B} & (a) \\ R_{1B} = 2R_{1A} \frac{S_A}{S_A + S_B} + R_{2B} \frac{S_A - S_B}{S_A + S_B} & (b) \end{cases} \qquad (9) $$

When causality cannot be on Riemann's variables, such as for a transition between a line and a volume, inductive behavior will be assumed for causality: pressure must then be imposed to the component as it appears in Figure 3.

### 3.5. Regenerator model

Because of temperature gradient, methods of characteristics is not easily applicable. Therefore the two first conservation equations (4.a) and (4.b) are simply discretized with Euler-type scheme and isothermal relation (5) is used. Temperature is considered to be a constant at any given position and has been obtained with a SAGE simulation.

Regenerators are supposed to be mostly resistive components so causality can be chosen. Pressure will be imposed on one side and volumetric flow on the other side.

## 4. Experimental validation

### 4.1. Experimental protocol

The compressor is powered with two identical sinusoidal voltages at 57.5 Hz of constant amplitude. During the whole experiment, time-averaged temperatures as well as cooling power are constant. Pressure is measured at five locations as shown in Figure 2 to validate the model. The present study focuses on sensor P1 to analyze the forces exerted by the thermodynamic system on the pistons. Sensor P1 is a Keller PAA-25S. Its measurement range is 0-50 bar and it has a bandwidth of 5 kHz. Its performances as well as the ones of its conditioner were certified. Data from those certificates enable to determine the relationship (10) between relative pressure P and sensor voltage signal U considering pressure measurement uncertainty $u_p$ and voltage conditioner uncertainty $u_U$.

$$ P = 263(U \pm u_U)^2 + 497551(U \pm u_U) + 2535 \pm u_p \qquad (10) $$
with $u_p \le 600 Pa$ and $u_U \le 1.1 mV$

As the considered uncertainties are uncorrelated, measurements are to be considered with a combined uncertainty of $u_c \le 814 Pa$.

### 4.2. Comparison

<table><caption>Table 1: Spectral content of pressure P1 in the compression chamber predicted by the model compared to experimental data and to the output of a simpler linear mass-spring-damper model. (Gray background highlights experimental data inside accuracy range)</caption><thead><tr><th>H#</th><th>Freq.<br>(Hz)</th><th>Experimental<br>P1 (bar)</th><th>Simulated<br>P1 (bar)</th><th>Difference<br>experiment/simulation<br>(%)</th><th>Simulated<br>H#/H1<br>(%)</th><th>Mass-spring-<br>damper system<br>P1H#/P1H1(%)</th></tr></thead><tbody><tr><td>1</td><td>57.5</td><td>2.732 ± 0.008</td><td>2.789</td><td>+2</td><td>100.00</td><td>100</td></tr><tr><td>2</td><td>115</td><td>0.083 ± 0.008</td><td>0.083</td><td>+1</td><td>2.98</td><td>0</td></tr><tr><td>3</td><td>172.5</td><td>0.064 ± 0.008</td><td>0.066</td><td>+3</td><td>2.37</td><td>0</td></tr><tr><td>4</td><td>230</td><td>0.012 ± 0.008</td><td>0.011</td><td>-10</td><td>0.40</td><td>0</td></tr><tr><td>5</td><td>287.5</td><td>0.013 ± 0.008</td><td>0.007</td><td>-41</td><td>0.27</td><td>0</td></tr><tr><td>6</td><td>345</td><td>0.006 ± 0.008</td><td>0.003</td><td>-43</td><td>0.12</td><td>0</td></tr><tr><td>7</td><td>402.5</td><td>0.003 ± 0.008</td><td>0.003</td><td>-6</td><td>0.10</td><td>0</td></tr><tr><td>8</td><td>460</td><td>0.001 ± 0.008</td><td>0.002</td><td>+47</td><td>0.07</td><td>0</td></tr></tbody></table>

The model was run considering the sinusoidal motions of the pistons as input. Amplitude of stroke was obtained using a numerical model of the compressor, by entering the effective coil voltage measured
during the test. These pistons strokes applied to the fluid model of Figure 3 give a pressure swing with harmonics in the compressor chamber. These harmonics amplitudes are reported and compared to experiment in Table 1.

Simulated harmonics are thus very close from measured ones despite measurements inside accuracy range for higher order harmonics. It shows that the mechanical behavior of thermodynamic system is not linear: when excited with a sinusoidal piston motion, it induces harmonics in the compression chamber's pressure. First and second harmonics amplitudes are about 3% of fundamental amplitude.

**5. Conclusion**

The new model proposed in this paper focuses on the mechanical behavior of the work-fluid in a pulse tube cryocooler. Thermodynamic, thermal and gas properties hypotheses are validated using SAGE software simulations. It uses conservation equations, a simplified Redlich Kwong state equation as well as polytropic equation. Whenever applicable, method of characteristics is used to solve conservation equations in order to lower simulation time. The model is built to respect integral causality and propagation phenomena.

The model reproduces the pressure fundamental and first 2 harmonics very accurately (3% or better, with regard to measurement). 4th and 5th harmonics are also given within the confidence interval. It highlights the non-linear mechanical behavior of the thermodynamic system. Excited with sinusoidal pistons positions, pressure in the compression chamber appears to be non-sinusoidal. Among other, first and second harmonics amplitudes are about 3% of fundamental amplitude.

Such model can now be easily integrated in a whole cryocooler model in order to predict power consumption, vibrations (especially to see how these non-linearities interact with those induced by the CPA), or electrical harmonics. Hopefully, the last two capabilities will be of great help to improve the present force cancellation strategies based on dissimilar voltages or currents.

**6. References**

[1] Dadd M W, Bailey P B, Davey G, Davis T and Thomlinson B J 2002 *Cryocoolers 11* (Springer, Boston, MA) pp 175–82

[2] Sobol S and Grossman G 2015 *IOP Conf. Ser. Mater. Sci. Eng. 101* 012093

[3] Gifford W E and Longsworth R C 1964 *J. Eng. Ind. 86* 264

[4] Storch P J, Radebaugh R and Zimmerman J E 1990 *Analytical model for the refrigeration power of the orifice pulse-tube refrigerator*. (National Inst. of Standards and Technology (NEL), Boulder, CO (United States). Chemical Engineering Science Div.)

[5] Kirkconnell C S and Colwell G T 1997 *Cryocoolers 9* 335–44

[6] Gedeon D 2016 *Sage User's Guide: Stirling, Pulse-Tube and Low-T Cooler Model Classes* (Athens, OH: Gedeon Associates)

[7] Zhang L, Dang H Z, Tan J, Bao D, Zhao Y B and Qian G Z 2015 *IOP Conf. Ser. Mater. Sci. Eng. 101* 012094

[8] Lopes D C 2011 *Développement d'un tube à gaz pulsé très haute fréquences* PhD thesis (Université de Grenoble)

[9] Redlich O and Kwong J N S 1949 *Chem. Rev.* **44** 233–44

[10] Dymond J H and Smith E B 1980 *The virial coefficients of pure gases and mixtures: a critical compilation* (Oxford : New York: Clarendon Press ; Oxford University Press)

[11] Bastin G and Coron J-M 2016 *Stability and Boundary Stabilization of 1-D Hyperbolic Systems* vol 88 (Cham: Springer International Publishing)

[12] Toro E F 2009 *Riemann solvers and numerical methods for fluid dynamics: a practical introduction* (Dordrecht ; New York: Springer)

[13] Zhao L, Guo L, Bai B, Hou Y and Zhang X 2003 *Int. J. Heat Mass Transf.* **46** 4779–88

[14] Krishna B S V S R 2012 *Int. J. Appl. Res. Mech. Eng.* **2** 31–6

[15] Gedeon D 1996 *NASA Contract. Rep.* 198442 74