Meshless kinetic upwind method for compressible, viscous rotating flows

Computers & Fluids 46 (2011) 325–332

Contents lists available at ScienceDirect

Computers & Fluids j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m p fl u i d

Meshless kinetic upwind method for compressible, viscous rotating flows A.K. Mahendra a,⇑, R.K. Singh b, G. Gouthaman a a b

Machine Dynamics Division, Bhabha Atomic Research Centre, Mumbai 400 085, India Reactor Safety Division, Bhabha Atomic Research Centre, Mumbai 400 085, India

a r t i c l e

i n f o

Article history: Received 30 April 2010 Received in revised form 15 October 2010 Accepted 24 October 2010 Available online 4 November 2010 Keywords: Meshless Rotating viscous flow Kinetic method KFVS Least squares Slip flow Chapman-Enskog distribution

a b s t r a c t This paper presents the split-stencil least square kinetic upwind method for Navier–Stokes (SLKNS) solver using kinetic flux vector splitting (KFVS) scheme with Chapman-Enskog distribution. SLKNS solver operates on an arbitrary distribution of points and uses a novel least squares method which differs from the normal equations approach as it generates the non-symmetric cross-product matrix by selective splitting of the set of neighbours to avoid ill-conditioning. SLKNS also uses the axi-symmetric formulation of the Boltzmann equation and kinetic slip boundary condition. SLKNS is capable of capturing weak secondary flows as well as features of strong rotation characterized by steep density gradient and thin boundary layers towards the peripheral region with a rarefied central core.  2010 Elsevier Ltd. All rights reserved.

1. Introduction For rapidly rotating viscous compressible flows the gas is held under rigid-body rotation characterized by an exponential density rise in the radial direction towards the periphery with thin boundary layers and a rarefied inner core [4,6]. Numerical studies performed by researchers Dickson and Jones [11], Park and Hyun [23] and Babarsky et al. [2] focused more on perturbative flows arising out of basic state of rigid body rotation. Most of the researchers have described such a rapidly rotating flows using continuum hydrodynamics eventhough a large portion of the volume is rendered rarefied near to the free molecular region. It was first shown by Johnson and Stopford [17] that flow predicted by kinetic theory differs from the results obtained using continuum hydrodynamics. Sharipov and Kremer [27] and Sharipov et al. [26] have investigated transport phenomena through a fluid confined between two coaxial cylinders over a wide range of gas rarefaction and observed that results based on kinetic approach do not follow from the Navier–Stokes and Fourier equations of continuum mechanics. Taheri and Struchtrup [29] have investigated the effects of rarefaction in microflows between two rotating coaxial cylinders using continuum approach and regularized 13-moment (R13) equations. In this paper we have carried out numerical simulation of rapidly rotating flows for practical engineering problems using kinetic ⇑ Corresponding author. Tel.: +91 22 25593611; fax: +91 22 25505151. E-mail address: [email protected] (A.K. Mahendra). 0045-7930/$ - see front matter  2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.compfluid.2010.10.015

flux vector scheme (KFVS) coupled with slip boundary conditions. KFVS was proposed by Deshpande [10], and Mandal and Deshpande [20] used it to solve Euler problem. KFVS is based on moment-method strategy which involves two levels, the Boltzmann level (upwind implementation) and the Euler level (level at which the state update operates). KFVS scheme for the viscous flows was proposed by Chou and Baganoff [9]. Prendergast and Xu [25] proposed a scheme based on the BGK simplification of the Boltzmann equation and May et al. [21] further modified and improved it. Ilgaz and Tuncer [16] used a gas-kinetic BGK scheme for the solution of 3-D viscous flows on unstructured hybrid grids. KFVS for viscous flows applies Courant splitting at the Boltzmann level, with moment-method strategy using Chapman-Enskog distribution to obtain split Navier–Stokes equations. Since the task of generating suitable grid for a complex multibody configuration can be very tedious hence we have used meshless scheme. The present approach is similar to least square kinetic upwind method (LSKUM) developed by Ghosh and Deshpande [13] for Euler flows. Mahendra [18] and Mahendra et al. [19] used KFVS for viscous flows using LSKUM. Ill-conditioning of the least square problem was observed when LSKUM operates on the stretched distribution of points near the boundary required for capturing viscous flow features. The prime motivation of the paper is to develop a meshless method using KFVS to solve a least square problem by avoiding the ill-conditioning and retaining the simplicity of normal equations approach. The present method split-stencil least square kinetic upwind method for Navier–Stokes (SLKNS) is a novel method which differs from normal equations

326

A.K. Mahendra et al. / Computers & Fluids 46 (2011) 325–332

