Cryogenics 65 (2015) 49–58

Contents lists available at ScienceDirect

Cryogenics
journal homepage: www.elsevier.com/locate/cryogenics

Transient simulation of a miniature Joule–Thomson (J–T) cryocooler

with and without the distributed J–T effect
R.M. Damle, M.D. Atrey ⇑
Refrigeration and Cryogenics Laboratory, Department of Mechanical Engineering, Indian Institute of Technology (IIT) Bombay, Mumbai 400076, Maharashtra, India

a r t i c l e i n f o a b s t r a c t

Article history: The aim of this work is to develop a transient program for the simulation of a miniature Joule–Thomson
Received 20 August 2014 (J–T) cryocooler to predict its cool-down characteristics. A one dimensional transient model is formulated
Received in revised form 18 October 2014 for the fluid streams and the solid elements of the recuperative heat exchanger. Variation of physical prop-
Accepted 27 October 2014
erties due to pressure and temperature is considered. In addition to the J–T expansion at the end of the
Available online 11 November 2014
finned tube, the distributed J–T effect along its length is also considered. It is observed that the distributed
J–T effect leads to additional cooling of the gas in the finned tube and that it cannot be neglected when the
Keywords:
pressure drop along the length of the finned tube is large. The mathematical model, method of resolution
J–T cryocooler
Transient
and the global transient algorithm, within a modular object-oriented framework, are detailed in this
Distributed J–T effect paper. As a part of verification and validation of the developed model, cases available in the literature
Heat exchanger are simulated and the results are compared with the corresponding numerical and experimental data.
Ó 2014 Elsevier Ltd. All rights reserved.

1. Introduction pressures (107 psia  689,475 bar) to liquefy nitrogen without a

recuperative heat exchanger [1]. Flynn [2] demonstrated that a
Miniature Joule–Thomson (J–T) cryocoolers are widely used for heat exchanger with 90% effectiveness is not sufficient to liquefy
cooling of infrared sensors, cryosurgery probes, thermal cameras, nitrogen for a cycle working between pressure limits of 1 bar and
etc. The main advantages of J–T type cryocoolers are: fast cool 100 bar. The effectiveness of the cryogenic heat exchangers needs
down time (of the order of seconds); simple design; no moving to be as high as 97%. Atrey [3] observed that a drop in heat exchan-
parts; high reliability; less maintenance; and low cost. The working ger effectiveness from 97% to 95% reduced the liquid helium yield
medium for such cryocoolers is usually Nitrogen or Argon gas by 12%.
which is easily available at a reasonable price. The main parts of Being the most critical part of a cryocooler, the recuperative
a J–T cryocooler are: storage vessel with high pressure gas; heat heat exchanger has been studied extensively by many authors.
exchanger; expansion valve; and evaporator/load as shown in Xue et al. [4] and Ng et al. [5] have carried out steady state analysis
Fig. 1. The high pressure gas from the vessel after flowing through of a miniature Hampson type heat exchanger with argon gas and
the heat exchanger expands through the expansion valve and cools compared it with their experimental data. An accurate geometrical
because of the J–T effect. The cooling effect is extracted in the model for the helical finned tube is included in the steady state
evaporator and thereafter, the gas flows through the low pressure thermodynamic model of the miniature heat exchanger by Chua
side of the heat exchanger. The low pressure gas cools the high et al. [6]. Hong et al. [7] used an effectiveness-NTU approach to
pressure gas on its way out to the atmosphere. This results in a predict the performance of the heat exchanger for pressures up
heat exchanger with a counter-flow arrangement and is called a to 500 bar with Argon and Nitrogen as working fluids. Recently,
recuperative heat exchanger as it recovers the ’cold’ from the Ardhapukar and Atrey [8] presented a steady state analysis for
return low pressure gas to cool the high pressure gas. the performance optimization of a miniature J–T cryocooler.
The working of the recuperative counter-flow heat exchanger Chou et al. [9] reported transient numerical analysis along with
directly affects the performance of the J–T cryocooler. The experimental data for a miniature J–T cryocooler. Chien et al. [10]
importance of a heat exchanger in cryogenic applications can be published a similar transient model where a self-regulating
emphasized by the fact that it would need extremely high bellows mechanism was added, to the model of Chou et al. [9], to
regulate the mass flow rate after reaching cryogenic temperatures.
⇑ Corresponding author. Tel.: +91 (22)2576 7522; fax: +91 (22)2572 6875. Hong et al. [11] carried out experiments to study the transient cool
E-mail address: [email protected] (M.D. Atrey). down characteristics of a miniature Joule–Thomson refrigerator

http://dx.doi.org/10.1016/j.cryogenics.2014.10.003
0011-2275/Ó 2014 Elsevier Ltd. All rights reserved.
50 R.M. Damle, M.D. Atrey / Cryogenics 65 (2015) 49–58

Nomenclature

A cross-sectional area, m2 lJT Joule–Thomson coefficient, K m2 =N

acf area correction factor q density, kg=m3
CV control volume / generic variable ðp; T; VÞ
4
Cp specific heat, J/kg K r Stefan–Boltzmann constant 5:67  108 ; W=m2 k
Dhel diameter of the helix, m sw wall shear stress, N=m2
dfi inner diameter of the finned tube, m
dx CV length, m Superscripts
ec kinetic energy, m2/s2 ⁄ guess value
er emissivity of shield 0 value at previous time step
f fanning friction factor
G mass velocity, kg=m2 s
Subscripts
h enthalpy, J/kg amb ambient
k thermal conductivity, W/m K c cold gas in the external annulus
l wetted perimeter, m h hot gas in the finned tube
L length of the finned tube/external annulus, m
ft finned tube
m mass, kg in inlet
m _ mass flow rate, kg/s m mandrel
Pr Prandtl number
out outlet
p pressure, N=m2 s shield
Re Reynolds number si inner surface of shield
SLPM standard litres per minute so outer surface of shield
t time, s
w finned tube wall
V velocity, m/s wi inner surface of finned tube
x distance in positive x-direction, m wo outer surface of finned tube
z vertical distance in z-direction, m () integral average over CV

Greek symbols
a heat transfer coefficient, W=m2 K
 precision for convergence

using nitrogen gas at inlet pressures up to 120 bar. Hong et al. [12] due to the distributed J–T effect can be significant when the pres-
also presented a numerical study of the operating characteristics of sure drop along the flow of high pressure gas is substantial. He also
a miniature Joule–Thomson refrigerator. commented on the aspects related to the modelling of this distrib-
Recently, Maytal [13] performed tests on four different con- uted J–T effect.
structions of Hampson type miniature J–T cryocoolers with mixed In general, transient numerical analysis has received less atten-
refrigerants. It was pointed out that the frictional pressure drop tion in the literature which is crucial for predicting the cool down
along the tube causes a distributed J–T expansion effect and any characteristics of a miniature J–T cryocooler. Although the model-
additional expansion device at the end of the tube can thus be ling of the time dependent accumulative (transient) terms in the
eliminated. It was also suggested that the temperatures change continuity, momentum and energy equations are inherent to a
transient formulation, Chou et al. [9] and Hong et al. [12] have
neglected the accumulative terms for the continuity and momen-
tum equations in their transient formulation. Also, the initial tem-
perature and pressure maps for all the components of the heat
exchanger at time t = 0, which are very important for the temporal
evolution of the physical phenomenon, are not clearly mentioned.
Moreover, previously presented numerical models, by Chou et al.
[9], Chien et al. [10], Ng et al. [5], and Hong et al. [12], have mod-
elled the J–T effect only at the outlet of the finned tube carrying
the high pressure gas. They have not considered the distributed
J–T effect along the length of the finned tube as pointed out by
Maytal [13]. Thus, changes in the enthalpy of the fluid due to tem-
perature alone (i.e., dh ¼ C p dT) have been considered along the
tube length. This assumption may not be valid for the finned tube
with a very large gradient of pressure along its length.
In this work, a one-dimensional transient model is developed
for the recuperative heat exchanger of the J–T cryocooler. Time
dependent terms are taken into account for the continuity,
momentum and energy equations. To simulate the entire cryogenic
cycle, as in a cryocooler, the J–T expansion process is also simu-
lated to calculate the inlet temperature of the return gas in the
external annulus. The so far neglected distributed J–T effect is also
Fig. 1. Parts of a J–T cryocooler. considered along the length of the finned tube where a huge
 R.M. Damle, M.D. Atrey / Cryogenics 65 (2015) 49–58 51

pressure drop (around 80–100 bar) is experienced by the fluid due figuration of the heat exchanger (i.e., finned tube geometry, heat
to wall friction. Both temperature and enthalpy formulations are exchanger length, etc.) and the working conditions drastically
worked out to explore the distributed J–T effect. Physical proper- affect the cool down characteristics. A numerical tool is therefore
ties are evaluated as a function of the local temperature and useful for the optimization of the geometrical and operating
pressure. For different working conditions, the numerically parameters of the heat exchanger.
obtained values of the outlet temperature of the fluid in the
external annulus are compared with experimental data. The
temperature profiles of the fluid streams over the heat exchanger 3. Mathematical model
length are also compared with available numerical data.
Also, an object-oriented and modular approach is employed for 3.1. Control volume arrangement
modelling the heat exchanger in this work. The heat exchanger is
modelled as a collection of different basic elements; i.e., fluid In this work, a one-dimensional transient model is developed
streams, separating finned tube, mandrel and shield. These ele- for the fluid streams and the solid elements (i.e., finned tube, man-
ments interact with each other only through boundary conditions. drel and shield) which together form the heat exchanger. For the
So, the method of resolution of each element can be different. This numerical resolution of the governing equations described in this
gives flexibility to choose models for individual elements as per section, it is first necessary to divide the different elements of
requirement. From the software point of view, the program is reus- the heat exchanger into a series of control volumes (CVs). The CV
able and easy to maintain. arrangements employed in this work for the finned tube, fluid in
the finned tube and the fluid in the external annulus are shown
in Fig. 3. For the finned tube, the nodes are placed at the centre
2. Miniature J–T cryocooler configuration of a CV while the nodes for the fluid streams are placed on the
CV faces. The CV arrangement for the mandrel and the shield is
The schematic diagram of a tube-in-tube counter-flow heat the same as that of the finned tube. As the finned tube is helically
exchanger for a miniature Joule–Thomson cryocooler is shown in wound on the mandrel, the total length of the tube is first calcu-
Fig. 2. The main parts of this heat exchanger are: a helical finned lated taking into account the pitch and diameter of the helix along
tube with high pressure gas, a mandrel over which the finned tube with the vertical height (L) of the heat exchanger. This length can
is wound, and a covering shield forming an external annular space then be divided into any number of control volumes. The length
for the returning low pressure gas. Normally, the high pressure gas of the fluid in the external annulus is also divided into a series of
enters the finned tube (hot side of the heat exchanger) at a pres- CVs which are of the form of annular rings.
sure in the range of 100–400 bar depending on the working fluid. In order to coincide the mesh pattern of the finned tube and
The temperature at the inlet is close to the ambient temperature inner hot fluid, same number of CVs are assigned to them. This is
(300 K). A very high pressure drop is experienced by the fluid in because the hot fluid and the finned tube, with same length, go
the helical finned tube along with heat transfer. The return gas in together around the mandrel in a helical way. It is therefore conve-
the external annulus (cold side of the heat exchanger) enters at nient to associate each CV of high pressure fluid with the corre-
low pressure around 1.5–2 bar and at an inlet temperature of about sponding CV of the finned tube. The outer annular space with the
80–110 K. The outer finned surface in the external annulus low pressure fluid is a straight vertical column of height L. It can
enhances the heat transfer to the return gas. The geometrical con- be divided into several CVs which can be different in number from
the CVs of the finned tube. The number of CVs for the shield and
the mandrel are the same as that for the outer annular fluid.
As the number of CVs of the finned tube and the inner fluid can
be different as compared to the outer annular fluid, every finned
tube CV is associated with its corresponding outer fluid CV. This
association of tubes CVs is based on their height with respect to
the height of outer fluid CVs. For example, as shown in the Fig. 4,
the tube CVs with height less than z1 are assigned to the first CV
of the outer fluid. Similarly, other tube CVs are also associated with
the outer fluid CVs. The surface area and perimeter of finned tube,
and flow area and hydraulic diameter of the external annulus are
calculated as in Ardhapurkar and Atrey [8].

