附录E有限体积法求解库塔流E.doc

-E. PAGE 16- 附录E 二维不可压缩黏性流体流动问题 的有限体积算法与计算程序 二维流动是一个不可压缩黏性流动的典型流动，并有解析解，可以用来检验数值算法计算精度和可靠性。对它采用有限体积算法一阶迎风型离散格式进行数值求解。同时，为了初学者入门和练习方便,这里给出了用语言和语言编写的计算二维不可压缩黏性流体流动问题计算程序，供大家学习参考。 E-1利用有限体积算法一阶迎风型离散格式求解二维 不可压缩黏性流体流动问题 1.二维不可压缩黏性流体流动问题的提法 二维不可压缩黏性流体流动：有两个无限长平板，间距为，两板之间充满密度为1、静止的不可压缩黏性流体。无限长平板组成的通道两端压力相等，下板固定不动，上板以量纲为一的速度1自左向右平移运动（图E.1）。 图E.1二维不可压缩黏性流体流动问题示意图 2. 基本方程组、初始条件和边界条件 设流体是黏性流体。二维流动问题在数学上可以由二维不可压缩黏性流动方程组来表示，把它写成通用变量的微分方程组形式，有： （E.1） 其中为变量在水平方向的流速，为在垂直方向的流速，为黏度，为源项。源项中不仅包含压力梯度项，也包含时间导数项。 初始条件：上板以量纲为一的速度1沿着上壁面方向自左向右运动。 边界条件：流动速度在上下壁面可采用无滑移边界条件，在左右两端采用自由输出边界条件；压强采用自由输出边界条件。 3. 计算网格划分和控制体单元与节点定义 采用交错网格。图E.2和图E.3是计算网格、控制体单元和节点示意图。 图E.3计算采用的交错网格示意图 图E.2方腔流动计算网格、 控制体单元和节点示意图 节点所在主控制体单元如图E.2中有阴影部分所示。在方向与节点相邻的节点为和，在方向与节点相邻的节点为和，主控制单元界面分别为。压力和速度分别在三套不同网格中，如图E.3中有阴影部分所示。 4. 有限体积算法离散格式 对方程（E.1）在图E.2所示节点主控制单元内积分，有： （E.2） 由于不可压缩黏性流体流动是二维问题，因此，控制体单元体积仅是面积，而它的边界是长度。设 ，利用定理，可将方程（E.2）改写成如下有限体积离散格式： （E.3） 对上式中采用一阶向前差分近似，则有： （E.4） 同时记： （E.5） （E.6） 则式（E.2）写成： （E.7） 式中都是主控制单元内节点上的已知量，如果利用差分计算得到主控制单元边界上的流通量，就可以求出节点上的未知量。、 为了便于讨论，现对一维对流扩散方程的一阶迎风型离散格式进行分析。首先讨论无源项一维对流扩散方程的一阶迎风型离散格式。 当流动为正向时，即，主控制单元界面取值： （E.8） 则方程（E.7）离散为： （E.9） 当流动为负向时，即，控制体单元界面取值为： （E.10） 则方程（E.7）离散为： （E.11） 将两种流动方向离散格式（E.9）和（E.11）合并后，得到统一的一维对流扩散方程一阶迎风型离散格式表达式： （E.12） 式中 （E.12a） 同理，可以得到带有源项的二维对流扩散方程的一阶迎风型离散格式： （E.13） 式中 （E.13a） 源项为： （E.14） 若把表示为时刻的动量，表示时刻的动量，则可以得到源项离散格式为： （E.15） 最后，得到带有源项的二维对流扩散方程有限体积算法一阶迎风型显式离散格式： （E.16） 式中系数为一阶迎风格式中各对应系数。 5. 计算结果分析 利用上述有限体积算法一阶迎风型离散格式和相应的初始条件和边界条件，采用MAC算法中压力耦合方程，求解二维不可压缩黏性流体流动问题。图E.4 给出了二维不可压缩黏性流体流动沿y方向的速度分布，并和精确解进行了比较，十分吻合。图E.5是二维不可压缩黏性流体流动水平x方向的速度云图。图E.6是二维不可压缩黏性流体流动速度矢量分布图。 图E.4 二维不可压缩黏性流体流动沿y方向的速度分布 图E.5 二维不可压缩黏性流体流动水平x方向速