Lattice_Boltzmann方法与粘滞流模拟解

文章编号:1000-2634(2003)02-0080-03

Lattice2Boltzmann方法与粘滞流模拟解Ξ

陈代　

(西南石油学院工业自动化系,四川南充637001)

摘要:应用Lattice2Boltzmann方法,由72Bit的FHP模型得到了描述粘滞流的N2S方程。在小Knudsen数条件下,给出了Poiseuille管流的密度分布函数和模拟解,确定了弛豫时间τ,模拟解与其解析解吻合。

关键词:Lattice2Boltzmann方法;N2S方程;密度分布函数;Poiseuille流

中图分类号:O37 文献标识码:A

引　言

Lattice2Boltzmann方法(LBM)是可用于模拟粘

滞流的一种数值计算方法,它由Lattice2G as方法

(L GM)发展而来[1]。在LBM中,流体场离散为由

格点组成,格点处粒子运动状态是由分布函数描述

的。由于粒子的运动和碰撞,其分布函数在离散的

时间步上也从一个格点移动到另一格点,且弛豫地

趋于局域平衡,在该过程中遵守质量和动量守恒。

在粒子的迁移和碰撞过程中,分布函数的变化服从

Boltzmann演化方程[2],这样可将流体运动的力学过

程转化为关于分布函数随时间变化的运动方程。因

此,LBM充分体现了宏观流体力学过程是微观运动

的统计平均。LBM作为一种全离散和并行运算的

局部动力学模型,现已被应用在流体、传热等领

域[3]。本文由72Bit的FHP模型运用LBM得到描

述粘滞流的N2S方程,给出了Poiseuille管流的LBM

模拟解。

1　LBM与描述粘滞流的N2S方程

本文采用二维72Bit的FHP(正六边形格子)为

LBM的模型。在每一格点处将速度空间离散为由

一个速率为零(静止)和六个速率为c的速度Cα组

成,速度Cα为

Cα=c[cos 2π(α-1)

6x1

+

sin2

π(α-1)

6x2]

