CN112214869B - 一种求解欧拉方程的改进型高阶非线性空间离散方法 - Google Patents

Abstract

本发明提供了一种求解欧拉方程的改进型高阶非线性空间离散方法，包括以下步骤：步骤1、读取初始流场数据，对欧拉方程计算时刻的各节点上的正负通量；步骤2、对各节点上的正负通量进行特征投影，得到特征通量，并根据各节点上的特征通量计算间断侦测因子；步骤3、根据间断侦测因子构造半点上数值通量的高阶混合计算方法，完成欧拉方程的空间离散；步骤4、采用三阶龙格库塔法对时间项进行离散；步骤5、将时间推进至指定t_N结束计算，得到t_N时刻的流场数据。本发明的改进型高阶非线性空间离散方法WENN‑LC格式在同等网格下较传统NND格式具有更高的流动结构分辨率；此外，hyWENN‑LC混合格式不仅具有更高的分辨率，同时具有更快的计算效率。

Description

一种求解欧拉方程的改进型高阶非线性空间离散方法

技术领域

本发明涉及计算流体力学中数值计算方法领域，特别涉及一种求解欧拉方程的改进型 高阶非线性空间离散方法。

背景技术

在高超声速飞行器研制过程中，常规商业软件数值模拟预测得到的气动力/热与实际飞 行数据差别较大。其中控制方程中对流项(简化为欧拉方程)的数值离散会直接影响无粘 流区域的激波的分辨，并间接影响边界层内流动的预测。

当前面向工程问题计算的软件和In-house代码广泛采用的是二阶TVD格式(如NND格式)，该类数值格式具有良好的数值稳定性。但是TVD类格式存在数值耗散误差过大问题。近年来，加权类非线性格式发展是当下流行的高阶格式之一，广泛应用于空气动力学问题的科学计算中。

发明内容

针对现有技术中存在的问题，本发明提供了一种与二阶NND格式相同模板点的新型高 阶非线性空间离散方法，通过在下风引入一个三点的子模板，并采用非线性加权策略和激 波侦测技术，提高空间数值离散方法的精度，改进非线性离散方法的分辨率，并提高算法 计算效率。

本发明采用的技术方案如下：一种求解欧拉方程的改进型高阶非线性空间离散方法，包 括以下步骤：

步骤1、读取初始流场数据，对欧拉方程计算时刻的各节点上的正负通量；

步骤2、对各节点上的正负通量进行特征投影，得到特征通量，并根据各节点上的特征 通量计算间断侦测因子；

步骤3、根据间断侦测因子构造半点上数值通量的高阶混合计算方法，完成欧拉方程的 空间离散；

步骤4、采用三阶龙格库塔法对时间项进行离散；

步骤5、将时间推进至指定t_N结束计算，得到t_N时刻的流场数据。

进一步的，所述步骤1的具体过程为：准备t₀时刻流场数据，计算时刻的各节点正负通 量

具体如下：

对三维曲线坐标系(ξ，η，ζ)下的无量纲形式欧拉方程求解：

其中，Q为守恒变量，E、F、G为直角坐标系(x，y，z)下的无粘矢通量，具体表达式为：

其中，ρ、u≡(u，v，w)、p分别表示密度、速度矢量、压力；e为总能，具体表达式为：

其中，y为比热比；ξ_x，ξ_y，ξ_z，η_x，η_y，η_z，ζ_x，ζ_y，ζ_z称为网格度量系数，具体求解方法为：

ξ_x＝J(y_ηz_ζ-z_ηy_ζ)，ξ_y＝J(z_ηx_ζ-x_ηz_ζ)，ξ_z＝J(x_ηy_ζ-y_ηx_ζ)，

η_x＝J(y_ζz_ξ-z_ζy_ξ)，η_y＝J(z_ζx_ξ-x_ζz_ξ)，η_z＝J(x_ζy_ξ-y_ζx_ξ)，

ζ_x＝J(y_ξz_η-z_ξy_η)，ζ_y＝J(z_ξx_η-x_ξz_η)，ζ_z＝J(x_ξy_η-y_ξx_η).

对网格度量系数进行逆变换的雅可比行列式为：

无粘矢通量的一般形式记为

