一种基于非结构化网格的无压力型投影法

李国杰，李 杨，王 勇，李建峰，张 晓，王憨鹰，黄 萌，周丹丹，胡广涛

(1.榆林学院 能源工程学院，陕西 榆林 719000；2.西安热工研究院有限公司，陕西 西安 710054；3.陕西四季春清洁热源股份有限公司，陕西 西安 710061；4.西安现代控制技术研究所，陕西 西安 710065)

自然界和工程应用中，流动问题常常出现在复杂的流动区域，如弯曲河流的泥沙沉积、人体血管的血液流动等。介于非结构化网格对复杂区域具有良好的适用性[1]，基于非结构化网格的数值计算方法成为复杂流动区域内流体力学分析的一种有效方法，受到越来越多学者的关注。路川藤等[2]发展了基于非结构化网格的LU-SGS隐格式算法，并将其推广到水动力学研究中，完成了大范围海域潮流泥沙运动分析；周文[3]将IDEAL算法推广应用到非结构三角网格系统中，实现了粘性-粘弹性两相流数值研究；李杰[4]基于非结构化网格对流项TVD离散格式，发展了r因子算法、SIMPLERR算法，提高了在高Re数及细密网格数值研究中的计算效率；本文作者前期发展了基于非结构化网格的非牛顿SIMPLE算法[5]、PISO算法[6]，并将其应用到了血液动力学研究。目前，基于非结构化网格的流体力学数值方法主要集中在SIMPLE类方法的研究方面，而SIMPLE类方法属于压力耦合求解方法，在计算非稳态问题时，需在每一时层内完成内迭代求解，较为耗时。

投影法(Projection Method)早期由Chorin提出[7]，属于一种速度、压力解耦算法。其基本思想是引入中间速度场，将速度场和压力场的求解从动量方程中解耦，大大降低速度场和压力场耦合求解所带来的计算量，在计算非定常不可压缩流动的流场求解中优势明显。介于此，投影法得到了国内外学者的不断研究发展。Choi和Moin[8]提出了一种基于结构化网格的四步投影法，对流项和扩散项军采用隐式格式，整体计算精度达到二阶；Ni和Mohamed在Choi研究的基础了，将四步投影法扩展到了多相流数值研究[9]；刘淼儿[10]通过给定合适的人工边界条件将投影法提高到了三阶精度；Jan和Sheu[11]提出了一种基于非结构的准隐式投影法，采用了动量差值方法消除求解中压力场的非物理震荡问题，采用了准隐式的时间离散格式，加速了收敛；王坪[12]提出了一种基于结构化网格的不可压缩流体投影法，在小尺度流动问题数值模拟中具有较好的可靠性；常怀见[13]构造了一种基于非结构网格的高精度投影算法，时间和空间的收敛精度达到了二阶以上；王文全等[14]根据投影浸入边界法分步投影求解的特点，提出了一种投影浸入边界法快速求解方法，保证了计算的稳定性。根据求解中间速度场时是否引入压力项，现有投影法类方法又分为了压力增量型投影法[8-11]和无压力型投影法[12-13]。

目前很多学者已将压力型投影法扩展到了非结构化网格中[11]，而基于非结构化网格的无压力投影法发展相对还不够成熟[13]。因此，本文拟发展一种基于非结构化的无压力型投影法，使其能够完成复杂流动区域内非定常流动问题的高效求解。

1 控制方程的离散

对于不可压缩流体的非稳态流动，其通用控制方程可以写为：

(1)

式中，u为速度矢量；φ为广义变量，可为速度矢量u三个不同方向分量u，v，w等待求变量；λφ为广义扩散系数，一般为流体动力粘度；Qφ为广义源项，外力、磁场力等包含于广义源项中，不包含压力梯度项。

在任意一个非结构化网格单元P0中，对通用控制方程进行时间和空间积分，整理得：

(2)

