计算流体驱动方腔程序

lbm 方腔自然对流程序自然对流是指在没有外力驱动的情况下，由于温度差异引起的流体运动。

在生活中，我们常常会遇到一些与自然对流有关的现象，比如水壶中的水会因热胀冷缩而产生热气上升的现象，热空气会从窗户朝外流出等等。

而在工程领域，我们需要对自然对流进行深入研究，以便更好地理解和应用这一现象。

为了更准确地模拟和预测自然对流的行为，科学家们开发了许多数学模型和计算方法。

其中，LB（Lattice Boltzmann）方法是一种非常常用的数值计算方法，特别适用于模拟复杂的流体运动。

LB方法将流体看作是由大量微小粒子组成的，通过在空间中构建一个网格来模拟流体的运动。

在这个网格中，每个微小粒子都有一定的速度和密度，并通过碰撞和散射的过程来描述流体的运动。

在研究自然对流的过程中，LB方法可以很好地模拟流体在容器内的温度分布、速度场以及压力场等。

通过调节初始条件和边界条件，我们可以模拟不同温度差异下的自然对流现象。

这对于很多工程问题来说具有非常重要的意义。

比如在建筑物的设计中，我们需要考虑到自然对流的影响，以便优化空调系统的设计，提高能源利用效率。

此外，在电子设备散热和核能工程等领域，自然对流也是一个重要的研究课题。

在实际应用中，LB方法的计算结果与实验结果和其他数值方法的结果进行对比，可以发现其准确性和可靠性。

同时，LB方法还具有较少的计算资源占用和较快的计算速度的优点。

这使得它在工程领域的应用非常广泛。

总之，LB方法是研究自然对流的重要工具之一。

通过模拟自然对流现象，我们可以更好地理解和预测流体的行为，并在工程实践中应用这些知识。

未来，随着计算技术的发展，LB方法还将进一步改进和应用于更多领域，从而促进工程科学的发展。

附录F 二维不可压缩黏性流体方腔流动问题的有限体积算法与计算程序二维方腔流动问题是一个不可压缩黏性流动中典型流动。

虽然目前尚不能求得它的解析解，但是它常被用来作为检验各种数值算法计算精度和可靠性的算例。

文献中几乎大多数算法都对它进行过计算。

在本算例中采用有限体积算法三阶迎风型QUICK 离散格式对它进行数值求解。

同时，为了初学者入门和练习方便,这里给出了用C 语言和Fortran77语言编写的计算二维不可压缩黏性方腔流动问题计算程序，供大家学习参考。

F-1利用有限体积算法三阶迎风型QUICK 离散格式求解二维不可压缩黏性流体方腔流动问题1.二维不可压缩黏性流体方腔流动问题二维不可压缩黏性流体方腔流动（cavity flow ）：有一正方形腔室，其量纲为一的宽度为1.0，里面充满静止的不可压缩黏性流体，方腔内初始时刻压力和密度为=1.0=1.0p ρ、它周围壁面（左右壁面和底面）固定不动，上壁面以量纲为一的速度1.0沿着上壁面方向自左向右运动（图F.1）。

2. 基本方程组、初始条件和边界条件设流体是黏性流体。

二维方腔流动问题在数学上可以由二维不可压缩黏性流动N - S()()u v x y x x φρφρφμ∂∂∂∂⎛⎫+= ⎪∂∂∂∂⎝⎭⎝⎭） 其中u 为变量φ在水平x 方向的流速，v 为φ在垂直y 方向的流速，μ为黏度，S φ为源项。

源项中不仅包含压力梯度项，也包含时间导数项。

初始条件：方腔上壁面以量纲为一的速度1.0沿着上壁面方向自左向右运动。

图F.1二维不可压缩黏性方腔流问题示意图边界条件：流动速度u v 、均可采用无滑移边界条件，压力p 采用自由输出边界条件。

3.计算网格划分和控制体单元与节点定义采用交错网格，图F.2和图F.3是计算网格、控制体单元和节点示意图。

节点P 所在主控制单元如图F.2中有阴影部分所示。

在x 方向与节点P 相邻的节点为W 和E ，在y 方向与节点P 相邻的节点为S 和N ，主控制单元界面分别为w s e n 、、、。

涡量流函数法中心差分格式的二维方腔顶盖驱动计算本题是关于粘性流体方腔顶盖驱动的问题。

采用涡量流函数法，在均分网格下，用中心差分格式进行计算，结果与文献中所采用的其他方法和格式进行比较，认为中心差分格式符合计算的精度，但是其明显缺点是计算过程不稳定。