approach as it generates the non-symmetric cross-product matrix by selective splitting of the set of neighbours to avoid ill-conditioning encountered while using highly stretched distribution of points in the boundary layers. SLKNS also uses the axi-symmetric formulation of the Boltzmann equation and kinetic slip boundary condition. Further, SLKNS is able to handle strong rotating flows while capturing the weak secondary flow features. SLKNS also captures features characterized by steep density gradient and thin boundary layers towards the peripheral region with a rarefied central core. 2. Kinetic flux vector splitting for rotating viscous flows The Boltzmann transport equation describes the transient molecular distribution function, f ð~ x; ~ v ;I;tÞ : RN  RN  Rþ  Rþ ! Rþ . For conservation of total energy instead of translational energy alone an additional internal energy variable I 2 Rþ is added as polyatomic gas consists of particles with additional degree of freedom. For a gas in absence of external force and without internal degrees of freedom, the Boltzmann equation is as follows:

@f þ r~x  ð~ v f Þ ¼ Jðf ; f Þ @t

ð1Þ

where ~ x, the position vector and velocity vector, ~ v of molecules are given in RN . The right hand side factor J(f, f) is the collision integral describing the binary collision amongst the molecules which vanishes in the Euler limit. The moment of a function, W ¼ Wð~ v ; I; tÞ : RN  Rþ  Rþ ! R is defined as the inner product

Wð~ x; tÞ ¼ hWð~ v ; I; tÞ; f i 

Z

Z

RN

Rþ

Wð~ v ; I; tÞf ð~x; ~ v ; I; tÞdI d~ v

ð2Þ



 Using moment vector function defined as W ¼ 1; ~ v ; I þ 12 ~ vT~ v T the moments of the Boltzmann equation are as follows:

½q; q~ u; qET ¼ hW; f i 

Z Rþ

RT c1

Z

@f1 v x þ jv x j @f1 v x  jv x j @f1 v y þ jv y j @f1 v y  jv y j þ þ þ þ 2 2 2 2 @t @x @x @y @f1 v z þ jv z j @f1 v z  jv z j @f1  þ þ ¼0 2 2 @y @z @z

ð6Þ

Taking W moment leads to split flux Navier–Stokes equations.

@U @GX  @GY  @GZ  þ þ þ ¼0 @t @x @y @z

ð7Þ

where state vector, U = [q, qux, quy, quz, qE]T and GX±, GY±, GZ± are the split KFVS fluxes. 2.1. Treatment for strong rotation The velocity distribution associated with this rigid body rotation [11], fRB is a special case of Chapman-Enskog distribution as shear stress terms disappear and distribution becomes a Maxwellian [19,27]. Boltzmann equation for any perturbation near this rigid body rotation with f1 as the Chapman-Enskog distribution can be written as

@ðf1  fRB Þ @ðf1  fRB Þ @ðf1  fRB Þ @ðf1  fRB Þ þ vx þ vy þ vz ¼0 @t @x @y @z

ð8Þ

The subscript RB represents the state of rigid body rotation. The steep gradients observed in stationary inertial frame now appear to be a weak perturbation. Taking W moment leads to Navier– Stokes equations in the following special form

@ðU  U RB Þ @ðGX  GX RB Þ @ðGY  GY RB Þ @ðGZ  GZ RB Þ þ þ þ ¼0 @t @x @y @z ð9Þ

RN

Wf ð~ x; ~ v ; I; tÞdI d~ v

ð3Þ

1 ~T ~ ðu uÞ. 2

The Chapman-Enskog expansion of the where E ¼ þ Boltzmann equation with Knudsen number as parameter will give higher order constitutive relationship for shear and heat transfer terms [8,3]. Scaling the temperature gradient term in the Chapman-Enskog expansion by the Prandtl number [21] leads to a first-order Chapman-Enskog term with the correct Prandtl number as

  /1;s /1;q f1 ¼ f0 1   p p

ð4Þ

where f0 is the Maxwellian distribution [20] and /1,s and /1,q are Chapman-Enskog polynomial associated with the shear stress tensor and heat flux vector given as follows:

"

# " # 5c  7 I 5c  7 I 2 2 /1;s ¼ þ bcx sxx þ þ bcy syy   2ð5  3cÞ 4bI2o 2ð5  3cÞ 4bI2o " # 5c  7 I þ þ bc2z szz þ 2bcx cy sxy  2ð5  3cÞ 4bI2o

where state vector URB and fluxes GXRB, GYRB and GZRB are calculated using the primitive variables at rigid body rotation. Dickson and Jones [11] used a small perturbation defined by Rossby number e ¼ Dxx to expand the solution as a Taylor series in e about rigid body rotation. The local Rossby number, e based on z-component of vorticity vector, Xz is defined as

e¼

Dx

x

  ¼ 1 



Xz  ðXz ÞRB 

ð10Þ

