LBM numerical study on oscillating flow and heat transfer in porous media

Applied Thermal Engineering 54 (2013) 16e25

Contents lists available at SciVerse ScienceDirect

Applied Thermal Engineering journal homepage: www.elsevier.com/locate/apthermeng

LBM numerical study on oscillating flow and heat transfer in porous media Qunte Dai a, b, Luwei Yang a, * a b

Technical Institute of Physics and Chemistry, CAS, 29 Zhongguancun East Road, Haidia, Beijing 100190, China Graduate University of Chinese Academy of Sciences, Beijing 100190, China

h i g h l i g h t s < We initially adopted a mesoscopic method to the calculation of cryocoolers. < A thermal LBM was applied to calculate the velocity and temperature coupled field. < We constructed porous media by programming according to the practical mesh screens. < We obtained detailed information in porous media, though impossible by classical CFD.

a r t i c l e i n f o

a b s t r a c t

Article history: Received 26 September 2012 Accepted 19 January 2013 Available online 28 January 2013

Oscillating flow and heat transfer in regenerative cryocoolers, are very important for the optimization of cryocooler’s performance. A numerical study was conducted in such a system by Lattice Boltzmann Method (LBM), which is an efficiently new way compared to that of the traditional continuum Navier eStokes method. The simulation work was firstly developed in a two-dimensional empty planar channel, and then followed by porous media. The oscillating flow is driven by a periodic pressure wave, with the isothermal or adiabatic planar wall. A coupled double distribution function model, which is one of the typical thermal LBMs, is used to investigate the thermo-hydrodynamics problem. As a key step, the extrapolation and bounce back schemes are used to treat the boundaries. The simulation results have shown a great effectiveness of the implementation of the LBM in the present study. LBM would serve as a promising method for predicting flow and heat transfer characteristics in regenerative cryocoolers, such as pulse tube cryocoolers. In this work, the effects of the characteristic system parameters on the fluid flow and the heat transfer are investigated. Detailed information for system behavior, especially in porous media, are also included in this paper. Ó 2013 Elsevier Ltd. All rights reserved.

Keywords: Lattice Boltzmann method Oscillating flow Heat transfer Numerical simulation Cryogenics

1. Introduction Oscillating flow is a fundamental driving force in Stirling engines and various regenerative cryocoolers. The pulse tube cryocooler (PTC) is a type of regenerative cryocooler with the potential of longtime operation in space, and it can cool space-borne infrared devices to about 80 K in flight vehicles [1]. The schematic of PTC is shown in Fig. 1, which mainly consists of a pressure wave generator, a regenerator, a pulse tube, a hot end heat exchanger, a cold end heat exchanger, an inertance tube and a reservoir [2,3]. The working medium performs a complicated, periodically oscillatory motion in the pulse tube, which is induced by the compressor at one end, and

* Corresponding author. Tel.: þ86 010 82543446. E-mail address: [email protected] (L. Yang). 1359-4311/$ e see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.applthermaleng.2013.01.020

reflected by the other closed end. Besides, the oscillating flow and the fluid heat transfer play an important role in the PTC’s refrigeration performance. In classical studies, the phase relationship and flow resistance of the oscillatory flow in pulse tubes and the complicated porous configurated regenerators have usually been studied experimentally [4,5]. However, the complex microscopic structure of the regenerator’s porous media makes the experimental investigation difficult. Hence, the microscopic thermohydrodynamic characteristics inside the irregular structure of the regenerator could not be easily obtained by experiments. In order to obtain the microscopic mechanism of the oscillating flow and the heat transfer, numerical simulations on such systems may be a promising choice, and can expose some new information. Lattice Boltzmann method (LBM), which is a new and promising numerical technique for simulating fluid flow and modeling physics, is being successfully applied to an increasing number of

Q. Dai, L. Yang / Applied Thermal Engineering 54 (2013) 16e25

research fields in science and engineering [6e9], including fluid flows in porous media [6], multi-phase flow [7], magnetohydrodynamics [8], and nano- and microscale fluid flow [9] problems. Different from conventional Computational Fluid Dynamics (CFD), LBM is based on gas kinetic theory. The simulation by LBM is developed through tracking the evolution of fluid particles’ distribution function. As a mesoscopic method, there are many advantages over the NaviereStokes approach, such as a clear physical background, an easy boundary processing method, involving microscopic interactions and so on [10e12]. The conventional CFD method is difficult to compute the complicated interfaces like porous media. However, since the simple boundary treatment process, LBM can be proposed as a promising method to study the thermal flow in the complicated structure of the cryocoolers, such as the porous medium in a regenerator. Numerous studies of LBM investigation on porous media flows have been reported [13e17]. Some references also studied the thermal flow effects in various kinds of applications, including PEM fuel cell [15], conduits filled with porous media [16], and fractal porous medium [17]. The differences of the present work from them are the oscillatory flow and the application background of cryocoolers. Some researchers adopted the LBM to investigate the oscillating flow and heat transfer in an open channel [18,19] or closed tube [20,21]. Y. Wang [18] and F.K. Meng [19] both investigated the oscillating flows in 2D parallel-plate channel, and obtained some flow field results. F.K. Meng [19] also investigated fluide solid conjugate heat transfer of oscillating flow in the isothermal parallel-plate. Y. Wang [20,21] resonant simulated the oscillations in a closed tube and constructed Implicit-explicit finite-difference LBM model. However, all these studies investigated an empty channel or tube, without specific applying object. The study of F.K. Meng [19] was based on isothermal boundary, and the conjugate heat transfer mainly aimed at the heat transfer between wall and the fluid. In the present work, the oscillatory flow and heat transfer in porous media are investigated, aiming at the study of flow and heat transfer problems in regenerator of regenerative cryocoolers. Since the

  8  i 1h   ðeqÞ > > d d fi r; t  fi r; t ; ¼  ðr þ e t; t þ tÞ  f f > i r; t i < i

