CFD investigation on characteristics of oscillating flow and heat transfer in 3D pulse tube

International Journal of Heat and Mass Transfer 84 (2015) 401–408

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt

CFD investigation on characteristics of oscillating flow and heat transfer in 3D pulse tube Qunte Dai a,b, Yanyan Chen a, Luwei Yang a,⇑ a b

Technical Institute of Physics and Chemistry, Chinese Academy of Sciences, Beijing 100190, China Key Laboratory for Thermal Science and Power Engineering of Ministry of Education, Tsinghua University, Beijing 100084, China

a r t i c l e

i n f o

Article history: Received 7 April 2013 Received in revised form 1 March 2014 Accepted 5 January 2015

Keywords: Pulse tube refrigerator CFD Phase angle Annular effect Mass flow Heat transfer

a b s t r a c t A CFD method is used to investigate the three dimensional oscillating flow and heat transfer in the pulse tube (PT) and heat exchangers (HEs) of a pulse tube refrigerator (PTR). Some interesting phenomena and characteristics have been found in the present work, and quantitative analyses in detail have been made. The transient temperature variation resembles sinusoidal in the central part of the pulse tube, while exhibits non-sinusoidal within 10 mm away from the two end HEs, respectively. The lowest periodically averaged temperature appears not in the cold end HE but on the cross-section about 1 mm away from the cold end HE, while the highest averaged temperature appears at about 3 mm away from the hot end HE. This would increase the conduction loss, but benefits the coaxial space arrangement for a PTR. The heat transfer efficiency of the two HEs has significant influence on the PV power loss. Poor heat transfer can lead to a higher loss, even reaching to 10% of the total PV power when the heat transfer coefficient decreases to 102 W/m2 K. The efficiency of a pulse tube has been newly proposed, which can be used to describe the loss of pulse tube quantitatively. In summary, the effects of a series of factors on the loss of the pulse tube have been discussed, including the phase angle, the heat transfer coefficient of the HEs, the mass flow rate, and the gravity. The present study would give a better understanding on the mechanism of losses in a PT. Ó 2015 Elsevier Ltd. All rights reserved.

1. Introduction Pulse tube refrigerator was first invented in 1960s by two Americans, Gifford and Longsworth [1], which is called the basic type of pulse tube refrigerator. There are lots of remarkable advantages of pulse tube refrigerator (PTR) compared to other types of refrigerators, such as lower vibration, longer performance life, higher reliability and so on. A series of significant improvements have been made afterwards by introducing the orifice [2], double inlet [3], and the inertance tube [5] as phase shifting methods. Multi-stage [4] pulse tube has also been developed to achieve lower temperature. All these efforts have contributed greatly to the rapidly development of PTR. PTR has been applied broadly in recent years, especially in some space applications [6]. There are a large number of experimental studies on pulse tube refrigerator, compared to limited theoretical investigations. The mechanism of losses in a PT has not yet been completely under⇑ Corresponding author at: Technical Institute of Physics and Chemistry, Chinese Academy of Sciences, 29 Zhongguancun East Road, Beijing 100190, China. Tel./fax: +86 1062627901. E-mail address: [email protected] (L. Yang). http://dx.doi.org/10.1016/j.ijheatmasstransfer.2015.01.027 0017-9310/Ó 2015 Elsevier Ltd. All rights reserved.

stood, owing to the complicated flow process within, which involve several subjects including thermodynamics, fluid dynamics, heat transfer and acoustics. Gifford and Longsworth [7] first explained the cooling effect in the base type of PTR by a theory named surface heat pumping. Then Radebaugh [8] proposed enthalpy flow analysis theory, which indicate that phase angle between the pressure and the mass flow is the major factor for PTR’s cooling performance. Based on this, Perter and Radebaugh [9] explained the cooling mechanism in more detail by using a vector analysis method. Methods for improving the refrigerator’s performance have also been given. Mikulin [2] studied the orifice PTR in terms of the classical thermodynamics theory, and presented the expression of the cooling capacity. In 1990s, Liang et al. [10] proposed thermodynamic non-symmetry effect to illustrate the refrigeration process in a PTR. In 2004, Luo et al. [11] proposed the mesoscopic thermodynamic theory to investigate the thermodynamic cycles in the regenerator. However, all these theories are based on some ideal assumptions and are incapable of describing complex flow and heat transfer phenomenon in the PTRs. Though useful for qualitative analysis, obvious deviation exists for quantitative comparison from the experimental results.

402

Q. Dai et al. / International Journal of Heat and Mass Transfer 84 (2015) 401–408

