随着飞机设计技术的不断发展,对气动数据的精度要求逐渐提高。边界层转捩现象的发生对表面摩擦阻力、流动分离等流场形态影响显著。因此在飞机的精细化设计中,能否准确模拟边界层转捩现象具有十分重要的工程意义。
在计算流体力学领域中,N-S(Navier-Stockes)方程及其各种简化形式占据主导地位,但对其直接数值求解(DNS方法)或是采用大涡模拟(LES)等方法,由于昂贵的计算代价而仅停留在对转捩机制和建模方法的研究上,无法应用于工程计算。目前,具有工程实用性的转捩预测方法主要有eN[1, 2, 3]和由Langtry-Menter[4, 5, 6]提出的 
转捩模型。eN是一种基于稳定性理论的半经验方法,认为转捩过程是由层流边界层内最初的小扰动向下游传播放大达到eN时完成,放大因子N依赖于风洞和自由来流的环境而并非是特定的。在三维流动中,同时存在迥异的流向速度型和横流速度型,加之eN方法要求沿流线追踪扰动放大率的增长率,均导致了其在三维流动问题应用中的局限性[3]。 转捩模型[4]是完全基于当地变量的转捩模型,该方法综合了经验关系式判断转捩起始位置,利用间歇因子来控制湍流的生成。它完全基于当地化变量,适用于计算流体力学并行计算,因而在实际工程中已得到广泛应用。
近年来,转捩模型由于在工程应用中的优势得到不断的发展。Krause和Behr等[7]对失稳雷诺数(Reθc)和表征转捩区长度的参数(Flength)的经验关系式进行了修改,将来流湍流度引入经验关系式中,取得了不错的效果。Coder和Maughmer[8]发展了该模型的简化形式,直接建立了输运雷诺数 
的求解关系式达到简化目的。Menter和Smirnov等人[9]则通过建立失稳雷诺数求解关系式的方法进行模型简化。本文在团队开发的RANS方程求解程序上,对Langtry-Menter模型[4, 5, 6]进行简化研究,进而得到转捩模型的一方程形式,并利用Schubanuer and Klebanoff平板调试程序参数,对S809翼型及DLR-F5机翼进行了算例验证。
计算采用有限体积法求解RANS方程。无黏通量通过Roe的FDS格式离散,黏性通量采用中心差分格式进行离散,时间推进采用隐式近似因子分解法。本文涉及的湍流模型有SST[10]两方程湍流模型和SA[11]湍流模型。
1.1转捩模型
转捩模型包含2个输运方程来求解转捩临界动量厚度雷诺数和间歇因子。利用涡量雷诺数与动量厚度雷诺数的关系,进一步将流场中传统的转捩判据解除了对非当地变量的依赖而改造为“当地判据”,判断转捩现象的发生。最终与SST模型耦合完成转捩模拟。 的输运方程如下
Pθt为源项,其中的Fθt是边界层指示因子,使其在边界层外符合经验关系式,在边界层内 来源于边界层外的扩散与上游的对流。
间歇因子γ输运方程如下
式中，Flength和Fonset是源项的关键参数。Flength控制源项的强度,进而影响间歇因子的增长速度,控制转捩区的长度。Fonset是源项的开关,控制间歇因子的增长即转捩的起始。
对于转捩模型其余各参数的详细含义可参考文献[4]。
因具有2个输运方程可视为两方程转捩模型,仔细分析方程间的联系可以发现,输运方程的主要作用是将边界层外通过经验关系式得到的Reθt输运至边界层内得到 ,进而利用经验关系式得到Flength与失稳雷诺数Reθc。Flength直接作用于γ输运方程源项中控制转捩区长度。失稳雷诺数Reθc则与涡量雷诺数ReV构成转捩判据影响Fonset,控制转捩的起始。实际上,各参数间可能并不是单一的起作用而是互相影响耦合作用的,例如Flength与Reθc间不是相互独立,形式也不是固定的[12]。
通过以上的分析可以看出,倘若建立了适用于边界层内的经验关系式,就有望节省掉 输运方程而发展为简化的一方程转捩模型。
γ输运方程与之前形式一致
式中
源项中对Pγ进行了部分简化。Pγ中的关键参数Flength在两方程转捩模型中由经验关系式 Flength=F( )得到。而在本文一方程算法中,由于没有 输运方程而直接将Flength赋值为常数。通过针对SK平板的数值模拟与实验的对比,最终将Flength数值标定为21.0,ce2=50,ca2=0.06,σγ=1.0。Fonset的求解公式参考Menter的方法详见文献[9]。
在简化两方程转捩模型的过程中,除了对Flength的处理,最大的问题就在于建立经验关系式直接求解边界层内失稳雷诺数Reθc。Coder和Maughmer的一方程模型没有通过输运方程而直接建立经验关系式求解 ,进而得到Reθc以达到简化目的。本文采用的是Menter直接建立求解 Reθt经验关系式的做法,即无需中间变量 。在建立这个新的经验关系式之前要引出以下2个适用于边界层内的参数——TuL,λθL。
TuL在传统的自由来流湍流度的基础上引入了ωdw,dw是壁面距离,ω是湍流耗散率(specific turbulence dissipation rate)。ωdw用以表征边界层内的速度尺度替代原有模型湍流度公式中的来流速度U,这就构造出适用于边界层内的湍流度公式。
在原转捩模型中,压力梯度因子求解公式为
该求法中 
是边界层外沿流线方向的速率导数,引入适用于边界层内垂直物面方向速度梯度 
。同样θ代表边界层动量厚度,若要建立适用于边界层内部的压力梯度因子公式,直接用θ并不能满足转捩模型的要求,因此引入dw代替θ。最终考虑众多因素后将公式确定为
至此,适用于边界层内部的所需参数已经建立,下面给出失稳雷诺数的经验公式[9]
FPG(λθL)是通过Falkner-Skan速度型进行标定的经验关系式,为了避免FPG出现负值,在程序中为其加限制FPG=max(FPG,0)[9]。
CPG1、CPG2、CPG3、CPG1lim、CPG2lim、CTU1、CTU2、CTU3均为模型常数,通过改变它们的数值用于程序调试。在本文中,根据SK平板的数值模拟与试验数据对比进行标定,最终确定的各参数如下
至此,一方程转捩模型就构建完毕。相较于转捩模型,现有的一方程模型并没有对分离流转捩进行特殊处理,可能会对模拟由分离泡引发的转捩过程产生不利影响。因此,本文在构造一方程转捩模型时直接引入了分离流转捩的修正方程。
为了模拟分离流转捩,本文在建立的一方程基础上引入了γ修正公式[4],最终得到γeff。
关于公式及参数的详细表达可参照文献[4]。 1.4 与湍流模型的耦合 1.4.1 与SST湍流模型的耦合
利用修正后的γeff与SST湍流模型进行结合。
Pk和Dk是原来湍流模型的源项和破坏项,
和 
是与转捩模型耦合过后的源项与破坏项。将此模型命名为One SST。针对Schubauer and Klebanoff平板进行参数标定时,便利用了此模型。
实现与SST的耦合并进行参数标定后,继续将该一方程转捩模型拓展到SA湍流模型。
在结合之前,给出所需的耗散率ω和湍动能k的修正公式[8],具体可参考文献[10, 13]。
与SA湍流模型结合,SA模型公式如下
为初始的生成项和破坏项,其中ft2=ct3exp(-ct4χ2),详细参数可参考文献[11]。参考文献[10]方法推导分析SST与SA湍流模型输运方程之间的联系,可以发现间歇因子对源项的影响类似于SA模型中ft2函数对源项的影响。通过分析原间歇因子对湍流模型的作用效果,本文将ft2函数加以改进,完成与SA湍流模型的耦合。
改进后的ft2函数为
在层流区域时,f′t2与原函数作用相似,转捩发生时,通过f′t2函数进而影响湍流模型中源项。下文中将此一方程转捩模型简称为One SA。 2 算例验证与结果分析
本文选取Schubauer and Klebanoff平板进行One SST模型参数标定。并针对S809翼型以及DLR-F5机翼,计算验证One SST、One SA模型对转捩位置及转捩过程的模拟精度。同时对于S809翼型还与Langtry-Menter转捩模型的计算数据进行了对比分析。
2.1 Schubanuer and Klebanoff平板Schubanuer and Klebanoff平板无厚度,长度2.0 m。速度入口边界位于平板前缘,用一段长度0.015 m的对称边界和平板相连,压力出口直接与平板相连。物面布置336个网格节点,法向第一层网格高度5×10-5m,增长率1.1,满足计算时y+小于1。
计算状态参数为:Tu∞=0.18％,u∞=50.1 m／s,ρ=1.225 kg／m3,μ∞=1.86·10-5
本文通过计算该平板调试One SST模型的各参数,确定参数后将一方程转捩模型拓展耦合到SA湍流模型。
调试时对SK平板转捩位置及转捩过程的模拟与实验数据如图 1所示。
 |
