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

摘要

本发明提供了一种求解欧拉方程的改进型高阶非线性空间离散方法，包括以下步骤：步骤1、读取初始流场数据，对欧拉方程计算时刻的各节点上的正负通量；步骤2、对各节点上的正负通量进行特征投影，得到特征通量，并根据各节点上的特征通量计算间断侦测因子；步骤3、根据间断侦测因子构造半点上数值通量的高阶混合计算方法，完成欧拉方程的空间离散；步骤4、采用三阶龙格库塔法对时间项进行离散；步骤5、将时间推进至指定t

说明书

技术领域

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

背景技术

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

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

发明内容

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

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

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

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

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

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

步骤5、将时间推进至指定t

进一步的，所述步骤1的具体过程为：准备t

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

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

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

其中，y为比热比；ξ

ξ

η

ζ

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

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

其中，θ＝k

其中，

其中，Λ

其中，

其中，

其中，a为当地声速，

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

各节点上的正负通量为：

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

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

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

计算其中规约化后

其中，“s”表示为5×1矩阵中某一元素，有s＝1，...，5；f

计算间断侦测因子

其中：σ

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

其中，h隐式地定义为：

具体的：

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

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

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

先对

其中，其中

其中半点上的R

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

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

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

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

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

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

进一步的，所述步骤5的具体过程为：对于空间离散后的欧拉方程，采用R-K时间推进 法得到t

与现有技术相比，采用上述技术方案的有益效果为：本发明的改进型高阶非线性空间 离散方法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为比热比；ξ

ξ

η

ζ

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

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

其中θ＝k

其中

Λ

其中，

Λ

L

其中

a为当地声速

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

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

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

其中

带入式(15)为

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

下面以ξ方向为

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

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

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

c.计算其中规约化后

“s”表示为5×1矩阵中某一元素，有s＝1，...，5。f

d.计算间断侦测因子

其中参数为：σ

步骤3中，

对于

其中h隐式地定义为：

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

其中

其中半点上的R

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

下面以正通量

其中R

ω

其中下标“s”表示为5×1矩阵中某一元素，d

基于式(34)，IS

新的全局光滑度量τ

可以看到，如果重构算子

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

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

故半点上通量

其中σ

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

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

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

步骤5中，对于空间离散后的欧拉方程，采用R-K时间推进法得到t

本实施例还给出了一二维欧拉方程下来流马赫数Ma＝3.0前台阶问题的数值验证，网格 为均匀网格，间距Δx＝Δy＝1/320，计算时间至t

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

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

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

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