一种基于四阶紧致差分格式的不可压缩流动模拟方法与流程

本发明属于离散信号计算领域，尤其涉及一种基于四阶紧致差分格式的不可压缩流动模拟方法。

背景技术：

在完全交错网格中基于四阶紧致差分格式的不可压缩流动模拟方法，一方面通过采用四阶紧致差分格式提高数值方法的精度，从而提高数值模拟的准确性；另一方面，在完全交错网格中的投影方法能够精确保证求解结果满足质量守恒方程，并且不会出现非交错网格中的压力振荡现象。然而，基于四阶紧致格式的不可压缩流动模拟方法中，压力泊松方程的求解是很大的难点，这是本发明要解决的关键问题。此方法属于通用计算流体力学方法(cfd)，能够应用于各种不可压缩流动问题的仿真模拟，比如管道流动、槽道流动、机翼绕流、建筑物绕流等。

综上所述，现有技术存在的技术问题是：传统不可压缩流动求解(包括各种商用cfd软件)一般采用的二阶格式方法，对于计算气动力来讲精度是足够的；但如今工业界要求对由湍流引起的气动噪声、压力脉动等现象有更高的预测精度，而二阶格式本身精度较低、并有较大数值耗散，在应用中有很大局限性；本发明采用四阶紧致差分格式，能够显著提高计算精度和降低误差，非常适用于湍流的数值模拟以及在气动噪声、压力脉动预测等领域的应用。

技术实现要素：

为解决现有技术存在的问题，本发明的目的在于提供一种基于四阶紧致差分格式的不可压缩流动模拟方法。

本发明是这样实现的，一种完全交错网格中基于四阶紧致差分格式的不可压缩流动模拟方法，该完全交错网格中基于四阶紧致差分格式的不可压缩流动模拟方法包括：

在时间方向，通过微分算子的近似分裂实现半隐式crank-nicolson格式离散粘性项，对流项离散则采用显式adams-bashforth格式；

对于精确投影步所产生的压力poisson方程，将压力poisson方程变换至fourier空间并得到四阶紧致差分算子的修正波数，通过代数运算获得fourier空间中的压力解后，再经过fourier反变换得到物理空间中精确满足速度散度为零的压力场。

进一步，adams-bashforth格式离散对流项和涡粘系数项，得到的离散化方程为

其中

进一步，通过微分算子的近似分裂实现半隐式crank-nicolson格式离散粘性项的方法,包括：

通过block-lu分解将与δp解耦后得到的方程为

或

进一步，待求解的x方向动量方程为

进一步，待求解的y方向动量方程为

进一步，待求解的z方向动量方程为

进一步，采用四阶紧致格式离散压力梯度项，则有

即

其中

是n3×n3矩阵，矩阵；通过求解上述线性方程组得到压力梯度项的四阶离散结果。

进一步，所述压力poisson方程为：

dgδp＝q

其中，右端项

本发明提供的基于四阶紧致差分格式的不可压缩流动模拟方法，与现有四阶精度的求解方法相比，本发明通过采用fft变换求解压力泊松方程从而大大提高了计算效率、而采用紧致差分格式使得边界条件的处理更为简单。对二维taylor-oseen涡算例的计算结果表明，本发明所提出的方法具有四阶空间精度，即随着网格尺寸的减小，离散误差以网格尺寸的四次方减小。对reτ＝180的槽道湍流dns模拟结果与谱方法dns结果吻合很好，最大误差不超过5％。

附图说明

图1是本发明实施例提供的完全交错网格中基于四阶紧致差分格式的不可压缩流动模拟方法流程图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

下面结合附图对本发明的应用原理作详细描述。

如图1所示，一种完全交错网格中基于四阶紧致差分格式的不可压缩流动模拟方法，包括：

s101：在时间方向，通过微分算子的近似分裂实现半隐式crank-nicolson格式离散粘性项，对流项离散则采用显式adams-bashforth格式；

s102：对于精确投影步所产生的压力poisson方程，将压力poisson方程变换至fourier空间并得到四阶紧致差分算子的修正波数，通过代数运算获得fourier空间中的压力解后，再经过fourier反变换得到物理空间中精确满足速度散度为零的压力场。

进一步，adams-bashforth格式离散对流项和涡粘系数项，得到的离散化方程为

其中

进一步，通过微分算子的近似分裂实现半隐式crank-nicolson格式离散粘性项的方法,包括：