| 图 1 SK平板摩擦阻力系数与实验值对比 |
S809翼型是一个为风力机设计的相对厚度21％的层流翼型。风洞实验在Delft大学的低湍流度风洞中完成。实验雷诺数为2×106,马赫数为0.1,湍流度0.2％,湍流黏性比10。本文分别使用了简化的One SST、One SA模型及转捩模型进行计算,并与实验数据进行了对比分析。
S809翼型在1°迎角时分别采用上述2种一方程转捩模型进行计算,结果显示One SST和One SA转捩模型均能捕捉到翼型的层流分离及湍流再附现象。图 2为两模型计算得到的压力分布与实验值的对比。由图可以看出,分离导致了约52％弦长处小压力平台的出现,本文构建的2种分别耦合了SST与SA湍流模型的一方程转捩模型均捕捉到了这一现象。
 |
| 图 2 1°迎角时的压力分布 |
图 3和图 4描述的是S809翼型0~20°迎角范围内,利用模型与简化后2种转捩模型的计算结果与实验值的对比。图 3是升力系数曲线,图 4是阻力系数曲线。图 5是转捩位置的对比。从升、阻力系数图中均可以看出,在线性段One SST与One SA的计算结果与实验数据吻合很好,8°迎角内误差基本控制在3％以内;而在失速段,只能模拟出随迎角变化的力系数增长趋势,在数值上均与实验结果有较大误差。
 |
