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

本发明涉及计算流体力学中数值计算方法领域，特别涉及一种求解欧拉方程的改进型高阶非线性空间离散方法。
背景技术：
：在高超声速飞行器研制过程中，常规商业软件数值模拟预测得到的气动力/热与实际飞行数据差别较大。其中控制方程中对流项(简化为欧拉方程)的数值离散会直接影响无粘流区域的激波的分辨，并间接影响边界层内流动的预测。当前面向工程问题计算的软件和in-house代码广泛采用的是二阶tvd格式(如nnd格式)，该类数值格式具有良好的数值稳定性。但是tvd类格式存在数值耗散误差过大问题。近年来，加权类非线性格式发展是当下流行的高阶格式之一，广泛应用于空气动力学问题的科学计算中。技术实现要素：针对现有技术中存在的问题，本发明提供了一种与二阶nnd格式相同模板点的新型高阶非线性空间离散方法，通过在下风引入一个三点的子模板，并采用非线性加权策略和激波侦测技术，提高空间数值离散方法的精度，改进非线性离散方法的分辨率，并提高算法计算效率。本发明采用的技术方案如下：一种求解欧拉方程的改进型高阶非线性空间离散方法，包括以下步骤：步骤1、读取初始流场数据，对欧拉方程计算时刻的各节点上的正负通量；步骤2、对各节点上的正负通量进行特征投影，得到特征通量，并根据各节点上的特征通量计算间断侦测因子；步骤3、根据间断侦测因子构造半点上数值通量的高阶混合计算方法，完成欧拉方程的空间离散；步骤4、采用三阶龙格库塔法对时间项进行离散；步骤5、将时间推进至指定tn结束计算，得到tn时刻的流场数据。进一步的，所述步骤1的具体过程为：准备t0时刻流场数据，计算时刻的各节点正负通量具体如下：对三维曲线坐标系(ξ，η，ζ)下的无量纲形式欧拉方程求解：其中，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η).对网格度量系数进行逆变换的雅可比行列式为：无粘矢通量的一般形式记为其中，θ＝kxu+kyv+kzw，k取ξ，η，ζ时分别对应满足：其中，为k方向无粘矢通量雅克比矩阵，经相似变换得到：其中，λk为矩阵的特征值对角阵，rk和lk分别为左、右特征矩阵，λk，rk和lk的具体表现形式分别为：其中，其中，其中，a为当地声速，根据特征传播方向，将无粘通量分为正、负两个部分，一般形式的表达式为：各节点上的正负通量为：进一步的，所述步骤2的具体过程为：对各节点的正负通量进行特征投影得到特征通量计算半点i+1/2附近相邻节点特征通量的差值的绝对大小：计算其中规约化后的最大值θs，以及最大值与最小值的比值其中，“s”表示为5×1矩阵中某一元素，有s＝1，...，5；fr为规约函数，hv(x)为heaviside函数，具体为：计算间断侦测因子其中：σa1＝0.022，σa2＝0.029，σr1＝3.2，σr2＝3.8，rp＝0.7。进一步的，所述步骤3过程如下：对于有：其中，h隐式地定义为：为h±的一个数值近似，对于五点格式有：具体的：步骤31、判断间断侦测因子是否大于σc，若是则对特征通量进行wenn-lc重构，再进行反特征投影，完成重构得到若否则直接进行线性格式重构，得到步骤32、重构之后对各方向空间导数离散项求和，完成欧拉方程的空间离散。进一步的，所述步骤3中，重构得到的具体过程为：先对进行特征投影，在特征空间对特征通量进行重构：其中，其中为重构算子，再进行反特征投影回物理空间：其中半点上的rk，i±1/2和lk，i±1/2采用该半点两侧节点值的算术平均或roe平均计算得到；完成重构：进一步的，所述步骤31中，线性重构的具体方法为：进一步的，所述步骤32具体过程为：根据重构后参数对进行差分离散，带入具体方向即可得到三个方向的空间离散求和并移至欧拉方程右端，得到：完成欧拉方程的空间离散。进一步的，所述步骤4的具体过程为：采用显式三阶tvd龙格-库塔法(r-k)进行时间导数项离散：其中，上标“n”表示第n时刻步的值，上标“n+1”表示第n+1时刻步的值；完成欧拉方程的时空离散。进一步的，所述步骤5的具体过程为：对于空间离散后的欧拉方程，采用r-k时间推进法得到tn+1使刻的流场数据；判断tn+1-tn是否大于ε，是则表示tn+1时刻的流场数据为最终时刻tn的流场数据；否则继续进行时间推进，计算流场数据。与现有技术相比，采用上述技术方案的有益效果为：本发明的改进型高阶非线性空间离散方法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η).逆变换的雅克比行列式如下：无粘矢通量的一般形式记为其中θ＝kxu+kyv+kzw，k取ξ，η，ζ时分别对应由于具有齐次性，故满足：其中为k方向无粘矢通量雅克比矩阵，经相似变换λk为矩阵的特征值对角阵，特征值代表了各类波的特征速度，rk和lk分别为左、右特征矩阵，λk，rk和lk具体形式如下其中，λk为矩阵的特征值对角阵，特征值代表了各类波的特征速度，具体形式为lk和rk分别为左、右特征矩阵，具体形式为其中a为当地声速在数值计算中，根据特征传播方向，通常将无粘通量分为正、负两个部分，以为例：上标“+”表示这部分的无粘矢通量的特征速度指向计算坐标轴的正方向，上标“-”则表示该部分的特征速度指向计算坐标轴的负方向传播。采用steger-warming矢通量分裂，对角阵中每一个对角元素表示为其中带入式(15)为即计算各节点上的正负通量为：下面以ξ方向为例：步骤2中，在ξ方向半点i+1/2处，本发明给出的间断侦测因子计算步骤如下：a.计算各节点上特征通量为：b.计算半点i+1/2附近相邻节点特征通量的差值的绝对大小，以正方向为例，上标省略：c.计算其中规约化后的最大值θs，以及最大值与最小值的比值“s”表示为5×1矩阵中某一元素，有s＝1，...，5。fr为规约函数，hv(x)为heaviside函数，具体为：d.计算间断侦测因子其中参数为：σa1＝0.022，σa2＝0.029，σr1＝3.2，σr2＝3.8，rp＝0.7。步骤3中，对于有其中h隐式地定义为：为h±的一个数值近似，对于五点格式有在数值计算中，为了使算法具有更好的鲁棒性并获得更光滑的结果，在重构时，通常先对进行特征投影，在特征空间对特征变量进行重构：其中为重构算子，如本发明发展的wenn-lc格式重构，最后再“反投影”回物理空间其中半点上的rξ，i±1/2和lξ，i±1/2采用该半点两侧节点值的算术平均或roe平均计算得到。为了捕捉间断，本发明发展的wenn-lc，其中l表示为光滑因子基于拉格朗日插值多项式，c表示中心型格式，格式具体如下：下面以正通量为例(为了方便，下面上标省略)，其数值通量为：其中ri+1/2为半点处的右特征矩阵，为在候选模板s0、s1和s2上(如图2、3)的特征通量，分别为：ωk为格式的非线性权系数，具体为：其中下标“s”表示为5×1矩阵中某一元素，dk为理想权值，并有可以看到，当取理想权情况下，式(31)为标准四阶中心差分格式；αk为非标准化的非线性权重；c为常数，取5.0；另外，ε为一个很小的数，其值为ε＝10-40，以避免在光滑区域分母变为零；isk是当地的光滑度量因子，在x＝xi附近的n点模板{xi，...xi+n-1}上，可以构造拉格朗日插值多项式，基于式(34)，isk的具体形式为：新的全局光滑度量τ4定义如下：可以看到，如果重构算子为线性算子，式(30)与式(28)完全等价。故在在光滑区域，可直接采用如下差分格式：具体实施过程中不需要特征投影-反投影操作。故半点上通量的高阶混合计算方法(hybridwenn-lc，记为hywenn-lc格式)为：其中σc取-0.0687，将式(38)代入式(26)即完成对差分离散，同理另外两方向的空间离散和可类似得到，求和并移至方程右端，得：至此欧拉方程的空间离散结束。步骤4中，为了欧拉方程离散的完整性，下面给出左端时间导数项的离散方法，采用显式三阶tvd龙格库塔方法，形式如下：其中上标“n”表示第n时刻步的值，上标“n+1”表示第n+1时刻步的值。至此，完成了对三维欧拉方程的高阶时空离散。步骤5中，对于空间离散后的欧拉方程，采用r-k时间推进法得到tn+1使刻的流场数据；判断tn+1-tn是否大于ε，是则表示tn+1时刻的流场数据为最终时刻tn的流场数据；否则继续进行时间推进，计算流场数据。本实施例还给出了一二维欧拉方程下来流马赫数ma＝3.0前台阶问题的数值验证，网格为均匀网格，间距δx＝δy＝1/320，计算时间至tn＝4.0。图4-图7和表1给出了相关结果。表1不同空间离散方法的计算耗时比较(单位：秒)计算格式nndwenn-lchywenn-lc耗时(集群并行)0.48683e+03(1.00)0.68592e+03(1.41)0.42024e+03(0.86)本发明并不局限于前述的具体实施方式。本发明扩展到任何在本说明书中披露的新特征或任何新的组合，以及披露的任一新的方法或过程的步骤或任何新的组合。如果本领域技术人员，在不脱离本发明的精神所做的非实质性改变或改进，都应该属于本发明权利要求保护的范围。本说明书中公开的所有特征，或公开的所有方法或过程中的步骤，除了互相排斥的特征和/或步骤以外，均可以以任何方式组合。本说明书中公开的任一特征，除非特别叙述，均可被其他等效或具有类似目的的替代特征加以替换。即，除非特别叙述，每个特征只是一系列等效或类似特征中的一个例子而已。当前第1页12