用于高超声速强激波流场气动热预测的湍流模型修正方法与流程

1.本发明涉及计算流体力学技术领域，尤其涉及一种用于高超声速强激波流场气动热预测的湍流模型修正方法。


背景技术：

2.当飞行器在较为稠密的大气内做高超声速飞行时，由于来流总温很高，将引起十分严重的气动加热问题，并且此时的流动状态往往为湍流，会进一步加剧气动加热的严重程度。因此，为保证飞行安全，必须对气动热环境进行准确预测，尤其需要对湍流进行准确模拟。
3.随着数值方法与计算能力的发展，计算流体力学方法逐渐成为气动热环境预测的重要手段，其中对湍流的数值模拟主要基于雷诺平均(rans)方法开展。然而，高超声速流场中强激波间断的出现将严重影响湍流模型的准确性：由于湍流模型变量输运方程中的生成项表达式中包含一系列速度导数项，而激波捕捉算法是以高梯度来近似激波间断，因此在流场中激波附近的高速度梯度值将是非物理的，导致在激波处湍流模型变量的生成项会呈现过度增加，从而引起湍流变量如湍动能、比耗散率、湍流涡粘系数等的非物理增长，进而会显著影响壁面热流密度的模拟精度。例如，“sinha k,mahesh k,candler g v.modeling shock unsteadiness in shock/turbulence interaction[j].physics of fluids,2003,15(8):2290
‑
2297”中对比了湍流模型与直接数值模拟(dns)计算得到的正激波前后湍动能变化，如图1所示，未修正的原始湍流模型(standard model)得到的激波后湍动能显著高于dns结果。
[0004]
针对上述问题，目前国内外所见研究较少，主要为sinha与zhang针对湍流模型变量输运方程生成项提出的修正方法。“sinha k.，mahesh k.，candler g.v.modeling the effect of shock unsteadiness in shock/turbulent boundary
‑
layer interactions[j].aiaa journal，2005,43(3):586
‑
594”重新构造了输运方程中生成项与耗散项在激波处的表达式，通过dns的数值模拟结果确定表达式中模型系数的形式，使得激波处湍流变量的变化情况同dns的结果一致，原有的非物理增长现象消失。“zhang z,gao z,jiang c,et al.arans model correction on unphysical over
‑
prediction of turbulent quantities across shock wave[j].international journal of heat and mass transfer,2017,106:1107
‑
1119”基于weno格式中模板光滑因子的思想，提出了以无量纲压强作为检测变量的适用于全流场的激波捕捉函数，并在激波处构造湍动能衰减函数，使激波处的湍动能生成项为0，同样能够消除激波处湍流变量的非物理增长。
[0005]
sinha等人的修正方法，主要存在两个问题：第一，修正生成项中模型系数由激波波前法向马赫数得到，该参数在三维复杂流动模拟中很难在cfd程序中获得；第二，所用的激波检测函数需要得到激波根部的边界层厚度，这使得该方法无法应用至脱体激波存在的流动。zhang等人的修正方法所用激波检测函数形式过于复杂，实际应用中会增加一定的计算量，并且，改进的激波检测函数对于激波检测的准确性在一定程度上依赖于函数中的系
数取值，对于存在复杂波系的流场，一旦将某些不存在激波的湍流区域误判为激波，则会令湍流变量生成项为零进而引起很大的误差。
[0006]
因此，需要发展一种通用性更强、更为可靠，同时不显著增加计算量的适用于强激波流场的湍流模型修正方法。


技术实现要素：

