附录F 二维不可压缩黏性流体方腔流动问题的
有限体积算法与计算程序
二维方腔流动问题是一个不可压缩黏性流动中典型流动。虽然目前尚不能求得它的解析解，但是它常被用来作为检验各种数值算法计算精度和可靠性的算例。文献中几乎大多数算法都对它进行过计算。在本算例中采用有限体积算法三阶迎风型QUICK离散格式对它进行数值求解。同时，为了初学者入门和练习方便,这里给出了用C语言和Fortran77语言编写的计算二维不可压缩黏性方腔流动问题计算程序，供大家学习参考。
F-1利用有限体积算法三阶迎风型QUICK离散格式求解
二维不可压缩黏性流体方腔流动问题
1.二维不可压缩黏性流体方腔流动问题
二维不可压缩黏性流体方腔流动（cavity flow）：有一正方形腔室，其量纲为一的宽度为1.0，里面充满静止的不可压缩黏性流体，方腔内初始时刻压力和密度
?=1.0它周围壁面（左右壁面和为p=1.0、底面）固定不动，上壁面以量纲为一的速度1.0沿着上壁面方向自左向右运动（图F.1）。
2. 基本方程组、初始条件和边界条件
设流体是黏性流体。二维方腔流动问题在数学上可以由二维不可压缩黏性流动
N - S方程组来表示，把它写成通用变量的微分方程组形式，有：
图F.1 二维不可压缩黏性方 ??腔流动问题示意图???图F.1二维不可压缩黏性
方腔流问题示意图
??????（F.1） ???u?????v???????????S? ?x?y?x??x??y??y?其中u为变量?在水平x方向的流速，v为?在垂直y方向的流速，?为黏度，S?为源项。源项中不仅包含压力梯度项，也包含时间导数项。
初始条件：方腔上壁面以量纲为一的速度1.0沿着上壁面方向自左向右运动。
-F.1-
边界条件：
v均可采用无滑移边界条件，压力p采用自由输出边界条件。 流动速度u、
3.计算网格划分和控制体单元与节点定义
采用交错网格，图F.2和图F.3是计算网格、控制体单元和节点示意图。
图F.2方腔流动计算网格、 控制体单元和节点示意图
图F.3计算采用的交错网格示意图
节点P所在主控制单元如图F.2中有阴影部分所示。在x方向与节点P相邻的节点为W和E，在y方向与节点P相邻的节点为S和N，主控制单元界面分
sen。压力p和速度u、v分别在三套不同网格中如图F.3中有阴影部分别为w、、、所示。
4.有限体积算法三阶迎风型QUICK离散格式
对方程（F.1）在图F.2所示节点P所在控制体单元内积分，有：
???????????? ?u?dV??v?dV??dV????????dV??S?dV（F.2）???????x?y?x??x??y??y??V?V?V?V?V由于二维不可压缩黏性流体方腔流动是二维问题，因此控制体单元体积?V仅是面积，而它的边界是长度。设 Aw?Ae??y,As?An??x，利用Gauss定理，可将方程（F.2）改写成如下有限体积算法离散格式：
-F.2-
????u?A?e???u?A?w???????v?A?n???v?A?s?? （F.3） ???????????????????A????A????A????A??S??x?y?x?e??x?w??y?n??y?s?????????????????对上式中??、??、??采用一阶向前差分近似，则有： ??、?x?x?y??e??w??n??y?s
?P??W?????E??P?????,??????x?x??x?e??x?w?????N??P?????P??S,???????y?y?y??n??y?sFe???u?eAe,Fw???u?wAwFn???v?nAn,Fs???v?sAs （F.4）
同时记：
（F.5）
?eAe?A,Dw?ww,?xPE?xPW （F.6）
?nAn?sAsDn?,Ds??yPN?yPSDe?则可由式（F.2）写成：
Fe?e?Fw?w?Fn?n?Fs?s?De??E??P??Dw??P??W??Dn??N??P??Ds??P??S??S??x?y （F.7）
式中?P、?E、?W、?N、?S、De、Dw、Dn、Ds都是控制体单元内节点上的已知量，如果利用差分计算得到控制体单元边界上的流通量Fe?e、Fw?w、Fn?n、Fs?s，就可以求出节点上未知量?P、?E、?W、?N、?S。
-F.3-
图F.4三阶迎风型QUICK离散格式示意图