耦合流动及网格分布的重构限制器进行激波捕捉的方法与流程

1.本发明涉及计算流体力学中的激波捕捉技术领域，更为具体的，涉及一种耦合流动及网格分布的重构限制器进行激波捕捉的方法。


背景技术：

2.当前计算流体力学（cfd, computational fluid dynamics）在流体力学研究中的地位愈发重要，基于二阶精度的cfd方法已在飞行器气动性能分析与优化设计、气动热分析与评估、发动机燃烧、复杂流动机理分析、气动噪声、流动控制等方面发挥越来越主要的作用。但是，二阶精度计算格式数值耗散、色散较大，在网格规模受计算机条件限制时通常难以给出旋涡、分离、湍流等复杂流动的精细结构，对于一些关键的气动参数（如阻力、力矩和热流等）不能给出准确的计算结果，无法满足飞行器精细化设计需求。而高精度方法数值耗散和色散小，分辨率高，能够更精细模拟流场，因此各类高精度方法得到了广泛的研究和应用。
3.现有技术中，高精度间断galerkin有限元方法（dg, discontinuous galerkin method）作为一种高精度方法，具有众多的优点。它在数值模拟强间断流场时，有效的激波捕捉技术包括限制器、人工粘性及滤波器。但是dg方法中的这三种激波捕捉技术在数值模拟强间断时存在鲁棒性差、主控方程右端项残差（后称残差）下降量级少，收敛滞止的问题。
4.hweno限制器作为一类重构限制器，用于重构网格单元内变量分布多项式，同样存在上述问题，如何有效降低残差，抑制数值振荡是计算流体力学中的一个重要研究内容。
5.hweno限制器在数值模拟中构造单元自身模板的变量分布多项式以及单元的面相邻单元模板的变量分布多项式，从而形成一系列候选模板的变量分布多项式（后称变量分布多项式），见图1。随后基于这些候选模板的变量分布多项式，通过hweno加权重构过程形成单元新的变量分布多项式，从而调整单元内密度、速度、能量变量的分布。hweno限制器的重构过程关系到整个方法的鲁棒性和残差收敛性，许多研究人员在这方面开展了大量的工作。
6.hweno限制器的关键技术是如何建立候选模板多项式和如何加权候选模板多项式，如何选择加权权重系数关系到该限制器限制过程计算鲁棒性以及残差收敛性。传统的激波捕捉技术用到的加权权重系数计算方法为固定自身模板变量分布多项式与面相邻模板变量分布多项式的比例，比如1000:1，这样突出自身模板在重构过程的比重，从而保证限制器限制过程的计算精度，但是在激波区域可能带来数值振荡。这种固定加权权重系数选取方法未考虑流动特征与网格分布的耦合，在重构过程中仅考虑了变量分布多项式的一阶导数项以及高阶导数项构成的光滑因子影响，见具体实施方式中的公式（4），导致数值模拟复杂流场及复杂网格时右端项残差的振荡、残差下降量少，甚至计算发散。
7.基于以上原因，需要探索hweno限制器中的加权权重系数选取方法，实现该限制器在包含强激波流场的数值模拟过程中能减少残差的振荡，同时降低残差的绝对值，实现激波稳定捕捉。


技术实现要素：

