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

利用有限体积算法三阶迎风型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 yx x y y φρφρφρφρφφφφφμμμμ⎡⎤⎡⎤-+-⎣⎦⎣⎦⎛⎫⎛⎫∂∂∂∂⎛⎫⎛⎫=-+-+∆∆ ⎪ ⎪ ⎪ ⎪∂∂∂∂⎝⎭⎝⎭⎝⎭⎝⎭ （F.3） 对上式中e w n sx x y y φφφφ⎛⎫⎛⎫∂∂∂∂⎛⎫⎛⎫ ⎪ ⎪ ⎪ ⎪∂∂∂∂⎝⎭⎝⎭⎝⎭⎝⎭、、、采用一阶向前差分近似，则有： 图F.2方腔流动计算网格、 控制体单元和节点示意图图F.3计算采用的交错网格示意图,,P W E P e wN P P Sn sx x x x y y y y φφφφφφφφφφφφ--∂∂⎛⎫⎛⎫== ⎪ ⎪∂∆∂∆⎝⎭⎝⎭⎛⎫⎛⎫--∂∂== ⎪ ⎪∂∆∂∆⎝⎭⎝⎭ （F.4）同时记：()()()(),,e e w we w n n s sn s F u A F u A F v A F v A ρρρρ==== （F.5）,,,e e w w e w PE PWn n s sn s PN PSA AD 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 φ计算公式为：636888ww wW wP wWW φφφφ=+- （F.8a ）当0e u >时，界面参数w φ计算公式为：636888e 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） 其中()31883818W w w eE e eWW wP 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） 其中()3,861,8818W w w E e e w WW eP W E EE e w a D F a D F F a F a a a a F F =+=--==+++- （F.10a ） 将两种流动方向离散方程（F.9）和（F.10）合并后，可得到统一的一维对流扩散方程三阶迎风型QUICK 离散格式：P P W W E E EE EE EE EE a a a a a φφφφφ=+++（F.11） 其中()()()()()61318881836111888118W w w w e e w wWW w wE e e e e e w w EE e eP W E WW EE e w a D F F F a F a D F F F a F a a a a a F F αααααααα=+++-=-=-----=-=++++- （F.11a ）式中01000100w w w w e e e e F F F F αααα>=<=>=<=当时，；当时，；当时，；当时，。

（F.11b ）同理，可以得到带有源项的二维对流扩散方程三阶迎风型QUICK 离散格式为：P P W W E E EE EE EE EE S S N N SS SS NN NN a a a a a a a a a S Vφφφφφφφφφφ=++++++++∆ （F.12）其中S φ为有限体积算法中源项平均值。

式中各个系数为：()()()61318881836111888W w w w e e w wWW w wE e e e e e w wa D F F F a F a D F F F ααααααα=+++-=-=----- ()()()()()11861318881836111888118EE e e S s s s n n s sSS s s N n n n n n s sNN n na F a D F F F a F a D F F F a F ααααααααα=-=+++-=-=-----=-()P W E WW EE S N SS NN e w n s a a a a a a a a a F F F F =++++++++-+- （F.12a ）式中0100,,,k k k k F F k w e s nαα>=<==当时，；当时，。

（F.12b ）源项S φ为：()u pS t xφρ∂∂=--∂∂ （F.13） 若把()nu ρ表示n t 时刻动量，()1n u ρ+表示1n t +时刻动量，则可以得到源项S φ离散格式为：()()()1n nPPe w Vu u S dV x y p p y tφρρ+∆-⋅=-∆∆--∆∆⎰ （F.14）最后，得到有限体积算法二维对流扩散方程三阶迎风型QUICK 离散格式：()()()1n n n n nP P W W E E S S N N n nnn PPewa u a u a u a u a u u u x y p pytρρ+=+++--∆∆--∆∆ （F.15）式中系数k a 为一阶迎风格式中各对应系数。

5.计算结果分析利用三阶迎风型QUICK 离散格式和相应的初始条件和边界条件，求解二维不可压缩黏性流体方腔流动问题。

图F.5是不同雷诺数Re 条件下采用三阶迎风型QUICK 离散格式得到的二维不可压缩黏性流体方腔流动的计算结果。

计算结果和文献中其他高精度算法得到的计算结果进行了比较，两者计算结果十分吻合，能把方腔下壁面两个底角附近二次小涡清晰地计算出来。

