CN113486453A - 一种可压缩流体的数值模拟方法 - Google Patents

一种可压缩流体的数值模拟方法 Download PDF

Info

Publication number
CN113486453A
CN113486453A CN202111045330.9A CN202111045330A CN113486453A CN 113486453 A CN113486453 A CN 113486453A CN 202111045330 A CN202111045330 A CN 202111045330A CN 113486453 A CN113486453 A CN 113486453A
Authority
CN
China
Prior art keywords
flux
threshold
positive
negative
format
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.)
Granted
Application number
CN202111045330.9A
Other languages
English (en)
Other versions
CN113486453B (zh
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.)
Low Speed Aerodynamics Institute of China Aerodynamics Research and Development Center
Original Assignee
Low Speed Aerodynamics Institute of China Aerodynamics Research and Development Center
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 Low Speed Aerodynamics Institute of China Aerodynamics Research and Development Center filed Critical Low Speed Aerodynamics Institute of China Aerodynamics Research and Development Center
Priority to CN202111045330.9A priority Critical patent/CN113486453B/zh
Publication of CN113486453A publication Critical patent/CN113486453A/zh
Application granted granted Critical
Publication of CN113486453B publication Critical patent/CN113486453B/zh
Withdrawn - After Issue legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • 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
    • 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
    • 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

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • Evolutionary Computation (AREA)
  • Computer Hardware Design (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (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

本发明适用于计算流体力学领域，提供了一种可压缩流体的数值模拟方法，包括：S10：读取计算网格、时
Figure 100004_DEST_PATH_IMAGE001
刻的物理变量分布、边界条件、试验时长；S20：在计算网格上，根据
Figure 662121DEST_PATH_IMAGE001
时刻的物理变量获得守恒变量
Figure 313682DEST_PATH_IMAGE002
和对流通量
Figure 100004_DEST_PATH_IMAGE003
、
Figure 483632DEST_PATH_IMAGE004
、
Figure 100004_DEST_PATH_IMAGE005
；S30：分别计算对流通量
Figure 331372DEST_PATH_IMAGE003
、
Figure 145744DEST_PATH_IMAGE004
、
Figure 437048DEST_PATH_IMAGE005
在物理空间上的数值通量
Figure 376185DEST_PATH_IMAGE006
、
Figure 100004_DEST_PATH_IMAGE007
、
Figure 558774DEST_PATH_IMAGE008
；S40：将守恒变量
Figure 821128DEST_PATH_IMAGE002
和数值通量
Figure 204704DEST_PATH_IMAGE006
、
Figure 959034DEST_PATH_IMAGE007
、
Figure 430466DEST_PATH_IMAGE008
代入离散欧拉方程得到n+1时刻的守恒变量
Figure 100004_DEST_PATH_IMAGE009
及其对应的物理变量；S50：n从1~k遍历，其中k为
Figure 422693DEST_PATH_IMAGE010
时刻试验时长的序号，将守恒变量
Figure 100004_DEST_PATH_IMAGE011
及
Figure 242751DEST_PATH_IMAGE010
时刻的物理变量作为最终的试验结果。通过本发明的方法提高了计算效率、保证了数值计算时的本质无振荡特性、实现了稳定的计算。

Description

一种可压缩流体的数值模拟方法
技术领域
本发明涉及计算流体力学领域 ，尤其是涉及一种可压缩流体的数值模拟方法。
背景技术
研究关于飞行器三维复杂流动的数值模拟问题是目前计算流体力学的前沿热点问题，同时也是工程应用中的重要需求，计算流体力学为低成本地研究真实流动现象提供了巨大支持。和传统的风洞试验相比较，计算流体力学技术具有计算成本低、飞行大气环境仿真度高的优点，可以以较低计算代价获得较为精确的气动特性，为后续飞行器的相关设计提供快速和准确的指导。
对流扩散方程是粘性流体力学及其他实际问题中常见的一类基本控制方程。对于这类方程，除了极少数简单情形，大部分问题目前还无法求得精确解，只能利用数值方法进行数值模拟。
流体控制方程组是指基于连续介质假设由质量守恒、动量守恒、能量守恒定律推导得到的连续方程、动量方程和能量方程组成的方程组，是常规流动问题的核心数学描述。对于不同流动问题，流动控制方程组可以简化省略部分项，或增加其他方程。特别地，对于可压缩无粘流动，流动控制方程组简化为欧拉方程组。
加权本质无振荡格式（Weighted Essentially Non-oscillatory， 简称WENO）方法是近年来广泛流行的一种高分辨率数值方法，用于解决以对流为主的对流扩散方程，特别是双曲守恒律问题。加权本质无振荡格式具有在光滑区可以得到任意高阶精度，在间断区保持稳定的、本质无振荡性质和陡峭的间断过渡等优点。
传统的计算流体力学在求解含激波、接触间断、相界面等间断流场结构的流动问题时，一阶格式会引起间断结构的过度抹平，高阶线性格式则会在间断附近产生非物理的数值振荡，均难以实现有效模拟。加权本质无振荡格式方法对子模板上低阶格式进行非线性加权，既可以实现对变量光滑区域的高精度计算，又可以保证间断附近的本质无振荡特性，捕捉到陡峭的间断过渡，被广泛应用于含激波等间断结构的数值计算中。但经典的加权本质无振荡格式在整个计算区域进行非线性加权，计算量较大；可采用一定的判定准则检测间断，仅在变量间断附近进行加权本质无振荡格式非线性加权计算，提高格式效率。
采用加权本质无振荡格式求解欧拉方程组时，可以采用分量型加权本质无振荡格式，即在物理空间对分量形式方程进行求解，但这样仍会产生数值振荡；也可以采用特征型加权本质无振荡格式，即在特征空间对特征形式方程进行求解，这样可更好地抑制振荡，但涉及特征变换等矩阵计算，计算量过大。
发明内容
本发明实施例的目的是提供一种可压缩流体的数值模拟方法，目的在于解决现有技术中存在的技术问题，以欧拉方程组为例展开说明，包括如下步骤：
步骤S10：读取计算网格、
Figure 353741DEST_PATH_IMAGE001
时刻的物理变量分布、边界条件、试验时长；
步骤S20：在计算网格上，根据
Figure 474144DEST_PATH_IMAGE001
时刻的物理变量获得守恒变量
Figure 785039DEST_PATH_IMAGE002
和对流通量
Figure 524456DEST_PATH_IMAGE003
、
Figure 338829DEST_PATH_IMAGE004
、
Figure 692449DEST_PATH_IMAGE005
，其中，
Figure 100428DEST_PATH_IMAGE003
为笛卡尔坐标系下x方向上的对流通量、
Figure 830487DEST_PATH_IMAGE004
为笛卡尔坐标系下y方向上的对流通量；
Figure 702628DEST_PATH_IMAGE005
为笛卡尔坐标系下z方向上的对流通量；
步骤S30：分别计算对流通量
Figure 102516DEST_PATH_IMAGE003
、
Figure 122425DEST_PATH_IMAGE004
、
Figure 859437DEST_PATH_IMAGE005
在物理空间上的数值通量
Figure 789347DEST_PATH_IMAGE006
、
Figure 422453DEST_PATH_IMAGE007
、
Figure 929658DEST_PATH_IMAGE008
；
步骤S40：将守恒变量
Figure 345727DEST_PATH_IMAGE002
和数值通量
Figure 130143DEST_PATH_IMAGE006
、
Figure 996468DEST_PATH_IMAGE007
、
Figure 928652DEST_PATH_IMAGE008
代入离散欧拉方程得到n+1时刻的守恒变量
Figure 160131DEST_PATH_IMAGE009
及其对应的物理变量；
步骤S50：n从1~k遍历，其中k为
Figure 189267DEST_PATH_IMAGE010
时刻试验时长的序号，将守恒变量
Figure 164176DEST_PATH_IMAGE011
及
Figure 114814DEST_PATH_IMAGE010
时刻的物理变量作为最终的试验结果。
进一步地，步骤S30包括如下步骤：
步骤S301：构建对流通量
Figure 403844DEST_PATH_IMAGE003
、
Figure 21908DEST_PATH_IMAGE004
、
Figure 370980DEST_PATH_IMAGE005
在物理空间上的正负混合通量
Figure 12177DEST_PATH_IMAGE012
、
Figure 963953DEST_PATH_IMAGE013
、
Figure 436522DEST_PATH_IMAGE014
步骤S302：对物理空间上的正负混合通量
Figure 690917DEST_PATH_IMAGE012
、
Figure 819410DEST_PATH_IMAGE013
、
Figure 574877DEST_PATH_IMAGE014
进行间断检测，若检测到正负混合通量
Figure 777319DEST_PATH_IMAGE012
、
Figure 61670DEST_PATH_IMAGE013
、
Figure 208618DEST_PATH_IMAGE014
均未存在间断，将对流通量
Figure 643141DEST_PATH_IMAGE003
、
Figure 90303DEST_PATH_IMAGE004
、
Figure 686500DEST_PATH_IMAGE005
进行分裂和逼近，得到物理空间上的数值通量
Figure 789586DEST_PATH_IMAGE015
、
Figure 621275DEST_PATH_IMAGE016
、
Figure 329468DEST_PATH_IMAGE017
；若检测到正负混合通量
Figure 893305DEST_PATH_IMAGE012
、
Figure 952528DEST_PATH_IMAGE013
、
Figure 853488DEST_PATH_IMAGE014
存在间断，将对流通量
Figure 744083DEST_PATH_IMAGE003
、
Figure 744400DEST_PATH_IMAGE004
、
Figure 494182DEST_PATH_IMAGE005
投影至特征空间，并分裂为正负通量
Figure 198833DEST_PATH_IMAGE018
、
Figure 412776DEST_PATH_IMAGE019
、
Figure 380732DEST_PATH_IMAGE020
，并根据特征空间上的正负通量
Figure 883389DEST_PATH_IMAGE018
、
Figure 860572DEST_PATH_IMAGE019
、
Figure 929022DEST_PATH_IMAGE020
计算得到物理空间上的数值通量
Figure 474404DEST_PATH_IMAGE015
、
Figure 57832DEST_PATH_IMAGE016
、
Figure 573127DEST_PATH_IMAGE017
。
进一步地，步骤S302中，根据特征空间上的正负通量
Figure 964925DEST_PATH_IMAGE018
、
Figure 477946DEST_PATH_IMAGE019
、
Figure 79829DEST_PATH_IMAGE020
，计算得到物理空间上的数值通量
Figure 274181DEST_PATH_IMAGE015
、
Figure 113961DEST_PATH_IMAGE016
、
Figure 125780DEST_PATH_IMAGE017
的步骤如下：
步骤S3021：对特征空间上的正负通量
Figure 621483DEST_PATH_IMAGE018
、
Figure 619526DEST_PATH_IMAGE019
、
Figure 844971DEST_PATH_IMAGE020
进行间断检测，若未检测到正负通量
Figure 496532DEST_PATH_IMAGE018
、
Figure 745111DEST_PATH_IMAGE019
、
Figure 546845DEST_PATH_IMAGE020
存在间断，则采用加权本质无振荡格式对应的高阶线性格式对正负通量
Figure 361217DEST_PATH_IMAGE018
、
Figure 855783DEST_PATH_IMAGE019
、
Figure 326079DEST_PATH_IMAGE020
进行逼近，得到特征空间上的数值通量；若检测到正负通量
Figure 524979DEST_PATH_IMAGE018
、
Figure 193858DEST_PATH_IMAGE019
、
Figure 859326DEST_PATH_IMAGE020
存在间断，则采用加权本质无振荡格式对正负通量
Figure 816917DEST_PATH_IMAGE018
、
Figure 85088DEST_PATH_IMAGE019
、
Figure 608473DEST_PATH_IMAGE020
进行逼近，得到特征空间上的数值通量；
步骤S3022：将特征空间上的数值通量投影回物理空间，得到物理空间上的数值通量
Figure 167544DEST_PATH_IMAGE015
、
Figure 674748DEST_PATH_IMAGE016
、
Figure 418714DEST_PATH_IMAGE017
。
进一步地，步骤S302中，若检测到正负混合通量
Figure 734288DEST_PATH_IMAGE012
、
Figure 272717DEST_PATH_IMAGE013
、
Figure 408163DEST_PATH_IMAGE014
均未存在间断，将对流通量
Figure 752557DEST_PATH_IMAGE021
、
Figure 657059DEST_PATH_IMAGE022
、进行分裂和逼近，得到物理空间上的数值通量
Figure 163127DEST_PATH_IMAGE023
、
Figure 582607DEST_PATH_IMAGE024
、
Figure 996271DEST_PATH_IMAGE025
，获得物理空间上数值通量
Figure 489700DEST_PATH_IMAGE023
、
Figure 963407DEST_PATH_IMAGE024
、
Figure 401341DEST_PATH_IMAGE025
的方法如下：
采用通量分裂格式将对流通量
Figure 25221DEST_PATH_IMAGE021
、
Figure 638736DEST_PATH_IMAGE022
、
Figure 17764DEST_PATH_IMAGE026
分裂为正负通量
Figure 146257DEST_PATH_IMAGE027
、
Figure 777090DEST_PATH_IMAGE028
、
Figure 104166DEST_PATH_IMAGE029
；对正负通量
Figure 654096DEST_PATH_IMAGE027
、
Figure 473148DEST_PATH_IMAGE028
、
Figure 438830DEST_PATH_IMAGE029
采用加权本质无振荡格式对应的高阶线性格式进行逼近，获得物理空间中的数值通量
Figure 885992DEST_PATH_IMAGE023
、
Figure 341244DEST_PATH_IMAGE024
、
Figure 647591DEST_PATH_IMAGE025
。
进一步地，步骤S302，将对流通量
Figure 416964DEST_PATH_IMAGE021
、
Figure 984212DEST_PATH_IMAGE022
、
Figure 751310DEST_PATH_IMAGE026
投影至特征空间，并分裂为正负通量
Figure 669588DEST_PATH_IMAGE030
、
Figure 773810DEST_PATH_IMAGE031
、
Figure 133247DEST_PATH_IMAGE032
的方法如下：
将对流通量
Figure 336827DEST_PATH_IMAGE021
、
Figure 211242DEST_PATH_IMAGE022
、
Figure 587996DEST_PATH_IMAGE026
由物理空间投影至特征空间，得到特征空间上对流通量
Figure 5202DEST_PATH_IMAGE033
、
Figure 973158DEST_PATH_IMAGE034
、
Figure 334870DEST_PATH_IMAGE035
；采用通量分裂格式对特征空间上的对流通量
Figure 187419DEST_PATH_IMAGE033
、
Figure 521449DEST_PATH_IMAGE034
、
Figure 191464DEST_PATH_IMAGE035
分裂得到正负通量
Figure 774892DEST_PATH_IMAGE030
、
Figure 165554DEST_PATH_IMAGE031
、
Figure 681986DEST_PATH_IMAGE032
。
进一步地，步骤S302中和步骤S3021中的通量间断检测的方法如下：
获取对流通量正负混合通量
Figure 195006DEST_PATH_IMAGE036
、
Figure 406676DEST_PATH_IMAGE037
、
Figure 991241DEST_PATH_IMAGE038
和特征空间上的正负通量
Figure 831021DEST_PATH_IMAGE030
、
Figure 514943DEST_PATH_IMAGE031
、
Figure 213909DEST_PATH_IMAGE032
中的任意一个通量的分量，并将所述分量作为目标分量g；将g在计算网格的网格单元i中心处的离散值记为g i ，遍历所有计算网格单元i，对目标分量g在网格单元i附近依次进行间断检测；i=1,2,...,N，N为网格单元的总个数；
获取目标分量g在网格单元i对应的加权无本质振荡格式大模板S whole中的最大斜率值S i,max；
获取目标分量g在网格单元i对应的加权本质无振荡格式阈值模板S threshold中的斜率阈值S i , threshold，其中加权无本质振荡格式阈值模板S threshold为加权无本质振荡格式大模板S whole的子模板中任一迎风子模板；
当S i,max<S i , threshold时，则检测到目标分量g在网格单元i附近不存在间断；当S i,max
Figure 336586DEST_PATH_IMAGE039
S i , threshold时，则检测到目标分量g在网格单元i附近存在间断。
进一步地，获取目标分量g在网格单元i对应的加权无本质振荡格式大模板S whole中的最大斜率值S i,max的计算方法如下：
在网格单元i对应的加权无本质振荡格式的大模板S whole上对目标分量g采用加权本质无振荡格式对应的高阶线性格式进行逼近，得到逼近多项式函数
Figure 437397DEST_PATH_IMAGE040
；
计算逼近多项式函数
Figure 88958DEST_PATH_IMAGE040
在大模板S whole上的斜率的绝对值的最大值S i,max。
进一步地，获取目标分量g在网格单元i对应的加权本质无振荡格式阈值模板S threshold中的斜率阈值S i , threshold的计算方法如下：
选取网格单元i对应的加权本质无振荡格式中的阈值模板S threshold；
在加权本质无振荡格式阈值模板S threshold中计算基于网格最小可分辨波长的目标分量g的斜率阈值S i , threshold，其中，网格最小可分辨波长为网格单元长度的两倍。
进一步地，计算逼近多项式函数
Figure 337537DEST_PATH_IMAGE040
在大模板S whole上的斜率绝对值的最大值S i,max的公式如下：
Figure 998326DEST_PATH_IMAGE041
其中，
Figure 953643DEST_PATH_IMAGE042
为逼近多项式函数
Figure 244947DEST_PATH_IMAGE040
在大模板S whole上的导数。
进一步地，计算斜率阈值S i , threshold的公式如下：
Figure 43139DEST_PATH_IMAGE043
其中，
Figure 976460DEST_PATH_IMAGE044
为网格单元长度，r和b为加权本质无振荡格式中阈值模板S threshold对应的整数，m为网格单元编号，g m+1 和g m 为网格单元的中心位置处目标分量g的离散值，
Figure 798003DEST_PATH_IMAGE045
为正数小量。
本发明产生的技术效果至少具有以下方面：
1）通过本发明提供的可压缩流体的数值模拟方法，该数值模拟方法在求解欧拉方程组的过程中进行了两次间断检测，在未检测到间断的区域采用加权本质无振荡格式对应的高阶线性格式进行计算，简化了非必须的加权本质无振荡格式非线性加权计算，实现了加权本质无振荡格式整体计算效率的提高；在检测到间断的区域采用加权本质无振荡格式计算，保证数值计算时的本质无振荡特性，实现稳定的计算。需要说明的是本发明的方法主要适用于有限差分加权本质无振荡格式，但易于扩展至有限体积加权本质无振荡格式。
2）本发明的可压缩流体的数值模拟方法中，在第一次对对流通量在物理空间上进行间断检测时，对于检测到的正负混合通量均未存在间断时，直接在物理空间中对对流通量进行分裂和逼近计算，不对对流通量进行特征投影等耗时巨大的矩阵运算，减小了试验时的计算量，提高了试验效率。
3）本发明的可压缩流体的数值模拟方法中，在物理空间和特征空间进行了两次间断检测后的计算过程中，仅针对在特征空间中检测出间断的通量进行加权本质无振荡格式的逼近，大大减小了加权本质无振荡非线性加权计算的计算量。
4）本发明的可压缩流体的数值模拟方法中，在物理空间进行间断检测时，构建了正负混合通量，进行间断检测时仅是对正负混合通量进行间断检测，而不是对原对流通量中的各个分量分别进行间断检测，减少了多次间断检测的计算量。
5）本发明的可压缩流体的数值模拟方法中，在进行通量间断检测的方法中将斜率阈值作为间断检测的标准，除了正数小量外，不包含任何经验参数，具有更好的普适性，同时，斜率阈值与所研究的网格单元所处的阈值模板内通量差值的绝对值正相关，对于不同的网格单元取值不同，具有更好的自适应性。
附图说明
图1是本发明可压缩流体的数值模拟方法的总体流程图 ；
图2是本发明中通量进行间断检测的流程图；
图3是本发明中 5阶加权本质无振荡格式大模板与高阶线性格式逼近多项式函数示意图；
图4是本发明5阶加权本质无振荡格式各子模板与斜率阈值示意图。
具体实施方式
以下的说明提供了许多不同的实施例、或是例子，用来实施本发明的不同特征。以下特定例子所描述的元件和排列方式，仅用来精简的表达本发明，其仅作为例子，而并非用以限制本发明。
如图1-4所示，本发明实施例的目的是提供一种可压缩流体的数值模拟方法，目的在于解决现有技术中存在的技术问题，包括如下步骤：
步骤S10：读取计算网格、时
Figure 260208DEST_PATH_IMAGE046
刻的物理变量分布、边界条件、试验时长；
步骤S20：在计算网格上，根据
Figure 280117DEST_PATH_IMAGE046
时刻的物理变量获得守恒变量
Figure 423653DEST_PATH_IMAGE047
和对流通量
Figure 947039DEST_PATH_IMAGE021
、
Figure 642462DEST_PATH_IMAGE022
、
Figure 87350DEST_PATH_IMAGE026
，其中，
Figure 34577DEST_PATH_IMAGE021
为笛卡尔坐标系下x方向上的对流通量、
Figure 678048DEST_PATH_IMAGE022
为笛卡尔坐标系下y方向上的对流通量；
Figure 482056DEST_PATH_IMAGE026
为笛卡尔坐标系下z方向上的对流通量；
步骤S30：分别计算对流通量
Figure 945399DEST_PATH_IMAGE021
、
Figure 696317DEST_PATH_IMAGE022
、
Figure 459874DEST_PATH_IMAGE026
在物理空间上的数值通量
Figure 372466DEST_PATH_IMAGE048
、
Figure 791946DEST_PATH_IMAGE049
、
Figure 674451DEST_PATH_IMAGE050
；
步骤S40：将守恒变量
Figure 558094DEST_PATH_IMAGE047
和数值通量
Figure 376008DEST_PATH_IMAGE048
、
Figure 345101DEST_PATH_IMAGE049
、
Figure 234560DEST_PATH_IMAGE050
代入离散欧拉方程得到n+1时刻的守恒变量
Figure 582496DEST_PATH_IMAGE051
及其对应的物理变量；
步骤S50：n从1~k遍历，其中k为
Figure 695945DEST_PATH_IMAGE052
时刻试验时长的序号，将守恒变量
Figure 152334DEST_PATH_IMAGE053
及
Figure 579905DEST_PATH_IMAGE052
时刻的物理变量作为最终的试验结果。
对流扩散方程是粘性流体力学及其他实际问题中常见的一类基本控制方程。对于这类方程，除了极少数简单情形，大部分问题目前还无法求得精确解，只能利用数值方法进行数值模拟。
流体控制方程组是常规流动问题的核心数学描述，具体是指基于连续介质假设由质量守恒、动量守恒、能量守恒定律推导得到的连续方程、动量方程和能量方程组成的方程组。对于不同流动问题，流动控制方程组可以简化省略部分项，或者增加其他方程。特别地，对于可压缩无粘流动，流动控制方程组简化为欧拉方程组，欧拉方程组的形式如下：
Figure 579085DEST_PATH_IMAGE054
Figure 269960DEST_PATH_IMAGE055
，
Figure 620170DEST_PATH_IMAGE056
，
Figure 179327DEST_PATH_IMAGE057
，
Figure 501855DEST_PATH_IMAGE058
，
其中：U为守恒变量；F、G、H分别为笛卡尔坐标系下x、y、z三个方向上的对流通量；S为源项；
Figure 160370DEST_PATH_IMAGE059
为密度；u、v、w分别为笛卡尔坐标系下x、y、z方向上的速度分量；E为单位质量流体总能；p为压力
在均匀结构网格上，欧拉方程组的半离散形式如下：
Figure DEST_PATH_IMAGE061A
其中，i、j、k为网格单元编号；
Figure 466717DEST_PATH_IMAGE062
、
Figure 32828DEST_PATH_IMAGE063
、
Figure 741021DEST_PATH_IMAGE064
为x、y、z方向上对流通量对应的离散后的数值通量，
Figure 570437DEST_PATH_IMAGE065
、
Figure 98501DEST_PATH_IMAGE066
、
Figure 468302DEST_PATH_IMAGE067
为x、y、z方向的网格单元长度，S i,j,k 为网格单元中心源项值。
计算网格具体为：根据设计时的任务要求，构建需要设计对象的模型，所述模型为不同种类飞行器的模型，并建立复杂流场计算所需要的网格，由给定的网格得到需要的网格单元信息，网格单元信息包括网格尺度、节点坐标等。
边界条件指的是计算区域边界的物理变量的分布和取值。
试验时长指的是本次试验需要的时间。
具体地，物理变量指的是
Figure 890057DEST_PATH_IMAGE059
、u、v、w、p，物理变量通过简单的计算后可以获得守恒变量
Figure 93636DEST_PATH_IMAGE047
和对流通量
Figure 905734DEST_PATH_IMAGE021
、
Figure 344806DEST_PATH_IMAGE022
、
Figure 496432DEST_PATH_IMAGE026
。
上述方案中，根据物理空间上的对流通量
Figure 464388DEST_PATH_IMAGE021
、
Figure 91679DEST_PATH_IMAGE022
、
Figure 6545DEST_PATH_IMAGE026
计算出物理空间上的数值数值通量
Figure 278258DEST_PATH_IMAGE023
、
Figure 948274DEST_PATH_IMAGE024
、
Figure 266122DEST_PATH_IMAGE025
，然后将守恒变量
Figure 984680DEST_PATH_IMAGE047
和数值通量
Figure 376478DEST_PATH_IMAGE023
、
Figure 951816DEST_PATH_IMAGE024
、
Figure 163485DEST_PATH_IMAGE025
代入离散欧拉方程得到n+1时刻的守恒变量
Figure 216892DEST_PATH_IMAGE068
及其对应的物理变量；n从1~k遍历，其中k为
Figure 259934DEST_PATH_IMAGE069
时刻试验时长的序号，将守恒变量
Figure 271753DEST_PATH_IMAGE070
及
Figure 959000DEST_PATH_IMAGE069
时刻的物理变量作为最终的试验结果。
物理空间上的对流通量
Figure 753780DEST_PATH_IMAGE021
、
Figure 713646DEST_PATH_IMAGE022
、
Figure 896366DEST_PATH_IMAGE026
计算出物理空间上的数值数值通量
Figure 817048DEST_PATH_IMAGE048
、
Figure 743416DEST_PATH_IMAGE049
、
Figure 761051DEST_PATH_IMAGE050
的过程中，采用了两次间断检测，包括物理空间间断检测和特征空间间断检测，在物理空间间断检测时，首先构建对流通量
Figure 990038DEST_PATH_IMAGE021
、
Figure 522650DEST_PATH_IMAGE022
、
Figure 721550DEST_PATH_IMAGE026
在物理空间上的正负混合通量
Figure 593692DEST_PATH_IMAGE036
、
Figure 993580DEST_PATH_IMAGE037
、
Figure 13489DEST_PATH_IMAGE038
，然后对正负混合通量
Figure 422604DEST_PATH_IMAGE036
、
Figure 149252DEST_PATH_IMAGE037
、
Figure 313517DEST_PATH_IMAGE038
进行间断检测，若检测到正负混合通量
Figure 86301DEST_PATH_IMAGE036
、
Figure 767949DEST_PATH_IMAGE037
、
Figure 676999DEST_PATH_IMAGE038
均未存在间断，则直接对
Figure 481007DEST_PATH_IMAGE021
、
Figure 85295DEST_PATH_IMAGE022
、
Figure 429689DEST_PATH_IMAGE026
进行分裂和逼近，得到物理空间上的数值通量
Figure 193245DEST_PATH_IMAGE023
、
Figure 168155DEST_PATH_IMAGE024
、
Figure 525318DEST_PATH_IMAGE025
，若检测到正负混合通量
Figure 673402DEST_PATH_IMAGE036
、
Figure 25886DEST_PATH_IMAGE037
、
Figure 171697DEST_PATH_IMAGE038
存在间断，将对流通量
Figure 16156DEST_PATH_IMAGE021
、
Figure 967932DEST_PATH_IMAGE022
、
Figure 581447DEST_PATH_IMAGE026
投影至特征空间，并分裂为正负通量
Figure 694896DEST_PATH_IMAGE030
、
Figure 354548DEST_PATH_IMAGE031
、
Figure 719801DEST_PATH_IMAGE032
，并根据特征空间上的正负通量
Figure 781298DEST_PATH_IMAGE030
、
Figure 800069DEST_PATH_IMAGE031
、
Figure 353542DEST_PATH_IMAGE032
，对特征空间上的正负通量
Figure 912699DEST_PATH_IMAGE030
、
Figure 297544DEST_PATH_IMAGE031
、
Figure 159321DEST_PATH_IMAGE032
进行间断检测，计算得到物理空间上的数值通量
Figure 59144DEST_PATH_IMAGE023
、
Figure 31779DEST_PATH_IMAGE024
、
Figure 536709DEST_PATH_IMAGE025
。
因此，通过本发明提供的可压缩流体的数值模拟方法，该数值模拟方法在求解欧拉方程组的过程中进行了两次间断检测，在未检测到间断的区域采用加权本质无振荡格式对应的高阶线性格式进行计算，简化了非必须的加权本质无振荡非线性加权计算，实现了加权本质无振荡格式整体计算效率的提高；在检测到间断的区域采用加权本质无振荡格式计算，保证数值计算时的本质无振荡特性，实现稳定的计算。需要说明的是本发明的方法主要适用于有限差分加权本质无振荡格式，但易于扩展至有限体积加权本质无振荡格式。
进一步地，步骤S30包括如下步骤：
步骤S301：构建对流通量
Figure 428442DEST_PATH_IMAGE021
、
Figure 956507DEST_PATH_IMAGE022
、
Figure 60729DEST_PATH_IMAGE026
在物理空间上的正负混合通量
Figure 889008DEST_PATH_IMAGE036
、
Figure 889325DEST_PATH_IMAGE037
、
Figure 904685DEST_PATH_IMAGE038
；
步骤S302：对物理空间上的正负混合通量
Figure 78177DEST_PATH_IMAGE036
、
Figure 292121DEST_PATH_IMAGE037
、
Figure 260077DEST_PATH_IMAGE038
进行间断检测，若检测到正负混合通量
Figure 762734DEST_PATH_IMAGE036
、
Figure 739917DEST_PATH_IMAGE037
、
Figure 23348DEST_PATH_IMAGE038
均未存在间断，将对流通量
Figure 365468DEST_PATH_IMAGE021
、
Figure 480054DEST_PATH_IMAGE022
、
Figure 995349DEST_PATH_IMAGE026
进行分裂和逼近，得到物理空间上的数值通量
Figure 855989DEST_PATH_IMAGE023
、
Figure 900169DEST_PATH_IMAGE024
、
Figure 111838DEST_PATH_IMAGE025
；若检测到正负混合通量
Figure 696403DEST_PATH_IMAGE036
、
Figure 145970DEST_PATH_IMAGE037
、
Figure 626630DEST_PATH_IMAGE038
存在间断，将对流通量
Figure 450230DEST_PATH_IMAGE021
、
Figure 448273DEST_PATH_IMAGE022
、
Figure 345822DEST_PATH_IMAGE026
投影至特征空间，并分裂为正负通量
Figure 528541DEST_PATH_IMAGE030
、
Figure 714803DEST_PATH_IMAGE031
、
Figure 313275DEST_PATH_IMAGE032
，并根据特征空间上的正负通量
Figure 658805DEST_PATH_IMAGE030
、
Figure 622213DEST_PATH_IMAGE031
、
Figure 889247DEST_PATH_IMAGE032
计算得到物理空间上的数值通量
Figure 619305DEST_PATH_IMAGE023
、
Figure 491446DEST_PATH_IMAGE024
、
Figure 891335DEST_PATH_IMAGE025
。
上述方案中，特征空间上的正负混合通量
Figure 114506DEST_PATH_IMAGE036
的计算方式如下：
Figure 258042DEST_PATH_IMAGE071
其中，
Figure 47007DEST_PATH_IMAGE072
，c i 为网格单元i中心处音速。
正负混合通量
Figure 883376DEST_PATH_IMAGE037
、
Figure 328263DEST_PATH_IMAGE038
的计算方式与正负混合通量
Figure 134545DEST_PATH_IMAGE036
的计算方式相似。
在间断检测时，检测到正负混合通量
Figure 918962DEST_PATH_IMAGE036
、
Figure 254128DEST_PATH_IMAGE037
、
Figure 186312DEST_PATH_IMAGE038
均未存在间断指的是正混合通量
Figure 796285DEST_PATH_IMAGE073
、
Figure 435208DEST_PATH_IMAGE074
、
Figure 738013DEST_PATH_IMAGE075
和负混合通量
Figure 891914DEST_PATH_IMAGE076
、
Figure 774419DEST_PATH_IMAGE077
、
Figure 533428DEST_PATH_IMAGE078
两者都没有检测到间断，此时直接对对流通量
Figure 741555DEST_PATH_IMAGE021
、
Figure 586015DEST_PATH_IMAGE022
、
Figure 209894DEST_PATH_IMAGE026
进行分裂和逼近，得到物理空间上的数值通量
Figure 682464DEST_PATH_IMAGE023
、
Figure 61492DEST_PATH_IMAGE024
、
Figure 393248DEST_PATH_IMAGE025
；当正混合通量
Figure 820818DEST_PATH_IMAGE073
、
Figure 413473DEST_PATH_IMAGE079
、
Figure 838770DEST_PATH_IMAGE080
和负混合通量
Figure 251296DEST_PATH_IMAGE076
、
Figure 482558DEST_PATH_IMAGE077
、
Figure 195299DEST_PATH_IMAGE078
中任何一个检测到间断，则说明正负混合通量
Figure 525917DEST_PATH_IMAGE036
、
Figure 956898DEST_PATH_IMAGE037
、
Figure 726271DEST_PATH_IMAGE038
存在间断，此时对流通量
Figure 27939DEST_PATH_IMAGE021
、
Figure 795038DEST_PATH_IMAGE022
、
Figure 447737DEST_PATH_IMAGE026
投影至特征空间进行进一步的间断检测。
上述方案中，在第一次对对流通量
Figure 9082DEST_PATH_IMAGE021
、
Figure 306202DEST_PATH_IMAGE022
、
Figure 368836DEST_PATH_IMAGE026
在物理空间上进行间断检测时，对于检测到的正负混合通量
Figure 977672DEST_PATH_IMAGE036
、
Figure 557689DEST_PATH_IMAGE037
、
Figure 771632DEST_PATH_IMAGE038
均未存在间断时，直接在物理空间中对对流通量
Figure 270747DEST_PATH_IMAGE021
、
Figure 507824DEST_PATH_IMAGE022
、
Figure 750587DEST_PATH_IMAGE026
进行分裂和逼近计算，不用将对流通量
Figure 615775DEST_PATH_IMAGE021
、
Figure 957894DEST_PATH_IMAGE022
、
Figure 682268DEST_PATH_IMAGE026
进行特征投影等耗时巨大的矩阵运算，减小了试验时的计算量，提高了试验效率。
在物理空间对对流通量
Figure 463142DEST_PATH_IMAGE021
、
Figure 651678DEST_PATH_IMAGE022
、
Figure 367961DEST_PATH_IMAGE026
进行间断检测时，构建了正负混合通量
Figure 438685DEST_PATH_IMAGE036
、
Figure 757671DEST_PATH_IMAGE037
、
Figure 3976DEST_PATH_IMAGE038
，进行间断检测时仅是对正负混合通量
Figure 687898DEST_PATH_IMAGE036
、
Figure 511498DEST_PATH_IMAGE037
、
Figure 368595DEST_PATH_IMAGE038
进行间断检测，而不是针对原对流通量
Figure 203827DEST_PATH_IMAGE021
、
Figure 58651DEST_PATH_IMAGE022
、
Figure 635125DEST_PATH_IMAGE026
中的各个分量分别进行间断检测，减少了多次间断检测的计算量。
进一步地，步骤S302中，根据特征空间上的正负通量
Figure 171280DEST_PATH_IMAGE030
、
Figure 720073DEST_PATH_IMAGE031
、
Figure 339273DEST_PATH_IMAGE032
，计算得到物理空间上的数值通量
Figure 809569DEST_PATH_IMAGE023
、
Figure 414994DEST_PATH_IMAGE024
、
Figure 83872DEST_PATH_IMAGE025
的步骤如下：
步骤S3021：对特征空间上的正负通量
Figure 546078DEST_PATH_IMAGE030
、
Figure 300407DEST_PATH_IMAGE031
、
Figure 443944DEST_PATH_IMAGE032
进行间断检测，若未检测到正负通量
Figure 232908DEST_PATH_IMAGE030
、
Figure 69277DEST_PATH_IMAGE031
、
Figure 248586DEST_PATH_IMAGE032
存在间断，则采用加权本质无振荡格式对应的高阶线性格式对正负通量
Figure 789289DEST_PATH_IMAGE030
、
Figure 963918DEST_PATH_IMAGE031
、
Figure 440030DEST_PATH_IMAGE032
进行逼近，得到特征空间上的数值通量；若检测到正负通量
Figure 700110DEST_PATH_IMAGE030
、
Figure 982187DEST_PATH_IMAGE031
、
Figure 621109DEST_PATH_IMAGE032
存在间断，则采用加权本质无振荡格式对正负通量
Figure 127177DEST_PATH_IMAGE030
、
Figure 874553DEST_PATH_IMAGE031
、
Figure 960321DEST_PATH_IMAGE032
进行逼近，得到特征空间上的数值通量；
步骤S3022：将特征空间上的数值通量投影回物理空间，得到物理空间上的数值通量
Figure 719330DEST_PATH_IMAGE023
、
Figure 927457DEST_PATH_IMAGE024
、
Figure 771916DEST_PATH_IMAGE025
。
上述方案中，在对特征空间上的正负通量
Figure 926954DEST_PATH_IMAGE030
、
Figure 540469DEST_PATH_IMAGE031
、
Figure 591602DEST_PATH_IMAGE032
进行间断检测时，当未检测到间断时采用加权本质无振荡格式对应的高阶线性格式对正负通量
Figure 47991DEST_PATH_IMAGE030
、
Figure 413244DEST_PATH_IMAGE031
、
Figure 740320DEST_PATH_IMAGE032
进行逼近，当检测到间断时，则采用加权本质无振荡格式对正负通量
Figure 555830DEST_PATH_IMAGE030
、
Figure 171619DEST_PATH_IMAGE031
、
Figure 340563DEST_PATH_IMAGE032
进行逼近，在此过程中，仅针对在特征空间上检测出间断的通量进行加权本质无振荡格式的逼近，大大减小了加权本质无振荡非线性加权计算的计算量。
进一步地，步骤S302中，若检测到正负混合通量
Figure 787725DEST_PATH_IMAGE036
、
Figure 457958DEST_PATH_IMAGE037
、
Figure 357781DEST_PATH_IMAGE038
均未存在间断，将对流通量
Figure 330416DEST_PATH_IMAGE021
、
Figure 632085DEST_PATH_IMAGE022
、
Figure 664763DEST_PATH_IMAGE026
进行分裂和逼近，得到物理空间上的数值通量
Figure 255144DEST_PATH_IMAGE023
、
Figure 359366DEST_PATH_IMAGE024
、
Figure 46699DEST_PATH_IMAGE025
，获得物理空间上数值通量
Figure 984700DEST_PATH_IMAGE023
、
Figure 796798DEST_PATH_IMAGE024
、
Figure 501449DEST_PATH_IMAGE025
的方法如下：
采用通量分裂格式将对流通量
Figure 653075DEST_PATH_IMAGE021
、
Figure 621031DEST_PATH_IMAGE022
、
Figure 248322DEST_PATH_IMAGE026
分裂为正负通量
Figure 897609DEST_PATH_IMAGE027
、
Figure 169321DEST_PATH_IMAGE028
、
Figure 839337DEST_PATH_IMAGE029
；对正负通量
Figure 626028DEST_PATH_IMAGE027
、
Figure 875743DEST_PATH_IMAGE028
、
Figure 267542DEST_PATH_IMAGE029
采用加权本质无振荡格式对应的高阶线性格式进行逼近，获得物理空间中的数值通量
Figure 842879DEST_PATH_IMAGE023
、
Figure 320128DEST_PATH_IMAGE024
、
Figure 576797DEST_PATH_IMAGE025
。
上述方案中，通量分裂格式优选为Lax–Friedrichs格式，也可以选用其他的通量分裂格式，在此不做限制。
进一步地，步骤S302，将对流通量
Figure 947736DEST_PATH_IMAGE021
、
Figure 428396DEST_PATH_IMAGE022
、
Figure 861782DEST_PATH_IMAGE026
投影至特征空间，并分裂为正负通量
Figure 922142DEST_PATH_IMAGE030
、
Figure 288532DEST_PATH_IMAGE031
、
Figure 471252DEST_PATH_IMAGE032
的方法如下：
将对流通量
Figure 719831DEST_PATH_IMAGE021
、
Figure 521565DEST_PATH_IMAGE022
、
Figure 335937DEST_PATH_IMAGE026
由物理空间投影至特征空间，得到特征空间上对流通量
Figure 830503DEST_PATH_IMAGE033
、
Figure 97537DEST_PATH_IMAGE034
、
Figure 499699DEST_PATH_IMAGE035
；采用通量分裂格式对特征空间上的对流通量
Figure 575103DEST_PATH_IMAGE033
、
Figure 568466DEST_PATH_IMAGE034
、
Figure 73528DEST_PATH_IMAGE035
分裂得到正负通量
Figure 607278DEST_PATH_IMAGE030
、
Figure 333925DEST_PATH_IMAGE031
、
Figure 904715DEST_PATH_IMAGE032
。
上述方案中，通量分裂格式优选为Lax–Friedrichs格式，也可以选用其他的通量分裂格式，在此不做限制。
进一步地，步骤S302中和步骤S3021中的通量间断检测的方法如下：
获取对流通量正负混合通量
Figure 411920DEST_PATH_IMAGE036
、
Figure 359147DEST_PATH_IMAGE037
、
Figure 940301DEST_PATH_IMAGE038
和特征空间上的正负通量
Figure 72205DEST_PATH_IMAGE030
、
Figure 676493DEST_PATH_IMAGE031
、
Figure 20887DEST_PATH_IMAGE032
中的任意一个通量的分量，并将所述分量作为目标分量g；将g在计算网格的网格单元i中心处的离散值记为g i ，遍历所有计算网格单元i，对目标分量g在网格单元i附近依次进行间断检测；i=1,2,...,N，N为网格单元的总个数；
获取目标分量g在网格单元i对应的加权无本质振荡格式大模板S whole中的最大斜率值S i,max；
获取目标分量g在网格单元i对应的加权本质无振荡格式阈值模板S threshold中的斜率阈值S i , threshold，其中加权无本质振荡格式阈值模板S threshold为加权无本质振荡格式大模板S whole的子模板中任一迎风子模板；
当S i,max<S i , threshold时，则检测到目标分量g在网格单元i附近不存在间断；当S i,max
Figure 722126DEST_PATH_IMAGE039
S i , threshold时，则检测到目标分量g在网格单元i附近存在间断。
进一步地，获取目标分量g在网格单元i对应的加权无本质振荡格式大模板S whole中的最大斜率值S i,max的计算方法如下：
在网格单元i对应的加权无本质振荡格式的大模板S whole上对目标分量g采用加权本质无振荡格式对应的高阶线性格式进行逼近，得到逼近多项式函数
Figure 759353DEST_PATH_IMAGE040
；
计算逼近多项式函数
Figure 116516DEST_PATH_IMAGE040
在大模板S whole上的斜率的绝对值的最大值S i,max。
进一步地，获取目标分量g在网格单元i对应的加权本质无振荡格式阈值模板S threshold中的斜率阈值S i , threshold的计算方法如下：
选取网格单元i对应的加权本质无振荡格式中的阈值模板S threshold；
在加权本质无振荡格式阈值模板S threshold中计算基于网格最小可分辨波长的目标分量g的斜率阈值S i , threshold，其中，网格最小可分辨波长为网格单元长度的两倍。
进一步地，计算逼近多项式函数
Figure 530180DEST_PATH_IMAGE040
在大模板S whole上的斜率绝对值的最大值S i,max的公式如下：
Figure 74207DEST_PATH_IMAGE041
其中，
Figure 751176DEST_PATH_IMAGE042
为逼近多项式函数
Figure 595635DEST_PATH_IMAGE040
在大模板S whole上的导数。
进一步地，计算斜率阈值S i , threshold的公式如下：
Figure 547411DEST_PATH_IMAGE043
其中，
Figure 160926DEST_PATH_IMAGE044
为网格单元长度，r和b为加权本质无振荡格式中阈值模板S threshold对应的整数，m为网格单元编号，g m+1 和g m 为网格单元的中心位置处目标分量g的离散值，
Figure 477638DEST_PATH_IMAGE045
为正数小量。
上述方案中，如图3-4所示，本发明实施例中以5阶的加权本质无振荡格式为例进行说明，其中5阶加权本质无振荡格式构造模板包括大模板
Figure 668448DEST_PATH_IMAGE081
和三个子模板S 0 、S 1 、S 2 ，则可取
Figure 892756DEST_PATH_IMAGE081
为大模板S whole，S 1 为阈值模板；S threshold其中x为横坐标，x i 为网格单元的中心，g i 为目标分量g在网格单元的中心x i 处的值。
最大斜率阈值S i,max的计算方法：在网格单元i应的加权无本质振荡格式的大模板S whole上对目标分量g采用加权本质无振荡格式对应的高阶线性格式进行逼近，得到逼近多项式函数
Figure 626356DEST_PATH_IMAGE040
，在附3图中，虚线表示的为逼近多项式函数
Figure 113969DEST_PATH_IMAGE040
；然后通过求导计算逼近多项式函数
Figure 792075DEST_PATH_IMAGE040
的斜率值
Figure 961020DEST_PATH_IMAGE042
，取斜率值
Figure 142602DEST_PATH_IMAGE042
的绝对值的最大值作为目标分量g在网格单元i对应的加权无本质振荡格式大模板S whole中的最大斜率值S i,max；
在加权本质无振荡格式所有子模板中，取模板中心距界面
Figure 129013DEST_PATH_IMAGE082
最近的迎风子模板为阈值模板S threshold，在S threshold上计算基于网格最小可分辨波长（即2倍的网格单元长度）的斜率阈值S i , threshold。记阈值模板S threshold为
Figure 232098DEST_PATH_IMAGE083
，
Figure 939154DEST_PATH_IMAGE084
为第i个网格单元，r和b为与所采用格式相关的整数，在附图4中，弯曲的虚线的斜率的绝对值的最大值表示的是斜率阈值S i , threshold。
本发明中间断检测方法的基本原理为：采用高阶线性格式对应的逼近多项式函数
Figure 506402DEST_PATH_IMAGE040
进行数值模拟时，在间断附近会产生非物理的数值振荡，这些振荡伴随有变量的大斜率变化，这使得可以用
Figure 132555DEST_PATH_IMAGE040
的斜率大小表征变量是否出现间断；而在构建的阈值模板S threshold上的斜率阈值S i , threshold时，S i , threshold为波长为网格最小可分辨波长（即2倍的网格单元长度）的正弦波在幅值为
Figure 457357DEST_PATH_IMAGE085
时的斜率最大值，用以表示S threshold中允许出现的最大斜率。这里幅值之所以取迎风的阈值模板S threshold中的相邻网格通量差值绝对值的最小值
Figure 233683DEST_PATH_IMAGE085
，是为了减小S i , threshold的值，即减小模板内所允许的最大斜率，使得计算更为稳定。最后，判断
Figure 655437DEST_PATH_IMAGE040
对应的斜率值S i,max是否大于S i , threshold，当S i,max<S i , threshold时，则检测到目标分量g在网格单元i附近不存在间断；当最大斜率值S i,max
Figure 859017DEST_PATH_IMAGE039
S i , threshold斜率阈值时，则检测到g在网格单元i附近存在间断。
在进行通量间断检测的方法中将斜率阈值作为间断检测的标准，除了正数小量外，不包含任何经验参数，具有更好的普适性，同时，斜率阈值与所研究的网格单元所处的阈值模板内通量差值的绝对值正相关，对于不同的网格单元取值不同，具有更好的自适应性。
以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种可压缩流体的数值模拟方法，其特征在于，包括如下步骤：
步骤S10：读取计算网格、时
Figure DEST_PATH_IMAGE001
刻的物理变量分布、边界条件、试验时长；
步骤S20：在计算网格上，根据
Figure 266801DEST_PATH_IMAGE001
时刻的物理变量获得守恒变量
Figure 214028DEST_PATH_IMAGE002
和对流通量
Figure DEST_PATH_IMAGE003
、
Figure 60762DEST_PATH_IMAGE004
、
Figure DEST_PATH_IMAGE005
，其中，
Figure 769829DEST_PATH_IMAGE003
为笛卡尔坐标系下x方向上的对流通量、
Figure 967593DEST_PATH_IMAGE004
为笛卡尔坐标系下y方向上的对流通量；
Figure 718511DEST_PATH_IMAGE005
为笛卡尔坐标系下z方向上的对流通量；
步骤S30：分别计算对流通量
Figure 685330DEST_PATH_IMAGE003
、
Figure 362037DEST_PATH_IMAGE004
、
Figure 781517DEST_PATH_IMAGE005
在物理空间上的数值通量
Figure 398443DEST_PATH_IMAGE006
、
Figure DEST_PATH_IMAGE007
、
Figure 626293DEST_PATH_IMAGE008
；
步骤S40：将守恒变量
Figure 37683DEST_PATH_IMAGE002
和数值通量
Figure 646256DEST_PATH_IMAGE006
、
Figure 535715DEST_PATH_IMAGE007
、
Figure 742705DEST_PATH_IMAGE008
代入离散欧拉方程得到n+1时刻的守恒变量
Figure DEST_PATH_IMAGE009
及其对应的物理变量；
步骤S50：n从1~k遍历，其中k为
Figure 465942DEST_PATH_IMAGE010
时刻试验时长的序号，将守恒变量
Figure DEST_PATH_IMAGE011
及
Figure 302092DEST_PATH_IMAGE010
时刻的物理变量作为最终的试验结果。
2.如权利要求1所述的可压缩流体的数值模拟方法，其特征在于，步骤S30包括如下步骤：
步骤S301：构建对流通量
Figure 995241DEST_PATH_IMAGE003
、
Figure 260000DEST_PATH_IMAGE004
、
Figure 216455DEST_PATH_IMAGE005
在物理空间上的正负混合通量
Figure 596359DEST_PATH_IMAGE012
、
Figure DEST_PATH_IMAGE013
、
Figure 499724DEST_PATH_IMAGE014
；
步骤S302：对物理空间上的正负混合通量
Figure 415727DEST_PATH_IMAGE012
、
Figure 244881DEST_PATH_IMAGE013
、
Figure 613545DEST_PATH_IMAGE014
进行间断检测，若检测到正负混合通量
Figure 914076DEST_PATH_IMAGE012
、
Figure 684586DEST_PATH_IMAGE013
、
Figure 717264DEST_PATH_IMAGE014
均未存在间断，将对流通量
Figure 307646DEST_PATH_IMAGE003
、
Figure 848086DEST_PATH_IMAGE004
、
Figure 207523DEST_PATH_IMAGE005
进行分裂和逼近，得到物理空间上的数值通量
Figure DEST_PATH_IMAGE015
、
Figure 879944DEST_PATH_IMAGE016
、
Figure DEST_PATH_IMAGE017
；若检测到正负混合通量
Figure 880260DEST_PATH_IMAGE012
、
Figure 522593DEST_PATH_IMAGE013
、
Figure 267696DEST_PATH_IMAGE014
存在间断，将对流通量
Figure 376597DEST_PATH_IMAGE003
、
Figure 440106DEST_PATH_IMAGE004
、
Figure 886131DEST_PATH_IMAGE005
投影至特征空间，并分裂为正负通量
Figure 423422DEST_PATH_IMAGE018
、
Figure DEST_PATH_IMAGE019
、
Figure 405023DEST_PATH_IMAGE020
，并根据特征空间上的正负通量
Figure 457292DEST_PATH_IMAGE018
、
Figure 441429DEST_PATH_IMAGE019
、
Figure 895544DEST_PATH_IMAGE020
计算得到物理空间上的数值通量
Figure 611827DEST_PATH_IMAGE015
、
Figure 151393DEST_PATH_IMAGE016
、
Figure 375438DEST_PATH_IMAGE017
。
3.如权利要求2所述的可压缩流体的数值模拟方法，其特征在于，步骤S302中，根据特征空间上的正负通量
Figure 684060DEST_PATH_IMAGE018
、
Figure 899141DEST_PATH_IMAGE019
、
Figure 332527DEST_PATH_IMAGE020
，计算得到物理空间上的数值通量
Figure 897282DEST_PATH_IMAGE015
、
Figure 263672DEST_PATH_IMAGE016
、
Figure 118496DEST_PATH_IMAGE017
的步骤如下：
步骤S3021：对特征空间上的正负通量
Figure 632654DEST_PATH_IMAGE018
、
Figure 699967DEST_PATH_IMAGE019
、
Figure 481716DEST_PATH_IMAGE020
进行间断检测，若未检测到正负通量
Figure 710703DEST_PATH_IMAGE018
、
Figure 180998DEST_PATH_IMAGE019
、
Figure 848740DEST_PATH_IMAGE020
存在间断，则采用加权本质无振荡格式对应的高阶线性格式对正负通量
Figure 422679DEST_PATH_IMAGE018
、
Figure 884884DEST_PATH_IMAGE019
、
Figure 311317DEST_PATH_IMAGE020
进行逼近，得到特征空间上的数值通量；若检测到正负通量
Figure 517171DEST_PATH_IMAGE018
、
Figure 774977DEST_PATH_IMAGE019
、
Figure 844302DEST_PATH_IMAGE020
存在间断，则采用加权本质无振荡格式对正负通量
Figure 289189DEST_PATH_IMAGE018
、
Figure 564313DEST_PATH_IMAGE019
、
Figure 879888DEST_PATH_IMAGE020
进行逼近，得到特征空间上的数值通量；
步骤S3022：将特征空间上的数值通量投影回物理空间，得到物理空间上的数值通量
Figure 379834DEST_PATH_IMAGE015
、
Figure 780859DEST_PATH_IMAGE016
、
Figure 328515DEST_PATH_IMAGE017
。
4.如权利要求2所述的可压缩流体的数值模拟方法，其特征在于，步骤S302中，若检测到正负混合通量
Figure 233017DEST_PATH_IMAGE012
、
Figure 207927DEST_PATH_IMAGE013
、
Figure 657100DEST_PATH_IMAGE014
均未存在间断，将对流通量
Figure 8447DEST_PATH_IMAGE003
、
Figure 33035DEST_PATH_IMAGE004
、
Figure 178846DEST_PATH_IMAGE005
进行分裂和逼近，得到物理空间上的数值通量
Figure 85622DEST_PATH_IMAGE015
、
Figure 676878DEST_PATH_IMAGE016
、
Figure 883868DEST_PATH_IMAGE017
，获得物理空间上数值通量
Figure 872684DEST_PATH_IMAGE015
、
Figure 266756DEST_PATH_IMAGE016
、
Figure 225485DEST_PATH_IMAGE017
的方法如下：
采用通量分裂格式将对流通量
Figure 660883DEST_PATH_IMAGE003
、
Figure 679655DEST_PATH_IMAGE004
、
Figure 764285DEST_PATH_IMAGE005
分裂为正负通量
Figure DEST_PATH_IMAGE021
、
Figure 703203DEST_PATH_IMAGE022
、
Figure DEST_PATH_IMAGE023
；对正负通量
Figure 25732DEST_PATH_IMAGE021
、
Figure 215404DEST_PATH_IMAGE022
、
Figure 256173DEST_PATH_IMAGE023
采用加权本质无振荡格式对应的高阶线性格式进行逼近，获得物理空间中的数值通量
Figure 55239DEST_PATH_IMAGE015
、
Figure 560170DEST_PATH_IMAGE016
、
Figure 858427DEST_PATH_IMAGE017
。
5.如权利要求2所述的可压缩流体的数值模拟方法，其特征在于，步骤S302，将对流通量
Figure 448808DEST_PATH_IMAGE003
、
Figure 21872DEST_PATH_IMAGE004
、
Figure 348686DEST_PATH_IMAGE005
投影至特征空间，并分裂为正负通量
Figure 880162DEST_PATH_IMAGE018
、
Figure 957839DEST_PATH_IMAGE019
、
Figure 272277DEST_PATH_IMAGE020
的方法如下：
将对流通量
Figure 17379DEST_PATH_IMAGE003
、
Figure 454177DEST_PATH_IMAGE004
、
Figure 720948DEST_PATH_IMAGE005
由物理空间投影至特征空间，得到特征空间上对流通量
Figure 635814DEST_PATH_IMAGE024
、
Figure DEST_PATH_IMAGE025
、
Figure 564585DEST_PATH_IMAGE026
；采用通量分裂格式对特征空间上的对流通量
Figure 172284DEST_PATH_IMAGE024
、
Figure 958974DEST_PATH_IMAGE025
、
Figure 880794DEST_PATH_IMAGE026
分裂得到正负通量
Figure 334909DEST_PATH_IMAGE018
、
Figure 877624DEST_PATH_IMAGE019
、
Figure 823714DEST_PATH_IMAGE020
。
6.如权利要求2-3任意一项所述的可压缩流体的数值模拟方法，其特征在于，步骤S302中和步骤S3021中正负混合通量和正负通量的间断检测的方法如下：
获取对流通量
Figure 345962DEST_PATH_IMAGE028
、
Figure 857846DEST_PATH_IMAGE030
、
Figure 305883DEST_PATH_IMAGE032
在物理空间上的正负混合通量
Figure DEST_PATH_IMAGE033
、
Figure 473690DEST_PATH_IMAGE034
、
Figure DEST_PATH_IMAGE035
和特征空间上的正负通量
Figure 534050DEST_PATH_IMAGE036
、
Figure DEST_PATH_IMAGE037
、
Figure 602238DEST_PATH_IMAGE038
中的任意一个通量的分量，并将所述分量作为目标分量
Figure DEST_PATH_IMAGE040A
；将
Figure DEST_PATH_IMAGE040AA
在计算网格的网格单元
Figure DEST_PATH_IMAGE043
中心处的离散值记为
Figure DEST_PATH_IMAGE045
，遍历所有计算网格单元
Figure DEST_PATH_IMAGE047
，对目标分量
Figure DEST_PATH_IMAGE040AAA
在网格单元
Figure DEST_PATH_IMAGE047A
附近依次进行间断检测；
Figure DEST_PATH_IMAGE047AA
=1,2,...,N，N为网格单元的总个数；
获取目标分量g在网格单元i对应的加权无本质振荡格式大模板S whole中的最大斜率值S i,max；
获取目标分量g在网格单元i对应的加权本质无振荡格式阈值模板S threshold中的斜率阈值S i , threshold，其中加权无本质振荡格式阈值模板S threshold为加权无本质振荡格式大模板S whole的子模板中任一迎风子模板；
当S i,max<S i , threshold时，则检测到目标分量g在网格单元i附近不存在间断；当S i,max
Figure 725570DEST_PATH_IMAGE052
S i , threshold时，则检测到目标分量g在网格单元i附近存在间断。
7.如权利要求6所述的可压缩流体的数值模拟方法，其特征在于，获取目标分量g在网格单元i对应的加权无本质振荡格式大模板S whole中的最大斜率值S i,max的计算方法如下：
在网格单元i对应的加权无本质振荡格式的大模板S whole上对目标分量g采用加权本质无振荡格式对应的高阶线性格式进行逼近，得到逼近多项式函数
Figure DEST_PATH_IMAGE053
；
计算逼近多项式函数
Figure 911832DEST_PATH_IMAGE053
在大模板S whole上的斜率的绝对值的最大值S i,max。
8.如权利要求6所述的可压缩流体的数值模拟方法，其特征在于，获取目标分量g在网格单元i对应的加权本质无振荡格式阈值模板S threshold中的斜率阈值S i , threshold的计算方法如下：
选取网格单元i对应的加权本质无振荡格式中的阈值模板S threshold；
在加权本质无振荡格式阈值模板S threshold中计算基于网格最小可分辨波长的目标分量g的斜率阈值S i , threshold，其中，网格最小可分辨波长为网格单元长度的两倍。
9.如权利要求7所述的可压缩流体的数值模拟方法，其特征在于，计算逼近多项式函数
Figure 41462DEST_PATH_IMAGE053
在大模板S whole上的斜率的绝对值的最大值S i,max的公式如下：
Figure 958297DEST_PATH_IMAGE054
其中，
Figure DEST_PATH_IMAGE055
为逼近多项式函数
Figure 187284DEST_PATH_IMAGE053
在大模板S whole上的导数。
10.如权利要求8所述的可压缩流体的数值模拟方法，其特征在于，计算斜率阈值S i , threshold的公式如下：
Figure 923159DEST_PATH_IMAGE056
其中，
Figure DEST_PATH_IMAGE057
为网格单元长度，r和b为加权本质无振荡格式中阈值模板S threshold对应的整数，m为网格单元编号，g m+1 和g m 为网格单元的中心位置处目标分量g的离散值，
Figure 230381DEST_PATH_IMAGE058
为正数小量。
CN202111045330.9A 2021-09-07 2021-09-07 一种可压缩流体的数值模拟方法 Withdrawn - After Issue CN113486453B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111045330.9A CN113486453B (zh) 2021-09-07 2021-09-07 一种可压缩流体的数值模拟方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111045330.9A CN113486453B (zh) 2021-09-07 2021-09-07 一种可压缩流体的数值模拟方法

Publications (2)

Publication Number Publication Date
CN113486453A true CN113486453A (zh) 2021-10-08
CN113486453B CN113486453B (zh) 2021-11-19

Family

ID=77947242

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111045330.9A Withdrawn - After Issue CN113486453B (zh) 2021-09-07 2021-09-07 一种可压缩流体的数值模拟方法

Country Status (1)

Country Link
CN (1) CN113486453B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116384288A (zh) * 2023-06-05 2023-07-04 中国空气动力研究与发展中心计算空气动力研究所 一种可压缩流动高分辨率数值模拟方法、介质及设备

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104408291A (zh) * 2014-10-27 2015-03-11 西北工业大学 一种统一气体动理论格式的抑制振荡方法
CN105045954A (zh) * 2015-06-09 2015-11-11 北京交通大学 一种陡坎河床洪水冲刷演变模拟方法及系统
CN110781626A (zh) * 2019-10-31 2020-02-11 南京航空航天大学 有限差分多重分辨三角函数weno格式的模拟方法
CN111046615A (zh) * 2019-12-27 2020-04-21 中国人民解放军国防科技大学 一种黎曼解算器激波不稳定性抑制方法及系统
CN111859819A (zh) * 2020-06-16 2020-10-30 空气动力学国家重点实验室 一类高阶weno格式的构造方法
US20210012046A1 (en) * 2019-07-12 2021-01-14 Beihang University Meshless method for solid mechanics simulation, electronic device, and storage medium

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104408291A (zh) * 2014-10-27 2015-03-11 西北工业大学 一种统一气体动理论格式的抑制振荡方法
CN105045954A (zh) * 2015-06-09 2015-11-11 北京交通大学 一种陡坎河床洪水冲刷演变模拟方法及系统
US20210012046A1 (en) * 2019-07-12 2021-01-14 Beihang University Meshless method for solid mechanics simulation, electronic device, and storage medium
CN110781626A (zh) * 2019-10-31 2020-02-11 南京航空航天大学 有限差分多重分辨三角函数weno格式的模拟方法
CN111046615A (zh) * 2019-12-27 2020-04-21 中国人民解放军国防科技大学 一种黎曼解算器激波不稳定性抑制方法及系统
CN111859819A (zh) * 2020-06-16 2020-10-30 空气动力学国家重点实验室 一类高阶weno格式的构造方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
CHAI DELIN等: "An efficient modified WENO scheme based on the identification of inflection points", 《COMPUTERS & FLUIDS》 *
李冬冬: "激波与伴随蒸发和燃烧的铝液滴相互作用的数值模拟研究", 《中国优秀博硕士学位论文全文数据库(博士)基础科学辑》 *
王强等: "立楔诱导高超声速分离流动的被动控制研究", 《航空动力学报》 *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116384288A (zh) * 2023-06-05 2023-07-04 中国空气动力研究与发展中心计算空气动力研究所 一种可压缩流动高分辨率数值模拟方法、介质及设备
CN116384288B (zh) * 2023-06-05 2023-08-25 中国空气动力研究与发展中心计算空气动力研究所 一种可压缩流动高分辨率数值模拟方法、介质及设备

Also Published As

Publication number Publication date
CN113486453B (zh) 2021-11-19

Similar Documents

Publication Publication Date Title
Martínez-Tossas et al. The aerodynamics of the curled wake: a simplified model in view of flow control
US11168667B2 (en) Method and device for calculating power generation of wind farm
CN101858713B (zh) 有限元分析中的爆炸模拟
CN108763683B (zh) 一种三角函数框架下新weno格式构造方法
CN112100835B (zh) 一种适用于复杂流动的高效高精度翼型绕流数值模拟方法
JP2016530585A (ja) 技術的な系の出力量のモデルを算出する方法
CN109992889B (zh) 风电场模型的建立方法及系统、尾流值计算方法及系统
CN113486453B (zh) 一种可压缩流体的数值模拟方法
CN116029219B (zh) 一种飞行器气动热预测方法、装置、设备及存储介质
CN116738891B (zh) 一种增强飞行器流场模拟稳定性的lu-sgs改进方法
CN108509718B (zh) 一种基于质量守恒的远场尾流二维解析模型
CN115496006A (zh) 一种适用于高超声速飞行器的高精度数值模拟方法
CN112733473A (zh) 一种基于cfd的船舶自由横摇衰减数值模拟自动化智能化方法
CN112001109A (zh) 再生核粒子算法实现结构冲击动力学仿真方法
CN113864112A (zh) 风力发电机组的尾流流场的确定方法、装置及系统
CN106682292A (zh) 一种降维模拟退火算法的叶根结构优化方法
CN116341421B (zh) 高超声速流场数值模拟方法、系统、电子设备及存储介质
JP2017072922A (ja) 軸流ファンの解析方法，解析装置及び解析プログラム
Yang et al. A novel global optimization algorithm and its application to airfoil optimization
CN107122533A (zh) 一种基于efdc程序改进的水面热交换数值模拟方法
Jin Numerical simulation of wind turbine wakes based on actuator line method in NEK5000
CN114676522B (zh) 融合gan和迁移学习的气动形状优化设计方法及系统及设备
JP2021088974A (ja) 風車後流演算装置、及び風車後流演算方法
CN112182771B (zh) 基于数值模拟的数据处理方法、存储介质、电子装置
Wilkinson et al. Approximating urban wind interference

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
GR01 Patent grant
GR01 Patent grant
AV01 Patent right actively abandoned
AV01 Patent right actively abandoned
AV01 Patent right actively abandoned

Granted publication date: 20211119

Effective date of abandoning: 20211201

AV01 Patent right actively abandoned

Granted publication date: 20211119

Effective date of abandoning: 20211201