In recent years, along with the development of actual refrigerator, the self-made programming and commercial software are used for the simulation investigation of PTR. Wang [12–16] and Ju [17] are the pioneers who developed the one dimensional numerical simulation to investigate PTR. Through a series of improvements on the numerical model, they obtained a good accordance between the simulation results and the experimental results. However, one dimensional numerical simulation, which based on lots of assumptions and simplifications, is of big difference to the real situation. CFD investigations on PTRs were conducted recently [18–21]. Savio Antao et al. [18] developed a detailed time-dependent axisymmetric computational fluid dynamic (CFD) model, where heat transfer in the porous media is achieved and the effects of wall thicknesses and natural convective losses are included. He et al. [19] developed an approach to predict the performance of a PTR, by using a one-dimensional model for the regenerator and a two-dimensional pulse tube. The complicated fluid flow and heat transfer phenomena such as DC flow, velocity and temperature annular effects could be gained. Zhang et al. [20] performed a twodimensional axisymmetric CFD simulation. The complicated phenomena of flow and heat transfer were studied and a swirling flow pattern was observed in the pulse tube. Chen et al. [21] studied the thermodynamic cycles in the PTR by means of CFD method. It shows that gas parcels working in different parts undergo different thermodynamic cycles. All these work are based on one dimensional or two dimensional models. Though these studies reveal some typical characteristics, an in-depth and quantitative analysis on the pulse tube where the temperature varied most extraordinarily is still lacking. 3D CFD studies of PTRs would give a deeper understanding of both the 3D oscillating flow and the associated loss mechanisms. We have carried out a lot of experiments on PTRs [22] and found that some problems have not been investigated completely until now, for example, the effect of tube diameter on the portion of pulse tube loss to the total refrigeration, the effect of tilting angle for different diameters, the losses at 4 K temperature scope, and the effect of non-uniform in-flow. These practical problems are of great importance to the cryocooler’s design and performance. It is essential to develop three-dimensional model to study these problems systematically. In this paper, based on the real sizes of a pulse tube refrigerator [22], a three dimensional model is constructed on the pulse tube and the two HEs. In the process of developing the 3D model, some characteristics have been found in the pulse tube which were seldom paid attention to by the former researchers. So detailed analyses on the pulse tube mechanisms have been made. Mainly the effects of various factors on the performance of the pulse tube have been studied. These factors including heat transfer coefficient of HEs, tilting angle, PV power level and the mass flow rate, etc.

All the sizes of each part are presented in the figure. The diameters of all the three parts of the model are the same, and the diameter shown in the figure represents the inner diameter. In the present work, the commercial CFD code FLUENT is used to execute the calculation, where the FDM (Finite Difference Method) scheme is adopted. The governing equations of Fluent used to solve the mass, momentum, and energy for the working medium in this computational model are [23]:

@q þ r  ðq~ mÞ ¼ 0 @t

q

ð1Þ

D~ m n  r~ pþr ¼ qg~ Dt

qc

 





2 3

l r~m þ r~mT  r  ~mI



  @T þ ð~ m  rÞT ¼ r  ðkrTÞ @t

ð2Þ

ð3Þ

where q; ~ m; ~ p; T represent the density, velocity, pressure and temperature of the local working medium. t is the time, c the heat capacity, and g the acceleration of gravity. The flow in this model is driven by an oscillating pressure wave, with a pressure inlet boundary condition: Pin = Pm + Pa sin (xt), and a mass flow outlet boundary condition: min = ma sin (xt  h). Pm = 2.5 MPa represents the periodically averaged pressure, i.e. the operating condition pressure; h indicates the leading phase angle of inlet pressure to outlet mass flow. The two HEs worked under constant temperature surroundings. The external surfaces of both the cold end HE and hot end HE are set at constant 80 and 300 K respectively. The pulse tube is adiabatic. In order to speed up the calculation, an initial gas temperature in the pulse tube has been set to vary linearly along the tube. Taking into account the turbulent flow in the pulse tube, the Spalart–Allmaras (one equation) [24] turbulent model was adopted in the simulation of the pulse tube. The laminar flow model was used in the two HEs. The working medium helium was treated as an ideal gas because the temperature was higher than 40 K. The thermo-physical properties of the working gas were set to be temperature dependent. The thermal conductivity (W/m K) was set as following [25]: 2

T < 300 K; k ¼ 0:0038T 0:65 þ 0:0172e½0:025ðT6Þ  W=m K T P 300 K;  k ¼ 4:74  105 T 2 þ 0:35044T þ 60:848  103 W=m K The viscosity (kg/m s) was set as following [25]:

T < 300 K;

l¼



29:0 T

1:27

 þ 0:52T 0:64  106 kg=m s

T P 300 K;  l ¼ 5:7  106 T 2 þ 4:4393  102 T þ 7:5665  106 kg=m s

2. Physical model A schematic for the three dimensional simulation model is shown in Fig. 1. This model only includes partial components of a linear PTR, including cold end HE, pulse tube and hot end HE.

Both of the HEs adopted porous medium model, with the porosity of 0.7. Non-equilibrium model was adopted in the heat transfer

Fig. 1. Schematic of simulation model of pulse tube with two exchangers.

403

Q. Dai et al. / International Journal of Heat and Mass Transfer 84 (2015) 401–408

2.0

3. Simulation results and analyses 3.1. Oscillating flow characteristics and PV power Fig. 2 a presents the axial velocity distribution at Z = 50 mm (middle of the pulse tube) and X = 0 mm (central of the pulse tube) at 8 moments in a period. Fig. 2 b shows the velocity profile near the walls with a higher resolution. The maximum velocity emerges not on the central position, but in the vicinity of the pulse tube wall. Shiraishi [27] et al. also found this phenomenon in a visual experiment of working gas inside a pulse tube refrigerator. The velocity distribute evenly on the most region of the cross section. In addition, the region of the maximum velocity in the neighborhood of tube wall is shaped annularly, which we call the ‘‘velocity

4/8T T

0.8 0.4 0.0 -0.4 -0.8

-1.6

ð4Þ

where l is the dynamic viscosity coefficient of the working gas, and f is the frequency of the pressure wave. We set the boundary mesh based on the magnitude of the viscous penetration depth. Taking into account that the minimum viscous penetration depth emerges in the cold end, which is 0.07 mm, the first layer of boundary mesh was set as 0.03 mm depth. In the present research, many cases with different boundary conditions were investigated. The conditions of the base case are: the phase angle of 68°, zero gravity, the mass flux rate min = 8.16 kg/m2 s, the frequency f = 40 Hz, and the heat transfer coefficient of h = 105 W/m2 K. Other cases with different conditions were simulated for comparison. All the results in the present work were based on the base case unless otherwise notified. In the simulation process, convergence criterion is very important. Convergence conditions for oscillating flow in pulse tube are different from steady flow. For steady flow, convergence condition is achieved when the relative or absolute error is less than a preset infinitesimal value. But for an oscillating flow in a PTR which works under periodically steady state conditions, one needs periodically convergence. In other words, after one period of the iterations, all the important parameters like pressure, mass flow and especially temperature, is very close to the value of the previous period, with the difference within a certain number (de/e < 105). Then the simulation can be viewed as a steady oscillation. In actual, the pressure and mass flow is easy to become steady, while the temperature is much slower to be steady. As a result of the fast convergence of pressure and velocity field, the PV power also reaches convergence quickly. In addition, one can refine time step to make the temperature field to be steady. After the simulation is judged as steady, we can obtain the information of the amplitude and the phase angle of various variables through the FFT (Fast Fourier Transformation) analysis method, including that of pressure wave, velocity wave, temperature, mass flow and so on. And finally we can obtain the PV power of distribution.

3/8T 7/8T

-1.2

-6

-5

-4

-3

-2

-1

0

1

2

3

4

5

6

Radial direction (mm)

(a) Complete figure 2.0

1/8T 5/8T

1.6

2/8T 6/8T

3/8T 7/8T

4/8T T

1.2

Axial velocity (m/s)

kl ¼

2/8T 6/8T

1.2

rffiffiffiffiffiffiffiffiffi

l qpf

1/8T 5/8T

1.6

Axial velocity (m/s)

calculation, with the heat transfer coefficient varying from 102 to 105 W/m2 K. The three dimensional mesh model was implemented by means of GAMBIT [26] preprocessing software, which was used for geometry modeling and meshing. Hexahedral cells were adopted in all the regions. Boundary layer meshes were used for the near wall boundary regions. The local grid refinement was used in the inlet and outlet region and the interfaces between the pulse tube and the two HEs. The total amount of computational cells is above 23,000. The expression of viscous penetration depth of an oscillating flow is

0.8 0.4 0.0 -0.4 -0.8 -1.2 -1.6 -6

-5

-4

3

4

5

6

Radial direction (mm)

(b) Partially enlarged figure Fig. 2. Velocity distributions at Z = 50 mm and X = 0 mm at 8 moments (in the PT, base case).

