附录F有限体积法计算方腔流(F).docVIP

实用标准 PAGE 精彩文档 附录F 二维不可压缩黏性流体方腔流动问题的 有限体积算法与计算程序 二维方腔流动问题是一个不可压缩黏性流动中典型流动。虽然目前尚不能求得它的解析解，但是它常被用来作为检验各种数值算法计算精度和可靠性的算例。文献中几乎大多数算法都对它进行过计算。在本算例中采用有限体积算法三阶迎风型离散格式对它进行数值求解。同时，为了初学者入门和练习方便,这里给出了用语言和语言编写的计算二维不可压缩黏性方腔流动问题计算程序，供大家学习参考。 F-1利用有限体积算法三阶迎风型离散格式求解 二维不可压缩黏性流体方腔流动问题 1.二维不可压缩黏性流体方腔流动问题 二维不可压缩黏性流体方腔流动（cavity flow图F.1二维不可压缩黏性方腔流问题示意图）：有一正方形腔室，其量纲为一的宽度为，里面充满静止的不可压缩黏性流体，方腔内初始时刻压力和密度为它周围壁面（左右壁面和底面）固定不动，上壁面以量纲为一的速度沿着上壁面方向自左向右运动（图F.1）。 图F.1二维不可压缩黏性 方腔流问题示意图 2. 基本方程组、初始条件和边界条件 图F.1 二维不可压缩黏性方腔流动问题示意图设流体是黏性流体。二维方腔流动问题在数学上可以由二维不可压缩黏性流动N - S方程组来表示，把它写成通用变量的微分方程组形式，有： 图F.1 二维不可压缩黏性方腔流动问题示意图 （F.1） 其中为变量在水平方向的流速，为在垂直方向的流速，为黏度，为源项。源项中不仅包含压力梯度项，也包含时间导数项。 初始条件：方腔上壁面以量纲为一的速度沿着上壁面方向自左向右运动。 边界条件： 流动速度均可采用无滑移边界条件，压力采用自由输出边界条件。 3.计算网格划分和控制体单元与节点定义 采用交错网格，图F.2和图F.3是计算网格、控制体单元和节点示意图。 图F.2方腔流动计算网格、 图F.2方腔流动计算网格、 控制体单元和节点示意图 图F.3计算采用的交错网格示意图 节点所在主控制单元如图F.2中有阴影部分所示。在方向与节点相邻的节点为和，在方向与节点相邻的节点为和，主控制单元界面分别为。压力和速度分别在三套不同网格中如图F.3中有阴影部分所示。 4.有限体积算法三阶迎风型离散格式 对方程（F.1）在图F.2所示节点所在控制体单元内积分，有： （F.2） 由于二维不可压缩黏性流体方腔流动是二维问题，因此控制体单元体积仅是面积，而它的边界是长度。设 ，利用定理，可将方程（F.2）改写成如下有限体积算法离散格式： （F.3） 对上式中采用一阶向前差分近似，则有： （F.4） 同时记： （F.5） （F.6） 则可由式（F.2）写成： （F.7） 式中都是控制体单元内节点上的已知量，如果利用差分计算得到控制体单元边界上的流通量，就可以求出节点上未知量。 图F.4三阶迎风型离散格式 图F.4三阶迎风型离散格式示意图 为了便于讨论，现对一维对流扩散方程的三阶迎风型离散格式进行分析：在三阶迎风型离散格式中，计算主控制单元界面上流动量需要取主控制单元界面两侧3个节点处的流动量值进行插值计算得到，其中两个节点位于界面紧邻的两侧，第三个节点位于迎风一侧较远邻点，如图F.4所示。 当时，通过、和三个节点值拟合曲线来计算主控制单元左侧界面参数。通过节点、和三个节点值拟合曲线来计算主控制单元右侧界面参数。当，则分别通过节点、、和、、三个节点值计算主控制单元左、右两侧界面参数和。根据上述计算原则，可以得到界面参数计算公式如下： 当时，界面参数计算公式为： （F.8a） 当时，界面参数计算公式为： （F.8b） 对于一维无源项一维对流扩散方程三阶迎风型离散格式： 当时，三阶迎风型离散格式为： （F.9） 其中 （F.9） 同理，若，三阶迎风型离散格式为： （F.10） 其中 （F.10a） 将两种流动方向离散方程（F.9）和（F.10）合并后，可得到统一的一维对流扩散方程三阶迎风型离散格式： （F.11） 其中