其中，θ＝k_xu+k_yv+k_zw，k取ξ，η，ζ时分别对应

满足：

其中，

为k方向无粘矢通量

雅克比矩阵，经相似变换得到：

其中，Λ_k为矩阵

的特征值对角阵，R_k和L_k分别为左、右特征矩阵，

Λ_k，R_k和L_k的 具体表现形式分别为：

其中，

其中，

其中，a为当地声速，

根据特征传播方向，将无粘通量分为正、负两个部分，一般形式的表达式为：

各节点上的正负通量为：

进一步的，所述步骤2的具体过程为：

对各节点的正负通量进行特征投影得到特征通量

计算半点i+1/2附近相邻节点特征通量的差值的绝对大小：

计算其中规约化后

的最大值θ_s，以及

最大值与最小值的比值

其中，“s”表示为5×1矩阵中某一元素，有s＝1，...，5；f_r为规约函数，H_v(x)为Heaviside 函数，具体为：

计算间断侦测因子

其中：σ_a1＝0.022，σ_a2＝0.029，σ_r1＝3.2，σ_r2＝3.8，r_p＝0.7。

进一步的，所述步骤3过程如下：对于

有：

其中，h隐式地定义为：

为h^±的一个数值近似，对于五点格式有：

具体的：

步骤31、判断间断侦测因子

是否大于σ_c，若是则对特征通量

进行WENN-LC重构， 再进行反特征投影，完成重构得到

若否则直接进行线性格式重构，得到

步骤32、重构之后对各方向空间导数离散项求和，完成欧拉方程的空间离散。

进一步的，所述步骤3中，重构得到

的具体过程为：

先对

进行特征投影，在特征空间对特征通量

进行重构：

其中，其中

为重构算子，再进行反特征投影回物理空间：

其中半点上的R_k，i±1/2和L_k，i±1/2采用该半点两侧节点值的算术平均或Roe平均计算得到； 完成重构：

进一步的，所述步骤31中，线性重构的具体方法为：

进一步的，所述步骤32具体过程为：根据重构后参数对

进行差分离散，带入具 体方向即可得到三个方向的空间离散

求和并移至欧拉方程右端， 得到：

完成欧拉方程的空间离散。

进一步的，所述步骤4的具体过程为：

采用显式三阶TVD龙格-库塔法(R-K)进行时间导数项离散：

其中，上标“n”表示第n时刻步的值，上标“n+1”表示第n+1时刻步的值；完成欧 拉方程的时空离散。

进一步的，所述步骤5的具体过程为：对于空间离散后的欧拉方程，采用R-K时间推进 法得到t_n+1使刻的流场数据；判断t_n+1-t_N是否大于ε，是则表示t_n+1时刻的流场数据为最终 时刻t_N的流场数据；否则继续进行时间推进，计算流场数据。

与现有技术相比，采用上述技术方案的有益效果为：本发明的改进型高阶非线性空间 离散方法WENN-LC格式在同等网格下较传统NND格式具有更高的流动结构分辨率；此外， hyWENN-LC混合格式不仅具有更高的分辨率，同时具有更快的计算效率。

附图说明

图1是本发明的求解欧拉方程改进型高阶非线性空间离散方法的流程图。

图2是本发明中的改进型非线性WENN-LC格式和hyWENN-LC格式的正通量模板示意图。

图3是本发明中的改进型非线性WENN-LC格式和hyWENN-LC格式的负通量模板示意图。

图4是本发明一实施例中标准算例数值验证：二维欧拉方程前台阶问题密度分布图(NND 格式)。

图5是本发明一实施例中标准算例数值验证：二维欧拉方程前台阶问题密度分布图 (WENN-LC格式)。

图6是本发明一实施例中标准算例数值验证：二维欧拉方程前台阶问题密度分布图 (hyWENN-LC格式)。

图7是本发明一实施例中前台阶问题间断侦测因子

分布图(黑色部分为

区域)。

具体实施方式

下面结合附图对本发明做进一步描述。

如图1所示，本发明提供了一种求解欧拉方程的改进型高阶非线性空间离散方法，包 括以下步骤：

步骤1、读取初始流场数据，对欧拉方程计算时刻的各节点上的正负通量；

