有限容积法求解-方腔流动与传热

利用有限体积算法三阶迎风型QUICK 离散格式求解

二维不可压缩黏性流体方腔流动问题

1.二维不可压缩黏性流体方腔流动问题

二维不可压缩黏性流体方腔流动（cavity flow ）：有一正方形腔室，其量纲为一的宽度为1.0，里面充满静止的不可压缩黏性流体，方腔内初始时刻压力和密度为=1.0=1.0p ρ、它周围壁面（左右壁面和底面）固定不动，上壁面以量纲为一的速度1.0沿着上壁面方向自左向右运动（图F.1）。

2. 基本方程组、初始条件和边界条件

设流体是黏性流体。二维方腔流动问题在数学上可以由二维不可压缩黏性流动

N - S

()()u v x y x x φρφρφμ∂∂∂∂⎛⎫

+= ⎪∂∂∂∂⎝⎭⎝⎭

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

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

流动速度u v 、均可采用无滑移边界条件，压力p 采用自由输出边界条件。

3.计算网格划分和控制体单元与节点定义

采用交错网格，图F.2和图F.3是计算网格、控制体单元和节点示意图。

图F.1二维不可压缩黏性

方腔流问题示意图

节点P 所在主控制单元如图F.2中有阴影部分所示。在x 方向与节点P 相邻的节点为W 和E ，在y 方向与节点P 相邻的节点为S 和N ，主控制单元界面分

别为w s e n 、、、

。压力p 和速度u v 、分别在三套不同网格中如图F.3中有阴影部分所示。

4.有限体积算法三阶迎风型QUICK 离散格式

对方程（F.1）在图F.2所示节点P 所在控制体单元内积分，有：

()()V V V V V u dV v dV dV dV S dV x y x x y y φφφρφρφμμ∆∆∆∆∆⎛⎫

∂∂∂∂∂∂⎛⎫+=++ ⎪ ⎪∂∂∂∂∂∂⎝⎭⎝⎭⎰⎰⎰⎰⎰（F.2） 由于二维不可压缩黏性流体方腔流动是二维问题，因此控制体单元体积V ∆仅是面积，而它的边界是长度。设 ,w e s n A A y A A x ==∆==∆，利用Gauss 定理，可将方程（F.2）改写成如下有限体积算法离散格式：

()()()()e w n s e w n s u A u A v A v A A A A A S x y

x x y y φρφρφρφρφφφφφμμμμ⎡⎤⎡⎤-+-⎣⎦⎣⎦

⎛⎫⎛⎫∂∂∂∂⎛

⎫⎛⎫=-+-+∆∆ ⎪ ⎪ ⎪ ⎪∂∂∂∂⎝

⎭⎝⎭⎝⎭⎝⎭ （F.3） 对上式中e w n s

x x y y φφφφ⎛⎫⎛⎫

∂∂∂∂⎛⎫⎛⎫ ⎪ ⎪ ⎪ ⎪∂∂∂∂⎝⎭⎝⎭⎝⎭⎝⎭、、

、采用一阶向前差分近似，则有： 图F.2方腔流动计算网格、 控制体单元和节点示意图

图F.3计算采用的交错网格示意图

,,P W E P e w

N P P S

n s

x x x x y y y y φφφφφφφφφφφφ--∂∂⎛⎫⎛⎫

== ⎪ ⎪

∂∆∂∆⎝⎭⎝⎭⎛⎫⎛⎫--∂∂=

= ⎪ ⎪∂∆∂∆⎝⎭⎝⎭ （F.4）

同时记：

()()()(),,e e w w

e w n n s s

n s F u A F u A F v A F v A ρρρρ==== （F.5）

,,,e e w w e w PE PW

n n s s

n s PN PS

A A

D D x x A A D D y y μμδδμμδδ=

===

（F.6）

则可由式（F.2）写成：

()()()()e e w w n n s s e E P w P W n N P s P S F F F F D D D D S x y

φφφφφφφφφφφφφ-+-=---+

---+∆∆ （F.7）

式中P E W N S e w n s D D D D φφφφφ、、、、、、、、都是控制体单元内节点上的已知量，如果利用差分计算得到控制体单元边界上的流通量e e w w n n s s F F F F φφφφ、、、，就可以求出节点上未知量P E W N S φφφφφ、、、、。

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

图F.4三阶迎风型QUICK 离散格式示意图

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

当0w u >时，界面参数w φ计算公式为：

636888

ww wW wP wWW φφφφ=+- （F.8a ）

当0e u >时，界面参数w φ计算公式为：

636

888

e P E W φφφφ=+- （F.8b ）

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

P P W W E E WW WW a a a a φφφφ=++ （F.9） 其中

()

31

883

81

8

W w w e

E e e

WW w

P W E WW e w a D F F a D F a F a a a a F F =++=-=-=+++- （F.9a ） 同理，若0,0e w u u <<，三阶迎风型QUICK 离散格式为：

P P W W E E EE EE a a a a φφφφ=++ （F.10） 其中