[0007]
有鉴于此，本发明提供了一种用于高超声速强激波流场气动热预测的湍流模型修正方法，用以消除湍流模型变量在激波后的非物理增长，提供一种通用性更强、更为可靠的能够用于复杂波系存在的流动的湍流模型修正方法。
[0008]
本发明提供的一种用于高超声速强激波流场气动热预测的湍流模型修正方法，包括如下步骤：
[0009]
s1：针对高超声速飞行器的外形，构建结构化计算网格；
[0010]
s2：采用沿时间推进的求解方法，在每个时间步内，计算所述结构化计算网格中各网格点上基于应变率的原始生成项；
[0011]
s3：根据当前时间步下各网格点的流体流动速度以及湍流涡粘系数，计算涡量大小，得到各网格点上基于涡量的修正生成项；
[0012]
s4：比较所述原始生成项与所述修正生成项的量级差距，判断所述原始生成项是否高于所述修正生成项一个量级以上；若是，则所述湍流模型中的生成项采用所述修正生成项，若否，则所述湍流模型中的生成项采用所述原始生成项。
[0013]
在一种可能的实现方式中，在本发明提供的上述用于高超声速强激波流场气动热预测的湍流模型修正方法中，步骤s2，具体包括：
[0014]
k
‑
ω湍流模型与k
‑
ωsst湍流模型中湍动能的原始生成项p
k
的表达式为：
[0015][0016]
其中，μ
t
表示湍流涡粘系数，u表示流体流动速度，x表示当地坐标，ρ表示流体密度，k
t
表示湍动能；上标
“‑”
表示雷诺平均量，上标“～”表示favr
é
平均量，下标i,j,k分别表示笛卡尔坐标系的三个空间分量；δ
ij
为kronecker符号。
[0017]
在一种可能的实现方式中，在本发明提供的上述用于高超声速强激波流场气动热预测的湍流模型修正方法中，步骤s3，具体包括：
[0018]
涡量ω的表达式为：
[0019][0020]
从涡量出发，考虑量纲一致的原则，重新构造k
‑
ω湍流模型与k
‑
ωsst湍流模型的修正生成项p
k_mod
的表达式为：
[0021]
p
k_mod
＝μ
t
ω2[0022]
(3)。
[0023]
在一种可能的实现方式中，在本发明提供的上述用于高超声速强激波流场气动热预测的湍流模型修正方法中，步骤s4中，所述湍流模型中的生成项的最终表达式为：
[0024][0025]
本发明提供的上述用于高超声速强激波流场气动热预测的湍流模型修正方法，通过将湍流模型变量输运方程中原始的基于应变率构造的生成项在激波处改为基于涡量构造的生成项，从而消除生成项在激波处的非物理增长，达到消除激波后湍流模型变量非物理增长的目的。与现有的修正方法相比，本发明不仅能够显著消除湍流模型变量在激波后的非物理增长，而且对附体激波与脱体激波均有良好的适用性，与原始湍流模型相比增加的计算量较小，并且，可靠性与容错率更高，能够用于复杂波系存在的流动。
附图说明
[0026]
图1为现有技术中正激波前后湍动能变化图；
[0027]
图2为本发明提供的一种用于高超声速强激波流场气动热预测的湍流模型修正方法的流程图；
[0028]
图3为本发明实施例1中钝头体外形与几何尺寸示意图；
[0029]
图4为本发明实施例1中钝头体流场对称面湍流涡粘系数云图；
[0030]
图5为本发明实施例1中钝头体流场对称面迎风面热流密度图；
[0031]
图6为本发明实施例1中钝头体流场对称面生成项分布图；
[0032]
图7为本发明实施例1中湍流模型修正前后钝头体流场对称面涡粘系数云图；
[0033]
图8为本发明实施例1中湍流模型修正前后钝头体流场对称面迎风面热流密度图。
具体实施方式
[0034]
下面将结合本发明实施方式中的附图，对本发明实施方式中的技术方案进行清楚、完整的描述，显然，所描述的实施方式仅仅是作为例示，并非用于限制本发明。
[0035]
本发明提供的一种用于高超声速强激波流场气动热预测的湍流模型修正方法，如图2所示，包括如下步骤：
[0036]
s1：针对高超声速飞行器的外形，构建结构化计算网格；
[0037]
s2：采用沿时间推进的求解方法，在每个时间步内，计算结构化计算网格中各网格点上基于应变率的原始生成项；
[0038]
s3：根据当前时间步下各网格点的流体流动速度以及湍流涡粘系数，计算涡量大小，得到各网格点上基于涡量的修正生成项；
[0039]
s4：比较原始生成项与修正生成项的量级差距，判断原始生成项是否高于修正生成项一个量级以上；若是，则湍流模型中的生成项采用修正生成项，若否，则湍流模型中的生成项采用原始生成项。
[0040]
下面通过一个具体的实施例对本发明提供的上述用于高超声速强激波流场气动热预测的湍流模型修正方法的具体实施进行详细说明。为了证明本发明提供的上述用于高超声速强激波流场气动热预测的湍流模型修正方法的有效性，选取sch
ü
lein等人所做的高超声速钝头体表面热流测量实验数据作为对比数据，具体请参见sch
ü
lein e.,krogmann p.,stanewsky e.documentation of two
‑
dimensional impinging shock/turbulent boundary layer interaction flow[r].german aerospace center，report ib 223
‑
96a49。
[0041]
实施例1：
[0042]
第一步，针对高超声速飞行器的外形，构建结构化计算网格。
[0043]
第二步，采用沿时间推进的求解方法，在每个时间步内，计算结构化计算网格中各网格点上基于应变率的原始生成项。
[0044]
首先分析湍流模型变量在激波间断附近的变化情况及其对气动热环境预测的影响。所选取的湍流模型为k
‑
ω湍流模型与k
‑
ωsst湍流模型，所选取的钝头体外形与几何尺寸如图3所示。图4为应用k
‑
ω湍流模型(如图4中的(a)所示)与k
‑
ωsst湍流模型(如图4中的(b)所示)计算得到的钝头体流场对称面湍流涡粘系数云图，并以来流的湍流涡粘系数μ
t∞
对钝头体流场对称面的湍流涡粘系数进行无量纲化，表示为μ
t
/μ
t∞
。由图4可以看出，在对应脱体激波的位置，k
‑
ω湍流与k
‑
ωsst湍流模型计算得到的湍流涡粘系数在激波后均呈现非物理急剧增长，最大湍流涡粘系数分别约为激波前的60000倍与10000倍。同时，湍流涡粘系数非物理增长的影响区域均波及边界层内部，尤其会使驻点区附近边界层内的湍流涡粘系数显著偏高，进而导致壁面热流密度q
w
的增大，如图5所示，图5中的纵坐标q
w
采用热流密度常用单位“千瓦每平方米”，横坐标以模型全长l对当地坐标x进行无量纲化。由图5可以看出，在驻点附近，k
‑
ω湍流模型与k
‑
ωsst湍流模型计算得到的热流密度均显著高于对比数据(exp.)，其中，k
‑
ω湍流模型偏高约300％，k
‑
ωsst湍流模型偏高约100％。
[0045]
然后，由于湍流涡粘系数在激波后的非物理增长是由输运方程的生成项在激波处不适用导致的，因此，可以针对输运方程的生成项进行修正。k
‑
ω湍流模型与k
‑
ωsst湍流模型中湍动能的原始生成项p
k
的表达式为：
[0046][0047]
其中，μ
t
表示湍流涡粘系数，u表示流体流动速度，x表示当地坐标，ρ表示流体密度，k
t
表示湍动能；上标
“‑”
表示雷诺平均量，上标“～”表示favr
é
平均量，下标i，j，k分别表示笛卡尔坐标系的三个空间分量；表示流体在笛卡尔坐标系下i方向的favr
é
平均流动速度分量，表示流体在笛卡尔坐标系下j方向的favr
é
平均流动速度分量，表示流体在笛卡尔坐标系下k方向的favr
é
平均流动速度分量；表示流体的雷诺平均密度，x
i
、x
j
、x
k
分别表示笛卡尔坐标系下i、j、k三个方向的坐标，δ
ij
为kronecker符号。
[0048]
第三步，根据当前时间步下各网格点的流体流动速度以及湍流涡粘系数，计算涡量大小，得到各网格点上基于涡量的修正生成项。
[0049]
从公式(1)可以看出，由于生成项基于线变形项构造，表达式中包含速度分量对该分量所在方向的导数项。而激波本身作为强间断结构，在数值计算中被视为高梯度区域，导致线变形项在此处急剧增加，生成项产生非物理的极大值。因此，本发明提出将生成项修正为基于涡量构造。涡量ω的表达式为：
[0050][0051]
由于涡量中不包含线变形项，仅包含角速度项，即速度分量对其他方向的导数，因
此，在激波处将不会产生非物理的极大值。从涡量出发，考虑量纲一致的原则，重新构造k
‑
ω湍流模型与k
‑
ωsst湍流模型的修正生成项p
k_mod
的表达式为：
[0052]
p
k_mod
＝μ
t
ω2(3)。
[0053]
尽管在激波处基于涡量的修正生成项会显著低于基于线变形的原始生成项，但在流场其他区域，诸如边界层或分离区等涡量较大的区域，修正生成项有可能高于原始生成项，进而改变原始模型对边界层与分离流动模拟的特性。因此，直接在全流场采用上述修正生成项是不合适的。基于此，可以执行下面的第四步。
[0054]
第四步，比较原始生成项与修正生成项的量级差距，判断原始生成项是否高于修正生成项一个量级以上。如果原始生成项高于修正生成项一个量级以上，那么湍流模型中的生成项采用修正生成项；如果原始生成项小于或等于修正生成项一个量级以上，那么湍流模型中的生成项采用原始生成项。
[0055]
为进一步定量探究这一问题，图6对比了钝头体流场中不同位置原始生成项与修正生成项的大小。图6中的(a)为钝头体流场中所选取的三个位置的示意图，依次记为1、2、3号线。图6中的(b)、(c)、(d)分别为1、2、3号线上原始生成项与修正生成项的分布情况，纵坐标中，生成项p
k*
的大小以来流参数无量纲化，p
k*
分别取值p
k
和p
k_mod
，ρ
∞
为来流密度，v
∞
为来流速度，l
∞
为参考长度，l
∞
在本发明实施例1中取1m，横坐标s表示距离壁面的法向距离，其单位为m。由图6中的(b)可以看出，两种生成项在激波处的大小存在量级上的显著差距，而在图6中的(c)与(d)中，两种生成项在边界层内呈现较小的差别。根据这个显著特征，可以以量级的差距来判断进行生成项修正的区域，即湍流模型中的生成项的最终表达式为：
[0056][0057]
下面给定来流条件如下：马赫数ma＝6.03，静压p
∞
＝58.6pa，静温t
∞
＝216.65k，壁面温度t
w
＝300k，迎角α＝40
°
。利用以上提出的针对k
‑
ω湍流模型与k
‑
ωsst湍流模型生成项的修正模型开展数值模拟。图7对比了修正模型与原始模型计算得到的对称面湍流涡粘系数云图，虚线以上部分为未修正结果，虚线以下部分为修正结果。由图7可以看出，修正后的k
‑
ω湍流模型(如图7中的(a)所示)与k
‑
ωsst湍流模型(如图7中的(b)所示)的湍流涡粘系数在激波处没有产生明显变化，原有的非物理增长被消除。图8对比了k
‑
ω湍流模型(如图8中的(a)所示)与k
‑
ωsst湍流模型(如图8中的(b)所示)修正前后钝头体迎风面对称面的热流密度分布，由图8可以看出，在驻点处k
‑
ω湍流模型与k
‑
ωsst湍流模型原本严重偏高的热流密度显著下降，均同对比数据(exp.)吻合良好。以上结果表明，本发明提供的上述用于高超声速强激波流场气动热预测的湍流模型修正方法，能够消除脱体激波后湍流变量非物理增长现象，提高气动热的预测精度。
[0058]
本发明提供的上述用于高超声速强激波流场气动热预测的湍流模型修正方法，通过将湍流模型变量输运方程中原始的基于应变率构造的生成项在激波处改为基于涡量构造的生成项，从而消除生成项在激波处的非物理增长，达到消除激波后湍流模型变量非物理增长的目的。本发明提出的修正方法不仅能够显著消除湍流模型变量在激波后的非物理增长，而且对附体激波与脱体激波均有良好的适用性，sinha等人的修正方法需要得到激波根部的边界层厚度，而钝头体超声速绕流产生的脱体激波不存在激波根部区域，故无法使
用sinha等人的修正方法，本发明提出的修正方法仅需要使用流场中的速度导数，不受激波类型的限制。本发明提出的修正方法与原始湍流模型相比增加的计算量较小，zhang等人构造的激波检测函数基于weno格式中模板光滑因子的思想，以无量纲压强作为检测变量，计算量很大，而本发明提出的修正方法通过比较流场中原始生成项与修正生成项的量级大小检测激波位置，计算方法简单，增加的计算量很小。本发明提出的修正方法的可靠性与容错率更高，能够用于复杂波系存在的流动，zhang等人的修正方法在检测到激波后，对生成项的处理方式是将生成项大小设为0，然而激波检测的准确性在一定程度上依赖于函数中的系数取值，如果激波检测有误，此处生成项大小为0将严重影响整个流场，导致错误的计算结果，而本发明中修正生成项由涡量构造，在流场中除激波以外的区域，修正生成项同原始生成项的大小没有显著差异，即使在部分区域错误地将原始生成项改为修正生成项，也不会使计算结果产生显著的错误。
[0059]
显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。