This local Rossby number is used as a measure of departure from the rigid body solution, for example if e < es solver switches to this special form of Navier–Stokes equations. Mahendra et al. [19] based on numerical experiments with 3-D Navier Stokes using LSKUM have observed accuracy with es = 0.1. 2.2. Kinetic schemes for axi-symmetric rotating flows The velocity discretization of the transport operator is no longer trivial in the cylindrical system as inertia terms are velocity derivatives of the distribution function. It should be noted that characteristic curves of this form of Boltzmann equation are curves of R4 and certainly more complex compared to the Cartesian form [22]. Sugimoto and Sone [28] and Mieussens [22] have dealt with completely conservative form of Boltzmann equation. The Boltzmann equation cast in the cylindrical coordinate system [27] can be written as:

þ 2bcx cz sxz þ 2bcy cz syz   c1 7 cx I /1;q ¼ 2bð Þ cx   bc3x  bcx c2y  bcx c2z qx Io c 2    c1 7 cy I cy   bc2x cy  bc3y  bcy c2z qy þ 2b 2 Io c    c1 7 cz I cz   bc2x cz  bc2y cz  bc3z qz þ 2b 2 Io c ð5Þ 53c RT 2ðc1Þ

cous flows upwinding is implemented using Courant split Boltzmann equation as follows:

1 where b ¼ 2RT ; Io ¼ and ca = va  ua where a = x, y, z with ua and va being the fluid and molecular velocity. In the KFVS for vis-

@f @f @f @ v r @f @ v h @f þ vz þ vr þ þ ¼ Jðf ; f Þ @t @z @r @t @ v r @t @ v h @f @f v r @rf f v r v 2h @f v r v h @f þ vz þ þ  ¼ Jðf ; f Þ () þ @t @z r @r r r @v r r @v h ð11Þ

327

A.K. Mahendra et al. / Computers & Fluids 46 (2011) 325–332

Taking W moments with Chapman-Enskog distribution leads to split Navier–Stokes equations as

3.1. Split-stencil least square kinetic upwind method



Let / be any function of x, y and z and further it is assumed that it is given at all nodes i.e. it is given at all points in a cloud. We are interested in finding the derivatives of / at all the nodes. The derivatives are obtained in the following way. Consider a point o surrounded by n points as shown in Fig. 1 then error at any point i in the neighbourhood of o gives us,

@f1 v z þ jv z j @f1 v z  jv z j @f1 v r þ jv r j @rf1 v r  jv r j @rf1 W; þ þ þ þ 2 2 2r 2r @t @z @z @r @r

f1 v r v 2h @f1 v r v h @f1 ¼0  þ W; þ r r @v r r @v h ()

@U @GZ þ @GZ  1 @rGRþ 1 @rGR þ þ þ þ þS¼0 @t r @r r @r @z @z



ð12Þ

where U = [q, quz, qur, quh, qE]T and GZ±, GR± represents h axial and raqu2 dial split fluxes with un-split source term S ¼ 0; 0;  Pr  r h  T shh qur uh srh ; r þ r ; 0 . r 2.3. Treatment of slip flow in the rarefied core The boundary conditions at the surface of the solid object define the distribution function of the reflected particles as a sum of diffuse and specular reflections [31]. With a diffuse reflection or accommodation coefficient r, distribution function can be written as

f1R ðv x ; v y ; v z ; IÞ ¼ f1



v ix ; v iy ; v iz ; I



þ ð1  rÞf1  þ rf0wc v dx ; v dy ; v dz ; I v z <0

v z >0



v sx ; v sy ; v sz ; I

v z >0

ð13Þ

 where f1R ðv x ; v y ; v z ; IÞ is the total, f1 v ix ; v iy ; v iz ; I is the incident,  f1 v sx ; v sy ; v sz ; I is the specularly reflected Chapman-Enskog distribu tion and f0wc v dx ; v dy ; v dz ; I is the diffuse reflected Maxwellian distribution evaluated at the wall conditions denoted by superscript wc such that no particles penetrate the wall. Once the total distribution f1R is obtained the slip velocity is

R R v f R dI d~ v þ N ~ ~ uslip ðux ; uy ; 0Þ ¼ RR RR R1 ~ f dI d v Rþ RN 1

ð14Þ

Ei ¼ D/i  Dxi /xo  Dyi /yo  Dzi /zo

i ¼ 1; . . . ; m

ð15Þ