步骤2、对各节点上的正负通量进行特征投影，得到特征通量，并根据各节点上的特征 通量计算间断侦测因子；

步骤3、根据间断侦测因子构造半点上数值通量的高阶混合计算方法，完成欧拉方程的 空间离散；

步骤4、采用三阶龙格库塔法对时间项进行离散；

步骤5、将时间推进至指定tN结束计算，得到tN时刻的流场数据。

具体的，

步骤1中，以三维曲线坐标系(ξ，η，ζ)下的无量纲形式欧拉方程为例：

其中Q为守恒变量，E、F、G为直角坐标系(x，y，z)下的无粘矢通量，具体表达式为：

其中ρ，u≡(u，v，w)和p分别表示为密度、速度矢量、压力；e为总能，有如下关系：

其中y为比热比；ξ_x，ξ_y，ξ_z，η_x，η_y，η_z，ζ_x，ζ_y，ζ_z称为网格度量系数，求法如下：

ξ_x＝J(y_ηz_ζ-z_ηy_ζ)，ξ_y＝J(z_ηx_ζ-x_ηz_ζ)，ξ_z＝J(x_ηy_ζ-y_ηx_ζ)， (4)

η_x＝J(y_ζz_ξ-z_ζy_ξ)，η_y＝J(z_ζx_ξ-x_ζz_ξ)，η_z＝J(x_ζy_ξ-y_ζx_ξ)，

ζ_x＝J(y_ξz_η-z_ξy_η)，ζ_y＝J(z_ξx_η-x_ξz_η)，ζ_z＝J(x_ξy_η-y_ξx_η).

逆变换的雅克比行列式如下：

无粘矢通量的一般形式记为

其中θ＝k_xu+k_yv+k_zw，k取ξ，η，ζ时分别对应

由于

具有齐次性，故满足：

其中

为k方向无粘矢通量

雅克比矩阵，经相似变换

Λ_k为矩阵

的特征值对角阵，特征值代表了各类波的特征速度，R_k和L_k分别为左、右特 征矩阵，

Λ_k，R_k和L_k具体形式如下

其中，

Λ_k为矩阵

的特征值对角阵，特征值代表了各类波的特征速度，具体形式为

L_k和R_k分别为左、右特征矩阵，具体形式为

其中

a为当地声速

在数值计算中，根据特征传播方向，通常将无粘通量分为正、负两个部分，以

为例：

上标“+”表示这部分的无粘矢通量的特征速度指向计算坐标轴的正方向，上标“-”则表示该部 分的特征速度指向计算坐标轴的负方向传播。

采用Steger-Warming矢通量分裂，对角阵中每一个对角元素表示为

其中

带入式(15)为

即计算各节点上的正负通量为：

下面以ξ方向为

例：

步骤2中，在ξ方向半点i+1/2处，本发明给出的间断侦测因子

计算步骤如下：

a.计算各节点上特征通量为：

b.计算半点i+1/2附近相邻节点特征通量的差值的绝对大小，以正方向为例，上标省略：

c.计算其中规约化后

的最大值θ_s，以及

最大值与最小值的比值

“s”表示为5×1矩阵中某一元素，有s＝1，...，5。f_r为规约函数，H_v(x)为Heaviside函数，具体 为：

d.计算间断侦测因子

其中参数为：σ_a1＝0.022，σ_a2＝0.029，σ_r1＝3.2，σ_r2＝3.8，r_p＝0.7。

步骤3中，

对于

有

其中h隐式地定义为：

为h^±的一个数值近似，对于五点格式有

在数值计算中，为了使算法具有更好的鲁棒性并获得更光滑的结果，在重构

时，通常先 对

进行特征投影，在特征空间对特征变量

进行重构：

其中

为重构算子，如本发明发展的WENN-LC格式重构，最后再“反投影”回物理空间

其中半点上的R_ξ，i±1/2和L_ξ，i±1/2采用该半点两侧节点值的算术平均或Roe平均计算得到。

为了捕捉间断，本发明发展的WENN-LC，其中L表示为光滑因子基于拉格朗日插值多 项式，C表示中心型格式，格式具体如下：

下面以正通量

为例(为了方便，下面上标省略)，其数值通量为：

其中R_i+1/2为半点处的右特征矩阵，

