程序

附录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二维不可压缩黏性

方腔流问题示意图

） u v x y x x

其中u为变量 在水平x方向的流速，v为 在垂直y方向的流速， 为黏度，S 为源项。源项中不仅包含压力梯度项，也包含时间导数项。

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

程序

边界条件：

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）改写成如下有限体积算法离散格式：

程序

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 s

Fe u eAe,Fw u wAwFn v nAn,Fs v sAs

（F.4）

同时记：

（F.5）

eAe A

,Dw ww, xPE xPW

（F.6）

nAn sAs

Dn ,Ds

yPN yPS

De

则可由式（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.4三阶迎风型QUICK离散格式示意图

程序

为了便于讨论，现对一维对流扩散方程的三阶迎风型QUICK离散格式进行分析：在三阶迎风型QUICK离散格式中，计算主控制单元界面上流动量 需要取主控制单元界面两侧3个节点处的流动量值进行插值计算得到，其中两个节点位于界面紧邻的两侧，第三个节点位于迎风一侧较远邻点，如图F.4所示。

当ue 0,uw 0时，通过WW、W和P三个节点值拟合曲线来计算主控制单元左侧界面参数 w。通过节点W、P和E三个节点值拟合曲线来计算主控制单元右侧界面参数 e。当ue 0,uw 0，则分别通过节点W、P、E和P、E、EE三个节点值计算主控制单元左、右两侧界面参数 w和 e。根据上述计算原则，可以得到界面参数 w计算公式如下：

当uw 0时，界面参数 w计算公式为：

ww wW wP wWW （F.8a）

当ue 0时，界面参数 w计算公式为：

636

e P E W （F.8b）

888

683868

对于一维无源项一维对流扩散方程三阶迎风型QUICK离散格式： 当ue 0,uw 0时，三阶迎风型QUICK离散格式为：

aP P aW W aE E aWW WW （F.9） 其中

31

aW Dw Fw Fe

883

aE De Fe

8 （F.9a） 1

aWW Fw

8

aP aW aE aWW Fe Fw

同理，若ue 0,uw 0，三阶迎风型QUICK离散格式为：

aP P aW W aE E aEE EE （F.10）

程序

其中

3

aW Dw Fw,

861

aE De Fe Fw,

88 （F.10a） 1aWW Fe

8

aP aW aE aEE Fe Fw

将两种流动方向离散方程（F.9）和（F.10）合并后，可得到统一的一维对流扩散方程三阶迎风型QUICK离散格式： …… 此处隐藏：11391字，全部文档内容请下载后查看。喜欢就下载吧 ……