3.2. Governing equations

3.2.1. Assumptions
The assumptions made in the derivation of the governing equa-
tions are:

(i) heat transfer and fluid flow is one dimensional along the
length of solid and fluid elements of the heat exchanger;
(ii) axial conduction in the fluid is neglected;
(iii) body forces and axial stresses are negligible;
(iv) the helical tube is assumed to be perfectly circular and
closely spaced;
(v) fin efficiency is assumed to be 100%;
(vi) diametrical clearance between fins and shield is neglected;
Fig. 2. Schematic of a miniature J–T cryocooler with the recuperative counter-flow
(vii) emissivity of the shield is assumed to be constant and
heat exchanger. receives outside radiation at ambient temperature.
52 R.M. Damle, M.D. Atrey / Cryogenics 65 (2015) 49–58

(a)

(b)

Fig. 3. Control volume arrangement: (a) finned tube and inner fluid and (b) outer fluid.

Multiplying the mass conservation equation by h,  subtracting it

from Eq. (3) and putting dh ¼ C p dT  lJT C p dp according to Maytal
[13] the energy equation in temperature form can be rearranged as:

@T @

 0 C 0p
Aq þ _ p ðT  TÞ ¼ a  lp  ðT w  TÞ
mC
@t @x  
@p @ðq ec Þ @ðme
_ cÞ
þA  A þ
@t @t @x
@p @  
þ Aq  0 C 0p l
 0JT þ _ p lJT ðp  p
mC Þ ð4Þ
@t @x
Here, q 0 is calculated at mean temperature T 0 and mean pressure p0
Fig. 4. Association of the finned tube CVs with annular fluid CVs.
while C 0p and l  0JT are evaluated at the mean temperature of
ðT þ T 0 Þ=2 and mean pressure of ðp þp0 Þ=2. The last two terms in
the energy equation with the temperature formulation (Eq. (4))
(viii) The initial pressure of the fluid streams, at time t ¼ 0, is represent the distributed J–T effect or the changes in temperature
assumed to be equal to the inlet pressure of the low pressure due to variation of pressure. These terms are considered only for
side (the pressure to which the high pressure gas expands). the fluid in the finned tube due the large pressure drop (around
80–100 bar) over its length.
3.2.2. Mass, momentum and energy conservation The energy equations for the solid elements are the following:
The basic equations of conservation of mass, momentum and Finned tube:
energy for the fluid elements and energy conservation equations  
@T w @ @T w
for solid elements are written in a differential form. The mean qw Aw C pw ¼ kw Aw þ ah lwi ðT h  T w Þ  ac lwo ðT w  T c Þ
properties of the fluid at the centre of the CV are calculated from @t @x @x
the values at the nodes and are indicated with a overbar. The vari- ð5Þ
ables at the previous time step are indicated with a superindex ‘0’.
Mandrel:
The conservation of mass over a fluid CV is:
 
@q
 @m _ @T m @ @T m
A þ ¼0 ð1Þ qm Am C pm ¼ km Am  ac lm ðT m  T c Þ ð6Þ
@t @x @t @x @x
The conservation of momentum is given by: Shield:
@ðq _
 VÞ @ðmVÞ @p    
A þ ¼   A  sw lp ð2Þ @T s @ @T s
@t @x @x qs As C ps ¼ ks As  ac lsi ðT s  T c Þ  rer lso T 4s  T 4amb
@t @x @x
where sw ¼ f qV 2 =2 is the wall shear stress, f is the Fanning friction ð7Þ
factor, A is the cross-sectional area and lp is the wetted perimeter.
Energy equation in terms of enthalpy is written as:
 3.3. Boundary conditions
@ðq
 hÞ _
@ðmhÞ @p
A þ ¼ a  lp  ðT w  TÞ þ A
@t @x @t The inlet temperature, pressure and mass flow rate are known
 
@ðq ec Þ @ðme
_ cÞ
for the gas in the finned tube. The same are known at the inlet of
 A þ ð3Þ
@t @x the return gas in the external annulus. All the solid elements (i.e.,
 R.M. Damle, M.D. Atrey / Cryogenics 65 (2015) 49–58 53

finned tube, mandrel and shield) are assumed to be adiabatic at

