基于RAD求解器的激波计算方法与流程

基于rad求解器的激波计算方法
技术领域
1.本发明涉及激波计算领域，更为具体的，涉及基于rad求解器的激波计算方法。


背景技术：

2.近似riemann求解器常用在激波计算方法中，并广泛用于超声速流动的模拟计算。但是，在多维计算中使用近似riemann求解器，可能会遇到激波计算的红玉现象（carbuncle phenomenon）。
3.现有的riemann算子在计算激波问题时，依然会出现数值不稳定的红玉现象，对于超声速问题的求解具有很大的局限性，且存在耗散大，不能识别剪切波和接触波，从而不能精准刻画流场的精细结构等问题。


技术实现要素：

4.本发明的目的在于消除数值计算中的激波不稳定性现象，因此提供一种基于rad求解器的激波计算方法，该方法采用混合方法对反扩散矩阵增加适当的耗散，再利用构造的特征矩阵来构造特征值，再利用构造的特征值来构造rad通量等，因此本发明提出的新型rad求解器具有自适应耗散的良好特性，从而能够保持激波计算稳定性。并且，本发明提供的rad求解器克服了传统近似riemann求解器计算激波时出现红玉现象的缺陷，既能精确捕捉接触间断和激波，又能大幅提高对剪切层等精细结构的分辨率，该求解器应用到高阶加权紧致格式时，能够保持原有插值格式一致的高精度，并且格式的稳定性也很好，克服了传统方案中应用近似riemann求解器的缺陷，既能精确捕捉接触间断和激波，又能大幅提高对剪切层等精细结构的分辨率。
5.本发明的目的是通过以下方案实现的：基于rad求解器的激波计算方法，包括：s1，将激波计算系统中用到的特征矩阵通过混合反扩散矩阵来构造实现，并利用构造得到的特征矩阵来构造得到特征值，再利用构造得到的特征值来构造rad通量，将构造得到的rad通量记为rad求解器；s2，利用步骤s1中得到的rad求解器，在激波计算系统中计算激波。
6.进一步地，在步骤s1中，按照如下公式构造混合反扩散矩阵：其中， 均为自适应混合因子，且取值范围均为[0,1]；均为反扩散系数，分别定义为：，，为耗散参数，
在强激波时，否则；diag表示对角矩阵；在反扩散系数中，为正信号速度，为负信号速度，且，；其中， 为数值信号速度，且，，其中，q表示在各个方向的分速度u，v，w；是声速的roe平均值，max表示最大值，min表示最小值。
[0007]
进一步地，在步骤s1中，利用构造的反扩散矩阵按照如下公式构造特征矩阵d：其中，i为单位矩阵，“~”表示roe平均，左特征矩阵，右特征矩阵，特征值分别为：左特征矩阵中，其中是第k个维数的左特征向量；右特征矩阵中，是第k个维数的右特征向量；特征值矩阵中，是第k个维数的特征值，其中q表示在各个方向的分速度u，v，w；是声速的roe平均值。
[0008]
进一步地，在步骤s1中，利用特征矩阵d按照如下公式构造特征值q：进一步地，在步骤s1中，利用特征矩阵d按照如下公式构造特征值q：为特征值。
[0009]
进一步地，在步骤s1中，利用特征值按照如下公式构造rad通量：其中，，是第k个维数的波强度，是第k个维数
的右特征向量，为右守恒变量，为左守恒变量，f为流通量。
[0010]
进一步地，自适应混合因子的取法包括：；其中，，定义，定义，定义；其中p为压强，j为第j个网格点处的指标。
[0011]
本发明的有益效果是：本发明实施例中，针对现有计算方法在计算激波时的数值不稳定性问题，通过对反扩散矩阵增加适当的耗散，再利用构造的特征矩阵来构造特征值，再利用构造的特征值来构造rad通量等，因此本发明提出的新型rad求解器具有自适应耗散的良好特性，从而能够保持激波计算稳定性。并且，本发明提供的rad求解器克服了传统近似riemann求解器计算激波时出现红玉现象的缺陷，既能精确捕捉接触间断和激波，又能大幅提高对剪切层等精细结构的分辨率，该新型rad求解器应用到高阶加权紧致格式时，能够保持与原有差分格式一致的高精度。并且格式的稳定性也很好，克服了传统近似riemann求解器的缺陷，既能精确捕捉接触间断和激波，又能大幅提高对剪切层等精细结构的分辨率。
附图说明
[0012]
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
[0013]
图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为本发明的步骤流程图。
具体实施方式
[0014]
本说明书中所有实施例公开的所有特征，或隐含公开的所有方法或过程中的步骤，除了互相排斥的特征和/或步骤以外，均可以以任何方式组合和/或扩展、替换。本发明中rad求解器是自定义术语，即riemann solver with anti-diffusion，简称rad求解器。
[0015]
如图1~10所示，基于rad求解器的激波计算方法，包括：s1，将激波计算系统中用到的特征矩阵通过混合反扩散矩阵来构造实现，并利用构造得到的特征矩阵来构造得到特征值，再利用构造得到的特征值来构造rad通量，将构造得到的rad通量记为rad求解器；s2，利用步骤s1中得到的rad求解器，在激波计算系统中计算激波。
[0016]
进一步地，在步骤s1中，按照如下公式构造混合反扩散矩阵 ：
其中， 均为自适应混合因子，且取值范围均为[0,1]；均为反扩散系数，分别定义为：，，为耗散参数，在强激波时，否则；diag表示对角矩阵；在反扩散系数中，为正信号速度，为负信号速度，且，；其中， 为数值信号速度，且，，其中，q表示在各个方向的分速度u，v，w；是声速的roe平均值，max表示最大值，min表示最小值。
[0017]
进一步地，在步骤s1中，利用构造的反扩散矩阵按照如下公式构造特征矩阵d：其中，i为单位矩阵，“~”表示roe平均，左特征矩阵，右特征矩阵，特征值分别为：左特征矩阵中，其中是第k个维数的左特征向量；右特征矩阵中，是第k个维数的右特征向量；特征值矩阵中，是第k个维数的特征值，其中q表示在各个方向的分速度u，v，w；是声速的roe平均值。
[0018]
进一步地，在步骤s1中，利用特征矩阵d按照如下公式构造特征值q：进一步地，在步骤s1中，利用特征矩阵d按照如下公式构造特征值q：为特征值。
[0019]
进一步地，在步骤s1中，利用特征值按照如下公式构造rad通量：其中，，是第k个维数的波强度，是第k个维数的右特征向量，为右守恒变量，为左守恒变量，f为流通量。
[0020]
进一步地，自适应混合因子的取法包括：；其中，，定义，定义，定义；其中p为压强，j为第j个网格点处的指标。
[0021]
在本发明的其他实施例中，在步骤s1中可具体包括如下步骤：s11，按照如下公式构造混合反扩散矩阵：s12，利用构造的反扩散矩阵按照如下公式构造特征矩阵d：s13，利用特征矩阵d按照如下公式构造特征值q：s14，利用特征值按照如下公式构造rad通量：
对本发明实施例的技术效果验证如下：1. 精度验证为了验证本发明方法的计算精度，本发明实施例中，选取较为典型的格式：时间离散采用三阶tvd runge-kutta格式，空间离散采用五阶honcs插值格式进行实施。
[0022]
本发明实施例中以一维euler方程的对流密度波问题为验证标准，其精确解为：周期性边界条件的计算域为：。计算的最终时刻为，初始的计算网格为15个点，cfl=0.5。随着网格点数增多，cfl数为：，其中下标“1”表示上一个时间层的值，下标“2”表示下一个时间层的值，n表示网格点数，n表示空间格式精度。利用本发明实施例计算得到的精度表如下所示。
[0023]
表1五阶honcs格式的数值精度验证由本发明的上述实施例可知，从计算的精度表发现：五阶honcs插值格式都达到了设计精度，验证了本发明方法对于插值类型格式是满足高阶精度需求的。
[0024]
2．频谱分析本发明实施例中，采用adr (approximate dispersion relation)方法来分析五阶weno重构格式以及五阶honcs插值格式的色散性和耗散性。其中修正波数的实部对应于格式的色散性（分辨率），而修正波数的虚部对应于格式的耗散性，计算结果如图1、图2所示。
[0025]
从图1-2中可以看出，利用本发明实施例方法得到的五阶honcs格式的分辨率高于五阶weno格式分辨率，并且五阶honcs格式的耗散远小于五阶weno格式的耗散。
[0026]
3. 数值实验与分析在本发明实施例中，空间离散采用五阶weno格式以及五阶honcs插值格式，时间离散采用三阶tvd runge-kutta格式，采用本发明实施例计算了如下案例场景。
[0027]
因为hllem求解器和roe求解器的区别仅仅在数值信号速度的选择上，所以只给出hllem求解器的计算结果。
[0028]
(1) lax激波管问题该问题的初始条件为：计算区域取为，计算网格采用100个点，初始间断位于x=0处，最终计算时刻取为t=1.3。为了比较本发明实施例方法和传统的应用hll、hllc和hllem型riemann求解器的优劣，图3和图4分别给出了采用五阶weno格式和五阶honcs格式的计算结果。可以看出，这类rad型riemann求解器的计算结果在激波、膨胀波及接触间断附近都能达到基本无波动，基于五阶weno格式以及五阶honcs插值格式的本发明实施例方法的计算结果与精确解符合的都很好。从计算的结果可以看出，本发明实施例方法的表现最好，而应用hll型riemann算子的计算方案因耗散最大而导致计算结果最差。
[0029]
(2) 激波双马赫反射问题双马赫反射问题描述的是mach数为10的强运动斜激波以与x轴方向呈60度角的方向入射，入射点在（1/6,0），计算区域取[0,4] x[0,1]。激波波前静止空气的密度为1.4，压力为1，波后按激波关系给定初始条件。下边界在[1/6,4]的区域给定壁面边界条件，其它边界区域按照激波运动所在的位置，分别给定波前或波后的值。初始物理参数为:本发明实施例中采用1920x480的均匀网格计算到无量纲时间t=0.2的时刻。
[0030]
图5-7给出了一阶迎风格式、5阶weno格式和5阶honcs格式的计算结果，取密度为1.731到20.92共30条等值线作图。本发明实施例中给出一阶迎风格式计算结果的目的是为了更清晰的显示激波计算的红玉现象。
[0031]
由于采用hllem算子时的方案都导致计算发散，所以只给出采用hll、hllc和rad算子的计算方案的计算结果。可以看出，应用hllc算子的方案会导致计算激波不稳定，对于一阶迎风格式尤其严重，这种红玉现象在weno5格式和honcs5格式的计算结果（右下角马赫杆处）中也比较明显。
[0032]
利用本发明实施例方法的计算结果在激波附件都能达到基本无波动，同时对于剪切层的分辨率比hll算子的分辨率更高。
[0033]
(3) 前台阶问题前台阶问题的初始条件为mach数为3的强运动激波沿着带有台阶的风洞向右传播。台阶高度为0.2，前缘位置在0.6，计算域为[0,3] x[0,1]。初始物理参数为：左、右边界为开边界，其余边界取壁面反射条件。采用960x320的均匀网格计算截止到无量纲时间t=4的时刻。
[0034]
图8-9给出了weno5格式和honcs5格式的计算结果，取密度为0.2568到6.607共60条等值线作图。
[0035]
采用hllem算子时都导致计算发散。对于weno5格式，采用hll和hllc算子时都导致数值震荡。对于honcs5格式，采用hll和hllc算子时都导致计算发散，所以图中只给出rad算子的计算结果。
[0036]
可以看出，本发明实施例方法的计算结果在激波附近都能达到基本无波动，对剪切层的分辨率也非常高。
[0037]
本发明未涉及部分均与现有技术相同或可采用现有技术加以实现。
[0038]
上述技术方案只是本发明的一种实施方式，对于本领域内的技术人员而言，在本发明公开了应用方法和原理的基础上，很容易做出各种类型的改进或变形，而不仅限于本发明上述具体实施方式所描述的方法，因此前面描述的方式只是优选的，而并不具有限制性的意义。
[0039]
除以上实例以外，本领域技术人员根据上述公开内容获得启示或利用相关领域的知识或技术进行改动获得其他实施例，各个实施例的特征可以互换或替换，本领域人员所进行的改动和变化不脱离本发明的精神和范围，则都应在本发明所附权利要求的保护范围内。