二阶精确格子Boltzmann非牛顿流体模型

二阶精确格子Boltzmann非牛顿流体模型JBoyd1,2,JBuick1,2andSGreen21PhysicsandElectronics,UniversityofNewEngland,Armidale,NSW,2351,Australia2CardiovascularResearchGroup,UniversityofNewEngland,Armidale,NSW,2351,Australia摘要：二阶精确的格子Boltzmann模型，提出了非牛顿流体。非牛顿流动性是使用幕法则模型实现的。它可以估算出模型的精确程度，同时不会限制这个模型。二阶精度由剪切变稀和剪切增稠液体的幕法则模型参数范围给出。这些结果与Gabbanelli等人的结果相比，精确度更高，并且得到了更快的计算效率。结果表明了格子Boltzmann方法适用于非牛顿流体模拟。导言对非牛顿流体及其性质的研究在很多领域都有应用，包括地球物理(Ashby和Verrall1977),水文学(Federico1998),材料科学(Bird和West1995)以及生理学建模(Nichols和O'Rourke2005)等等。在很多情况下，由于流体复杂的几何性质以及非牛顿性，分析方法并不适用。因此，数值模拟成为一种有效的方法。本文我们将考虑用格子Boltzmann方法对非牛顿流体进行模拟。格子Boltzmann方法使用了一个简化的动力学方程，是一个有二阶精度的流体模拟方法。它已经应用于很多一般性的问题上，包括湍流(Cosgrove等2003),磁流体力学(Chen等1991)，多孔介质流量(Manz等1999)，多相流(Swift等1996)和血流量(Fang等2002,Tamagawa和matsuo2004,Artoli等2004,Boyd等2005,Yi等2005)，非牛顿流动(Gabbanelli等2006)等问题。其基础格子结构和局部计算性使其非常适用于并行实现(Kandhai等1998)。本文用格子Boltzmann方法模拟了非牛顿流体流动模型，并且在二维刚性管流中证明了其精确性。用幕法则模型来描述非牛顿流体的性质。这个模型可以比较结果与解析解的差别，同时不会限制这个技术的使用。我们发现，格子Boltzmann方法保持了非牛顿流的二阶精确性，表明其适用于对剪切依赖非牛顿流的模拟，包括那些涉及较复杂结构的流体。理论背景2.1格子Boltzmann方法格子Boltzmann方法(Chen和Doolen1998)最近已经发展成为对流体流动模拟的一种替代方法。在格子Boltzmann方法粒子分布函数中，f.(x，t)在点x时间t处只限于i在正则格上同步移动。分布函数在保持质量，动量格，各项同性和伽利略不变性的同时，对格子有所影响。在这里，i表示格子是否与分布函数相链接。在本文中使用的格子是D2Q9，如图1所示。图1.D2Q9格子。黑圆为结点，线为连接方向，编号从1到811)ii格子上分布函数的变化是由离散Boltzmann方程表示的(Chen和Doolen1998):f(x+eAx,t+At)二f(x,t)+0(x,t),(i二0,1,...M)如格子D2Q9，见图1。其中e=(0,0),(i=0),兀兀22e=(cos((i-1)),sm(i-1)),(i=1,2,3,4),22i22，0，0为碰撞算子。在每个节点处，由P=工f和Pu=工fe(3)，流体密度P和速iiiiii度u可以直接由分布函数计算得到。假设分布函数/可以围绕局部均衡分布正式展开，使得f-feq+£fneq(4)，其

iiii中，e是一个很小的参数，通常取Knudson数，feq为一个均衡分布函数，fneq为非均ii衡分布函数，选择feq使其满足P=工fe且Pu=工f妣(5)，并且假设非均衡分iiiiii布函数fne河以进一步展开为fneq二f⑴+f⑵+o(82)(6)，其中(k(k)e二0,k=1,2(7)。碰撞算子0由Bhatnagar-Gross-Krook近似给出(Bhatnagar等1954,Chen和Doolen-11998)，即0二[f(x,t)—feq(x,t)](8)，其中T为松弛时间。D2Q9格子在二维中TOC\o"1-5"\h\ziTii93分布函数的均衡形式由式feq(x,t)二wP(1+3e-u+(e-u)2—u2)(9)给出。其iii2i2中当i=1，2，3，4时w二4/9,w二1/9，i=5，6，7，8时，w二1/36，松弛时间T0ii2T-1与运动粘度v有关，且满足v—(10)。格子Boltzmann方法在几乎不可压缩极限条件下与NavierStokes方程相同，并且在流体的内部满足二阶精确性(Chen和Dollen1998)。压强为p的不可压缩流体的应力张量由式◎°=-p3+2qS(11)给出。其中5apapapap1为Kronecker函数且为应变率张量并满足S-(▽u+▽u)(12)。ap2paap可以证明在格子Boltzmann方法中S可以在每个节点的局部区域中计算得出，即ap3ap2TS=-Yf⑴ee(ap2T在格子Boltzmann方法算法中，f(1)项通常作为速度计算的一部分计算。因此由于i这种方法的计算剪切无需计算速度衍生，它是有效的。此外，在局部计算剪切时，如果并行使用格子Boltzmann方法将会非常便捷。2.2幂法则模型在随后的讨论中，我们定义应变率张量的第二个变量为D=£SS(14)IIapapa,P=1其中在二维情况下l=2。剪切速率定义为*=2\气(15)。幂法则模型是对非牛顿流体最简单的概括之一。在此模型中，表观粘度由v(Y)=m?|n一1(16)给出(Quarteroni等2000)，其中m和n通常是由拟合方程(16)n亠n亠1)n](17)，G丄nLn£1的解析解，如下(Robson2003)：u(y)=()n()[()n2mn+12其中L为管径，G=-为驱动流动的压力梯度。dx参数n的值决定了在剪切速率变化时流体的响应，当n<1时，流体为剪切变稀,当n=1时，流体为牛顿流体，当n>1时，流体为剪切增稠。我们注意到，对剪切变稀流体(即n<1)，有limm沖-1=。同时m的量u2-nLn=—0—u2-nLn=—0—m方程(16)可以无量纲地得到下列类似于Reynold数的无量纲数Re曲(18)，其中m和n为幕法则参数，u°为在管宽L中的最大流量。方法和结果通过联结方程(10)和方程(16)，在格子Boltzmann方法中使用了幕法则非牛顿流，从而在每个节点处给出了依赖于剪切的松弛时间工。通过方程(13)可以计算出剪应力，从而计算出剪应变。为了测试在非牛顿情形下格子Boltzmann方法的精确性，在RepL=100，n分别为0.25,0.5,0.75时进行模拟。子网格的二阶精确边界是用来实现管道的几何形状。为了避免不可压缩，应确保马赫数小于0.03，此时模拟中的G，L和m可以取很多值。当满足如下的标准时，模拟才可进行，即满足£||u(x,t)-u(x,t-1)||<£(19)，x其中£很小，可取其为£=1X10-10。得到的结果与方程(17)给出的解析解相比较，计算出全局误差(20),其中下标a指精确解析解，下标b指Boltzmann模拟值。全局误差结果见图值。全局误差结果见图2。1^05.0」0.0Jg0.001O.CXHJPipewidth(Latticeunits)R«R«R«1.0图2.幂法则模型流的全局误差1^05.0」0.0Jg0.001O.CXHJPipewidth(Latticeunits)R«R«R«1.0图2.幂法则模型流的全局误差-J-0.75-0.5-0.2505JNormaliseilradialposition图4.正常幂法则模型流的流速剖面。粗黑线代表标准牛顿流的抛物线分布图2中的黑线代表了斜率为-2的线，表明了二阶性。可以看出，不同n参数代表的数据与线的斜率相符合。当n较小时，误差较大。图3(a)-(d)描述了在(a)n=0.25,(b)n=0.50,(c)n=0.75和(d)n=1.25情形下标准化流速剖面分析结果(实线)与格子Boltzmann方法所得结果(圈)的比较。可以看出，格子Boltzmann方法精确地模拟出剪切变稀(n〈l)和剪切增稠(n〉l)流体的流速剖面。当n〈l时，我们看到流速剖面整体比较平，且n越小,流速剖面越平。相反，当n=1.25时，在中间剖面的高峰值处有较大的曲率。在图4中可以更明显地看出这点，它描述了当n不同时，正常流体剖面间的比较。讨论图2表明在二维幕法则流通过一个刚性管的情形下，格子Boltzmann方法保持了二阶精确性。这种方法是对Gabbanelli等方法的改进，他们用一阶差分方法(Aharanov和Rothman1993)来估计剪切，只取得了一阶精确地结果。这里的结果提升到二阶的精确度，并且加快了算法的计算效率。在图3中，对所有的n,模拟得到的流速剖面与解析方法得到的剖面基本相符。当n=0.25时，误差最大，此时格子Boltzmann方法与解析方法相比，在中间高峰值处不够平。结论提出了对剪切依赖非牛顿流体的二阶精确格子Boltzmann方法。这种方法避免了计算剪切时耗时推到导速度数据。结果表明格子Boltzmann方法对剪切依赖的非牛顿流体模拟是适用的，包括具有更为复杂结构方程的流体。致谢这份工作由SigmaXiGrantno10040015和theAustralianPostgraduateAward部分支持，在此深表感谢。参考书目AharanovEandRothmanDH1993Non-Newtonianflow(throughporousmedia):alattice-BoltzmannmethodGeophys.Res.Lett.20679-82ArtoliA2003MesoscopicComputationalHaemodynamics(Wageningen:Ponsen&Looijen)ArtoliAM,KandhaiD,HoefslootHCJ,HoekstraAGandSloopPMA2004LatticeBGKsimulationsofflowinasymmetriebifurcationFutureGener.Comput.Syst20909-16AshbyMFandVerrallRA1977Michromechanismsofflowandfracture,andtheirrelevancetotherheologyoftheuppermantlePhil.Trans.R.Soc.A28859-95BhatnagarPL,GrossEPandKrookM1954Amodelforcollisionprocessesingases:I.Smallamplitudeprocessesinchargedandneutralone-componentsystemPhys.Rev.94511-25BirdRBandWestJM1995ConstitutiveequationsforpolymericliquiAdsnnu.Rev.Fluid.Mech.27169-93BoydJ,BuickJ,CosgroveJAandStansellP2005ApplicationofthelatticeBoltzmannmodeltosimulatedstenosisgrowthinatwo-dimensionalcarotidarteryPhys.Med.Biol504783-96ChenS,ChenH,MartinezDandMatthauesW1991LatticeBoltzmannmodelforsimulationofmagnetohydrodynamicsPhys.Rev.Lett.673776-9ChenSandDoolenGD1998LatticeBoltzmannmethodforfluidflowsAnnu.Rev.FluidMech.30329-64CosgroveJA,BuickJM,TongeSJ,MunroCG,GreatedCAandCampbellDM2003ApplicationofthelatticeBoltzmannmethodtotransitioninoscillatorychannelflowJ.Phys.A:Math.Gen.362609-20FangHP,WangZW,LinZFandLiuMR2002LatticeBoltzmannmethodforsimulatingtheviscousflowinlargedistensiblebloodvesselsPhys.Rev.E65051925FedericoVD1998Non-NewtonianFlowinavariableaperturefractureTransportPorousMed.3075-86GabbanelliSG,DrazerGandKoplikJ2006LatticeBoltzmannmethodfornon-NewtonianfluidflowsPhys.Rev.E72046312GuoZ,ZhengCandShiB2002AnextrapolationmethodforboundaryconditionsinlatticeBoltzmannmethodPhys.Fluids142007-10KandhaiD,KoponenA,KoekstraAG,KatajaM,TimonenJandSlootPMA1998Lattice-BoltzmannhydrodynamicsonparallelsystemsComput.Phys.Commun.11114-26ManzB,GladdenLFandWarrenPB1999Flowanddispersioninporousmedia:latice-BoltzmannandNMRstudiesAIChEJ.451845-54NicholsWMandO'RourkeMF2005McDonald'sBloodFlowinArteriesThetoirceal,ExperimentalandClinicalPrinciples5thedn(NewYork:OxfordUniversityPress)QuarteroniA,TuveriMandVeneziani