通过block-lu分解将与δp解耦后得到的方程为

或

进一步，待求解的x方向动量方程为

进一步，待求解的y方向动量方程为

进一步，待求解的z方向动量方程为

进一步，采用四阶紧致格式离散压力梯度项，则有

即

其中

是n3×n3矩阵，矩阵；通过求解上述线性方程组得到压力梯度项的四阶离散结果。

进一步，所述压力poisson方程为：

dgδp＝q

其中，右端项

下面结合实施例对本发明的应用原理作进一步描述。

1、控制方程为三维非定常不可压缩navier-stokes方程

以及连续方程

在三维直角坐标系中，navier-stokes方程和连续方程的完整形式为：

2、本发明变量定义

空间离散采用交错网格。三个方向的速度以及压力定义为u(n1,n2-1,n3)、υ(n1,n2,n3)、w(n1,n2-1,n3)、p(n1-1,n2-1,n3-1)。x、y和z三个方向的网格间距分别为δx、δy和δz，其中x和z方向采用均匀网格，而y方向采用非均匀网格。

定义hj为垂向相邻网格单元中心的间距，即

定义与流向和展向速度垂向二阶导数离散相关的系数

以及(akselvoll&moin，1995)

与垂向速度沿垂向二阶导数离散相关的系数为

定义空间离散算子h、l、g、d。其中，h为非线性对流项的离散算子：

l是粘性项的离散拉普拉斯算子：

g是梯度算子，

d是散度算子：

3、本发明半隐式离散(四阶紧致差分)

在时间方向采用adams-bashforth格式离散对流项、crank-nicolson格式离散粘性项，adams-bashforth格式离散涡粘系数项，得到的离散化方程为

其中

与全隐式离散类似，通过block-lu分解将与δp解耦后得到的方程为

或

进一步展开可得

p^n+1/2＝p^n-1/2+δp,

并且

则满足的方程为

最终得到

进一步，将l分解为l＝l1+l2+l3，分别代表对x、y和z的导数离散项，上述方程可得到如下具有二阶时间精度的近似

本发明待求解的x方向动量方程为

分解之后为：

其中，

(1)左端项离散结果

首先离散方程

其中的二阶中心差分离散结果为：

因此

j＝1：

j＝n2-1：

其次离散方程

采用二阶导数的四阶紧致格式

即

其中，

和是n1×n1矩阵。在方程两端乘以可得到

而是循环三对角矩阵，可直接求解从而得到u^***。

最后离散方程

同样采用四阶紧致格式离散

即

其中，

和是n3×n3矩阵。在方程两端乘以可得到

也是循环三对角矩阵，可直接求解得到u^*。

(2)右端项离散结果

首先对

进行离散。

首先采用二阶中心差分格式离散通过线性插值可得到

则的离散结果为：

对于首先插值得到离散点上的y方向速度v，即

再采用中心差分得到

对于通过四阶紧致插值格式得到

即

并且

是n1×n1矩阵，是n1×n1矩阵。通过求解上面的线性方程组得到

然后采用四阶紧致格式离散

即

其中

f1^x和p1^x是n1×n1矩阵。求解上述线性方程组，得到的值。

对于采用四阶紧致差分离散

采用四阶紧致格式离散

即

其中

f1^x和p1^x是n1×n1矩阵。求解上述线性方程组，得到的值。

最后离散首先通过四阶紧致插值格式得到

即

并且

是n3×n3矩阵，是n1×n1矩阵，是n3×n3矩阵，是n1×n1矩阵。通过求解上面的线性方程组得到

而后采用四阶紧致格式离散

即

其中

f1^z是n3×n3矩阵，p1^z是n3×n3矩阵。求解上述线性方程组，得到的值。

对于首先利用四阶紧致格式插值得到离散点处的w，即

以及

离散由涡粘系数产生的导数项

首先离散

得出离散点a和b上的采用四阶紧致格式

即

其中

f1^x和p1^x是n1×n1矩阵。求解上述线性方程组，得到

同样的网格离散离散点上的和νt已知，采用四阶紧致格式离散的一阶导数：

即

其中

f1^x和p1^x是n1×n1矩阵。求解上述线性方程组，得到

离散

首先插值得到a和b点的νt，得到

j＝1:

j＝n2m:

同样的网格，离散

j＝n2m

j＝1