ends. Thus,
dT w dT m dT s
at x ¼ 0 & t > 0; T ¼ T h;in ; p ¼ ph;in ; ¼ 0; ¼ 0; ¼ 0 ð8Þ
dx dx dx
dT w dT m dT s
at x ¼ L & t > 0; T ¼ T a;e ; p ¼ pc;in ; ¼ 0; ¼ 0; ¼ 0 ð9Þ
dx dx dx
T a;e is the temperature of the gas after isenthalpic expansion, from
its state at the exit of the finned tube, to the pressure pc;in in the
external annulus. This temperature goes on reducing from its initial
value (ambient temperature) to cryogenic temperature T c;in at
steady state. In addition, the initial temperature map for all the
solid and fluid elements and pressure map for the fluid elements
are specified as:
Fig. 5. Base and derived classes for the fluid elements.
at t ¼ 0; T 0 ¼ T amb ; p0 ¼ pc;in for 06x6L ð10Þ

functions and variables of the base class can be reused by the

3.4. Heat transfer and friction factor correlations
derived class and at the same time the internal method of resolu-
tion for both can be different. Now, from the derived class Fluid2,
The Fanning friction factor (f) for the flow through the helical
an object f 2 is created for the return fluid in the external annulus.
finned tube, with the correlation of Timmerhaus and Flynn [14],
On the same lines, a finned tube object t1 is created from the base
is calculated as:
class Tube. Objects s1 and m1 are created from a base class
 
dfi Annular  tube representing the shield and the mandrel. These
f ¼ 0:184 1 þ 3:5 Re0:2 ð11Þ objects (elements) are then interlinked to form a thermal system
Dhel
(heat exchanger in this case) as shown in Fig. 6a.
The convective heat transfer coefficient (ah ) for the turbulent flow All these objects solve themselves for given boundary condi-
in the finned tube, calculated according to Timmerhaus and Flynn tions. These boundary conditions come from the neighbouring
[14], is given by: objects to which they are connected. For example, the f 1 and f 2
  objects get the wall temperature from the object t1, while the t1
dfi
ah ¼ 0:023C ph Gh Re0:2 Pr2=3 1 þ 3:5 ð12Þ object receives fluid temperatures and heat transfer coefficients
Dhel
from the fluid objects f 1 and f 2. The advantage is that the method
The cold gas, on its way back to the atmosphere, exchanges heat by
convection with the finned surface, outer surface of the mandrel
and the inner surface of the shield. Reduction of flow area of the
external annular region due to the finned tube is also considered (a)
according to Gupta et al. [15]. The convective heat transfer coeffi-
cient (ac ) for the return flow, with the correlation from Timmerhaus
and Flynn [14], is estimated as:

ac ¼ 0:26C pc Gc Re0:4 Pr2=3 ð13Þ

The Fanning friction factor (f) for the flow through the external
annulus is calculated according to:

f ¼ 16=Re Re < 2300 ð14Þ

0:25
f ¼ 0:079Re Re > 2300 ð15Þ

4. Numerical resolution

4.1. Modular object-oriented methodology (b)

In this work, a modular object-oriented approach is employed
for the numerical simulation of the counter-flow heat exchanger
which forms the most of the cryocooler assembly. From the pro-
gramming point of view in C++ basic classes for fluid and solid
parts are identified and defined. Instantiations, called objects, can
be derived from these basic classes to create individual elements
of the heat exchanger, namely, the high pressure fluid, the low
pressure fluid, the finned tube, the mandrel and the shield. For
example, from a base ‘Class’ called Fluid an object f 1 can be created
for the high pressure fluid in the finned tube as shown in Fig. 5. For
the low pressure fluid in the external annular space a new ‘Class’
with name Fluid2 is derived from the base class Fluid. The derived
class inherits all the properties of the base class and in addition Fig. 6. Modular methodology: (a) modular heat exchanger and (b) different levels
new features can be implemented. The advantage is that all the of modelling different elements in the same simulation.
54 R.M. Damle, M.D. Atrey / Cryogenics 65 (2015) 49–58

of resolution of each element can be different (f1 can be solved 4.3. Global algorithm
with SIMPLEC algorithm [16] while f2 can resolved with a step-
by-step method). Only the necessary boundary information The global resolution algorithm for the transient simulation of
(changing during the transient evolution) are fed to the concerned the cryocooler, modelled as a collection of different elements, is
object. Thus, each element of the system can be modelled with dif- shown in Fig. 7. Geometrical parameters, initial map (at time
ferent level of detail (see Fig. 6b) in the same simulation which t = 0) and boundary conditions are first set for each element of
gives flexibility to choose individual element models according to the heat exchanger. An iteration at a given time step consists of
the desired accuracy. Moreover, new elements (classes in C++) going to each element once and resolving the corresponding gov-
and features can be added to already existing elements without erning equations to get new value of variables (e.g., temperature,
changing basic framework. The program is thus reusable, easy to pressure and velocity). During an iteration, latest variable values
debug and also, easy to maintain. from the linked elements are available to an element whose calcu-
lation is under process. For example, if the finned tube calculation
4.2. Resolution of fluid and solid elements is done after resolving both the fluid streams then the latest values
of the fluid temperatures will be available for the heat transfer cal-
Step-by-step method is employed for the resolution of the fluid culation of the finned tube wall.
streams in the finned tube and the external annulus. In the step-by After each iteration, at each CV, absolute or relative differences
_ at a
step method, starting with the values of variables (e.g. p, T, m) are evaluated between the new variable values and the variable
given cross-section i the respective values at the next cross-section values of the previous iteration. Each element calculates these dif-
i þ 1 are calculated. This method is suitable here because the pres- ferences at each of its CV and a maximum value from these calcu-
sure, temperature and mass flow rate at inlet cross-sections of both lated values is returned as the difference value for that element. A
the finned tube and the external annulus are known and marching global maximum value (it ) is then found out from these individual
in the flow direction at each time step is possible to obtain the var- element maximum values. This value, it , represents the maximum
iable values at subsequent locations. The fluid CVs receive wall absolute or relative difference for the entire heat exchanger for
temperature from the surrounding CVs of the solid wall. that iteration. If it > ts (the convergence criteria set for every time
For the solid elements, integration of Eqs. (5)–(7) over a CV step) iterations at the same time step are continued with updated
results in a system of linear algebraic equations. TDMA (Tri-Diago- variable values / ¼ /. Time step calculation is over when the con-
nal Matrix Algorithm) method is used for solving this system of vergence criteria (i.e., it < ts ) is achieved. Thereafter, the next
equations. Heat transfer coefficients and fluid temperatures are time step calculation begins by updating the temporal map
received by the solid wall CVs from the fluid CVs in contact with /0 ¼ / of all variables. The process continues until the set steady
them. state criteria, steady , is reached.