where Dxi = xi  xo, Dyi = yi  yo, Dzi = zi  zo and D/i = /i  /o. Finding the derivative at point o is a least squares problem where error norm kEk2 is to be minimized with respect to /xo, /yo and /zo using stencil N(Po). N(Po) = {"Pi:d(Po, Pi) < h} is also defined as connectivity for node Po, the term d(Pi, Pj) is the Euclidean distance between nodes Pi and Pj such that it is within the ball of radius h. LSKUM uses the normal equations approach to find Uo ¼ ½/xo ; /yo ; /zo T 2 Rn such that kANUo  D/Nk2 is minimized where data matrix AN 2 Rmn and observation D/N ¼ ½D/1 ; D/2 ; . . . ; D/m T 2 Rm . The normal equations approach uses smaller cross-product matrix C ¼ ATN AN 2 Rnn . The most reliable approach is reduction of matrix AN to various canonical forms via orthogonal transformations [15] and QR approach is one of the way to compute orthonormal basis for a set of vectors. Not all ill effects inherent in the normal equations approach can be avoided by using orthogonal transformation as condition number is still relevant to some extent [14]. Approach of normal equations as well as QR produces inaccurate results when applied to ill-conditioned problem [15]. The cross-product matrix, C in LSKUM becomes ill-conditioned when points in the connectivity lies in a thin disc or a thin cylinder [24]. Thus, we require a formulation which leads to well conditioned matrix C. SLKNS uses a different approach to solve a least squares problem as it generates the non-symmetric cross-product matrix C – ATN AN by suitable selection of sub-stencils such that the matrix is diagonally dominant and well conditioned. The present method like normal equations approach requires (m + n/3)n2 flops as against 2(m  n/3)n2 flops required in QR approach for a full rank least squares problem. For large m  n the SLKNS provides stable results involving about half the arithmetic. 3.2. SLKNS and Selection of sub-stencils

3. Split-stencil least square kinetic upwind method for Navier– Stokes Split-stencil least square kinetic upwind method for Navier–Stokes, SLKNS requires cloud of points and its connectivity. The pre-processing step is similar to LSKUM [24] where the points are generated around each component of the multibody configuration using simple grid generator and then the points around each components are merged to form the cloud of points.

As described earlier, LSKUM follows normal equations approach where kANUo  D/Nk2 is minimized for stencil N(Po). Whereas, SLKNS minimizes kANx Uo  D/Nx k2 ; kANy Uo  D/Ny k2 and kANz Uo  D/Nz k2 with respect to /xo, /yo and /zo, respectively for each carefully selected sub-stencils N x ðPo Þ 2 Rmx ; N y ðP o Þ 2 Rmy ; N z ðPo Þ 2 Rmz and D/Na ¼ ½D/1 ; D/2 ; . . . ; D/ma T 2 Rma where a = x, y and z. For example matrix ANx 2 Rmx n AN contains entries of only those nodes which have Dy and Dz tending to zero such that sub-stencil Nx(Po) N(Po)

Fig. 1. Point connectivity for: (a) LSKUM; and (b) SLKNS with its split six sub-stencils.

328

A.K. Mahendra et al. / Computers & Fluids 46 (2011) 325–332

contains mx nodes. Similarly sub-stencil Ny(Po) N(Po) contains my nodes forms matrix ANy 2 Rmy n AN and sub-stencil Nz(Po) N(Po) contains mz nodes forms matrix ANz 2 Rmz n AN . Minimization of kANa Uo  D/Na k2 with respect to its gradient /ao where a = x, y and z leads to

a

!

Ka ATNa ANa Uo ¼

X a

Ka ATNa D/Na () C Uo ¼ D/Nxyz

ð16Þ

where C ¼ Kx ATNx ANx þ Ky ATNy ANy þ Kz ATNz ANz and D/Nxyz ¼ Kx ATNx D/Nx þ Ky ATNy D/Ny þ Kz ATNz D/Nz with diagonal matrices Kx, Ky and Kz 2 Rnn . The diagonal matrices are subsets of Identity matrix, I such that Kx + Ky + Kz = I with Kx = diag(1, 0, 0), Ky = diag(0, 1, 0) and Kz = diag(0, 0, 1). With the proper choice of sub-stencils Nx, Ny and Nz we can make the matrix C diagonally dominant and off-diagonal terms very small. In SLKNS the basic stencil N(Po) gets divided into  þ six different kind of upwind sub-stencils: N þ x ðP o Þ; N x ðP o Þ; N y ðP o Þ;  þ  N y ðP o Þ; N z ðPo Þ and N z ðP o Þ. For example connectivity parameters  S  yx and Szx define   conical sub-stencils N x ðPo Þ such that Dyi  Dzi  Dxi  ¼ syx i 6 Syx and Dxi  ¼ szxi 6 Szx as follows:

N þx ðPo Þ ¼ f8Pi : Pi 2 NðP o Þ9Dxi 6 0; syx i 6 Syx ; szxi 6 Szx g

ð17Þ

N x ðPo Þ ¼ f8Pi : Pi 2 NðP o Þ9Dxi P 0; syx i 6 Syx ; szxi 6 Szx g

Similarly, connectivity parameters Sxy and Szy define conical sub   Dxi  Dzi  stencils N y ðP o Þ such that Dyi  ¼ sxy i 6 Sxy and Dyi  ¼ szy i 6 Szy . Connectivity parameters Sxz and Syz define conical sub-stencils N   z ðP o Þ     such that DDxzi  ¼ sxzi 6 Sxz and DDyzi  ¼ syz i 6 Syz . Minimization of the i i kANa Uo  D/Na k2 with respect to its gradient /ao where a = x, y and z leads to C Uo ¼ D/Nxyz as follows: P P 3 2 ð Dx2i syx i ÞN ð Dx2i szx i ÞN x x P P 1 6 ð Dx2i ÞN ð Dx2i ÞN 7 3 72 6 3 2 P x ð D/i Dxi ÞNx 6 P 2 P 2 x 7 /xo 6 ð Dyi sxy i ÞN Dyi szy i Þ  7 ð P 7 6 N 6 7 y 76 P ð D/ Dy Þ  7 1 Drxyz 6 P Dy2 y 4 /yo 5 ¼ 6 4 P i i Ny 5 ð Dy2i ÞN 7 7 6 ð i ÞN  y y 7 /zo 6 P ð D/i Dzi ÞNz 7 |fflfflffl{zfflfflffl} 6 ð Dz2 s Þ ðP Dz2 s Þ 5 4 |fflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflffl} i xz i N i yz i N  P 2 z P 2 z Uo 1 D/Nxyz ð Dzi ÞN ð Dzi ÞN z z |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl ffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl ffl} C