插值出网格点a、b上的νt和的值，再进一步离散可得

离散

首先对νt进行两次四阶插值，得出a,b两点的

即

并且

是n1×n1矩阵，是n3×n3矩阵，是n1×n1矩阵，是n3×n3阵，通过求解上面的线性方程组得到

同样的四阶格式离散a、b点

即

其中

f1^z和p1^z是n3×n3矩阵。求解上述线性方程组，得到

最终的四阶差分离散公式为：

即

其中

f1^z和p1^z是n3×n3矩阵。求解上述线性方程组，得

离散

首先插值a、b点的νt和

即

和是n1×n1矩阵。求解上述线性方程组，得到在a、b点的值。

涡粘系数插值公式为

则最终离散结果为：

最后离散

插值出a、b点上的νt和四阶紧致格式离散

四阶紧致格式插值νt

则最终离散结果为

采用四阶紧致格式离散压力梯度项，即

即

其中

f1^m是n1×n1矩阵，p1^mn1×n1矩阵。求解上述线性方程组得到压力梯度项的四阶离散结果。

4、压力poisson方程

压力poisson方程为

dgδp＝q

其中，右端项

其四阶紧致差分离散结果为

p′和q沿x和z方向的离散傅里叶变换结果为

以及逆变换

首先采用四阶紧致格式离散x方向的压力梯度，即

将用离散傅里叶展开，得到

则在傅里叶空间中，对x方向压力梯度的四阶紧致离散结果可表示为

再采用四阶紧致格式离散x方向压力梯度的散度，即

则在傅里叶空间中，x方向压力梯度的散度的四阶紧致离散结果可表示为

其中，

类似的，z方向压力梯度的散度的四阶紧致离散结果可表示为

其中，

将上述各项傅里叶变换展开式代入压力poisson方程，得到

其中，对的二阶中心差分离散结果为

j＝1(neumann条件)：

j＝n2-1(neumann条件)：

则傅里叶空间内的压力possion方程可以离散为

其线性方程组的系数矩阵为三对角矩阵，采用thomas追赶法求解。求得后，便可通过离散傅里叶逆变换得到p′i,j,k。

5、本发明提供的三对角线性方程组求解方法

(1)三对角方程

对于系数矩阵为三对角矩阵的线性方程组，比如

采用thomas追赶法(tdma算法)求解，本质上仍然是gauss消去法。算法由forwardelimination和backwardsubstitution两步组成。

在forwardelimination步，通过gauss消元将系数矩阵变为一个上三角矩阵，具体步骤为：

在backwardsubstitution步，求解得到

(2)循环三对角方程

在周期边界条件下，会出现系数矩阵为循环三对角矩阵的线性方程组，即

这里给出两种循环三对角方程组的求解方法。

a.解法一

首先将循环三对角矩阵做lu分解，得到

通过等式左右两端对比可知

β1＝b1,γ1＝c1/β1,

βi＝bi-aiγi-1,γi＝ci/βi,i＝2,...,n-1

p1＝cn,q1＝a1/β1

pi＝-pi-1γi-1,qi＝-aiqi-1/βi，i＝2,...,n-2

pn-1＝an-pn-2γn-2,

pn＝bn-(p1q1+p2q2+...+pn-1qn-1).

forwardsubstitution:

lz＝d

z1＝d1/β1,

zi＝(di-aizi-1)/βi,i＝2,...,n-1

backwardsubstitution:

ux＝z,

xn＝zn,

xn-1＝zn-1-qn-1xn,

xi＝zi-γixi+1-qixn,i＝1,...,n-2.

b.解法二

根据sherman-morrison公式，令

则有

a＝a′+uυ^t.

采用thomas追赶法求解下列三对角方程组

a′y＝d,a′q＝u,

则原方程组的解为

6、保证流量守恒的处理方法

假设压力沿流向的变化可以表示为

p＝p′+x·px,

px是流向平均压力梯度。p′满足与p相同的possion方程(边界条件不同)，即

以及

δp＝δp′+x·δpx,

gp＝gp′+px.

流量守恒条件可以表示为

∫vρδu·dv＝0,

表明不同时刻的体积流量相同。

利用

δu＝δu^*-δtgδp,

代入到流量守恒条件中，可得

从而

得到后，n+1时刻的速度更新为

压力更新为

保持展向流量守恒的方法与流向相同。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。