8.本发明的目的在于克服现有技术的不足，提供一种耦合流动及网格分布的重构限制器进行激波捕捉的方法，通过耦合局部区域密度梯度以及局部区域的网格分布，构建重构限制器的自适应加权权重系数，以此自适应标记重构限制器中每个重构模板的激波强度，降低激波强度大且模板方向与密度梯度方向角度大的重构模板的重构权重，实现激波稳定捕捉，从而抑制激波带来的数值振荡，能够给出关键气动参数（如阻力、力矩和热流等）可靠的计算结果，获取高分辨率的流场结构，满足飞行器精细化设计需求，适用于包含强间断流场以及在网格规模受计算机条件限制时的场景。
9.本发明的目的是通过以下方案实现的：一种耦合流动及网格分布的重构限制器进行激波捕捉的方法，包括步骤：s1，以密度作为重构限制器中定义自适应加权权重系数的物理量，在第n时间迭代步计算所有网格单元密度分布多项式，并探测出包含激波的问题单元，n为正整数；s2，定义重构限制器中重构模板构造方式；s3，按照步骤s2定义的重构模板构造方式，重构获得每个重构模板内的密度分布多项式，该多项式为重构限制器中的重构分量，用于重构单元最终的密度分布多项式实现激波捕捉；s4，计算步骤s3获得的每个重构模板的密度分布多项式的光滑因子，同时计算重构模板的密度分布多项式的密度梯度，该密度梯度将作为重构模板内是否包含激波的基本依据，并计算问题单元到面相邻单元之间的体心矢量；s5，将重构模板的密度梯度投影到该重构模板的体心矢量，实现当地流动特征与网格分布的耦合，以投影后的密度梯度计算自适应加权权重系数，然后依据计算的自适应加权权重系数以及步骤s4计算的光滑因子计算每个重构模板密度分布多项式的权重因子，该权重因子代表每个重构模板中激波强度；s6，根据步骤s3获得的每个重构模板的密度分布多项式以及步骤s5计算的每个重构模板密度分布多项式的权重因子，加权计算得到问题单元内密度变量最终的密度分布多项式，实现重构限制器对对单元密度变量的重构，实现激波捕捉。
10.进一步地，在步骤s1中，设通过问题单元探测器探测到0单元为问题单元。
11.进一步地，在步骤s2中，所述定义重构模板构造方式包括子步骤：将0,1,2,3表示单元体心编号，a、b、c、d、e、f表示单元的节点编号，在时间迭代步第n步时获得0单元及其面相邻单元1、2、3单元内的密度变量分布多项式分布表示为：。
12.进一步地，在步骤s3中，包括子步骤：构造4个候选模板，分别是0-0单元模板、0-1单元模板、0-2单元模板、0-3单元模板，并将模板名分别定义为，重构0单元内的密度变量，则四个模板对应
的密度分布多项式分别如下：四个模板的密度分布多项式满足以下条件：式中，表示0单元的体积，代表笛卡尔坐标系的三个方向。
13.进一步地，在步骤s4中，采用如下公式计算步骤s3获得的每个候选模板的密度分布多项式的光滑因子：式中，表示光滑因子。
14.进一步地，计算重构模板的密度分布多项式的梯度以及问题单元到面相邻单元之间的体心矢量包括子步骤：采用如下公式计算模板的密度分布多项式的梯度以及0单元到面相邻单元1,2,3之间的体心矢量：；。
15.进一步地，在步骤s5中，依据体心矢量以及密度梯度计算得到自适应加权权重系数，包括子步骤：依据体心矢量以及梯度，并设定如下自适应加权权
重系数关系：式中，考虑0单元密度梯度与四个重构模板的方向，增大与梯度方向平行的重构模板权重，减小与梯度方向垂直的重构模板权重。
16.进一步地，在步骤s5中，采用如下公式计算所述计算每个重构模板多项式的权重因子：；式中，表示一个防止分母为零的极小数，，表示步骤s5中计算出的自适应加权权重系数，表示未归一化的重构模板多项式的权重因子，表示归一化后的重构模板多项式的权重因子。
17.进一步地，在步骤s6中，采用如下公式计算0单元内密度变量的分布多项式：。
18.本发明的有益效果是：本发明方法通过耦合问题单元变量的梯度以及问题单元与其面相邻单元的位置关系，增大平行于变量梯度方向的模板的重构权重，减小垂直于变量梯度方向的模板的重构权重，能够有效抑制hweno限制器限制过程中残差的数值振荡，减小残差的数值大小，为hweno限制器的多项式重构过程提供一种新的方法。新的方法根据流场特征和网格关系自适应调整每个重构模板的激波强度在重构限制器中的重构权重，降低激波强度大且模板方向与密度梯度方向角度大的重构模板的重构权重，实现激波稳定捕捉，从而抑制激波带来的数值振荡，能够对关键的气动参数（如阻力、力矩和热流等）给出准确的计算结果，满足飞行器精细化设计需求，适用于包含强间断的流场以及在网格规模受计算机条件限制时的场景，给出旋涡、分离等复杂流动的高分辨率精细结构，提高计算的鲁棒性。
附图说明
19.为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
20.图1为网格单元及面相邻网格；图2为网格位置关系及密度梯度方向；图3为onera-m6机翼的压力云图；图4为hweno限制器有无自适应加权权重系数的密度残差随迭代步数变化曲线；图5为hweno限制器有无自适应加权权重系数的升力系数随迭代步数变化曲线。
具体实施方式
21.本说明书中所有实施例公开的所有特征，或隐含公开的所有方法或过程中的步骤，除了互相排斥的特征和/或步骤以外，均可以以任何方式组合和/或扩展、替换。
22.下面根据附图1~图5，对本发明的技术构思、工作原理、功效和工作过程作进一步详细说明。
23.本发明的设计思想是这样的：本发明的目的在于提供一种耦合流动及网格分布的重构限制器进行激波捕捉的方案。该方案中，首先给定线性加权权重系数，然后在每一时间迭代步后获得单元内变量分布多项式，同时计算单元变量梯度以及单元到单元的面相邻单元的体心方向；随后将单元变量梯度投影到体心方向并将其归一化；最后以此归一化的梯度投影乘以给定的线性加权权重系数获得候选模板变量分布多项式的自适应加权权重系数；该自适应加权权重系数耦合了流动特征和网格分布，避免了hweno限制器重构过程只考虑候选模板变量分布多项式的光滑程度，而不考虑当地的流动特征和网格分布，从而能够抑制时间迭代过程残差的数值振荡并降低残差的值。
24.在具体实施时，本发明实施例的技术方案包括：一种种耦合流动及网格分布的重构限制器进行激波捕捉的方法，用于激波捕捉中，旨在解决背景中提出的问题，具体实施时，可以包括以下步骤：s1，以密度作为重构限制器中定义自适应加权权重系数的物理量，在第n时间迭代步计算所有网格单元密度分布多项式，并探测出包含激波的问题单元，n为正整数；s2，定义重构限制器中重构模板构造方式；s3，按照步骤s2定义的重构模板构造方式，重构获得每个重构模板内的密度分布多项式，该多项式为重构限制器中的重构分量，用于重构单元最终的密度分布多项式实现激波捕捉；s4，计算步骤s3获得的每个重构模板的密度分布多项式的光滑因子，同时计算重构模板的密度分布多项式的密度梯度，该密度梯度将作为重构模板内是否包含激波的基本依据，并计算问题单元到面相邻单元之间的体心矢量；s5，将重构模板的密度梯度投影到该重构模板的体心矢量，实现当地流动特征与网格分布的耦合，以投影后的密度梯度计算自适应加权权重系数，然后依据计算的自适应加权权重系数以及步骤s4计算的光滑因子计算每个重构模板密度分布多项式的权重因子，该权重因子代表每个重构模板中激波强度；s6，根据步骤s3获得的每个重构模板的密度分布多项式以及步骤s5计算的每个重构模板密度分布多项式的权重因子，加权计算得到问题单元内密度变量最终的密度分布多项式，实现重构限制器对对单元密度变量的重构，实现激波捕捉。在实际应用时，包括如下步骤：
步骤一：如图1所示，0,1,2,3表示单元体心编号，a、b、c、d、e、f表示单元的节点编号。假如cfd数值模拟方法为二阶精度dg方法，在时间迭代步第n步时获得0单元及其面相邻单元1、2、3单元内的密度变量分布多项式为：（1）步骤二：如果0单元判定为“问题单元”，此时需要采用hweno限制器重构0单元内的守恒变量分布。按照hweno限制器的重构过程，首先构造4个候选模板，分别是0-0单元模板、0-1单元模板、0-2单元模板、0-3单元模板，将模板名定义为，以重构0单元内的密度变量为例，四个模板对应的密度分布多项式分别如下：（2）四个模板的密度分布多项式满足以下条件：（3）步骤三：计算光滑因子，光滑因子用来计算每个重构模板的密度分布多项式的权重，采用如下方法计算光滑因子：（4）该光滑因子度量模板中密度分布多项式的梯度量级，从公式（4）看到光滑因子只考虑了重构多项式梯度大小，未考虑梯度的方向。
25.步骤四：计算模板密度分布多项式的梯度以及0单元到面相邻单元1,
2,3之间的体心矢量，每个矢量与梯度都存在一定角度，如图2所示：（5）（6）步骤五：计算自适应加权权重系数，设定线性权，为体现网格分布在重构过程中的影响，计算新的：（7）考虑了0单元密度梯度与4个重构模板的方向，增大与梯度方向平行的重构模板权重，减小与梯度方向垂直的重构模板权重。
26.步骤六：计算重构多项式的权重因子：：（8）步骤七：依据计算的重构多项式权重因子，重构0单元内新的密度分布多项式：（9）至此完成0单元内密度的重构，基于重构后的0单元内密度可实现激波更加稳定捕捉。
27.选定的几何构型如图3所示，该构型为onera-m6机翼跨音速绕流，其有简单的外形及复杂的流场结构，常用来考核程序及方法的正确性。计算来流马赫数ma=0.8395，来流攻角3.06
°
，来流静压p=101325pa、来流静温t=280.0k。使用nnw-flowstar软件中uho非结构高阶模块计算流场，得到流场变量在网格单元上的分布，图3给出了采用基于自适应加权权重系数的hweno限制器获得的表面压力云图。
28.根据步骤一，首先在第n时间迭代步获得计算网格中所有网格的密度变量分布多项式。假设通过问题单元探测器探测到0单元为“问题单元”，按照步骤二定义的候选模板构造方式，重构每个候选模板的密度分布多项式，获得四个模板的密度候选多项式，。
29.根据计算的候选模板的密度分布多项式，计算每个多项式的光滑因子；同时计
算模板的密度分布多项式的梯度以及0单元到面相邻单元1,2,3之间的体心矢量。
30.依据计算的体心矢量以及梯度，计算得到自适应加权权重系数，最后依据计算的加权权重系数以及光滑因子计算每个重构模板多项式的权重因子。
31.根据步骤六计算得到的权重因子以及步骤二的重构模板多项式，获得0单元内密度变量新的分布多项式，从而完成hweno限制器对密度变量的限制过程，修改第n时间迭代步获得的0单元密度变量的自由度。图4给出了数值模拟onera-m6机翼跨音速绕流时hweno限制器有无自适应加权权重系数的密度残差随迭代步数变化曲线，从图4看到本发明实施例方法能够减少密度残差振荡，提高计算的鲁棒性，在原方法基础上进一步降低右端项残差的大小。图5给出了数值模拟onera-m6机翼跨音速绕流时hweno限制器有无自适应加权权重系数的升力随迭代步数变化曲线，可以看到本发明实施例方法获得的升力系数与原方法一致。
32.本发明功能如果以软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可读存储介质中。基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分或者该技术方案的部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中，在一台计算机设备（可以是个人计算机，服务器，或者网络设备等）以及相应的软件中执行本发明各个实施例所述方法的全部或部分步骤。而前述的存储介质包括：u盘、移动硬盘、或者光盘等各种可以存储程序代码的介质，进行测试或者实际的数据在程序实现中存在于只读存储器（random access memory，ram）、随机存取存储器（random access memory，ram）等。