su

Another outstanding predominance of this method, due to the kinetic nature, is that the micro-pore structure’s detailed local information of the flow and heat transfer can be obtained, which can be used to investigate macroscopic relations. In order to study the microscopic characteristics of the oscillating flow and heat transfer in the pulse tube, the pulsating flow in a two-dimensional planar channel is simulated by double distribution function LBM, a typical thermal LBM. Then, the oscillating flow and heat transfer characteristics on different flow conditions are presented, including both the differences and similarities. 2. Thermal lattice Boltzmann model 2.1. LBM model Originally, the lattice Boltzmann method can only be applied to isothermal or athermal flows, but as a result of the significance of thermal effects in many engineering applications, energy conservation was taken into account gradually. The double distribution function (DDF) model is one of the widely used thermal LBM models, and its reliability has been confirmed for modeling thermal flows. The basic idea is that the fluid field and the temperature field are described by two independent distribution functions, the single particle density distribution function and the temperature distribution function, respectively. DDF model shows good simulation stability, and can simulate large temperature variation problems with a simple lattice construction. In the present study, the D2Q9 Double Distribution Function (DDF) LBM model developed by Guo et al. [22] is adopted. This DDF LBM model is based on Boussinesq hypothesis. The model applies the improved density equilibrium distribution function (see Eq. (3) below in detail) compared to the traditional D2Q9 model, since still based on the standard D2Q9 lattice of discrete velocity. In addition, it can be applied to both steady and unsteady flows and has a smaller simulation error. For this coupled thermal LBM model, the velocity and temperature evolution equations are respectively defined as following [22]

i ¼ 0; 1; 2; 3; 4; 5; 6; 7; 8

   i > 1h   > ðeqÞ > : Ti ðr þ ei dt; t þ dtÞ  Ti r; t ¼  Ti r; t  Ti r; t ;

sT

previous works are not carried out for regenerators, they cannot explain the flow and heat transfer characteristics inside them. The temporal variation of pressure has either not been described in their works, which is a critical factor for the performance of the cryocooler. In the present work, the influences of dimensionless characteristic parameters like the Womersley and Reynolds numbers on the flow and heat transfer are analyzed systematically.

8 < ð0; 0Þ; ðcos½ði  1Þp=2; sin½ði  1Þp=2Þ$c; ei ¼ : pffiffiffi 2ðcos½ði  5Þp=2 þ p=4; sin½ði  5Þp=2 þ p=4Þ$c;

17

(1) i ¼ 1; 2; 3; 4

where fi is the single particle density distribution function, and Ti is the temperature distribution function. su and sT represent the dimensionless characteristic collision time for velocity and temperature, respectively, which denote the relaxation time to a local equilibrium distribution feq or Teq. In this LBM model, a two-dimensional and nine directional speed square lattice model is used, where ei is the discrete microscopic velocity of a single particle, which can be given by:

i ¼ 0 i ¼ 1; 2; 3; 4 ; i ¼ 5; 6; 7; 8

(2)

18

Q. Dai, L. Yang / Applied Thermal Engineering 54 (2013) 16e25

Fig. 1. Schematic of a pulse tube cryocooler.

where c is the basic speed of the lattice. In this work, we set c ¼ 1 as usual. ðeqÞ is defined as The density equilibrium distribution function fi [22]

ðeqÞ fi

8   p > r0  4s 2 þ r0 s0 u ; i ¼ 0 > > > c > > < p   l 2 þ r0 si u ; i ¼ 1  4 ¼ > c > > > >   > : g p þ r0 si u ; i ¼ 5  8 2 c

P 8 ei fi u ¼ > > > i > > > " # > > >   < c2 X p ¼ r0 fi þ s0 u 4s > is0 > > > > > 4 > P > > Ti :T ¼

(6)

i¼1

(3) Through the ChapmaneEnskog procedure, a multi-scaling expansion, the incompressible NaviereStokes and energy equations can be derived from the lattice Boltzmann model (see Eq. (1)) as follow:

with

" #   ei $u ðei $uÞ2 u2 si u ¼ ui þ  2 c2s 2c4s 2cs