Fig. 7. Global resolution algorithm.

 R.M. Damle, M.D. Atrey / Cryogenics 65 (2015) 49–58 55

5. Results and discussion Table 2

Operating parameters of the test cases ([5]).

In this work, cases published in the literature by Xue et al. [4] Case Flow rate (SLPM) ph,in (bar) Th,in (K) pc,in (bar) Tc,in (K)
and Ng et al. [5] have been simulated and compared with the A 10.145 140.47 291.94 1.3426 108.70
numerical model developed. The dimensions of a recuperative heat B 11.943 160.10 291.25 1.6362 109.90
exchanger [5] and its components (i.e., finned tube, mandrel, etc.) C 13.927 179.12 291.49 1.7272 110.36
are listed in Table 1. The operating parameters for these cases are
listed in Table 2. Argon is used as the working fluid in all the cases.
The properties of Argon as a function of temperature and pressure 1
are obtained from the commercial software AspenONE [17]. At any T

Maximum absolute differences (φ-φr)

p
given instant, for each CV, the physical properties like density, 0.01
V
viscosity, specific heat, thermal conductivity and enthalpy are
0.0001
calculated as a function of temperature and pressure. The values
of thermal conductivity, as a function of temperature, for the solid
1e-06
elements like finned tube, mandrel and shield are evaluated from
the functions given by Ng et al. [5]. 1e-08

5.1. Convergence criteria and mesh requirement 1e-10

1e-12
The convergence at each time step is declared if the absolute
differences of pressure, temperature and velocities are below
1e-14
1  104. Steady state is declared when the absolute differences 0 10 20 30 40 50
for temperature and velocity, and relative differences of pressure Time, (s)
are below to 1  103. It is observed that at least 500 CVs are
Fig. 8. Maximum differences at each time step for temperature, pressure and
required on the high pressure fluid in the finned tube to meet
velocity (mesh: 500/100).
the convergence criteria at each time step. Fig. 8 shows the evolu-
tion of the maximum absolute differences of temperature, pressure
and velocity during a transient simulation. It can be seen that all Table 3
these values are below 1  104 as per the set criteria. By keeping Outlet parameters for case A with different meshes.
500 CVs on the high pressure side, the number of CVs on the low Sr. no. Fluid CVs inner/ Th,out (K) ph,out Tc,out (K) pc,out
pressure side are varied to see the effect of mesh refinement on dif- outer (bar) (bar)
ferent output parameters. Table 3 demonstrates that the tempera- 1 500/100 180.390 68.3582 285.785 1.16836
ture and pressure values at the exit of finned tube and external 2 500/250 180.387 68.3549 285.789 1.16821
annulus do not vary with the number of CVs on the outer side. 3 500/500 180.388 68.3549 285.789 1.16822
An additional case with 600 CVs for both inner and outer fluids 4 600/600 180.384 68.3618 285.791 1.16824

shows that the variation with further mesh refinement is not

significant.
the distributed J–T effect is significant, a tube without wall heat
5.2. The distributed J–T effect transfer is simulated. The case is simulated with the temperature
formulation (Eq. (4)) with and without the terms for distributed
In the literature, researchers [9,10,5,12] have modelled the J–T effect. Also, the same case is simulated with the enthalpy for-
changes in enthalpy of the high pressure fluid in the finned tube mulation (Eq. (3)) as it is expected to take care of the distributed
as a function of temperature only (i.e., dh ¼ C p dT is assumed). J–T effect intrinsically. This is done to compare the results using
Maytal [13] commented that the distributed J–T effect should be both formulations and to cross check the way of implementation
considered according to: of the extra terms, representing the distributed J–T effect, in the
temperature formulation. The inlet condition of pressure, tempera-
dh ¼ C p dT  lJT C p dp ð16Þ ture and mass flow rate for this case is the same as that for case A in
Table 2. The tube diameter is taken equal to the inner diameter of
Eq. (16) attributes the changes in enthalpy of the fluid not only to
the finned tube as given in Table 1. The steady state temperature
the temperature changes but also to the changes of pressure
profiles of the gas along the tube length are shown in Fig. 9. The fric-
through the Joule–Thomson coefficient lJT . To explore whether
tional pressure drop over the length of the finned tube in this case is
of the order of 100 bar. Without the term for distributed J–T effect,
Table 1 i.e., by assuming dh ¼ C p dT, temperature profile of the gas in the
Dimensions of the recuperative heat exchanger ([5]). tube is a horizontal line of constant temperature. There is no change
Geometrical parameters Size in temperature of the gas as there is no heat transfer. When the
additional terms for distributed J–T effect (lJT C p dp) are considered
Inside diameter of finned tube, dfi 0.3 mm
Outside diameter of finned tube, dfo 0.5 mm in the temperature formulation, there is a substantial drop in tem-
Inside diameter of mandrel, dmi 2.3 mm perature of the gas (around 30 K) even without wall heat transfer.
Outside diameter of mandrel, dmo 2.5 mm The enthalpy formulation with no such extra term to account for
Inside diameter of shield, dsi 4.5 mm distributed J–T effect also shows exactly the same decrease in tem-
Outside diameter of shield, dso 4.8 mm
Length of recuperative heat exchanger, L 50 mm
perature. This is because the distributed cooling effect has been
Fin height, hf 0.25 mm inherent to the enthalpy formulation as expected. It is observed that
Fin thickness, tf 0.1 mm the enthalpy remains constant during this process and therefore the
Fin density, fd 3.3 fins/mm term distributed J–T effect is appropriate for the cooling effect
Helix diameter, Dhel 3.5 mm
produced over the length of the tube. The temperature reduction
Helix pitch, Phel 1.0 mm
for argon is large and demonstrates the distributed J–T effect in a
56 R.M. Damle, M.D. Atrey / Cryogenics 65 (2015) 49–58