ð18Þ P

 P 2  P 2  where Dr xyz ¼ Dx2i N Dyi N Dzi N . The expressions ðÞNx ; x y z ðÞNy and ðÞNz requires values evaluated using split stencils N  x ðP o Þ;  N y ðP o Þ and N z ðP o Þ. When these connectivity parameters Syx ; Szx ; Sxy ; Szy ; Sxz and Syz all tend to zero then the matrix C as shown in





v x þ jv x j @f1 þ W; 2 @x Nþx ðPo Þ



v x  jv x j @f1 v y þ jv y j @f1 þ W; þ W; 2 2 @x Nx ðPo Þ @y Nþy ðPo Þ



v y  jv y j @f1 v z þ jv z j @f1 þ W; þ W; 2 2 @y Ny ðPo Þ @z Nþz ðPo Þ

v z  jv z j @f1 þ W; ¼0 2 @z Nz ðPo Þ

W;

@f1 @t



U kþ1

3 2  k @GX  k @GX þ þ @x @x N  6 x ðP o Þ 7 Nþ x ðP o Þ 7 6   k 7 6  þ k 7 6 þ @GY k @GY þ @y ¼ U  Dt 6 7 @y  Nþ N ðP Þ ðP Þ 6 o o 7 y y 7 6  5 4 @GZ k þ k þ @GZ þ  @z @z N ðP o Þ þ Nz ðPo Þ

ð20Þ

z

where GX±, GY± and GZ± represents the split inviscid and viscous fluxes [18]. Viscous fluxes are upwinded and treated similar to inviscid fluxes. For example viscous part of the mass flux contains terms of shear stress tensor as well as heat flux vector because the Chapman-Enskog polynomials /1,s and /1,q satisfy the Onsager reciprocal relation [32]. The upwind treatment of the viscous term is required for effective capture of cross phenomena of thermal transpiration or thermal creep and the mechanocalorific effect. Second order accuracy in SLKNS is achieved through a two-step procedure similar to one used in LSKUM described in [1], this when coupled with inner iterations leads to second order of accuracy [1].

0.6

Y NACA0012 airfoil at M =0.8 , AOA =10 deg. Re=500

0.5

X Z M 1.08 1.01 0.93 0.86 0.79 0.72 0.64 0.57 0.50 0.43 0.35 0.28 0.21 0.14 0.06

SLKNS solver Fortunato and Magi Catalano et al.

0.4 0.3 0.2 0.1 0 -0.1 -0.2

0

0.25

0.5

0.75

1

X

(a)

ð19Þ

Taking W moment leads to split flux Navier–Stokes state update equations for iteration k + 1 [18,19]

Cf

X

Eq. (18) is always diagonally dominant and well conditioned as off diagonal terms are very small. When all the connectivity parameters tend to zero SLKNS becomes less dissipative as it approaches the Taylor series based finite difference method. SLKNS while retaining the simplicity of LSKUM avoids the ill-conditioning of the matrix C which is the weakness of normal equations approach. The KFVS for Navier–Stokes for 3-D geometries [18] can be derived by using Courant splitting at the Boltzmann level as

(b)

Fig. 2. Transonic flow past NACA0012 airfoil at Mach = 0.8, AOA = 10, Re = 500: (a) Mach contours and flow separation; and (b) skin friction plot.

329

A.K. Mahendra et al. / Computers & Fluids 46 (2011) 325–332

0.2

