基于k-两方程模式的湍流绕流流场数值模拟

基于k-两方程模式的湍流绕流流场数值模拟

1大攻角特性分析现代战争对战术导弹的射程和机动性提出了很高的要求。为了适应这种要求,先进战术导弹多采用小升力面甚至无翼布局,以减小巡航段的零升阻力,在需要过载时提高飞行攻角,充分利用大攻角时弹体的非线性涡升力输出所需的过载。飞行攻角增大到一定程度,由于三维边界层的分离,细长体背风面将形成分离涡,卷起形成特定的空间涡结构。背风面分离流动的特性对飞行器的气动性能有着重要的影响,是影响飞行器大攻角稳定性和操纵性的关键因素,给飞行器气动外形和控制系统的设计带来许多新问题。因此尽可能准确地预测前体大攻角特性具有很大的工程应用价值。分离涡是三维边界层分离形成的,这一物理本质决定了对细长体背风面分离涡的准确模拟必须考虑粘性。此外,真实飞行器飞行雷诺数很高,可达到107～109量级,其绕流绝大部分呈湍流流态。因此细长体背风面分离流动模拟的准确度在很大程度上又与所采用的湍流模式有关。近几年来,国内外学者在湍流模式研究中发展了一系列的非线性涡粘性模式,对研究湍流的近壁特性和流场曲率影响效果显著。为探讨不同湍流模式对细长体分离流的模拟能力,作者曾选用不同湍流模式,针对尖头细长旋成体,研究过超音速中等大攻角下的背风面分离流。为了深入分析细长体绕流的物理特征以及湍流模式的适用性,本文对前阶段工作做了进一步的深化和细化,选用Shih根据可实现性条件提出的k-ε模式对不同速度下旋成体背风面分离流进行了数值模拟,从背风面分离涡的强度和位置、物面压力分布、集中力和力矩等方面与已有试验结果进行了对比,并分析了分离涡内的湍流特性。2数值模拟方法2.1u3000u3000e+对密度和压力采用Renolds平均,速度、能量和温度采用Favre平均,结合Boussinesque涡粘性假设,可压缩N-S方程可以写为∂ρ∂t+∂ρuj∂xj=0∂ρui∂t+∂ρujui∂xj=-∂p∂xi+∂σij∂xj∂ρE∂t+∂ρujΗ∂xj=∂∂xj(uiσij-qj)}(1)∂ρ∂t+∂ρuj∂xj=0∂ρui∂t+∂ρujui∂xj=−∂p∂xi+∂σij∂xj∂ρE∂t+∂ρujH∂xj=∂∂xj(uiσij−qj)⎫⎭⎬⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪(1)其中σij=(μ+μt)(2Sij-23Sllδij)-23ρkδijqj=-γ(μΡr+μtΡrt)∂e∂xjE=e+12uiui+kΗ=h+12uiui+k}(2)σij=(μ+μt)(2Sij−23Sllδij)−23ρkδijqj=−γ(μPr+μtPrt)∂e∂xjE=e+12uiui+kH=h+12uiui+k⎫⎭⎬⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪(2)式中e=CvT,h=CpT,p=(γ-1)ρe。用有限体积法对方程(1)进行空间离散,多步Runge-Kutta显式方法进行时间推进。2.2可实现性条件是构造通用湍流模式的一个极端重要和有效的约束条件,要求模拟后的量不应当产生物理上不可能的值,如负的正应力或湍动能,关联系数大于1等。Shih等人根据可实现性条件提出的k-ε模式为∂∂t(ρk)+∂∂xi(ρkui)=∂∂xi[(μ+μtσk)∂k∂xi]+Ρk-ρε-YΜ(3)∂∂t(ρk)+∂∂xi(ρkui)=∂∂xi[(μ+μtσk)∂k∂xi]+Pk−ρε−YM(3)∂∂t(ρε)+∂∂xi(ρεui)=∂∂xi[(μ+μtσε)∂ε∂xi]+ρC1Sε-ρC2ε2k+√νε(4)∂∂t(ρε)+∂∂xi(ρεui)=∂∂xi[(μ+μtσε)∂ε∂xi]+ρC1Sε−ρC2ε2k+νε√(4)式中C2是常数,σk和σε是对应于k和ε的Prandtl数,Pk和YM分别是湍动能生成项和压缩性影响项。湍流粘性定义为μt=ρCμk2εμt=ρCμk2ε其中Ρk=-ρ¯u´iu´jSij‚YΜ=2ρεΜ2tC1=max(0.43‚ηη+5)‚η=SkεC2=1.9‚σk=1.0‚σε=1.2Cμ=1A0+ASU(*)k/ε2.3k-模式下的k-模式定义湍流雷诺数ReyRey=ρ√ky/μ根据Rey的大小将壁面附近的计算域分为两层:Rey>200时认为是完全湍流区,采用两方程k-ε模式;而当Rey<200时,采用Wolfstein提出的近壁处理方法,只求解k方程,分别用关系式μt=ρCμ√klμ‚ε=k32/lε计算湍流粘性系数和湍流耗散率ε,其中的长度尺度取为lμ=C1y[1-e(-ReyAμ)]‚lε=C1y[1-e(-ReyAε)]式中的常数取为C1=κC-3/4μ,Au=70,Aε=2C13流场的对称面文献提供了一个尖拱形旋成体的试验结果。试验模型是13倍长细比的尖拱形旋成体,头部长细比等于3,柱段长细比等于10。本文选用其中两种状态进行数值模拟:(1)M=2.5,α=14°,Red=1.23×106(2)M=0.7,α=14°,Red=0.667×106利用流场的对称性特点,计算采用半模。O型结构网格,周向70点,径向65点,轴向131点,总共约60万网格单元。周向网格在侧面和背面进行了加密;径向网格在靠近物面处进行了加密,最靠近物面的一层网格间距等于0.0001d。物面按绝热边界处理,对称面上使用对称边界条件,计算域远场采用无反射边界条件。来流湍流度取0.5%,湍流粘性比等于1。4结果分析4.1次分离锚和极限流线图1给出两个计算马赫数下三个不同轴向位置背风面截面流线。在X=3.5d截面,两个马赫数均出现了小范围的分离。相比之下,M2.5的分离区更大,涡核位置也更靠外。在X=6.5d截面,分离区同时向高度方向和展向扩张,形成很强的主分离涡。M2.5情况还出现了清晰的二次分离涡。X=11.5d截面,主分离涡展向位置稍有外移,但高度方向上移动很大。这时M2.5的二次分离涡仍然很强、很清晰,M0.7也分辨出了二次分离。整体来看,三个截面上M2.5的主分离点都更靠前,分离区也明显比M0.7的大。图2是物面极限流线。M2.5的分离区明显大于M0.7,主分离线更靠近迎风侧。M2.5的极限流线中还可以分辨出三次分离线,而M0.7情况则只能看出二次分离。此外,背风面分离流的空间流线表明,在这种中等大攻角下,分离涡飘起的高度不大,其空间走向基本上能保持与弹身轴线平行。4.2压心长度偏差表1列出了M=2.5条件下力和力矩系数结果。同试验结果相比,法向力偏低约2%,压心稍微靠后,相对弹身长度的误差为0.5%。阻力系数偏大2.5%。4.3相对总压的变化曲线图3是X=11.5d截面上沿穿过涡核水平方向(z方向)和垂直方向(y方向)上的相对总压的变化曲线。总的来看,数值计算很好地抓住了相对总压的最小值及其位置,意味着很好地模拟了主涡的强度及涡核位置。表2给出了X=11.5d截面上主涡涡核位置。计算结果与试验结果吻合得很好。相比之下,M2.5的涡核位置更靠外、更高。4.4x=6.5截面的分离锚三个不同横截面上物面压力分布如图4所示。X=3.5d截面由于分离区小、涡的强度弱,分离对压力分布的影响有限,背风面压力单调升高。M0.7的分离涡位置偏后,约在160°,而M2.5的分离涡位置稍靠前,在150°以前。X=6.5d截面该截面分离涡的影响较为明显、复杂。M2.5的压力曲线清楚地反映出了主分离涡和二次分离涡造成的吸力峰,分别位于155°和130°。M0.7分离涡影响区域明显较小,位置也靠后(165°)。试验数据中主分离涡的吸力峰值比计算结果低很多,在135°位置存在一个很弱的压力凹坑,预示着该位置可能存在一个二次涡。这两个迹象表明在该位置数值模拟计算出的涡强度偏弱。X=11.5d截面M2.5的分离区进一步扩大,受分离涡的影响,整个背风面压力基本不变,分辨不出主分离涡和二次分离涡。M0.7主分离涡位置仍在162°附近,且形成了一个压力低谷。4.5边界层纤相变化的纵向内湍流特性,主要分为5.从三个横截面上湍流粘性比可以分析得出,沿弹身轴线方向,由于不断地有边界层涡量的输入,分离区湍流粘性比在增大。相比之下,M=0.7时湍流粘性比更大,几倍于M=2.5时的值,表明此时湍流发展更快、更充分。5数值结果与试验结果的比较本文用一种符合可实现性条件的k-ε模