CN114676378B - 基于rad求解器的激波计算方法 - Google Patents

基于rad求解器的激波计算方法 Download PDF

Info

Publication number
CN114676378B
CN114676378B CN202111390993.4A CN202111390993A CN114676378B CN 114676378 B CN114676378 B CN 114676378B CN 202111390993 A CN202111390993 A CN 202111390993A CN 114676378 B CN114676378 B CN 114676378B
Authority
CN
China
Prior art keywords
rad
matrix
solver
constructed
shock wave
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.)
Active
Application number
CN202111390993.4A
Other languages
English (en)
Other versions
CN114676378A (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.)
Computational Aerodynamics Institute of China Aerodynamics Research and Development Center
Original Assignee
Computational 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 Computational Aerodynamics Institute of China Aerodynamics Research and Development Center filed Critical Computational Aerodynamics Institute of China Aerodynamics Research and Development Center
Priority to CN202111390993.4A priority Critical patent/CN114676378B/zh
Publication of CN114676378A publication Critical patent/CN114676378A/zh
Application granted granted Critical
Publication of CN114676378B publication Critical patent/CN114676378B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T3/00—Geometric image transformation in the plane of the image
    • G06T3/40—Scaling the whole image or part thereof
    • G06T3/4007—Interpolation-based scaling, e.g. bilinear interpolation
    • 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)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Data Mining & Analysis (AREA)
  • Computational Mathematics (AREA)
  • Operations Research (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Complex Calculations (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了基于RAD求解器的激波计算方法，包括步骤：S1，将激波计算系统中用到的特征矩阵通过混合反扩散矩阵来构造实现，并利用构造得到的特征矩阵来构造得到特征值，再利用构造得到的特征值来构造RAD通量，将构造得到的RAD通量记为RAD求解器；S2，利用步骤S1中得到的RAD求解器，在激波计算系统中计算激波。本发明应用到高阶加权紧致格式时，能够保持原有插值格式一致的高精度，并且格式的稳定性也很好，既能精确捕捉接触间断和激波，又能大幅提高对剪切层等精细结构的分辨率。

Description

基于RAD求解器的激波计算方法
技术领域
本发明涉及激波计算领域，更为具体的，涉及基于RAD求解器的激波计算方法。
背景技术
近似Riemann求解器常用在激波计算方法中，并广泛用于超声速流动的模拟计算。但是，在多维计算中使用近似Riemann求解器，可能会遇到激波计算的红玉现象（carbunclephenomenon）。
现有的Riemann算子在计算激波问题时，依然会出现数值不稳定的红玉现象，对于超声速问题的求解具有很大的局限性，且存在耗散大，不能识别剪切波和接触波，从而不能精准刻画流场的精细结构等问题。
发明内容
本发明的目的在于消除数值计算中的激波不稳定性现象，因此提供一种基于RAD求解器的激波计算方法，该方法采用混合方法对反扩散矩阵增加适当的耗散，再利用构造的特征矩阵来构造特征值，再利用构造的特征值来构造RAD通量等，因此本发明提出的新型RAD求解器具有自适应耗散的良好特性，从而能够保持激波计算稳定性。并且，本发明提供的RAD求解器克服了传统近似Riemann求解器计算激波时出现红玉现象的缺陷，既能精确捕捉接触间断和激波，又能大幅提高对剪切层等精细结构的分辨率，该求解器应用到高阶加权紧致格式时，能够保持原有插值格式一致的高精度，并且格式的稳定性也很好，克服了传统方案中应用近似Riemann求解器的缺陷，既能精确捕捉接触间断和激波，又能大幅提高对剪切层等精细结构的分辨率。
本发明的目的是通过以下方案实现的：
基于RAD求解器的激波计算方法，包括：
S1，将激波计算系统中用到的特征矩阵通过混合反扩散矩阵来构造实现，并利用构造得到的特征矩阵来构造得到特征值，再利用构造得到的特征值来构造RAD通量，将构造得到的RAD通量记为RAD求解器；
S2，利用步骤S1中得到的RAD求解器，在激波计算系统中计算激波。
进一步地，在步骤S1中，按照如下公式构造混合反扩散矩阵
Figure 579601DEST_PATH_IMAGE001
：
Figure 604375DEST_PATH_IMAGE003
其中，
Figure 796322DEST_PATH_IMAGE004
均为自适应混合因子，且取值范围均为[0,1]；
Figure 141853DEST_PATH_IMAGE005
均为反扩散系数，分别定义为：
Figure 495474DEST_PATH_IMAGE006
，
Figure 289303DEST_PATH_IMAGE007
，
Figure 489340DEST_PATH_IMAGE008
为耗散参数， 在强激波时
Figure 13863DEST_PATH_IMAGE009
，否则
Figure 570789DEST_PATH_IMAGE010
；diag表示对角矩阵；
在反扩散系数中，
Figure 838959DEST_PATH_IMAGE011
为正信号速度，
Figure 893503DEST_PATH_IMAGE012
为负信号速度，且
Figure 854506DEST_PATH_IMAGE013
，
Figure 627290DEST_PATH_IMAGE014
；
其中，
Figure 699151DEST_PATH_IMAGE015
为数值信号速度，且
Figure 873780DEST_PATH_IMAGE016
，
Figure 740105DEST_PATH_IMAGE017
，其中，q表示在各个方向的分速度u，v，w；
Figure 185DEST_PATH_IMAGE018
是声速的Roe平 均值，Max表示最大值，Min表示最小值。
进一步地，在步骤S1中，利用构造的反扩散矩阵
Figure 875737DEST_PATH_IMAGE019
按照如下公式构造特征矩阵D：
Figure 904873DEST_PATH_IMAGE020
其中，I为单位矩阵，“~”表示Roe平均，左特征矩阵
Figure 207679DEST_PATH_IMAGE021
，右特征矩阵
Figure 689476DEST_PATH_IMAGE022
，特征值
Figure 368719DEST_PATH_IMAGE023
分 别为：
Figure 252361DEST_PATH_IMAGE024
左特征矩阵
Figure 726068DEST_PATH_IMAGE025
中，其中
Figure 695161DEST_PATH_IMAGE026
是第k个维数的左特征向量；右特征矩阵
Figure 912515DEST_PATH_IMAGE027
中，
Figure 916243DEST_PATH_IMAGE028
是第k个维数的右特征向量；特征值矩阵
Figure 560851DEST_PATH_IMAGE023
中，
Figure 17241DEST_PATH_IMAGE029
是第k个维数的特征值，其中q表示在各个方 向的分速度u，v，w；
Figure 38286DEST_PATH_IMAGE030
是声速的Roe平均值。
进一步地，在步骤S1中，利用特征矩阵D按照如下公式构造特征值Q：
Figure 630942DEST_PATH_IMAGE031
Figure 446451DEST_PATH_IMAGE032
为特征值。
进一步地，在步骤S1中，利用特征值
Figure 655715DEST_PATH_IMAGE032
按照如下公式构造RAD通量
Figure 949293DEST_PATH_IMAGE033
：
Figure 927614DEST_PATH_IMAGE034
其中
Figure 648445DEST_PATH_IMAGE035
，
Figure 73567DEST_PATH_IMAGE036
，
Figure 170836DEST_PATH_IMAGE037
是第k个维数的波强度，
Figure 3663DEST_PATH_IMAGE038
是第k个维数 的右特征向量，
Figure 160975DEST_PATH_IMAGE039
为右守恒变量，
Figure 79252DEST_PATH_IMAGE040
为左守恒变量，F为流通量。
进一步地，自适应混合因子的取法包括：
Figure 980212DEST_PATH_IMAGE041
；
其中，
Figure 401966DEST_PATH_IMAGE042
，
定义
Figure 730179DEST_PATH_IMAGE043
，
定义
Figure 870174DEST_PATH_IMAGE044
，
定义
Figure 574824DEST_PATH_IMAGE045
；
其中p为压强，j为第j个网格点处的指标。
本发明的有益效果是：
本发明实施例中，针对现有计算方法在计算激波时的数值不稳定性问题，通过对反扩散矩阵增加适当的耗散，再利用构造的特征矩阵来构造特征值，再利用构造的特征值来构造RAD通量等，因此本发明提出的新型RAD求解器具有自适应耗散的良好特性，从而能够保持激波计算稳定性。并且，本发明提供的RAD求解器克服了传统近似Riemann求解器计算激波时出现红玉现象的缺陷，既能精确捕捉接触间断和激波，又能大幅提高对剪切层等精细结构的分辨率，该新型RAD求解器应用到高阶加权紧致格式时，能够保持与原有差分格式一致的高精度。并且格式的稳定性也很好，克服了传统近似Riemann求解器的缺陷，既能精确捕捉接触间断和激波，又能大幅提高对剪切层等精细结构的分辨率。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
图1为格式的色散特性比较图；
图2为格式的耗散特性比较图；
图3中（a）为WENO5格式求解Lax问题中采用不同Riemann求解器的计算方法的比较图；
图3中（b）为WENO5格式求解Lax问题中的局部放大图；
图4中（a）为HONCS5格式求解Lax问题中采用不同Riemann求解器的的计算方法比较图；
图4中（b）为HONCS5格式求解Lax问题中的局部放大图；
图5中（a）为一阶迎风格式求解双马赫反射问题中采用HLL求解器的计算方法的计算结果图；
图5中（b）为一阶迎风格式求解双马赫反射问题中采用HLLC求解器的计算方法的计算结果图；
图5中（c）为一阶迎风格式求解双马赫反射问题中利用本发明实施例的方法的计算结果图；
图6中（a）为WENO5格式求解双马赫反射问题中采用HLL求解器的计算方法的计算结果图；
图6中（b）为WENO5格式求解双马赫反射问题中采样HLLC求解器的计算方法的计算结果图；
图6中（c）为WENO5格式求解双马赫反射问题中采用本发明实施例的方法的计算结果图；
图7中（a）为HONCS5格式求解双马赫反射问题中采用本发明实施例的方法的计算结果图
图7中（b）为HONCS5格式求解双马赫反射问题中采用HLLC求解器的计算方法的计算结果图；
图7中（c）为HONCS5格式求解双马赫反射问题中采用本发明实施例的方法的计算结果图；
图8中（a）为WENO5格式求解前台阶问题中采用HLL求解器的计算方法的计算结果图；
图8中（b）为WENO5格式求解前台阶问题中采用HLLC求解器的计算方法的计算结果图；
图8中（c）为WENO5格式求解前台阶问题中利用本发明实施例的方法的计算结果图；
图9为HONCS5格式求解前台阶问题中利用本发明实施例的方法的计算结果图；
图10为本发明的步骤流程图。
具体实施方式
本说明书中所有实施例公开的所有特征，或隐含公开的所有方法或过程中的步骤，除了互相排斥的特征和/或步骤以外，均可以以任何方式组合和/或扩展、替换。本发明中RAD求解器是自定义术语，即Riemann solver with Anti-Diffusion，简称RAD求解器。
如图1~10所示，基于RAD求解器的激波计算方法，包括：
S1，将激波计算系统中用到的特征矩阵通过混合反扩散矩阵来构造实现，并利用构造得到的特征矩阵来构造得到特征值，再利用构造得到的特征值来构造RAD通量，将构造得到的RAD通量记为RAD求解器；
S2，利用步骤S1中得到的RAD求解器，在激波计算系统中计算激波。
进一步地，在步骤S1中，按照如下公式构造混合反扩散矩阵
Figure 382243DEST_PATH_IMAGE001
：
Figure 508648DEST_PATH_IMAGE003
其中，
Figure 751411DEST_PATH_IMAGE004
均为自适应混合因子，且取值范围均为[0,1]；
Figure 147757DEST_PATH_IMAGE005
均为反扩散系数，分别定义为：
Figure 817773DEST_PATH_IMAGE006
，
Figure 978813DEST_PATH_IMAGE007
，
Figure 495245DEST_PATH_IMAGE008
为耗散参数， 在强激波时
Figure 336162DEST_PATH_IMAGE009
，否则
Figure 938044DEST_PATH_IMAGE010
；diag表示对角矩阵；
在反扩散系数中，
Figure 522609DEST_PATH_IMAGE011
为正信号速度，
Figure 159127DEST_PATH_IMAGE012
为负信号速度，且
Figure 170946DEST_PATH_IMAGE013
，
Figure 994545DEST_PATH_IMAGE014
；
其中，
Figure 382801DEST_PATH_IMAGE046
为数值信号速度，且
Figure 608246DEST_PATH_IMAGE016
，
Figure 790966DEST_PATH_IMAGE017
，其中，q表示在各个方向的分速度u，v，w；
Figure 367441DEST_PATH_IMAGE018
是声速的Roe平 均值，Max表示最大值，Min表示最小值。
进一步地，在步骤S1中，利用构造的反扩散矩阵
Figure 559387DEST_PATH_IMAGE019
按照如下公式构造特征矩阵D：
Figure 645198DEST_PATH_IMAGE020
其中，I为单位矩阵，“~”表示Roe平均，左特征矩阵
Figure 264398DEST_PATH_IMAGE021
，右特征矩阵
Figure 62590DEST_PATH_IMAGE022
，特征值
Figure 792649DEST_PATH_IMAGE023
分 别为：
Figure 727107DEST_PATH_IMAGE024
左特征矩阵
Figure 517208DEST_PATH_IMAGE025
中，其中
Figure 802696DEST_PATH_IMAGE026
是第k个维数的左特征向量；右特征矩阵
Figure 70866DEST_PATH_IMAGE027
中，
Figure 125410DEST_PATH_IMAGE028
是第k个维数的右特征向量；特征值矩阵
Figure 86413DEST_PATH_IMAGE023
中，
Figure 859197DEST_PATH_IMAGE029
是第k个维数的特征值，其中q表示在各个方 向的分速度u，v，w；
Figure 931058DEST_PATH_IMAGE030
是声速的Roe平均值。
进一步地，在步骤S1中，利用特征矩阵D按照如下公式构造特征值Q：
Figure 105687DEST_PATH_IMAGE031
Figure 972012DEST_PATH_IMAGE032
为特征值。
进一步地，在步骤S1中，利用特征值
Figure 232092DEST_PATH_IMAGE032
按照如下公式构造RAD通量
Figure 842065DEST_PATH_IMAGE033
：
Figure 871201DEST_PATH_IMAGE034
其中
Figure 908427DEST_PATH_IMAGE035
，
Figure 655803DEST_PATH_IMAGE036
，
Figure 69467DEST_PATH_IMAGE037
是第k个维数的波强度，
Figure 953109DEST_PATH_IMAGE038
是第k个维数 的右特征向量，
Figure 692395DEST_PATH_IMAGE039
为右守恒变量，
Figure 661488DEST_PATH_IMAGE040
为左守恒变量，F为流通量。
进一步地，自适应混合因子的取法包括：
Figure 613264DEST_PATH_IMAGE041
；
其中，
Figure 616992DEST_PATH_IMAGE047
，
定义
Figure 996021DEST_PATH_IMAGE043
，
定义
Figure 452410DEST_PATH_IMAGE044
，
定义
Figure 207876DEST_PATH_IMAGE045
；
其中p为压强，j为第j个网格点处的指标。
在本发明的其他实施例中，在步骤S1中可具体包括如下步骤：
S11，按照如下公式构造混合反扩散矩阵
Figure 66111DEST_PATH_IMAGE048
：
Figure 881620DEST_PATH_IMAGE049
S12，利用构造的反扩散矩阵
Figure 108463DEST_PATH_IMAGE050
按照如下公式构造特征矩阵D：
Figure 402041DEST_PATH_IMAGE051
S13，利用特征矩阵D按照如下公式构造特征值Q：
Figure 114782DEST_PATH_IMAGE052
S14，利用特征值
Figure 101192DEST_PATH_IMAGE053
按照如下公式构造RAD通量
Figure 532174DEST_PATH_IMAGE054
：
Figure 629443DEST_PATH_IMAGE055
对本发明实施例的技术效果验证如下：
1. 精度验证
为了验证本发明方法的计算精度，本发明实施例中，选取较为典型的格式：时间离散采用三阶TVD Runge-Kutta格式，空间离散采用五阶HONCS插值格式进行实施。
本发明实施例中以一维Euler方程的对流密度波问题为验证标准，其精确解为：
Figure 196690DEST_PATH_IMAGE056
周期性边界条件的计算域为：
Figure 354002DEST_PATH_IMAGE057
。计算的最终时刻为
Figure 272280DEST_PATH_IMAGE058
，初始的计算网格 为15个点，CFL=0.5。随着网格点数增多，CFL数为：
Figure 173239DEST_PATH_IMAGE059
，其中下标 “1”表示上一个时间层的值，下标“2”表示下一个时间层的值，N表示网格点数，n表示空间格 式精度。利用本发明实施例计算得到的精度表如下所示。
表1五阶HONCS格式的数值精度验证
Figure 860573DEST_PATH_IMAGE060
由本发明的上述实施例可知，从计算的精度表发现：五阶HONCS插值格式都达到了设计精度，验证了本发明方法对于插值类型格式是满足高阶精度需求的。
2．频谱分析
本发明实施例中，采用ADR (approximate dispersion relation)方法来分析五阶WENO重构格式以及五阶HONCS插值格式的色散性和耗散性。其中修正波数的实部对应于格式的色散性（分辨率），而修正波数的虚部对应于格式的耗散性，计算结果如图1、图2所示。
从图1-2中可以看出，利用本发明实施例方法得到的五阶HONCS格式的分辨率高于五阶WENO格式分辨率，并且五阶HONCS格式的耗散远小于五阶WENO格式的耗散。
3. 数值实验与分析
在本发明实施例中，空间离散采用五阶WENO格式以及五阶HONCS插值格式，时间离散采用三阶TVD Runge-Kutta格式，采用本发明实施例计算了如下案例场景。
因为HLLEM求解器和Roe求解器的区别仅仅在数值信号速度的选择上，所以只给出HLLEM求解器的计算结果。
(1) Lax激波管问题
该问题的初始条件为：
Figure 188786DEST_PATH_IMAGE061
计算区域取为
Figure 328780DEST_PATH_IMAGE062
，计算网格采用100个点，初始间断位于x=0处，最终计算 时刻取为t=1.3。为了比较本发明实施例方法和传统的应用HLL、HLLC和HLLEM型Riemann求 解器的优劣，图3和图4分别给出了采用五阶WENO格式和五阶HONCS格式的计算结果。可以看 出，这类RAD型Riemann求解器的计算结果在激波、膨胀波及接触间断附近都能达到基本无 波动，基于五阶WENO格式以及五阶HONCS插值格式的本发明实施例方法的计算结果与精确 解符合的都很好。从计算的结果可以看出，本发明实施例方法的表现最好，而应用HLL型 Riemann算子的计算方案因耗散最大而导致计算结果最差。
(2) 激波双马赫反射问题
双马赫反射问题描述的是Mach数为10的强运动斜激波以与x轴方向呈60度角的方向入射，入射点在（1/6,0），计算区域取[0,4] X[0,1]。激波波前静止空气的密度为1.4，压力为1，波后按激波关系给定初始条件。下边界在[1/6,4]的区域给定壁面边界条件，其它边界区域按照激波运动所在的位置，分别给定波前或波后的值。初始物理参数为:
Figure 33431DEST_PATH_IMAGE063
本发明实施例中采用1920x480的均匀网格计算到无量纲时间t=0.2的时刻。
图5-7给出了一阶迎风格式、5阶WENO格式和5阶HONCS格式的计算结果，取密度为1.731到20.92共30条等值线作图。本发明实施例中给出一阶迎风格式计算结果的目的是为了更清晰的显示激波计算的红玉现象。
由于采用HLLEM算子时的方案都导致计算发散，所以只给出采用HLL、HLLC和RAD算子的计算方案的计算结果。可以看出，应用HLLC算子的方案会导致计算激波不稳定，对于一阶迎风格式尤其严重，这种红玉现象在WENO5格式和HONCS5格式的计算结果（右下角马赫杆处）中也比较明显。
利用本发明实施例方法的计算结果在激波附件都能达到基本无波动，同时对于剪切层的分辨率比HLL算子的分辨率更高。
(3) 前台阶问题
前台阶问题的初始条件为Mach数为3的强运动激波沿着带有台阶的风洞向右传播。台阶高度为0.2，前缘位置在0.6，计算域为[0,3] X[0,1]。初始物理参数为：
Figure 575271DEST_PATH_IMAGE064
左、右边界为开边界，其余边界取壁面反射条件。采用960x320的均匀网格计算截止到无量纲时间t=4的时刻。
图8-9给出了WENO5格式和HONCS5格式的计算结果，取密度为0.2568到6.607共60条等值线作图。
采用HLLEM算子时都导致计算发散。对于WENO5格式，采用HLL和HLLC算子时都导致数值震荡。对于HONCS5格式，采用HLL和HLLC算子时都导致计算发散，所以图中只给出RAD算子的计算结果。
可以看出，本发明实施例方法的计算结果在激波附近都能达到基本无波动，对剪切层的分辨率也非常高。
本发明未涉及部分均与现有技术相同或可采用现有技术加以实现。
上述技术方案只是本发明的一种实施方式，对于本领域内的技术人员而言，在本发明公开了应用方法和原理的基础上，很容易做出各种类型的改进或变形，而不仅限于本发明上述具体实施方式所描述的方法，因此前面描述的方式只是优选的，而并不具有限制性的意义。
除以上实例以外，本领域技术人员根据上述公开内容获得启示或利用相关领域的知识或技术进行改动获得其他实施例，各个实施例的特征可以互换或替换，本领域人员所进行的改动和变化不脱离本发明的精神和范围，则都应在本发明所附权利要求的保护范围内。

Claims (1)

1.基于RAD求解器的激波计算方法，其特征在于，包括：
S1，将激波计算系统中用到的特征矩阵通过混合反扩散矩阵来构造实现，并利用构造得到的特征矩阵来构造得到特征值，再利用构造得到的特征值来构造RAD通量，将构造得到的RAD通量记为RAD求解器；
在步骤S1中，按照如下公式构造混合反扩散矩阵
Figure FDA0004106355640000011
Figure FDA0004106355640000012
其中，
Figure FDA0004106355640000013
均为自适应混合因子，且取值范围均为[0,1]；
Figure FDA00041063556400000115
δ1、δ4均为反扩散系数，分别定义为：
Figure FDA0004106355640000015
ε0为耗散参数，在强激波时ε0＝0.1，否则ε0＝0；diag表示对角矩阵；
在反扩散系数中，S+为正信号速度，S-为负信号速度，且S-＝Min(SL，0)，S+＝Max(SR，0)；
其中，SL、SR为数值信号速度，且
Figure FDA0004106355640000016
Figure FDA0004106355640000017
其中，q表示在各个方向的分速度u，v，w，/>
Figure FDA0004106355640000018
表示在各个方向的分速度u，v，w的Roe平均值；/>
Figure FDA0004106355640000019
是声速的Roe平均值，Max表示最大值，Min表示最小值；
在步骤S1中，利用构造的反扩散矩阵
Figure FDA00041063556400000110
按照如下公式构造特征矩阵D：
Figure FDA00041063556400000111
其中，I为单位矩阵，“～”表示Roe平均，左特征矩阵
Figure FDA00041063556400000112
右特征矩阵/>
Figure FDA00041063556400000113
特征值/>
Figure FDA00041063556400000114
分别为：
Figure FDA0004106355640000021
左特征矩阵
Figure FDA0004106355640000022
中，其中/>
Figure FDA0004106355640000023
是第k个维数的左特征向量；右特征矩阵/>
Figure FDA0004106355640000024
中，/>
Figure FDA0004106355640000025
是第k个维数的右特征向量；特征值矩阵/>
Figure FDA0004106355640000026
中，/>
Figure FDA0004106355640000027
是第k个维数的特征值，其中q表示在各个方向的分速度u，v，w；/>
Figure FDA0004106355640000028
是声速的Roe平均值；
在步骤S1中，利用特征矩阵D按照如下公式构造特征值Q：
Figure FDA0004106355640000029
Qk为特征值；
在步骤S1中，利用特征值Qk按照如下公式构造RAD通量FRAD：
Figure FDA00041063556400000210
其中FR＝F(UR)，FL＝F(UL)，
Figure FDA00041063556400000211
是第k个维数的波强度，/>
Figure FDA00041063556400000212
是第k个维数的右特征向量，UR为右守恒变量，UL为左守恒变量，F为流通量；
S2，利用步骤S1中得到的RAD求解器，在激波计算系统中计算激波。
CN202111390993.4A 2021-11-23 2021-11-23 基于rad求解器的激波计算方法 Active CN114676378B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111390993.4A CN114676378B (zh) 2021-11-23 2021-11-23 基于rad求解器的激波计算方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111390993.4A CN114676378B (zh) 2021-11-23 2021-11-23 基于rad求解器的激波计算方法

Publications (2)

Publication Number Publication Date
CN114676378A CN114676378A (zh) 2022-06-28
CN114676378B true CN114676378B (zh) 2023-05-26

Family

ID=82070171

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111390993.4A Active CN114676378B (zh) 2021-11-23 2021-11-23 基于rad求解器的激波计算方法

Country Status (1)

Country Link
CN (1) CN114676378B (zh)

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5322242A (en) * 1991-07-08 1994-06-21 Tracy Richard R High efficiency, supersonic aircraft
CN107679319A (zh) * 2017-09-29 2018-02-09 北京航空航天大学 一种叶轮机通流模型中周向脉动应力项的代数建模方法
CN110516310A (zh) * 2019-07-31 2019-11-29 中国空气动力研究与发展中心 旋转爆震反压的非定常数值模拟方法
CN112100835A (zh) * 2020-09-06 2020-12-18 西北工业大学 一种适用于复杂流动的高效高精度数值模拟方法
CN112214869A (zh) * 2020-09-03 2021-01-12 空气动力学国家重点实验室 一种求解欧拉方程的改进型高阶非线性空间离散方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5322242A (en) * 1991-07-08 1994-06-21 Tracy Richard R High efficiency, supersonic aircraft
CN107679319A (zh) * 2017-09-29 2018-02-09 北京航空航天大学 一种叶轮机通流模型中周向脉动应力项的代数建模方法
CN110516310A (zh) * 2019-07-31 2019-11-29 中国空气动力研究与发展中心 旋转爆震反压的非定常数值模拟方法
CN112214869A (zh) * 2020-09-03 2021-01-12 空气动力学国家重点实验室 一种求解欧拉方程的改进型高阶非线性空间离散方法
CN112100835A (zh) * 2020-09-06 2020-12-18 西北工业大学 一种适用于复杂流动的高效高精度数值模拟方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
闫文辉 等.应用GAO-YONG湍流模式数值模拟三维激波/湍流边界层干扰.《航空动力学报》.2009,第24卷(第10期),第2193-2200页. *

Also Published As

Publication number Publication date
CN114676378A (zh) 2022-06-28

Similar Documents

Publication Publication Date Title
CN106487358B (zh) 一种机动目标转弯跟踪方法
CN105301666B (zh) 一种航磁干扰补偿系数的自适应调整方法
US9160280B1 (en) Memory polynomial based digital predistorter
CN109307855A (zh) 基于网格误差模型的无网格稀疏近似最小方差doa估计方法
CN110632555B (zh) 一种基于矩阵特征值扰动的tdoa直接定位方法
CN103500455B (zh) 一种基于无偏有限冲击响应滤波器（ufir）的改进机动目标跟踪方法
Kay et al. Improvement of TDOA position fixing using the likelihood curvature
CN111220222B (zh) 一种超声波燃气表流量的测定算法
CN114676378B (zh) 基于rad求解器的激波计算方法
CN112665820B (zh) 基于变量差及相对位移的r型网格自适应移动方法及设备
CN109582915B (zh) 应用于纯方位跟踪的改进非线性可观测度自适应滤波方法
CN105117537A (zh) 一种基于权值比较的粒子滤波系统重采样方法
CN109375160B (zh) 纯方位无源定位中一种测角误差估计方法
KR102453305B1 (ko) 비대칭도를 이용한 시간이력풍하중 생성 방법
CN111487440A (zh) 一种五孔探针的标定方法
CN114465851B (zh) 优化核宽最大箕舌线准则的簇稀疏水声信道估计方法
CN116559579A (zh) 基于改进的VMD和Teager能量算子故障定位方法
CN114637956B (zh) 一种基于双卡尔曼滤波器实现目标位置预测的方法
JP7321612B2 (ja) ゲイン位相誤差が存在する場合、スパース再構成に基づく到来方向推定方法
CN114791994A (zh) 一种引入法向量优化的ransac点云平面拟合方法
CN113452350A (zh) 一种变步长块稀疏仿射投影自适应滤波器
CN110602635A (zh) 一种室内地图匹配增强定位方法、设备及存储设备
CN105717527A (zh) 一种利用边坡变形数据快速确定监测点移动轨迹的方法
Yu et al. Determination of vocal-tract shapes from formant frequencies based on perturbation theory and interpolation method
CN112540342B (zh) 一种移动信号源定位方法

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