annular effect’’. In detail, the maximum velocity emerges at a position which is about 0.5 mm distant from the tube wall, and the situations on other sections are quite similar. In other words, the influencing area of the ‘‘velocity annular effect’’ is within 1 mm from the tube wall which can be treated as the viscous boundary layer. When pulse tube diameter becomes smaller, its portion to total flow area will be larger and the effect on the system performance will be greater. Shown in Fig. 3 is the velocity distribution at Z = 10 mm (middle of the cold end HE) and X = 0 mm at 8 moments. It is evidently shown that there is no velocity annular effect in the HE compared to that in the pulse tube, the velocity distribution is more uniform. Situations on other sections in the two HEs are also very similar. This is due to the porous media model used in the HEs compared to an empty tube model for the pulse tube. The practical HEs are composed by uniformly distributed small ducts to enhance the heat transfer. So in the porous media model settings, the flow resistance is manually set and with its value much higher than that in the pulse tube. For the model in calculation, the cases with different phase angle and different heat transfer coefficient were simulated. Fig. 4 shows the comparisons of cross-sectional mean velocity amplitude along the Z axis in different conditions. When the phase

Q. Dai et al. / International Journal of Heat and Mass Transfer 84 (2015) 401–408

1.6 1.4 1.2 1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 -0.6 -0.8 -1.0 -1.2 -1.4 -6

U ¼ jU j sin ðxt þ hÞ 1/8T 5/8T

2/8T 6/8T

3/8T 7/8T

4/8T T

W ¼ 0:5jPjjU j cos h

-5

-4

-3

-2

-1

0

1

2

3

4

5

6

Radial direction (mm) Fig. 3. Velocity distributions at Z = 10 mm and X = 0 mm at 8 moments (in the HE, base case).

2.6

Axial velocity amplitude (m/s)

2.4 2.2 2.0 1.8 1.6 1.4 1.2 1.0

0

68 0 0

0.8

0

0

60 30 0 2 68 & h=100W/m .K

0.6 0

10

20

30

40

50

60

70

80

90

Axial direction (mm) Fig. 4. Comparisons of sectional velocity amplitude along the Z axis in different conditions.

angle h is 0°, which represents a pure resistance impedance at the hot end of the pulse tube as the case in an orifice type of pulse tube, the velocity amplitude drops gradually from the cold end to the hot end. The velocity amplitude distributes evenly along the tube when h is 30°, while increases gradually along the tube when h is bigger. Further observation shows that the heat transfer coefficient has little effect on the velocity field. The maximum within a period of the working gas in the pulse tube can be calculated by 2|u|/x. |u| represents the axial velocity amplitude, x represents the angular frequency. For the base case with 68° of phase angle, the maximum displacement is 15.07 mm in the hot end, 25.11% of the pulse tube length, and 9.09 mm in the cold end respectively, equals to 15.14% of the pulse tube length. The working gas within and near the HEs oscillate in such a large displacement would increase the cooling losses through the pulse tube. This is why the length of the pulse tube is normally several times larger than the working gas displacement to separate the hot and the cold ends. The definition of PV power is:

W¼

1

s

Z

ð7Þ

Then

ð8Þ

where |P| indicates the pressure wave amplitude on a cross section, and |U| indicates the axial volume flow rate amplitude on the same cross section, h is the corresponding angle between pressure wave and velocity wave. Fig. 5 shows the PV diagram on Z = 21, 53 and 77 mm sections, respectively. The volume on the abscissa is the volume that the sectional fluid deviate from the equilibrium position. Each curve is an ellipse, and the area of the ellipse indicates the PV power of this section. Actually, the areas of the three ellipses are very close to each other, which validates the following results shown in Figs. 6–8. Fig. 6 shows the comparison of sectional PV power along the Z axis with different phase angles. The six curves exhibit the same variation tendency. It is reasonable that for the same mass flow boundary condition, a smaller phase angle results in a higher PV power as shown in Eq. (8). PV power reduction can be observed in the HEs for all conditions especially for higher power cases, although the mass flow rate within remains similar for all these cases. Fig. 7 presents the PV power variation along the pulse tube with different HE heat transfer coefficient varying from 102 to 105 W/m2 K. The phase angle of all the cases here is 68°, the same with the base case. For heat transfer coefficient equals 105 or 104 W/m2 K, the minimum PV power level in Fig. 7 means the HE can be treated as ideal in heat transfer for the working fluid and the heat source. Along the pulse tube, the PV power keeps nearly constant at an approximate value of 11.0 W for all heat transfer coefficients. Nevertheless, there is loss near the two interfaces of about 2%. This is because of the local inhomogeneous flow induced by the sudden change in velocity at these interfaces which lead to an increase in local viscous drag and consequently to a PV power loss. The sudden velocity change occurs because of the sudden change of the flow cross-section area as the porosity is 1 for the pulse tube while 0.7 for the HEs. In addition, with the decreasing of the coefficient, the differences of PV power drop in the two HEs increase obviously. There is an obvious PV power reduction along the two HEs, especially when the heat transfer coefficient decreases to 102 W/m2 K. For the 102 W/m2 K case, PV power of 1.2 W in the two HEs is lost, which is more than 10% of the total PV power.

