[发明专利]针对高精度格式的鲁棒高效隐式时间推进方法

权利要求书

1.一种针对高精度格式的鲁棒高效隐式时间推进方法，其特征在于：包括如下步骤：

步骤一、首先给定计算网格，计算每个网格单元的体积；

步骤二、求解修正的内层迭代的虚拟局部时间步长Δτ；

步骤三、将得到的局部时间步长，代入内循环进行迭代求解，直至方程平均残差降两个量级以完成一个物理时间步的推进；

步骤四、重复步骤一至步骤三，直至完成设定的物理时间推进。

2.根据权利要求1所述的针对高精度格式的鲁棒高效隐式时间推进方法，其特征在于：步骤二所述求解修正的内层迭代的虚拟局部时间步长Δτ的方法为：采用如下公式对虚拟局部时间步长Δτ进行修正：

Δτ^*＝max[(Vol_min/Vol_local)^1.0/power,threshold]*Δτ

其中Δτ^*是修正后的虚拟时间步长，Vol_min是最小网格单元体积，Vol_local是当地的网格单元体积，修正指数power的取值范围为3.0～6.0，修正阀值threshold的取值范围是0.01～1.0。

3.根据权利要求2所述的针对高精度格式的鲁棒高效隐式时间推进方法，其特征在于：对于控制单元体积比或长宽比相对较大的网格，取相对较小的power值。

4.根据权利要求2所述的针对高精度格式的鲁棒高效隐式时间推进方法，其特征在于：对于控制单元体积比或长宽比相对较大的网格，取相对较小的threshold值。

5.根据权利要求1所述的针对高精度格式的鲁棒高效隐式时间推进方法，其特征在于：步骤三所述将得到的局部时间步长，代入内循环进行迭代求解的方法为：

第一步、将一般曲线坐标系下的Navier-Stockes方程简写为如下形式：

N=∂Q∂t+RHS(Q)=0]]>

RHS(Q)=-Σi=j,k,l[(Hi+12-Hi-12)+(Hv,i+12-Hv,i-12)]]]>

其中，Q为守恒变量，H和H_v分别为无粘通量和粘性通量，记Navier-Stockes方程的所有余项为R(Q)；

第二步、对Navier-Stockes方程的时间项进行二阶精度三点后差离散，并在每一物理时间步中引入虚拟时间导数，得到如下离散方程：

J3Qn+1,m+1-4Qn+Qn-12Δt+JQn+1,m+1-Qn+1,mΔτ+RHS(Qn+1,m+1)=0]]>

式中，n表示物理时间步，Δt表示物理时间步长，m表示虚拟时间推进步数；

第三步、分别采用迎风和中心格式计算离散方程左端项的无粘通量和粘性通量；

第四步、将离散方程展开得到：

[J·IΔτ+3J·I2Δt+δξ-A++δξ+A-+δη-B++δη+B-+δζ-C++δζ+C--Av,j-1+2Av,j-Av,j+1-Bv,k-1+2Bv,k-Bv,k+1-Cv,l-1+2Cv,l-Cv,l+1]mδQm=RHS(Qm)]]>

其中，A_v,B_v,C_v表示由粘性通量求导得到的系数矩阵，δ⁺和δ^-是迎风算子；令D,分别是对角线矩阵、下三角矩阵和上三角矩阵，

L~=-Aj-1+-Bk-1+-Cl-1+-Av,j-1-Bv,k-1-Cv,l-1]]>

D=J.IΔτ+3J·I2Δt+Aj+-Aj-+Bk+-Bk-+Cl+-Cl-+2Av,j+2Bv,k+2Cv,l]]>

U~=Aj+1-+Bk+1-+Cl+1--Av,j+1-Bv,k+1-Cv,l+1]]>

第五步、令则得到如下近似方程：

LD^-1UδQ^m＝RHS(Q^m)

第六步、取

A±=12(A±ρ~(A)I)]]>

B±=12(B±ρ~(B)I)]]>

C±=12(C±ρ~(C)I)]]>

ρ~=ρ~(A)+ρ~(B)+ρ~(C)]]>

则得到：

L=D-Aj-1+-Bk-1+-Cl-1+-Av,j-1-Bv,k-1-Cv,l-1]]>

D=ρ~I+J·IΔτ+2Av,j+2Bv,k+2Cv,l]]>

U=D+Aj+1-+Bk+1-+Cl+1--Av,j+1-Bv,k+1-Cv,l+1]]>

式中：表示A,B,C的谱半径,|A_v,j|,|B_v,k|,|C_v,l|表示这三个粘性矩阵的范数；此时的对角块矩阵简化为一个标量与单位矩阵的乘积；

第七步、将第五步所述近似方程分解成三个算子的形式：

L(δQ₁)＝RHS(Q)

D^-1(δQ₂)＝δQ₁

U(δQ)＝δQ₂

其中，第一个和第三个算子用追赶法求解，第二个算子的求解采用标量求逆。