一种可压/不可压混合流动的高效扰动域推进方法

1.本公开涉及计算流体力学领域技术领域，尤其涉及一种可压/不可压混合流动的高效扰动域推进方法。


背景技术：

2.在航空航天领域，许多重要的问题都是可压缩/不可压缩的混合流动，而直接采用可压缩计算方法求解可压/不可压混合流动时，特征值的差异导致了严重的病态条件数，造成了效率降低和计算困难。传统的低速流动通常基于压力耦合方程组的半隐式方法进行扩展，这种扩展在计算低速问题时效率更高，但分离式求解压力和速度守恒性与密度基相比较差，且难以直接推广现有的激波捕捉格式。而以时间推进法和激波捕捉格式为特点的可压缩流动求解方法已经得到了长足的发展，但在直接求解低速流动时，由于声速与流速的过大差异，系统刚性过大，收敛速度过慢。预处理方法的出现使得应用可压缩方程求解可压/不可压混合流场成为了一条最优途径。
3.然而，预处理方法虽然提高了可压缩方程在求解可压/不可压流场的收敛速度，但在工程实际的应用中，对其计算效率提出了更高的要求。扰动域更新方法是针对可压缩流动求解提出的一种新型数值模拟加速方法，有效提高了数值模拟的计算效率。但传统的扰动域更新方法无法应用于预处理方法对可压/不可压流动的求解，其不足之处在于以下两点：
4.一是，预处理方法改变了对流项雅克比矩阵的特征值，而这特征值即为传统扰动域更新方法中作为扰动传播判据的当地速度与当地声速，这一改变影响了扰动域更新法中的诸多判据；二是，可压/不可压混合流动往往绝大多数区域为低马赫数流动，密度和温度变化较小，将总温或者总焓作为粘性分区判据会导致判断失效。
5.因此，目前缺少一种有效提高可压/不可压流动的计算效率，服务于相关领域的数值模拟计算的方法。


技术实现要素：