Displacement (mm) 0

2

4

6

8

10

12

14

300 Z=21mm Z=53mm Z=77mm

200

Pressure amplitude P (kPa)

Axial velocity (m/s)

404

100 0 -100 -200

T

P  Udt

ð5Þ

200

400

600

800 1000 1200 1400 1600 1800 3

As

P ¼ P0 þ jPj sin ðxtÞ

-300 0

0

Volume (mm )

ð6Þ

Fig. 5. PV diagram on Z = 21, 53 and 77 mm sections (base case).

Q. Dai et al. / International Journal of Heat and Mass Transfer 84 (2015) 401–408

30 28 26

PV power (W)

24 22 o

75 o 60

20 18

o

65 o 0

50

60

68 o 30

o

16 14 12 10 8 0

10

20

30

40

70

80

90

Axial direction (mm) Fig. 6. Comparisons of PV power along the Z axis with different phase angles.

12.0

PV power (W)

11.8

2

100000 W/m .K 2 10000 W/m .K 2 1000 W/m .K 2 100 W/m .K

11.6 11.4 11.2

405

102 W/m2 K. Obviously, the efficiency of the HEs influences the efficiency of pulse tube a lot. When the efficiency is low, it requires higher input PV power into the pulse tube. When the gravity is considered, due to the different densities of the hot and cold fluid, natural convection may happen as heavy fluid sink down while lighter fluid goes up. This can cause gas mixing in the pulse tube between the hot and cold fluids, which can lead to additional losses. In order to investigate the influence of the gravity, a typical 3-D simulation result is shown in Fig. 8, in which the sectional PV power variation along the Z axis is presented for different tilting angles. The tilting angle is defined as the axial direction of the pulse tube with respect to the direction of the gravity vector. The tilting angle is defined as 0° when the pulse tube is vertical with the cold end facing downwards. The tilting angle of 90° means the pulse tube is horizontal. As shown in Fig. 8, the PV power along the tube is nearly the same for the condition of inclination angle 0° or 45°, which is also similar to zero gravity cases. The PV power loss begins to increase at a tilting angle of 90°, and attains a maximum at 135°. This phenomenon is accordance with the experimental results. In addition, the PV power curve with gravity influence fluctuates along the tube, while for zero gravity the PV power curve almost keeps static along the tube. In general, the tilting angle of 135° is the worst condition where the PV power drop is the biggest, which should be avoided in practical applications. In fact, for the present condition with PV power of 11 W in the PT, the influence of natural convection driven by gravity is slightly small. This is because the forced convection effect is much larger than the natural convection. But if the input PV power decreases, the forced convection effect becomes weaker in the PT, then the influence of the natural convection would be much more significant.

11.0

3.2. Temperature characteristics 10.8 0

10

20

30

40

50

60

70

80

90

Axial direction (mm) Fig. 7. Comparison of PV power varying along the Z axis with different HEs.

11.4 zero gravity 0 tilting angle=0 0 tilting angle=45 0 tilting angle=90 0 tilting angle=135 0 tilting angle=180

11.3

PV power (W)

11.2 11.1 11.0 10.9 10.8 0

10

20

30

40

50

60

70

80

90

Axial direction (mm) Fig. 8. Comparisons of PV power varying along the Z axis with different tilting angles.

Here we propose a definition of the efficiency of the pulse tube with HEs, which is the hot end PV power divided by the cold end PV power. The efficiency of the pulse tube is 0.96 when the HEs are ideal (for the heat transfer coefficient P104 W/m2 K), while it drops to 0.90 when the heat transfer coefficient decreases to

Shown in Fig. 9a are the comparisons of periodically averaged temperature (gas) distribution along the pulse tube in different conditions. In the region of the pulse tube near the cold end HE, the periodically averaged temperature is below 80 K, which means refrigeration effect due to gas expansion. Accordingly, it is above 300 K in the region of the pulse tube near the hot end HE due to compression. Fig. 9 shows the periodically averaged cross-sectional mean temperature below 80 K near the cold end and above 300 K near the hot end, with a higher resolution. Also in Fig. 9 a, we can find that the mass flux plays a great role to the temperature field. When the mass flux is smaller, there is a dramatic temperature change near the HEs, while the temperature changes gently along the pulse tube. In addition, the phase angle has only a slightly small influence on the temperature distribution, and its effect becomes slightly bigger only when the angle is small. For the base case, the length of space above 300 K in the pulse tube is about 7 mm, while below 80 K about 4 mm long. As a typical 3-D result for the temperature considering the gravity, Fig. 10 shows the periodically averaged cross-sectional mean temperature distribution for different tilting angles. The temperature distribution along the tube for tilting angle of 0° and 45° is also nearly the same as that for zero gravity. The temperature curve for 135° deviates the farthest from the zero gravity condition. Once again, it can be concluded that the tilting angle of 135° is the worst condition where the natural convection is the strongest. The cross-sectional mean temperature distributions on various sections in the pulse tube for base case are presented in Fig. 11 as a function of time. It can be seen that the temperature in most sections in the middle region varies like a sinusoidal wave. In addition, there is a small offset in the phase angle of the temperature curves, of about 10° from the cold end to the hot end. Further

