基于投影法的三维流场数值求解方法
来源:用户上传
作者:骆佳玲 李伟忠 王文全
摘 要：本文利用投影法和浸入边界法使用的固定笛卡尔网格技术，三维泊松方程（poisson equation）采用有限七点差分格式离散，结合迭代法和直接法将三维的数值求解问题简化为二维问题，然后利用二维的直接数值解法进行快速的迭代求解。采用VC++编写数值计算代码，通过三维数值算例验证了三维压力Poisson方程求解数值方法的有效性和求解精度，并以三维槽道流场为基准数值算例，验证利用投影法的三维流场数值计算结果的可靠性和适用性。
关键词：投影法; 三维Poisson方程; 三维流场
中图分类号：O35
文献标志码：A
投影法是数值求解N-S（navier stokes）方程的一种十分有效的方法，最初是由CHORIN[1]在1968年提出的，在之后的发展过程中不断得到改进，得到一系列的改进后的投影格式，使其求解精度得到了极大的提高，例如，常怀见[2]构造了一种基于非结构网格的求解非稳态流动的高精度数值算法，时间和空间的收敛精度达到了二A以上。胡锐锋等[3]提出了具有四阶空间精度的不可压缩流动的投影格式。浸入边界法中采用固定的笛卡尔网格，且在求解中不需要更新网格，在变形大的流固耦合问题计算中具有很大的优势，所以投影法和浸入边界法的结合受到了很多研究人员的关注。王文全等[4]在传统的投影法基础上发展了一种投影浸入边界法，提高了计算的精度。TIAN等[5]利用有限元法和投影浸入边界法三维数值模拟了大变形的流固耦合问题，研究了其在生物流体力学中的具体应用。李宇航[6]在投影法与直接力浸入边界法结合的基础上修正体积力，研究了仿飞蛇滑翔平板的空气动力学机理。
在N-S方程求解中，压力泊松方程（Poisson equation）的数值求解是难点之一，常用的数值求解方法有迭代法和直接法。二维Poisson方程具有五点差分格式[7]，可利用直接法进行求解，三维Poisson方程的差分格式在有些文献[8]中进行了一般性的介绍，文献[9]中也给出了具体的三维Poisson方程的七点差分格式， 但是其精度及其实用性没有得到验证。二维Poisson方程可利用直接法进行快速稳定的求解[10]，但是三维问题得到的线性方程组相对于二维情况变得很大，这对求解方程组的数值方法和计算机的存储都具有很大的挑战。
为此，本文运用七点差分格式并结合迭代法原理，将三维的数值求解简化为二维的数值求解，结合方程组求解的迭代求解法和直接求解法，建立投影法的三维数学模型和数值求解策略，对三维问题进行快速稳定的求解。采用VC++编写投影法的数值计算代码，通过求解三维数值算例验证了三维压力Poisson方程数值方法的求解精度和准确性，并以三维槽道流场为基准数值算例，验证了三维压力Poisson方程数值计算结果的可靠性和适用性。
1 数值计算方法
11 控制方程
流体为黏性不可压缩流体时，张量形式的控制方程表示如下[11]，
ui=0（1）
uit+uj・ui=-p+1ReΔui（2）
式中，ui，uj为速度矢量;p为压强;Re为流动雷诺数。
12 时间推进和对流扩散项的离散方法
利用投影法对方程（1），（2）进行分步求解，基本步骤如下：
1）速度预算步。不考虑压强对时间离散求解得到预估速度u*i，
u*i=uni-Δt（C+V）（3）
式中，C表示对流项;V表示黏性扩散项。在求解方程（3）时，对C和V的数值离散化，参照文献[10]中二维的数值离散化，将其直接推广至三维问题。
2）速度校正步。考虑压强项和瞬态时间项离散求解得校正速度u′i，
u′i=u*i-Δtpn+1xi（4）
式中，pn+1表示下一时间步的压力值，可由压力泊松方程求出，
2pn+1x2i=xiu*iΔt=fi，j，k（5）
求解（5）式时使用Neumann边界条件，
pn+1xi=-u′i-u*iΔt=Ci（6）
13 三维泊松方程的数值求解方法
在整个数值求解计算过程中，求解三维压力Poisson方程是关键也是难点。我们考虑一长方体区域，在x、y、z方向上分别划分为M、N、P等份，三个方向上的网格间距分别为Δx、Δy、Δz。利用标准的七点差分格式空间离散三维Poisson方程，即：
pi+1，j，k-2pi，j，k+pi-1，j，kΔx2+pi，j+1，k-2pi，j，k+pi，j-1，kΔy2+
pi，j，k+1-2pi，j，k+pi，j，k-1Δz2=fi，j，k（7）
将包含z方向变化的数值项视为已知数，移至等式右边整理得，
pi-1，j，kΔx2+pi，j-1，kΔy2-21Δx2+1Δy2+1Δz2pi，j，k+pi，j+1，kΔy2+pi+1，j，kΔx2=fi，j，k-pi，j，k+1+pi，j，k-1（Δz）2=Ci，j，k（8）
第一时间步内分别在各个xy平面上进行直接迭代求解，后续时间步内不用迭代法，右端项的压力项采用上一时间步的值用直接法直接求解，边界内部采用中心差分格式离散，详细的边界处理步骤见文献[10，12-13]。
离散化后，通过整理可以得到与二维同阶的大型稀疏方程组，表示如下：
AP=D （9）
式中，
P=[P0，j，k，P1，j，k，…，Pi，j，k，…，PM-1，j，k，PM，j，k]T
Pi，j，k=[pi，0，k，pi，1，k，…，pi，j，k，…，pi，N-1，k，pi，N，k]T，i=0，1，…，M
转载注明来源:https://www.xzbu.com/1/view-15423418.htm
