计算流体大作业-基于C++的顶部驱动流模拟

（3）
4 3 u2 9 (1 2 c 2 ) i=0 C u 9 (Ci u )2 3 u 2 1 Fi (1 3 i 2 ) i=1,2,3,4 c 2 c4 2 c2 9 1 Ci u 9 (Ci u )2 3 u 2 ) i=5,6,7,8 (1 3 2 c 2 c4 2 c2 36
语言编程，求解了顶部驱动流模型，运用 tecplot360 对所得数据进行处理，所 得流线图显示了该流动模型的典型现象，验证了方法的正确性。
关键词：C/C++ 顶部驱动流 D2Q9 格子波尔兹曼方法 1.引言
格子气自动机(lattice gas automata, LGA)以及从它发展而来的格子 Boltzmann 方法(lattice boltzmann method, LBM)是使用简单的微观模型来模拟流体的宏观行 为的一种新方法。LBM 的基本理论来源于格子气自动机 LBM 与以宏观连续方程的离散 化为基础的传统数值方法不同，其基础是微观模型与微观运动论方法。它保留了分子 运动学的许多优点，如物理条件清晰，边界条件易于实现等。LBM 的基本思想是对系 统中的微观粒子给出一个简化的动力模型，然后通过统计平均，得到宏观物理量 (如 密度、速度)的流体动力学方程。 顶部驱动流是非常传统的流动问题， 很多学者都对其进行过研究。 本文借助 D2Q9 格子波尔兹曼方法，在 Visual Studio 2012 平台上运用 C/C++语言编程，对上述传统 问题进行了求解，并借助 Tecplot360 工具，对所得数据进行后处理。
7.结论
本文借助 D2Q9 格子波尔兹曼方法， 在 Visual Studio 2012 平台上运用 C/C++ 语言编程，求解了顶部驱动流模型，运用 tecplot360 对所得数据进行处理，得 到了图 6 至图 8 三幅对应于不同雷诺数的流线图， 显示了该流动模型的典型现象， 验证了本文所述方法的正确性。
式中，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
在 Boltzmann 方程的基础上，LBM 引入 Chaffman Enskog 展开和 Taylor 级 数展开,用 Bhatnagar Gross Krook(BGK)格式替代碰撞过程，得到常用的 BGK 方 程。本文使用二维九速(D2Q9)不可压 LB 模型，见图 2。其演化方程为:
Fi ( x ci t , t t ) Fi ( x, t ) t
2.物理模型
u=1, v=0 u=0 v=0 y u=0 v=0 x u=0, v=0 l
图 1 物理模型示意图
l
如图 1 所示，为本文所述物理模型。考虑其为一长和宽分别为 1、1 的水平 放置的正方形。其上边界以一个固定的水平速度（U=0.1）移动，而其它三个边 界保持静止。这种流动的基本特征是，在腔的中央出现一个较大的漩涡，而在角
密度和动量可由分布函数求得：
（4）
(x,t)= Fi ( x, t )
i
（5） （6）
j(x,t)= (x,t)u(x,t)= Ci Fi ( x, t )
i
应用 Chapman-Enskog 展开方法， 可以得到不可压缩条件下的纳威-斯托克斯 方程：
u=0 i u (u )u 1 p v 2 u K p K bT / m KT 1 v b ( ) t m t 2
-9-
-5-
计算流体力学大作业
第6页 共9页
理方法得到的精度有限。如果让固壁节点参与碰撞,因为 F2 不一定总是等于 F4, 所以并不能保证壁面切向速度为零。为了满足无滑移条件,对部分分布函数 Fi 进 行重新分配,对具有切向分量的流入流体内部的分布函数 Fi 作如下分配:
Fi Fi 0 Fi ' Fi '0
（8）
其中 i 是流入的方向,I’表示其相反方向。在给定压力和速度的边界上,平衡态 函数由式(3)直接得到。在下一时刻,流出边界方向上的分布函数是已知的,再结 合式(7),使含未知分布函数 Fi 的方程组封闭,从而得到流入方向上的分布函数 值。碰撞过程在边界格点和流体格点上同时进行。
6.数值计算结果及分析
-2-
计算流体力学大作业
第3页 共9页
落处出现两个次级漩涡。 在本文中，假定所述研究对象满足以下假设： 1）定常流动。腔体中空气密度均匀，故所述流动为定常流动。 2）二维流动。由于物体的对称性，故在研究时只考虑 XY 方向。 3）刚性固体壁边界条件。 在本文中假定所述对象满足刚性固体壁边界条件。
3.控制方程 3.1 D2Q9 格子波尔兹曼方法
Re
ul Re ，其中顶盖流速 u = 0.1，l = 1
GAM 1 Re
4.程序组织
本文借助 D2Q9 格子波尔兹曼方法，在 Visual Studio 2012 平台上运用 C/C++语 言编程，运用有限差分方法，对上述传统问题进行了求解，并借助 Tecplot360 工具， 对所得数据进行后处理。所述程序组织如下图所示（完整程序参见所提交作业的文件 夹） ：
设定本例的初始条件为：顶盖驱动速度为 0.1，雷诺数分别为 2000、5000 以及 10000，计算本例，运行界面如图 4 所示。
-6-
计算流体力学大作业
第7页 共9页
图 4 程序运行界面图
对图 5 所示的数据运用 tecplot360 做后处理。
图 5 数据存储文件图
对各雷诺数情形下的结果分别进行后处理，其流线图如下面图形所示：
南 京 理 工 大 学
能 源 与 动 力 工 程 学 院
计算流体力学大作业
题 专
姓名 电话
目 业
基于 C/C++的顶部驱动圆柱 绕流模拟
动力工程
学号 512081733 成绩
教师评语：
年பைடு நூலகம்
月
日
第1页 共9页
计算流体力学大作业
第2页 共9页
基于 C/C++的顶部驱动流模拟
摘要：借助 D2Q9 格子波尔兹曼方法，在 Visual Studio 2012 平台上运用 C/C++
8.参考文献
[1] Ghia U, Ghia K N and Shin C T. High-Re solutions for incompressible flow using the Navier-Stokes equations and a multi-grid method[J]. J Comput Phys, 1982. 48: 387-411 [2] 陶文铨编著. 数值传热学（第二版），西安：西安交通大学出版社，2001 [3] 传热与流体流动的数值计算课程教学程序，EXAMPLES
Fi (0) ( x, t )
（ 2）
为平衡态分布函数，仅与当地的密度和动量有关，可以由熵增加原
理推出:
-3-
计算流体力学大作业
第4页 共9页
Fi (0) ( x, t ) Fi (0) ( , j )
更直接地表示为：
Wi
0
{
m m m ci j [ (ci j )2 j 2 ]} K BT 2 K BT K BT

[ Fi ( x, t ) Fi (0) ( x, t )]
（1）
式中，Fi 为沿 i 方向的粒子速度分布函数， x 和 t 分别为空间步长和时间步长，
为量纲一松弛时间。
图 2 二维九速(D2Q9)不可压 LB 模型
Ci 表示粒子的运动方向，课表示为
C1,3 , C2,4 (c, 0), (0, c) C5,6 , C7,8 (c, c), (c, c) C0 (0, 0)

-4-
计算流体力学大作业
第5页 共9页
uu uv 1 p 2u 2u 2 2 x y x y x uv vv 1 p 2v 2v 2 2 x y y y x
-7-
计算流体力学大作业
第8页 共9页
图 6
Re=2000 时顶盖驱动流的流线图
图 7
Re=5000 时的顶盖驱动流流线图
-8-
计算流体力学大作业
第9页 共9页
图 8
Re=10000 时顶盖驱动流流线图
如图 6 至图 8 三幅图所示，基于 C/C++顶盖驱动流的流线图中，具有如下特 征： a) 方腔中心具有一个大漩涡； b) 方腔的左上角、左下角以及右下角具有次级漩涡； c) 随着雷诺数的增大，右下角的次级漩涡逐渐变小，左上角的次级漩涡逐 渐变大。
图 3 程序组织示意图
5.边界条件
在本文中，边界条件采用格子波尔兹曼方法中的刚性固壁边界条件。以图 4 所示的左壁面为例，对于固壁上的节点 A,由 A 邻近的 B、C、D、E、F 点可以分 别确定格点 A 的 5 个方向的分布函数 F2、F6、F3、F7、F4。根据壁面反弹的思 想可知 F1=F3=F5=F7=F8=F6。则点 A 所有方向的分布函数可以确定。然而这种处