406

Q. Dai et al. / International Journal of Heat and Mass Transfer 84 (2015) 401–408

320

350

240 200 160

0

0

75 68 0 0 65 60 0 0 30 0 o 2 75 & ma=3.6 kg/m .s

120 80 20

o

40

50

150 100

68 & ma=3.6 kg/m .s

30

200

2

60

70

50 0.000

80

0.005

310

Mean temperature (K)

300 290 280 0

0

75 68 0 0 65 60 0 0 30 0 o 2 75 & ma=3.6 kg/m .s

100 90 80

o

2

68 & ma=3.6 kg/m .s

70 22

24

26

74

76

78

80

Axial direction (mm)

(b) Partially enlarged figure

270

80

Tcold (K)

82

240 210 180

60 30

40

50

60

0.025

70

78 Tcold

76 312

zero gravity 0 tilting angle=0 0 tilting angle=45 0 tilting angle=90 0 tilting angle=135 0 tilting angle=180

90

0.020

observation shows that the temperature variation changes from sinusoidal to non-sinusoidal near the Z = 68 mm section. Correspondingly, the temperature variation shows deviations from sinusoidal behavior near the cold end at about Z = 29 mm. The transient temperature variation in one period on Z = 20.5 and 79.5 mm sections are presented in Fig. 12, respectively. The temperature of the cold end remains around 80 K for nearly half a period, shown in the right part of the figure. On the other half period, the temperature is lower than 80 K, which means cooling effect, shown in the left side. For the hot end temperature, the fluid keeps around 300 K for almost half a period, as shown between the two vertical lines. On the other half period, the temperature is higher than 300 K, which means releasing heat to the environment. But the division between in-flow and out-flow is not so clear if only surveying the temperature oscillation. This is different from the one dimension theoretical analysis or calculation in former researches: one period was simply divided into two parts, and there is a clear boundary point between them. The in-flow keeps 300 or 80 K while out-flow maintains below 80 K or higher than 300 K. Our simulation results show much more complex variation than that of the simple model. Even very close to the interface, the temperature variation not simple, especially for the small part higher than 80 K at the cold end and the small part lower than 300 K at the hot end. Fig. 13 shows the temperature variations near the cold end of the pulse tube at sections Z = 20.1, 20.3, 20.5, 21 and 22 mm. The

300

120

0.015

Time (s)

308

Thot (K)

Mean temperature (K)

Fig. 9. Periodically averaged sectional mean temperature distributions along the axis (in the PT).

150

0.010

Fig. 11. Mean temperature distributions on various sections in the PT (base case).

(a) Complete figure

20

Z=43 mm Z=68 mm

250

Axial direction (mm)

20

Z=33 mm Z=63 mm

300

Temperature (K)

Mean temperature (K)

280

Z=29 mm Z=53 mm

80

Axial direction (mm)

304

Thot

300 296 0.000

0.005

0.010

0.015

0.020

0.025

Time (s) Fig. 10. Periodically averaged sectional mean temperature distributions along the axis (in the PT).

Fig. 12. Periodical temperature distributions at Z = 20.5 and 79.5 mm (base case).

407

Q. Dai et al. / International Journal of Heat and Mass Transfer 84 (2015) 401–408

316

84 83

312

82

Temperature (K)

Temperature (K)

81 80 79 78

Z=20.1 mm Z=20.3 mm Z=20.5 mm Z=21 mm Z=22 mm

77 76 75 74 0.000

0.005

0.010

0.015

308 304 300 Z=79.5 mm Z=79 mm Z=78 mm Z=77 mm

296 292

0.020

0.025

0.000

0.005

0.010

Time (s)

0.015

0.020

0.025

Time (s)

Fig. 13. Periodical temperature distributions around the cold end (base case).

Fig. 15. Periodical temperature distributions around the hot end (base case).

