顶盖驱动流的格子Boltzmann模拟代码
基于格子Boltzmann方法的气动声学计算
基于格子Boltzmann方法的气动声学计算司海青;石岩;王兵;吴晓军【摘要】研究了应用于气动声学计算中的格子Boltzmann方法(Lattice Boltzmann method,LBM),采用非平衡外推格式处理壁面边界条件,远场无反射边界条件采用吸收层边界条件.首先,将LBM应用于顶盖驱动的空腔流动模拟进行程序验证;然后,数值研究气动声学中的几个典型问题,特别讨论了黏性对LBM数值解的影响,并与传统的四阶精度低色散保持格式(Dispersion-relation-preserving,DRP)比较,检验了该方法模拟气动声学基本问题的能力,为进一步运用该方法模拟复杂物体产生噪声奠定基础.研究表明,尽管标准LBM方法在时间、空间上仅有二阶精度,但是,LBM方法计算得到的结果能和分析解保持一致,它是有效且可行的气动声学计算方法.【期刊名称】《南京航空航天大学学报》【年(卷),期】2013(045)005【总页数】5页(P616-620)【关键词】格子Boltzmann方法;壁面边界条件;吸收层边界条件;计算气动声学;计算流体力学【作者】司海青;石岩;王兵;吴晓军【作者单位】南京航空航天大学民航/飞行学院,南京,210016;南京航空航天大学民航/飞行学院,南京,210016;南京航空航天大学民航/飞行学院,南京,210016;中国空气动力研究与发展中心,绵阳,621000【正文语种】中文【中图分类】V211.3格子Boltzmann方法（Lattice Boltzmann method，LBM）的产生与发展，不仅在计算流体力学领域中产生了深远的影响，它所使用的处理方法和观点对其他学科也是富有启发性的。
尽管LBM方法［1］在计算流体力学领域中已得到较多应用，但它在计算气动声学领域中的研究［2－5］与应用相对较晚。
近几年，在计算动力学（Computational aeroacoustics，CAA）研究领域，LBM正逐渐受到国外研究者们的足够重视，Buick等［5］运用LBM解决无黏性的声传播问题，然后，Dellar等［6］运用该方法求解含有黏性的声传播问题。
用于界面对流模拟的格子Boltzmann方法
摘 要 ：通过 引入外力修正 系数和界 面对流强度 ，在介观尺度上运用格 子Ｂ ｌｍａ ｎ ｏｔ ｎ 方法模拟 乙酸 乙酯 吸收二 ｚ
氧化碳传质过程 的界ＮＲ ｙｅ ｈ ａｌｉ 对流现 象。基 于实验 测定 的流 场平均速 度值 来确 定外力修 正系数，在外力修 正 ｇ 系数确定 的条件下 ， 进行 了平均速度分 析 ，其结果 与实验相符 ；并且考查 了平均浓度和界面对流强度随 时间变
Ｆ Ｕ Ｂｏ， ＹＵ ＡＮ ｉ ａ ｇ， Ｃｈ ｎＷ ｅ ， Ｃｈ ｎ Ｓ ｕ ｎｇ Ｘ ｇｎ ｅ ｉ ｅ ｈ ｙｏ
（ｔｔＫ ｙＬ ｂ ｒｔ ｙｏ Ｃ ｅ ｉ ｌ ｎ ｉｅｒ ｇ Ｔ ｎｉ Ｕ ｉ ｒｉ， ｉ ｊｎ３ ０ ７， ｈ ａ Ｓ ｅ ｅ ａ ｏａ ｒ ｈ ｍ ｃ ｇｎｅｉ ， ｉ ｊ ｎｖ ｓｙ Ｔａ ｉ ００ ２ ｉ ） ａ ｏ ｆ ａＥ ｎ ａ ｎ ｅ ｔ ｎ Ｃ ｎ
化关系。
关键 词 ：化学工程 ；外力修 正系数 ；格子Ｂ ｌ ｍａ ｎ ｏｔ ｎ 方法 ；界面对流强度 ｚ 中图分类号 ：Ｔ ２ ． Ｑ０ １ ４ 文献标识码 ：Ａ 文章编号 ：１７ —７ ８ （０ ８ １ —０ ４ —６ ６ ３ １０２ ０ ）２ ９ １
Ａ ｔｉ ｅ Ｂｏ ｔ ｍ ａ ｎ ｍ ｅ ｈｏ ｏ ｉ ｕｌ ｔｏ Ｌａ ｔｃ — ｌｚ ｎ ｔ ｄ ｆ ｒｓｍ ａ ｉ ｎ ｏ ｔ ｒ ａ ｉ ｌ ｏ ｖ ｃ ｉ ｎ ｐｈ ｎ ｍ ｅ ｏ ｆ ｎ ｅ ｆ ｃ ａ ｎ ｅ ｔｏ ｅ ｏ ｉ ｃ ｎ ｎ
Ａｂｓｒ ｃ ：Ｂｙ ｉ ｔｏ ｕ ｉｇ ｃ ｒｅａｉ ｎ ｃ ｅ ｃｅ ｔｏ ｘ ｅ ｎ ｌｆ ｒ ｅ ａ ｄ ｉｔ ｒａ ｉ ｌｃ ｎ ｅ ｔｏ ｎ ｅ ｓｔ ，ａ １ｔ ｃ ｔａ ｔ ｎ ｒ ｄ ｃｎ ｏ ｒ ｌｔｏ ｏ ｆ ｉｎ ｆ ｅ ｔ ｒ ａ ｏ ｃ ｎ ｎ ｅ ｆ ｃａ ｏ ｖ ｃｉｎ ｉ ｔｎ ｉ ｉ ｙ ａ ｔ ｅ ｉ Ｂｏ ｔｍａ ｎ ｍ ｅｈ ｄ ｉ ｓ ｄ ｆ ｒｓｍ ｕ ａｉ ｎ ｏ ｈ ｎ ｅ ｆ ｃａ ｙ ｅｇ ｏ ｖ ｃｉ ｎ ｐ ｅ ｏ ｅ ａ ｉ ｈｅ ｍ ａ ｓｔａ ｓｅ ｌ ｚ ｎ ｔ ｏ Ｓｕ ｅ ｏ ｉ ｌｔｏ ｆｔ ｅ ｉ ｔｒａ ｉｌＲａ ｌｉ ｈ ｃ ｎ ｅ ｔｏ ｈ ｎ ｍ ｎ ｎ ｔ ｓ ｒ ｎ ｆ ｒ
格子Boltzmann方法模拟二维轴对称狭窄血管内的脉动流
201020446(2) 北京师范大学学报(自然科学版)Journal of Beijing Normal University (Natural Science ) 139 格子Boltzmann 方法模拟二维轴对称狭窄血管内的脉动流3张立换 康秀英 吉驭嫔(北京师范大学物理学系,100875,北京)摘要 将格子Boltzmann 方法应用到二维轴对称余弦狭窄血管模型,模拟比较加入脉动后流场速度、压强和剪切应力分布,并详细分析了不同狭窄模型、Reynolds 数和Womersley 数对血液流动规律的影响,从而为研究血管壁病变和动脉硬化形成机制提供了有用的理论参考.关键词 格子Boltzmann 方法;Reynolds 数;Womersley 数;脉动流;动脉狭窄3北京师范大学青年科学基金资助项目通信作者收稿日期:2009205219 格子Boltzmann 方法(lattice Boltzmann met hod ,简称LBM )是20世纪80年代迅速发展起来的一种新的流体动力学数值模拟方法[122].与以宏观连续方程的离散化为基础的传统数值方法不同,LBM 从微观层次出发,采用统计物理方法得出流体的宏观特性,而且在可操作性方面,它计算方便,编程易于实现,边界易于处理等优点已经得到广泛地证实.由于心血管疾病多集中于具有复杂几何形状和具有复杂流动特性的区域,流动区域和剪切应力的分布对理解、诊断和治疗这种疾病有很重要的作用.近年来,LBM 在血液动力学方面的应用越来越受到重视[326].本文的主要工作是用格子Boltzmann 方法模拟二维轴对称狭窄血管内脉动流的流动特性.首先对狭窄血管内定常流特性进行了研究,模拟比较不同狭窄模型和不同Reynolds 数对管壁切应力、压强和压力梯度分布的影响.然后对二维轴对称狭窄血管内脉动流的流动特性进行了研究,模拟比较在改变Reynolds 数、Womersley 数时动脉血流的流动特性,找到动脉血流的非定常性对狭窄血管中流场速度、压强和剪切应力分布的影响,从而对常见的心血管疾病发展机制给出物理解释,为进一步分析动脉粥样硬化的形成、发展及其影响提供新的研究方法和理论参考.1 二维轴对称狭窄血管内定常流特性的研究111 管壁几何模型 假定血管的狭窄处为轴对称,如图1所示,狭窄形状采用常用的余弦形状,即y =h2[1+co sπL(x -x 0)],(1)图1 二维轴对称余弦狭窄模型 其中h 是狭窄的最大高度,对应于x =x 0处,L 是狭窄总长度的一半,L x 是血管段的长度,L y 是狭窄发生前的血管宽度.112 数值计算 模拟中,计算网格选为N x ×N y =300×40,狭窄中心处为x 0=121,通过调整h 和L 来控制血管狭窄程度.血管出入口采用压强边界条件[7],管壁边界采用Mei 改进的曲线边界条件[8].为了研究不同狭窄情况下管壁的切应力、压强和压强梯度的变化规律,我们选择3个不同的狭窄模型,如表1.表1 不同的狭窄模型狭窄模型M1M2M3狭窄高度h L y /8L y /4L y /4狭窄长度2L16h 8h 16h 在保证Reynolds 数(Re =ρUL y μ=UL yν,ν=μ/ρ为流体运动学黏滞系数,U 为入口附近的平均速度)一定时,计算得3种模型管壁切应力、压强和压强梯度见 140 北京师范大学学报(自然科学版)第46卷 图2～4.Re =114,狭窄中心x 0=121.图2 3种狭窄模型下管壁切应力分布 从图2中可以看出,管壁切应力振荡的负峰值在靠近狭窄中心(x 0=121)的上游,这个峰值达到一定值后,该部位血管内皮组织易发生机械应力损伤.当狭窄长度一定时,狭窄高度越大,切应力的负峰值越大,如图2中的M1和M2;当狭窄高度一定时,狭窄长度越短,切应力的负峰值越大,如图2中的M2和M3.同时也可以看出在狭窄处的下游切应力变小,特别是M2,血液容易在此处发生流体分离.模拟得到狭窄区域的压强和压强梯度分布如图3和4所示.在相同狭窄长度下,狭窄高度越大,血管狭窄上游压强下降越大,下游压强上升越大,同时狭窄区域前后的压强落差越大,如图3中的M1和M2.另一方面,在相同狭窄高度下,狭窄长度越长,血管狭窄上游压强下降越大,同时狭窄区域前后的压强落差越大,如图3中的M2和M3.压强梯度在狭窄区域波动加图3 管壁上压强分布(Re =114),p 0是狭窄发生前的压强,u 0是x =20处的中心流速 图4 管壁上的压强梯度分布(Re =114) 剧,压强梯度波动最大的是狭窄模型M2(图4),其对应的切应力负峰值也为最大值,狭窄部位管壁切应力与压强梯度的变化规律具有相似性.选择模型M2,比较管壁切应力和狭窄附近的流场分布随Re 的变化规律,如图5和6.从图5中可以看出,狭窄模型一定时,随着Re 的增加,管壁切应力增大,在狭窄区域的下游,切应力的增加相对减小,这是由于出现了流体分离,如图6的流场分布.图6显示了模型M2在不同Re 下狭窄附近的流场分布,可以看出,随着Re 的增大,在狭窄下游管壁处出现流动分离区,且Re 越大,流动分离区越大.113 分析与结论 通过改变参数,我们获得了大量有关狭窄血管中的流场的信息.模拟结果表明,血管局部图5 管壁切应力随Reynolds 数的变化曲线(狭窄模型M2) 第2期张立换等:格子Boltzmann方法模拟二维轴对称狭窄血管内的脉动流141图6 不同Re下的流场分布(M2,Re=114、215、318)狭窄会对血液的流动状态产生明显的影响,从而带来一系列的生理和病理方面的复杂变化.例如,动脉硬化斑块主要发生在几何形状急剧变化和高Re流动状态的血管内.在动脉硬化斑块发展的初期,血管狭窄度比较小,对于黏度是常数的血液流体,其Re比较小,无流动分离,管壁切应力可能达到临界应力值,对狭窄上游血管壁内皮细胞造成损伤,使壁面进一步异常增生,导致血管狭窄度增加,进而导致此处流动Re的增加.当血管狭窄增大到一定值时,在狭窄下游管壁附近就会有流动分离区形成,在该区域内血液会发生滞留,血液中的血小板和纤维蛋白就会沉积,并在血管壁处形成网络结构致使血液中的脂质颗粒沉积,而最终导致动脉粥样硬化现象的出现.同时,狭窄度较大时,对应的压力梯度的值也会较大,也可以反映病变血管的异常血液流动情况.2 二维轴对称狭窄血管内脉动流的流动特性选择模型M2为研究对象,模拟中选取周期T=10000,流动的Womersley数(α=L y2ων,ω=2πf=2π/T是脉动的角频率)为α=31357,入口压强随时间周期性变化,即p(0,t)=Δp cosωt+p out,Δp为一常量,出口压强pout设为定值,图7显示一个周期8个不同时刻的脉动流管道中心中轴线上的压强分布.从图7中可以看出,中轴线上的压强不是线性变化,在靠近狭窄部位压强下降幅度明显增加,在最大狭窄处附近压强出现极小值,狭窄下游压强又逐渐回升,远离狭窄后,压强变化逐渐恢复类直管变化趋势,并且压强随时间的波动存在一定的滞后,如图中1/8T和7/8T,2/8T和6/8T以及3/8T和5/8T不完全重合.狭窄中心x0=121,狭窄长度为78.图7 iT/8时刻中轴线上的压强分布 142 北京师范大学学报(自然科学版)第46卷 脉动流前半周期的流场分布如图8所示.从图中可以看出,在T/4时刻,在狭窄下游管壁附近开始出现流动分离区,且分离区逐渐扩大,如3T/8时刻,接着又缓慢消失,如T/2时刻,流体平滑地流过凸包.图8 脉动流在前半周期内不同时刻的流场分布 需要注意的是心脏的周期性泵血作用使动脉中的血液以脉动的形式流动,动脉中血液流动的参量———压强、流量等流动参数也会随时间变化,虽然动脉中血液的流动是脉动流而不是定常流,但动脉中血流的方向平均来说却是始终不变的,即总是从动脉流向毛细血管,再流向静脉.因此,可以把由心脏收缩和舒张所引起的动脉中的脉动流看作是一定常流分量与一振荡分量的叠加,即在图8所示的流场分布中叠加上一个定常流,最终倒流的出现时间将非常短暂,且流速很小.对应于一个周期中的不同时刻,我们发现,管壁切应力的随时间的波动也存在一定的滞后.如图9给出前半周期的切应力分布.3 结束语我们讨论了二维余弦狭窄血管中血液流动的切应力、流场速度、压强和压强梯度在不同狭窄模型和不同图9 前半周期内管壁切应力的变化曲线Re下的分布规律,所得结论与用其他实验,理论和数值模拟得到的结论相同[9211],但用LBM方法编程简单,参数易于选择,从分布函数就可以得到所有主要宏 第2期张立换等:格子Boltzmann方法模拟二维轴对称狭窄血管内的脉动流143观量,证实了LBM在此模型下的适用性.考虑到血液流动的脉动性,研究了一个脉动周期中流场的变化特点,并与定常流动比较,分析其差异.由于Womersley数的选择在血流参数范围内,故认为上述结论具有参考性.值得注意的是,流动分离区并不同于定常流动所述那样在管壁处停留,而是随着时间的演化,流动分离区间歇性的出现,如对α=710797的流场分布模拟显示,与α=31357的不同点是流动分离区在管壁附近产生后,随着时间的推移,又会向管轴附近发展.与定常流情况下在Re达到300后才出现明显的分离区不同,对于脉动流,在Re较小时,就已经可以观察到明显的流动分离区了.4 参考文献[1] Qian Y H,d’Humieres D,lallemand ttice B GKmodels for Navier2Stokes equation[J].Europhys Lett, 1992,17:479[2] Chen H,Chen S,Matthaeus W H.Recovery of theNavier2Stokes using a lattice2gas Boltzmann method[J].Phys Rev A,1992,45:R5339[3] Artoli A M,Kandhai D,Hoef sloot H ttice B GKsimulations of flow in a symmetric bif urcation[J].FutureG eneration Computer Systems,2004,20:909[4] Boyd J,Buick J,Cosgrove J A,et al.Application of thelattice Boltzmann model to simulated stenosis growth in a two2dimensional carotid artery[J].Phys Med Biol,2005, 50:4783[5] Li H B,Fang H P,Lin Z ttice Boltzmannsimulation on particle suspensions in a two2dimensional symmetric stenotic artery[J].Phys Rev E,2004,69: 031919[6] 康秀英,刘大禾,周静,等.用格子Boltzmann方法模拟动脉分叉流场[J].北京师范大学学报:自然科学版,2005, 41(4):364[7] Z ou Q,He X.On pressure amd velocity boundaryconditions for the lattice Boltzmann B GK model[J].Phys Fluids,1997,9(6):1591[8] Mei R,L uo L S,L uo Shyy W.An accurate curvedboundary treatment in the lattice Boltzmann method[J].J Comput Phys,1999,155:307[9] 姚力,李大治.刚性轴对称狭窄血管内压强及其梯度的研究[J].应用数学和力学,2006,27(3):311[10] 刘国涛,王先菊,艾保全,等.狭窄动脉血管中Poiseuille流动对管壁切应力的影响[J].中山大学学报:自然科学版,2004,4(6):29[11] 秦杰,刘辉,孙利众,等.刚性狭窄管内血流压力分布的研究[J].生物力学,1989,4(6);57SIMU LATING B LOOD FLOW IN A TWO2DIMENSIONALSYMMETRIC STENOTIC ARTER Y BYTHE LATTICE BOL TZMANN METH ODZHAN G Lihuan KAN G Xiuying J I Yupin(Depart ment of Physics,Beijing Normal University,100875,Beijing,China)Abstract In t his st udy t he lattice Boltzmann met hod has been applied to a two2dimensional symmet ric stenotic artery.The velocity,p ressure and shear st ress distribution of blood flow were simulated and compared when p ulsatio n over t he blood was added.We have observed t he impact of blood flow when changing t he steno sis struct ure,Reynolds number and Womersley number.These data provide a p hysical explanation for blood vessel lesions and arterio sclero sis.K ey w ords lattice Boltzmann met hod;Reynolds number;Womersley number;p ulsating blood;steno sed artery。
格子boltzmann方法
格子boltzmann方法格子玻尔兹曼方法是一种常用的数值计算方法，它主要用于模拟稀薄气体等流体力学问题。
下面我将从方法原理、模拟过程和应用领域三个方面详细介绍格子玻尔兹曼方法。
首先，格子玻尔兹曼方法基于玻尔兹曼方程和格子Boltzmann方程，通过将连续的物理系统离散化为网格系统进行模拟。
网格系统中的每个格子代表一个微观粒子的状态，而碰撞、传输和外部力的作用通过计算和更新这些格子的状态来实现。
该方法主要包含两个步骤：碰撞和传输。
在碰撞过程中，格子中的粒子通过相互作用和碰撞来改变其速度和方向，从而模拟了分子之间的碰撞过程。
在传输过程中，碰撞后的粒子根据流体的速度场进行移动，从而模拟了背景流场对粒子运动的影响。
其次，在格子玻尔兹曼方法中，模拟的过程可以简化为两个部分：演化和碰撞。
在每个时间步长内，系统首先根据粒子速度和位置的信息计算出相应格点上的分布函数，然后通过碰撞步骤更新这些分布函数以模拟粒子之间的碰撞效应。
通过迭代演化和碰撞步骤，系统的宏观行为可以得到。
格子玻尔兹曼方法中最常用的碰撞操作是BGK碰撞算子，它根据粒子的速度和位置信息计算出新的分布函数，并用该新分布函数代替原来的分布函数。
而在传输过程中，粒子通过碰撞后得到的新速度和方向进行移动。
最后，格子玻尔兹曼方法在流体力学领域具有广泛的应用，特别是在稀薄气体流动、微纳尺度流动和多相流等问题中。
由于其适用于模拟分子尺度和介观尺度流动问题，因此在利用普通的Navier-Stokes方程难以模拟的问题中表现出了良好的效果。
此外，格子玻尔兹曼方法还可以用于模拟流动中的热传导问题、气体分子在多孔介质中的传输问题以及颗粒与流体相互作用等多种复杂流动现象。
近年来，随着计算机性能的不断提高，格子玻尔兹曼方法也得到了快速发展，在模拟大规模真实流体问题方面取得了不错的结果。
总结来说，格子玻尔兹曼方法通过将连续的物理系统离散化为网格系统，模拟粒子碰撞和传输过程，实现了对流体力学问题的数值模拟。
格子Boltzmann方法三种边界格式的对比分析
中分别使 用半步长反 弹、 非平衡反 弹、 以及 非平衡 外推 三种边界处理格 式 ， 并得 到 了不同格 式对应 的流 线分
布 ， 函数最 小值、 流 涡心 坐标、 几何 中心线速度分布等。通过将所得结果与基 准解进行 比较 ， 三种边界格 式 就
的计算效率 ， 计算精 度、 以及计算稳 定性等方 面进行 了讨论和分析 ， Ｌ Ｍ 计算 中边界格 式 的选择 提供 了有 为 Ｂ
Ａｂ ｔ ａ ｔ ｎ ｔ ｉ ｐ ｐ ｒ ｗｏ ｄ ｍｅ ｓｏ ａ ｉ － ｒ ｅ ａ ｉ ｏ ｓ ｓ ｌａｅ ｙ ｌｔ ｃ ｌ ｍａ ｎ ｍｅｈ ｄ Ｔ ｒ ｅ ｂ ｕ ｄ ｒ ｓｒ ｃ ：Ｉ ｈ ｓ ａ ｅ ，ｔ － ｉ ｎ ｉｎ ｌｌ ｄ ｉ ｎ ｃ ｖｔ ｆ ｗ ｉ ｉ ｄ ｖ ｙｌ ｍｕ ｔｄ ｂ ａｔ ｅ Ｂｏｔ ｎ ｔｏ ． ｈ ｅ ｏ ｎ ａ ｙ ｉ ｚ
益 的参考。 关键 词 ： 格子 Ｂ ｈｍａｎ方法； ｏｚ ｎ 边界处理格 式； 半步长反弹格式 ； 非平衡 反弹格 式； 非平衡 外推格式 中图分类号 ： ７７ １ ０ ５ ． 文献标识码 ： Ａ 文章编号 ：０７ ４ １ （０ ２ ０ — ０ ８ ０ １０ － ４ ４ ２ １ ） １ ０ １ — ５ Ｔ ｅｉ ｍｐ ｒ ｔ ｅａ ａｙ ｉ ｆｔ ｒ ｅｂ ｕ ｄ ｒ ｃ ｅ ａ ｅ ｉ Ｉｔ ｅ ｌ ｍａ ｎ ｍｅｈ ｄ ｈ ｏ ａ ａｉ ｎ ｌｓ ｏ ｅ ｏ ｎ ａ ｙ ｓ ｈ ｍｅｂ ｓｄ Ｏｌ ａ ｆ ｅＢｏｔ ： ｖ ｓ ｈ ｉ ｚ ｎ ｔｏ
ｂ ｕ ｄ ｒ ｃ ｅ ｅ ｅ ｔ ｇ ｉ ＢＭ ｉ ｌ ｔ ｎ ｏ ｎ ａｙ ｓ ｈ ｍｅ ｓ ｌ ｃｉ Ｌ ｎ ｎ ｓ ｍｕａｉ ． ｏ
Ｋｅ ｒ ｓ：ｌｔ ｃ ｏｔｍａ ｎ ｍｅｈ ｄ ｙ ｗｏ ｄ ａｔ ｅ Ｂ ｌ ｎ ｔ ｏ ；ｂ ｕ ｄ ｒ ｃ ｅ ；ｈ ｌ — ｙ ｂ ｕ ｃ —ｂ ｃ ｃ ｅ ；ｎ ｎ ｑ ｉ ｂ ｉ ｍ ｏ ｃ —ｂ ｃ ｉ ｚ ｏ ｎ ａ ｙ ｓ ｈ ｍｅ ａ ｗａ ｏ ｎ ｅ ａ ｋ ｓ ｈ ｍｅ ｏ ｅ ｕ ｌ ｒｕ ｂ ｕ ｅ ａ ｋ ｆ ｉ ｎ ｓ ｈ ｍｅ；ｎ ｎ ｑ ｉ ｂ ｉ ｍ ｘｒ ｐ ｌｔ ｎ ｓ ｈ ｍｅ ｃｅ ｏ ｅ ｕｌ ｒｕ ｅ ￣ ｏ ｏ ｕａ ｏ ｏ ｉｅｅ ｔｓｈ ｍｅ ａｅ ｄｓｕｓ ｄ ｔｂＨ ｆｃｍｐ ｔｔ ｎ ｆｒｄｆ ｒｎ ｃ ｅ ｒ ｉｃ ｓｅ ．ｗｈｃ ｒ ｘ ｅｔｄ ｔ ｍｄｄ ａｕ ｂｅ ｒｆｒｎ ｅ ｔ ｉ ｆ ｉｈ ａｅ ｅｐ ｃｅ ｏ ｐ ｅ ｖ ａ ｌ ｅｅｅ ｃ ｏ ｌ
顶盖驱动流计算程序
!此?程ì序ò用?QUICK格?式?求ó解a顶￥盖?驱y动ˉ流ⅰ?!对?压1力求ó解a利?用?人?工¤压1缩?法ぁ?求ó解aprogram QUICK_cavityimplicit none!mx和ímy分?别纄表括?示?网?格?节ú点?数簓integer :: i,j,mx,my,numcharacter(len=15) :: name,name1!c2表括?示?人?工¤压1缩?系μ数簓的?平?方?,dt非?稳è态?前°后ó时骸?间?间?隔?!dx和ídy表括?示?网?格?间?距à,x1和íy1表括?示?边?长¤,err表括?示?判D断?人?工¤压1缩?法ぁ?求ó解a收?敛?的?标括?准?!u0表括?示?顶￥盖?运?动ˉ的?速ù度è,relax表括?示?人?工¤压1缩?系μ数簓real(kind=8) :: c2,relax,dt,dx,dy,x1,y1,err,abc,u0!density表括?示?流ⅰ?体?密ü度è,Viscosity表括?示?流ⅰ?体?粘3度è,Re表括?示?雷ぁ?诺μ数簓real(kind=8) :: density,Viscosity,Re!表括?示?t时骸?刻ì的?速ù度è和í压1力值μreal(kind=8),allocatable :: u(:,:),v(:,:),p(:,:)!表括?示?t+dt时骸?刻ì的?速ù度è和í压1力值μreal(kind=8),allocatable :: un(:,:),vn(:,:),pn(:,:)!最?终?各÷节ú点?上?的?速ù度è和í压1力值μreal(kind=8),allocatable :: uc(:,:),vc(:,:),pc(:,:)!定¨义?涡D量?和í流ⅰ?函ˉ数簓real(kind=8),allocatable:: flow_function_u(:,:),flow_function_v(:,:),vortex_function(:,:)!给?各÷节ú点?定¨义?坐?标括?矩?阵óreal(kind=8),allocatable :: x(:),y(:)open(3,file="parameter.dat")read(3,*) mx,my,x1,y1 !读á入?网?格?节ú点?数簓和í边?长¤read(3,*) relax,dt,u0 !读á入?压1缩?系μ数簓,时骸?间?间?隔?和í顶￥盖?移?动ˉ速ù度èread(3,*) density,Viscosity !读á入?密ü度è和í粘3度èclose(3)allocate(u(mx,my+1),v(mx+1,my),p(mx+1,my+1))allocate(un(mx,my+1),vn(mx+1,my),pn(mx+1,my+1))allocate(uc(mx,my),vc(mx,my),pc(mx,my))allocate(flow_function_u(mx,my),flow_function_v(mx,my),vortex_function(mx,my))allocate(x(mx),y(my))dx=x1/(mx-1)dy=y1/(my-1)num=0err=100.0c2=relax**2Re=u0*density*x1/Viscosity!对?流ⅰ?场?进?行D初?始?化ˉdo i=1,mx+1,1do j=1,my+1,1p(i,j)=1.0end doend dodo i=1,mx,1do j=1,my+1,1u(i,j)=0.0if(j==my+1) u(i,j)=3.0*u0/2.0if(j==my) u(i,j)=u0/2.0end doend dodo i=1,mx+1,1do j=1,my,1v(i,j)=0.0end doend do!利?用?人?工¤压1缩?法ぁ?求ó解a流ⅰ?场?值μdo while(err>0.00001.and.num<1000000)err=0.0!调獭?用?quick子哩?程ì序ò中D的?QUICK离?散Α?格?式?计?算?un和ívn的?值μcall quick(u0,u,v,p,mx,my,dx,dy,dt,density,Viscosity,un,vn)!调獭?用?calp子哩?程ì序ò求ó解a压1强?pn值μcall calp(un,vn,pn,mx,my,dx,dy,dt,c2,p)!校￡验é流ⅰ?场?信?息￠，?得?到?收?敛?判D断?准?则òerr，?同?时骸?更ü新?u、￠v、￠p !利?用?单蹋?位?时骸?间?间?隔?前°后ó时骸?刻ì之?间?的?差?值μ的?绝?对?值μ作痢?为a 判D据Ydo i=1,mx,1do j=1,my+1,1abc=abs(un(i,j)-u(i,j))/dtif(abc>err) err=abcu(i,j)=un(i,j)end doend dodo i=1,mx+1,1do j=1,my,1abc=abs(vn(i,j)-v(i,j))/dtif(abc>err) err=abcv(i,j)=vn(i,j)end doend dodo j=1,my+1,1abc=(abs(pn(i,j)-p(i,j))/c2)/dtif(abc>err) err=abcp(i,j)=pn(i,j)end doend dowrite(*,*) "error=",errnum=num+1write(*,*) numend do!---------循-环·结á束?-------------------------------------------------------do i=1,mx,1x(i)=(i-1)*dxend dodo j=1,my,1y(j)=(j-1)*dyend do!计?算?节ú点?上?的?涡D量?do i=2,mx-1,1do j=2,my-1,1vortex_function(i,j)=(un(i,j+1)-un(i,j))/dy-(vn(i+1,j)-vn(i,j))/dx end doend dodo j=1,my,1vortex_function(1,j)=(un(1,j+1)-un(1,j))/dy-(vn(2,j)-vn(1,j))/dxvortex_function(mx,j)=(un(mx,j+1)-un(mx,j))/dy-(vn(mx+1,j)-vn(mx,j))/dx end dodo i=2,mx-1,1vortex_function(i,1)=(un(i,2)-un(i,1))/dy-(vn(i+1,1)-vn(i,1))/dxvortex_function(i,my)=(un(i,my+1)-un(i,my))/dy-(vn(i+1,my)-vn(i,my))/dx end do!计?算?节ú点?上?的?速ù度è值μdo i=1,mx,1do j=1,my,1uc(i,j)=0.5*(u(i,j)+u(i,j+1))vc(i,j)=0.5*(v(i,j)+v(i+1,j))pc(i,j)=0.25*(p(i,j)+p(i,j+1)+p(i+1,j)+p(i+1,j+1))end doend do!计?算?节ú点?上?的?流ⅰ?函ˉ数簓flow_function_u=0.0flow_function_v=0.0do i=2,mx-1,1flow_function_u(i,j)=dy*uc(i,j)+flow_function_u(i,j-1)flow_function_v(i,j)=-dx*vc(i,j)+flow_function_u(i-1,j)end doend do!输?出?数簓据Y到?文?件t夹D中Dwrite(name,"(f10.4)") Reopen(10,file='Reout'//trim(adjustl(name))//'.dat')write(10,*) 'TITLE = "result"'write(10,*) 'VARIABLES = "X","Y","U","V","P"'write(10,*) 'ZONE I=',mx,'J=',my,'F=POINT'write(10,"(5(f15.8,5x))") ((x(i),y(j),uc(i,j),vc(i,j),pc(i,j),i=1,mx),j=1,my)close(10)write(name1,"(f10.4)") Reopen(20,file='Refunction'//trim(adjustl(name1))//'.dat')write(20,*) 'TITLE = "result"'write(20,*) 'VARIABLES = "X","Y","FLOW_FUNCTION_U","FLOW_FUNCTION_V","VORTEX_FUNCTION"'write(20,*) 'ZONE I=',mx,'J=',my,'F=POINT'write(20,"(5(f15.8,5x))")((x(i),y(j),flow_function_u(i,j),flow_function_v(i,j),vortex_function(i,j),i=1,mx),j=1, my)close(20)stopend!利?用?一?阶×迎?风?格?式?和íQUICK格?式?求ó解a速ù度è值μsubroutine quick(u0,u,v,p,mx,my,dx,dy,dt,density,Viscosity,un,vn)implicit noneinteger :: i,j,mx,myreal(kind=8) :: dx,dy,dt,Viscosity,density,u0real(kind=8) :: u(mx,my+1),v(mx+1,my),p(mx+1,my+1)real(kind=8) :: un(mx,my+1),vn(mx+1,my)real(kind=8) :: fw,fe,fs,fn,df,dw,de,ds,dnreal(kind=8) :: aw,aww,ae,aee,as,ass,an,ann,apreal(kind=8),external :: alfa!求ó解ax方?向ò的?速ù度èundo i=3,mx-2,1do j=3,my-1,1!求ó解aQUICK格?式?中D的?系μ数簓fw=0.5*density*(u(i-1,j)+u(i,j))*dyfe=0.5*density*(u(i,j)+u(i+1,j))*dyfs=0.5*density*(v(i,j-1)+v(i+1,j-1))*dxfn=0.5*density*(v(i,j)+v(i+1,j))*dxdf=fe-fw+fn-fsdw=Viscosity*dy/dxde=Viscosity*dy/dxds=Viscosity*dx/dydn=Viscosity*dx/dyaw=dw+0.75*alfa(fw)*fw+0.125*alfa(fe)*fe+0.375*(1.0-alfa(fw))*fwaww=-0.125*alfa(fw)*fwae=de-0.375*alfa(fe)*fe-0.75*(1.0-alfa(fe))*fe-0.125*(1.0-alfa(fw))*fwaee=0.125*(1.0-alfa(fe))*feas=ds+0.75*alfa(fs)*fs+0.125*alfa(fn)*fn+0.375*(1.0-alfa(fs))*fsass=-0.125*alfa(fs)*fsan=dn-0.375*alfa(fn)*fn-0.75*(1.0-alfa(fn))*fn-0.125*(1.0-alfa(fs))*fsann=0.125*(1.0-alfa(fn))*fnap=aw+aww+ae+aee+as+ass+an+ann+dfun(i,j)=u(i,j)+dt/(dx*dy)*(-ap*u(i,j)+aw*u(i-1,j)+ae*u(i+1,j)+as*u(i,j-1)+an*u(i,j+1)+& aww*u(i-2,j)+aee*u(i+2,j)+ass*u(i,j-2)+ann*u(i,j+2))-dt*(p(i+1,j)-p(i,j))/dx end doend do!----------------------------------------------------------------------!内ú层?边?界?由?一?阶×迎?风?离?散Α?格?式?计?算?得?到?j=2do i=3,mx-2,1call upbound_u(u,v,p,mx,my,dx,dy,dt,density,Viscosity,i,j,un)end doj=mydo i=3,mx-2,1call upbound_u(u,v,p,mx,my,dx,dy,dt,density,Viscosity,i,j,un)end doi=2do j=2,my,1call upbound_u(u,v,p,mx,my,dx,dy,dt,density,Viscosity,i,j,un)end doi=mx-1do j=2,my,1call upbound_u(u,v,p,mx,my,dx,dy,dt,density,Viscosity,i,j,un)end do!--------------------------------------------------------------------!计?算?外猘层?边?界?do i=2,mx-1,1un(i,1)=-un(i,2)un(i,my+1)=2.0*u0-un(i,my)end dodo j=1,my+1,1un(1,j)=0.0un(mx,j)=0.0end do!--------------------------------------------------------------------!求ó解ay方?向ò的?速ù度èvndo i=3,mx-1,1do j=3,my-2,1!求ó解aQUICK格?式?中D的?系μ数簓fw=0.5*density*(u(i-1,j)+u(i-1,j+1))*dyfe=0.5*density*(u(i,j)+u(i,j+1))*dyfs=0.5*density*(v(i,j-1)+v(i,j))*dxfn=0.5*density*(v(i,j)+v(i,j+1))*dxdf=fe-fw+fn-fsdw=Viscosity*dy/dxde=Viscosity*dy/dxds=Viscosity*dx/dydn=Viscosity*dx/dyaw=dw+0.75*alfa(fw)*fw+0.125*alfa(fe)*fe+0.375*(1.0-alfa(fw))*fwaww=-0.125*alfa(fw)*fwae=de-0.375*alfa(fe)*fe-0.75*(1.0-alfa(fe))*fe-0.125*(1.0-alfa(fw))*fwaee=0.125*(1.0-alfa(fe))*feas=ds+0.75*alfa(fs)*fs+0.125*alfa(fn)*fn+0.375*(1.0-alfa(fs))*fsass=-0.125*alfa(fs)*fsan=dn-0.375*alfa(fn)*fn-0.75*(1.0-alfa(fn))*fn-0.125*(1.0-alfa(fs))*fsann=0.125*(1.0-alfa(fn))*fnap=aw+aww+ae+aee+as+ass+an+ann+dfvn(i,j)=v(i,j)+dt/(dx*dy)*(-ap*v(i,j)+aw*v(i-1,j)+ae*v(i+1,j)+as*v(i,j-1)+an*v(i,j+1)+& aww*v(i-2,j)+aee*v(i+2,j)+ass*v(i,j-2)+ann*v(i,j+2))-dt*(p(i,j+1)-p(i,j))/dy end doend do!------------------------------------------------------------------------------!内ú层?边?界?由?一?阶×迎?风?离?散Α?格?式?计?算?得?到?j=2do i=3,mx-1,1call upbound_v(u,v,p,mx,my,dx,dy,dt,density,Viscosity,i,j,vn)end doj=my-1do i=3,mx-1,1call upbound_v(u,v,p,mx,my,dx,dy,dt,density,Viscosity,i,j,vn)end doi=2do j=2,my-1,1call upbound_v(u,v,p,mx,my,dx,dy,dt,density,Viscosity,i,j,vn) end doi=mxdo j=2,my-1,1call upbound_v(u,v,p,mx,my,dx,dy,dt,density,Viscosity,i,j,vn) end do!----------------------------------------------------------------------!计?算?外猘层?边?界?do i=2,mx,1vn(i,1)=0.0vn(i,my)=0.0end dodo j=1,my,1vn(1,j)=-vn(2,j)vn(mx+1,j)=-vn(mx,j)end doreturnend subroutine!-----------------------------------------------------------------!定¨义?一?个?函ˉ数簓function alfa(x)implicit nonereal(kind=8) :: alfa,xif(x>0.0) thenalfa=1.0elsealfa=0.0end ifreturnend!利?用?一?阶×迎?风?格?式?求ó解a内ú层?边?界?上?的?速ù度è值μsubroutine upbound_u(u,v,p,mx,my,dx,dy,dt,density,Viscosity,i,j,un) implicit noneinteger :: i,j,mx,myreal(kind=8) :: dx,dy,dt,density,Viscosityreal(kind=8) :: u(mx,my+1),v(mx+1,my),p(mx+1,my+1)real(kind=8) :: un(mx,my+1)real(kind=8) :: df,dw,de,ds,dnreal(kind=8) :: aw,ae,as,an,apdw=Viscosity*dy/dxde=Viscosity*dy/dxds=Viscosity*dx/dydn=Viscosity*dx/dyaw=dw+max(0.5*density*(u(i-1,j)+u(i,j))*dy,0.0)ae=de+max(0.0,-0.5*density*(u(i,j)+u(i+1,j))*dy)as=ds+max(0.5*density*(v(i,j-1)+v(i+1,j-1))*dx,0.0)an=dn+max(0.0,-0.5*density*(v(i,j)+v(i+1,j))*dx)df=0.5*density*(u(i+1,j)-u(i-1,j))*dy+0.5*density*(v(i,j)+v(i+1,j)-v(i,j-1)-v(i+1,j-1)) *dxap=aw+ae+as+an+dfun(i,j)=u(i,j)+(dt/(dy*dx))*(-ap*u(i,j)+aw*u(i-1,j)+ae*u(i+1,j)+as*u(i,j-1)+an*u(i,j+1) )-dt*(p(i+1,j)-p(i,j))/dxreturnend subroutine!利?用?一?阶×迎?风?格?式?求ó解a内ú层?边?界?上?的?速ù度è值μsubroutine upbound_v(u,v,p,mx,my,dx,dy,dt,density,Viscosity,i,j,vn)implicit noneinteger :: i,j,mx,myreal(kind=8) :: dx,dy,dt,density,Viscosityreal(kind=8) :: u(mx,my+1),v(mx+1,my),p(mx+1,my+1)real(kind=8) :: vn(mx+1,my)real(kind=8) :: df,dw,de,ds,dnreal(kind=8) :: aw,ae,as,an,apdw=Viscosity*dy/dxde=Viscosity*dy/dxds=Viscosity*dx/dydn=Viscosity*dx/dyaw=dw+max(0.5*density*(u(i-1,j)+u(i-1,j+1))*dy,0.0)ae=de+max(0.0,-0.5*density*(u(i,j)+u(i,j+1))*dy)as=ds+max(0.5*density*(v(i,j-1)+v(i,j))*dx,0.0)an=dn+max(0.0,-0.5*density*(v(i,j)+v(i,j+1))*dx)df=0.5*density*(u(i,j)+u(i,j+1)-u(i-1,j)-u(i-1,j+1))*dy+0.5*density*(v(i,j+1)-v(i,j-1)) *dxap=aw+ae+as+an+dfvn(i,j)=v(i,j)+(dt/(dy*dx))*(-ap*v(i,j)+aw*v(i-1,j)+ae*v(i+1,j)+as*v(i,j-1)+an*v(i,j+1) )-dt*(p(i,j+1)-p(i,j))/dyreturnend subroutine!人?工¤压1缩?法ぁ?求ó解a下?一?时骸?刻ì的?压1强?值μsubroutine calp(un,vn,pn,mx,my,dx,dy,dt,c2,p)implicit noneinteger :: i,j,mx,myreal(kind=8) :: dx,dy,dt,c2real(kind=8) :: un(mx,my+1),vn(mx+1,my),pn(mx+1,my+1)real(kind=8) :: p(mx+1,my+1)!根ù据Yun和ívn求ó解a压1强?pn的?值μ!利?用?连?续?性?方?程ì求ó解a压1力do i=2,mx,1do j=2,my,1pn(i,j)=p(i,j)-dt*c2*((un(i,j)-un(i-1,j))/dx+(vn(i,j)-vn(i,j-1))/dy) end doend do!利?用?虚é拟a边?界?条?件t求ó解a压1力值μdo i=2,mx,1pn(i,1)=pn(i,2)pn(i,my+1)=pn(i,my)end dodo j=1,my+1,1pn(1,j)=pn(2,j)pn(mx+1,j)=pn(mx,j)end doreturnend subroutine。
流动沸腾的格子boltzmann方法模拟
流动沸腾的格子boltzmann方法模拟随着科学技术的不断发展，对复杂流体动力学过程的研究以及对快速过程的模拟变得更加重要，由于大多数过程都是多尺度耦合的，而且有一个高度不确定的初始条件，因此传统的办法往往很难分析以及模拟。
格子Boltzmann方法（LBM）是一种计算流体机械系统的模拟方法，它允许将多物理学过程直接耦合在一起，同时也满足计算机能力的限制。
在过去的几十年里，格子Boltzmann方法已经成功应用于各种复杂系统中，并且得到了广泛的应用。
在本文中，我们将重点关注格子Boltzmann方法模拟流动沸腾过程。
流动沸腾过程是指流体作为一个定常动态系统，当流体间的温度波动超过沸点时，出现蒸发和汽化的共同现象。
在实际的工程应用中，这种现象可以被用来模拟大气中的蒸发和汽化过程，也可以用来模拟熔炼过程中的熔融气体的性质，当物体在极端条件下的表面上出现沸腾时，也可以用来模拟这种情况。
格子Boltzmann方法是一种经典的集总体方法，通常是将空间区域划分为多个格子大小的单元格，然后模拟每个单元格内流体粒子的运动方式。
首先，根据Boltzmann方程，定义系统中每个单元格内流体粒子的分布函数，然后根据撞击原理，建立系统的碰撞模型，最后，利用快速Fourier变换（FFT）计算出每个单元格的空间和时间变化，最终得到系统中流体粒子分布函数的空间和时间变化，从而模拟出系统的动力学过程。
格子Boltzmann方法模拟流动沸腾过程的优势在于它可以解决复杂的动力学问题，在流动沸腾过程中，可以使用格子Boltzmann方法考虑到热传导和蒸汽压力的影响，同时考虑到蒸发和汽化等热力学效应，从而得到较准确的结果。
此外，采用格子Boltzmann方法模拟流动沸腾过程，具有计算效率高、可靠性高、精度高的优势。
由于上述优点，格子Boltzmann方法在模拟流动沸腾过程方面已经取得了许多成功的应用。
例如，Zhao等对水汽-气体混合流动的二相沸腾过程进行了模拟，他们利用格子Boltzmann方法模拟了水汽-气体混合的沸腾过程，比较了侧向温度波及其在慢一点的过程中的水汽沉积速率，并表明其模拟结果良好，与实际的实验结果相一致。
枝晶和气孔演化规律的格子boltzmann方法数值模拟
枝晶和气孔演化规律的格子boltzmann方法
数值模拟
格子Boltzmann方法是一种新兴的、用于模拟二维和三维
物质在复杂流变行为中的方法，在过去几年中，由于晶体和气
孔演化机制已经被深入研究，格子Boltzmann方法也受到越来
越多的关注。
首先，格子Boltzmann方法通过对一系列物理量进行条件
分布来复杂模拟晶体和气孔演化。
这一方法可以为建模和模拟
气孔和晶体形成机制提供有益依据，它能够模拟系统的温度依
赖性和宏观动力学行为的变化。
同时，格子Boltzmann方法也
能够有效地提取气孔演化和晶体结构演化的行为特征。
此外，格子Boltzmann方法可以将复杂的物理过程抽象化，从而提取出系统中比较容易理解的量化模型。
基于定属性原理，该方法可以用来模拟数据的动态变化，从而获得晶体的形貌、
尺寸和气孔的分布状况等信息，为研究晶体和气孔演化提供了
重要的参考资料。
最后，格子Boltzmann方法可以有效地解决抽象的高维物
理系统问题，其优势在于能够精准地模拟晶体和气孔演化的行为，它已经广泛应用于材料科学、工程学和生物学等领域，成
为现代互联网科学领域中解决复杂过程问题的基石。
计算流体大作业-基于C++的顶部驱动流模拟
-7-
计算流体力学大作业
第8页 共9页
图 6
Re=2000 时顶盖驱动流的流线图
图 7
Re=ห้องสมุดไป่ตู้000 时的顶盖驱动流流线图
-8-
计算流体力学大作业
第9页 共9页
图 8
Re=10000 时顶盖驱动流流线图
如图 6 至图 8 三幅图所示，基于 C/C++顶盖驱动流的流线图中，具有如下特 征： a) 方腔中心具有一个大漩涡； b) 方腔的左上角、左下角以及右下角具有次级漩涡； c) 随着雷诺数的增大，右下角的次级漩涡逐渐变小，左上角的次级漩涡逐 渐变大。
Re
ul Re ，其中顶盖流速 u = 0.1，l = 1
GAM 1 Re
4.程序组织
本文借助 D2Q9 格子波尔兹曼方法，在 Visual Studio 2012 平台上运用 C/C++语 言编程，运用有限差分方法，对上述传统问题进行了求解，并借助 Tecplot360 工具， 对所得数据进行后处理。所述程序组织如下图所示（完整程序参见所提交作业的文件 夹） ：
式中，p 为压力，v 为动力粘性系数。 根据上述分析，可以方便地设计出格子波尔兹曼算法。
（7）
3.2 顶部驱动流控制方程
方腔顶盖驱动流是考核程序的经典算例之一，其满足如下控制方程：
uu uv 2u 2u p 2 2 y x y x x uv vv 2v 2v p x 2 y 2 y y x
7.结论
本文借助 D2Q9 格子波尔兹曼方法， 在 Visual Studio 2012 平台上运用 C/C++ 语言编程，求解了顶部驱动流模型，运用 tecplot360 对所得数据进行处理，得 到了图 6 至图 8 三幅对应于不同雷诺数的流线图， 显示了该流动模型的典型现象， 验证了本文所述方法的正确性。
格子-boltzmann 顶盖流热力学模拟程序 VB版
顶盖流热力学模拟程序（VB版）概述天大热能李扬一、程序介绍本程序用格子-Boltzmann法模拟顶盖流运动的速度场和温度场分布，编写语言为VB。
程序可以实时绘制速度场及温度场的分布情况并在任意时刻输出图片，可以自由设定格子数、迭代精度等各项参数，可以随时暂停循环并定义每循环多少次时自动输出数据。
二、几个代表时刻的温度场、速度场分布（100*100格）温度场：速度场（注意最终稳定时会出现中心、左下角、右下角的漩涡）：三、程序注意事项1. 实时显示速度场及温度场大大降低了计算的速度，64*64格迭代30000次需要约2个小时，因而建议在不随时观察的情况下将绘制频率加大2. 以密度作为迭代精度判定条件时收敛性不好，因而在程序中加入不考虑终止精度的选择。
以速度作为迭代精度较为合适，程序可以选择采用密度或速度作为判定条件（对于速度场，50*50格大概在迭代15000次之后稳定，100*100格大概在迭代30000次后稳定）对于用密度作为迭代残差不收敛的问题如果有同学知道原因请不吝赐教！万分感谢！本程序还有许多许多不完善、不科学的地方，尤其计算速度很慢，望诸位高人海涵！QQ：513682028四、程序代码声明：本代码仅供学习参考之用，请有选择性、学习性的借鉴，切勿不加思考简单修改就用于作业或考试！谢谢！（我在编写计算部分时也是很头疼参考了很多程序并和同学商量，借鉴没关系，但一定要有自己的收获哦~）控件代号：程序部分(红字部分为主要计算程序)：Dim w(8), cs2, midu0, c(8, 1), tf, tg '定义全局变量Dim u0, re, pr, npDim midu, f, u, miduold, fp, ip, jp, m, n, start, g, gp, t, t0, color, xx, vPrivate Sub Combo4_click() '切换误差种类Label5.Caption = "两轮" & Right(Combo4.Text, Len(Combo4.Text) - 1) & "：未开始"End SubPrivate Sub Command1_Click() '初始化'取得雷诺数re和普朗特数prre = Val(Text1.Text)pr = Val(Text2.Text)'取得顶盖流速度u0u0 = Val(Text3.Text)'取得每行格子数npnp = Val(Text4.Text)'计算tf、tgtf = u0 * np * 3 / re + 0.5tg = 0.5 + (tf - 0.5) / pr'定义各点uReDim u(np, np, 1)ReDim t(np, np)'定义第一行水平速度为u0For i = 0 To npu(i, np, 0) = u0t(i, np) = t0Next'定义其他点x、y方向速度均为0 For i = 0 To npFor j = 0 To np - 1u(i, j, 0) = 0u(i, j, 1) = 0t(i, j) = 0NextNext'定义各点初始miduReDim midu(np, np)ReDim v(i, j)For j = 0 To npmidu(i, j) = midu0v(i, j) = (u(i, j, 0) ^ 2 + u(i, j, 1) ^ 2) ^ 0.5NextNext'计算初始feqReDim f(np, np, 8)ReDim g(np, np, 8)For i = 0 To npFor j = 0 To npFor k = 0 To 8f(i, j, k) = feq(k, midu(i, j), u(i, j, 0), u(i, j, 1))g(i, j, k) = feq(k, midu(i, j), u(i, j, 0), u(i, j, 1)) * t(i, j) NextNextNext'改变桌面控件属性Text1.Enabled = FalseText2.Enabled = FalseText3.Enabled = FalseText4.Enabled = Falsestart = 1Command1.Caption = "暂停"ElseIf start = 1 Thenstart = 2Command1.Caption = "继续"Command2.Enabled = TrueCommand3.Enabled = TrueElseIf start = 2 ThenCommand1.Caption = "暂停" Command2.Enabled = False Command3.Enabled = FalseEnd IfEnd SubPrivate Sub Command2_Click() '输出数据Call output(n - 1, 1) '调用输出数据End SubPrivate Sub Form_Load() '定义初始变量'定义w在9个方向上的数值w(0) = 4 / 9For i = 1 To 4w(i) = 1 / 9NextFor i = 5 To 8w(i) = 1 / 36Next'定义cs2cs2 = 1 / 3'定义初始密度midu0 = 1t0 = 1'定义c在9个方向上的数值c(0, 0) = 0c(0, 1) = 0c(1, 0) = 1c(1, 1) = 0c(2, 1) = 1c(3, 0) = -1c(3, 1) = 0c(4, 0) = 0c(4, 1) = -1c(5, 0) = 1c(5, 1) = 1c(6, 0) = -1c(6, 1) = 1c(7, 0) = -1c(7, 1) = -1c(8, 0) = 1c(8, 1) = -1n = 1start = 0'定义10个温度区间ReDim color(10)color(0) = RGB(255, 13, 0) color(1) = RGB(255, 77, 0) color(2) = RGB(255, 176, 0) color(3) = RGB(247, 248, 0) color(4) = RGB(176, 255, 0) color(5) = RGB(75, 255, 0) color(6) = RGB(0, 255, 176) color(7) = RGB(0, 248, 247) color(8) = RGB(0, 176, 255) color(9) = RGB(0, 75, 255) '画温度示意图For i = 0 To 9Picture2.Line (1, i * 360 + 1)-(360, i * 360 + 360), color(i), BFNext'定义控件Label5.Caption = "两轮" & Right(Combo4.Text, Len(Combo4.Text) - 1) & "：未开始" Combo2.AddItem "1e-6"Combo2.AddItem "1e-4"Combo2.AddItem "1e-5"Combo2.AddItem "1e-7"Combo3.AddItem "5"Combo3.AddItem "1"Combo3.AddItem "20"Combo3.AddItem "50"Combo4.AddItem "1各点密度和之差"Combo4.AddItem "2各点速度和之差"Combo5.AddItem "1000"Combo5.AddItem "500"Combo5.AddItem "2000"Combo5.AddItem "5000"End SubFunction feq(k, midu, u1, u2) 'feq函数Dim cu, uu2cu = c(k, 0) * u1 + c(k, 1) * u2uu2 = u1 ^ 2 + u2 ^ 2feq = w(k) * midu * (1 + cu / cs2 + cu ^ 2 / 2 / (cs2 ^ 2) - uu2 / 2 / cs2)End FunctionPrivate Sub Label10_Click() '唔MsgBox "O(∩_∩)O~祝大家开心~"End SubPrivate Sub Text3_Change() '检查速度If Val(Text3.Text) >= 1 ThenDim uoveruover = MsgBox("速度大于或等于1可能会导致计算出错，依然继续？", 52, "请考虑速度是否合适")If uover = vbNo Then Text3.Text = "0.1"End IfEnd SubPrivate Sub Timer1_Timer() '主程序写在timer里是为了方便暂停If start = 1 ThenReDim fp(np, np, 8)ReDim gp(np, np, 8)Dim midu_all, miduold_all, v_all, vold_all'主过程'取得上一次数据Select Case Left(Combo4.Text, 1)Case 1miduold = 0For i = 0 To npFor j = 0 To npmiduold_all = miduold_all + midu(i, j)NextNextCase 2vold_all = 0For i = 0 To npFor j = 0 To npvold_all = vold_all + v(i, j)NextNextEnd Select'碰撞For i = 1 To np - 1For j = 1 To np - 1For k = 0 To 8ip = i - c(k, 0)jp = j - c(k, 1)fp(i, j, k) = f(ip, jp, k) + (feq(k, midu(ip, jp), u(ip, jp, 0), u(ip, jp, 1)) - f(ip, jp, k)) / tfgp(i, j, k) = g(ip, jp, k) + (feq(k, midu(ip, jp), u(ip, jp, 0), u(ip, jp, 1)) * t(ip, jp) - g(ip, jp, k)) / tg NextNextNext'内部For i = 1 To np - 1For j = 1 To np - 1midu(i, j) = 0u(i, j, 0) = 0u(i, j, 1) = 0t(i, j) = 0For k = 0 To 8f(i, j, k) = fp(i, j, k)g(i, j, k) = gp(i, j, k)midu(i, j) = midu(i, j) + f(i, j, k)u(i, j, 0) = u(i, j, 0) + c(k, 0) * f(i, j, k)u(i, j, 1) = u(i, j, 1) + c(k, 1) * f(i, j, k)t(i, j) = t(i, j) + g(i, j, k)Nextu(i, j, 0) = u(i, j, 0) / midu(i, j)u(i, j, 1) = u(i, j, 1) / midu(i, j)t(i, j) = t(i, j) / midu(i, j)NextNext'水平边界For j = 1 To np - 1For k = 0 To 8midu(np, j) = midu(np - 1, j)f(np, j, k) = feq(k, midu(np, j), u(np, j, 0), u(np, j, 1)) + (f(np - 1, j, k) - feq(k, midu(np - 1, j), u(np - 1, j,0), u(np - 1, j, 1)))midu(0, j) = midu(1, j)f(0, j, k) = feq(k, midu(0, j), u(0, j, 0), u(0, j, 1)) + (f(1, j, k) - feq(k, midu(1, j), u(1, j, 0), u(1, j, 1))) NextNext'垂直边界For i = 0 To npFor k = 0 To 8midu(i, 0) = midu(i, 1)f(i, 0, k) = feq(k, midu(i, 0), u(i, 0, 0), u(i, 0, 1)) + (f(i, 1, k) - feq(k, midu(i, 1), u(i, 1, 0), u(i, 1, 1))) midu(i, np) = midu(i, np - 1)u(i, np, 0) = u0f(i, np, k) = feq(k, midu(i, np), u(i, np, 0), u(i, np, 1)) + (f(i, np - 1, k) - feq(k, midu(i, np - 1), u(i, np - 1,0), u(i, np - 1, 1)))NextNext'画温度场、速度场If n Mod Combo3.Text = 0 ThenPicture3.Line (1, 1)-(Picture3.Width, Picture3.Height), RGB(255, 255, 255), BFDim addpxIf np > 20 Thenaddpx = Int(xx / 2)Else: addpx = 0End Ifxx = Int(6000 / (np + 1))Picture1.Width = xx * (np + 1) + addpxPicture1.Height = xx * (np + 1) + addpxPicture3.Width = xx * (np + 1) + addpxPicture3.Height = xx * (np + 1) + addpxFor i = 0 To npFor j = 0 To npPicture1.Line (i * xx + 1, (np - j) * xx + 1)-(i * xx + xx, (np - j) * xx + xx), drawt(t(i, j)), BF Call drawu(i, j, u(i, j, 0), u(i, j, 1))NextNextEnd If'显示残差、迭代次数Dim errorSelect Case Left(Combo4.Text, 1)Case 1midu_all = 0For j = 0 To npFor i = 0 To npmidu_all = midu_all + midu(i, j)error = Abs(midu_all - miduold_all)NextNextCase 2v_all = 0For j = 0 To npFor i = 0 To npv(i, j) = (u(i, j, 0) ^ 2 + u(i, j, 1) ^ 2) ^ 0.5v_all = v_all + v(i, j)error = Abs(v_all - vold_all)NextNextEnd SelectLabel5.Caption = "两轮" & Right(Combo4.Text, Len(Combo4.Text) - 1) & "：" & error Label6.Caption = "迭代次数：" & n'迭代终止条件Dim jingduSelect Case Combo2.TextCase "1e-4"jingdu = 0.01Case "1e-5"jingdu = 0.00001Case "1e-6"jingdu = 0.000001Case "1e-7"jingdu = 0.0000001End SelectIf Check2.Value = 1 ThenIf n > 1 ThenIf error < jingdu Thenstart = 2Command1.Caption = "继续"Label5.Caption = "两轮" & Right(Combo4.Text, Len(Combo4.Text) - 1) & "：" & error Call output(-1, 1)End IfEnd IfEnd IfDoEvents '防止死机n = n + 1If Check1.Value = 1 Then '自动暂停If n - 1 = Text5.Text Thenstart = 2Command1.Caption = "继续"Command2.Enabled = TrueCommand3.Enabled = TrueMsgBox "已自动暂停！", 64, "自动暂停"End IfEnd IfIf Check3.Value = 1 Then '自动输出If (n - 1) Mod Combo5.Text = 0 ThenCall outputimgCall output(n - 1, 0)Label7.Caption = "上次自动输出数据及两场图：第" & n - 1 & "轮"End IfEnd IfEnd IfEnd SubFunction drawt(t) '返回温度对应的颜色Select Case tCase 0 To 0.05drawt = color(9)Case 0.05 To 0.1drawt = color(8)Case 0.1 To 0.25drawt = color(7)Case 0.25 To 0.4drawt = color(6)Case 0.4 To 0.5drawt = color(5)Case 0.5 To 0.6drawt = color(4)Case 0.6 To 0.7drawt = color(3)Case 0.7 To 0.8drawt = color(2)Case 0.8 To 0.9drawt = color(1)Case 0.9 To 1drawt = color(0)Case Elsedrawt = color(9)End SelectEnd FunctionFunction drawu(i, j, v0, v1) '画速度场If n Mod Combo3.Text = 0 ThenDim ux, uy, x0, y0, colort, v2x0 = i * xx + xx / 2y0 = (np - j) * xx + 1 + xx / 2v2 = (v0 ^ 2 + v1 ^ 2) ^ 0.5If v2 <> 0 ThenSelect Case v2 / u0Case 0 To 0.05colort = color(9)Case 0.05 To 0.1colort = color(8)Case 0.1 To 0.25colort = color(7)Case 0.25 To 0.4colort = color(6)Case 0.4 To 0.5colort = color(5)Case 0.5 To 0.6colort = color(4)Case 0.6 To 0.7colort = color(3)Case 0.7 To 0.8colort = color(2)Case 0.8 To 0.9colort = color(1)Case 0.9 To 2colort = color(0)Case Elsecolort = color(9)End Selectux = xx * v0 / v2 / 2uy = xx * v1 / v2 / 2Picture3.Line (x0 - ux, y0 + uy)-(x0 + ux, y0 - uy), colortPicture3.Circle (x0 + ux, y0 - uy), 2, RGB(0, 0, 0)End IfEnd IfEnd FunctionFunction output(isend, from) '输出数据Dim boxtitle, boxaddIf isend = -1 Thenisend = "final"boxtitle = "迭代结束"boxadd = vbCrLf & "若明显未达到稳态请继续"Elseboxtitle = "输出数据"End IfOpen App.Path & "/ux_" & np & "_" & n - 1 & ".txt" For Output As #1 Open App.Path & "/uy_" & np & "_" & n - 1 & ".txt" For Output As #2 Open App.Path & "/t_" & np & "_" & n - 1 & ".txt" For Output As #3 For j = np To 0 Step -1For i = 0 To npPrint #1, u(i, j, 0) & " ";Print #2, u(i, j, 1) & " ";Print #3, t(i, j) & " ";NextPrint #1, ""Print #2, ""Print #3, ""NextClose #1Close #2Close #3If from = 1 Then MsgBox "结尾为" & np & "_" & n - 1 & "的温度及速度数据成功输出至程序所在目录下" & boxadd, 64, boxtitleEnd FunctionPrivate Sub Form_QueryUnload(Cancel As Integer, UnloadMode As Integer) '确认关闭Dim res As Longres = MsgBox("确定退出吗？", 36, "确认退出")If res = vbNo ThenCancel = 1 '取消ElseFor Each f In FormsUnload fNextEnd '完全关闭End IfEnd SubPrivate Sub Command3_Click() '输出两场图Call outputimgMsgBox "结尾为" & np & "_" & n - 1 & "的温度场及速度场成功输出至程序所在目录下", 64, "输出两场图"End SubFunction outputimg() '输出两场图SavePicture Picture1.Image, App.Path & "/t_field_" & np & "_" & n - 1 & ".bmp"SavePicture Picture3.Image, App.Path & "/v_field_" & np & "_" & n - 1 & ".bmp"End Function。
应用Lattice Boltzmann方法数值模拟顶盖驱动流
ＴＡＮ ｕ ｎ Ｉ Ｙ ａ ， ，Ｐｉｇ ， Ｎ ， ｎ Ｌ４ Ｋｅ ＨＵ Ｐａｎ一 ｎ ｇ
（ ｈｎ ｑｎ ｎｖｒ ｔ ｏ ｃ ｎｅ＆ Ｔ ｃ ｎｌ ｙ ｈｎ ｑｎ ０ ３ １ １Ｃ ｏｇ ｉｇＵ ｉ ｓｙ ｆ ｉ ｃ ｅ ｉ Ｓ ｅ ｅｈ ｏ ｇ ，Ｃ ｏｇ ｉｇ １３ ； ｏ ４
４Ｇ ｓＰｏ ｕｔ ｎａ ｄＴａ ｓ ｉ ｉ ｐ ｒｔ ｎＡ ｅ ， ｈ ｎｑｎ ０ ０ ０， ｈｎ ） ａ ｒ ｃ ｏ ｎ ｒｎｍ ｓ ｏ Ｏ ｅａｏ ｒａ Ｃ ｏ ｇ ｉ ４ ０ Ｃ ｉａ ｄ ｉ ｓｎ ｉ ｇ４
Ａｂ ｔ ａ ｔ ａｔ ｅ Ｂ ｈ ｍａ ｎ ｍｅｈ ｄ ｉ ａ ｎ ｗ ｔ ｅ ｒ ｅ ｅｏ ｅ ｎ ｒ ｃ ｎ ｅ ｒ ．Ｔ ｉ ｍｅｈ ｄ ｉ ｂ ｓ ｄ ｏ ｅ Ｋｉ ｅｉ ｓｒ ｃ ：Ｌ ｔ ｃ ｏ ｚ ｎ ｔ ｏ ｓ ｅ ｏｙ ｄ ｖ ｌｐ ｄ ｉ ｅ ｅ ｔ ａ ｓ ｈ ｓ ｔ ｏ ａ ｅ ｎ ｔ ｎ ｔ ｉ ｈ ｙ ｓ ｈ ｃ
关 键 词 ： ｔ ｅ ｏｚａｎ ａ ｉ Ｂ Ｌ ｔ ｈｍ ｎ 方法 ； ｃ 数值模拟 ； 顶盖驱动；Ｇ Ｂ Ｋ碰撞模型 ；２ ９ Ｄ Ｑ 模型
Ｔｈ ａ ｆＳ ｍ ｕ ａ ｉ ｇ Ｃａ —ｄｒｖ ｎ Ｆｌ ｗ ｔ ｔｉｅ Ｂｏｔ ｍ ａ ｎ Ｅｑ ａ ｉ ｎ ｅ Ｗ ｙ ｏ ｉ ｌ ｔｎ ｐ — ｉ ｅ ｏ ｗｉｈ Ｌａ ｔｃ ｌｚ ｎ ｕ ｔｏ
非混相驱替过程的格子Boltzmann模拟
非混相驱替过程的格子Boltzmann模拟指进是一种相界面不稳定现象。
通常情况下,其诱发的主要原因是流体间的黏度差,故有时也称黏性指进。
黏性指进现象是多相渗流的重要特征之一,在很多工程应用中起着至关重要的作用,如化工、石油开采以及地下水污染的治理等领域。
以石油开采为例,该现象的过早出现会减少波及效率,造成驱替流体的窜流,降低石油采收率。
非混相驱油过程中的黏性指进现象是一个涉及到相界面捕捉、多组分多相流理论以及复杂几何流场内流体流动的复杂流体动力学问题。
其影响因素众多,涉及到多孔介质的孔隙结构、壁面润湿性以及流体的物性等。
该问题的内在产生机理至今尚不明晰,因此对于黏性指进现象的研究有着重要的学术意义和应用价值。
本文以黏性指进现象为主要研究对象,针对采用Shan-Chen (SC)模型模拟具有不同密度比和黏度比流体系统时存在的问题,进行了改进,提出了一种新的LBM多组分模型,并利用格子Boltzmann方法(Lattice Boltzmann Method, LBM)模拟单孔隙和多孔介质内的指进现象,综合考虑不同因素对于相界面形态以及驱替效率的影响,具体内容包括：(1)基于SC模型,提出一个新的可以模拟具有不同密度比和黏度比流体系统的LBM多组分模型。
新模型中,对SC模型施加作用力的方式进行改进,消除了原模型中由于作用力的施加方式引入的离散误差。
不同组分粒子的格子速度和格子声速均不同,因此,在一个时间步内,不同组分粒子的迁移距离也不同,相互之间的比例关系由流体间的密度比确定。
迁移后格子节点上的粒子分布函数由双线性插值格式确定。
给出考虑固体边界的特殊处理方法。
利用新的LBM多组分模型模拟了在不同管径的圆管内、不同体积流量下,水驱癸烷的非混相驱替过程,并与实验结果进行比较。
模拟中的平均驱替速度和驱替完成时间与实验结果相比,相对误差在3.2%以内,证明了新模型的正确性。
新模型具有良好的守恒性。
通过适当地选取边界条件,并对流场进行适当地处理,该模型可以用于模拟流体速度远小于模型本身伪速度的低速驱替过程。
格子Boltzmann
格子Boltzmann 方法模拟C/C 复合材料颗粒沉积过程罗思璇()Particle Deposition Process Simulation in C/C Composites by Lattice-Boltzmann MethodLuo Sixuan()Abstract: Lattice Boltzmann method is used here to study the particle deposition process on C/C composites surface. This method considered the boudary condition change during particle deposition. Finally, the deposition pattern is obtained. Keywords: LB Method; flow-particle coupling; C/C composites; deposition摘要：本文使用格子Boltzmann 方法研究了固体火箭发动机中C/C 复合材料表面上颗粒的沉积模态。
该方法考虑了沉积过程中边界形貌的变化对流场的影响，最终得到了颗粒在碳纤维表面的沉积形态。
关键词：LB 方法；流固耦合；C/C 复合材料；沉积0 引言C/C 复合材料是目前新材料领域重点研究和开发的一种新型超高温热结构材料，具有密度小，比强度大、热膨胀系数低、热导率高等特点，是理想的航空航天高温材料[1, 2]。
C/C 复合材料在工作过程中其表面流过的工质为高温燃气。
高温燃气中通常带有燃烧产生的固体颗粒，如选用较高比冲的含铝推进剂时会产生一定量的凝聚相(Al2O3颗粒)。
固体颗粒在C/C 复合材料表面的沉积、冲刷及烧蚀会造成材料内型面的破坏，甚至影响气动性能。
本文使用格子Boltzmann 方法模拟C/C 复合材料中碳纤维上颗粒沉积过程及形态。
格子_Boltzmann方法模拟多孔介质内自然对流蓄热过程
收稿日期:2012－01－18基金项目:江西省科技支撑计划项目(2011BBE50031);江西省科学院引进博士项目(2011－YYB －02);江西省科学院国家级预研项目(2011－YGY －01)格子－Boltzmann 方法模拟多孔介质内自然对流蓄热过程万斌，罗成龙(江西省科学院能源研究所，江西南昌330029)摘要:以在高温热量存储、太阳能利用和建筑节能中有着广泛应用的多孔蓄热装置为研究对象，基于格子－Boltzmann 方法(LBM )的基本原理，建立表征体元尺度(REV )上多孔介质自然对流的热流耦合方程，对多孔介质区域内的定温加热过程进行数值计算，探索了蓄热装置工作效果与多孔介质材料和内部流体特性的关系。
分析获得了多孔介质渗透率和工质热膨胀率增大对多孔介质蓄热的强化作用，以及强自然对流作用下温度场分布所出现的不均匀特性。
关键词:多孔介质;表征体元尺度;格子－Boltzmann ;自然对流中图分类号:O359.1文献标志码:ANumerical Simulation of Heat Storage Process for NaturalConvection in Porous Media with Lattice －Boltzmann MethodWAN Bin ，LUO Cheng-long(Energy Institute of Jiangxi Academy of Sciences ，Nanchang 330029，China )Abstract :In this paper ，heat storage devices with porous media were investigated．These devices were applied in many areas ，such as high-temperature heat storage ，solar energy and HVAC．Based on the basic principles of Lat-tice －Boltzmann Method ，the equations of heat natural convection in porous media was established on REV scale．The process of constant temperature heating in porous media was calculated by LBM．The results of simulation showed that the effect of heat transfer was strengthened by permeability and thermal expansion porous media ，uneven distribution of temperature field was appeared in strong natural convection．Key Words :porous media ;REV ;Lattice －Boltzmann ;natural convection 随着人类对能源的利用和消耗急剧增加能源危机的出现，迫使人们重视节能技术和开发新能源，能源储存和蓄热技术也在此背景下快速发展。
格子-Boltzmann方法模拟霜结晶生长
格子-Boltzmann方法模拟霜结晶生长龚建英;孙金绢;李国君【摘要】为了研究冷壁面上霜结晶的生长过程以及各参数的变化规律,利用格子-Boltzmann方法(LBM)方法,建立一个二维介观模型,实验验证了模型的可靠性,分析了霜层温度变化和密度增长规律,并直观模拟出霜结晶凝聚变化过程.结果表明:霜层表面温度在早期阶段迅速增加,但增加率随着结霜时间增加而减小;霜层内部温度随霜层厚度的增加呈线性增长;霜层平均密度随结霜时间增加呈现出先慢后快的规律;随着结霜过程的进行,由于越往上,霜晶体积分数越小,导致霜层内部密度随霜层厚度的上升而减小.%In order to calculate the frost crystal growth process on a cold fiat surface and study the variation laws of parameters,a two-dimensional mesoscopic model was proposed based on LBM.The predicted results were found to be in good agreement with the data reported in the literatures.This simulation was concerned with temperature and density variation of the frost.The results show that the surface temperature increases rapidly in the early-stage,but the increase rate decreases with frosting.In addition,the interior temperatures rise in a linear curve with increase of the frost layer thickness.The average frost layer density increases with increasing frosting time.It was also found that the nearer the frost layer gets to the cold flat surfaces the higher the frost density is.【期刊名称】《低温工程》【年(卷),期】2013(000)004【总页数】5页(P10-13,28)【关键词】格子-Boltzmann方法;冷壁面;霜结晶【作者】龚建英;孙金绢;李国君【作者单位】西安交通大学能源与动力工程学院热流科学与工程教育部重点实验室西安710049;西安交通大学能源与动力工程学院热流科学与工程教育部重点实验室西安710049;西安交通大学能源与动力工程学院热流科学与工程教育部重点实验室西安710049【正文语种】中文【中图分类】TB611;TK121结霜现象的研究是涉及制冷、空调、航空、航天，农业等多领域中具有重要意义的一项基础性研究。
用格子boltzmann方法模拟mkdv方程
用格子boltzmann方法模拟mkdv方
程
广义双参数MKDV方程是用于二维模拟超深海浪的重要模型，不仅对深海的表
观水动力学提出了更为准确的数学描述，而且能够更为形象地表示不同浪流间的关系。
由于这一模型的复杂程度，往往需要使用计算数值方法来求解，其中最常用的方法就是格子Boltzmann方法（Lattice Boltzmann Method，LBM）。
格子Boltzmann方法是一种基于分子动力学（Molecular Dynamics，MD）思想，以离散时间步长作为微观模型，使用众多拟合子发展而成的数值求解方法。
它采用经典的“扰动—演化”过程，通过将精确且简单的条件和分布函数引入离散的网格空间，将难以直接解决的连续型方程变换为由离散分组定义的有限元算法，从而实现实时的动力学模拟。
与传统的数值法对求解常微分方程的精度较低不同，格子Boltzmann方法可以获得更高的准确性，而且大大降低了求解复杂方程的复杂度。
此外，格子Boltzmann方法还具有计算复杂度低、可并行计算（Parallel computation）、简洁的编程难度、可扩展性强等优点，可以有效解决广义双参数MKDV方程求解的问题。
因此，在互联网安全领域中，格子Boltzmann方法不断发展、应用广泛，是一种有效解决深海表观水动力学机理的技术手段。
顶盖驱动流的格子Boltzmann模拟代码
for(i=1;i<NX;i++) for(j=1;j<NY;j++) for(k=0;k<Q;k++) { ip=i-e[k][0]; jp=j-e[k][1]; F[i][j][k]=f[ip][jp][k]+(feq(k,rho[ip][jp], u[ip][jp])-f[ip][jp][k])/tau_f; } for(i=1;i<NX;i++) for(j=1;j<NY;j++) { u0[i][j][0]=u[i][j][0]; u0[i][j][1]=u[i][j][1]; rho[i][j]=0; u[i][j][0]=0; u[i][j][1]=0; for(k=0;k<Q;k++) { f[i][j][k]=F[i][j][k]; rho[i][j]+=f[i][j][k]; u[i][j][0]+=e[k][0]*f[i][j][k]; u[i][j][1]+=e[k][1]*f[i][j][k]; } u[i][j][0]/=rho[i][j]; u[i][j][1]/=rho[i][j]; } for(j=1;j<NY;j++) for(k=0;k<Q;k++) { rho[NX][j]=rho[NX-1][j]; f[NX][j][k]=feq(k,rho[NX][j],u[NX][j])+(f[NX-1] [j][k]-feq(k,rho[NX-1][j],u[NX-1][j]))*(1-1/tau_f); rho[0][j]=rho[1][j]; f[0][j][k]=feq(k,rho[0][j],u[0][j])+(f[1][j][k] -feq(k,rho[1][j],u[1][j]))*(1-1/tau_f); } for(i=0;i<=NX;i++) for(k=0;k<Q;k++) { rho[i][0]=rho[i][1]; f[i][0][k]=feq(k,rho[i][0],u[i][0])+(f[i][1][k] -feq(k,rho[i][1],u[i][1]))*(1-1/tau_f); rho[i][NY]=rho[i][NY-1]; u[i][NY][0]=U;
变高宽比空腔顶盖驱动流的格子Boltzmann方法模拟
变高宽比空腔顶盖驱动流的格子Boltzmann方法模拟曹先齐;杲东彦;文先太;蔡宁;王亮【摘要】采用格子Boltzmann方法构建空腔顶盖驱动流的数学模型,进行模型验证,并重点分析不同雷诺数Re下高宽比K对腔体流函数分布的影响.结果表明:Re=400,K=0.25时,二级涡由两个变成了三个,除了左下角和右下角的二级涡之外,在一级涡的内部还出现了第三个二级涡,位于腔体的中上部;Re=2000,K=0.5时,腔体出现了两个一级涡,一个位于腔体左侧,一个位于腔体右侧,基本呈对称分布;当腔体高度增加时,腔体内只有两个一级涡,一个一级涡在上方,另一个一级涡在下方,二级涡则未出现;Re=2000时腔体的流函数分布与Re=400时基本相同;腔体宽度对流函数分布的影响大于腔体高度.【期刊名称】《南京工程学院学报（自然科学版）》【年(卷),期】2018(016)003【总页数】6页(P1-6)【关键词】顶盖驱动流;高宽比;格子Boltzmann方法【作者】曹先齐;杲东彦;文先太;蔡宁;王亮【作者单位】南京工程学院能源与动力工程学院,江苏南京211167;南京工程学院能源与动力工程学院,江苏南京211167;南京工程学院能源与动力工程学院,江苏南京211167;南京工程学院能源与动力工程学院,江苏南京211167;南京工程学院能源与动力工程学院,江苏南京211167【正文语种】中文【中图分类】TK121顶盖空腔驱动流的研究在现代飞机、汽车的噪声控制、城市街道峡谷内污染物扩散机理、建筑墙体节能保温等诸多领域中均有着广泛的应用价值.如城市街道峡谷内污染物扩散问题,街道峡谷是指两旁都有连续的高大建筑物的相对狭长的街道空间,是城市的重要组成部分,也是居民活动中最为活跃的场所[1-3]，这些高大建筑物虽然可以在夏天为行人遮阴,但也存在不利的因素,例如阻碍交通污染物的扩散,对人类健康造成严重危害.随着社会的不断发展,城市化进程不断加块,城市大气污染这一问题也越来越严重,影响经济社会的可持续发展[4-6].因此,为了研究城市街道中污染物的分布状况,减缓城市大气污染对人体及周围环境的负面影响,需要深入研究城市街道峡谷内空气的流场特性.流场中蕴含有诸如多级涡、不稳定层流和紊流等复杂物理现象,因此顶盖空腔驱动流的研究成为复杂流体流动流场最为经典的物理模型,也成为验证各类流体流动数值计算方法准确度和计算效率的理想研究对象[7].除控制体积法、有限元法、有限差分法等传统流体流动数值计算方法外,格子Boltzmann方法由于清晰的物理背景、高效的并行效率等优异特征,近年来得到广泛的关注[8].国内外许多学者采用格子Boltzmann 方法对空腔顶盖驱动流进行了数值模拟研究[9-13]：文献[9]模拟了方腔顶盖驱动流,给出了详细的求解过程,提出了一种简化的固壁边界条件处理方法,讨论了雷诺数Re的影响；文献[10]模拟周期性流体驱动的方腔顶盖驱动流,研究了不同雷诺数和振荡频率的影响，结果表明,腔内涡的发展很大程度上依赖于Re和ω,在低Re和ω时,顶盖运动的影响可以有效地传递给腔内流体，然而,在较高的Re和ω时,只有在靠近顶盖狭窄区域的流体会受到顶盖驱动的影响；文献[11]研究了十字形腔体内的双顶盖驱动流问题,考察了雷诺数和摆动频率的影响，研究发现摆动的双壁面比均匀壁面驱动的流体混合得更加均匀；文献[12]模拟具有纳米流体的方腔顶盖驱动流,腔内含有一个薄的水平加热装置，研究表明纳米流体可以增强顶盖驱动流腔体的传热，可以通过升高加热器的位置来减少腔体冷壁面的热损失；文献[13]模拟了具有磁流体等粘塑性流体在空腔内的流动,壁面是非等温的,左壁面为高温,右壁面为低温,考察了雷诺数、Hartmann数、Bingham数、高宽比对传热、流动的影响.上述研究对象大多都是方腔,对非方腔体中驱动流的研究却相对较少.以城市街道峡谷内污染物扩散问题为例,其研究对象大多是非方腔体,而为长方体.本文以变高宽比情况下空腔顶盖驱动流为研究对象,采用格子Boltzmann方法,探讨高宽比对顶盖驱动流流场特性的影响.1 问题描述及控制方程图1 顶盖驱动流示意图顶盖驱动流示意图如图1所示.图1中,腔体的上边界以一个恒定速度u水平移动,其它三个边界则保持静止状态；腔体的水平宽度为L,垂直高度为H,定义高宽比K=H/L.宏观控制方程为：·u=0(1)u=-p+·τ(2)式中：u、ρ、p、τ分别为宏观速度、密度、压力、无量纲松弛时间.2 格子Boltzmann方程模型格子Boltzmann方程为(3)式中:r为空间矢量位置;t为时间;ei为格子离散速度;Δt为时间步长;τ为无量纲松弛时间;为平衡态密度分布函数,其中,wi为权系数,ρ为密度,格子速度c=Δx/Δt,Δx为格子步长.宏观密度和速度公式为(4)(5)式中fi指某一点上i个分量或方向上的分布函数.松弛时间系数公式为(6)式中v为流体动力粘性系数.为了保证数值精度及运行的稳定性,松弛时间系数通常在2≥τ>0.5范围即可.使用平衡态时分布函数赋值给初始分布函数(7)边界条件为:固体壁面无滑移，即ux=uy=0;顶盖速度ux=0.1,uy=0.采用非平衡外推方法[14]进行固体壁面速度边界的处理,其整体精度为2阶,具有较好的数值稳定性.(8)式中:rw为边界上的点；rn为离边界点最近一个内部点.3 模拟结果讨论3.1 模型验证图2 雷诺数Re=400时方腔内流函数图图2为方腔内流函数图,雷诺数Re=400,网格数256×256.从图2中可以看出,腔体内部出现了三个涡:一个大的一级涡和两个小的二级涡.大的一级涡位于腔体的中央,两个二级涡分别位于腔体的左下角和右下角，与文献[7]、[10]、[15]的流函数图基本一致.为了进一步验证模型的准确性,测量方腔中央的一级涡以及左右下角附近的二级涡的中心位置,结果见表1.从表1可以看出,本文程序的模拟结果与文献[7]、[10]、[15]中的涡中心位置的数值基本一致,结果吻合很好.表1 顶盖驱动流的一级涡和二级涡的位置Re=400一级涡左下涡右下涡xyxyxy文献[7]0.556 70.602 90.049 60.046 20.889 10.120 1文献[10]0.556 00.606 90.054 00.047 60.883 70.123 4文献[15]0.566 80.606 90.047 10.048 20.885 70.123 1本文0.558 60.609 10.046 30.045 80.885 20.125 53.2 腔体宽度的影响考察雷诺数Re=400和Re=2 000时腔体宽度对流函数分布的影响,见图3和图4. 从图3可以看出,K=1时,大的一级涡位于腔体的中央,两个二级涡分别位于腔体的左下角和右下角;K=0.5时,一级涡不再是基本对称的圆形,而是横向拉伸的椭圆，并且涡中心的水平位置发生了移动,位于腔体偏右的地方，两个二级涡的位置没变,但是涡的大小发生了变化，左下角的二级涡变大,右下角的二级涡变小；K=0.25时,一级涡继续在水平方向拉伸,涡中心的水平位置继续右移.二级涡由两个变成了三个,除左下角和右下角的二级涡之外,在一级涡的内部还出现了第三个二级涡,位于腔体的中上部.(a) K=1(b) K=0.5(c) K=0.25图3 Re=400时宽度对流函数分布的影响从图4可以看出，当K=1时,腔体流线分布与Re=400时基本一致；K=0.5时,腔体出现了两个一级涡(一个位于腔体左侧,一个位于腔体右侧,基本呈对称分布,只是右侧一级涡稍大一些)和两个二级涡(基本未发生变化)；K=0.25时,两个一级涡都在横向进行了拉伸,横向尺寸变大，二级涡则发生了变化,一个位于右下角,另一个位于腔体中上部,与Re=400时类似，而原先位于左下角的二级涡未出现.(a) K=1(b) K=0.5(c) K=0.25图4 Re=2 000时宽度对流函数分布的影响以上结果表明,随着腔体宽度的增加,顶盖对腔体内流体流动的影响越来越显著. 3.3 腔体高度的影响考察Re=400和Re=2 000时腔体高度对流函数分布的影响,见图5和图6.当Re=400,高度增加时,腔体内只有两个一级涡,一个一级涡在上方,另一个一级涡在下方,二级涡则消失了,腔体的下部只有一些单向的流线.这是由于随着腔体高度的增加,顶盖的移动对腔体下部的影响越来越弱,因此腔体下部只有一个单向的流线,而未形成漩涡.当Re=2 000时,腔体的流函数分布与Re=400时基本相同,腔体仅有上下两个一级涡.有所区别的是,下方的一级涡的位置往左有微弱的偏移.4 结论本文建立了顶盖驱动流的格子Boltzmann方程模型并进行了数值模拟,重点分析了不同雷诺数Re下高宽比K对腔体流函数分布的影响,结论为：1) 腔体内顶盖驱动流的格子Boltzmann方法模拟结果与文献结果一致,模型得到了验证；2) Re=400时,随着腔体宽度的增加(K=0.25),一级涡继续在水平方向拉伸,涡中心的水平位置继续右移，二级涡由两个变成了三个,除了左下角和右下角的二级涡之外,在一级涡的内部还出现了第三个二级涡,位于腔体的中上部，当Re=2 000时,随着腔体宽度的增加(K=0.5),腔体出现了两个一级涡,一个位于腔体左侧,一个位于腔体右侧,基本呈对称分布.当腔体的宽度进一步增加时(K=0.25),二级涡则发生了变化,一个位于右下角,另一个则位于腔体中上部,原先位于左下角的二级涡则未出现.(a) K=1(b) K=2(c) K=4图5 Re=400时高度对流函数分布的影响(a) K=1(b) K=2(c) K=4图6 Re=2 000时高度对流函数分布的影响3)当腔体高度增加时,腔体内只有两个一级涡,一个一级涡在上方,另一个一级涡在下方,二级涡则未出现.Re=2 000时,腔体的流函数分布与Re=400时基本相同.腔体宽度对流函数分布的影响大于腔体高度.参考文献:【相关文献】[1] LIU C H, WONG C C. On the pollutant removal, dispersion, and entrainment over two-dimensional idealized street canyons[J]. Atmospheric Research, 2014, 135-136: 128-142.[2] PARK S J, KIM J J, KIM M J, et al. Characteristics of flow and reactive pollutant dispersion in urban street canyons[J]. Atmospheric Environment, 2015, 108: 20-31.[3] PRASHANT T, SHARAD G. The impact of traffic-flow patterns on airquality in urban street canyons[J]. Environmental Pollution, 2016, 208:161-169.[4] FU X W, LIU J F, BAN-WEISS G A, et al. Effects of canyon geometry on the distribution of traffic-related air pollution in a large urban area: Implications of a multi-canyon air pollution dispersion model[J]. Atmospheric Environment, 2017, 165: 111-121.[5] KUBILAY A, NEOPHYTOU M K, MATSENTIDES S, et al. The pollutant removal capacity of urban street canyons as quantified bythe pollutant exchange velocity[J]. Urban Climate, 2017, 21: 136-153.[6] NOSEK S, KUKACKA L, JURCAKOVA K, et al. Impact of roof height non-uniformity onpollutant transport between astreet canyon and intersections[J]. Environmental Pollution, 2017, 227: 125-138.[7] ARUN S, SATHEESH A. Analysis of flow behaviour in a two sided lid driven cavity using lattice boltzmann technique[J]. Alexandria Engineering Journal, 2015, 54 (4): 795-806. [8] 何雅玲,李庆,王勇,等. 格子Boltzmann方法的工程热物理应用[J]. 科学通报,2009,54(18):2638-2656.[9] 韩善灵,朱平,林忠钦. 基于格子Boltzmann方法模拟方腔顶盖驱动流[J]. 中国机械工程,2005,16(1):64-66, 73.[10] MENDU S S, DAS P K. Fluid flow in a cavity driven by an oscillating lid-A simulation by lattice Boltzmann method[J]. European Journal of Mechanics B/Fluids, 2013(39): 59-70. [11] STHAVISHTHA B R, PERUMAL D A, YADAV A K. Computation of fluid flow in double sided cross-shaped lid-driven cavities using Lattice Boltzmann method[J]. European Journal of Mechanics/B Fluids, 2018(70): 46-72.[12] DAHANI Y, HASNAOUI M, AMAHMID A, et al. Lattice-Boltzmann modeling of forced convection in a lid-driven square cavity filled with a nanofluid and containing a horizontal thin heater[J]. Energy Procedia, 2017, 139: 134-139.[13] KEFAYATI G R, TANG H. MHD mixed convection of viscoplastic fluids in different aspect ratios of a lid-driven cavity using LBM[J]. International Journal of Heat and Mass Transfer, 2018, 124: 344-367.[14] 郭照立,郑楚光. 格子Boltzmann方法的原理及应用[M]. 北京:科学出版社,2009.[15] 何雅玲,王勇,李庆. 格子Boltzmann方法的理论及应用[M]. 北京:科学出版社,2009.。
Matlab实现格子玻尔兹曼方法(Lattice Boltzmann Method,LBM)模拟
Matlab实现格子玻尔兹曼方法(Lattice Boltzmann Method,LBM)模拟%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% cylinder.m: Flow around a cyliner, using LBM %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% This program is free software; you can redistribute it and/or% modify it under the terms of the GNU General Public License% as published by the Free Software Foundation; either version 2% of the License, or (at your option) any later version.% This program is distributed in the hope that it will be useful,% but WITHOUT ANY WARRANTY; without even the implied warranty of % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the % GNU General Public License for more details.% You should have received a copy of the GNU General Public% License along with this program; if not, write to the Free% Software Foundation, Inc., 51 Franklin Street, Fifth Floor,% Boston, MA 02110-1301, USA. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clear% GENERAL FLOW CONSTANTSlx = 250;ly = 51;obst_x = lx/5+1; % position of the cylinder; (exactobst_y = ly/2+1; % y-symmetry is avoided)obst_r = ly/10+1; % radius of the cylinderuMax = 0.02; % maximum velocity of Poiseuille inflowRe = 100; % Reynolds numbernu = uMax * 2.*obst_r / Re; % kinematic viscosityomega = 1. / (3*nu+1./2.); % relaxation parametermaxT = 400000; % total number of iterationstPlot = 5; % cycles% D2Q9 LATTICE CONSTANTSt = [4/9, 1/9,1/9,1/9,1/9, 1/36,1/36,1/36,1/36];cx = [ 0, 1, 0, -1, 0, 1, -1, -1, 1];cy = [ 0, 0, 1, 0, -1, 1, 1, -1, -1];opp = [ 1, 4, 5, 2, 3, 8, 9, 6, 7];col = [2:(ly-1)];[y,x] = meshgrid(1:ly,1:lx);obst = (x-obst_x).^2 + (y-obst_y).^2 <= obst_r.^2;obst(:,[1,ly]) = 1;bbRegion = find(obst);% INITIAL CONDITION: (rho=0, u=0) ==> fIn(i) = t(i)fIn = reshape( t' * ones(1,lx*ly), 9, lx, ly);% MAIN LOOP (TIME CYCLES)for cycle = 1:maxT% MACROSCOPIC VARIABLESrho = sum(fIn);ux = reshape ( ...(cx * reshape(fIn,9,lx*ly)), 1,lx,ly) ./rho;uy = reshape ( ...(cy * reshape(fIn,9,lx*ly)), 1,lx,ly) ./rho;% MACROSCOPIC (DIRICHLET) BOUNDARY CONDITIONS% Inlet: Poiseuille profileL = ly-2; y = col-1.5;ux(:,1,col) = 4 * uMax / (L*L) * (y.*L-y.*y);uy(:,1,col) = 0;rho(:,1,col) = 1 ./ (1-ux(:,1,col)) .* ( ...sum(fIn([1,3,5],1,col)) + ...2*sum(fIn([4,7,8],1,col)) );% Outlet: Zero gradient on rho/uxrho(:,lx,col) = rho(:,lx-1,col);uy(:,lx,col) = 0;ux(:,lx,col) = ux(:,lx-1,col);% COLLISION STEPfor i=1:9cu = 3*(cx(i)*ux+cy(i)*uy);fEq(i,:,:) = rho .* t(i) .* ...( 1 + cu + 1/2*(cu.*cu) ...- 3/2*(ux.^2+uy.^2) );fOut(i,:,:) = fIn(i,:,:) - ...omega .* (fIn(i,:,:)-fEq(i,:,:));end% MICROSCOPIC BOUNDARY CONDITIONSfor i=1:9% Left boundaryfOut(i,1,col) = fEq(i,1,col) + ...18*t(i)*cx(i)*cy(i)* ( fIn(8,1,col) - ...fIn(7,1,col)-fEq(8,1,col)+fEq(7,1,col) );% Right boundaryfOut(i,lx,col) = fEq(i,lx,col) + ...18*t(i)*cx(i)*cy(i)* ( fIn(6,lx,col) - ...fIn(9,lx,col)-fEq(6,lx,col)+fEq(9,lx,col) );% Bounce back regionfOut(i,bbRegion) = fIn(opp(i),bbRegion);end% STREAMING STEPfor i=1:9fIn(i,:,:) = ...circshift(fOut(i,:,:), [0,cx(i),cy(i)]);end% VISUALIZATIONif (mod(cycle,tPlot)==0)u = reshape(sqrt(ux.^2+uy.^2),lx,ly);u(bbRegion) = nan;imagesc(u');axis equal off; drawnowendend。
基于格子Boltzmann方法模拟方腔顶盖驱动流
基于格子Boltzmann方法模拟方腔顶盖驱动流
韩善灵;朱平;林忠钦
【期刊名称】《中国机械工程》
【年(卷),期】2005(016)001
【摘要】格子Boltzmann方法是使用时间和空间完全离散的细观模型来模拟流体宏观行为的一种新方法,模型的平均行为符合宏观的Navier-Stokes方程.给出了格子Boltzmann方法详细的求解过程,提出了一种简化的固壁边界条件处理方法.用格子Boltzmann方法对方腔顶盖驱动流进行了数值模拟,并对模拟结果进行了分析和讨论.通过与前人已有的试验及数值研究结果对比,验证了格子Boltzmann方法的正确性.
【总页数】4页(P64-66,73)
【作者】韩善灵;朱平;林忠钦
【作者单位】上海交通大学,上海,200030;上海交通大学,上海,200030;上海交通大学,上海,200030
【正文语种】中文
【中图分类】U462
【相关文献】
1.格子Boltzmann方程模拟高雷诺数三维方腔流 [J], 康海贵;李承功
2.变高宽比空腔顶盖驱动流的格子Boltzmann方法模拟 [J], 曹先齐;杲东彦;文先太;蔡宁;王亮
3.基于格子Boltzmann方法的倾斜方腔自然对流模拟 [J], 朱建奇;侯佳煜;杲东彦;林福建;陈玮玮;鹿世化
4.基于时间的驱动方腔流的高精度多重网格方法数值模拟 [J], 葛永斌;田振夫;吴文权
5.用Gao-Yong湍流方程组数值模拟高雷诺数顶盖驱动方腔流 [J], 闫文辉;张常贤;陈宁宁;高歌
因版权原因，仅展示原文概要，查看原文内容请购买。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
if(n%1000==0) output(n); if(error<1.0e-6) break; } } } return 0; } void init() { dx=1.0; dy=1.0; Lx=dx*double(NY); Ly=dy*double(NX); dt=dx; c=dx/dt; rho0=1.0; Re=100; niu=U*Lx/Re; tau_f=3.0*niu+0.5; std::cout<<"tau_f= "<<tau_f<<endl; for(i=0;i<=NX;i++) for(j=0;j<=NY;j++) { u[i][j][0]=0; u[i][j][1]=0; rho[i][j]=rho0; u[i][NY][0]=U; for(k=0;k<Q;k++) { f[i][j][k]=feq(k,rho[i][j],u[i][j]); } } } double feq(int k,double rho,double u[2]) { double eu,uv,feq; eu=(e[k][0]*u[0]+e[k][1]*u[1]); uv=(u[0]*u[0]+u[1]*u[1]); feq=w[k]*rho*(1.0+3.0*eu+4.5*eu*eu-1.5*uv); return feq; } void evolution() {
for(i=1;i<NX;i++) for(j=1;j<NY;j++) for(k=0;k<Q;k++) { ip=i-e[k][0]; jp=j-e[k][1]; F[i][j][k]=f[ip][jp][k]+(feq(k,rho[ip][jp], u[ip][jp])-f[ip][jp][k])/tau_f; } for(i=1;i<NX;i++) for(j=1;j<NY;j++) { u0[i][j][0]=u[i][j][0]; u0[i][j][1]=u[i][j][1]; rho[i][j]=0; u[i][j][0]=0; u[i][j][1]=0; for(k=0;k<Q;k++) { f[i][j][k]=F[i][j][k]; rho[i][j]+=f[i][j][k]; u[i][j][0]+=e[k][0]*f[i][j][k]; u[i][j][1]+=e[k][1]*f[i][j][k]; } u[i][j][0]/=rho[i][j]; u[i][j][1]/=rho[i][j]; } for(j=1;j<NY;j++) for(k=0;k<Q;k++) { rho[NX][j]=rho[NX-1][j]; f[NX][j][k]=feq(k,rho[NX][j],u[NX][j])+(f[NX-1] [j][k]-feq(k,rho[NX-1][j],u[NX-1][j]))*(1-1/tau_f); rho[0][j]=rho[1][j]; f[0][j][k]=feq(k,rho[0][j],u[0][j])+(f[1][j][k] -feq(k,rho[1][j],u[1][j]))*(1-1/tau_f); } for(i=0;i<=NX;i++) for(k=0;k<Q;k++) { rho[i][0]=rho[i][1]; f[i][0][k]=feq(k,rho[i][0],u[i][0])+(f[i][1][k] -feq(k,rho[i][1],u[i][1]))*(1-1/tau_f); rho[i][NY]=rho[i][NY-1]; u[i][NY][0]=U;
f[i][NY][k]=feq(k,rho[i][NY],u[i][NY])+(f[i][NY-1] [k]-feq(k,rho[i][NY-1],u[i][NY-1]))*(1-1/tau_f); } } void output(int m) { ostringstream name; name<<"cavity_100_"<<m<<".dat"; ofstream out(name.str().c_str()); out<<"Title= \"LBM Lid Driven Flow\"\n"<<"VARIABLES=\"X\",\"Y\",\"U\",\"V\"\n"<<"ZONE T=\"BOX\",I=" <<NX+1<<",J="<<NY+1<<",F=POINT"<<endl; for(j=0;j<=NY;j++) for(i=0;i<=NX;i++) { out<<double(i)/Lx<<" " <<double(j)/Ly<<" "<<u[i][j][0]<<" "<< u[i][j][1]<<endl; } } void Error() { double temp1,temp2; temp1=0; temp2=0; for(i=1;i<NX;i++) for(j=1;j<NY;j++) { temp1 +=( (u[i][j][0]-u0[i][j][0])*(u[i][j][0]-u0[i][j][0]) +(u[i][j][1]-u0[i][j][1])*(u[i][j][1]-u0[i][j][1])); temp2 +=(u[i][j][0]*u[i][j][0]+u[i][j][1]*u[i][j][1]); } temp1=sqrt(temp1); temp2=sqrt(te0); }