为在候选模板S₀、S₁和S₂上(如图2、3)的特征通量， 分别为：

ω_k为格式的非线性权系数，具体为：

其中下标“s”表示为5×1矩阵中某一元素，d_k为理想权值，并有

可以看到， 当取理想权情况下，式(31)为标准四阶中心差分格式；α_k为非标准化的非线性权重；C为 常数，取5.0；另外，ε为一个很小的数，其值为ε＝10^-40，以避免在光滑区域分母变为零；IS_k是当地的光滑度量因子，在x＝x_i附近的n点模板{x_i，...x_i+n-1}上，可以构造拉格朗日插值多项式，

基于式(34)，IS_k的具体形式为：

新的全局光滑度量τ₄定义如下：

可以看到，如果重构算子

为线性算子，式(30)与式(28)完全等价。

故在在光滑区域，可直接采用如下差分格式：

具体实施过程中不需要特征投影-反投影操作。

故半点上通量

的高阶混合计算方法(hybridWENN-LC，记为hyWENN-LC格式)为：

其中σ_c取-0.0687，将式(38)代入式(26)即完成对

差分离散，同理另外两方向的空间离 散

和

可类似得到，求和并移至方程右端，得：

至此欧拉方程的空间离散结束。

步骤4中，为了欧拉方程离散的完整性，下面给出左端时间导数项的离散方法，采用显 式三阶TVD龙格库塔方法，形式如下：

其中上标“n”表示第n时刻步的值，上标“n+1”表示第n+1时刻步的值。至此，完 成了对三维欧拉方程的高阶时空离散。

步骤5中，对于空间离散后的欧拉方程，采用R-K时间推进法得到t_n+1使刻的流场数据； 判断t_n+1-t_N是否大于ε，是则表示t_n+1时刻的流场数据为最终时刻t_N的流场数据；否则继续 进行时间推进，计算流场数据。

本实施例还给出了一二维欧拉方程下来流马赫数Ma＝3.0前台阶问题的数值验证，网格 为均匀网格，间距Δx＝Δy＝1/320，计算时间至t_N＝4.0。图4-图7和表1给出了相关结果。

表1不同空间离散方法的计算耗时比较(单位：秒)

计算格式	NND	WENN-LC	hyWENN-LC
				耗时(集群并行)	0.48683E+03(1.00)	0.68592E+03(1.41)	0.42024E+03(0.86)

本发明并不局限于前述的具体实施方式。本发明扩展到任何在本说明书中披露的新特 征或任何新的组合，以及披露的任一新的方法或过程的步骤或任何新的组合。如果本领域 技术人员，在不脱离本发明的精神所做的非实质性改变或改进，都应该属于本发明权利要 求保护的范围。

本说明书中公开的所有特征，或公开的所有方法或过程中的步骤，除了互相排斥的特 征和/或步骤以外，均可以以任何方式组合。

本说明书中公开的任一特征，除非特别叙述，均可被其他等效或具有类似目的的替代 特征加以替换。即，除非特别叙述，每个特征只是一系列等效或类似特征中的一个例子而 已。

Claims (7)

1.一种求解欧拉方程的改进型高阶非线性空间离散方法，其特征在于，包括以下步骤：

步骤1、读取初始流场数据，对欧拉方程计算时刻的各节点上的正负通量；

步骤2、对各节点上的正负通量进行特征投影，得到特征通量，并根据各节点上的特征通量计算间断侦测因子；

步骤3、根据间断侦测因子构造半点上数值通量的高阶混合计算方法，完成欧拉方程的空间离散；

步骤4、采用三阶龙格库塔法对时间项进行离散；

步骤5、将时间推进至指定t_N结束计算，得到t_N时刻的流场数据；

所述步骤1的具体过程为：准备t₀时刻流场数据，计算时刻的各节点正负通量

具体如下：

对三维曲线坐标系ξ,η,ζ下的无量纲形式欧拉方程求解：

其中，Q为守恒变量，E、F、G为直角坐标系x,y,z下的无粘矢通量，具体表达式为：

其中，ρ、u≡(u,v,w)、p分别表示密度、速度矢量、压力；e为总能，具体表达式为:

其中，γ为比热比；ξ_x,ξ_y,ξ_z,η_x,η_y,η_z,ζ_x,ζ_y,ζ_z称为网格度量系数，具体求解方法为：