1、参数的无量纲化令 X =x H; Y =y H; U =u u t ; ⁄⁄⁄V =v u t ; ⁄Ω=w w 0 ; ⁄Ψ=φφ0⁄ 其中，由涡量的定义式：u v y xω∂∂=-∂∂ 将速度的导数项无量纲化，并将上述定义式代入，即ω=∂u ∂y −∂v ∂x =u t ∂U H ∂Y −u t ∂V H ∂X =u t H (∂U ∂Y −∂V ∂X)令ut H ⁄=w 0 则涡量的无量纲形式可以写为： Ω=∂U ∂Y−∂V ∂X=ωω0再由流函数的定义式：u =∂ψ/∂y 同理，令 ψ0=u t ∗H可以得到流函数的无量纲表达式为：Ψ=ψψ0流函数边界条件(X,Y )∈AB Ψ=0(X,Y )∈BC Ψ=0 (X,Y )∈CD Ψ=0 (X,Y )∈AD Ψ=0涡量的边界条件（Thom 公式）(X,Y )∈AB Ω1，j =2(Ψ2,j −Ψ1,j )δX 2 (X,Y )∈BC Ωi ，1=2(Ψi,2−Ψi,1)δY 2 (X,Y )∈CD ΩL1,j =2(ΨL2,j −ΨL1,j )δX 2 (X,Y )∈AD Ωi ，M1=2(Ψi,M2−Ψi,M1)δY 2+ 2δY2、方程的无量纲处理过程2.1、非守恒型方程的处理（1）将以上假设的各式代入到非守恒型方程组：2222()u v x y x y ωωωωρρμ∂∂∂∂+=+∂∂∂∂22220x y ψψω∂∂+-=∂∂得到以下无量纲形式的方程： ①涡量控制方程∂2Ω∂X 2+∂2Ω∂Y 2=R e (U ∂Ω∂X +V ∂Ω∂Y) ②流函数控制方程∂2Ψ∂X 2+∂2Ψ∂Y 2−Ω=0 （2）采用有限差分法离散非守恒型涡量的无量纲控制方程：R e (U P ΩE −ΩW 2ΔX +V P ΩN −ΩS 2ΔY )=ΩW −2ΩP +ΩE ∆X 2+ΩS −2ΩP +ΩN∆Y 2整理得到：(2∆X 2+2∆Y 2)ΩP =(1∆X 2+R e U P2ΔX)ΩW +(1∆X 2−R e U P 2ΔX )ΩE+(1∆Y 2+R e V P 2ΔY)ΩS +(1∆Y 2−R e V P2ΔY)ΩN简化为：4ΩP =(1+R e U P ΔX2)ΩW +(1−R e U P ΔX2)ΩE+(1+R e V P ΔY2)ΩS +(1−R e V P ΔY2)ΩN写成一般形式是: a p Ωp =a W ΩW +a E ΩE +a S ΩS +a N ΩN 其中 a p =4 ；a W =(1+R e U P ΔX2) ；a E =(1−R e U P ΔX2)a S =(1+R e V P ΔY2) ； a N =(1−R e V P ΔY2)2.2、守恒型方程的处理（1）将假设的各个无量纲量代入守恒型涡量方程 得到其无量纲形式为：Re[∂(UΩ)∂X +∂(VΩ)∂Y ]=∂2Ω∂X 2+∂2Ω∂Y2（2）采用有限容积法离散守恒型涡量的无量纲控制方程：对于扩散项，有∫∫[ ∂∂X ewn s (∂Ω∂X )+∂∂Y (∂Ω∂Y )]dXdY=∂Ω∂X |w e ∆Y +∂Ω∂Y |sn∆X =[ΩE −ΩP (δX )e −ΩP −Ωw (δX)w ]∆Y +[ΩN −ΩP (δY )n −ΩP −ΩS(δY)s]∆X对于对流项，有∫∫[∂(UΩ)∂X +∂(VΩ)∂Y ]ewn s dXdY =[(UΩ)e −(UΩ)w ]∆Y +[(VΩ)n −(VΩ)s ] ∆X对于界面上的取值，采用如下形式：U e ≥0 时，Ωe =ΩP +(f e +−ΩP )； U e ≤0 时，Ωe =ΩE +(f e −−ΩE ) U w ≥0时，Ωw =ΩW +(f e +−ΩW )； U w ≤0 时，Ωw =ΩP +(f e −−ΩP ) V n ≥0时，Ωn =ΩP +(f e +−ΩP )； V n ≤0 时，Ωn =ΩN +(f e −−ΩN ) V s ≥0时，Ωs =ΩS +(f e +−ΩS )； V s ≤0 时，Ωs =ΩP +(f e −−ΩP )其中，界面上涡量的插值采用中心差分格式，令f e +=f e −=ΩP +ΩE2 ；f w +=f w −=ΩP +Ωw2；f n +=f n −=ΩP +ΩN2；f s +=f s −=ΩP +ΩS2所以，Ωe =ΩP +ΩE2Ωw =ΩP +ΩW2Ωn =ΩP +ΩN2Ωs =ΩP +ΩS2关于界面上速度的插值，可以用上一层次的流函数来表示，形式如下： U e =ΨN +ΨNE −ΨS −ΨSE4∆Y U w =ΨNW +ΨN −ΨS −ΨSW4∆YV n =ΨNW +ΨW −ΨE −ΨNE 4∆XV s =ΨW +ΨSW −ΨSE −ΨE4∆X由以上各式，整理得到守恒型涡量控制方程无量纲形式的离散方程为：R e [ΨN +ΨNE −ΨS −ΨSE 4∆YΩP +ΩE2−ΨNW +ΨN −ΨS −ΨSW 4∆YΩP +ΩW2]∆Y+R e [ΨNW +ΨW −ΨE −ΨNE 4∆XΩP +ΩN2−ΨW +ΨSW −ΨSE −ΨE 4∆XΩP +ΩS2] ∆X= [ΩN −ΩP(δY )n−ΩP −ΩS (δY)s]∆X +[ΩE −ΩP (δX )e−ΩP −Ωw (δX)w]∆Y综上，可整理得到（同位的均分网格）：4ΩP = (1−R e 8(ΨN +ΨNE −ΨS −ΨSE ))ΩE +(1+R e 8(ΨNW +ΨN −ΨS −ΨSW ))ΩW+(1−R e8(ΨNW+ΨW−ΨE−ΨNE))ΩN+(1+R e8(ΨW+ΨSW−ΨSE−ΨE))ΩS写成一般形式为：b pΩp=b EΩE+b WΩW+b NΩN+b SΩS 其中，b p=4b E=1−R e8(ΨN+ΨNE−ΨS−ΨSE)b W=1+R e8(ΨNW+ΨN−ΨS−ΨSW)b N=1−R e8(ΨNW+ΨW−ΨE−ΨNE)b S=1+R e8(ΨW+ΨSW−ΨSE−ΨE)2.3、对于流函数方程，离散可得到：(2∆X2+2∆Y2)ΨP=1∆X2ΨW+1∆X2ΨE+1∆Y2ΨS+1∆Y2ΨN−ΩP整理并写成一般形式是：c pΨp=c WΨW+c EΨE+c SΨS+c NΨN−d其中，c p=4 ；c W=1 ；c E=1 ；c S=1 ；c N=1；d=ΩP∙∆X2并且，采用有限容积法离散的结果与采用这里的方法得到的结果是相同的。