6.有鉴于此，本公开提供了一种可压/不可压混合流动的高效扰动域推进方法，以解决现有技术在不影响计算域有限性带来的误差的同时，大幅提高计算效率的技术问题。
7.本公开提供了一种可压/不可压混合流动的高效扰动域推进方法，包括：
8.s1数据读入与流场初始化；
9.s2根据所述流场初始化的方式，建立初始动态计算域，其中，动态计算域包括对流动态域和粘性动态域；
10.s3在所述动态计算域内，求解可压缩纳维-斯托克斯方程当前时间步的残差，并进行时间积分迭代求解网格参数；
11.s4基于s3得到的残差和网格参数，对所述动态计算域进行更新；
12.s5判断全局残差是否小于给定收敛残差阈值，若是，则计算完成并执行s6，若否，
则返回s3；
13.s6输出计算所得网格参数的值。
14.进一步地，所述s1，包括：
15.读入计算网格、边界条件、网格参数和计算设置；
16.基于所述网格参数，根据来流条件或给定流场进行所述流场初始化，确定计算前还未定义边界和域的物理量初值以及所选物理模型中公式参数，其中，物理模型中公式参数包括气体常数，粘性系数，湍流模型系数。
17.进一步地，所述s2，包括：
18.当根据所述来流条件初始化时，所有网格参数设定为来流参数，依据壁面条件建立所述对流动态域，同时仅取紧邻壁面的1层网格单元作为初始粘性动态域；
19.当根据给定流场初始化时，所有网格参数按照给定流场设置，依据所述给定流场与来流的差异建立所述粘性动态域，与来流流动特性不符的网格单元作为初始对流动态域，并选取对流动态域中的粘性主导网格单元作为粘性动态域。
20.进一步地，所述s3，包括：
21.基于流场初始化给定的所有网格参数，求解可压缩纳维-斯托克斯方程，并选择预处理后的roe格式，计算对流通量，获得当前时间步的残差；
22.基于所述残差结果，进行时间积分迭代求解，更新获得当前时间步的网格参数。
23.进一步地，所述s4，包括：
24.基于当前时间步的网格参数，根据预处理后的速度和声速，判断对流动态域是否需要增大，再结合残差的收敛情况，判断对流动态域是否需要减少；
25.判断当前网格单元是否受粘性效应主导，若受粘性单元主导，则将所述当前网格单元所有紧邻网格单元增加至粘性动态域，若不受粘性单元主导，则将当前网格单元从粘性动态域中删除。
26.进一步地，所述s4中所述对流动态域的增大与减小，具体包括：
27.通过衡量s3中计算得到的网格参数更新量，判断当前网格单元是否受到无粘扰动，若受到无粘扰动，确定会受影响的紧邻网格单元，对于处于亚声速流动的单元，所有紧邻网格单元都会受到扰动的影响，将所有紧邻网格单元加入到对流动态域中；对于处于超声速流动的单元，则根据对流通量的特征值确定扰动的传播方向，将传播方向下游的网格单元加入对流动态域中；
28.判断当前网格单元是否已收敛，若收敛，则判断所述当前网格单元是否位于所述对流动态域的最上游，若已位于最上游，则判断所述对流动态域中其他网格单元是否不再对所述当前网格单元产生影响，若其他所述网格单元不再对所述当前网格单元产生影响，则从所述对流动态域中删除当前网格单元。
29.进一步地，所述预处理后的速度和声速是基于当地速度和声速获得的，具体包括：
30.基于所述当地速度u与声速c获得预处理后的速度和声速；
31.基于所述预处理后速度和声速，获得预处理后的特征值，
32.其中，预处理后的速度和声速的计算式分别如下：
33.λ＝[u,u,u,u
′
c
′
,u
′‑c′
]
[0034][0035]ur2
＝min(c2,max(|u|2,k|u|2))
[0036]
其中，λ为计算对流通量时用到的特征值，u为当地速度，u
′
为预处理后的速度，c为声速，c
′
为预处理声速，α、β都是求解预处理特征值所用的参数，ur为预处理参数，由预处理方法给出，ρ,ρ
p
,ρ
t
分别为密度、密度对压力的导数、密度对温度的导数，由流场参数计算得到，c
p
为定压比热容常数，k是一个小值(10-5
)。
[0037]
进一步地，所述s4中所述粘性动态域的增大与减小，具体包括：
[0038]
将总焓变化量的阈值除以马赫数平方作为粘性判断依据，进行所述粘性动态域的增大与缩小，具体方式如下：
[0039]
若|1-h
0,∞
/h0|/ma2≥ε
d,v
，则将当前网格单元的所有紧邻网格单元增加至粘性动态域；
[0040]
若|1-h
0,∞
/h0|/ma2《ε
d,v
，则将当前网格单元从粘性动态域中删除；
[0041]
其中，h
0,∞
为来流总焓，h0为当前网格单元总焓，ma为当前网格单元马赫数，ε
d,v
为给定的粘性域判断阈值。
[0042]
进一步地，所述网格参数包括：
[0043]
密度、速度、压力和温度。
[0044]
本公开与现有技术相比存在的有益效果是：
[0045]
1.本发明基于预处理后的当地声速与当地速度，建立可压/不可压扰动域更新方法的扰动传播计算公式，从而建立可压/不可压扰动域传播和影响的判据，解决了雅克比矩阵特征值改变和低马赫数区域判据不适用的问题，收缩计算域；
[0046]
2.提出了适用于预处理求解可压/不可压流场的粘性主导判断依据，将总焓变化量的阈值除以马赫数平方作为粘性分区的判断阈值作为粘性分区的判断阈值，避免了采用总焓变化量作为判据时粘性分区对马赫数的依赖性；
[0047]
3.避免了大量无效计算，大幅提高计算效率。
附图说明
[0048]
为了更清楚地说明本公开中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其它的附图。
[0049]
图1为本发明提供的一种可压/不可压混合流动的高效扰动域推进方法的示意图；
[0050]
图2是本发明提供的对流动态域与粘性主导区域占比大小的示意图；
[0051]
图3a是本发明实施例提供的未采用可压/不可压流动的扰动域更新方法计算结果的对比示意图；图3b是本发明实施例提供的采用可压/不可压流动的扰动域更新方法计算结果的对比示意图；
[0052]
图4是本发明实施例提供的是否采用可压/不可压流动的扰动域更新方法收敛情况的对比示意图。
具体实施方式
[0053]
以下描述中，为了说明而不是为了限定，提出了诸如特定系统结构、技术之类的具体细节，以便透彻理解本公开实施例。然而，本领域的技术人员应当清楚，在没有这些具体细节的其它实施例中也可以实现本公开。在其它情况中，省略对众所周知的系统、装置、电路以及方法的详细说明，以免不必要的细节妨碍本公开的描述。
[0054]
下面将结合附图详细说明根据本公开的一种可压/不可压混合流动的高效扰动域推进方法。
[0055]
图1为本发明提供的一种可压/不可压混合流动的高效扰动域推进方法的示意图。
[0056]
如图1所示，该高效扰动域更新方法包括：
[0057]
s1，数据读入与流场初始化。
[0058]
s1包括：
[0059]
s11，读入计算网格、边界条件、计算设置。
[0060]
s12，基于网格参数，根据来流条件或给定流场进行流场初始化，确定计算前还未定义边界和域的物理量初值以及所选物理模型中公式参数。
[0061]
基于网格参数，进行流场初始化。流场初始化能够确定计算前还未定义边界和域的物理参数初值以及所选物理模型中公式参数，其中，物理模型中公式参数包括气体常数，粘性系数，湍流模型系数。
[0062]
s2，根据流场初始化的方式，建立初始动态计算域，其中，动态计算域包括对流动态域和粘性动态域。
[0063]
s21，当根据来流条件初始化时，所有网格参数设定为来流参数，依据壁面条件建立对流动态域，同时仅取紧邻壁面的1层网格单元作为初始粘性动态域。
[0064]
当根据来流条件初始化时，将所有网格参数设定为来流参数，依据壁面条件建立对流动态域，通常取与壁面边界相邻的10层网格单元作为对流动态域的初始网格单元，粘性动态域取与壁面边界相邻的第1层网格单元作为初始粘性动态域。
[0065]
s22，当根据给定流场初始化时，所有网格参数按照给定流场设置，依据给定流场与来流的差异建立粘性动态域，与来流流动特性不符的网格单元作为初始对流动态域，并选取对流动态域中的粘性主导网格单元作为粘性动态域。
[0066]
s3，在动态计算域内，求解可压缩纳维-斯托克斯方程当前时间步的残差，并进行时间积分迭代求解网格参数。
[0067]
求解可压缩纳维-斯托克斯方程当前时间步的残差，并进行时间推进，更新动态计算域的网格参数。
[0068]
s3包括：
[0069]
s31，基于流场初始化给定的所有网格参数，求解可压缩纳维-斯托克斯方程，并选择预处理后的roe格式，计算对流通量，获得当前时间步的残差。
[0070]
基于流场初始化给定的所有网格单元的参数，进行残差估计，残差估计得到的是更新的网格参数。首先给定了网格参数的初始值，在残差估计后得到网格参数的的更新量。残差的获得可以理解为采用如下表达式：
[0071]
残差＝初始值-下一步的值
[0072]
s32，基于残差结果，进行时间积分迭代求解，更新获得当前时间步的网格参数。
[0073]
s4，基于s3得到的残差和网格参数，对动态计算域进行更新。
[0074]
基于当前时间步的网格参数，根据预处理后的速度和声速，判断对流动态域是否需要增大，再结合残差的收敛情况，判断对流动态域是否需要减少；
[0075]
判断当前网格单元是否受粘性效应主导，若受粘性单元主导，则将当前网格单元所有紧邻网格单元增加至粘性动态域，若不受粘性单元主导，则将当前网格单元从粘性动态域中删除。
[0076]
s4中对流动态域的增大与减小，具体包括：
[0077]
通过衡量s3中计算得到的网格参数更新量，判断当前网格单元是否受到无粘扰动，若受到无粘扰动，确定会受影响的紧邻网格单元，对于处于亚声速流动的单元，所有紧邻网格单元都会受到扰动的影响，将所有紧邻网格单元加入到对流动态域中；对于处于超声速流动的单元，则根据对流通量的特征值确定扰动的传播方向，将传播方向下游的网格单元加入对流动态域中；
[0078]
判断当前网格单元是否已收敛，若收敛，则判断当前网格单元是否位于对流动态域的最上游，若已位于最上游，则判断对流动态域中其他网格单元是否不再对当前网格单元产生影响，若其他网格单元不再对当前网格单元产生影响，则从对流动态域中删除当前网格单元。
[0079]
预处理后的速度和声速是基于当地速度和声速获得的，具体包括：
[0080]
基于当地速度u与声速c获得预处理后的速度和声速；
[0081]
基于预处理后速度和声速，获得预处理后的特征值，
[0082]
其中，预处理后的速度和声速的计算式分别如下：
[0083]
λ＝[u,u,u,u
′
c
′
,u
′‑c′
]
[0084][0085]ur2
＝min(c2,max(|u|2,k|u|2))
[0086]
其中，λ为计算对流通量时用到的特征值，u为当地速度，u
′
为预处理后的速度，c为声速，c
′
为预处理声速，α、β都是求解预处理特征值所用的参数，ur为预处理参数，由预处理方法给出，ρ,ρ
p
,ρ
t
分别为密度、密度对压力的导数、密度对温度的导数，由流场参数计算得到，c
p
为定压比热容常数，k是一个小值(10-5
)。
[0087]
s4中粘性动态域的增大与减小，具体包括：
[0088]
将总焓变化量的阈值除以马赫数平方作为粘性判断依据，进行粘性动态域的增大与缩小，具体方式如下：
[0089]
若|1-h
0,∞
/h0|/ma2≥ε
d,v
，则将当前网格单元的所有紧邻网格单元增加至粘性动态域；
[0090]
若|1-h
0,∞
/h0|/ma2《ε
d,v
，则将当前网格单元从粘性动态域中删除；
[0091]
其中，h
0,∞
为来流总焓，h0为当前网格单元总焓，ma为当前网格单元马赫数，ε
d,v
为给定的粘性域判断阈值。
[0092]
网格参数包括：
[0093]
密度、速度、压力和温度。
[0094]
示例性地，在本实施例中粘性域判断阈值取0.001，该数值适用于扰动域更新法针对可压缩流动计算的要求。针对本实施例，对比是否除以马赫数平方进行粘性分区判断的计算结果。以靠近壁面的一个网格单元为例，该网格必定处于粘性主导区域。但无量纲化计算得到的来流总焓为25000.5，当前网格单元总焓为25002.0064556779，原判断依据|1-h
0,∞
/h0|计算所得的结果为6.025
×
10-5
，远小于给定的粘性判断阈值。在除以马赫数平方后，该网格单元判断依据|1-h
0,∞
/h0|/ma2＝0.6025，即可判断得到该网格单元处于粘性主导区域，证明了本项技术手段能够有效地判断不可压流动的粘性分区，而在更高速度的可压缩流动中，马赫数平方趋近于1，也未改变原判断依据对可压缩流动的适用性，即本发明采用的技术手段可以适用于可压/不可压混合流动的计算。
[0095]
从总体来看，在采用除以马赫数平方之前，本实施例的粘性主导区域始终保持在1.4％(即初始粘性域)，失去了计算的正确性。而采用除以马赫数平方之后，粘性主导区域占比大小如图2所示，十分明显，更加符合物理规律，并可得到准确的计算结果。
[0096]
判断并执行粘性动态域的缩小。粘性动态域的缩小与上一步粘性域的扩大对应，当粘性效应不占据主导因素时，则可将该网格单元从粘性动态域中删除，即缩小了粘性动态域的范围。
[0097]
s5判断全局残差是否小于给定收敛残差阈值，若是，则计算完成并执行s6，若否，则返回s3。
[0098]
通过全局残差与给定收敛残差阈值的大小关系进行判定计算是否完成，即全局残差是否小于给定收敛残差阈值，若是，则进行s6，否则，跳转至s3。
[0099]
(1)若全局残差小于给定阈值，则计算完成，进入s6。
[0100]
r《ε
lim
[0101]
其中，r为全局残差，ε
lim
为给定收敛残差阈值。
[0102]
(2)当全局残差大于等于给定的收敛残差阈值时，则跳转至s3。
[0103]
r≥ε
lim
[0104]
s6，输出计算所得网络参数的值。
[0105]
在计算完成后，输出计算所得的各物理量的数值。
[0106]
实施例1
[0107]
图1为本发明提供的一种可压/不可压混合流动的高效扰动域推进方法的示意图，下面进行具体介绍。本发明提供了一种适用于可压/不可压流动的扰动域更新方法，以二维翼型naca0012流场求解为例，将预处理方法应用于可压缩纳维-斯托克斯方程求解的实施方式包括如下步骤：
[0108]
第一步：数据读入与流场初始化。
[0109]
(1)读入naca0012翼型的具体计算参数，本实施例的具体参数为：马赫数0.01，雷诺数为10000，层流有粘流动；
[0110]
(2)读入naca0012翼型的计算网格和边界条件的设置；
[0111]
(3)本实施例采用根据来流条件初始化的方法，采用来流边界条件初始化所有计算网格上的流场参数。
[0112]
第二步：建立初始动态计算域。
[0113]
(1)由于采用根据来流条件初始化的方法，仅根据壁面边界建立对流、粘性动态
域；
[0114]
第三步：残差估计与时间积分迭代求解。
[0115]
(1)根据不同的计算条件，选择合适的重构格式与对流格式，求解可压缩纳维-斯托克斯方程当前时间步的残差。在本实施例中，选择预处理后的roe格式，计算对流通量，获得当前时间步的残差；其中，计算对流通量时用到特征值，特征值包括当地速度和声速。
[0116]
(2)得到残差后，采用显式或隐式时间推进求解的方法进行时间推进，更新得到当前时间步计算域的物理量。在本实施例中，选择隐式lu-sgs格式进行时间推进求解，获得当前时间步的物理量更新量。
[0117]
第四步：动态域更新
[0118]
(1)判断并执行对流域的增大。首先通过衡量第三步中计算得到的物理量更新量，判断当前网格单元是否受到无粘扰动。若受到无粘扰动，则需确定会受影响的紧邻单元。若单元处于亚声速流动，所有紧邻单元都会受到扰动的影响，若单元处于超声速流动，则根据对流通量的特征值确定扰动的传播方向，从而确定受扰单元，将传播方向下游的网格单元加入对流动态域中。
[0119]
计算对流通量时使用到的特征值需采用预处理后的特征值，该特征值计算采用的预处理后的速度和声速是基于当地速度和声速获得的。
[0120]
在这一步中，计算对流通量时使用到的特征值需采用预处理后的特征值，其求解方法为：
[0121]
λ＝[u,u,u,u
′
c
′
,u
′‑c′
]
[0122]
预处理后的速度和声速的计算式如下：
[0123][0124]
其中，λ为计算对流通量时用到的特征值，u为当地速度，u
′
为预处理后的速度，c为声速，c
′
为预处理声速，ur为预处理参数，由预处理方法给出，ρ,ρ
p
,ρ
t
分别为密度、密度对压力的导数、密度对温度的导数，由流场参数计算得到，c
p
为定压比热容常数。
[0125]
(2)判断并执行对流域的缩小。首先判断当前单元是否已经收敛，若收敛，则判断当前网格单元是否位于对流动态计算域的最上游，若已位于最上游，则需判断对流动态域中其他单元是否不再对该单元产生影响，则从对流动态域中删除当前网格单元。
[0126]
若其他单元不再对该单元产生影响，则可以进行对流域的缩小。
[0127]
(3)判断并执行粘性域的增大。判断该单元是否受粘性效应主导，若受粘性单元主导，则将所有紧邻单元增加至粘性域。本发明提出了适用于预处理求解可压/不可压流场的粘性主导判断依据：
[0128]
|1-h
0,∞
/h0|
×
ma2≥ε
d,v
[0129]
其中，h
0,∞
为来流总焓，h0为当前单元总焓，ma为当前单元马赫数，ε
d,v
为给定的粘性域判断阈值，在本实施例中取0.001，该数值适用于扰动域更新法针对可压缩流动计算的要求。
[0130]
针对本实施例，对比是否除以马赫数平方进行粘性分区判断的计算结果。以靠近壁面的一个网格单元为例，该网格必定处于粘性主导区域。但无量纲化计算得到的来流总
焓为25000.5，当前单元总焓为25002.0064556779，原判断依据|1-h
0,∞
/h0|计算所得的结果为6.025
×
10-5
，远小于给定的粘性判断阈值。在除以马赫数平方后，该单元判断依据|1-h
0,∞
/h0|
×
ma2＝0.6025，即可判断得到该单元处于粘性主导区域，证明了本项技术手段能够有效地判断不可压流动的粘性分区，而在更高速度的可压缩流动中，马赫数平方趋近于1，也未改变原判断依据对可压缩流动的适用性，即本性技术手段可以适用于可压不可压流动的计算。
[0131]
图2是本发明提供的对流动态域与粘性主导区域占比大小的示意图。
[0132]
从总体来看，若无除以马赫数平方，本实施例的粘性主导区域始终保持在1.4％(即初始粘性域)，失去了计算的正确性。而除以马赫数平方后，粘性主导区域占比大小如图2所示，明显更加符合物理规律，并可得到正确的计算结果。
[0133]
如图2所示，size_a和size_v分别指对流动态域和粘性动态域在所有计算域中的占比。
[0134]
(4)判断并执行粘性域的缩小。粘性域的缩小与上一步粘性域的扩大对应，当粘性效应不占据主导因素时，则可将该单元从对流域中删除。
[0135]
第五步：判断计算是否完成。
[0136]
(1)判断全局残差与给定收敛残差阈值的大小关系，若全局残差小于给定阈值，则计算完成，进入第六步。
[0137]
r《ε
lim
[0138]
其中，r为全局残差，ε
lim
为给定收敛残差阈值，即若全局残差小于给定阈值，则计算完成，进入第六步。
[0139]
(2)反之，若全局残差仍大于给定阈值，则循环进行至第三步。
[0140]
第六步：输出计算结果。
[0141]
在计算完成后，输出计算所得的各物理量的数值。
[0142]
图3a是本发明实施例提供的未采用可压/不可压流动的扰动域更新方法计算结果的对比示意图；图3b是本发明实施例提供的采用可压/不可压流动的扰动域更新方法计算结果的对比示意图。
[0143]
图4是本发明实施例提供的是否采用可压/不可压流动的扰动域更新方法收敛情况的对比示意图。
[0144]
如图4所示，drum指本发明的方法，gum指未采用本发明的方法，横坐标为迭代步数，纵坐标为最大残差收敛情况。
[0145]
在本实施例完成计算后，是否采用本发明提出的可压/不可压流动的扰动域更新方法计算结果和计算收敛情况对比如图3a、图3b和图4所示，可以看出，翼型表面的压力分布二者完全一致，总体残差收敛情况也基本一致，即可说明采用本发明提出的方法，保证了计算的准确性。而应用本发明提出计算方法后，本实施例的总计算时间由40183.3s减少到了26197.6s，节省了35％的时间，证明了本发明可以大幅提高可压/不可压混合流动的计算效率。
[0146]
根据本公开实施例提供的技术方案，基于预处理后的当地声速与当地速度，建立可压/不可压扰动域更新方法的扰动传播计算公式，从而建立可压/不可压扰动域传播和影响的判据，解决了雅克比矩阵特征值改变和低马赫数区域判据不适用的问题，收缩计算域；
提出了适用于预处理求解可压/不可压流场的粘性主导判断依据，将总焓变化量的阈值除以马赫数平方作为粘性分区的判断阈值作为粘性分区的判断阈值，避免了采用总焓变化量作为判据时粘性分区对马赫数的依赖性，大幅提高计算效率。
[0147]
上述所有可选技术方案，可以采用任意结合形成本技术的可选实施例，在此不再一一赘述。
[0148]
应理解，上述实施例中各步骤的序号的大小并不意味着执行顺序的先后，各过程的执行顺序应以其功能和内在逻辑确定，而不应对本公开实施例的实施过程构成任何限定。
[0149]
以上实施例仅用以说明本公开的技术方案，而非对其限制；尽管参照前述实施例对本公开进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本公开各实施例技术方案的精神和范围，均应包含在本公开的保护范围之内。