295
(a) 300

280 time
290
T
260
h(T,p)=Constant
285 T (μJT,dist.) 240
Temperature (K)

Th-0 s

Temperature, (K)
flow direction

220 Th-2 s
280
ph,in=140 bar Th-8 s
flow direction 200
. Th-16 s
275
m=0.282 g/s
180 Th-20 s
Th-24 s
270 160 Th-28 s ph,in=140 bar
Th-32 s
140 pc,in=1.34 bar
265 Th-36 s .
Th-steady m=0.282 g/s
120
260 100
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.2 0.4 0.6 0.8 1
Distance (x/L) Distance (x/L)

Fig. 9. Gas temperature profiles along the adiabatic tube with temperature
formulation (T), enthalpy formulation (h(T, p) = Constant), and temperature formu- (b) 300
  
280
lation with distributed J–T effect T lJT;dist . time
260

240

Temperature, (K)
Tc-0 s
better way. This behaviour can be explained from the fact that for Tc-2 s flow direction
220
an ideal gas dh ¼ C p dT as lJT ¼ 0. Thus, there is no cooling of the Tc-8 s
200
gas due to changes in pressure. On the other hand, for a real gas, Tc-16 s
below the inversion temperature, lJT > 0 and cooling takes place 180 Tc-20 s
Tc-24 s
due to the drop in pressure. Therefore, if pressure changes are large 160 Tc-28 s pc,in=1.34 bar
over the length of a tube, the cooling effect due to the distributed J– Tc-32 s
140 ph,in=140 bar
T effect cannot be neglected as suggested by Maytal [13]. The Tc-36 s .
Tc-steady m=0.282 g/s
120
enthalpy formulation takes care of this effect intrinsically. In this
work, the effect of this distributed cooling is studied for all the cases 100
0 0.2 0.4 0.6 0.8 1
mentioned in Table 2.
Distance (x/L)

Fig. 10. Transient evolution of temperatures for case A: (a) gas in the finned tube
5.3. Comparison of temperature profiles and (b) gas in the annulus.

Cases presented in Table 2 are simulated to test the numerical

model developed in this work. The geometrical parameters of the fluid do not change with time, i.e., when steady state is reached.
heat exchanger are listed in Table 1. For these cases, Ng et al. [5] It can also be seen, that after 36 s, the temperature profiles are
have published the numerically obtained temperature profiles for practically unchanged and after 37.64 s the program terminates
cases A and C under steady state conditions. These cases are simu- as the steady state criteria is achieved. Fig. 10 also shows that
lated with the numerical model developed in this work. The simu- the criteria for declaring steady state is satisfactory.
lations are carried out with and without the distributed J–T effect The steady state temperature distributions for the case A,
terms in the energy Eq. (4) with temperature formulation. The obtained as a consequence of transient evolution, are shown in
effect of applying an area correction factor to the outer surface area Fig. 11. The temperature profiles of the hot and the cold fluid
of the finned tube, as presented by Ardhapurkar and Atrey [8], is streams without the distributed J–T effect pass well over the
also studied. numerical profiles predicted by Ng et al. [5]. As the predicted cold
As the program developed in this work is of transient nature, side temperatures are on the higher side, the hot side temperatures
the steady state is obtained as consequence of transient evolution are also higher. Chua et al. [6] have mentioned that area correction
from initial conditions for these cases. The initial temperature map factors were used by Xue et al. [4] and Ng et al. [5] to account for
for all the elements is equal to the inlet temperature of the high the reduction in the outer heat transfer area of the finned tube as it
pressure gas in the finned tube and the initial pressure is set as is spirally wound over the mandrel.
inlet pressure of the gas in the external annulus (i.e., pc;in ). At time Ardhapurkar and Atrey [8], in their steady state analysis, also
t > 0, mass flow rate is imposed at the inlet cross sections of the predicted higher temperatures on both hot and cold sides. Their
finned tube and external annulus which are then resolved in the profiles matched well with the profiles of Ng et al. [5] for an area
direction of flow with the step-by-step method considering the correction factor (acf) of 0.7 for the outer surface area of the finned
heat transfer interactions aforementioned. tube. They concluded that realistic area of the outer finned surface
The transient evolution of the temperature profiles of the hot is necessary for correct numerical prediction. In this work also, the
and the cold fluid streams for case A are shown in Fig. 10. Initially, temperature profiles agree well for an acf of 0.7 without consider-
all the temperatures along the length of the inner and outer tubes ing the distributed J–T effect. This in a way validates the present
are at the inlet temperature of the high pressure gas, i.e., at 291 K. numerical analysis.
It can be seen from Fig. 10 that the high pressure gas starts cooling When the distributed J–T effect is taken into account, even
due to the J–T effect and therefore the cold side inlet temperature without any correction factor to the outer surface of the finned
starts decreasing. This, in turn, cools the gas in the finned tube. tube, the hot side temperatures drop considerably and the temper-
Subsequent expansions and heat exchange result in further lower- ature curve falls below the temperature profiles predicted by Ng
ing of temperatures of both the high pressure gas in the finned et al. [5] and that with acf = 1.0. Lower temperatures on the hot
tube and the low pressure gas in the external annulus. This process side due to the distributed J–T effect also reduce the temperatures
continues until the temperature profiles of the hot and the cold on the cold side due to lesser heat dissipation to the cold side
 R.M. Damle, M.D. Atrey / Cryogenics 65 (2015) 49–58 57