第1篇1. 方腔流动概述方腔流动是指在一个封闭的矩形腔体中，由于外部驱动力的作用，腔体内流体产生的流动现象。

常见的驱动方式有顶盖驱动、双边驱动等。

2. 流动特性分析2.1 顶盖驱动方腔流动- 雷诺数（Re）：方腔尺寸为1m×1m，顶盖以1m/s的速度向右移动，流体粘性系数为0.001m²/s，因此Re = UL/ν = 1000，属于过渡流区域。

- 网格划分：采用100×100的网格，贴壁第一层网格高度为1.3mm，保证壁面处y 足够小。

- 求解器：使用OpenFOAM中的icoFoam和pisoFoam求解器，分别对应层流和湍流模型。

- 结果分析：两种求解器获得的最终流场完全相同，说明在过渡流区域，层流和湍流模型对结果的影响不大。

2.2 双边反向驱动方腔流动- 研究目的：探究双边驱动方腔内流流场的过渡流临界特性，捕捉各种流动分岔点，分析其对流场特性带来的改变。

- 方法：基于格子玻尔兹曼方法的数值模拟，计算各流动状态发生变化时的临界雷诺数，绘制Hopf和Neimark-Sacker流动分岔点以及湍流临界点随速度比的函数图像。

- 结果分析：- 双边驱动内流流场的稳定性较强，说明第二条边的驱动条件可以有效提高流场的稳定性。

- 流场稳定性最强时，双边驱动条件相同时可以更好地提高流场稳定性。

3. 数值模拟方法3.1 OpenFOAM- OpenFOAM是一款开源的CFD（计算流体力学）软件，广泛应用于各种流体流动问题的模拟。

- icoFoam和pisoFoam是OpenFOAM中的不可压缩求解器，分别适用于层流和湍流模拟。

3.2 格子玻尔兹曼方法- 格子玻尔兹曼方法是一种基于统计物理的数值模拟方法，适用于模拟复杂流体流动问题。

- 该方法具有计算效率高、并行性好等优点。

4. 总结方腔流动是一个典型的流体力学问题，通过数值模拟方法可以研究其流动特性。