ber of 500 solved for 10,143 nodes. Fig. 2b shows plot of skin friction coefficient for SLKNS solver and a finite volume laminar viscous solver of Fortunato and Magi [12]. The plot also shows comparison with results of Catalano et al. [7] obtained using the grid size of 257  65 = 16,705 nodes. This paper presents following four test cases to demonstrate the results of SLKNS solver for compressible, rotating viscous flows. Test cases-I,II,III use 3D axi-symmetric formulation of KFVS while test case-IV uses the 3D Cartesian formulation of KFVS.

0.1

4.1. Test Case-I

0.5

0.4

0.3

0

0

0.25

0.5

0.75

1

(R-Ri)/(Ro-Ri) Fig. 3. Comparison of the non-dimensional tangential velocity for SLKNS and DSMC.

4. Results and discussions SLKNS code was validated with non-rotational flow over NACA0012 aerofoil. Fig. 2a shows transonic flow at Mach = 0.8 past a NACA0012 aerofoil at 10 degree angle of attack at Reynolds num-

0.1

0.1

0.09

0.09 Temperature (deg K)

0.08

310.00 308.57 307.14 305.71 304.29 302.86 301.43 300.00 298.57 297.14 295.71 294.29 292.86 291.43 290.00

0.07 0.06 0.05 0.04 0.03 0.02 0.01 0 0.1

0.12

0.14

0.16

0.18

0.08

Temperature (deg K)

310.00 308.57 307.14 305.71 304.29 302.86 301.43 300.00 298.57 297.14 295.71 294.29 292.86 291.43 290.00

0.07 0.06 0.05 0.04 0.03 0.02 0.01 0 0.1

0.2

0.12

0.14

0.16

r

r

(a)

(b)

0.18

0.2

0.00025 0.0002

Kinetic Solver, SLKNS DSMC

0.00015

Axial Mass Flux

z

We first present a test case of rarefied cylindrical Couette flow of argon gas at STP. Inner and outer cylinder have tangential momentum accommodation coefficient, r of 0.1 and radii of 3k and 5k, respectively, where mean free path, k = 6.25  108 m. Inner cylinder rotates with angular speed, x = 5.17  108 rad/s and outer cylinder is held stationary [30]. SLKNS solver is able to predict the anomalous behaviour of increase in gas velocity with the radial distance from inner rotating cylinder. Fig. 3 shows the comparison of the non-dimensional tangential velocity with respect to non-dimensional radial distance using SLKNS and Discrete Simulation Monte Carlo (DSMC) [5]. SLKNS took around 15 times less computational time with fraction of memory requirement as compared to DSMC.

z

Vt/Vsb

SLKNS with σ =0.1 DSMC with σ =0.1

0.0001 5E-05 0 -5E-05 -0.0001 -0.00015 -0.0002 0.1

0.125

0.15

0.175

0.2

Radius

(c) Fig. 4. Comparison of SLKNS and DSMC results: (a) DSMC; (b) SLKNS; and (c) axial mass flux for DSMC and SLKNS solver.

330

A.K. Mahendra et al. / Computers & Fluids 46 (2011) 325–332

0.5

Axial velocity

0.4 0.2

z 0.1

0.3

z

600.83 573.34 545.85 518.36 490.87 463.38 435.90 408.41 380.92 353.43 325.94 298.45 270.96 243.47 215.98

0.15

0.2

0.1

0.05

0 0.75

2.76 2.47 2.18 1.88 1.59 1.30 1.01 0.72 0.43 0.14 -0.16 -0.45 -0.74 -1.03 -1.32

Azimuthal velocity

0.8

0.85

0.9

0.95

0

1

0.7

0.8

0.9

r

r

(a)

(b)

1

2.5 SLKNS solver Analytical (Eigenvalue problem) Radial velocity

1.02 0.77 0.52 0.27 0.03 -0.22 -0.47 -0.72 -0.97 -1.22 -1.47 -1.72 -1.97 -2.22 -2.47

z

0.03

0.02

0.01

0 0.95

Axial Velocity (m/s)

0.04 2

1.5

1

0.5

0.96

0.97

0.98

0.99

1

0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75

r

Axial length, z (m)

(c)

(d)

Fig. 5. Rotating cylidrical annulus with a slowly rotating ring: (a) Cloud of points; (b) recirculation in Stewartson layer; (c) velocity vector in Ekman and Stewartson layers; and (d) axial decay of axial flow at radius = 0.77 m.

4.2. Test Case-II In this test case we present result for argon gas with a stream density of 1.4  1021 mol/m3 in a rotating cylindrical annulus of inner radius 0.1 m and outer radius of 0.2 m and height of 0.1 m with top lid held at 310 K and bottom at 290 K. Fig. 4 shows the comparison between Discrete Simulation Monte Carlo (DSMC) with 100  100 mesh and SLKNS results with cloud of 14,400 nodes. Apart from capturing flow features associated with strong rotation SLKNS solver was also able to capture a weak axial flow of the order 2 ms1 with maximum wall velocity of around 628 ms1. Performance of SLKNS solver depends on speed memory trade-off whereas DSMC’s performance is governed by time-step, number of particles per cell and sampling criteria. Typically, SLKNS takes around 10 to 20 times less computational time with fraction of memory requirement as compared to DSMC. Memory requirement of a DSMC problem is a product of solver specific parameter (memory per particle) and problem specific parameter (particles per cell). In this particular case DSMC solver takes 3.1 K Bytes/particle and the present test case uses 15 particles/cell. This equates to

