一种边界层分布式粗糙元诱导转捩预测方法

1.本发明属于转捩预测技术领域，具体涉及一种可预测不同壁面分布式粗糙度下的边界层转捩预测方法。


背景技术：

2.围绕壁面粗糙度作用的转捩问题，国内外相关研究工作已经取得了一系列较为显著的研究成果。但受限于问题本身的复杂性，对于壁面粗糙度的转捩影响机理并未形成统一完整的解释，这使得当前对粗糙元诱导转捩的预测不得不引入一系列的经验拟合，且只对特定外形具有较高的预测精度。
3.这其中，通过建立转捩参量和扰动参量的拟合式发展而来的转捩准则，便是最具代表性的方法之一。转捩参量y
tr
和扰动参量x的关系式可表示为y
tr
＝axn(a》0,n《0)，进而通过大量飞行和地面风洞实验数据对a和n的值进行标定。目前，在高超声速飞行器转捩研究中，公布的经验公式已达数十种，扰动参数与转捩参数的相关关系也不断优化，并逐步确立了re
θ
/me准则和rek准则的主导地位。比较著名的如：美国兰利研究中心建立的基于孤立粗糙元的re
θ
/me准则；reda等人针对具有边界层分布式粗糙元的钝头体飞行器发展的分段式rek准则等。
4.由经验拟合的转捩准则方法严重依赖于拟合的数据库，不考虑转捩过程的物理机制，因此不具备普适性，但其依旧为基于rans的粗糙元转捩模型的构建和发展带来相当多的启发。具体而言，尽管转捩模型对边界层转捩过程的模化基于间歇因子输运方程，但其对粗糙壁面下转捩起始位置的判据，则主要借鉴于转捩准则方法及其研究基础。特别是转捩准则方法中发展的re
θ
/me、rek以及粗糙元等效高度ks等概念，在粗糙元转捩模型中依旧具有重要意义。
5.在转捩模型中，考虑粗糙元对边界层转捩影响的方式主要有两种。第一种方法是stripf等提出的所谓离散粗糙元法，该方法通过在输运方程中引入一些附加项来描述粗糙元后流动的阻塞效应。这一方法虽然避免了ks的局限性，但却增加了计算量。类似的，易淼荣和赵慧勇等通过增大γ-re
θ
模型中分离诱导转捩参数s1的值，来提高γ-re
θ
模型对分离诱导转捩的敏感性，从而提高模型对粗糙元强制诱导转捩的预测能力，但普适性较差。第二种基于等效粗糙元高度(ks)的概念，而忽略粗糙元之间的相互作用。stripf等和elsner等利用等效粗糙元高度ks和边界层位移厚度δ*对转捩起始位置判据进行了修正，以使其反映壁面粗糙效应的影响。但边界层位移厚度的引入使得该方法不严格基于当地变量，不适用于现代cfd方法。此外，根据代数拟合建立的转捩判据无法反映壁面粗糙元的非当地效应和历史效应。
6.以上分析可知，相较于转捩准则，基于rans的转捩模型方法考虑了一定的粗糙元诱导转捩的物理机理，从而一定程度上提高了模型的预测性能和普适性。但距离实际三维外形的粗糙元诱导转捩预测仍有一定距离。


技术实现要素：