这表明有限体积算法三阶迎风型QUICK 离散格式具有相当高的计算精度。

图F.5不同雷诺数Re 条件下采用三阶迎风型QUICK离散格式计算二维不可压缩黏性方腔流动的计算结果Re =100 Re =1 000Re =5 000 Re =10 000从图F.5中可以看出：二维不可压缩黏性流体方腔流动的中心大涡并不在中心位置，方腔内流动也并不对称。

这是因为，方腔上壁面以量纲为一的速度1.0沿着上壁面方向自左向右运动时，在方腔上壁面两侧的两个顶角处不再满足边界条件，这是一个带奇性的方腔流动。

计算结果表明，方腔流动和雷诺数有关，随着雷诺数Re增加，计算精度在降低。

当雷诺数Re较低时，方腔下壁面的两个底角附近的二次小涡十分清晰，随着雷诺数Re的增加，二次小涡变得越来越模糊。

由于三阶迎风型QUICK离散格式计算精度较高，因此三阶迎风型QUICK离散格式计算效果一般要比一阶迎风型离散格式相对来说好些。

using System;using System.Collections.Generic;using ponentModel;using System.Data;using System.Drawing;using System.Linq;using System.Text;using System.IO;using System.Windows.Forms;using ;using System.Collections;namespace FVM_计算方腔流场{public partial class Form1 : Form{public Form1(){InitializeComponent();}public int MX= 30;public int MY =30;public double Re =100.00;public double Ha = 10.00;public double Ri = 1.0;public double Pr = 0.733;public double dt =0.0005;public double c2 =2.25;public double[,] u, v, p, t,uo, vo, po,to, un, vn, pn,tn,streamline;//public double dx, dy;FileStream fileResult, file1,f4;StreamWriter fileResultWriter, fileW1,fw4;public double dx, dy, err, value,x,y;private void Form1_Load(object sender, EventArgs e){u = new double[MX + 1, MY + 2];v = new double[MX + 2, MY + 1];p = new double[MX + 2, MY + 2];t = new double[MX + 2, MY + 1];uo = new double[MX + 1, MY + 2];vo = new double[MX + 2, MY + 1];po = new double[MX + 2, MY + 2];to = new double[MX + 2, MY + 1];un = new double[MX + 1, MY + 2];vn = new double[MX + 2, MY + 1];pn = new double[MX + 2, MY + 2];tn = new double[MX + 2, MY + 1];streamline = new double[MX + 1, MY + 1];dx = 1.0 / MX;dy = 1.0 / MY;fileResult = new FileStream(" X_Matrix.txt", FileMode.OpenOrCreate); fileResultWriter = new StreamWriter(fileResult);file1 = new FileStream(" file1.txt", FileMode.OpenOrCreate);fileW1 = new StreamWriter(file1);f4 = new FileStream(" f4.txt", FileMode.OpenOrCreate);fw4 = new StreamWriter(f4);//double dx1, dy1,//init(u, v, p, dx, dy);//MessageBox.Show(u[5, MY + 1].ToString());}public double max(double a, double b){if (a < b)return b;elsereturn a;}public double alfa(double x){if (x >= 0)return 1.0;elsereturn 0.0;}void init(double[,] u, double[,] v, double[,] p, double[,] t, double dx, double dy) {//u = new double[MX + 1, MY + 2];//v = new double[MX + 2, MY + 1];//p = new double[MX + 2, MY + 2];int i,j;//dx=1.0/MX;//dy=1.0/MY;for(i=0;i<=MX;i++){for(j=0;j<=MY+1;j++){u[i,j]=0.0;if(j==MY+1) u[i,j]=4.0/3.0;if(j==MY) u[i,j]=2.0/3.0;}}for(i=0;i<=MX+1;i++)for(j=0;j<=MY;j++)v[i,j]=0.0;for(i=0;i<=MX+1;i++)for(j=0;j<=MY+1;j++)p[i,j]=1.0;for (i = 0; i <= MX+1; i++){for (j = 0; j <= MY; j++){t[i, j] = 0.0;if (i == MX + 1) t[i, j] = 4.0 / 3.0;if (i == MX) t[i, j] = 2.0 / 3.0;}}}//// 一阶迎风型离散格式// 二维的三阶迎风型离散格式为9点格式，因此有两层边界网格需要//处理，本程序采用一阶迎风型离散格式处理内层，用物理边界条件处理外层。