temperature near the cold end fluctuates around 80 K. Obviously, in the half period of out-flow, the temperatures in different sections vary similarly and approximately sinusoidal. For comparison, in the other half period of in-flow, the temperature oscillations show a quite different behavior. Closer to the interface, the temperature is closer to 80 K. It can also be found that the maximum temperature in the small area within 2 mm varied in a range of almost 4 K. And there is a small part of area where the temperature is higher than 80 K. This phenomenon leads to a heat transfer of gas mixing, which is bad for reaching a lower refrigeration temperature. For a fixed position, the gas temperature varied periodically and the wall temperature is nearly the average value. There is heat transfer between the gas and wall, which also lead to a loss. One could calculate the heat exchange amount between the gas and the wall. This heat transfer is quite like the shuttle loss in a Stirling cooler. It is useful to gain this heat transfer amount in simulation and try to decrease it. Fig. 14 represents the maximum, minimum and average temperature at different sections around the cold end. The lowest value of periodically averaged temperature, namely 78.57 K, appears at the Z = 21 mm section, while the lowest value of temperature in a period, namely 74.47 K, appears at the Z = 25 mm section. Fig. 15 shows the periodical temperature variation near the hot end of the pulse tube at sections Z = 77, 78, 79 and 79.5 mm. The temperature close to the hot end in the pulse tube fluctuates near 300 K. In the half period of out-flow, the temperature at different

sections perform similarly and exhibit an approximately sinusoidal variation. Whereas in the half period of in-flow, the temperatures show a rather different behavior. Closer to the interface, the temperature is closer to 300 K. Fig. 16 illustrates the maximum, minimum and average temperature at different sections around the hot end. The highest value of periodically averaged temperature appears on the Z = 77 mm section, 304.92 K, while the highest value of temperature in a period appears on the Z = 68 mm section, 321.74 K. As seen from Figs. 13–16, there exists a high temperature platform and a low temperature platform at both pulse tube ends. These platforms evidently shorten the thermal baffle length between the hot temperature Th and the cold temperature Tc, which will result in a higher conduction loss and may bring hot and cold fluid mixing if the length is less than the fluid displacement. However, this effect is usually not considered in pulse tube design, from the present simulation it seems that actual conduction loss is at least 10% higher than calculated from the length of the pulse tube. Therefore, the pulse tube should be prolonged to decrease the conduction loss and reduce the possibility of mixing. However, the existence of these two temperature platforms is beneficial to the choice of pulse tube and regenerator length in a coaxial refrigerator design. Nowadays, by use of different lengths of pulse tube and regenerator tube it is easy to arrange their relative position in a coaxial configuration. The present simulation results may be used as a guide for a proper arrangement of the tubes.

100

Mean Temp. Highest Temp. Lowest Temp.

310

92

Temperature (K)

Temperature (K)

96

320

88 84 80

300 290 280 270 Mean Temp. Highest Temp. Lowest Temp.

260

76 250

72 20

21

22

23

24

25

26

27

28

29

Axial direction (mm) Fig. 14. Three kinds of temperature distributions around the cold end (base case).

66

68

70

72

74

76

78

80

Axial direction (mm) Fig. 16. Three kinds of statistical temperature around the hot end (base case).

408

Q. Dai et al. / International Journal of Heat and Mass Transfer 84 (2015) 401–408

4. Conclusions

References

A three dimensional simulation investigation is conducted on a high frequency pulse tube by use of FLUENT software. The simulation results make possible a quantitative analysis of the ‘‘velocity annular effect’’ and the temperature distribution characteristics. We have also investigated the effects of a series of factors on the performance of a pulse tube, including phase angle, heat transfer coefficient of the HEs, mass flow, and gravity. The influence area of the ‘‘velocity annular effect’’ is within about 0.5 mm distance from the tube wall and the influence of the effect on the performance is stronger when the tube diameter is smaller. We have also found that the ‘‘annular effect’’ is evidently small in porous media or at higher flow resistance. In the pulse tube, the heat transfer between gas and tube wall is strong, which will contribute to the pulse tube loss. There exist a high temperature and a low temperature platform at the two ends of the pulse tube, which will increase the conduction loss by at least 10%. On the other hand, these two temperature platforms can be helpful for the proper arrangement and choice of pulse tube and regenerator length in a coaxial refrigerator design. The phase angle plays an important role in the velocity distribution, whereas the influence on the temperature distribution is small. Poor heat transfer in the HEs will cause a large PV power loss. When the heat transfer coefficient decreases to 102, a big drop of PV power of 1.2 W is calculated, which is more than 10% of the total PV power. A definition of the efficiency of a pulse tube is proposed, which is mainly affected by the quality of the HEs and the annular flow. Gravity also influences the performance of an inclined pulse tube by natural convection. A tilting angle of 135° is the worst condition for which the PV power loss is the highest. For tilting angles up to 90°, the influence of natural convection by gravity is weak. The effect of natural convection would be more significant when the input PV power is reduced. The flow and heat transfer characteristics found in this study will be helpful for a better understanding of the losses in a pulse tube, and can also provide theoretical guidance for system improvement and component design.