7.本发明提供一种边界层分布式粗糙元诱导转捩预测方法，引入了粗糙元放大因子(ar)输运方程来描述壁面粗糙度对边界层转捩的影响机理和作用效果以体现壁面粗糙，ar与壁面粗糙程度直接相关，由边界条件给定，方程的对流和扩散项允许ar被输运到流场中，以此反映壁面粗糙的历史效应和非当地效应对转捩的影响，并且实现了完全当地化。
8.为实现上述技术目的，本发明采用如下技术方案：
9.一种边界层分布式粗糙元诱导转捩预测方法，用于对飞行器的表面边界层进行分离诱导转捩预测，包括以下步骤：
10.s1：根据飞行器外形壁面的粗糙程度确定等效粗糙度ks，根据壁面摩擦速度u
τ
和等效粗糙度ks，按表达式(1)计算飞行器外形壁面的等效沙粒粗糙元高度k
+
；
[0011][0012]
式中，下角标w表示壁面量，ρw，τw，μw，νw分别为壁面密度、壁面剪切应力、壁面分子动力粘性系数和壁面运动粘性系数；
[0013]
s2：建立有关于粗糙元放大因子ar的输运方程如表达式(2)，并按表达式(3)基于等效沙粒粗糙元高度k
+
给定粗糙元放大因子ar的边界条件ar|
wall
：
[0014][0015][0016]
式中，ρ为密度，t为时间，x为坐标，μ为分子粘性系数，μ
eff
为有效涡粘性粘性系数；u为速度，j＝1,2,3；为模式封闭参数；k
+
是飞行器外形壁面无量纲的等效沙粒粗糙元高度；
[0017]
s3：以原始γ-re
θ
转捩模型为基本框架，向其动量厚度输运方程中增加粗糙元源项实现粗糙元诱导转捩的预测；其中粗糙元源项与粗糙元放大因子ar之间的关系式为：
[0018][0019]
其中，c
θt
和为模式封闭参数，是采用γ-re
θ
模型中的输运方程求得的动量厚度雷诺数，f
θt
为原始γ-re
θ
转捩模型中的边界层识别参数，re
scf
为横流雷诺数，f(ar)为粗糙度控制参数。
[0020]
在更优的技术方案中，将原始γ-re
θ
转捩模型中湍动能比耗散率ω的输运方程的壁面边界条件修正为：
[0021][0022]
其中：
[0023][0024]
在更优的技术方案中，将粗糙元源项增加在动量厚度雷诺数输运方程中，得边界层分布式粗糙元转捩预测模型，如表达式(6)所示：
[0025][0026]
式中，p
θt
表示该输运方程的源项，σ
θt
为模型封闭参数，μ
t
表示涡粘性系数。
[0027]
在更优的技术方案中，步骤1中等效粗糙度ks的确定方法为：取粗糙元均方根高度ra的2.8～4.8倍。
[0028]
在更优的技术方案中，粗糙度控制参数的表达式为：
[0029][0030]
其中，c
ar1
和c
ar2
为模式封闭参数。
[0031]
在更优的技术方案中，所述横流雷诺数re
scf
，基于现有风洞实验数据，通过cfd层流解得到不同壁面粗糙度的临界动量厚度雷诺数，采用最小二乘法对关系系数进行求解，得到判据关系为：
[0032][0033]
式中，θ表示边界层动量厚度，因同时出现在等式两边，可通过牛顿迭代法求解。
[0034]
有益效果
[0035]
与传统粗糙元诱导转捩预测技术相比，本发明具有以下优势：
[0036]
1、传统粗糙元诱导转捩预测方法较少反映壁面粗糙度的非当地效应和历史效应，本发明引入一个粗糙元放大因子ar输运方程，方程的对流和扩散项允许ar被输运到流场中，从而避免了该缺陷；
[0037]
2、传统方法大多基于经验拟合，本发明构造的粗糙元放大因子ar与壁面粗糙程度直接相关，其由边界条件给定，从而可反映粗糙元诱导转捩的物理机理；
[0038]
3、本发明采用牛顿迭代法求解边界层动量厚度，实现了转捩判据的完全当地化，适用于现代基于大规模并行和非结构网格的计算方法。
附图说明
[0039]
图1是本发明粗糙元诱导转捩预测模式总体框架；
[0040]
图2是本发明等效粗糙度与粗糙元均方根高度ra之间的关系曲线；
[0041]
图3是本发明粗糙元放大因子的壁面边界条件ar|
wall
随k
+
的变化曲线；
[0042]
图4是本发明方案(记为γ-re
θ
_rough)预测的粗糙壁面球头壁面热流分布云图与风洞实验结果(记为exp.)的对比；
[0043]
图5是本发明方案(记为γ-re
θ
_rough)预测的壁面子午线热流值与风洞实验结果
(记为exp.)以及原始γ-re
θ
转捩模型(记为γ-re
θ
_orig)计算结果的对比。
具体实施方式
[0044]
下面对本发明的实施例作详细说明，本实施例以本发明的技术方案为依据开展，给出了详细的实施方式和具体的操作过程，对本发明的技术方案作进一步解释说明。
[0045]
本实施例公开了一种边界层分布式粗糙元诱导转捩预测方法，用于对飞行器的表面边界层进行分离诱导转捩预测，参考图1所示，包括以下步骤：
[0046]
s1：根据飞行器外形壁面的粗糙程度确定等效粗糙度ks，根据壁面摩擦速度u
τ
和等效粗糙度ks，按表达式(9)计算飞行器外形壁面无量纲的等效沙粒粗糙元高度k
+
；
[0047][0048]
式中，下角标w表示壁面量，ρw，τw，μw，νw分别为壁面密度、壁面剪切应力、壁面分子动力粘性系数和壁面运动粘性系数；
[0049]
已有研究发现，等效粗糙度ks与粗糙元均方根高度ra间存在线性关系，如图2所示，ks为ra的2.8～4.8倍。因此，取ks＝[2.8ra,4.8ra]。
[0050]
s2：建立有关于粗糙元放大因子ar的输运方程如表达式(10)；并以sigmoid函数，且基于等效沙粒粗糙元高度k
+
，对粗糙元放大因子ar给定边界条件ar|
wall
如表达式(11)：
[0051][0052][0053]
式中，ρ为密度，t为时间，x为坐标，μ为分子粘性系数，μ
eff
为有效涡粘性粘性系数；u为速度，j＝1,2,3；为模式封闭参数，用以表征非直接影响因素对边界层转捩过程的影响，其取值与飞行器气动外形以及来流条件相关；
[0054]
表达式(11)的等号右边的常数项，用以保证在光滑壁面条件下ar＝0。
[0055]
上述构造的粗糙元放大因子ar输运方程，其本质目的是放大转捩区粗糙元对边界层内小扰动的影响。已有研究表明，当k
+
《10，粗糙元对转捩的作用可忽略，此时可视为光滑壁面；当10《k
+
《70，壁面粗糙元开始影响转捩；当k
+
》70，壁面粗糙元对转捩其主导作用。图3所示为本文建立的ar壁面边界条件随k
+
的变化曲线，显然，当k
+
《10时，ar＝0；当10《k
+
《70，ar开始增长；当k
+
》70，ar迅速增长并达到其峰值。
[0056]
s3：以原始γ-re
θ
转捩模型为基本框架，向其动量厚度雷诺数输运方程中增加粗糙元源项得到改进的边界层分布式粗糙元转捩预测模型，以实现粗糙元诱导转捩的预测。
[0057]
其中，粗糙元源项与粗糙元放大因子ar之间的关系式为：
[0058][0059]
其中，c
θt
和为模式封闭参数，是采用γ-re
θ
模型中的输运方程求得的动量
厚度雷诺数，f
θt
为原始γ-re
θ
转捩模型中的边界层识别参数；
[0060]
re
scf
为横流雷诺数，定义如下：
[0061][0062]
式中θ
t
为边界层动量厚度，其存同时存在于等式两侧，可需通过牛顿迭代法求解，从而避免了积分求解动量厚度，实现了当地化处理。
[0063]
f(ar)为粗糙度控制参数，基于现有风洞实验数据，通过cfd层流解得到不同壁面粗糙度的临界动量厚度雷诺数，采用最小二乘法对关系系数进行求解，得到判据关系式如下：
[0064][0065]
其中，c
ar1
和c
ar2
为模式封闭参数。
[0066]
增加粗糙元源项后得到的动量厚度雷诺数输运方程，如表达式(15)所示：
[0067][0068]
式中，p
θt
表示该输运方程的源项，σ
θt
为模型封闭参数，μ
t
表示涡粘性系数。
[0069]
s4：当转捩到全湍流后，γ-re
θ
转捩模型退化为sst湍流模型。原始sst湍流模型中比耗散率ω的壁面条件为：
[0070][0071]
其中η为壁面法向距离。
[0072]
当壁面存在粗糙度时，上述边界条件得到的湍流区热流和比摩擦系数偏小。因此，对于全湍流区域的粗糙度壁面，为避免在全湍流区域预测得到的热流和比摩擦系数较实验值偏低，本实施例对原始sst湍动模型中的湍动能比耗散率(ω)的输运方程的壁面边界条件进行修正，修改后的表达式为：
[0073][0074]
其中：
[0075][0076]
总体而言，壁面粗糙度的增加将导致ar增大，从而使得f(ar)函数增大，进而促进输运方程中耗散项增大，导致转捩起始位置处减小，转捩位置提前，促进湍动能的发展，能量由层流动能向湍动能输运，从而促进转捩的发生。
[0077]
图4给出了采用本发明方案(记为γ-re
θ
_rough)预测的粗糙壁面球头壁面热流分布云图，及其与风洞实验结果(记为exp.)的对比；图5给出了本发明方案(记为γ-re
θ
_
rough)预测的壁面子午线热流值，及其与风洞实验结果(记为exp.)以及原始γ-re
θ
转捩模型(记为γ-re
θ
_orig)计算结果的对比。
[0078]
以上实施例为本技术的优选实施例，本领域的普通技术人员还可以在此基础上进行各种变换或改进，在不脱离本技术总的构思的前提下，这些变换或改进都应当属于本技术要求保护的范围之内。