式中，n为非结构化网格边界数(亦是非结构化网格相邻网格数)。离散中用到的几何参数如1所示：Sj是界面面积矢量，方向垂直界面向外；dj是控制容积中心节点P0到相邻网格中心节点Pj的距离矢量；rj-rpj是中心节点Pj到界面交点j的距离矢量。本文算法中，速度和压力变量储存在网格中心节点上。

图1 非结构化网格中的变量设置

离散中，时变项采用一阶迎风离散格式进行离散，对流项由通过界面的流量和界面处的广义变量值相乘获得，扩散项离散为界面处的广义变量梯度与界面面积矢量的内积，广义源项的处理采用局部线性化方法[1, 15]。其中界面上的广义变量采用混合格式，变量值界面处的广义变量梯度采用仿动量插值法进行离散[15]。为消除求解过程中不合理的压力场，界面流速uj采用动量插值的方法获得，即

(3)

而界面广义变量平均梯度采用线性插值求解，即

(4)

经过上述离散后，可得到动量离散方程为：

(5)

其中,

2 不可压缩流体无压力型投影法的实施

无压力型投影法，求解中间速度场时不引入压力项，速度和压力解耦求解，计算量小[16]。算法实施中，首先基于预测方程完成中间速度求解，然后根据修正方程完成压力场求解，最后基于压力场对速度场进行修正。边界条件的处理可以参考前期研究[5]。本文算法具体实施步骤如下：

(1)中间速度预测。中间速度场时不引入压力项，具体预测步方程如下：

(6)

从而可以得到中间速度场变量为：

(7)

(2)压力场求解。基于修正方程，完成压力场求解，具体中间速度修正方程如下：

(8)

由此可以得到下一时层节点速度场变量为：

(9)

(10)

将其带入到连续性方程，可得节点压力求解方程：

(11)

其中，

(3)速度场更新。基于中间速度修正方程完成速度场更新。

3 算法验证

本文采用方腔顶盖驱动流动典型算例对算法实施的准确性进行验证。三维空间方腔驱动流示意图如图2所示。初始时刻腔内流体静止，从t=0时刻顶盖开始以速度u沿x正方向运动，腔内流体在顶盖的驱动作用下运动，最终达到稳定状态。流动中雷诺数Re的定义基于顶盖速度u和腔体深度l。采用Re=400的方腔顶盖驱动流验证流场稳定时刻结果的准确性，其方腔尺寸为1∶1∶1(x∶y∶z)，网格数为55 205，如图3(a)所示；然后采用Re=1 000的方腔顶盖驱动流验证不同时刻流场结果的准确性，其方腔尺寸为1∶2∶1(x∶y∶z)，网格数为112 803，如图3(b)所示。

图2 顶盖驱动示意图

(a) Re=400 (b) Re=1000图3 顶盖驱动方腔网格划分

图4给出了Re=400时流场稳定后中心界面(y=0)上中心线速度分布，从图中可以看出，本文计算结果与前人计算结果[17-19]吻合良好，验证本文算法实施的正确性。

图4 中心截面(y=0)上速度分布(Re=400)

图5给出了Re=1000时不同时刻下中心界面(y=0)处中心线速度分布，从图中可以看出，在时间尺度上本文算法的计算结果均与Guermond等人[20]的实验结果较好吻合，验证了本文算法在时间维度计算结果的准确性。

(a)t=4 (b)t=8

(c)t=10 (d)t=12图5 不同时刻中心线上速度分布对比(Re=1000)

4 结论

考虑到流动区域的复杂性、投影法求解非定常问题的高效性，本文发展了一种基于非结构化网格的非定常无压力型投影法。无压力型投影法在中间速度求解中不引入压力项，速度和压力完全解耦求解，计算量小，效率高。算法实施中，首先基于不包含压力项的预测方程完成中间速度求解，然后根据包含压力项的修正方程完成压力场求解，最后基于压力场对速度场进行修正。通过对方腔顶盖驱动流典型算例进行模拟分析，与前人结果进行对比，验证了本文算法实施的准确性。