下载文档原格式
基于格子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］运用该方法求解含有黏性的声传播问题。
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方程，通过将连续的物理系统离散化为网格系统进行模拟。
网格系统中的每个格子代表一个微观粒子的状态，而碰撞、传输和外部力的作用通过计算和更新这些格子的状态来实现。
该方法主要包含两个步骤：碰撞和传输。
在碰撞过程中，格子中的粒子通过相互作用和碰撞来改变其速度和方向，从而模拟了分子之间的碰撞过程。
在传输过程中，碰撞后的粒子根据流体的速度场进行移动，从而模拟了背景流场对粒子运动的影响。
其次，在格子玻尔兹曼方法中，模拟的过程可以简化为两个部分：演化和碰撞。
在每个时间步长内，系统首先根据粒子速度和位置的信息计算出相应格点上的分布函数，然后通过碰撞步骤更新这些分布函数以模拟粒子之间的碰撞效应。
通过迭代演化和碰撞步骤，系统的宏观行为可以得到。
格子玻尔兹曼方法中最常用的碰撞操作是BGK碰撞算子，它根据粒子的速度和位置信息计算出新的分布函数，并用该新分布函数代替原来的分布函数。
而在传输过程中，粒子通过碰撞后得到的新速度和方向进行移动。
最后，格子玻尔兹曼方法在流体力学领域具有广泛的应用，特别是在稀薄气体流动、微纳尺度流动和多相流等问题中。
由于其适用于模拟分子尺度和介观尺度流动问题，因此在利用普通的Navier-Stokes方程难以模拟的问题中表现出了良好的效果。
此外，格子玻尔兹曼方法还可以用于模拟流动中的热传导问题、气体分子在多孔介质中的传输问题以及颗粒与流体相互作用等多种复杂流动现象。
近年来，随着计算机性能的不断提高，格子玻尔兹曼方法也得到了快速发展，在模拟大规模真实流体问题方面取得了不错的结果。
总结来说，格子玻尔兹曼方法通过将连续的物理系统离散化为网格系统，模拟粒子碰撞和传输过程，实现了对流体力学问题的数值模拟。
!此?程ì序ò用?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方法（LBM）是一种计算流体机械系统的模拟方法，它允许将多物理学过程直接耦合在一起，同时也满足计算机能力的限制。
在过去的几十年里，格子Boltzmann方法已经成功应用于各种复杂系统中，并且得到了广泛的应用。
在本文中，我们将重点关注格子Boltzmann方法模拟流动沸腾过程。
流动沸腾过程是指流体作为一个定常动态系统，当流体间的温度波动超过沸点时，出现蒸发和汽化的共同现象。
在实际的工程应用中，这种现象可以被用来模拟大气中的蒸发和汽化过程，也可以用来模拟熔炼过程中的熔融气体的性质，当物体在极端条件下的表面上出现沸腾时，也可以用来模拟这种情况。
格子Boltzmann方法是一种经典的集总体方法，通常是将空间区域划分为多个格子大小的单元格，然后模拟每个单元格内流体粒子的运动方式。
首先，根据Boltzmann方程，定义系统中每个单元格内流体粒子的分布函数，然后根据撞击原理，建立系统的碰撞模型，最后，利用快速Fourier变换（FFT）计算出每个单元格的空间和时间变化，最终得到系统中流体粒子分布函数的空间和时间变化，从而模拟出系统的动力学过程。
格子Boltzmann方法模拟流动沸腾过程的优势在于它可以解决复杂的动力学问题，在流动沸腾过程中，可以使用格子Boltzmann方法考虑到热传导和蒸汽压力的影响，同时考虑到蒸发和汽化等热力学效应，从而得到较准确的结果。
此外，采用格子Boltzmann方法模拟流动沸腾过程，具有计算效率高、可靠性高、精度高的优势。
由于上述优点，格子Boltzmann方法在模拟流动沸腾过程方面已经取得了许多成功的应用。
例如，Zhao等对水汽-气体混合流动的二相沸腾过程进行了模拟，他们利用格子Boltzmann方法模拟了水汽-气体混合的沸腾过程，比较了侧向温度波及其在慢一点的过程中的水汽沉积速率，并表明其模拟结果良好，与实际的实验结果相一致。
枝晶和气孔演化规律的格子boltzmann方法
数值模拟
格子Boltzmann方法是一种新兴的、用于模拟二维和三维
物质在复杂流变行为中的方法，在过去几年中，由于晶体和气
孔演化机制已经被深入研究，格子Boltzmann方法也受到越来
越多的关注。
首先，格子Boltzmann方法通过对一系列物理量进行条件
分布来复杂模拟晶体和气孔演化。
这一方法可以为建模和模拟
气孔和晶体形成机制提供有益依据，它能够模拟系统的温度依
赖性和宏观动力学行为的变化。
同时，格子Boltzmann方法也
能够有效地提取气孔演化和晶体结构演化的行为特征。
此外，格子Boltzmann方法可以将复杂的物理过程抽象化，从而提取出系统中比较容易理解的量化模型。
基于定属性原理，该方法可以用来模拟数据的动态变化，从而获得晶体的形貌、
尺寸和气孔的分布状况等信息，为研究晶体和气孔演化提供了
重要的参考资料。
最后，格子Boltzmann方法可以有效地解决抽象的高维物
理系统问题，其优势在于能够精准地模拟晶体和气孔演化的行为，它已经广泛应用于材料科学、工程学和生物学等领域，成
为现代互联网科学领域中解决复杂过程问题的基石。
顶盖流热力学模拟程序（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。