本文分析了顶盖驱动和双边反向驱动方腔流动的数值模拟结果，揭示了流场稳定性和过渡流临界特性等方面的规律。

一、问题描述方腔顶盖驱动流动如图1所示的一个简化两维方腔（高，宽都等于L)，内部充满水分。

上表面为移动墙，非维化速度为ｕ／ｕ0 =1。

其他三面为固定墙。

试求方腔内水分流动状态。

u=1, v=0u=0, v=0 u=0,v=0u=0, v=0图1常微分方程理论只能求解极少一类常微分方程；实际中给定的问题不一定是解析表达式，而是函数表，无法用解析解法.二、离散格式数值解法:求解所有的常微分方程 计算解函数 y(x) 在一系列节点a = x 0< x 1<…<x n = b 处的近似值),...,1()(n i x y y i i =≈节点间距为步长，通常采用等距节点，即取 hi = h (常数)。

步进式：根据已知的或已求出的节点上的函数值计算当前节点上的函数值，一步一步向前推进。

因此只需建立由已知的或已求出的节点上的函数值求当前节点函数值的递推公式即可。

欧拉方法1(,) 0,1,...n n n n y y h f x y n +=+=几何意义在假设 y n = y (x n )，即第 n 步计算是精确的前提下，考虑公式或方法本身带来的误差: R n = y (x n +1) y n +1 , 称为局部截断误差.截断误差: 实际上，y (x n ) ? y n ， y n 也有误差，它对y n +1的误差也有影响，见下图。

但这里不考虑此误差的影响，仅考虑方法或公式本身带来的误差，因此称为方法误差或截断误差。

局部截断误差的分析：由于假设y n = y (x n ) ，即y n 准确，因此分析局部截断误差时将y (x n +1) 和 y n +1都用点x n 上的信息来表示，工具：Taylor 展开。

显式欧拉公式一阶向前差商近似一阶导数推导如下：223111232()[()()()()][ (,)] ()()h n n n n n n n n n h n R y x y y x hy x y x O h y hf x y y x O h +++'''=-=+++-+''=+1()()()n n n y x y x y x h+-'≈111()()() ()()(,)n n n n nn n n n n y x y x hy x y x y y x y y h f x y +++'≈+↑≈≈=+隐式欧拉公式xn +1点向后差商近似导数 推导如下：几何意义设已知曲线上一点 P n (x n , y n ),过该点作弦线，斜率为(x n +1 , y n +1 ) 点的方向场f (x ,y )方向,若步长h 充分小，可用弦线和垂线x =x n +1的交点近似曲线与垂线的交点。

一、二、问题描述方腔顶盖驱动流动如图1所示的一个简化两维方腔（高，宽都等于L)，内部充满水分。

上表面为移动墙，非维化速度为ｕ／ｕ0 =1。

其他三面为固定墙。

试求方腔内水分流动状态。

u=1, v=0u=0, v=0u=0,v=0u=0, v=0图1常微分方程理论只能求解极少一类常微分方程；实际中给定的问题不一定是解析表达式，而是函数表，无法用解析解法.二、离散格式数值解法:求解所有的常微分方程 计算解函数 y(x) 在一系列节点a = x 0< x 1<…<x n = b 处的近似值),...,1()(n i x y y i i =≈节点间距为步长，通常采用等距节点，即取 hi = h (常数)。

步进式：根据已知的或已求出的节点上的函数值计算当前节点上的函数值，一步一步向前推进。

因此只需建立由已知的或已求出的节点上的函数值求当前节点函数值的递推公式即可。

欧拉方法1(,) 0,1,...n n n n y y h f x y n +=+=几何意义在假设 y n = y (x n )，即第 n 步计算是精确的前提下，考虑公式或方法本身带来的误差: R n = y (x n +1)y n +1 , 称为局部截断误差.显式欧拉公式一阶向前差商近似一阶导数223111232()[()()()()][ (,)] ()()h n n n n n n n n n h n R y x y y x hy x y x O h y hf x y y x O h +++'''=-=+++-+''=+推导如下：隐式欧拉公式x n +1点向后差商近似导数 推导如下：1()()()n n n y x y x y x h+-'≈111()()() ()()(,)n n n n nn n n n n y x y x hy x y x y y x y y h f x y +++'≈+↑≈≈=+11()()()n n n y x y x y x h++-'≈11()()()()n n n n ny x y x hy x y x y ++'≈+↑≈几何意义设已知曲线上一点 P n (x n , y n ),过该点作弦线，斜率为(x n +1 , y n +1 ) 点的方向场f (x ,y )方向,若步长h 充分小，可用弦线和垂线x =x n +1的交点近似曲线与垂线的交点。