[发明专利]一种适用于复杂流动的高效高精度翼型绕流数值模拟方法

权利要求书

1.一种适用于复杂流动的高效高精度翼型绕流数值模拟方法，其特征在于：包括以下步骤：

步骤1：根据设计任务要求，构建翼型模型，并建立计算网格，获得需要的网格单元信息；

步骤2：在步骤1建立的网格空间中构造空间离散形式的半离散控制方程：

当采用欧拉方程时，其微分形式的控制方程为

其中，

q是守恒形式的流场变量，f和g表示x和y方向的通量，ρ，u，v，p，E分别表示流体密度、x方向速度，y方向速度，压强和能量；

对通量项进行空间离散得到：

其中i，j是单元节点编号；f_i+1/2，j和g_i，j+1/2分别为x方向和y方向的界面数值通量；

步骤3、采用HWENO重构方案构造重构多项式

首先对控制方程分别沿x和y进行求导，求导后的控制方程转化为：

其中，f_x(q，r)＝f′(q)r，g_x(q，r)＝g′(q)r，f_y(q，s)＝f′(q)s，g_y(q，s)＝g′(q)s，r和s分别是变量q关于x和y的导数，

选择含有V₀～V₈在内共9个单元的模板，并将其分为8个子模板：

S₁＝{V₀，V₁，V₂，V₈}，S₂＝{V₀，V₂，V₃，V₄}，S₃＝{V₀，V₄，V₅，V₆}，S₄＝{V₀，V₆，V₇，V₈}S₅＝{V₀，V₁，V₂，V₃，V₇，V₈}，S₆＝{V₀，V₁，V₂，V₃，V₄，V₅}，S₇＝{V₀，V₃，V₄，V₅，V₆，V₇}，S₈＝{V₀，V₁，V₅，V₆，V₇，V₈}

其中V₀～V₈分别指网格单元V_i，j、V_i+1，j+1、V_i，j+1、V_i-1，j+1、V_i-1，j、V_i-1，j-1、V_i，j-1、V_i+1，j-1、V_i+1，j；

在子模板S₁，S₂，S₃，S₄中，关于任意变量q的插值多项式p_n需要满足如下约束条件：

其中，

n＝1，k＝0，1，2，8，k_x＝8，k_y＝2；n＝2，k＝0，2，3，4，k_x＝4，k_y＝2；

n＝3，k＝0，4，5，6，k_x＝4，k_y＝6；n＝4，k＝0，6，7，8，k_x＝8，k_y＝6.

在子模板S₅，S₆，S₇，S₈中，关于任意变量q的插值多项式p_n需要满足如下约束条件：

其中，

n＝5，k＝0，1，2，3，7，8；n＝6，k＝0，1，2，3，4，5；

n＝7，k＝0，3，4，5，6，7；n＝8，k＝0，1，5，6，7，8.

在每个子模板中，具有三阶精度的插值多项式写为：

p_n(x，y)＝a₀+a₁(x-x₀)+a₂(y-y₀)+a₃(x-x₀)(y-y₀)+a₄(x-x₀)²+a₅(y-y₀)²，n＝1，2，3，4，5，6，7，8

将插值多项式代入约束条件，从而在每个子模板上得到一组关于多项式系数a_k(k＝0，1，2，3，4，5)的线性代数方程组，求解该方程组得到各子模板中插值多项式p_n的系数a_k；

在得到各子模板上的插值多项式后，采用WENO限制器的方法，通过光滑指示因子求得9个多项式的权重，加权组合成最终的重构多项式P_i，j(x，y)：

光滑指示因子定义如下：

其中|α|＝α₁+α₂，根据光滑指示因子得到每个多项式的权重如下：

最后通过加权组合的方式得到单元V_i,j上的空间重构多项式：

步骤4、根据步骤3得到的最终的重构多项式求解多维黎曼求解器所需的重构状态量；

对于单元V_i,j，首先通过最终的重构多项式求得变量在包含界面中点和角点在内的8个插值点(x_i+1,j，y_i+1/2,j)，(x_i,j，y_i+1,j)，(x_i-1,j，y_i,j)，(x_i,j，y_i-1,j)，(x_i+1,j，y_i+1/2,j)，(x_i-1/2,j，y_i+1/2,j)，(x_i-1/2,j，y_i-1/2,j)，(x_i+1/2,j，y_i-1/2,j)处的重构值，继而求得界面i+1/2处中点以及上下角点处的多维黎曼求解器所需的状态量其中上标“R”和“L”分别表示界面中点两侧的重构变量值，上标“RU”,“LU”,“LD”和“RD”表示角点四周的重构变量值，下标“i+1/2，j”表示界面中点，下标“i+1/2，j+1/2”和“i+1/2，j-1/2”分别表示界面的上下角点；

步骤5、采用多维黎曼求解器进行界面通量求解：

在步骤2的半离散控制方程中，界面通量的具体求解公式如下：

其中，ω₁＝1/6,ω₂＝4/6,ω₃＝1/6为权重系数，分别为x方向和y方向上的辛普森插值点，在x方向上，y_j1＝y_j+1/2,分别表示了界面i+1/2的上角点、中点和下角点；辛普森插值点处的数值通量为：

其中，通过经典的一维HLLE格式求得：

其中，上表“m”表示与界面中点相关的物理量，下标“R”和“L”分别表示界面两侧的重构变量值，由步骤2中求解得到；和分别表示左右传播的最大波速，采用如下公式进行计算：

a是声速，上标“～”表示Roe平均；

界面角点处的通量和则通过Balsara的真正二维HLLE格式求得；

步骤6、根据界面通量求解残差，并将半离散有限体积格式转化为时空全离散有限体积格式，全流场进行时间推进求解，得到最终的流场解。