(a) 300
(a) 300
280 280
260 260
flow direction
Temperature, (K)

240 flow direction

240

Temperature, (K)
220 220
200 200
Ng et al. (2002) Ng et al. (2002)
180 acf=1.0 180 acf=1.0
acf=0.7 acf=0.7
160 acf=1.0, μJT ph,in=140 bar 160 acf=1.0 with μJT ph,in=179 bar
140 pc,in=1.34 bar 140 pc,in=1.72 bar
. .
120
m=0.282 g/s m=0.386 g/s
120
100 100
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
Distance x/L Distance x/L

(b) 300 (b) 300

280 280
260 260
240 240
Temperature, (K)

Temperature, (K)
220 flow direction flow direction
220
200 200
Ng et al. (2002) Ng et al. (2002)
180 acf=1.0 180 acf=1.0
acf=0.7 acf=0.7
160 acf=1.0, μJT pc,in=1.34 bar 160 acf=1.0 with μJT pc,in=1.72 bar
140 ph,in=140 bar 140 ph,in=179 bar
. .
m=0.282 g/s m=0.386 g/s
120 120
100 100
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
Distance x/L Distance x/L

Fig. 11. Temperature distribution for case A: (a) gas in the finned tube and (b) gas Fig. 12. Temperature distribution for case C: (a) gas in the finned tube and (b) gas in
in the annulus. the annulus.

values are less than 2.6% for all the cases. The prediction of the out-
which further decrease the hot side temperatures. The tempera- let temperature of the return gas is in good agreement with the
ture at the end of the finned tube is also lower with the distributed experimental data with acf = 0.7 and distributed J–T effect with
J–T effect. Similar observations are made for case C as shown in no area reduction i.e., acf = 1.0. The relative differences are more
Fig. 12. for the cases without the distributed effect and acf = 1.0. Also, sim-
Area reduction with a correction factor leads to lesser heat dis- ilar to the experimental observation, the outlet temperature is les-
sipation to the cold side and in turn lower temperatures on both ser for the highest inlet pressure with largest mass flow rate. Thus
hot and cold sides. Although the effect of employing an acf is sim- the model developed in this work is able to capture the trends and
ilar to considering the distributed J–T effect, the physical heat physics of the heat transfer and fluid flow phenomenon.
transfer interactions are different. It is quite possible that both The cool down time for the different cases are shown in Table 5.
these effects could be present in the actual working of the J–T cryo- It can be seen that the cool down time reduces with increase in
cooler with different respective weightage. Also, the steady state inlet pressure as observed by Chou et al. [9]. Moreover, this
temperature profiles given by Ng et al. [5] are their numerical pre- decrease in the cool down time is observed for each column i.e.,
dictions without the distributed J–T effect and with area correction for acf = 1.0 and acf = 0.7 without the distributed J–T effect and
factor. Comparison of the temperature profiles with experimental acf = 1.0 with distributed J–T effect. This shows the consistent
values along the length of the heat exchanger is necessary to behaviour of the model developed. When the distributed effect is
understand the heat exchange process and whether any area cor- not considered, the cool down time is higher for acf = 0.7 than that
rection factor is indeed required with the distributed J–T effect. It for acf = 1.0. This is because with acf = 1.0 more heat is dissipated
is though clear that with and without the distributed J–T effect, to the cold fluid as compared to acf = 0.7 due to larger area and
the temperature profiles are different and that the contribution due to lower temperature at the exit of the finned tube. Finally,
of pressure drop in lowering the temperatures along the high pres- with distributed J–T effect, the cool down time reduces further in
sure side cannot be ignored. all the three cases. This is due to the lower temperatures along
the finned tube length and lower temperatures on the cold side
5.4. Model validation with all the outer finned surface area available for heat transfer.
The literature lacks to give more information about the geometrical
Ng et al. [5] have also published the experimental values of the configuration of the heat exchangers for miniature J–T cryocoolers.
outlet temperature (T c;out ) of the gas in the external annulus for all Cool down times are also not specified explicitly for a given config-
the cases in Table 2. In this work the numerically obtained values uration and operating conditions. Ng et al. [5] have mentioned in
of outlet temperature of the return gas in the external annulus are their work that steady state conditions were achieved in less than
compared with the experimental data for all the three cases afore- a minute for all the cases. Thus, the cool down time predicted in
mentioned. The comparison is shown in Table 4. The percent rela- this work are realistic although individual case comparisons are
tive differences (% r.d.) between the numerical and experimental necessary to fully validate the model.
58 R.M. Damle, M.D. Atrey / Cryogenics 65 (2015) 49–58

Table 4
Numerical and experimental values of Tc,out (working fluid: Argon).