(4)

8 V$u ¼ 0 > >   > > > < vu þ V$ uu ¼ Vp þ nV2 u r0 vt ; >   > > > vT > : þ V$ uT ¼ DV$T vt

(7)

(

4=9; i ¼ 0 1=9; i ¼ 1; 2; 3; 4 . s, l, g are 1=36; i ¼ 5; 6; 7; 8 the model parameters of the density equilibrium distribution ( lþg ¼ s . function, which should meet the requirements: l þ 2g ¼ 1=2 and the weight coefficients wi ¼

ðeqÞ

The temperature equilibrium distribution function Ti given by Ref. [22]

ðeqÞ

Ti

¼

  T e $u 1þ2 i2 ; 4 c

i ¼ 1; 2; 3; 4:

can be

(5)

The macroscopic fluid velocity u, the pressure p, and the temperature T can be obtained by their respective distribution functions as following [22]:

where the fluid kinematic viscosity is defined as n ¼ 1/3c2(su  1/2) dt and the thermal diffusivity is defined as D ¼ 1/2c2(sT  1/2)dt. The corresponding Prandtl number is given by Pr ¼ n/D ¼ 2/3(2su  1)/ (2sT  1). 2.2. Boundary conditions used in the simulations In the lattice Boltzmann method, the basic evolution variables of fi and Ti are not given directly on the boundaries. They can only be obtained by appropriate boundary treatment schemes. In this simulation, the boundary treatment schemes used in the present research are non-equilibrium extrapolation scheme [23] and bounce back scheme [24]. They are used to respectively treat the simulation domain boundaries and the boundaries of complicated porous media wall, as shown in Figs. 2 and 3. Bounce back scheme is a simple scheme very suitable for complicated boundary

(i-1,2)

g6

(i+1,2)

(i,2)

g8

g7

g4 g2

g5

(i,1)

(a) Non- equilibrium extrapolation scheme Fig. 2. Boundary treatment schemes.

(b) Bounce back scheme

Q. Dai, L. Yang / Applied Thermal Engineering 54 (2013) 16e25

19

The analytical solution of the flow velocity field is expressed as [12]:

     A cos½lð2y=W  1Þ iut 1 e u y; t ¼ Re i u cos l

Fig. 3. Physical model of the oscillating flow in porous media.

treatment like porous media, which is also the reason of LBM’s significant advantage of easy boundary processing. The accuracy of the non-equilibrium extrapolation scheme is second order in space, matching the accuracy of the lattice Boltzmann model. The distribution function of boundary nodes O (see Fig. 2(a)) can be given by Ref. [23]

       

ga O; t ¼ gaeq O; t þ ga B; t  gaeq B; t

(8)

where g represents the density or temperature distribution function (f or T), a indicate the lattice direction, varying from 0 to 8. O and B both represent the mesh nodes, while t represents the time. And analogously, we can obtain the distribution functions of all the other boundary nodes, including velocity and temperature variables. For the complicated porous media wall, the bounce back scheme is adopted in the calculation of flow field. In the bounce back scheme, the particle’s distribution function is assumed to be reversal after the particle collide the boundary. As shown in Fig. 2(b), the unknown distribution function on the boundary nodes can be given by Ref. [24]

    g2;5;6 i; 1 ¼ g4;7;8 i; 1

(9)

where g represents the density distribution function f, i represent the mesh node, 1 represent the current step, while 2 represent the next time step. Analogously, we can obtain the distribution functions of all the other boundary nodes. 3. Numerical simulations and discussion 3.1. Physical model The unsteady oscillating flow in a two-dimensional channel is defined in the region ½ðX; YÞj0  X  100; 0  Y  50. The physical model and boundary conditions are shown in Fig. 3, where L and W are the geometry length and width of the channel respectively. We first investigated the oscillating flow and heat transfer in an empty planar channel, and then in a channel placed with some porous media. For the empty channel model, the boundary plates are at constant temperature Tw, and the temperature of working gas in the inlet and outlet are at uniform temperature Ti. The flow is driven by a periodic pressure wave between the entrance and the exit of the channel, which is given by P ¼ Pm þ Pasin ut. The pressure gradient is given by dP/dx ¼ Re[Aeiwt] ¼ Acos(ut), with an amplitude A and an angular frequency u. The mean pressure at the outlet of the region is Pout, and that at inlet is Pm ¼ Pout þ ALcos(ut). The lattice size is Nx  Ny ¼ 100  50. The values of relaxation times su and sT, are based on the limitations of su ¼ 3nþ0.5 and Pr ¼ 2=3ð2su  1Þ=ð2sT  1Þ, where n is the fluid kinematic viscosity.

(10)

pffiffiffiffiffiffiffiffi with l2 ¼ ia2, a ¼ W=2 u=n, and a represents the Womersley number of the pulsating flow. Womersley number is a dimensionless expression about the pulsating flow frequency divided by the dynamic viscosity. Re[c] (Real) is the real part of the complex number c. n and m are the dynamic and the kinetic viscosity respectively. The maximum velocity in the center of the channel is umax ¼ AW2/8n, and other parameters like the Reynolds number and the Mach number are defined by Re ¼ Wumax/n, Ma ¼ umax =cs ¼ AW 2 =8ncs . The spatial variation of pressure (or density) should not be too large, in order to make sure the propagation of pressure wave or density fluctuation in the fluid is regarded as instantaneous. To archive this, the period of oscillating flow must be greater than L/cs, which is the time taken by the sound wave to travel the distance L. So the driving frequency should not be too high [12]. The initial state of velocity field u is always set to be zero in the whole system. The initial state of temperature field is set to Ti, and the initial pressure is Pout. In the two-dimensional channel, the inlet and outlet of the fluid flow (left and right boundaries), as shown in Fig. 3, are flow boundaries, and the top and bottom parallel walls are solid boundaries. Here we adopt the non-equilibrium extrapolation scheme to process all these boundaries [23]. Based the above model, an improved model was constructed. According to the structure of mesh screen packing in the regenerator, an innovative porous media model was constructed. A modified simulation model was developed by placing a certain amount of porous media in the center of the planar plate, as shown in Fig. 3, with the plate wall boundary adiabatic or isothermal. It is aimed at investigating the oscillating flow and heat transfer features through porous media in the regenerator. The porous media was initially constructed by self-programming, and the porosity of the porous media is 0.7, which is very common in regenerator of cryocooler. We adopt the bounce back [24] to treat the complicated boundaries of the interfaces between solid and fluid in the porous media, as elated in the upper part. For the flow field, the solid mesh of solid porous media was not calculated. But for temperature field, both the solid and fluid mesh were calculated. In the process of calculation, the solid and fluid mesh adopted the different relaxation time in the evolution equation (Eq. (1)), sT_s and sT_f respectively. 3.2. Results and discussion In the simulation of the empty channel, with the amplitude A fixed at 1.0  104, we calculated the velocity and temperature field for different Womersley number, from less than 1 to 10, and analyzed their dependence on it. Three important parameters, including period, Reynolds number, and angular frequency, are shown in Table 1. The computational curves are shown from Figs. 4e9, and the results will be discussed in the following part in detail. Table 1 Corresponding parameters at three different periods. Cases

Period (Tk)

Re number (Re)

Angular frequency (u ¼ 2p/Tk)

Womersley pnumber ffiffiffiffiffiffiffiffi (a ¼ W=2 u=n)

a b c d

8000 1000 200 100

5 10 8 10

0.000785 0.00628 0.0314 0.0628

0.937 < 1 3.15 > 1 6.66 > 1 9.97 > 1

20

Q. Dai, L. Yang / Applied Thermal Engineering 54 (2013) 16e25

1.4

0.4

1.2

1/8T 5/8T

1.0

2/8T 6/8T

3/8T 7/8T

4/8T 8/8T

0.8 0.4

U/Umax

U/Umax

2/8T 6/8T

0.2 0.0 -0.2

30

40

0.1 0.0

-0.2

-0.6 -0.8

-0.3 0

10

20

30

40

50

0

10

20

(a)

α = 0.937

(b) α = 3.15

0.20 1/8T 5/8T

0.15

2/8T 6/8T

3/8T 7/8T

50

Y

Y

4/8T 8/8T

0.16 0.14 0.12 0.10 0.08 0.06 0.04 0.02 0.00 -0.02 -0.04 -0.06 -0.08 -0.10 -0.12

1/8T 5/8T

2/8T 6/8T

3/8T 7/8T

4/8T 8/8T

U/Umax

0.10

U/Umax

4/8T 8/8T

-0.1

-0.4

0.05 0.00

-0.05 -0.10 -0.15

3/8T 7/8T

0.2

0.6

-1.0

1/8T 5/8T

0.3

0

10

20

Y

30

40

50

(c) α = 6.66

0

10

20

(d)

Y

30

40

50

α = 9.97

Fig. 4. Velocity distribution at X ¼ 0.5L profile in a whole period for different a.

Fig. 4 shows the velocity fluctuation at the middle of the channel (X ¼ 0.5L) in a whole period for different Womersley numbers a. When the Womersley number is close or less than 1, the frequency of pulsation is sufficiently slight so that a parabolic velocity profile has enough time to be fully developed during each cycle. In this case the flow will be nearly in phase with the pressure wave, and can be approximated by Poiseuille’s law at each time step. The results agree well with the theoretical analyses (see Eq. (10)). As seen in Fig. 4(a), the velocity along the centerline is closer to the maximum velocity and always higher than the velocity near the wall at any point. However, when the Womersley number is much higher than 1, as shown in Fig. 4(b)e(d), the velocity annular effect show up more and more distinctly. The velocity annular effect is formed as a result of both the wall viscous phenomenon and the fluid flow inertial force in the middle of the channel. This is consistent with experimental results measured by Richardson [25], who first found the annular effect of oscillating flow in the pipe. It should be noted that the velocity distribution profiles are relatively flat or plug-like, especially in Fig. 4(c) and (d), where the maximum velocity appears near the wall, but not at the center of the channel. Sometimes we can also find the interesting phenomenon that the velocity direction is not the same along the same cross section at certain moments, which means that a velocity reversal is taking place. Moreover, the results of the velocity field agree well with the analytical solution. In summary, the velocity profile and shape have

close relationship with the Womersley number and the period of driving pressure. It can be also noted that the velocity profiles are periodically symmetrical. The temperature distribution at the middle cross section of the channel (X ¼ 0.5L) in a whole period and for several Womersley numbers is presented in Fig. 5. The comparison of the temperature distribution at different periods, shows that when the period Tk is 100 and a > 1, the temperature distribution is annular, and the maximum value of temperature emerged nearby the wall. While the period Tk is 1000 and a > 1, the temperature distribution is annular or parabolic. When the period Tk is 8000 and a < 1, the temperature distribution is parabolic and the fluid temperature is lower than Tw. That is to say, the annular effects can also be found in temperature distributions when a > 1. Fig. 6 presents the periodical pressure wave along the centerline in a whole period, which is the driving force of the flow field. We can see that the pressure changes periodically with time, and the pressure difference between the inlet and outlet is positive for half the period and negative for the other half. In addition, the pressure is nearly linear along the X-axis direction at each moment. These phenomena agree well with the previous numerical simulations [18] and theoretical results. In fact, the pressure always keeps almost constant along the same vertical profile. Fig. 7 illustrates the pressure fluctuation in a whole period at the middle cross section of the channel on the

Q. Dai, L. Yang / Applied Thermal Engineering 54 (2013) 16e25

21

1.4 1.05

1/8T 5/8T

3/8T 7/8T

4/8T 8/8T

1.2

0.95 0.90 0.85

3/8T 7/8T

4/8T 8/8T

1.1 1.0 0.9

0.7 0

10

20

30

40

0.6

50

0

Y

(a)

( 1/8T 5/8T

1.4

2/8T 6/8T

3/8T 7/8T

10

20

(b)

)

1.6

1.6

4/8T 8/8T

Y

30

( 1/8T 5/8T

2/8T 6/8T

40

50

) 3/8T 7/8T

4/8T 8/8T

1.4

(T-Ti)/(TW-Ti)

(T-Ti)/(TW-Ti)

2/8T 6/8T

0.8

0.80 0.75

1/8T 5/8T

1.3

(T-Ti)/(TW-Ti)

(T-Ti)/(TW-Ti)

1.00

2/8T 6/8T

1.2

1.2

1.0

1.0

0.8

0.8

0.6

0.6 0

10

20

30

Y

(c)

40

50

0

10

20

30

40

50

Y

(

)

(d)

9.97 ( Tk

100 )

Fig. 5. Temperature distribution at X ¼ 0.5L profile in a whole period for different a.

vertical profile (X ¼ 0.5L) for a ¼ 0.937 and 6.66. We can see that the pressure maintains nearly the same value along the vertical profile. When the period Tk is large and a < 1 (see Fig. 7(a)), the pressure profiles are extremely flat, and the pressure variation

1/8T 5/8T

(P-Pout)/(Pin-Pout)

1.0

2/8T 6/8T

3/8T 7/8T

4/8T 8/8T

0.5 0.0

-0.5 -1.0 0

20

40

X

60

80

100

Fig. 6. Pressure distributions in a whole period along the centerline at Y ¼ 0.5W profile.

along the vertical profile at one moment is near to 0. However, when the period Tk is small and a > 1, there is a little fluctuation on the pressure profiles (see Fig. 7(b)), that means there exists a pressure gradient along the vertical direction in the channel, which is the reason for the velocity annular effects. In other words, the high frequency of the pressure wave will results in vertical section pressure fluctuation. Besides, the pressure variation on the vertical profile in a period will enlarge with increasing Womersley number. Periodically variations of velocity, pressure and temperature at three positions are presented in Fig. 8. The three positions are all on the vertical profile of X ¼ 0.5L, with Y ¼ 0.1,0.5,0.9W respectively. For the velocity distribution, as shown in Fig. 8(a), the velocity changes periodically with time. Due to the spatially symmetry of the top and the bottom of the channel flow, the velocity and temperature variation curves for Y ¼ 0.1W and Y ¼ 0.9W are coincided with each other. The amplitude of the velocity at Y ¼ 0.5W is larger than that at Y ¼ 0.1W or 0.9W, which again indicate that the maximum velocity of the channel is in the center when Womersley number is small. We can also see that the pressure variation curves are nearly the same for the three positions, again denoting the phenomenon shown in Fig. 7. Moreover, the velocity and temperature in different locations along the same vertical profile are in phase with each other.

22

Q. Dai, L. Yang / Applied Thermal Engineering 54 (2013) 16e25

0.8 1/8T 5/8T

2/8T 6/8T

3/8T 7/8T

1.0

4/8T 8/8T

0.5

0.4 0.2

0.0

0.0 -0.5

-0.2 -0.4

-1.0

0

10

(a)

Y

0.937 ( Tk 1/8T 5/8T

1.0

20

2/8T 6/8T

30

40

0

50

8000 )

3/8T 7/8T

50

100

150

200

t( º )

250

300

350

0.6 X=0.5 Y=0.1 Y=0.5 Y=0.9

0.4

4/8T 8/8T

(P-Pout)/(Pin-Pout)

-0.6

(P-Pout)/(Pin-Pout)

X=0.5 Y=0.1 Y=0.5 Y=0.9

U/Umax

(P-Pout)/(Pin-Pout)

0.6

0.5

0.2 0.0

-0.2

0.0

-0.4 -0.5

-0.6 0

10

20

Y

30

40

0

50

100

50

150

200

t( º )

250

300

350

1.00

(b)

6.66 ( Tk

200 ) 0.95

Fig. 9 shows the periodically variations of velocity, pressure and temperature at three different locations along the centerline of the channel, where Y ¼ 0.5W, X ¼ 0.1,0.5,0.9L respectively. The velocity variations along the centerline at different positions are equal with each other, which means that the flow is the same along the X-axis. The three pressure variations along the centerline keep changing with the same phase, while the amplitude gradually decreases from the inlet to the outlet. The temperature variations at the three positions along the centerline have different phases, and the one at X ¼ 0.1L has approximately the opposite phase with that at X ¼ 0.9L. In addition, the average temperature for X ¼ 0.5L is higher than at the other two places. We also investigated the flow and heat transfer in porous media, which we constructed by programming according to the features of practical mesh screens in cryocoolers. Fig. 10(a) shows the streamline chart of oscillating flow through porous media at a time in a period. The fluid in the channel flow out from the two ends, and the zero velocity profile moves from left to right. It is noted that there are small vortexes near the surface of solid-phase porous media. In actual, the streamline keep changing in a period. This phenomenon is due to the constantly changing pressure boundary condition. The combined action of the periodically changed inertial

(T-Ti)/(TW-Ti)

Fig. 7. Pressure distribution at X ¼ 0.5L profile for different a.

0.90 0.85 X=0.5 Y=0.1 Y=0.5 Y=0.9

0.80 0.75

0

50

100

150

200

t( º )

250

300

350

Fig. 8. Periodically variations of U, P, and T at X ¼ 0.5L, Y ¼ 0.1,0.5,0.9W (a ¼ 0.937).

motion and the viscous force between fluid and solid caused the typical characteristics of oscillating flow in porous media. The temperature contour of the working medium and the porous media at a time in a period was presented in Fig. 10(b). In a period of oscillating flow, the solid porous media continuously absorb heat from or release heat to the working medium. Fig. 10(b) shows a representative time when solid porous media absorb heat from the fluid.

Q. Dai, L. Yang / Applied Thermal Engineering 54 (2013) 16e25

23

Fig. 10. (a) Streamline chart of oscillating flow through porous media at a time 10. (b) Temperature contour in the porous media at a time.

respectively. The dimensionless velocity and temperature are defined as U/Umax and (T  Ti)/(Tw  Ti), respectively. The plate wall is isothermal in Fig. 12. In actual, there is almost no difference for the flow between isothermal and adiabatic boundary conditions. From Fig. 12, we can see that the velocity always keep zero in some places, corresponding to the solid-phase of porous media. The flow in the pores keeps changing periodically with time, which is socalled oscillating flow. And the oscillating flow in each pore can be seen as a microscopic oscillating flow in planar channel, as shown in Fig. 4. It means that the characteristics of oscillating flow in two parallel plane appear in the pore of the porous media. For the temperature distribution with isothermal boundary in Fig. 13, the temperature at any position keeps changing periodically with the time. The fluctuations in the pores are more volatile than that of solid-phase, which is as a result of the higher heat capacity of solid-

Fig. 9. Variations of U, P, and T at Y ¼ 0.5W, X ¼ 0.1,0.5,0.9L (a ¼ 0.937).

We also compute the regenerator connected by a pulse tube, as shown in Fig. 11. We can see that the flow in the entrance part to the regenerator in the pulse tube is very nonuniform, which is a flow loss to the cryocooler. Figs. 12 and 13 illustrate the velocity and temperature distribution on central cross section (X ¼ 0.5L profile) in a whole period,

Fig. 11. Streamline chart around regenerator and pulse tube at a time.

24

Q. Dai, L. Yang / Applied Thermal Engineering 54 (2013) 16e25

0.14 0.12

1/8T 5/8T

0.10

2/8T 6/8T

3/8T 7/8T

1.0

4/8T 8/8T

Dimensionless temperature

Dimensionless velocity

0.08 0.06 0.04 0.02 0.00 -0.02 -0.04 -0.06 -0.08 -0.10 0

20

40

60

Y axis

80

as/af=1 as/af=5

0.8

as/af=10 as/af=17.213

0.6

0.4

0.2

0.0

100

as/af=0.656

0

20

40

60

80

100

Y axis

Fig. 12. X-directional velocity’s periodical distribution on X ¼ 0.5L profile.

Fig. 15. Temperature distribution on X ¼ 0.5L profile, at 1/8T, with different thermal diffusivity ratio (isothermal boundary).

5.0

1/8T 5/8T

4.0

2/8T 6/8T

3/8T 7/8T

4/8T 8/8T

3.5 3.0

Dimensional Temperature

Dimensionless temperature T

4.5

2.5 2.0 1.5 1.0 0.5 0.0

-0.5 0

20

40

Y axis

60

80

100

Fig. 13. X-directional temperature’s periodical distribution on X ¼ 0.5L profile (isothermal boundary).

phase than that of the fluid. The working fluid periodically absorb or release heat with the solid-phase porous media. The flow keeps the same either the plates is isothermal or adiabatic. Fig. 14 shows the temperature distribution with adiabatic boundary. It can be easily found that the temperature distribute differently in the area near the plates from that of isothermal.

1.0020

1/8T 5/8T

1.0016

2/8T 6/8T

3/8T 7/8T

4/8T 8/8T

1.0012

T/Tin

1.0008 1.0004 1.0000 0.9996 0.9992 0.9988 0

20

40

60

80

100

Y axis Fig. 14. X-directional temperature’s periodical distribution on X ¼ 0.5L profile (adiabatic boundary).

1.0014 1.0012 1.0010 1.0008 1.0006 1.0004 1.0002 1.0000 0.9998 0.9996 0.9994 0.9992 0.9990 0.9988 0.9986

a /a =0.656 a /a =1 a /a =5 a /a =10 a /a =17.213 0

20

40

60

80

100

Y axis Fig. 16. Temperature distribution on X ¼ 0.5L profile, at 1/8T, with different thermal diffusivity ratio (adiabatic boundary).

The material’s thermal diffusivity is defined as a ¼ k/rcp, where k is the conductivity, and rcp represents the heat capacity. Thermal diffusivity ratio between solid-phase and fluid is defined as as/af, where as and af respectively represent the thermal diffusivity of solid and fluid. For stainless-steel woven wire screens, the ratio equals to 0.656, while for copper woven wire screens, the ratio equals to 17.213. Figs. 15 and 16 show the temperature distributions at 1/8T, with different solidefluid thermal diffusivity ratio, on X ¼ 0.5L profile. The two figures respectively represent the situation of isothermal boundary and adiabatic boundary. From Fig. 15, we can see the temperature of fluid and porous media are lower than the two parallel planes, which is because of the given constant temperature boundary condition. The heat is transferred from the solid planes to the channel. With the increasing of thermal diffusivity ratio, the temperature curves are rising constantly, which means the higher the ratio is, the better the heat transfer performance is. From Fig. 16, the temperature fluctuates more and more feebly with the increasing of thermal diffusivity ratio. The ratio influences a lot on the solid temperature, but not on the fluid temperature. In other words, with the increasing of the thermal diffusivity of solid, the solid packing can transfer heat more quickly with the ambient fluid. On the whole, it can be concluded that the thermal diffusivity ratio plays a big part in the heat transfer between fluid and solid.

Q. Dai, L. Yang / Applied Thermal Engineering 54 (2013) 16e25

4. Conclusion In the present study, a typical double distribution function thermal Lattice Boltzmann model is proposed to simulate the thermal oscillating flow both in planar planes and in porous media. Since its emergence, the LBM is applied to an increasing number of research fields in science and engineering, and has obtained some good achievements in mesoscopic scale, especially in some special fields where conventional CFD is not appropriate. Porous media flow is a typical example. Since we cannot obtain the real flow field of porous media by traditional NaviereStokes method, LBM is a better choice. The introduction and implementation of LBMs to the investigation of regenerative cryocoolers is a totally new trial. The simulation results show that: via the LBM simulation, the microscopic characteristics in planar planes and porous media can be obtained. We presented the dependence of flow field and temperature field at different Womersley numbers. It can be noted that the velocity and temperature profiles are mainly influenced by the Womersley number. The velocity annular effect emerges when Womersley number is more than 1, and becomes more and more distinct with the increase of the Womersley number. In the channel, the pressure along the same vertical profile always stays nearly constant. The results agree well with the previous experimental investigation results. For the flow in porous media, we found that there are some little vortexes near the surface of solid-phase porous media at some time. The working gas periodically absorb or release heat with the solid-phase porous media. The thermal diffusivity ratio between solid and fluid also have great influences on heat transfer. The successful implementation of LBM shows that this method is useful for simulating oscillating flow and heat transfer in regenerative cryocooler systems, and may be further developed as a powerful tool for investigating and improving the design of pulse tube cryocoolers. Acknowledgements This research is supported by National Science Foundation of China (No. 50890181, No. 51076160). References [1] J.T. Liang, Y. Zhou, W.X. Zhu, Study on miniature pulse tube cryocooler for space application, Cryogenics 40 (2000) 229e233. [2] L.W. Yang, G. Thummes, High frequency two-stage pulse tube cryocooler with base temperature below 20 K, Cryogenics 45 (2005) 155e159.

25

[3] J.T. Liang, A. Ravex, P. Rolland, Study on pulse tube refrigeration part 1: thermodynamic nonsymmetry effect, part 2: theoretical Modelling, part 3: experimental verification, Cryogenics 36 (1996) 87e106. [4] H.L. Chen, J.L. Yang, L.W. Yang, Dynamic characteristics of oscillating flow in pulse tube coolers, in: The Proc. of the CEC/ICMC 2007, pp. 715e719. [5] X.L. Wang, M.G. Zhao, J.H. Cai, Experimental analysis of the oscillating flow characteristics for high frequency regenerators, in: The 20th Int. CEC, Beijing, China 2004, pp. 201e204. [6] Z.L. Guo, T.S. Zhao, Lattice Boltzmann model for incompressible flows through porous media, Phys. Rev. E 66 (2002) 036301e036304. [7] X. He, S. Chen, R. Zhang, A lattice Boltzmann scheme for incompressible multiphase flow and its application in simulation of RayleigheTaylor instability, J. Stat. Phys. 152 (1999) 642e663. [8] H. Yamaguchi, X.R. Zhang, X.D. Niu, K. Yoshikawa, Thermomagnetic natural convection of thermo-sensitive magnetic fluids in cubic cavity with heat generating object inside, J. Magn. Magn. Mater. 322 (2010) 698e704. [9] D. Raabe, Overview of the lattice Boltzmann method for nano- and microscale fluid dynamics in materials science and engineering, Model. Simulat. Mater. Sci. Eng. 12 (2004) 13e46. [10] S. Succi, Lattice Boltzmann Equation for Fluid Dynamics and Beyond, first ed., Clarendon, Oxford, 2001, pp. 40e43. [11] S.Y. Chen, G.D. Doolen, Lattice Boltzmann method for fluid flows, Annu. Rev. Fluid Mech. 30 (1998) 329e364. [12] X.Y. He, L.S. Luo, Lattice Boltzmann model for the incompressible Naviere Stokes equation, J. Stat. Phys. 88 (1997) 927e988. [13] Y.L. He, Q. Li, Y. Wang, Lattice Boltzmann method and its applications in engineering thermophysics, Chin. Sci. Bull. 54 (2009) 4117e4134. [14] N. Jeong, D.H. Choi, C.L. Lin, Estimation of thermal and mass diffusivity in a porous medium of complex structure using a lattice Boltzmann method, Int. J. Heat Mass. Transfer 51 (2008) 3913e3923. [15] J. Park, M. Matsubara, X. Li, Application of lattice Boltzmann method to a micro-scale flow simulation in the porous electrode of a PEM fuel cell, J. Power Sources 173 (2007) 404e414. [16] H. Shokouhmand, F. Jam, M.R. Salimpour, Simulation of laminar flow and convective heat transfer in conduits filled with porous media using Lattice Boltzmann Method, Int. Commun. Heat Mass. Transfer 36 (2009) 378e384. [17] J. Cai, X.L. Huai, Study on fluid-solid coupling heat transfer in fractal porous medium by lattice Boltzmann method, Appl. Therm. Eng. 30 (2010) 715e723. [18] Y. Wang, Y.L. He, G.H. Tang, Simulation of two-dimensional oscillating flow using the lattice Boltzmann method, Int. J. Mod. Phys. C 17 (2006) 615e630. [19] F.K. Meng, M. Wang, Z.X. Li, Lattice Boltzmann simulations of conjugate heat transfer in high-frequency oscillating flows, Int. J. Heat Fluid Flow 29 (2008) 1203e1210. [20] Y. Wang, Y.L. He, Q. Li, Numerical simulations of gas resonant oscillations in a closed tube using lattice Boltzmann method, Int. J. Heat Mass. Transfer 51 (2008) 3082e3090. [21] Y. Wang, Y.L. He, J. Huang, Impliciteexplicit finite-difference lattice Boltzmann method with viscid compressible model for gas oscillating patterns in a resonator, Int. J. Numer. Methods Fluids 59 (2009) 853e872. [22] Z.L. Guo, B.C. Shi, C.G. Zheng, A coupled lattice BGK model for the Boussinesq equations, Int. J. Numer. Methods Fluids 39 (2002) 325e342. [23] Z.L. Guo, C.G. Zheng, B.C. Shi, Non-equilibrium extrapolation method for velocity and boundary conditions in the lattice Boltzmann method, Chin. Phys. 11 (2002) 0366e0374. [24] S.Y. Chen, D. Martinez, R.W. Mei, On boundary conditions in lattice Boltzmann methods, Phys. Fluids 8 (9) (1996) 2257e2536. [25] E.G. Richardson, E. Tyler, Transverse velocity gradient near the mouths of pipes in which an alternating or continues flow air is established, Proc. Phys. Soc. London 42 (1929) 1.