46.5 K Bytes/cell for DSMC compared to 1.8 K Bytes/node for SLKNS solver. 4.3. Test Case-III We present a third case of argon gas in a cylindrical annulus of inner radius of 0.3 m, outer radius of 1.0 m, height of 1.0 m with a ring of radius 0.01 m located at a radius of 0.85 m and axial location of 0.1 m from bottom. Initial density at the wall being 0.0528 Kg/m3 is held at 300 K and rotates with 100 revolutions per second. Ring which is also held at 300 K rotates at angular speed of 95 revolutions per second. The total number of nodes in the cloud was 15,010. With this highly stretched connectivity at the wall the LSKUM solver with normal equations approach as well as QR approach failed but SLKNS solver on the other hand was able to capture weak secondary flow. Fig. 5a–c shows the cloud of points, streamline plot and vectors in Ekman and Stewartson layers. Variation of axial flow profile can be cast as an eigenvalue problem for unknown decay rates ki. Thus the axial flow variation can be expressed as series of decaying exponentials e.g.

331

A.K. Mahendra et al. / Computers & Fluids 46 (2011) 325–332

Z Z

Y Y

X

O U T F L O W

Mach 3.46 3.29 3.11 2.94 2.77 2.59 2.42 2.25 2.08 1.90 1.73 1.56 1.38 1.21 1.04 0.86 0.69 0.52 0.35 0.17

0.1 0.09 -0.02 0.08

x

0.07

-0.01

X

y

0.06 0.08 0

0.05 y

(a)

(b)

100

Density

SLKNS Solver LSKUM Solver Analytical Rigid Body solution

10-1

10

-2

10

-3

0.05

0.06

0.07

0.08

0.09

0.1

Radius

(c) Fig. 6. Hemisphere facing a rotating flow: (a) Cloud of points generated using hemisphere and cylindrical geometry; (b) Mach contours; and (c) density variation at the bottom lid of the outflow based on SLKNS, LSKUM and analytical rigid body solution.

vz ¼

P

ki z [11,2]. The decay of axial flow at radial location i fi ðrÞe r = 0.77 m predicted by SLKNS solver compares well with the theoretical result as shown in Fig. 5d.

4.4. Test Case-IV We present a fourth test case of an hemisphere of radius 0.005 m placed mid way at a radius of 0.075 m in a annular cylindrical sector of outer radius 0.1 m and inner radius of 0.05 m and height of 0.035 m rotating at 2000 revolutions per second with density of air at wall being 1.271 Kg/m3. The top lid and the rear side of the hemisphere have an outflow boundary as shown in Fig. 6. The total number of nodes in the cloud was 122,921. Result shows a curved shock with strong radially inward flow, the exponential rise in the density compared well with the results of LSKUM solver [19] and analytical result [11] given as h v2  v2 i ð h Þr ð h Þw 2RT . qr ¼ qw e 5. Concluding remarks SLKNS captures weak secondary flow even under strong rotation and is able to simulate the slip flow in the rarefied core. The highly stretched distribution in the boundary region for viscous

problems leads to degeneracy while solving the least squares problem using the normal equations approach. SLKNS resolves this problem by careful selection of sub-stencil by splitting the stencils along the coordinate directions to make the cross-product matrix diagonally dominant. Future scope of work will include detailed research and investigation to further improve this method. Acknowledgements We thank anonymous referees whose remarks helped to improve the paper significantly. The first author is grateful to Professor S.M. Deshpande, Mr. T.K. Bera and Mr. A. Sanyal for encouragement and their support. References [1] Arora K, Rajan NKS, Deshpande SM. On the robustness and accuracy of least squares kinetic upwind method (LSKUM). In: Proceedings of the 12th Asian congress of fluid mechanics, Daejeon, Korea, 2008. [2] Babarsky RJ, Herbst IW, Wood III HG. A new variational approach to gas flow in a rotating system. Phys Fluid 2002;14(10):3624–40. [3] Balakrishnan R. An approach to entropy consistency in second-order hydrodynamic equations. J Fluid Mech 2004;503:201–45. [4] Bark FH, Bark T H. On vertical boundary layers in a rapidly rotating gas. J Fluid Mech 1976;78:815–25. [5] Bird GA. Molecular gas dynamics and the direct simulation of gas flows. Oxford: Clarendon; 1994.