[1] W. Gifford, R. Longsworth, Pulse tube refrigeration, Trans. ASME J. Eng. Ind. (Ser. B) 86 (1964) 264–268. [2] E.A. Mikulin, A.A. Tarasov, M.P. Shkrebyonock, Low-temperature expansion pulse tubes, Adv. Cryo Eng. 29 (1984) 629–637. [3] S. Zhu, P. Wu, Z. Chen, A single stage double inlet pulse tube refrigerator capable of reaching 42 K, in: ICEC 13 Proceedings of Cryogenics, vol. 30, Beijing, 1990, pp. 256–261. [4] J. Cai, J. Wang, Y. Zhou, Experimental analysis of double-inlet principle in pulse tube refrigerator, Cryogenics 33 (1993) 522–525. [5] K. Kanao, N. Watanabe, Y. Kanazawa, A miniature pulse tube refrigerator for temperature below 100 K, Cryogenics 34 (suppl. 1) (1994) 167–170. [6] R.G. Ross, R.F. Boyle, An overview of NASA space cryocooler programs-2006, in: 14th International Cryocooler Conference, pp. 1. [7] W. Gifford, R. Longsworth, Surface heat pumping, Adv. Cryo Eng. 11 (1966) 171–181. [8] R. Radebaugh, S. Herrmann, Refrigeration efficiency of pulse-tube refrigerator, in: 4th International Cryocoolers, Conference, 1986, pp. 119. [9] S. Peter, R. Radebaugh, Development and experimental test of an analysis model of the orifice pulse tube refrigerator, Adv. Cryo Eng. 33 (1987) 851–859. [10] J. Liang, A. Ravex, P. Rolland, Study on pulse tube refrigeration Part 1: thermodynamic non-symmetry effect, Cryogenics 36 (2) (1996) 87–93. [11] E. Luo, W. Dai, Z. Wu, et al., Meso-scope thermodynamic theory for cyclic flow engines, Cryogenics 1 (2004) 1–11 [in Chinese]. [12] C. Wang, P.Y. Wu, Z.Q. Chen, Numerical modeling of an orifice pulse tube refrigerator, Cryogenics 32 (1992) 785. [13] C. Wang, P.Y. Wu, Z.Q. Chen, Numerical analysis of a double-inlet pulse tube refrigerator, Cryogenics 33 (1993) 526. [14] C. Wang, P. Wu, Z. Chen, Theoretical and experimental study of a double-in let r eversible pulse tube refrigerator, Cryogenics 33 (1993) 648. [15] C. Wang, Numerical analysis of 4 K pulse tube coolers: Part I. Numerical simulation, Cryogenics 37 (4) (1997) P207–P213. 04. [16] C. Wang, Numerical analysis of 4 K pulse tube coolers: Part II. Performances and internal processes, Cryogenics 37 (4) (1997) 215–220. [17] Y.L. Ju, C. Wang, Y. Zhou, Numerical simulation and experimental verification of the oscillating flow in pulse tube refrigerator, Cryogenics 38 (2) (1988) P169–P176. 02. [18] Dion Savio Antao, Bakhtier Farouk, Experimental and numerical investigations of an orifice type cryogenic pulse tube refrigerator, Appl. Therm. Eng. 50 (1) (2013) 112–123. [19] Y.L. He, Y.B. Tao, F. Gao, A new computational model for entire pulse tube refrigerators: model description and numerical validation, Cryogenics 49 (2009) 84–93. [20] X.B. Zhang, L.M. Qiu, Z.H. Gan, et al., CFD study of a simple orifice pulse tube cooler, Cryogenics 47 (2007) 315–321. [21] L. Chen, Y. Zhang, E.C. Luo, T. Li, X.L. Wei, CFD analysis of thermodynamic cycles in a pulse tube refrigerator, Cryogenics 50 (11–12) (2010) 743–749. [22] L.W. Yang, H.L. Chen, M.G. Zhao, Development of high frequency pulse tube coolers at CAS, in: Cryogenics and Refrigeration, Proceedings, 2008, pp. 159– 162. [23] Fluent INC. Fluent 14.0 user manual. Fluent INC, 2011. [24] P.R. Spalart, S.R. Allmaras, A One-Equation Turbulence Model for Aerodynamic Flows, Recherche Aerospatiale, No. 1, 1994, pp. 5–21. [25] . [26] Gambit 2.4.6. Gambit user manual. [27] M. Shiraishi et al., Visualization of velocity profiles and displacements of working gas inside a pulse tube refrigerator, Cryocoolers 9 (1997) 355–364. Plenum Press.

Conflict of interest None declared. Acknowledgment This work is supported by the National Science Foundation of China (Grant Nos. 50890181, 51076160).