| 图 3不同模型升力系数与实验对比结果 |
 |
| 图 4 不同模型阻力系数与实验的对比结果 |
从图 5可以看出,简化的一方程转捩模型虽然减少了 输运方程,但其计算效果并不逊色于两方程转捩模型,较准确地捕捉到了转捩位置。然而,翼型上翼面转捩位置在5°到8°迎角间变化剧烈,One SST、One SA 以及模型均没有实现准确模拟的目的,需要进一步的分析与研究。但通过修正ft2函数耦合SA湍流模型的方法得到了证明,后期可对其继续研究。
 |
| 图 5 转捩位置预测结果与实验对比图 |
DLR-F5机翼前缘后掠角为20°,翼型采用超临界对称翼型。风洞实验由Sobieczky在1994年完成,实验中将机翼直接安装在实验段侧壁上,马赫数为0.82,迎角2°,基于参考弦长0.15 m的雷诺数为1.5×106,参考面积0.16 m2。计算网格量为600万。网格附面层第一层高度保证y+小于1,计算状态按实验条件设置,取来流湍流度为0.5％,黏性比为10。
使用本文中转捩模型One SST及One SA进行计算,实验结果与数值模拟结果如图 6所示。
 |
| 图 6 实验结果与计算结果转捩位置对比 |
该机翼表面存在激波,正是激波的逆压梯度触发边界层转捩的发生。由图 6可见,本文2种转捩模型的计算结果在中翼段和外翼段与实验吻合良好,转捩位置较接近。2种转捩模型在翼根处区别较明显,但其翼根处转捩是由横流不稳定导致的,两者尚不具备预测横流转捩的能力,故均未能正确模拟。将横流转捩预测引入本文的转捩模型可作为下一步的研究计划。
3 结 论本文的主要研究结论有:
1) 在建立一方程转捩模型的过程中,求解失稳雷诺数关系式时,引入了边界层内湍流度和压力梯度因子的计算公式,使得适用于求解边界层内失稳雷诺数的经验关系式顺利建立。这也是一方程转捩模型在简化过程中可以省去 输运方程的关键所在。
2) 在团队开发的数值求解程序上将模型简化,建立了一方程转捩模型,选取Schubanuer and Klebanoff平板进行参数标定。实现了其与2种湍流模型的耦合,并进行算例验证。证明该模型在仅有一个输运方程的基础上实现了对边界层转捩现象的判定。
3) 在得到的一方程转捩模型基础上,为模拟分离诱导转捩现象,引入γeff进行修正,使该模型得以正确模拟分离引起的转捩现象。
4) 在引入SA湍流模型时,未直接类比耦合SST的做法,而是通过分析SST与SA输运方程之间的联系,对ft2函数进行修正来模拟间歇因子对源项的影响进行耦合。
本文构建的一方程转捩模型经算例验证与实验数据对比后发现吻合较好,与模型类似。但还是存在一些问题需要进一步的分析与研究。
1) 该转捩模型在对S809翼型转捩位置的模拟中,对于变化剧烈的个别迎角没有准确模拟。此外,对于失速段的力系数的模拟只能模拟变化趋势,目前包括转捩模型以及使用SST,SA等全湍计算均无法正确模拟数值大小。需要更为深入的研究。
2) 模型本身无法模拟横流转捩,所以在其基础上进行简化得到的一方程模型也不能捕捉横流转捩现象。在后续的研究中可以尝试将横流转捩引入其中,使其具备该预测能力。
| [1] | Smith A M O, Gamberoni N. Transition, Pressure Gradient and Stability Theory[R]. Douglas Aircraft Company, Long Beach, Calif Rep No. ES 263881956 |
| Click to display the text | |
| [2] | Van Ingen J L. A Suggested Semi-Empirical Method for the Calculation of the Boundary Layer Transition Region[R]. Univ of Delft, Dept Aerospace Engineering, Delft, The Netherlands, Rep VTH-74,1956 |
| Click to display the text | |
| [3] | Andreas Krumbein. eN Transition Prediction for 3D Wing Configurations Using Database Methods and a Local, Linear Stability Code[J]. Aerospace Science and Technology, 2008(12): 592-598 |
| Click to display the text | |
| [4] | Langtry R B, Menter F R. Correlation-Based Transition Modeling for Unstructured Parallelized Computational Fluid Dynamics Codes[J]. AIAA Journal, 2009, 47(12): 2894-2906 |
| Click to display the text | |
| [5] | Menter F R, Langtry R B, Likki S R, et al. A Correlation-Based Transition Model Using Local Variables-PartⅠ: Model Formulation[J]. Journal of Turbomachinery, 2006, 128(3): 413-422 |
| Click to display the text | |
| [6] | Langtry R B, Menter F R, Likki S R, et al. A Correlation-Based Transition Model Using Local Variables-PartⅡ: Test Cases and Industrial Applications[J]. Journal of Turbomachinery, 2006, 128(3): 423-434 |
| Click to display the text | |
| [7] | Krause M, Barek M, Ballmann J. Modeling of Transition Effects in Hypersonic Intake Flows Using a Correlation-Based Intermittency Model[R]. AIAA-2008-2598 |
| Click to display the text | |
| [8] | Coder J G, Maughmer M D. One-Equation Transition Closure for Eddy-Viscosityturbulence Models in CFD[R]. AIAA-2012-0672 |
| Click to display the text | |
| [9] | Menter F R,Smirnov P E, Liu Tao. Avancha Ravikanth. A One-Equation Local Correlation-Based Transition Model[J]. J Flow Turbulence Combust, 2015, 95(4): 1-37 |
| [10] | Menter F R. Two-Equation Eddy-Viscosity Turbulence Models for Enginnering Application[J]. AIAA Journal, 1994, 32(8): 1598-1605 |
| Click to display the text | |
| [11] | Spalart P R, Allmaras S R. A One-Equation Turbulence Model for Aerodynamic Flows[J]. Recherche Aerospatiale, 1994(1): 5-21 |
| Click to display the text | |
| [12] |
张玉伦, 王光学, 孟德虹,等. 转捩模型的标定研究[J]. 空气动力学学报,2011, 29(3): 295-301 Zhang Yulun, Wang Guangxue, Meng Dehong, et al. Calibration of Transition Model[J]. Acta Aerodynamica Sinica, 2011, 29(3): 295-301 (in Chinese)
|
| Cited By in Cnki (15) s | |
| [13] | Menter F R. Eddy Viscosity Transport Equations and Their Relation to the k-ε Model[J]. Journal of Fluid Engineering, 1997, 119(4): 876-884 |
| Click to display the text |