332

A.K. Mahendra et al. / Computers & Fluids 46 (2011) 325–332

[6] Brouwers JJH. On compressible flow in a rotating cylinder. J Eng Math 1978;12:265–85. [7] Catalano LA, De Palma P, Naplitano M, Pascazio G. Genuinely multidimensional upwind methods for accurate and efficient solutions of compressible flows. Notes on Numerical Fluid Mechanism Vieweg-Verlag, Braunschweig 1997;221–50. [8] Chapman S, Cowling TG. Mathematical theory of non-uniform gases. 3rd ed. Cambridge Press; 1970. [9] Chou SY, Baganoff D. Kinetic flux-vector splitting for the Navier–Stokes equations. J Comput Phys 1996;130:217–30. [10] Deshpande SM. A second order accurate, kinetic-theory based method for inviscid compressible flows. NASA Technical Paper No. 2613; 1986. [11] Dickinson GJ, Jones IP. Numerical solutions for the compressible flow in a rapidly rotating cylinder. J Fluid Mech 1981;107:89–107. [12] Fortunato B, Magi V. An implicit Lambda method for 2-D viscous compressible flows. In: Deshpande SM, Desai SS, Narasimha R, editors. Fourteenth international conference on numerical methods in fluid dynamics. Lecture notes in physics. Springer; 1995. p. 259–64. [13] Ghosh AK, Deshpande SM. Least squares kinetic upwind method for inviscid compressible flows. AIAA paper 95-1735; 1995. [14] Golub GH, Saunders MA. Linear least squares and quadratic programming. Technical report no. CS 134, Stanford University; 1969. [15] Golub GH, Van Loan CF. Matrix computations. 3rd ed. Baltimore and London: The John Hopkins University Press; 1996. [16] Ilgaz M, Tuncer IH. A gas-kinetic BGK scheme for parallel solution of 3-D viscous flows on unstructured hybrid grids, AIAA paper 2007-4582; 2007. [17] Johnson EA, Stopford PJ. Shear flow in the presence of ’strong rotation’: II. Approximations for continuum-plus-rarefied flow. J Phys D: Appl Phys 1983;16:1207–15. [18] Mahendra AK. Application of least squares kinetic upwind method to strongly rotating viscous flows. M.Sc.(Engg.) Thesis. Bangalore: Indian Institute of Science; 2003. [19] Mahendra AK, Sanyal A, Gouthaman G. Parallel meshless solver for strongly rotating flows. In: Zeng S, editor. Proceedings of the 9th international

[20] [21]

[22]

[23]

[24]

[25] [26] [27]

[28]

[29] [30] [31] [32]

workshop on separation phenomena in liquids and gases. Beijing: Tsinghua University Press; 2006. p. 127–31. Mandal JC, Deshpande SM. Kinetic flux vector splitting for euler equations. Comput Fluids 1994;23(2):447–78. May G, Srinivasan B, Jameson A. An improved gas-kinetic BGK finite-volume method for three-dimensional transonic flow. J Comput Phys 2007;220:856–78. Mieussens L. Discrete-velocity models and numerical schemes for the Boltzmann–BGK equation in plane and axisymmetric geometries. J Comput Phys 2000;162:429–66. Park JS, Hyun JM. Flow of a compressible uid in a rapidly rotating pipe with azimuthally varying wall thermal condition. J Fluid Mech 2004;518: 125–45. Praveen C, Ghosh AK, Deshpande SM. Positivity preservation, stencil selection and applications of LSKUM to 3-D inviscid flows. Comput Fluids 2009;38(8):1481–94. Prendergast K, Xu K. Numerical hydrodynamics from gas-kinetic theory. J Comput Phys 1993;109:53–66. Sharipov FM, Cumin LMG, Kremer GM. Transport phenomena in rotating rarefied gases. Phys Fluid 2001;13(1):335–46. Sharipov FM, Kremer GM. Non-isothermal Couette flow of a rarefied gas between two rotating cylinders. Eur J Mech B/Fluids 1999;18(1): 121–30. Sugimoto H, Sone Y. Numerical analysis of steady flows of a gas evaporating from its cylindrical condensed phase on the basis of kinetic theory. Phys Fluid A 1992;4(2):419–40. Taheri P, Struchtrup H. Effects of rarefaction in microflows between coaxial cylinders. Phys Rev E 2009;80:066317-1–066317-16. Tibbs KW, Baras F, Garcia AL. Anomalous flow profile due to the curvature effect on slip length. Phys Rev E 1997;56(2):2282–3. Xu K, Li Z. Microchannel flow in the slip regime: gas-kinetic BGK-Burnett solutions. J Fluid Mech 2004;513:87–110. Zhdanov VM, Roldughin VI. Non-equilibrium thermodynamics and kinetic theory of rarefied gases. Physics-Uspekhi 1998;41(4):349–78.