附录F-有限体积法计算方腔流(F).doc

《附录F-有限体积法计算方腔流(F).doc》由会员上传分享，免费在线阅读，更多相关内容在教育资源-天天文库。

1、附录F二维不可压缩黏性流体方腔流动问题的有限体积算法与计算程序二维方腔流动问题是一个不可压缩黏性流动中典型流动。虽然目前尚不能求得它的解析解，但是它常被用来作为检验各种数值算法计算精度和可靠性的算例。文献中几乎大多数算法都对它进行过计算。在本算例中采用有限体积算法三阶迎风型离散格式对它进行数值求解。同时，为了初学者入门和练习方便,这里给出了用语言和语言编写的计算二维不可压缩黏性方腔流动问题计算程序，供大家学习参考。F-1利用有限体积算法三阶迎风型离散格式求解二维不可压缩黏性流体方腔流动问题1.二维不可压缩黏性流体方腔流动问题二维不可压缩黏性流体方腔流动（cavit

2、yflow图F.1二维不可压缩黏性方腔流问题示意图）：有一正方形腔室，其量纲为一的宽度为，里面充满静止的不可压缩黏性流体，方腔内初始时刻压力和密度为它周围壁面（左右壁面和底面）固定不动，上壁面以量纲为一的速度沿着上壁面方向自左向右运动（图F.1）。2.基本方程组、初始条件和边界条件图F.1二维不可压缩黏性方腔流动问题示意图设流体是黏性流体。二维方腔流动问题在数学上可以由二维不可压缩黏性流动N-S方程组来表示，把它写成通用变量的微分方程组形式，有：（F.1）其中为变量在水平方向的流速，为在垂直方向的流速，为黏度，为源项。源项中不仅包含压力梯度项，也包含时间导数项。初

3、始条件：方腔上壁面以量纲为一的速度沿着上壁面方向自左向右运动。边界条件：流动速度均可采用无滑移边界条件，压力采用自由输出边界条件。3.计算网格划分和控制体单元与节点定义采用交错网格，图F.2和图F.3是计算网格、控制体单元和节点示意图。图F.2方腔流动计算网格、控制体单元和节点示意图图F.3计算采用的交错网格示意图节点所在主控制单元如图F.2中有阴影部分所示。在方向与节点相邻的节点为和，在方向与节点相邻的节点为和，主控制单元界面分别为。压力和速度分别在三套不同网格中如图F.3中有阴影部分所示。4.有限体积算法三阶迎风型离散格式对方程（F.1）在图F.2所示节点所在

4、控制体单元内积分，有：（F.2）由于二维不可压缩黏性流体方腔流动是二维问题，因此控制体单元体积仅是面积，而它的边界是长度。设，利用定理，可将方程（F.2）改写成如下有限体积算法离散格式：（F.3）对上式中采用一阶向前差分近似，则有：（F.4）同时记：（F.5）（F.6）则可由式（F.2）写成：（F.7）式中都是控制体单元内节点上的已知量，如果利用差分计算得到控制体单元边界上的流通量，就可以求出节点上未知量。图F.4三阶迎风型离散格式示意图为了便于讨论，现对一维对流扩散方程的三阶迎风型离散格式进行分析：在三阶迎风型离散格式中，计算主控制单元界面上流动量需要取主控制单

5、元界面两侧3个节点处的流动量值进行插值计算得到，其中两个节点位于界面紧邻的两侧，第三个节点位于迎风一侧较远邻点，如图F.4所示。当时，通过、和三个节点值拟合曲线来计算主控制单元左侧界面参数。通过节点、和三个节点值拟合曲线来计算主控制单元右侧界面参数。当，则分别通过节点、、和、、三个节点值计算主控制单元左、右两侧界面参数和。根据上述计算原则，可以得到界面参数计算公式如下：当时，界面参数计算公式为：（F.8a）当时，界面参数计算公式为：（F.8b）对于一维无源项一维对流扩散方程三阶迎风型离散格式：当时，三阶迎风型离散格式为：（F.9）其中（F.9）同理，若，三阶迎风型

6、离散格式为：（F.10）其中（F.10a）将两种流动方向离散方程（F.9）和（F.10）合并后，可得到统一的一维对流扩散方程三阶迎风型离散格式：（F.11）其中（F.11a）式中（F.11b）同理，可以得到带有源项的二维对流扩散方程三阶迎风型离散格式为：（F.12）其中为有限体积算法中源项平均值。式中各个系数为：（F.12a）式中（F.12b）源项为：（F.13）若把表示时刻动量，表示时刻动量，则可以得到源项离散格式为：（F.14）最后，得到有限体积算法二维对流扩散方程三阶迎风型离散格式：（F.15）式中系数为一阶迎风格式中各对应系数。5.计算结果分析利用三阶迎风

7、型离散格式和相应的初始条件和边界条件，求解二维不可压缩黏性流体方腔流动问题。图F.5是不同雷诺数条件下采用三阶迎风型离散格式得到的二维不可压缩黏性流体方腔流动的计算结果。计算结果和文献中其他高精度算法得到的计算结果进行了比较，两者计算结果十分吻合，能把方腔下壁面两个底角附近二次小涡清晰地计算出来。这表明有限体积算法三阶迎风型离散格式具有相当高的计算精度。=1000=100=10000=5000图F.5不同雷诺数条件下采用三阶迎风型离散格式计算二维不可压缩黏性方腔流动的计算结果从图F.5中可以看出：二维不可压缩黏性流体方腔流动的中心大涡并不在中心位置，方腔内流动也并

8、不对称。这