ξ_x＝J(y_ηz_ζ-z_ηy_ζ),ξ_y＝J(z_ηx_ζ-x_ηz_ζ),ξ_z＝J(x_ηy_ζ-y_ηx_ζ),

η_x＝J(y_ζz_ξ-z_ζy_ξ),η_y＝J(z_ζx_ξ-x_ζz_ξ),η_z＝J(x_ζy_ξ-y_ζx_ξ),

ζ_x＝J(y_ξz_η-z_ξy_η),ζ_y＝J(z_ξx_η-x_ξz_η),ζ_z＝J(x_ξy_η-y_ξx_η).

对网格度量系数进行逆变换的雅可比行列式为：

无粘矢通量的一般形式记为

其中，θ＝k_xu+k_yv+k_zw，k取ξ,η,ζ时分别对应

满足：

其中，

为k方向无粘矢通量

雅克比矩阵，经相似变换得到：

其中，Λ_k为矩阵

的特征值对角阵，R_k和L_k分别为左、右特征矩阵，

Λ_k,R_k和L_k的具体表现形式分别为：

其中，

其中，

其中，a为当地声速，

根据特征传播方向，将无粘通量分为正、负两个部分，一般形式的表达式为：

各节点上的正负通量为：

所述步骤3过程如下：对于

有:

其中h隐式地定义为：

为h^±的一个数值近似，对于五点格式有:

具体的：

步骤31、判断间断侦测因子

是否大于σ_c，若是则对特征通量

进行WENN-LC重构，再进行反特征投影，完成重构得到

否则直接进行线性格式重构，得到

步骤32、重构之后对各方向空间导数离散项求和，完成欧拉方程的空间离散。

2.根据权利要求1所述的求解欧拉方程的改进型高阶非线性空间离散方法，其特征在于，所述步骤2的具体过程为：

对各节点的正负通量进行特征投影得到特征通量

计算半点i+1/2附近相邻节点特征通量的差值的绝对大小：

计算其中规约化后

的最大值θ_s，其中l＝1,2,3，以及

最大值与最小值的比值θ_s：

其中，“s”表示为5×1矩阵中某一元素，有s＝1,…,5；f_r为规约函数，H_v(x)为Heaviside函数，具体为：

计算间断侦测因子

其中：σ_a1＝0.022,σ_a2＝0.029,σ_r1＝3.2,σ_r2＝3.8,r_p＝0.7。

3.根据权利要求1所述的求解欧拉方程的改进型高阶非线性空间离散方法，其特征在于，所述步骤31中，重构得到

的具体过程为：

先对

进行特征投影，在特征空间对特征通量

进行重构，其中，l＝i-2,...,i+2：

其中，其中

为重构算子，再进行反特征投影回物理空间：

其中半点上的R_k,i±1/2和L_k,i±1/2采用该半点两侧节点值的算术平均或Roe平均计算得到；完成重构：

4.根据权利要求1或3所述的求解欧拉方程的改进型高阶非线性空间离散方法，其特征在于，所述步骤31中，线性重构的具体方法为：

5.根据权利要求1所述的求解欧拉方程的改进型高阶非线性空间离散方法，其特征在于，所述步骤32具体过程为：根据重构后参数对

进行差分离散，带入具体方向即可得到三个方向的空间离散

求和并移至欧拉方程右端，得到：

完成欧拉方程的空间离散。

6.根据权利要求5所述的求解欧拉方程的改进型高阶非线性空间离散方法，其特征在于，所述步骤4的具体过程为：

采用显式三阶TVD龙格库塔法进行时间导数项离散：

其中，上标“n”表示第n时刻步的值，上标“n+1”表示第n+1时刻步的值；完成欧拉方程的时空离散。

7.根据权利要求6所述的求解欧拉方程的改进型高阶非线性空间离散方法，其特征在于，所述步骤5的具体过程为：对于空间离散后的欧拉方程，采用R-K时间推进法得到t_n+1使刻的流场数据；判断t_n+1-t_N是否大于ε，是则表示t_n+1时刻的流场数据为最终时刻t_N的流场数据；否则继续进行时间推进，计算流场数据。

Images