Case ph,in (bar) Flow rate (SLPM) Numerical values of Tc,out Tc,out Experimental (K)
acf = 1.0 (K) acf = 0.7 (K) acf = 1.0, lJT (K)
A % r.d. 140.47 10.14 289.865 (1.71) 285.791 (0.28) 286.063 (0.38) 284.98
B % r.d. 160.10 11.94 290.389 (1.97) 286.173 (0.49) 286.456 (0.59) 284.77
C % r.d. 179.12 13.93 289.734 (2.53) 285.461 (1.02) 285.825 (1.15) 282.57

Table 5 time is also lower with the distributed J–T effect. The numerical
Cool down time with and without distributed J–T effect (working fluid: argon). values of the outlet temperature of the gas in the external annulus
Case ph,in (bar) Cool down time (numerical) are in good agreement with the experimental values of Ng et al. [5]
acf = 1.0 (s) acf = 0.7 (s) acf = 1.0, lJT (s) and the cool down times are also realistic. The model developed in
this work is able to capture the heat transfer and fluid flow charac-
A 140.47 37.64 49.96 23.78
B 160.10 31.60 41.76 19.62
teristics with physically realistic and consistent results.
C 179.12 23.54 34.30 16.11
Acknowledgements

6. Conclusions The authors thank Mr. P.M. Ardharpukar, of S. S. G. M. College of

Engineering, Shegaon – 444 203, India, for the techincal discus-
A transient model for the simulation of miniature J–T cryocool- sions about the work presented in this paper.
ers has been developed with a modular and object oriented
approach for modelling the counter-flow recuperative heat References
exchanger. The heat exchanger is modelled as collection of basic
elements linked with each other. Step-by-step algorithm is imple- [1] Barron RF. Cryogenics systems. New York: McGraw-Hill book company; 1966.
[2] Flynn TM. Cryogenic engineering. 2nd rev. ed. New York: Marcel Dekker; 2005.
mented for resolving mass, momentum and energy equations of [3] Atrey MD. Thermodynamic analysis of Collins helium liquefaction cycle.
the fluid streams. Realistic initial boundary conditions have been Cryogenics 1998;38:1199–206.
used for simulating the transient evolution. Axial heat conduction [4] Xue H, Ng KC, Wang JB. Performance evaluation of the recuperative heat
exchanger in on a miniature Joule–Thomson cooler. Appl Therm Eng
and convection at surfaces is considered for the solid elements like 2001;21:1829–44.
finned tube, mandrel and covering shield. Radiative heat transfer [5] Ng KC, Xue H, Wang JB. Experimental and numerical study on a miniature
from the external surface of the shield is taken into account. Isen- Joule–Thomson cooler for steady-state characteristics. Int J Heat Mass Transfer
2002;45:609–18.
thalpic expansion process is also simulated to complete the cryo- [6] Chua HT, Wang X, Teo HY. A numerical study of the Hampson-type miniature
genic cycle. Physical properties are evaluated as a function of Joule–Thomson cryocooler. Int J Heat Mass Transfer 2006;49:582–93.
local temperature and pressure. [7] Hong YJ, Park SJ, Choi YD. A numerical study of the performance of a
heat exchanger for a miniature Joule–Thomson refrigerator. International
Traditionally, the numerical models on miniature J–T heat
Cryocooler Conference Cryocoolers, vol. 15; 2009. p. 379–86.
exchangers have attributed the changes in enthalpy of the fluid [8] Ardhapurkar PM, Atrey M. Performance optimization of a miniature Joule–
in the finned tube to changes in temperature only. The enthalpy Thomson cryocooler using numerical model. Cryogenics 2014;63:94–101.
changes due to variation of pressure, the distributed effect as sug- [9] Chou FC, Pai CF, Chien SB, Chen JS. Preliminary experimental and numerical
study of transient characteristics for Joule–Thomson cryocooler. Cryogenics
gested by Maytal [13], have been neglected. By carrying out simu- 1995;35:311–6.
lations with enthalpy formulation and temperature formulation [10] Chien SB, Chen JS, Chou FC. A study on the transient characteristics of a self-
(with and without distributed J–T effect), it is shown that the regulating Joule–Thomson cryocooler. Cryogenics 1996;36:979–84.
[11] Hong YJ, Park SJ, Kim HB, Choi YD. The cool-down characteristics of a
changes of enthalpy, and therefore of temperature, due to variation miniature Joule-Thomson refrigerator. Cryogenics 2006;46:391–5.
in pressure along the length of the finned tube are significant and [12] Hong YJ, Park SJ, Choi YD. A numerical study on operating characteristics of a
cannot be neglected when the pressure variations are large. miniature Joule–Thomson refigerator. Supercond Cryogenics 2010;12:41–5.
[13] Maytal B-Z. Hampson’s type cryocoolers with distributed Joule–Thomson
Steady state profiles, obtained as consequence of transient evo- effect for mixed refrigerants closed cycle. Cryogenics 2014;61:92–6.
lution, with area correction factor and without the distributed J–T [14] Timmerhaus KD, Flynn TM. Cryogenic process engineering. New York,
effect match well with the numerical values of Ng et al. [5]. Higher USA: Plenum press; 1989.
[15] Gupta PK, Kush PK, Tiwari A. Design and optimization of coil finned-tube heat
temperature values on both hot and cold sides are predicted with- exchangers for cryogenic applications. Cryogenics 2007;47:322–32.
out the distributed J–T effect and no area correction. Considering [16] Patankar SV. Numerical heat transfer. Hemisphere Publishing Corporation;
the distributed J–T effect, even without area correction, lower tem- 1980.
[17] AspenONE V 7.1. Aspen Technology Inc., Burlington, MA 01803, USA; 2009.
peratures of the hot and the cold sides are observed. The cool down