(α=1,2, (6)

Cα=0(1)

这样,流体流动的空间也相应被离散为由正六

边形网格组成(见图1),粒子从一格点到另一格点

的离散时间步长为δt

。

图1　72Bit的FHP模型网格

在任一格点处,流体的质量密度ρ和宏观速度

u定义为

ρ(x,t)=Σ6

α=0

fα(x,t)

u i(x,t)=

1

ρ

Σ6

α=1

cαi fα(x,t)(2)

式中的fα(x,t)表示在t时刻、x处α方向具有

单位质量的粒子密度分布函数,在t+δt时刻粒子

移动到空间x+Cαδt处,其分布函数为fα(x+

Cαδt,t+δt)。考虑到在流动过程中除了粒子的迁

移和碰撞,还有外压力场存在,对于低速的粘滞流,

外压力克服粘滞力,描述fα(x+Cαδt,t+δt)和

fα(x,t)间演化关系的Boltsmann方程[2]可表为

第25卷　第2期 西南石油学院学报 Vol.25　No.2　2003年　4月 Journal of S outhwest Petroleum Institute Apr　2003　Ξ收稿日期:2002-06-14

基金项目:“油气藏地质及开发工程”国家重点实验室开放基金(PLN9733)。

作者简介:陈代　(1949-),女(汉族),重庆长寿人,副教授,理学硕士,主要从事物理教学、渗流物理机理及计算物理研究。

f α(x +C αδt ,t +δt )-f α(x ,t )=

Ωαc +B α

δt (3)其中碰撞项Ωαc =(5t f α)c δt 、

外压力项B αδt 分别表示在δt 时间内碰撞、外压力场的存在对α方向粒子密度分布变化的影响。

对f α(x +C αδt ,t +δt )作泰勒级数展开有f α(x +C αδt ,t +δt )=f α(x ,t )+δt [5t f α+

c αi 5i f α]+(δt )2[12

52

t f α+c αi 5t 5i f α+ 12

c αi c αj 52ij f α]+O (δt 3

max )(4)对(4)式取到2阶,考虑低速缓变的情况,忽略(4)式右端中含对时间的求导项,代入(3)式,两边乘以C α,并对粒子运动的所有α方向求和。因为碰撞

过程不引起系统的总动量改变,有∑6

α=1

C αΩα=0。设

外压力场沿x i 方向,压力梯度为5i p ,由在时间δt

内外压力作用对动量增量δ(ρu )i 的关系有

-5i p δt =δ(ρu )i =δt ∑6

α

=0c αi B α

得外压力项B α表示为[4]

B α=0

α=0-1

4c 5i p

α=1,2,61

4c 5i

p

α=3,4,5 (5)

由此可得用张量表述的方程为

5t ρu i +5j

Πij +1252jk T ijk

=-5i p (6)

其中(6)式左端的有关项为

Σ6

α=0c αi f α=ρu i Σ6α=0

c αi c αj f α=Πij

δt Σ6α

=0

c αi c αj c αk f α=T ijk

要由(6)式得到流体运动方程,必须求得f α(x ,t )。由时间弛豫法[2]知(3)式中的(5t f α)c 为

(5t f α)

c

=

f α-f (0)

α

τ

(7)

其中τ表示弛豫时间,f (0)

α

(x ,t )表示在(x ,t )处α方向的局域平衡分布函数。在低Much 数,即

|u |值很小时,由费米分布f (0)α可表为

[3,5]

f (0)

α

=ρ[ω0+ω1(C α?u )+

ω20(C α?u )(C α?u )+ω21(u ?u )](α=1,2, (6)

f (0)

α=ρ(z 0+z 21|u |2) (α=0)

其中ω0、ω1、ω20、ω21、z 0和z 21为待定参数。由非平衡统计力学中Chapman 2Enskog 截断[6]

可得f α与f (0)

α的关系为

f α=f (0)α-τc αi 5i f (0)

α+τ

B α(8)

由正六边形格子的对称性、对流项满足的伽利略不变式条件、以及为保证不出现非物理压力因素,要消除粒子密度对宏观流速u 2的依赖等条件,可求

得局域平衡分布函数f (0)

α的各参数为

ω0=z 0=17 ω1=13c

2

ω20=13c 4 ω21=16c 2

z 21=-1

c

2

这样f (0)

α完全确定,由(8)式f α可确定,且(6)

式表为

5t ρu i +5j ρu i u j =-5i ρ+

ξ5i 5j ρu j +η5j 5j ρu i

(9)

(9)式是描述粘滞流的N 2S 方程,其中ξ=

c 2

2

(τ

-δt 2),η=c 2

4(τ-δt 2)。由此知,由已知的流体宏观量—流体密度ρ、流速u 的分布和外压力场,利用

(2)、(3)、(8)式可确定下一时刻流体宏观量。以下针

对Poiseuille 流,给出其LBM 解。

2　由LBM 得到的Poiseuille 管流解

本文所谈的毛细管流满足Knudsen 数K n <0.01条件,管壁处u =0。在圆柱坐标中设外压力场

和流体的流向均沿圆柱毛管的轴向z ,毛管半径为R ,由(9)式,其Poiseuille 流满足的方程为ηρ1r 55r (r 5u 5r )-5p 5z

=0

(11)

(11)式的解析解为

u =u 0(1-r 2

R

2)(12)

其中u 0=-R

2

4ηρ5p

z

。由(8)式以及Poiseuille 管

流的特点,知分布函数f α是r 的函数,表述为

f 0(r )=ρ[

17-1c 2u 2

]f α(r )=ρ[17±13c u +1

2c

2u 2]

τ14c 5p 5z

(α=1,4) 1

8第2期 陈代　:　Lattice 2Boltzmann 方法与粘滞流模拟解　

f α(r )=ρ[17+1

6c u ]±τ3ρu 06R 2

r +

τ14c 5p 5z

(α=2,6) f 3(r )=ρ[

17-1

6c u ] τ3ρu 06R 2

r -

τ14c 5p 5z

(α=3,5) (13)由η的表达式知弛豫时间τ应大于δt/

2,以下在满足原边界条件下,采用在毛管外虚拟反向流场的方法确定τ。

图2　毛管边界示意图

原管壁位置用w 标示,A 是毛管内临近管壁w

的格点(见图2)。

设在毛管流场外存在有一反向的压力场,其压力梯度与原毛管内的大小相等,其中的流体流速u ′与原毛管内的流速u 等大反向,仍为Poiseuille 流,即流速分布呈抛物线状,由此知在原管壁位置w 处有u |±R =0,满足原题的边界条件。设

A ′是与A 相邻的w 外的格点之一,它们的坐标分别

为r A =R -3c δt/4和r A ′=R +3c δt/4,在A ′和

A 格点的分布函数间满足Boltsmann 演化方程(3)

f 2(A ′)-f 2(A )=-δt

τ[f 2(A )- f (0)

2(A )]+B 2δt (14)对于Poiseuille 流,在A 、A ′两点的位置处流速

大小相等,流动方向相反,即u A =-u A ′,由(2)式知

应有f 5(A )-f 2(A ′),因此由(13)、(14)两式及η可解得

8τ2-12τδt +3(δt )2

=0

(15)

关于弛豫时间τ的解为

τ=3±34

δt =

3+3

4δt 3-3

4

δt 舍去不适的解τ=

3-34δt ,得τ=3+3

4

δt ≈1.175δt 。δt 代表了粒子在格点间的迁移时间,即

粒子间发生一次碰撞的时间,计算得到的τ与

δt 的

大小关系与要达到局域平衡的条件是吻合的[6]。由

LBM 算出的任一时刻的速率u 由(2)和(13)式给出。图3给出了(12)式表示的Poiseuille 流解析解曲线与LBM 计算结果的比较,二者完全相符。

图3　u/u 0-r 2/R 2曲线

3　结　论

本文应用Lattice 2Boltzmann 方法,由72Bit 的FHP 模型得到了描述粘滞流的N 2S 方程。在小Knudsen 数条件下,给出了Poiseuille 管流的密度分

布函数f α(x ,t )和模拟解,确定了弛豫时间τ,模拟

解与其解析解吻合。

参考文献:

[1]　Chen H ,Chen S ,Matthaeus W H.Recovery of the

Navier 2Stokes equations using a lattice 2gas Boltzmann method[J ].Physical review A ,1992,8(45):R5339-5342.

[2]　王竹溪.统计物理学导轮[M ].北京:人民教育出版社,

1965.

[3]　Benzi R ,Succi S ,Vergassola M.The Lattice Boltzmann

equation :theory and applications [J ].North 2Holland ,Physics Reports ,1992,222(3):145-197.

[4]　Noble D R ,Chen S ,G eorgiadis J G ,et al.A consistent

hydrodynamic condition for the lattice Boltzmann method [J ].Phys Fliuds ,1995,7(1):203-209.

[5]　李元香,康立山,陈毓屏.格子气自动机[M ].北京:清

华大学出版社,广西:科学技术出版社,1994.

[6]　李如生.平衡与非平衡统计力学[M ].北京:清华大学

出版社,1995.

(编辑　罗先碧)

28西南石油学院学报 2003年

and the strain factor k2is reasonably constant and is not affected by the specimen ligament thickness and loading rate,i.e.,n k is considered as a material constant affecting the ductile crack propagation resistance.According to the Necking Z one Model, the increase in fracture resistance with res pect to loading rate is mainly attributed to the strainΟrate strengthening effect of the material.The Necking Z one Model not only suggests a simple engineering approach of calculating the energy dissipated by the ductile fracture during impact case,but also suggests a simple experimental method of measuring dynamic crack resistance us2 ing dynamic tensile test and static tensile tearing test.

K ey w ords:ductile fracture;crack propagation;necking zone model;mild steel;double cantilever beam s pecimen;im2 pact loading

ANA LYSIS OF THE FACT ORS OF THE CUTTING EFFI2 CIENC Y AB OUT SING LEΟCONE BIT

ZHAO Y an(Scientific and Technology University of Huazhong,Wuhan Hubei430074,China),ZHAN G FuΟrun, GUO Wei.JOU RNAL O F S OU T HW ES T PET ROL EU M IN2 S TI TU T E,V OL.25,N O.2,69-71,2003(IS S N1000-2634,IN CHIN ES E)

Through analyzing the cutting mechanism of the single cone bit,the bottomhole coverage and cutting uniformity of the compacts were researched from three as pects:compact expo2 sure,shape and setting pattern.The cutting structure for the singleΟcone bit was designed,the compact exposure should be higher than the height of rock ridges at bottom hole,phase dif2 ference of the compact in the longitude and latitude directions of the spherical cone should be optimized,and when fitted com2 pacts onto the cone,compacts of different shapes which suit the movement features of s pecific areas on the cone should be select2 ed.

K ey w ords:singleΟcone bit;cutting structure;compact ex2 posure;compact shape

DECISIONΟMAKING ANA LYSIS ON THE SCHEME T O DE2 SIGN AN D MANUFACTURE THE HERRINGB ONE GEAR WHEE L ON FΟ1300MU D PUMP

L I Qin(S outhwest Petroleum Institute,Nanchong Sichuan 637001,China),HUAN G ZhiΟqiang,L IU DeΟming,et al. JOU RNAL O F S OU T HW ES T PET ROL EU M IN S TI TU T E, V OL.25,N O.2,72-75,2003(IS S N1000-2634,IN CHI2 N ES E)

On the basis of analyzing and comparing three kind of schemes including existing runout slot,no runout slot and spliced assembly to design and manufacture the herringbone gear of FΟ1300mud pump which have three carboys and single ef2fect,the three kind of schemes were analyzed by evaluating the strength of gear,the woke life and the cost.The conclusion in2 dicates that the scheme is the excellent to manufacture the her2 ringbone gear which has no runout slot by shaving gear,and the indicator of the scheme was the highest when the ratio of work life of gear and the cost price to manufacture gear was regarded as the weighing indicator.Finally,considering three factors of other of influencing scheme evaluation,by the method of weighted scoring in the value engineering theory,it has been carried on the synthetical analysis,and the mentioned above conclusion has been confirmed by the result further.

K ey w ords:drilling pump;gear;scheme;analysis

CONSTRUCTION SCHEME OF DAZHANGTU O UN DER2 GR OUN D G AS ST ORAGE

CU I LiΟhong(Dagang Oil Field of CNPC,Beijing101149, China),SHU ZhuangΟzhi,Y AN G ShuΟhe,et al.JOU RNAL O F S OU T HW ES T PET ROL EU M IN S TI TU T E,V OL.25, N O.2,76-79,2003(IS S N1000-2634,IN CHIN ES E) Based on the gas reservoir geology comprehensive study, for the design project of Dazhangtuo Underground G as Storage as an example,many engineering methods were discussed by the different reservoirs,and the designed gas storage running basic index according to the using gas regulation of Beijing city, and at the same time making use of the test well and the gas well nodal analysis to determine the gas well proper productivi2 ty,it will provide the theoretical evidence for gas storage project design,it has proved that the scheme of the gas storage was true.

K ey w ords:gas storage;proper productivity;running pro2 ject;storage volume

SOL UTIONS OF LATTICEΟB OLTZMANN METH OD AN D VISCOSIT Y F LOW

CHEN DaiΟxun(S outhwest Petroleum Institute,Nanchong Sichuan637001,China).JOU RNAL O F S OU T HW ES T PET ROL EU M IN S TI TU T E,V OL.25,N O.2,80-82,2003 (IS S N1000-2634,IN CHIN ES E)

NavierΟStokes(NS)equations of viscosity flow were ex2 pressed with LatticeΟBoltzmann method in low Mach number range and a7BitΟFHP model.The population distribution,the function of the density distribution of Poiseuille pipeline flow and the simulating solution was presented.the relaxation timeτwas determined under the small Knudsen number.The solution of the simulation was in agreement with the analytical solution.

K ey w ords:LatticeΟBoltzmann Method;NavierΟStokes e2 quations;density distribution;poiseuille flow

ⅥAbstract 2003