CN115496006A - 一种适用于高超声速飞行器的高精度数值模拟方法 - Google Patents

一种适用于高超声速飞行器的高精度数值模拟方法 Download PDF

Info

Publication number
CN115496006A
CN115496006A CN202211089524.3A CN202211089524A CN115496006A CN 115496006 A CN115496006 A CN 115496006A CN 202211089524 A CN202211089524 A CN 202211089524A CN 115496006 A CN115496006 A CN 115496006A
Authority
CN
China
Prior art keywords
reconstruction
interface
grid
value
gaussian
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN202211089524.3A
Other languages
English (en)
Inventor
屈峰
刘青松
周伯霄
孙迪
白俊强
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Northwestern Polytechnical University
Original Assignee
Northwestern Polytechnical University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Northwestern Polytechnical University filed Critical Northwestern Polytechnical University
Priority to CN202211089524.3A priority Critical patent/CN115496006A/zh
Publication of CN115496006A publication Critical patent/CN115496006A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/28—Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/15—Vehicle, aircraft or watercraft design
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2113/00—Details relating to the application field
    • G06F2113/08—Fluids
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Geometry (AREA)
  • General Physics & Mathematics (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • General Engineering & Computer Science (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Pure & Applied Mathematics (AREA)
  • Aviation & Aerospace Engineering (AREA)
  • Computational Mathematics (AREA)
  • Automation & Control Theory (AREA)
  • Algebra (AREA)
  • Computing Systems (AREA)
  • Fluid Mechanics (AREA)
  • Mathematical Physics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明提出一种适用于高超声速飞行器的高精度数值模拟方法，首先建立高超声速飞行器几何模型与计算域，并划分和读取高超声速飞行器计算网格，获得网格信息并根据来流条件赋予初值；再对控制方程进行离散，得到有限体积形式的半离散格式；之后根据高超声速飞行器计算域中网格单元均值获得网格单元界面两侧的重构值；利用所有网格单元界面两侧的重构值，获得高超声速流场中所有网格单元界面高斯点处的重构值用于多维黎曼求解器；利用得到的多维黎曼求解器求得界面通量；根据界面通量确定残差，并进行时间推进求解，得到最终的高超声速飞行器流场。本发明能够为更加精准的高超声速数值模拟任务和高超声速飞行器设计工作提供技术支撑。

Description

一种适用于高超声速飞行器的高精度数值模拟方法
技术领域
本发明涉及计算流体力学技术领域，具体为一种适用于高超声速飞行器的高精度数值模拟方法。
背景技术
高超声速飞行器因其飞行速度快、突防能力强等特点，具有极其重要的军事和民用价值。例如，高超声速巡航武器，相比于常规武器，其杀伤力与命中精度更高，被拦截的可能性更低；高超声速战斗机攻击速度快，突防能力强，能够在短时间内高效地完成全球打击任务；高超声速防空导弹能够快速到达拦截空域，缩短拦截反应时间，提高拦截效率。因此，美国、俄罗斯、欧盟、日本等都投入了大量精力，致力于高超声速相关技术的研究，制定了一系列研究计划，如Hyper-X、HyFly、HIFiRE、HAWC、HSSW等。
从相关的研究可以发现，高超声速飞行器的气动环境非常复杂，其流场中存在激波、旋涡、粘性层、剪切层及其相互作用和干扰，既会增加飞行阻力，影响飞行器气动效率，又会产生强烈的气动加热，影响飞行安全。例如，1967年10月，X-15在进行马赫数为6.72的破纪录飞行时，发动机舱产生的激波与支架激波发生强烈干扰，导致飞行器多处严重受损。高超声速流动中的复杂相互作用与干扰使得高超声速飞行器的研究与发展面临着巨大的挑战。
目前，高超声速飞行器气动问题的研究主要通过数值模拟、风洞实验和飞行试验开展。而风洞实验和飞行试验成本巨大，实验测量技术也有限，无法作为复杂流动的常规研究手段。数值模拟因其便利性、灵活性和较低的成本成为了近年来研究复杂流动问题的重要手段。
然而传统数值模拟方法在高超声速飞行器的模拟中存在着模拟精度较低、鲁棒性不高的缺陷，仍有较大的改进空间。一方面，对于通量格式而言，传统的一维黎曼求解器在模拟多维复杂流动时会出现计算精度不足和稳定性降低的缺陷，Balsara提出了一种基于“角点框架”模式的真正多维黎曼求解器，该求解器通过在网格单元角点处求解二维黎曼通量，简单、高效地实现了格式的多维效应。然而该格式只具备一阶精度，无法满足复杂流动的求解。另一方面，重构格式也是提高计算精度的重要方法。传统的MUSCL格式、TVD格式、WENO格式等经过多年的发展，已广泛应用于各类流动求解器。但多维黎曼求解器具有四个间断初值，需要完成角点处四个状态变量的重构，传统适用于结构网格的高阶重构方法无法实现单元角点处的变量重构，从而限制了其向更高阶精度的推广。WENO格式由于其高精度和基本无振荡特性，是重构过程的理想方案，但除Balsara发展的ADER-WENO格式之外鲜有研究。因此研究多维高阶重构格式有利于进一步提高多维黎曼求解器的优势并提高对于复杂流动的计算精度。
发明内容
为解决现有技术存在的问题，提高针对高超声速飞行器中激波与复杂流动干扰的模拟精度，并将多维高阶重构方案应用于多维黎曼求解器，本发明提出一种适用于高超声速飞行器的，包含多维黎曼求解器、多维五阶WENO格式的高精度数值模拟方法，为更加精准的高超声速数值模拟任务和高超声速飞行器设计工作提供技术支撑。
本发明首先针对传统CFD求解器计算精度低，稳定CFL数小的弊端，采用多维黎曼求解器计算界面通量。其次针对传统高阶重构格式难以应用于多维黎曼求解器的不足，发展了基于维度分裂的五阶WENO格式以完成多维重构。最后，针对计算效率不够高、高阶格式仍存在振荡等问题，采用间断探测技术、保精度的限制器以提高重构效率与鲁棒性。
本发明的技术方案为：
所述一种适用于高超声速飞行器的高精度数值模拟方法，包括以下步骤：
步骤1：建立高超声速飞行器几何模型与计算域，并划分和读取高超声速飞行器计算网格，获得网格信息并根据来流条件赋予初值；
步骤2：对控制方程进行离散，得到有限体积形式的半离散格式；
步骤3：根据高超声速飞行器计算域中网格单元均值获得网格单元界面两侧的重构值；
步骤4：利用步骤3中所有网格单元界面两侧的重构值，获得高超声速流场中所有网格单元界面高斯点处的重构值用于多维黎曼求解器；
步骤5：利用步骤4得到的多维黎曼求解器求得界面通量；
步骤6：根据步骤5得到的界面通量确定残差，并进行时间推进求解，得到最终的高超声速飞行器流场。
进一步的，所述步骤2包括以下内容：
飞行器绕流微分形式控制方程如下：
Figure BDA0003836461890000021
其中，
Figure BDA0003836461890000022
t表示时间，x,y分别表示高超声速飞行器计算域中的横坐标和纵坐标，q是守恒变量，f和g表示x和y方向的通量，ρ,u,v,p,E分别表示高超声速流场密度、x方向速度，y方向速度，压强和能量；
对通量项进行空间离散得到：
Figure BDA0003836461890000031
其中i，j是网格单元节点编号，Δx,Δy分别表示计算域中单个网格在x方向和y方向的宽度；fi+1/2,j和gi,j+1/2分别为x方向和y方向的界面通量，由高阶重构格式和多维黎曼求解器求解得到。
进一步的，步骤3中，对于结点编号为i,j的网格单元Vi,j，其单元均值为
Figure BDA0003836461890000032
则利用单元均值沿x方向重构获得界面i+1/2左右两侧的均值
Figure BDA0003836461890000033
和
Figure BDA0003836461890000034
利用单元均值沿y方向重构获得界面j+1/2两侧的均值
Figure BDA0003836461890000035
和
Figure BDA0003836461890000036
进一步的，步骤3中，对于结点编号为i,j的网格单元Vi,j，计算
Figure BDA0003836461890000037
时，沿x方向建立三个子模板如下：
S0＝{Vi-2,j,Vi-1,j,Vi,j},S1＝{Vi-1,j,Vi,j,Vi+1,j},S2＝{Vi,j,Vi+1,j,Vi+2,j}
将计算域的全局坐标转换为单元内的局部坐标，记插值点横坐标为xG＝xi-1/2+αΔx，其中xi-1/2,j表示单元Vi,j左侧界面的横坐标，Δx表示单元Vi,j在x方向的宽度；α在单元Vi,j中的取值范围为α∈[0,1]，转换坐标后每个子模板的插值多项式记为
Figure BDA0003836461890000038
对于
Figure BDA0003836461890000039
取α＝1，各插值多项式如下：
Figure BDA00038364618900000310
Figure BDA00038364618900000311
Figure BDA00038364618900000312
计算每个多项式对应的光滑因子IS0,IS1,IS2：
Figure BDA00038364618900000313
Figure BDA00038364618900000314
Figure BDA00038364618900000315
计算理想权重γ0,γ1,γ2，和无振荡权重
Figure BDA00038364618900000316
Figure BDA0003836461890000041
Figure BDA0003836461890000042
Figure BDA0003836461890000043
Figure BDA0003836461890000044
Figure BDA0003836461890000045
τ5＝|IS0-IS2|
其中，
Figure BDA0003836461890000046
均表示中间计算变量，无特别物理意义，ε表示一个小量常数，用于避免分母为零，最后得到重构值
Figure BDA0003836461890000047
Figure BDA0003836461890000048
计算
Figure BDA0003836461890000049
时，沿x方向建立以Vi+1,j为中心的模板，并令α＝0；计算
Figure BDA00038364618900000410
时，沿y方向建立以Vi,j为中心的模板，并令α＝1；计算
Figure BDA00038364618900000411
时，沿y方向建立以Vi,j+1为中心的模板，并α＝0。
进一步的，所述步骤4包括以下步骤：
步骤4.1：利用多维五阶WENO方法获得单元界面高斯点处初步重构值；
步骤4.2：利用限制器限制步骤4.1得到的初步重构值，从而得到单元界面高斯点处重构值，并用于多维黎曼求解器。
进一步的，步骤4.1中，对于网格界面i+1/2，界面左侧具有均值为
Figure BDA00038364618900000412
的单元为
Figure BDA00038364618900000413
界面右侧具有均值为
Figure BDA00038364618900000414
的单元为
Figure BDA00038364618900000415
对于五阶格式，界面具有三个高斯点G1、G2、G3，纵坐标分别为
Figure BDA00038364618900000416
yj2＝yj+1/2，
Figure BDA00038364618900000417
对于第二个高斯点G2，应用二维黎曼求解器计算通量，令
Figure BDA00038364618900000418
Figure BDA00038364618900000419
其中
Figure BDA00038364618900000420
分别表示高斯点G2处左下角、左上角、右下角和右上角的初步重构值；
对于第一个高斯点G1和第三个高斯点G3，分别计算四个区域的初步重构值
对于高斯点G1左下角重构变量
Figure BDA00038364618900000421
和高斯点G3处左上角的重构变量
Figure BDA00038364618900000422
建立三个子模板如下：
Figure BDA0003836461890000051
Figure BDA0003836461890000052
将计算域的全局坐标转换为单元内的局部坐标，记插值点坐标为yG＝yj-1/2+αΔy，其中yj-1/2表示单元
Figure BDA0003836461890000053
最下侧的纵坐标，Δy表示单元
Figure BDA0003836461890000054
在y方向的宽度；α在单元
Figure BDA0003836461890000055
中的取值范围为α∈[0,1]，转换坐标后每个子模板的插值多项式记为
Figure BDA00038364618900000520
对于
Figure BDA0003836461890000057
Figure BDA0003836461890000058
对于
Figure BDA0003836461890000059
Figure BDA00038364618900000510
子模板上建立的重构多项式如下：
Figure BDA00038364618900000511
Figure BDA00038364618900000512
Figure BDA00038364618900000513
计算每个子模板对应的光滑因子IS0,IS1,IS2：
Figure BDA00038364618900000514
Figure BDA00038364618900000515
Figure BDA00038364618900000516
计算理想权重γ0,γ1,γ2：
Figure BDA00038364618900000517
Figure BDA00038364618900000518
Figure BDA00038364618900000519
根据理想权重的正负判断是否需要针对负权做特殊处理，如果min(γ0,γ1,γ2)＞0，则直接计算无振荡权重：
Figure BDA0003836461890000061
Figure BDA0003836461890000062
τ5＝|IS0-IS2|
其中，
Figure BDA0003836461890000063
均为中间计算变量，ε表示一个小量常数，用于避免分母为零；根据无振荡权重和插值多项式计算重构值
Figure BDA0003836461890000064
和
Figure BDA0003836461890000065
Figure BDA0003836461890000066
如果min(γ0,γ1,γ2)＜0，则将理想权重分为正负两个部分
Figure BDA0003836461890000067
Figure BDA0003836461890000068
θ表示加权参数，据此计算新的正负理想权重
Figure BDA0003836461890000069
Figure BDA00038364618900000610
其中，σ±表示原理想权重正负两部分之和，然后计算正负两部分的无振荡权重
Figure BDA00038364618900000611
Figure BDA00038364618900000612
其中，
Figure BDA00038364618900000613
为中间计算参数，利用正负两部分权重形成正负分离的重构多项式Q±(α)：
Figure BDA00038364618900000614
最后重新组装为插值多项式计算G1和G3处的重构值：
Figure BDA00038364618900000615
对于高斯点G1处左上角的重构变量
Figure BDA00038364618900000616
和高斯点G3处左下角的重构变量
Figure BDA00038364618900000617
分别建立以
Figure BDA00038364618900000618
和
Figure BDA00038364618900000619
为中心的模板进行重构，α的取值分别为
Figure BDA00038364618900000620
界面右侧所需的重构值利用单元右侧的均值进行重构。
进一步的，步骤4.2中，对于
Figure BDA00038364618900000621
该值位于网格单元Vi,j内，首先以网格单元Vi,j为中心划定3×3的网格单元范围，计算激波探测器指示因子Detector：
Figure BDA0003836461890000071
Figure BDA0003836461890000072
其中，
Figure BDA0003836461890000073
表示9个单元内的流场变量的均值，σ表示9个单元流场变量的方差，用流场中的密度来计算激波探测器指示因子；根据Detector的大小判断是否应用限制器：
Figure BDA0003836461890000074
Figure BDA0003836461890000075
Figure BDA0003836461890000076
Figure BDA0003836461890000077
其中
Figure BDA0003836461890000078
表示限制器指示因子，当
Figure BDA0003836461890000079
时表示不使用限制器，当
Figure BDA00038364618900000710
时表示激活限制器，MV和mV表示重构值的最大值和最小值，M和m表示9个单元内原始变量的最大值和最小值，最终经过限制后的重构值为：
Figure BDA00038364618900000711
进一步的，步骤5中，根据一维和二维黎曼求解器计算所有高斯点处的通量并组合，获得界面通量。
进一步的，步骤5中，对于界面i+1/2，先求界面i+1/2第一个高斯点G1处的通量；作如下变量替换：
Figure BDA00038364618900000712
其中，qLD,qLU,qRD,qRU为替换后高斯点处左下角、左上角、右下角和右上角的重构变量，根据重构变量计算波速：
Figure BDA00038364618900000713
其中，SR和SL表示多维波传播模型中沿x方向的最大波速和最小波速，SU和SD表示多维波传播模型中沿y方向传播的最大波速和最小波速；
Figure BDA0003836461890000081
表示状态变量qRU处x方向的最大波速，
Figure BDA0003836461890000082
其中uRU表示x方向速度，cRU表示声速；
Figure BDA0003836461890000083
表示状态qRU处x方向的最小波速，
Figure BDA0003836461890000084
Figure BDA0003836461890000085
表示(qLU,qRU)之间Roe平均状态沿x方向的最大波速；
Figure BDA0003836461890000086
表示x方向Roe平均速度，
Figure BDA0003836461890000087
表示Roe平均声速，则
Figure BDA0003836461890000088
Figure BDA0003836461890000089
表示(qLU,qRU)之间Roe平均状态沿x方向的最小波速，
Figure BDA00038364618900000810
多维波传播模型中，，二维黎曼通量表示为：
Figure BDA00038364618900000811
其中，fLD,fLU,fRD,fRU分别多维波传播模型表示左下角、左上角、右下角、右上角x方向的通量，gLD,gLU,gRD,gRU分别表示多维波传播模型左下角、左上角、右下角、右上角y方向的通量；
对于第三个高斯点G3处的通量
Figure BDA00038364618900000812
采用与高斯点G1相同的方法得到；
对于第二个高斯点G2，利用一维黎曼求解器计算通量：
Figure BDA00038364618900000813
其中，fL,fR分别表示左右两侧的通量，
Figure BDA00038364618900000814
分别表示高斯点处经过限制后左右两侧的重构值，上标“m”代表与界面中点处相关的物理量，
Figure BDA00038364618900000815
和
Figure BDA00038364618900000816
分别波传播模型中表示向左和向右传播的最大波速，其值定义如下：
Figure BDA00038364618900000817
式中a表示当地声速，上标“～”表示Roe平均；
最后，应用高斯积分得到界面通量fi+1/2,j：
Figure BDA0003836461890000091
进一步的，步骤6包括以下步骤：
对时间变量采用三阶龙格-库塔离散公式将半离散有限体积格式转化为时空全离散有限体积格式：
q(1)＝qn+ΔtL(qn)
Figure BDA0003836461890000092
Figure BDA0003836461890000093
L表示残差，上标“n”表示时间步，利用时空全离散有限体积格式求解下一时间步上的流场变量值，依次推进，得到高超声速飞行器全流场稳定的数值模拟结果。
有益效果
本发明与现有技术相比，具有如下优点：
1.本发明利用二维黎曼求解器计算界面通量，相对于一维黎曼求解器，在网格量相同时对高超声速飞行器流场中的激波间断和接触间断的分辨率更高。
2.本发明利用多维重构策略将五阶WENO重构方法应用于二维黎曼求解器，相对于目前的低阶重构方法，其捕捉激波、流动干扰等复杂流动结构的精度更高，能在网格量较少时仍具有较高精度，能够更精细地模拟高超声速飞行器的流场。
本发明的附加方面和优点将在下面的描述中部分给出，部分将从下面的描述中变得明显，或通过本发明的实践了解到。
附图说明
本发明的上述和/或附加的方面和优点从结合下面附图对实施例的描述中将变得明显和容易理解，其中：
图1是本发明的实现流程图。
图2是重构网格界面左右均值的示意图。
图3是计算激波探测因子是选定的网格范围。
图4是多维重构过程示意图。
图5是多维黎曼求解器多维波传播模型的示意图。
图6是一维黎曼求解器波传播模型示意图。
图7是实施例一中，基于低阶格式和本发明方案计算的二维黎曼问题密度等值线图；(a)二阶重构方案，(b)五阶重构方案。
图8是实施例二中，采用一维黎曼求解器和本发明方案求得的径向黎曼问题密度等值线图；(a)一维黎曼求解器和二阶重构方案，(b)二维黎曼求解器和五阶重构方案。
图9是实施例三中，采用低阶重构方案和本发明方案求得的NACA0012翼型超声速绕流压力云图；(a)二阶重构方案，(b)五阶重构方案。
图10是实施例四中，采用低阶重构方案和本发明方案求得的高超声速双椭圆绕流马赫数云图；(a)二阶重构方案，(b)五阶重构方案。
图11是实施例四中，采用低阶重构方案和本发明方案求得的高超声速双椭圆绕流激波探测器云图。
具体实施方式
下面以高超声速飞行器流动模拟为例结合附图对本发明作进一步详细的说明，所述实施例是示例性的，旨在用于解释本发明，而不能理解为对本发明的限制。
步骤1、建立高超声速飞行器几何模型与计算域，并划分和读取高超声速飞行器计算网格，获得网格信息并根据来流条件赋予初值；
根据待分析的高超声速飞行器流动问题生成高超声速计算网格，读取网格信息，获得飞行器网格尺度、壁面节点坐标等。然后根据来流条件给计算域赋予一定的初值，例如网格单元Vi,j的单元均值为
Figure BDA0003836461890000101
步骤2、构造有限体积形式的半离散格式
微分形式的控制方程如下：
Figure BDA0003836461890000102
其中，
Figure BDA0003836461890000103
其中，t表示时间，x,y分别表示横坐标和纵坐标，q是守恒变量，f和g表示x和y方向的通量，ρ,u,v,p,E分别表示高超声速飞行器流场的密度、x方向速度，y方向速度，压强和能量。
对通量项进行空间离散可以得到：
Figure BDA0003836461890000104
其中i，j是单元节点编号；fi+1/2,j和gi,j+1/2分别为x方向和y方向的界面数值通量，由高阶重构格式和多维黎曼求解器求解得到，具体求解过程详见下述步骤。
步骤3：根据高超声速飞行器计算域中网格单元均值获得网格单元界面两侧的重构值。
单元两侧界面均值以及多维重构均基于一般性的五阶WENO方法，一般性的五阶WENO方法如下。
假设某一维均匀网格单元为Ii＝(xi-1/2,xi+1/2)，网格间距为Δx，网格单元上存在未知函数分布f(x)，已知f(x)在这些单元上的积分平均值为
Figure BDA0003836461890000111
采用以单元Ii为中心的5个相邻单元{Ii-2,Ii-1,Ii,Ii+1,Ii+2}作为插值模板，重构出f(x)在单元V内任意插值点处的五阶近似，并将该重构多项式记为Q(x)。
在一般性的WENO方法中，具备五阶精度的四次重构多项式Q(x)应该表示为3个二次多项式的组合，以此来实现光滑区域的高精度和抑制间断附近的振荡。因此首先需要将上述五单元插值模板划分为如下三个子模板：
S0＝{Ii-2,Ii-1,Ii},S1＝{Ii-1,Ii,Ii+1},S2＝{Ii,Ii+1,Ii+2}
然后在每个子模板上构造具备三阶精度的二次多项式pj(x),j＝0,1,2，每个多项式具备相近的表达式pj(x)＝ajx2+bjx+cj，三个插值多项式需要满足以下条件：
Figure BDA0003836461890000112
Figure BDA0003836461890000113
Figure BDA0003836461890000114
其中
Figure BDA0003836461890000115
表示单元Ij的均值，将每个多项式的表达形式带入上式即可求得所有系数，这里为了便于计算，将计算域的全局坐标转换为单元内的局部坐标，记插值点坐标为xG＝xi-1/2+αΔx，其中xi-1/2表示单元Ii最左侧的坐标，α在单元Ii中的取值范围为α∈[0,1]，转换坐标后的插值多项式记为
Figure BDA0003836461890000116
通过计算转换，求得的点xG处的五阶近似为：
Figure BDA0003836461890000121
Figure BDA0003836461890000122
Figure BDA0003836461890000123
同理，可以在光滑连续的五单元模板上构造具备五阶精度的四次插值多项式Q(x)，其也需满足以下条件：
Figure BDA0003836461890000124
同样地将坐标转换为局部坐标后可以求得该多项式，记为
Figure BDA0003836461890000125
理想情况下，三个二次多项式
Figure BDA0003836461890000126
的加权组合应等于
Figure BDA0003836461890000127
因此引入理想权重γ0,γ1,γ2，使得插值多项式在任意网格单元内满足：
Figure BDA0003836461890000128
对比等号两端的多项式系数可以求得三个理想权重系数：
Figure BDA0003836461890000129
Figure BDA00038364618900001210
Figure BDA00038364618900001211
在非理想情况下(间断附近)，包含间断的子模板的权重应该非常小以避免非物理振荡的产生。因此就需要一定的方法判断插值子模板的光滑程度并根据理想权重和光滑程度计算每个模板新的权重。新的权重可以避免加权后的插值多项式产生非物理振荡，因此可以称为无振荡权重。表示子模板光滑程度的参数称为光滑因子，三个子模板的光滑因子分别用IS0,IS1,IS2表示，其计算方法为：
Figure BDA00038364618900001212
Figure BDA00038364618900001213
Figure BDA00038364618900001214
然后根据光滑因子和理想权重可以求得无振荡权重：
Figure BDA0003836461890000131
Figure BDA0003836461890000132
τ5＝|IS0-IS2|
Figure BDA0003836461890000133
均为中间计算变量，无需赋予实际物理意义，ε表示一个小量常数，能够避免分母为零，可取ε＝1×10-16。至此，任意情况下的插值多项式可以表示为：
Figure BDA0003836461890000134
下面阐述
Figure BDA0003836461890000135
的具体实施过程，其余位置的重构变量用相同的方法计算。如图2所示，令均值为
Figure BDA0003836461890000136
的单元为Vi,j，利用单元均值沿x方向重构获得界面i+1/2左右两侧的均值
Figure BDA0003836461890000137
和
Figure BDA0003836461890000138
再利用单元均值沿y方向重构获得界面j+1/2两侧的均值
Figure BDA0003836461890000139
和
Figure BDA00038364618900001310
求
Figure BDA00038364618900001311
时，沿x方向建立三个子模板如下：
S0＝{Vi-2,j,Vi-1,j,Vi,j},S1＝{Vi-1,j,Vi,j,Vi+1,j},S2＝{Vi,j,Vi+1,j,Vi+2,j}
将计算域的全局坐标转换为单元内的局部坐标，记插值点坐标为xG＝xi-1/2+αΔx，其中xi-1/2,j表示单元Vi,j左侧界面的横坐标，Δx表示单元Vi,j在x方向的宽度。α在单元Vi,j中的取值范围为α∈[0,1]，转换坐标后每个子模板的插值多项式记为
Figure BDA00038364618900001312
对于
Figure BDA00038364618900001313
取α＝1，三个子模板建立的插值多项式如下：
Figure BDA00038364618900001314
Figure BDA00038364618900001315
Figure BDA00038364618900001316
计算每个多项式对应的光滑因子IS0,IS1,IS2：
Figure BDA00038364618900001317
Figure BDA00038364618900001318
Figure BDA00038364618900001319
计算理想权重和无振荡权重：
Figure BDA0003836461890000141
Figure BDA0003836461890000142
Figure BDA0003836461890000143
Figure BDA0003836461890000144
Figure BDA0003836461890000145
τ5＝|IS0-IS2|
最后计算重构值
Figure BDA0003836461890000146
Figure BDA0003836461890000147
对于
Figure BDA0003836461890000148
沿x方向建立以Vi+1,j为中心的模板，令α＝0，用同样的方法可以求得
Figure BDA0003836461890000149
对于
Figure BDA00038364618900001410
和
Figure BDA00038364618900001411
沿y方向分别建立以Vi,j和Vi,j+1为中心的模板，分别令α＝1和0即可求解。这样可以求得所有网格交界面左右两侧的均值。
步骤4、计算多维黎曼求解器所需重构值
多维重构仍然基于一般性的五阶WENO方法，但由于重构的位置为高斯点，其理想权重可能为负，因为在一维情况下，插值点的α一般取0或1，理想权重γ0,γ1,γ2均为正，但当α取其他值时理想权重可能存在负值。而本发明需要进行多维重构，插值点处的α取值不为0或1，因此需要对其进行特殊处理。此外尽管基于WENO方法，仍然需要一定的限制措施防止数值振荡。
当理想权重为负时，可以采用简单的分离方法改进插值多项式。当min(γ0,γ1,γ2)＜0时，可以将理想权重分解为
Figure BDA00038364618900001412
和
Figure BDA00038364618900001413
正负两部分：
Figure BDA00038364618900001414
θ为加权参数，一般取3，然后构造出正负两部分的理想权重
Figure BDA00038364618900001415
Figure BDA00038364618900001416
其中，σ±为上一步分解的正负两部分之和，由此可以得到两个分离的多项式Q±(α)：
Figure BDA0003836461890000151
这两个多项式满足以下条件：
Q(α)＝σ+Q+(α)-σ-Q-(α)
同样的，需要考虑光滑因子和无振荡权重，分别利用正负两部分的理想权重计算：
Figure BDA0003836461890000152
Figure BDA0003836461890000153
即为正负两部分的无振荡权重，
Figure BDA0003836461890000154
表示中间参数，ε表示一个小量常数，可取ε＝1×10-16，能够避免分母为零。然后可以得到正负分裂的多项式Q±(α)：
Figure BDA0003836461890000155
最终的WENO重构公式为：
Q(α)＝σ+Q+(α)-σ-Q-(α)
根据五阶WENO方法，当重构i+1/2界面上所需要的变量时，可以先沿x方向应用五阶WENO公式进行一维重构，得到界面左右两侧的均值
Figure BDA0003836461890000156
和
Figure BDA0003836461890000157
然后利用上一步重构得到界面两侧均值沿y方向再次进行一维重构即可获得多维黎曼求解器所需的4个输入，即
Figure BDA0003836461890000158
理论上来说，上述策略可以完成重构并满足基本无振荡的需求，但在复杂流动中，间断附近依然存在数值振荡。因此本发明引入了一个限制器对重构值进行了二次限制以避免出现新的极值。
假定间断附近的流场变量q满足q∈[m,M]，则高斯点处重构的变量也应满足qGauss∈[m,M]。求x方向的通量时，qGauss包括
Figure BDA0003836461890000159
直接利用上述五阶WENO方法所得的重构一般不能满足qGauss∈[m,M]，因此需要一定的方法进行限制。
假定单元Vi,j内具备五阶精度的插值多项式为pi,j(x,y)，限制之后的多项式表示为
Figure BDA00038364618900001510
利用以下方法可以进行限制：
Figure BDA0003836461890000161
其中，
Figure BDA0003836461890000162
Figure BDA0003836461890000163
上式中，M和m表示流场解的极值，定义为目标单元与邻居单元的最大值与最小值。MV和mV表示重构多项式极值，单元内多项式分布一般比较难获取，因此利用重构得到的值计算：
Figure BDA0003836461890000164
Figure BDA0003836461890000165
限制器只需要在间断附近开启，在光滑区域并不需要，因此需要引入一个激波探测器确定计算域中的间断位置，间断探测器定义如下，
Figure BDA0003836461890000166
Figure BDA0003836461890000167
上式
Figure BDA0003836461890000168
和σ表示一定范围内流场变量的均值和标准差，一般情况下选择3×3的范围计算间断探测器，如图3所示，故n＝9，h表示网格尺度。因此结合激波探测器的限制器可以表示为：
Figure BDA0003836461890000169
Figure BDA00038364618900001610
根据一般性的WENO重构方法结合激波探测器和限制器，采用维度分裂的策略，即可形成完整的多维五阶重构策略。计算fi+1/2,j时，多维重构策略及通量求解方法如下：
用二维网格的积分平均值
Figure BDA00038364618900001611
在x方向进行一维重构，得到界面i+1/2处变量q在x方向上的高阶近似。基于偏左和偏右两个模板可重构出两个值，分别记作
Figure BDA00038364618900001612
和
Figure BDA0003836461890000171
令网格积分均值为
Figure BDA0003836461890000172
的网格单元为Vi,j，基于偏左的模板即以Vi,j为中心的模板，基于偏右的模板即以Vi+1,j为中心的模板。对于一维黎曼求解器而言，这就是其计算所需的输入状态量。
分别用x方向上一维重构得到的网格积分平均值
Figure BDA0003836461890000173
和
Figure BDA0003836461890000174
沿y方向求出高斯点yjk上的高阶近似值。这里假定界面左侧平均值为
Figure BDA0003836461890000175
的单元为
Figure BDA0003836461890000176
界面右侧平均值为
Figure BDA0003836461890000177
的单元为
Figure BDA0003836461890000178
对于五阶格式，共有3个高斯点，分别为
Figure BDA0003836461890000179
yj2＝yj+1/2，
Figure BDA00038364618900001710
其中第一个和第三个点为二维黎曼问题，第二个点实际上为一维黎曼问题，可以用一维黎曼求解器的公式，也可直接应用二维黎曼求解器。
对于第一个和第三个点，
Figure BDA00038364618900001711
和
Figure BDA00038364618900001712
利用以
Figure BDA00038364618900001713
为中心的模板重构，α的取值为
Figure BDA00038364618900001714
Figure BDA00038364618900001715
和
Figure BDA00038364618900001716
分别利用以
Figure BDA00038364618900001717
和
Figure BDA00038364618900001718
为中心的模板进行重构，α的取值为
Figure BDA00038364618900001719
对于第二个点，如果用一维黎曼求解器，则仅需
Figure BDA00038364618900001720
和
Figure BDA00038364618900001721
两个重构值，如果用二维黎曼求解器，可以令
Figure BDA00038364618900001722
同理，界面i+1/2右侧所需的重构值则利用界面右侧的单元进行重构，α的取值相同。
经过上文的重构，每个高斯点处已经得到了初步的重构值，记作
Figure BDA00038364618900001723
Figure BDA00038364618900001724
然后对于每一个值应用限制器以避免数值振荡。例如，对于
Figure BDA00038364618900001725
限制后的值为
Figure BDA00038364618900001726
Figure BDA00038364618900001727
的取值参照上文限制器计算方法。
最后利用已经限制之后重构值计算界面通量，对于任意一个高斯点，利用上述二维黎曼求解器
Figure BDA00038364618900001728
求解该点通量。
对于界面i+1/2，详细实施过程如下，其余位置网格界面重构方法相同。
如图4所示，假定交界面i+1/2左侧具有均值为
Figure BDA00038364618900001729
的单元为
Figure BDA00038364618900001730
界面右侧具有均值为
Figure BDA00038364618900001731
的单元为
Figure BDA00038364618900001732
对于五阶格式，界面具有三个高斯点G1、G2、G3，纵坐标分别为
Figure BDA00038364618900001734
第二个高斯点G2考虑为一维黎曼问题或特殊的二维黎曼问题，若希望应用二维黎曼求解器计算通量，则令
Figure BDA0003836461890000181
分别表示高斯点G2处左下角、左上角、右下角和右上角的初步重构值。然后，G1和G3为二维黎曼问题，四个区域的重构变量需分别计算。
对于高斯点G1左下角重构变量
Figure BDA0003836461890000182
和高斯点G3处左上角的重构变量
Figure BDA0003836461890000183
建立三个子模板如下：
Figure BDA0003836461890000184
Figure BDA0003836461890000185
将计算域的全局坐标转换为单元内的局部坐标，记插值点坐标为yG＝yj-1/2+αΔy，其中yj-1/2表示单元
Figure BDA0003836461890000186
最下侧的纵坐标，Δy表示单元
Figure BDA0003836461890000187
在y方向的宽度。α在单元
Figure BDA0003836461890000188
中的取值范围为α∈[0,1，转换坐标后每个子模板的插值多项式记为
Figure BDA0003836461890000189
对于
Figure BDA00038364618900001810
对于
Figure BDA00038364618900001811
子模板上建立的重构多项式如下：
Figure BDA00038364618900001812
Figure BDA00038364618900001813
Figure BDA00038364618900001814
计算每个子模板对应的光滑因子IS0,IS1,IS2：
Figure BDA00038364618900001815
Figure BDA00038364618900001816
Figure BDA00038364618900001817
计算理想权重γ0,γ1,γ2：
Figure BDA00038364618900001818
Figure BDA00038364618900001819
Figure BDA00038364618900001820
根据理想权重的正负判断是否需要针对负权做特殊处理，如果min(γ0,γ1,γ2)＞0，则直接计算无振荡权重：
Figure BDA0003836461890000191
Figure BDA0003836461890000192
τ5＝|IS0-IS2|
其中，
Figure BDA0003836461890000193
均为中间计算变量，ε表示一个小量常数，可取ε＝1×10-16，能够避免分母为零，根据无振荡权重和子模板插值多项式直接计算重构变量：
Figure BDA0003836461890000194
Figure BDA0003836461890000195
如果min(γ0,γ1,γ2)＜0，则需要将理想权重分为正负两个部分
Figure BDA0003836461890000196
Figure BDA0003836461890000197
θ表示加权参数，一般取3，据此计算新的正负理想权重
Figure BDA0003836461890000198
Figure BDA0003836461890000199
其中，σ±表示原理想权重正负两部分之和，然后计算正负两部分的无振荡权重
Figure BDA00038364618900001910
光滑因子与前文相同：
Figure BDA00038364618900001911
其中，
Figure BDA00038364618900001912
为中间计算参数，ε表示一个小量常数，可取ε＝1×10-16，能够避免分母为零，利用两部分权重形成正负分离的重构多项式Q±(α)：
Figure BDA00038364618900001913
最后重新组装为插值多项式计算G1和G3处的重构值：
Figure BDA00038364618900001914
Figure BDA00038364618900001915
此外，对于G1处左上角的重构变量
Figure BDA00038364618900001916
和G3处左下角的重构变量
Figure BDA00038364618900001917
需要分别建立以
Figure BDA0003836461890000201
和
Figure BDA0003836461890000202
为中心的模板进行重构，α的取值为
Figure BDA0003836461890000203
Figure BDA0003836461890000204
界面右侧所需的重构值利用单元右侧的均值进行重构。同理，单元界面j+1/2所有高斯点处的重构值利用界面j+1/2左右两侧的均值沿x方向进行重构。这样，能获得所有界面高斯点处的初步重构值。
利用限制器限制初步重构变量：
对于
Figure BDA0003836461890000205
来说，该值位于单元Vi,j内，首先以Vi,j为中心划定3×3的网格单元范围，计算激波探测器的指示因子Detector：
Figure BDA0003836461890000206
Figure BDA0003836461890000207
其中，
Figure BDA0003836461890000208
表示9个单元内的流场变量的均值，σ表示9个单元变量的方差，根据Detector的大小判断是否应用限制器：
Figure BDA0003836461890000209
Figure BDA00038364618900002010
Figure BDA00038364618900002011
Figure BDA00038364618900002012
其中
Figure BDA00038364618900002013
表示限制器指示因子，当
Figure BDA00038364618900002014
时表示不使用限制器，当
Figure BDA00038364618900002015
时表示应用限制器，MV和mV表示重构值的最大值和最小值，M和m表示9个单元内原变量的最大值和最小值。最终经过限制后的重构值为：
Figure BDA00038364618900002016
据此方法对所有高斯点处的初步重构值进行判断和限制。
步骤5、计算界面通量
求解界面通量时为了充分考虑沿界面法向和横向传播的信息，提高计算精度，利用二维黎曼求解器计算界面通量。
多维波传播模型示意图如图5所示，根据多维波传播模型，网格角点处含有qLD,qLU,qRD,qRU四个间断初值构成二维黎曼问题，据此可以求得x方向通量为：
Figure BDA0003836461890000211
上式中，fLD,fLU,fRD,fRU分别多维波传播模型表示左下角、左上角、右下角、右上角x方向的通量，gLD,gLU,gRD,gRU分别表示多维波传播模型左下角、左上角、右下角、右上角y方向的通量。SR和SL表示多维波传播模型中沿x方向的最大波速和最小波速，SU和SD表示多维波传播模型中沿y方向传播的最大波速和最小波速。这四个极限波速的表达式如下：
Figure BDA0003836461890000212
其中，
Figure BDA0003836461890000213
表示状态qRU处x方向的最大波速；
Figure BDA0003836461890000214
表示状态qRU处x方向的最小波速；
Figure BDA0003836461890000215
表示(qLU,qRU)之间Roe平均状态沿x方向的最大波速；
Figure BDA0003836461890000216
表示(qLU,qRU)之间Roe平均状态沿x方向的最小波速。
同理，y方向的通量可以表示为：
Figure BDA0003836461890000217
常规的界面通量求解方法考虑界面中点和角点的黎曼问题，界面中点视为一维黎曼问题，角点视为二维黎曼问题，利用辛普森公式进行加权，界面i+1/2处的通量fi+1/2,j计算公式如下：
Figure BDA0003836461890000218
其中，
Figure BDA0003836461890000221
表示角点通量，利用二维黎曼求解器计算，
Figure BDA0003836461890000222
表示界面中点通量，利用一维黎曼求解器计算。一维黎曼求解器的波传播模型如图5所示，其计算方法如下：
Figure BDA0003836461890000223
其中，fL,fR表示左右两侧的通量，qL,qR表示左右两侧的流场变量。上标“m”代表与界面中点处相关的物理量，
Figure BDA0003836461890000224
和
Figure BDA0003836461890000225
分别表示向左向右传播的最大波速，其定义如下：
Figure BDA0003836461890000226
上式中，a表示当地声速，上标“～”表示Roe平均。
但辛普森积分仅具有三阶精度，无法匹配达到五阶精度的重构格式。因此本发明对界面通量计算方法进行了一定的改进，利用三点高斯积分达到五阶精度，其中fi+1/2,j和gi,j+1/2分别用下面的两个公式计算：
Figure BDA0003836461890000227
对于本发明的五阶格式，有β1＝4/9,β2＝5/18,β3＝4/9，对于界面i+1/2，三个高斯点分别为
Figure BDA0003836461890000228
对于界面j+1/2，三个高斯点分别为
Figure BDA0003836461890000229
其中第二个点为界面中点，该点由于与两侧插值模板距离一样，可以视为一维黎曼问题，用一维黎曼求解器计算，但上式为了公式的简洁和形式一致，仍然用二维公式表示，因为当qLD＝qLU,qRD＝qRU时，该二维公式便退化为一维黎曼求解器。
下面详细描述界面i+1/2处通量的详细计算过程，其余位置计算方法相同。
先求i+1/2界面第一个高斯点G1处的通量，经过第4步多维重构，已经获得了高斯点处的重构变量，为了简化公式，作如下变量替换：
Figure BDA00038364618900002210
首先计算波速：
Figure BDA0003836461890000231
其中，
Figure BDA0003836461890000232
表示状态qRU处x方向的最大波速，
Figure BDA0003836461890000233
其中uRU表示x方向速度，cRU表示声速；
Figure BDA0003836461890000234
表示状态qRU处x方向的最小波速，
Figure BDA0003836461890000235
Figure BDA0003836461890000236
表示(qLU,qRU)之间Roe平均状态沿x方向的最大波速。
Figure BDA0003836461890000237
表示x方向Roe平均速度，
Figure BDA0003836461890000238
表示Roe平均声速，则
Figure BDA0003836461890000239
Figure BDA00038364618900002310
表示(qLU,qRU)之间Roe平均状态沿x方向的最小波速，
Figure BDA00038364618900002311
则二维黎曼通量为：
Figure BDA00038364618900002312
第三个高斯点方法同上，表示为
Figure BDA00038364618900002313
对于第二个高斯点，利用一维黎曼求解器计算通量。
Figure BDA00038364618900002314
其中，上标“m”代表与界面中点处相关的物理量，
Figure BDA00038364618900002315
和
Figure BDA00038364618900002316
分别表示向左向右传播的最大波速，其值定义如下：
Figure BDA00038364618900002317
上式中，a表示当地声速，上标“～”表示Roe平均。
最后，应用高斯积分得到界面通量：
Figure BDA0003836461890000241
然后，对于沿y方向穿过界面j+1/2的通量，可以用同样的方法进行计算，但沿y方向的多维黎曼通量计算方法与x方向略有不同，同样根据图5的多维波传播模型，界面j+1/2上第一个高斯点通量
Figure BDA0003836461890000242
的计算方法如下，y方向其余点利用相同的方法计算：
Figure BDA0003836461890000243
其中公式中各符号的意义与上文相同，但qLD,qLU,qRD,qRU需要利用界面y+1/2上第一个高斯点处的重构值。
步骤6、时间推进求解
对时间变量采用三阶龙格-库塔离散公式将半离散有限体积格式转化为时空全离散有限体积格式：
q(1)＝qn+ΔtL(qn)
Figure BDA0003836461890000244
Figure BDA0003836461890000245
L表示残差，上标“n”表示时间步。利用时空全离散有限体积格式求解下一时间步上的流场变量值，依次推进，得到高超声速飞行器全流场稳定的数值模拟结果。
下面给出4个实施算例作为本发明公开方法的具体实施案例。
实施例一、二维黎曼问题。
该问题描述了2段激波和2段接触间断之间的相互作用，初始条件为：
Figure BDA0003836461890000246
计算网格为1000×1000。该算例虽然仅为数值算例，但包含较为复杂的流动结构，能够用于证明本发明的较高的数值精度。图7给出了t＝0.25时分别利用二阶重构格式和本发明所述五阶重构方案获得的密度分布，共30条轮廓线，范围为0.54-1.70。从图中可以看出本发明所述方案相对于低阶方案激波分辨率更高，捕捉的干扰区波系也更加精细。
实施例二、径向黎曼问题
该算例计算区域为[0,1]×[0,1]，CFL数为0.6，计算网格为200×200，初始条件为
Figure BDA0003836461890000251
图8给出了t＝0.13时刻分别采用一维黎曼求解器方案和本发明方案求得的密度等值线图。从图中可以看出在网格量相同时，本发明方案的二维黎曼求解器配合高阶重构方案后，对于激波间断和接触间断的捕捉精度更高。
实施例三、NACA0012超声速绕流问题
该算例为翼型超声速绕流，在设计超声速和高超声速飞行器时需要进行翼型的计算。其中来流马赫数为2.0，迎角为10°，网格为300×300。图9为分别采用二阶重构方案和本发明方案所求的压力云图，从图中可以看出当网格量较少时，本发明方案的高阶格式对头部和尾部两道激波的捕捉精度更高。
实施例四、高超声速双椭圆绕流
该算例为高超声速双椭圆绕流，一定程度上能够代表高超声速飞行器的流动。来流马赫数为8.15，迎角为30°，来流密度为0.0231Kg/m3，来流静压为370.7Pa，来流静温为56K。图10为分别利用二阶重构方案和本发明的五阶重构方案求得的马赫数云图，从图中可以看出对于飞行器头部脱体激波的模拟二者精度相当，但在座舱前对于分离激波和二次激波的模拟，五阶格式的精度更高。从图11的激波探测器也可以看出五阶格式对于分离激波、二次激波以及附近的流动干扰求解更加精细。
尽管上面已经示出和描述了本发明的实施例，可以理解的是，上述实施例是示例性的，不能理解为对本发明的限制，本领域的普通技术人员在不脱离本发明的原理和宗旨的情况下在本发明的范围内可以对上述实施例进行变化、修改、替换和变型。

Claims (10)

1.一种适用于高超声速飞行器的高精度数值模拟方法，其特征在于：包括以下步骤：
步骤1：建立高超声速飞行器几何模型与计算域，并划分和读取高超声速飞行器计算网格，获得网格信息并根据来流条件赋予初值；
步骤2：对控制方程进行离散，得到有限体积形式的半离散格式；
步骤3：根据高超声速飞行器计算域中网格单元均值获得网格单元界面两侧的重构值；
步骤4：利用步骤3中所有网格单元界面两侧的重构值，获得高超声速流场中所有网格单元界面高斯点处的重构值用于多维黎曼求解器；
步骤5：利用步骤4得到的多维黎曼求解器求得界面通量；
步骤6：根据步骤5得到的界面通量确定残差，并进行时间推进求解，得到最终的高超声速飞行器流场。
2.根据权利要求1所述一种适用于高超声速飞行器的高精度数值模拟方法，其特征在于：所述步骤2包括以下内容：
飞行器绕流微分形式控制方程如下：
Figure FDA0003836461880000011
其中，
Figure FDA0003836461880000012
t表示时间，x,y分别表示高超声速飞行器计算域中的横坐标和纵坐标，q是守恒变量，f和g表示x和y方向的通量，ρ,u,v,p,E分别表示高超声速流场密度、x方向速度，y方向速度，压强和能量；
对通量项进行空间离散得到：
Figure FDA0003836461880000013
其中i，j是网格单元节点编号，Δx,Δy分别表示计算域中单个网格在x方向和y方向的宽度；fi+1/2,j和gi,j+1/2分别为x方向和y方向的界面通量，由高阶重构格式和多维黎曼求解器求解得到。
3.根据权利要求1所述一种适用于高超声速飞行器的高精度数值模拟方法，其特征在于：步骤3中，对于结点编号为i,j的网格单元Vi,j，其单元均值为
Figure FDA0003836461880000021
则利用单元均值沿x方向重构获得界面i+1/2左右两侧的均值
Figure FDA0003836461880000022
和
Figure FDA0003836461880000023
利用单元均值沿y方向重构获得界面j+1/2两侧的均值
Figure FDA0003836461880000024
和
Figure FDA0003836461880000025
4.根据权利要求3所述一种适用于高超声速飞行器的高精度数值模拟方法，其特征在于：步骤3中，对于结点编号为i,j的网格单元Vi,j，计算
Figure FDA0003836461880000026
时，沿x方向建立三个子模板如下：
S0＝{Vi-2,j,Vi-1,j,Vi,j},S1＝{Vi-1,j,Vi,j,Vi+1,j},S2＝{Vi,j,Vi+1,j,Vi+2,j}
将计算域的全局坐标转换为单元内的局部坐标，记插值点横坐标为xG＝xi-1/2+αΔx，其中xi-1/2,j表示单元Vi,j左侧界面的横坐标，Δx表示单元Vi,j在x方向的宽度；α在单元Vi,j中的取值范围为α∈[0,1]，转换坐标后每个子模板的插值多项式记为
Figure FDA0003836461880000027
对于
Figure FDA0003836461880000028
取α＝1，各插值多项式如下：
Figure FDA0003836461880000029
Figure FDA00038364618800000210
Figure FDA00038364618800000211
计算每个多项式对应的光滑因子IS0,IS1,IS2：
Figure FDA00038364618800000212
Figure FDA00038364618800000213
Figure FDA00038364618800000214
计算理想权重γ0,γ1,γ2，和无振荡权重
Figure FDA00038364618800000215
Figure FDA0003836461880000031
Figure FDA0003836461880000032
Figure FDA0003836461880000033
Figure FDA0003836461880000034
Figure FDA0003836461880000035
τ5＝|IS0-IS2|
其中，
Figure FDA0003836461880000036
τ5均表示中间计算变量，无特别物理意义，ε表示一个小量常数，用于避免分母为零，最后得到重构值
Figure FDA0003836461880000037
Figure FDA0003836461880000038
计算
Figure FDA0003836461880000039
时，沿x方向建立以Vi+1,j为中心的模板，并令α＝0；计算
Figure FDA00038364618800000310
时，沿y方向建立以Vi,j为中心的模板，并令α＝1；计算
Figure FDA00038364618800000311
时，沿y方向建立以Vi,j+1为中心的模板，并令α＝0。
5.根据权利要求1所述一种适用于高超声速飞行器的高精度数值模拟方法，其特征在于：所述步骤4包括以下步骤：
步骤4.1：利用多维五阶WENO方法获得单元界面高斯点处初步重构值；
步骤4.2：利用限制器限制步骤4.1得到的初步重构值，从而得到单元界面高斯点处重构值，并用于多维黎曼求解器。
6.根据权利要求5所述一种适用于高超声速飞行器的高精度数值模拟方法，其特征在于：步骤4.1中，对于网格界面i+1/2，界面左侧具有均值为
Figure FDA00038364618800000312
的单元为
Figure FDA00038364618800000313
界面右侧具有均值为
Figure FDA00038364618800000314
的单元为
Figure FDA00038364618800000315
对于五阶格式，界面具有三个高斯点G1、G2、G3，纵坐标分别为
Figure FDA00038364618800000316
对于第二个高斯点G2，应用二维黎曼求解器计算通量，令
Figure FDA00038364618800000317
Figure FDA00038364618800000318
其中
Figure FDA00038364618800000319
分别表示高斯点G2处左下角、左上角、右下角和右上角的初步重构值；
对于第一个高斯点G1和第三个高斯点G3，分别计算四个区域的初步重构值：
对于高斯点G1左下角重构变量
Figure FDA0003836461880000041
和高斯点G3处左上角的重构变量
Figure FDA0003836461880000042
建立三个子模板如下：
Figure FDA0003836461880000043
Figure FDA0003836461880000044
将计算域的全局坐标转换为单元内的局部坐标，记插值点坐标为yG＝yj-1/2+αΔy，其中yj-1/2表示单元
Figure FDA0003836461880000045
最下侧的纵坐标，Δy表示单元
Figure FDA0003836461880000046
在y方向的宽度；α在单元
Figure FDA0003836461880000047
中的取值范围为α∈[0,1]，转换坐标后每个子模板的插值多项式记为
Figure FDA0003836461880000048
对于
Figure FDA0003836461880000049
对于
Figure FDA00038364618800000410
子模板上建立的重构多项式如下：
Figure FDA00038364618800000411
Figure FDA00038364618800000412
Figure FDA00038364618800000413
计算每个子模板对应的光滑因子IS0,IS1,IS2：
Figure FDA00038364618800000414
Figure FDA00038364618800000415
Figure FDA00038364618800000416
计算理想权重γ0,γ1,γ2：
Figure FDA00038364618800000417
Figure FDA00038364618800000418
Figure FDA00038364618800000419
根据理想权重的正负判断是否需要针对负权做特殊处理，如果min(γ0,γ1,γ2)＞0，则直接计算无振荡权重：
Figure FDA0003836461880000051
Figure FDA0003836461880000052
τ5＝|IS0-IS2|
其中，
Figure FDA0003836461880000053
τ5均为中间计算变量，ε表示一个小量常数，用于避免分母为零；根据无振荡权重和插值多项式计算重构值
Figure FDA0003836461880000054
和
Figure FDA0003836461880000055
Figure FDA0003836461880000056
如果min(γ0,γ1,γ2)＜0，则将理想权重分为正负两个部分
Figure FDA0003836461880000057
Figure FDA0003836461880000058
θ表示加权参数，据此计算新的正负理想权重
Figure FDA0003836461880000059
Figure FDA00038364618800000510
其中，σ±表示原理想权重正负两部分之和，然后计算正负两部分的无振荡权重
Figure FDA00038364618800000511
Figure FDA00038364618800000512
其中，
Figure FDA00038364618800000513
为中间计算参数，利用正负两部分权重形成正负分离的重构多项式Q±(α)：
Figure FDA00038364618800000514
最后重新组装为插值多项式计算G1和G3处的重构值：
Figure FDA00038364618800000515
对于高斯点G1处左上角的重构变量
Figure FDA00038364618800000516
和高斯点G3处左下角的重构变量
Figure FDA00038364618800000517
分别建立以
Figure FDA00038364618800000518
和
Figure FDA00038364618800000519
为中心的模板进行重构，α的取值分别为
Figure FDA00038364618800000520
界面右侧所需的重构值利用单元右侧的均值进行重构。
7.根据权利要求5所述一种适用于高超声速飞行器的高精度数值模拟方法，其特征在于：步骤4.2中，对于
Figure FDA0003836461880000061
该值位于网格单元Vi,j内，首先以网格单元Vi,j为中心划定3×3的网格单元范围，计算激波探测器指示因子Detector：
Figure FDA0003836461880000062
Figure FDA0003836461880000063
其中，
Figure FDA0003836461880000064
表示9个单元内的流场变量的均值，σ表示9个单元流场变量的方差，用流场中的密度来计算激波探测器指示因子；根据Detector的大小判断是否应用限制器：
Figure FDA0003836461880000065
Figure FDA0003836461880000066
Figure FDA0003836461880000067
Figure FDA0003836461880000068
其中
Figure FDA0003836461880000069
表示限制器指示因子，当
Figure FDA00038364618800000610
时表示不使用限制器，当
Figure FDA00038364618800000611
时表示激活限制器，MV和mV表示重构值的最大值和最小值，M和m表示9个单元内原始变量的最大值和最小值，最终经过限制后的重构值为：
Figure FDA00038364618800000612
8.根据权利要求1所述一种适用于高超声速飞行器的高精度数值模拟方法，其特征在于：步骤5中，根据一维和二维黎曼求解器计算所有高斯点处的通量并组合，获得界面通量。
9.根据权利要求8所述一种适用于高超声速飞行器的高精度数值模拟方法，其特征在于：步骤5中，对于界面i+1/2，先求界面i+1/2第一个高斯点G1处的通量；作如下变量替换：
Figure FDA00038364618800000613
其中，qLD,qLU,qRD,qRU为替换后高斯点处左下角、左上角、右下角和右上角的重构变量，根据重构变量计算波速：
Figure FDA0003836461880000071
其中，SR和SL表示多维波传播模型中沿x方向的最大波速和最小波速，SU和SD表示多维波传播模型中沿y方向传播的最大波速和最小波速；
Figure FDA0003836461880000072
表示状态变量qRU处x方向的最大波速，
Figure FDA0003836461880000073
其中uRU表示x方向速度，cRU表示声速；
Figure FDA0003836461880000074
表示状态qRU处x方向的最小波速，
Figure FDA0003836461880000075
Figure FDA0003836461880000076
表示(qLU,qRU)之间Roe平均状态沿x方向的最大波速；
Figure FDA0003836461880000077
表示x方向Roe平均速度，
Figure FDA0003836461880000078
表示Roe平均声速，则
Figure FDA0003836461880000079
Figure FDA00038364618800000710
表示(qLU,qRU)之间Roe平均状态沿x方向的最小波速，
Figure FDA00038364618800000711
多维波传播模型中，，二维黎曼通量表示为：
Figure FDA00038364618800000712
其中，fLD,fLU,fRD,fRU分别多维波传播模型表示左下角、左上角、右下角、右上角x方向的通量，gLD,gLU,gRD,gRU分别表示多维波传播模型左下角、左上角、右下角、右上角y方向的通量；
对于第三个高斯点G3处的通量
Figure FDA00038364618800000713
采用与高斯点G1相同的方法得到；
对于第二个高斯点G2，利用一维黎曼求解器计算通量：
Figure FDA00038364618800000714
其中，fL,fR分别表示左右两侧的通量，
Figure FDA00038364618800000715
分别表示高斯点处经过限制后左右两侧的重构值，上标“m”代表与界面中点处相关的物理量，
Figure FDA00038364618800000716
和
Figure FDA00038364618800000717
分别波传播模型中表示向左和向右传播的最大波速，其值定义如下：
Figure FDA0003836461880000081
式中a表示当地声速，上标“～”表示Roe平均；
最后，应用高斯积分得到界面通量fi+1/2,j：
Figure FDA0003836461880000082
10.根据权利要求1所述一种适用于高超声速飞行器的高精度数值模拟方法，其特征在于：步骤6包括以下步骤：
对时间变量采用三阶龙格-库塔离散公式将半离散有限体积格式转化为时空全离散有限体积格式：
q(1)＝qn+ΔtL(qn)
Figure FDA0003836461880000083
Figure FDA0003836461880000084
L表示残差，上标“n”表示时间步，利用时空全离散有限体积格式求解下一时间步上的流场变量值，依次推进，得到高超声速飞行器全流场稳定的数值模拟结果。
CN202211089524.3A 2022-09-07 2022-09-07 一种适用于高超声速飞行器的高精度数值模拟方法 Pending CN115496006A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202211089524.3A CN115496006A (zh) 2022-09-07 2022-09-07 一种适用于高超声速飞行器的高精度数值模拟方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202211089524.3A CN115496006A (zh) 2022-09-07 2022-09-07 一种适用于高超声速飞行器的高精度数值模拟方法

Publications (1)

Publication Number Publication Date
CN115496006A true CN115496006A (zh) 2022-12-20

Family

ID=84468680

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202211089524.3A Pending CN115496006A (zh) 2022-09-07 2022-09-07 一种适用于高超声速飞行器的高精度数值模拟方法

Country Status (1)

Country Link
CN (1) CN115496006A (zh)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116070554A (zh) * 2023-04-06 2023-05-05 中国人民解放军国防科技大学 高超声速飞行器气动热载荷计算方法、装置及设备
CN116227388A (zh) * 2023-04-21 2023-06-06 中国空气动力研究与发展中心计算空气动力研究所 高超流动模拟cfl数动态调整方法、系统、设备及介质
CN116341421A (zh) * 2023-05-22 2023-06-27 中国空气动力研究与发展中心计算空气动力研究所 高超声速流场数值模拟方法、系统、电子设备及存储介质
CN116522827A (zh) * 2023-07-04 2023-08-01 北京凌云智擎软件有限公司 一种网格单元边界面的流动变量重构方法、设备及装置

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116070554A (zh) * 2023-04-06 2023-05-05 中国人民解放军国防科技大学 高超声速飞行器气动热载荷计算方法、装置及设备
CN116227388A (zh) * 2023-04-21 2023-06-06 中国空气动力研究与发展中心计算空气动力研究所 高超流动模拟cfl数动态调整方法、系统、设备及介质
CN116227388B (zh) * 2023-04-21 2023-07-07 中国空气动力研究与发展中心计算空气动力研究所 高超流动模拟cfl数动态调整方法、系统、设备及介质
CN116341421A (zh) * 2023-05-22 2023-06-27 中国空气动力研究与发展中心计算空气动力研究所 高超声速流场数值模拟方法、系统、电子设备及存储介质
CN116341421B (zh) * 2023-05-22 2023-08-25 中国空气动力研究与发展中心计算空气动力研究所 高超声速流场数值模拟方法、系统、电子设备及存储介质
CN116522827A (zh) * 2023-07-04 2023-08-01 北京凌云智擎软件有限公司 一种网格单元边界面的流动变量重构方法、设备及装置
CN116522827B (zh) * 2023-07-04 2023-10-20 北京凌云智擎软件有限公司 一种气动热环境计算的流动变量重构方法、设备及装置

Similar Documents

Publication Publication Date Title
CN115496006A (zh) 一种适用于高超声速飞行器的高精度数值模拟方法
Ghoreyshi et al. Reduced order unsteady aerodynamic modeling for stability and control analysis using computational fluid dynamics
Aeschliman et al. Experimental methodology for computational fluid dynamics code validation
Rizzi et al. Computation of flow around wings based on the Euler equations
CN108763683B (zh) 一种三角函数框架下新weno格式构造方法
CN115238397B (zh) 高超声速飞行器热环境计算方法、装置和计算机设备
Frink et al. Computational aerodynamic modeling tools for aircraft loss of control
CN111859645B (zh) 冲击波求解的改进musl格式物质点法
Morelli et al. Frequency-domain method for automated simulation updates based on flight data
Nabawy et al. Aerodynamic shape optimisation, wind tunnel measurements and CFD analysis of a MAV wing
Kishi et al. Supersonic forward-swept wing design using multifidelity efficient global optimization
Stradtner et al. Introduction to AVT-351: Enhanced Computational Performance and Stability & Control Prediction for NATO Military Vehicles
Millidere Optimal input design and system identification for an agile aircraft
Dean et al. Aircraft stability and control characteristics determined by system identification of CFD simulations
Li Feasibility of supersonic aircraft concepts for low-boom and flight trim constraints
Moran et al. Wind-Tunnel based Free-Flight Testing of a Viscous Optimised Hypersonic Waverider
CN115730533A (zh) 一种侧向喷流干扰流场模拟计算的网格自适应混合判据
Barrett et al. Airfoil shape design and optimization using multifidelity analysis and embedded inverse design
Jirásek et al. Reduced order modeling of X-31 wind tunnel model aerodynamic loads
Nelson et al. Experimental and numerical investigation of flight dynamics of a generic lambda wing configuration
Silva et al. An overview of the semi-span super-sonic transport (S4T) wind-tunnel model program
Giblette Rapid prediction of low-boom and aerodynamic performance of supersonic transport aircraft using panel methods
Righi et al. Uncertainties Quantification of CFD-Based Flutter Prediction
Xu et al. Expert's experience-informed hierarchical kriging method for aerodynamic data modeling
Morton et al. CFD-Based Reduced-Order Models of F-16XL Static and Dynamic Loads Using Kestrel

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination