一种适用于高超声速飞行器的高精度数值模拟方法
1.本发明涉及计算流体力学技术领域，具体为一种适用于高超声速飞行器的高精度数值模拟方法。
背景技术：
2.高超声速飞行器因其飞行速度快、突防能力强等特点，具有极其重要的军事和民用价值。例如，高超声速巡航武器，相比于常规武器，其杀伤力与命中精度更高，被拦截的可能性更低；高超声速战斗机攻击速度快，突防能力强，能够在短时间内高效地完成全球打击任务；高超声速防空导弹能够快速到达拦截空域，缩短拦截反应时间，提高拦截效率。因此，美国、俄罗斯、欧盟、日本等都投入了大量精力，致力于高超声速相关技术的研究，制定了一系列研究计划，如hyper-x、hyfly、hifire、hawc、hssw等。
3.从相关的研究可以发现，高超声速飞行器的气动环境非常复杂，其流场中存在激波、旋涡、粘性层、剪切层及其相互作用和干扰，既会增加飞行阻力，影响飞行器气动效率，又会产生强烈的气动加热，影响飞行安全。例如，1967年10月，x-15在进行马赫数为6.72的破纪录飞行时，发动机舱产生的激波与支架激波发生强烈干扰，导致飞行器多处严重受损。高超声速流动中的复杂相互作用与干扰使得高超声速飞行器的研究与发展面临着巨大的挑战。
4.目前，高超声速飞行器气动问题的研究主要通过数值模拟、风洞实验和飞行试验开展。而风洞实验和飞行试验成本巨大，实验测量技术也有限，无法作为复杂流动的常规研究手段。数值模拟因其便利性、灵活性和较低的成本成为了近年来研究复杂流动问题的重要手段。
5.然而传统数值模拟方法在高超声速飞行器的模拟中存在着模拟精度较低、鲁棒性不高的缺陷，仍有较大的改进空间。一方面，对于通量格式而言，传统的一维黎曼求解器在模拟多维复杂流动时会出现计算精度不足和稳定性降低的缺陷，balsara提出了一种基于“角点框架”模式的真正多维黎曼求解器，该求解器通过在网格单元角点处求解二维黎曼通量，简单、高效地实现了格式的多维效应。然而该格式只具备一阶精度，无法满足复杂流动的求解。另一方面，重构格式也是提高计算精度的重要方法。传统的muscl格式、tvd格式、weno格式等经过多年的发展，已广泛应用于各类流动求解器。但多维黎曼求解器具有四个间断初值，需要完成角点处四个状态变量的重构，传统适用于结构网格的高阶重构方法无法实现单元角点处的变量重构，从而限制了其向更高阶精度的推广。weno格式由于其高精度和基本无振荡特性，是重构过程的理想方案，但除balsara发展的ader-weno格式之外鲜有研究。因此研究多维高阶重构格式有利于进一步提高多维黎曼求解器的优势并提高对于复杂流动的计算精度。
技术实现要素：
6.为解决现有技术存在的问题，提高针对高超声速飞行器中激波与复杂流动干扰的
模拟精度，并将多维高阶重构方案应用于多维黎曼求解器，本发明提出一种适用于高超声速飞行器的，包含多维黎曼求解器、多维五阶weno格式的高精度数值模拟方法，为更加精准的高超声速数值模拟任务和高超声速飞行器设计工作提供技术支撑。
7.本发明首先针对传统cfd求解器计算精度低，稳定cfl数小的弊端，采用多维黎曼求解器计算界面通量。其次针对传统高阶重构格式难以应用于多维黎曼求解器的不足，发展了基于维度分裂的五阶weno格式以完成多维重构。最后，针对计算效率不够高、高阶格式仍存在振荡等问题，采用间断探测技术、保精度的限制器以提高重构效率与鲁棒性。
8.本发明的技术方案为：
9.所述一种适用于高超声速飞行器的高精度数值模拟方法，包括以下步骤：
10.步骤1：建立高超声速飞行器几何模型与计算域，并划分和读取高超声速飞行器计算网格，获得网格信息并根据来流条件赋予初值；
11.步骤2：对控制方程进行离散，得到有限体积形式的半离散格式；
12.步骤3：根据高超声速飞行器计算域中网格单元均值获得网格单元界面两侧的重构值；
13.步骤4：利用步骤3中所有网格单元界面两侧的重构值，获得高超声速流场中所有网格单元界面高斯点处的重构值用于多维黎曼求解器；
14.步骤5：利用步骤4得到的多维黎曼求解器求得界面通量；
15.步骤6：根据步骤5得到的界面通量确定残差，并进行时间推进求解，得到最终的高超声速飞行器流场。
16.进一步的，所述步骤2包括以下内容：
17.飞行器绕流微分形式控制方程如下：
[0018][0019]
其中，
[0020][0021]
t表示时间，x,y分别表示高超声速飞行器计算域中的横坐标和纵坐标，q是守恒变量，f和g表示x和y方向的通量，ρ,u,v,p,e分别表示高超声速流场密度、x方向速度，y方向速度，压强和能量；
[0022]
对通量项进行空间离散得到：
[0023][0024]
其中i，j是网格单元节点编号，δx,δy分别表示计算域中单个网格在x方向和y方向的宽度；f
i+1/2,j
和g
i,j+1/2
分别为x方向和y方向的界面通量，由高阶重构格式和多维黎曼求解器求解得到。
[0025]
进一步的，步骤3中，对于结点编号为i,j的网格单元v
i,j
，其单元均值为则利用单元均值沿x方向重构获得界面i+1/2左右两侧的均值和利用单元均值沿y方向重构获得界面j+1/2两侧的均值和
[0026]
进一步的，步骤3中，对于结点编号为i,j的网格单元v
i,j
，计算时，沿x方向建立三个子模板如下：
[0027]
s0＝{v
i-2,j
,v
i-1,j
,v
i,j
},s1＝{v
i-1,j
,v
i,j
,v
i+1,j
},s2＝{v
i,j
,v
i+1,j
,v
i+2,j
}
[0028]
将计算域的全局坐标转换为单元内的局部坐标，记插值点横坐标为xg＝x
i-1/2
+αδx，其中x
i-1/2,j
表示单元v
i,j
左侧界面的横坐标，δx表示单元v
i,j
在x方向的宽度；α在单元v
i,j
中的取值范围为α∈[0,1]，转换坐标后每个子模板的插值多项式记为对于取α＝1，各插值多项式如下：
[0029][0030][0031][0032]
计算每个多项式对应的光滑因子is0,is1,is2：
[0033][0034][0035][0036]
计算理想权重γ0,γ1,γ2，和无振荡权重
[0037][0038][0039][0040][0041][0042]
τ5＝|is
0-is2|
[0043]
其中，均表示中间计算变量，无特别物理意义，ε表示一个小量常数，用于避免分母为零，最后得到重构值
[0044][0045]
计算时，沿x方向建立以v
i+1,j
为中心的模板，并令α＝0；计算时，沿y方向建立以v
i,j
为中心的模板，并令α＝1；计算时，沿y方向建立以v
i,j+1
为中心的模板，并α＝0。
[0046]
进一步的，所述步骤4包括以下步骤：
[0047]
步骤4.1：利用多维五阶weno方法获得单元界面高斯点处初步重构值；
[0048]
步骤4.2：利用限制器限制步骤4.1得到的初步重构值，从而得到单元界面高斯点处重构值，并用于多维黎曼求解器。
[0049]
进一步的，步骤4.1中，对于网格界面i+1/2，界面左侧具有均值为的单元为界面右侧具有均值为的单元为对于五阶格式，界面具有三个高斯点g1、g2、g3，纵坐标分别为y
j2
＝y
j+1/2
，
[0050]
对于第二个高斯点g2，应用二维黎曼求解器计算通量，令对于第二个高斯点g2，应用二维黎曼求解器计算通量，令其中分别表示高斯点g2处左下角、左上角、右下角和右上角的初步重构值；
[0051]
对于第一个高斯点g1和第三个高斯点g3，分别计算四个区域的初步重构值
[0052]
对于高斯点g1左下角重构变量和高斯点g3处左上角的重构变量建立三个子模板如下：
[0053][0054][0055]
将计算域的全局坐标转换为单元内的局部坐标，记插值点坐标为yg＝y
j-1/2
+αδy，其中y
j-1/2
表示单元最下侧的纵坐标，δy表示单元在y方向的宽度；α在单元中的取值范围为α∈[0,1]，转换坐标后每个子模板的插值多项式记为对于对于对于对于子模板上建立的重构多项式如下：
[0056][0057]
[0058][0059]
计算每个子模板对应的光滑因子is0,is1,is2：
[0060][0061][0062][0063]
计算理想权重γ0,γ1,γ2：
[0064][0065][0066][0067]
根据理想权重的正负判断是否需要针对负权做特殊处理，如果min(γ0,γ1,γ2)＞0，则直接计算无振荡权重：
[0068][0069][0070]
τ5＝|is
0-is2|
[0071]
其中，均为中间计算变量，ε表示一个小量常数，用于避免分母为零；根据无振荡权重和插值多项式计算重构值和
[0072][0073]
如果min(γ0,γ1,γ2)＜0，则将理想权重分为正负两个部分
[0074][0075]
θ表示加权参数，据此计算新的正负理想权重
[0076][0077]
其中，σ
±
表示原理想权重正负两部分之和，然后计算正负两部分的无振荡权重
[0078][0079]
其中，为中间计算参数，利用正负两部分权重形成正负分离的重构多项式q
±
(α)：
[0080][0081]
最后重新组装为插值多项式计算g1和g3处的重构值：
[0082][0083]
对于高斯点g1处左上角的重构变量和高斯点g3处左下角的重构变量分别建立以和为中心的模板进行重构，α的取值分别为
[0084]
界面右侧所需的重构值利用单元右侧的均值进行重构。
[0085]
进一步的，步骤4.2中，对于该值位于网格单元v
i,j
内，首先以网格单元v
i,j
为中心划定3
×
3的网格单元范围，计算激波探测器指示因子detector：
[0086][0087][0088]
其中，表示9个单元内的流场变量的均值，σ表示9个单元流场变量的方差，用流场中的密度来计算激波探测器指示因子；根据detector的大小判断是否应用限制器：
[0089][0090][0091][0092][0093]
其中表示限制器指示因子，当时表示不使用限制器，当时表示激活限制器，mv和mv表示重构值的最大值和最小值，m和m表示9个单元内原始变量的最大值和最小值，最终经过限制后的重构值为：
[0094]
[0095]
进一步的，步骤5中，根据一维和二维黎曼求解器计算所有高斯点处的通量并组合，获得界面通量。
[0096]
进一步的，步骤5中，对于界面i+1/2，先求界面i+1/2第一个高斯点g1处的通量；作如下变量替换：
[0097][0098]
其中，q
ld
,q
lu
,q
rd
,q
ru
为替换后高斯点处左下角、左上角、右下角和右上角的重构变量，根据重构变量计算波速：
[0099][0100]
其中，sr和s
l
表示多维波传播模型中沿x方向的最大波速和最小波速，su和sd表示多维波传播模型中沿y方向传播的最大波速和最小波速；表示状态变量q
ru
处x方向的最大波速，其中u
ru
表示x方向速度，c
ru
表示声速；表示状态q
ru
处x方向的最小波速，处x方向的最小波速，表示(q
lu
,q
ru
)之间roe平均状态沿x方向的最大波速；表示x方向roe平均速度，表示roe平均声速，则表示roe平均声速，则表示(q
lu
,q
ru
)之间roe平均状态沿x方向的最小波速，
[0101]
多维波传播模型中，，二维黎曼通量表示为：
[0102][0103]
其中，f
ld
,f
lu
,f
rd
,f
ru
分别多维波传播模型表示左下角、左上角、右下角、右上角x方向的通量，g
ld
,g
lu
,g
rd
,g
ru
分别表示多维波传播模型左下角、左上角、右下角、右上角y方向的通量；
[0104]
对于第三个高斯点g3处的通量采用与高斯点g1相同的方法得到；
[0105]
对于第二个高斯点g2，利用一维黎曼求解器计算通量：
[0106]
[0107]
其中，f
l
,fr分别表示左右两侧的通量，分别表示高斯点处经过限制后左右两侧的重构值，上标“m”代表与界面中点处相关的物理量，和分别波传播模型中表示向左和向右传播的最大波速，其值定义如下：
[0108][0109]
式中a表示当地声速，上标“～”表示roe平均；
[0110]
最后，应用高斯积分得到界面通量f
i+1/2,j
：
[0111][0112]
进一步的，步骤6包括以下步骤：
[0113]
对时间变量采用三阶龙格-库塔离散公式将半离散有限体积格式转化为时空全离散有限体积格式：
[0114]q(1)
＝qn+δtl(qn)
[0115][0116][0117]
l表示残差，上标“n”表示时间步，利用时空全离散有限体积格式求解下一时间步上的流场变量值，依次推进，得到高超声速飞行器全流场稳定的数值模拟结果。
[0118]
有益效果
[0119]
本发明与现有技术相比，具有如下优点：
[0120]
1.本发明利用二维黎曼求解器计算界面通量，相对于一维黎曼求解器，在网格量相同时对高超声速飞行器流场中的激波间断和接触间断的分辨率更高。
[0121]
2.本发明利用多维重构策略将五阶weno重构方法应用于二维黎曼求解器，相对于目前的低阶重构方法，其捕捉激波、流动干扰等复杂流动结构的精度更高，能在网格量较少时仍具有较高精度，能够更精细地模拟高超声速飞行器的流场。
[0122]
本发明的附加方面和优点将在下面的描述中部分给出，部分将从下面的描述中变得明显，或通过本发明的实践了解到。
附图说明
[0123]
本发明的上述和/或附加的方面和优点从结合下面附图对实施例的描述中将变得明显和容易理解，其中：
[0124]
图1是本发明的实现流程图。
[0125]
图2是重构网格界面左右均值的示意图。
[0126]
图3是计算激波探测因子是选定的网格范围。
[0127]
图4是多维重构过程示意图。
[0128]
图5是多维黎曼求解器多维波传播模型的示意图。
[0129]
图6是一维黎曼求解器波传播模型示意图。
[0130]
图7是实施例一中，基于低阶格式和本发明方案计算的二维黎曼问题密度等值线图；(a)二阶重构方案，(b)五阶重构方案。
[0131]
图8是实施例二中，采用一维黎曼求解器和本发明方案求得的径向黎曼问题密度等值线图；(a)一维黎曼求解器和二阶重构方案，(b)二维黎曼求解器和五阶重构方案。
[0132]
图9是实施例三中，采用低阶重构方案和本发明方案求得的naca0012翼型超声速绕流压力云图；(a)二阶重构方案，(b)五阶重构方案。
[0133]
图10是实施例四中，采用低阶重构方案和本发明方案求得的高超声速双椭圆绕流马赫数云图；(a)二阶重构方案，(b)五阶重构方案。
[0134]
图11是实施例四中，采用低阶重构方案和本发明方案求得的高超声速双椭圆绕流激波探测器云图。
具体实施方式
[0135]
下面以高超声速飞行器流动模拟为例结合附图对本发明作进一步详细的说明，所述实施例是示例性的，旨在用于解释本发明，而不能理解为对本发明的限制。
[0136]
步骤1、建立高超声速飞行器几何模型与计算域，并划分和读取高超声速飞行器计算网格，获得网格信息并根据来流条件赋予初值；
[0137]
根据待分析的高超声速飞行器流动问题生成高超声速计算网格，读取网格信息，获得飞行器网格尺度、壁面节点坐标等。然后根据来流条件给计算域赋予一定的初值，例如网格单元v
i,j
的单元均值为
[0138]
步骤2、构造有限体积形式的半离散格式
[0139]
微分形式的控制方程如下：
[0140][0141]
其中，
[0142][0143]
其中，t表示时间，x,y分别表示横坐标和纵坐标，q是守恒变量，f和g表示x和y方向的通量，ρ,u,v,p,e分别表示高超声速飞行器流场的密度、x方向速度，y方向速度，压强和能量。
[0144]
对通量项进行空间离散可以得到：
[0145][0146]
其中i，j是单元节点编号；f
i+1/2,j
和g
i,j+1/2
分别为x方向和y方向的界面数值通量，由高阶重构格式和多维黎曼求解器求解得到，具体求解过程详见下述步骤。
[0147]
步骤3：根据高超声速飞行器计算域中网格单元均值获得网格单元界面两侧的重
构值。
[0148]
单元两侧界面均值以及多维重构均基于一般性的五阶weno方法，一般性的五阶weno方法如下。
[0149]
假设某一维均匀网格单元为ii＝(x
i-1/2
,x
i+1/2
)，网格间距为δx，网格单元上存在未知函数分布f(x)，已知f(x)在这些单元上的积分平均值为采用以单元ii为中心的5个相邻单元{i
i-2
,i
i-1
,ii,i
i+1
,i
i+2
}作为插值模板，重构出f(x)在单元v内任意插值点处的五阶近似，并将该重构多项式记为q(x)。
[0150]
在一般性的weno方法中，具备五阶精度的四次重构多项式q(x)应该表示为3个二次多项式的组合，以此来实现光滑区域的高精度和抑制间断附近的振荡。因此首先需要将上述五单元插值模板划分为如下三个子模板：
[0151]
s0＝{i
i-2
,i
i-1
,ii},s1＝{i
i-1
,ii,i
i+1
},s2＝{ii,i
i+1
,i
i+2
}
[0152]
然后在每个子模板上构造具备三阶精度的二次多项式pj(x),j＝0,1,2，每个多项式具备相近的表达式pj(x)＝ajx2+bjx+cj，三个插值多项式需要满足以下条件：
[0153][0154][0155][0156]
其中表示单元ij的均值，将每个多项式的表达形式带入上式即可求得所有系数，这里为了便于计算，将计算域的全局坐标转换为单元内的局部坐标，记插值点坐标为xg＝x
i-1/2
+αδx，其中x
i-1/2
表示单元ii最左侧的坐标，α在单元ii中的取值范围为α∈[0,1]，转换坐标后的插值多项式记为通过计算转换，求得的点xg处的五阶近似为：
[0157][0158][0159][0160]
同理，可以在光滑连续的五单元模板上构造具备五阶精度的四次插值多项式q(x)，其也需满足以下条件：
[0161][0162]
同样地将坐标转换为局部坐标后可以求得该多项式，记为
[0163]
理想情况下，三个二次多项式的加权组合应等于因此引入理想权重γ0,γ1,γ2，使得插值多项式在任意网格单元内满足：
[0164][0165]
对比等号两端的多项式系数可以求得三个理想权重系数：
[0166][0167][0168][0169]
在非理想情况下(间断附近)，包含间断的子模板的权重应该非常小以避免非物理振荡的产生。因此就需要一定的方法判断插值子模板的光滑程度并根据理想权重和光滑程度计算每个模板新的权重。新的权重可以避免加权后的插值多项式产生非物理振荡，因此可以称为无振荡权重。表示子模板光滑程度的参数称为光滑因子，三个子模板的光滑因子分别用is0,is1,is2表示，其计算方法为：
[0170][0171][0172][0173]
然后根据光滑因子和理想权重可以求得无振荡权重：
[0174][0175][0176]
τ5＝|is
0-is2|
[0177]
均为中间计算变量，无需赋予实际物理意义，ε表示一个小量常数，能够避免分母为零，可取ε＝1
×
10-16
。至此，任意情况下的插值多项式可以表示为：
[0178][0179]
下面阐述的具体实施过程，其余位置的重构变量用相同的方法计算。如图2所示，令均值为的单元为v
i,j
，利用单元均值沿x方向重构获得界面i+1/2左右两侧的均值和再利用单元均值沿y方向重构获得界面j+1/2两侧的均值和求时，沿x方向建立三个子模板如下：
[0180]
s0＝{v
i-2,j
,v
i-1,j
,v
i,j
},s1＝{v
i-1,j
,v
i,j
,v
i+1,j
},s2＝{v
i,j
,v
i+1,j
,v
i+2,j
}
[0181]
将计算域的全局坐标转换为单元内的局部坐标，记插值点坐标为xg＝x
i-1/2
+αδx，其中x
i-1/2,j
表示单元v
i,j
左侧界面的横坐标，δx表示单元v
i,j
在x方向的宽度。α在单元v
i,j
中的取值范围为α∈[0,1]，转换坐标后每个子模板的插值多项式记为对于取α＝1，三个子模板建立的插值多项式如下：
[0182][0183][0184][0185]
计算每个多项式对应的光滑因子is0,is1,is2：
[0186][0187][0188][0189]
计算理想权重和无振荡权重：
[0190][0191][0192][0193][0194][0195]
τ5＝|is
0-is2|
[0196]
最后计算重构值
[0197][0198]
对于沿x方向建立以v
i+1,j
为中心的模板，令α＝0，用同样的方法可以求得对于和沿y方向分别建立以v
i,j
和v
i,j+1
为中心的模板，分别令α＝1和0即可求解。这样可以求得所有网格交界面左右两侧的均值。
[0199]
步骤4、计算多维黎曼求解器所需重构值
[0200]
多维重构仍然基于一般性的五阶weno方法，但由于重构的位置为高斯点，其理想权重可能为负，因为在一维情况下，插值点的α一般取0或1，理想权重γ0,γ1,γ2均为正，但
当α取其他值时理想权重可能存在负值。而本发明需要进行多维重构，插值点处的α取值不为0或1，因此需要对其进行特殊处理。此外尽管基于weno方法，仍然需要一定的限制措施防止数值振荡。
[0201]
当理想权重为负时，可以采用简单的分离方法改进插值多项式。当min(γ0,γ1,γ2)＜0时，可以将理想权重分解为和正负两部分：
[0202][0203]
θ为加权参数，一般取3，然后构造出正负两部分的理想权重
[0204][0205]
其中，σ
±
为上一步分解的正负两部分之和，由此可以得到两个分离的多项式q
±
(α)：
[0206][0207]
这两个多项式满足以下条件：
[0208]
q(α)＝σ
+q+
(α)-σ-q-(α)
[0209]
同样的，需要考虑光滑因子和无振荡权重，分别利用正负两部分的理想权重计算：
[0210][0211]
即为正负两部分的无振荡权重，表示中间参数，ε表示一个小量常数，可取ε＝1
×
10-16
，能够避免分母为零。然后可以得到正负分裂的多项式q
±
(α)：
[0212][0213]
最终的weno重构公式为：
[0214]
q(α)＝σ
+q+
(α)-σ-q-(α)
[0215]
根据五阶weno方法，当重构i+1/2界面上所需要的变量时，可以先沿x方向应用五阶weno公式进行一维重构，得到界面左右两侧的均值和然后利用上一步重构得到界面两侧均值沿y方向再次进行一维重构即可获得多维黎曼求解器所需的4个输入，即
[0216]
理论上来说，上述策略可以完成重构并满足基本无振荡的需求，但在复杂流动中，间断附近依然存在数值振荡。因此本发明引入了一个限制器对重构值进行了二次限制以避免出现新的极值。
[0217]
假定间断附近的流场变量q满足q∈[m,m]，则高斯点处重构的变量也应满足q
gauss
∈[m,m]。求x方向的通量时，q
gauss
包括直接利用上述五阶weno方法所得的重构一般不能满足q
gauss
∈[m,m]，因此需要一定的方法进行限制。
[0218]
假定单元v
i,j
内具备五阶精度的插值多项式为p
i,j
(x,y)，限制之后的多项式表示为利用以下方法可以进行限制：
[0219][0220]
其中，
[0221][0222][0223]
上式中，m和m表示流场解的极值，定义为目标单元与邻居单元的最大值与最小值。mv和mv表示重构多项式极值，单元内多项式分布一般比较难获取，因此利用重构得到的值计算：
[0224][0225][0226]
限制器只需要在间断附近开启，在光滑区域并不需要，因此需要引入一个激波探测器确定计算域中的间断位置，间断探测器定义如下，
[0227][0228][0229]
上式和σ表示一定范围内流场变量的均值和标准差，一般情况下选择3
×
3的范围计算间断探测器，如图3所示，故n＝9，h表示网格尺度。因此结合激波探测器的限制器可以表示为：
[0230][0231][0232]
根据一般性的weno重构方法结合激波探测器和限制器，采用维度分裂的策略，即可形成完整的多维五阶重构策略。计算f
i+1/2,j
时，多维重构策略及通量求解方法如下：
[0233]
用二维网格的积分平均值在x方向进行一维重构，得到界面i+1/2处变量q在x方向上的高阶近似。基于偏左和偏右两个模板可重构出两个值，分别记作和令网格积分均值为的网格单元为v
i,j
，基于偏左的模板即以v
i,j
为中心的模板，基于偏右的模板即以v
i+1,j
为中心的模板。对于一维黎曼求解器而言，这就是其计算所需的输入状态量。
[0234]
分别用x方向上一维重构得到的网格积分平均值和沿y方向求出高斯
点y
jk
上的高阶近似值。这里假定界面左侧平均值为的单元为界面右侧平均值为的单元为对于五阶格式，共有3个高斯点，分别为y
j2
＝y
j+1/2
，
[0235]
其中第一个和第三个点为二维黎曼问题，第二个点实际上为一维黎曼问题，可以用一维黎曼求解器的公式，也可直接应用二维黎曼求解器。
[0236]
对于第一个和第三个点，和利用以为中心的模板重构，α的取值为为和分别利用以和为中心的模板进行重构，α的取值为
[0237]
对于第二个点，如果用一维黎曼求解器，则仅需和两个重构值，如果用二维黎曼求解器，可以令同理，界面i+1/2右侧所需的重构值则利用界面右侧的单元进行重构，α的取值相同。
[0238]
经过上文的重构，每个高斯点处已经得到了初步的重构值，记作经过上文的重构，每个高斯点处已经得到了初步的重构值，记作然后对于每一个值应用限制器以避免数值振荡。例如，对于限制后的值为限制后的值为的取值参照上文限制器计算方法。
[0239]
最后利用已经限制之后重构值计算界面通量，对于任意一个高斯点，利用上述二维黎曼求解器求解该点通量。
[0240]
对于界面i+1/2，详细实施过程如下，其余位置网格界面重构方法相同。
[0241]
如图4所示，假定交界面i+1/2左侧具有均值为的单元为界面右侧具有均值为的单元为对于五阶格式，界面具有三个高斯点g1、g2、g3，纵坐标分别为第二个高斯点g2考虑为一维黎曼问题或特殊的二维黎曼问题，若希望应用二维黎曼求解器计算通量，则令分别表示高斯点g2处左下角、左上角、右下角和右上角的初步重构值。然后，g1和g3为二维黎曼问题，四个区域的重构变量需分别计算。
[0242]
对于高斯点g1左下角重构变量和高斯点g3处左上角的重构变量建立三个子模板如下：
[0243][0244][0245]
将计算域的全局坐标转换为单元内的局部坐标，记插值点坐标为yg＝y
j-1/2
+αδy，
其中y
j-1/2
表示单元最下侧的纵坐标，δy表示单元在y方向的宽度。α在单元中的取值范围为α∈[0,1，转换坐标后每个子模板的插值多项式记为对于对于子模板上建立的重构多项式如下：
[0246][0247][0248][0249]
计算每个子模板对应的光滑因子is0,is1,is2：
[0250][0251][0252][0253]
计算理想权重γ0,γ1,γ2：
[0254][0255][0256][0257]
根据理想权重的正负判断是否需要针对负权做特殊处理，如果min(γ0,γ1,γ2)＞0，则直接计算无振荡权重：
[0258][0259][0260]
τ5＝|is
0-is2|
[0261]
其中，均为中间计算变量，ε表示一个小量常数，可取ε＝1
×
10-16
，能够避免分母为零，根据无振荡权重和子模板插值多项式直接计算重构变量：
[0262]
[0263][0264]
如果min(γ0,γ1,γ2)＜0，则需要将理想权重分为正负两个部分
[0265][0266]
θ表示加权参数，一般取3，据此计算新的正负理想权重
[0267][0268]
其中，σ
±
表示原理想权重正负两部分之和，然后计算正负两部分的无振荡权重光滑因子与前文相同：
[0269][0270]
其中，为中间计算参数，ε表示一个小量常数，可取ε＝1
×
10-16
，能够避免分母为零，利用两部分权重形成正负分离的重构多项式q
±
(α)：
[0271][0272]
最后重新组装为插值多项式计算g1和g3处的重构值：
[0273][0274][0275]
此外，对于g1处左上角的重构变量和g3处左下角的重构变量需要分别建立以和为中心的模板进行重构，α的取值为为中心的模板进行重构，α的取值为界面右侧所需的重构值利用单元右侧的均值进行重构。同理，单元界面j+1/2所有高斯点处的重构值利用界面j+1/2左右两侧的均值沿x方向进行重构。这样，能获得所有界面高斯点处的初步重构值。
[0276]
利用限制器限制初步重构变量：
[0277]
对于来说，该值位于单元v
i,j
内，首先以v
i,j
为中心划定3
×
3的网格单元范围，计算激波探测器的指示因子detector：
[0278][0279][0280]
其中，表示9个单元内的流场变量的均值，σ表示9个单元变量的方差，根据detector的大小判断是否应用限制器：
[0281][0282][0283][0284][0285]
其中表示限制器指示因子，当时表示不使用限制器，当时表示应用限制器，mv和mv表示重构值的最大值和最小值，m和m表示9个单元内原变量的最大值和最小值。最终经过限制后的重构值为：
[0286][0287]
据此方法对所有高斯点处的初步重构值进行判断和限制。
[0288]
步骤5、计算界面通量
[0289]
求解界面通量时为了充分考虑沿界面法向和横向传播的信息，提高计算精度，利用二维黎曼求解器计算界面通量。
[0290]
多维波传播模型示意图如图5所示，根据多维波传播模型，网格角点处含有q
ld
,q
lu
,q
rd
,q
ru
四个间断初值构成二维黎曼问题，据此可以求得x方向通量为：
[0291][0292]
上式中，f
ld
,f
lu
,f
rd
,f
ru
分别多维波传播模型表示左下角、左上角、右下角、右上角x方向的通量，g
ld
,g
lu
,g
rd
,g
ru
分别表示多维波传播模型左下角、左上角、右下角、右上角y方向的通量。sr和s
l
表示多维波传播模型中沿x方向的最大波速和最小波速，su和sd表示多维波传播模型中沿y方向传播的最大波速和最小波速。这四个极限波速的表达式如下：
[0293][0294]
其中，表示状态q
ru
处x方向的最大波速；表示状态q
ru
处x方向
的最小波速；表示(q
lu
,q
ru
)之间roe平均状态沿x方向的最大波速；表示(q
lu
,q
ru
)之间roe平均状态沿x方向的最小波速。
[0295]
同理，y方向的通量可以表示为：
[0296][0297]
常规的界面通量求解方法考虑界面中点和角点的黎曼问题，界面中点视为一维黎曼问题，角点视为二维黎曼问题，利用辛普森公式进行加权，界面i+1/2处的通量f
i+1/2,j
计算公式如下：
[0298][0299]
其中，表示角点通量，利用二维黎曼求解器计算，表示界面中点通量，利用一维黎曼求解器计算。一维黎曼求解器的波传播模型如图5所示，其计算方法如下：
[0300][0301]
其中，f
l
,fr表示左右两侧的通量，q
l
,qr表示左右两侧的流场变量。上标“m”代表与界面中点处相关的物理量，和分别表示向左向右传播的最大波速，其定义如下：
[0302][0303]
上式中，a表示当地声速，上标“～”表示roe平均。
[0304]
但辛普森积分仅具有三阶精度，无法匹配达到五阶精度的重构格式。因此本发明对界面通量计算方法进行了一定的改进，利用三点高斯积分达到五阶精度，其中f
i+1/2,j
和g
i,j+1/2
分别用下面的两个公式计算：
[0305][0306]
对于本发明的五阶格式，有β1＝4/9,β2＝5/18,β3＝4/9，对于界面i+1/2，三个高斯点分别为对于界面j+1/2，三个高斯点分别为其中第二个点
为界面中点，该点由于与两侧插值模板距离一样，可以视为一维黎曼问题，用一维黎曼求解器计算，但上式为了公式的简洁和形式一致，仍然用二维公式表示，因为当q
ld
＝q
lu
,q
rd
＝q
ru
时，该二维公式便退化为一维黎曼求解器。
[0307]
下面详细描述界面i+1/2处通量的详细计算过程，其余位置计算方法相同。
[0308]
先求i+1/2界面第一个高斯点g1处的通量，经过第4步多维重构，已经获得了高斯点处的重构变量，为了简化公式，作如下变量替换：
[0309][0310]
首先计算波速：
[0311][0312]
其中，表示状态q
ru
处x方向的最大波速，其中u
ru
表示x方向速度，c
ru
表示声速；表示状态q
ru
处x方向的最小波速，处x方向的最小波速，表示(q
lu
,q
ru
)之间roe平均状态沿x方向的最大波速。表示x方向roe平均速度，表示roe平均声速，则表示roe平均声速，则表示(q
lu
,q
ru
)之间roe平均状态沿x方向的最小波速，
[0313]
则二维黎曼通量为：
[0314][0315]
第三个高斯点方法同上，表示为
[0316]
对于第二个高斯点，利用一维黎曼求解器计算通量。
[0317][0318]
其中，上标“m”代表与界面中点处相关的物理量，和分别表示向左向右传播的最大波速，其值定义如下：
[0319]
[0320]
上式中，a表示当地声速，上标“～”表示roe平均。
[0321]
最后，应用高斯积分得到界面通量：
[0322][0323]
然后，对于沿y方向穿过界面j+1/2的通量，可以用同样的方法进行计算，但沿y方向的多维黎曼通量计算方法与x方向略有不同，同样根据图5的多维波传播模型，界面j+1/2上第一个高斯点通量的计算方法如下，y方向其余点利用相同的方法计算：
[0324][0325]
其中公式中各符号的意义与上文相同，但q
ld
,q
lu
,q
rd
,q
ru
需要利用界面y+1/2上第一个高斯点处的重构值。
[0326]
步骤6、时间推进求解
[0327]
对时间变量采用三阶龙格-库塔离散公式将半离散有限体积格式转化为时空全离散有限体积格式：
[0328]q(1)
＝qn+δtl(qn)
[0329][0330][0331]
l表示残差，上标“n”表示时间步。利用时空全离散有限体积格式求解下一时间步上的流场变量值，依次推进，得到高超声速飞行器全流场稳定的数值模拟结果。
[0332]
下面给出4个实施算例作为本发明公开方法的具体实施案例。
[0333]
实施例一、二维黎曼问题。
[0334]
该问题描述了2段激波和2段接触间断之间的相互作用，初始条件为：
[0335][0336]
计算网格为1000
×
1000。该算例虽然仅为数值算例，但包含较为复杂的流动结构，能够用于证明本发明的较高的数值精度。图7给出了t＝0.25时分别利用二阶重构格式和本发明所述五阶重构方案获得的密度分布，共30条轮廓线，范围为0.54-1.70。从图中可以看出本发明所述方案相对于低阶方案激波分辨率更高，捕捉的干扰区波系也更加精细。
[0337]
实施例二、径向黎曼问题
[0338]
该算例计算区域为[0,1]
×
[0,1]，cfl数为0.6，计算网格为200
×
200，初始条件为
[0339][0340]
图8给出了t＝0.13时刻分别采用一维黎曼求解器方案和本发明方案求得的密度等值线图。从图中可以看出在网格量相同时，本发明方案的二维黎曼求解器配合高阶重构方案后，对于激波间断和接触间断的捕捉精度更高。
[0341]
实施例三、naca0012超声速绕流问题
[0342]
该算例为翼型超声速绕流，在设计超声速和高超声速飞行器时需要进行翼型的计算。其中来流马赫数为2.0，迎角为10
°
，网格为300
×
300。图9为分别采用二阶重构方案和本发明方案所求的压力云图，从图中可以看出当网格量较少时，本发明方案的高阶格式对头部和尾部两道激波的捕捉精度更高。
[0343]
实施例四、高超声速双椭圆绕流
[0344]
该算例为高超声速双椭圆绕流，一定程度上能够代表高超声速飞行器的流动。来流马赫数为8.15，迎角为30
°
，来流密度为0.0231kg/m3，来流静压为370.7pa，来流静温为56k。图10为分别利用二阶重构方案和本发明的五阶重构方案求得的马赫数云图，从图中可以看出对于飞行器头部脱体激波的模拟二者精度相当，但在座舱前对于分离激波和二次激波的模拟，五阶格式的精度更高。从图11的激波探测器也可以看出五阶格式对于分离激波、二次激波以及附近的流动干扰求解更加精细。
[0345]
尽管上面已经示出和描述了本发明的实施例，可以理解的是，上述实施例是示例性的，不能理解为对本发明的限制，本领域的普通技术人员在不脱离本发明的原理和宗旨的情况下在本发明的范围内可以对上述实施例进行变化、修改、替换和变型。
技术特征：
1.一种适用于高超声速飞行器的高精度数值模拟方法，其特征在于：包括以下步骤：步骤1：建立高超声速飞行器几何模型与计算域，并划分和读取高超声速飞行器计算网格，获得网格信息并根据来流条件赋予初值；步骤2：对控制方程进行离散，得到有限体积形式的半离散格式；步骤3：根据高超声速飞行器计算域中网格单元均值获得网格单元界面两侧的重构值；步骤4：利用步骤3中所有网格单元界面两侧的重构值，获得高超声速流场中所有网格单元界面高斯点处的重构值用于多维黎曼求解器；步骤5：利用步骤4得到的多维黎曼求解器求得界面通量；步骤6：根据步骤5得到的界面通量确定残差，并进行时间推进求解，得到最终的高超声速飞行器流场。2.根据权利要求1所述一种适用于高超声速飞行器的高精度数值模拟方法，其特征在于：所述步骤2包括以下内容：飞行器绕流微分形式控制方程如下：其中，t表示时间，x,y分别表示高超声速飞行器计算域中的横坐标和纵坐标，q是守恒变量，f和g表示x和y方向的通量，ρ,u,v,p,e分别表示高超声速流场密度、x方向速度，y方向速度，压强和能量；对通量项进行空间离散得到：其中i，j是网格单元节点编号，δx,δy分别表示计算域中单个网格在x方向和y方向的宽度；f
i+1/2,j
和g
i,j+1/2
分别为x方向和y方向的界面通量，由高阶重构格式和多维黎曼求解器求解得到。3.根据权利要求1所述一种适用于高超声速飞行器的高精度数值模拟方法，其特征在于：步骤3中，对于结点编号为i,j的网格单元v
i,j
，其单元均值为则利用单元均值沿x方向重构获得界面i+1/2左右两侧的均值和利用单元均值沿y方向重构获得界面j+1/2两侧的均值和4.根据权利要求3所述一种适用于高超声速飞行器的高精度数值模拟方法，其特征在于：步骤3中，对于结点编号为i,j的网格单元v
i,j
，计算时，沿x方向建立三个子模板如下：
s0＝{v
i-2,j
,v
i-1,j
,v
i,j
},s1＝{v
i-1,j
,v
i,j
,v
i+1,j
},s2＝{v
i,j
,v
i+1,j
,v
i+2,j
}将计算域的全局坐标转换为单元内的局部坐标，记插值点横坐标为x
g
＝x
i-1/2
+αδx，其中x
i-1/2,j
表示单元v
i,j
左侧界面的横坐标，δx表示单元v
i,j
在x方向的宽度；α在单元v
i,j
中的取值范围为α∈[0,1]，转换坐标后每个子模板的插值多项式记为对于取α＝1，各插值多项式如下：取α＝1，各插值多项式如下：取α＝1，各插值多项式如下：计算每个多项式对应的光滑因子is0,is1,is2：：：计算理想权重γ0,γ1,γ2，和无振荡权重，和无振荡权重，和无振荡权重，和无振荡权重，和无振荡权重，和无振荡权重τ5＝|is
0-is2|其中，τ5均表示中间计算变量，无特别物理意义，ε表示一个小量常数，用于避免分母为零，最后得到重构值分母为零，最后得到重构值计算时，沿x方向建立以v
i+1,j
为中心的模板，并令α＝0；计算时，沿y方向建立以v
i,j
为中心的模板，并令α＝1；计算时，沿y方向建立以v
i,j+1
为中心的模板，并令α
＝0。5.根据权利要求1所述一种适用于高超声速飞行器的高精度数值模拟方法，其特征在于：所述步骤4包括以下步骤：步骤4.1：利用多维五阶weno方法获得单元界面高斯点处初步重构值；步骤4.2：利用限制器限制步骤4.1得到的初步重构值，从而得到单元界面高斯点处重构值，并用于多维黎曼求解器。6.根据权利要求5所述一种适用于高超声速飞行器的高精度数值模拟方法，其特征在于：步骤4.1中，对于网格界面i+1/2，界面左侧具有均值为的单元为界面右侧具有均值为的单元为对于五阶格式，界面具有三个高斯点g1、g2、g3，纵坐标分别为对于第二个高斯点g2，应用二维黎曼求解器计算通量，令对于第二个高斯点g2，应用二维黎曼求解器计算通量，令其中分别表示高斯点g2处左下角、左上角、右下角和右上角的初步重构值；对于第一个高斯点g1和第三个高斯点g3，分别计算四个区域的初步重构值：对于高斯点g1左下角重构变量和高斯点g3处左上角的重构变量建立三个子模板如下：子模板如下：将计算域的全局坐标转换为单元内的局部坐标，记插值点坐标为y
g
＝y
j-1/2
+αδy，其中y
j-1/2
表示单元最下侧的纵坐标，δy表示单元在y方向的宽度；α在单元中的取值范围为α∈[0,1]，转换坐标后每个子模板的插值多项式记为对于对于子模板上建立的重构多项式如下：下：下：计算每个子模板对应的光滑因子is0,is1,is2：
计算理想权重γ0,γ1,γ2：：：根据理想权重的正负判断是否需要针对负权做特殊处理，如果min(γ0,γ1,γ2)＞0，则直接计算无振荡权重：重：τ5＝|is
0-is2|其中，τ5均为中间计算变量，ε表示一个小量常数，用于避免分母为零；根据无振荡权重和插值多项式计算重构值和和如果min(γ0,γ1,γ2)＜0，则将理想权重分为正负两个部分)＜0，则将理想权重分为正负两个部分θ表示加权参数，据此计算新的正负理想权重θ表示加权参数，据此计算新的正负理想权重其中，σ
±
表示原理想权重正负两部分之和，然后计算正负两部分的无振荡权重表示原理想权重正负两部分之和，然后计算正负两部分的无振荡权重其中，为中间计算参数，利用正负两部分权重形成正负分离的重构多项式q
±
(α)：
最后重新组装为插值多项式计算g1和g3处的重构值：对于高斯点g1处左上角的重构变量和高斯点g3处左下角的重构变量分别建立以和为中心的模板进行重构，α的取值分别为界面右侧所需的重构值利用单元右侧的均值进行重构。7.根据权利要求5所述一种适用于高超声速飞行器的高精度数值模拟方法，其特征在于：步骤4.2中，对于该值位于网格单元v
i,j
内，首先以网格单元v
i,j
为中心划定3
×
3的网格单元范围，计算激波探测器指示因子detector：范围，计算激波探测器指示因子detector：其中，表示9个单元内的流场变量的均值，σ表示9个单元流场变量的方差，用流场中的密度来计算激波探测器指示因子；根据detector的大小判断是否应用限制器：密度来计算激波探测器指示因子；根据detector的大小判断是否应用限制器：密度来计算激波探测器指示因子；根据detector的大小判断是否应用限制器：密度来计算激波探测器指示因子；根据detector的大小判断是否应用限制器：其中表示限制器指示因子，当时表示不使用限制器，当时表示激活限制器，mv和mv表示重构值的最大值和最小值，m和m表示9个单元内原始变量的最大值和最小值，最终经过限制后的重构值为：8.根据权利要求1所述一种适用于高超声速飞行器的高精度数值模拟方法，其特征在于：步骤5中，根据一维和二维黎曼求解器计算所有高斯点处的通量并组合，获得界面通量。9.根据权利要求8所述一种适用于高超声速飞行器的高精度数值模拟方法，其特征在于：步骤5中，对于界面i+1/2，先求界面i+1/2第一个高斯点g1处的通量；作如下变量替换：其中，q
ld
,q
lu
,q
rd
,q
ru
为替换后高斯点处左下角、左上角、右下角和右上角的重构变量，
根据重构变量计算波速：其中，s
r
和s
l
表示多维波传播模型中沿x方向的最大波速和最小波速，s
u
和s
d
表示多维波传播模型中沿y方向传播的最大波速和最小波速；表示状态变量q
ru
处x方向的最大波速，其中u
ru
表示x方向速度，c
ru
表示声速；表示状态q
ru
处x方向的最小波速，处x方向的最小波速，表示(q
lu
,q
ru
)之间roe平均状态沿x方向的最大波速；表示x方向roe平均速度，表示roe平均声速，则表示roe平均声速，则表示(q
lu
,q
ru
)之间roe平均状态沿x方向的最小波速，多维波传播模型中，，二维黎曼通量表示为：其中，f
ld
,f
lu
,f
rd
,f
ru
分别多维波传播模型表示左下角、左上角、右下角、右上角x方向的通量，g
ld
,g
lu
,g
rd
,g
ru
分别表示多维波传播模型左下角、左上角、右下角、右上角y方向的通量；对于第三个高斯点g3处的通量采用与高斯点g1相同的方法得到；对于第二个高斯点g2，利用一维黎曼求解器计算通量：其中，f
l
,f
r
分别表示左右两侧的通量，分别表示高斯点处经过限制后左右两侧的重构值，上标“m”代表与界面中点处相关的物理量，和分别波传播模型中表示向左和向右传播的最大波速，其值定义如下：式中a表示当地声速，上标“～”表示roe平均；
最后，应用高斯积分得到界面通量f
i+1/2,j
：10.根据权利要求1所述一种适用于高超声速飞行器的高精度数值模拟方法，其特征在于：步骤6包括以下步骤：对时间变量采用三阶龙格-库塔离散公式将半离散有限体积格式转化为时空全离散有限体积格式：q
(1)
＝q
n
+δtl(q
n
))l表示残差，上标“n”表示时间步，利用时空全离散有限体积格式求解下一时间步上的流场变量值，依次推进，得到高超声速飞行器全流场稳定的数值模拟结果。
技术总结
本发明提出一种适用于高超声速飞行器的高精度数值模拟方法，首先建立高超声速飞行器几何模型与计算域，并划分和读取高超声速飞行器计算网格，获得网格信息并根据来流条件赋予初值；再对控制方程进行离散，得到有限体积形式的半离散格式；之后根据高超声速飞行器计算域中网格单元均值获得网格单元界面两侧的重构值；利用所有网格单元界面两侧的重构值，获得高超声速流场中所有网格单元界面高斯点处的重构值用于多维黎曼求解器；利用得到的多维黎曼求解器求得界面通量；根据界面通量确定残差，并进行时间推进求解，得到最终的高超声速飞行器流场。本发明能够为更加精准的高超声速数值模拟任务和高超声速飞行器设计工作提供技术支撑。技术支撑。技术支撑。
