附录F 有限体积法计算方腔流(F) 下载本文

附录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离散格式示意图