湍流大涡破碎装置的减阻研究_陈强

第13卷第4期 弹　道　学　报 Vol.13No.4 2001年12月 Journal of Ballistics Dec.2001湍流大涡破碎装置的减阻研究＊

陈　强　唐登斌　曹起鹏

(南京航空航天大学空气动力学系,南京210016)

摘要　采用在边界层中引入控制装置的方法进行湍流减阻问题研究.着重研究了对

于正确模拟减阻过程至关重要的湍流模型,不同的模型(包括Johnson-King涡粘/雷

诺应力封闭模型和对近壁区的低雷诺数流动作修正的k-ε模型)用于NS方程的数值

模拟.计算从控制装置的上游开始,其结果更为准确可靠.为得到最佳减阻效果,分析

研究了装置的排列和参数的影响,给出了单个装置和并排双装置的典型减阻结果,并

和实验数据进行了比较.

关键词　减阻,湍流边界层,大涡破碎,湍流模型

中图分类号　V211,O357.5

飞行器(飞机、导弹等)的减阻技术研究是其设计工作的重要内容,减阻与提高飞行器的性能紧密相关.以湍流为主的表面摩擦阻力在飞行器总阻力中占有很大比重,因此,湍流减阻问题的研究具有特别重要的意义.通常,湍流减阻是通过改变或制约湍流结构和流动状态来达到减少表面摩阻的目的,其中采用大涡破碎(Large Eddy Break-Up,“LEBU”)装置的方法减阻和通过纵向表面微槽的方法减阻是目前各种湍流减阻方法中最为有效的[1-3].

采用LEBU减阻是把一种装置(这里为一很小的平板,简称“平片”)沿流向置入到边界层内,通过其切割作用使大涡破裂,改变湍流的生成模式,修正边界层的内层和外层的传输特性,以达到减少表面摩阻的目的.尽管湍流问题极为复杂,给机理研究带来很大困难,但研究结果表明[3],在装置后下游的一定范围内,平均的和脉动的速度型都有明显变化,从而会引起表面摩阻的变化.为避免许多计算是从装置后缘开始而带来的问题[4],本文是从装置上游开始进行全流场计算,并采用了不同的湍流模型.为得到最佳的减阻效果,对LEBU装置的各种不同情况进行的大量计算和详细分析将为该湍流减阻方法的实际应用提供基础.

1　控制方程

采用时间相关二维不可压Navier-Stokes方程作为控制方程,这里写成用涡量(ω)和流函数(ψ)表达的形式

ω/ t=-V·( ω)+v e2ω(1)

2ψ=-ω(2)

收稿日期:2001-09-28

＊　航空科学基金资助项目(98A52008)

式中

ω= v / x - u / y

(3)V =(u ,v )T =( ψ/ y ,- ψ/ x )T (4)

有效粘性系数νe =ν+νt ,ν为运动粘性系数,νt 为湍流涡粘性系数,是湍流脉动应力项-u ′v ′=νt

u / y 中的速度梯度的系数,将由湍流模型得到.边界条件(包括在所取矩形计算区域中的上游、下游边界,上边界以及主平面和平片的壁面边界):

(1)上游边界,给定流向速度型

u (y ),ω1,j =( u / y )1,j ,ψ1,j =

∫y 0u (y )d y (5.1)(2)下游边界

( ω/ x )I ,j =0,( ψ/ x )I ,j =0

(5.2)(3)上边界

ωi ,J =0,　ψi ,J =ψ1,J

(5.3)(4)主平面壁面

ψi ,1=0,　ωi ,1=-2ψi ,2/(Δy )2(5.4)

(5)平片壁面,u =0、v =0ψ=ψd ,ωd (u )=-2(ψd +1-ψd )/(Δy )2,ωd (l )=2(ψd -1-ψd )

/(Δy )2(5.5)式中,下标中(u )代表上表面;(l )代表下表面.

2　湍流模型

2.1　J -K 模型

该模型兼有涡粘性模型和雷诺应力模型的一些特性,通过调整所满足的关于最大雷诺应力沿流向变化的常微分方程,能反映一定的上游影响

[5],可写成如下形式νt =νt0〔1-ex p (-νt i /νt0)〕(6)νt i =D 2ky (-u ′v ′m )1/2(7)

νt0=σ(x )(0.0168u e δ＊i γ

)(8)式中:D =1-exp 〔1-u ′v ′m )1/2y /ν

A +〕;A +=15,γ为Klebanoff 间隙因子,γ=1.0/〔1+5.5(y /δ)6〕;σ(x )通过解下列最大雷诺应力的常微分方程确定

(-u ′v ′m )1/2=(-u ′v ′m ,eq )1/2-〔L m u m /a 1(-u ′v ′m )〕d -(u ′v ′m )/d x -L m D m /-(u ′v ′m )

(9)

这里的下标m 表示在该站上的-u ′v ′取最大值,a 1、L m 为模型参数,-u ′v ′m ,e q 为σ(x )

=1时所得到的最大雷诺应力,D m 值由下式给出

D m =C dif (-u ′v ′m )1/2|1-σ(x )1/2

|/a 1δ〔0.7-(y /δ)m 〕29第4期 陈　强等　湍流大涡破碎装置的减阻研究

2.2　k -ε模型

考虑到在近壁区的低雷诺数流动的影响,改进的k -ε两方程湍流模型

[6]可写成如下形

式

k / t + uk / x + vk / y = 〔(νt +ν) k / y 〕/ y +νt (

u / y )2-ε-2ν( k / y )2(10) ε/ t + u ε/ x + v ε/ y = 〔(νt +ν) ε/ y 〕/ y +C e1νt (

ε/k )( u / y )2-C e2ε2/k -2νt μ

( 2u / y 2)2(11)式中k 为湍动能,ε为湍能耗散率,νt =C μ

k 2/ε.有关系数可写成:C μ=0.09exp 〔2.5/(1+R T /50)〕;C e 1=1.44;C e 2=2〔1-0.3exp (-R 2T )〕;R T =k 2/νε.3　数值方法

在所取矩形计算区域中,其内部网格点随(i +j )为偶数或奇数分成二类,给出初始猜值,进行偶—奇循环迭代,包括二大步.

(1)辅助步(偶数网格点)ωn +1/2i ,j +1=(ωn i ,j +ωn i +1,j +1+ωn i -1,j +1+ωn i ,j +2)/4-(Δt /8Δx Δy )(ψn i ,j +2-ψn i ,j )

(ωn i +1,j +1-ωn i -1,j +1+(Δt /8Δy Δx )(ψn i +1,j +1-ψn i -1,j +1)(ωn i ,j +1-ωn i ,j )

-〔Δt /2(Δx 2+Δy 2)〕νe (ωn i -1,j +ωn i -1,j +2+ωn i +1,j +ωn i +1,j +2-4ωn i ,j +1)

(12)(ψn +1/2i -1,j +ψn +1/2i -1,j +2+ψn +1/2i +1,j +ψn +1/2i +1,j +2)-4ψn +1/2i ,j +1=(Δx 2+Δy 2)ωn +1/2i ,j +1

(13)(2)主步(奇数网格点)ωn +1i ,j =ωn i ,j -(Δt /4Δx Δy )(ψn +1/2i ,j +1-ψn +1/2i ,j -1)(ωn +1/2i +1,j -ωn +1/2i -1,j )

+(Δt /4Δy Δx )(ψn +1/2i +1,j -ψn +1/2i -1,j )(ωn +1/2i ,j +1-ωn +1/2i ,j -1)-〔Δt /(Δx 2+Δy 2)〕νe (ωn +1/2i +1,j +1+

ωn +1/2i -1,j +1+ωn +1/2i +1,j -1+ωn +1/2i -1,j -1-4ωn +1/2i ,j )

(14)(ψn +1i +1,j +1+ψn +1i -1,j +1+ψn +1i +1,j -1+ψn +1i -1,j -1)-4ψn +1i ,j =(Δx 2+Δy 2)ωn +1i ,j (15)这样,可以用高斯消去法求解流函数ψ的矩阵微分方程,进而得到流场的速度分布以及摩阻系数分布.

对于湍流模型方程的计算,采用有限差分方法求解k -ε偏微分方程(10)、(11).J -K 模型中的常微分方程(9)用稳式Euler 法求解,为解决该法需要知道上游开始处的雷诺应力分

布,可用简单的Cebeci -Smith 代数模型算出该处的u ′v ′m ,再由J -

K 模型向下游推进.图1　典型的LEBU 装置及其几何参数

4　结果和分析

分别研究了单平片和双平片二种典型装置的

减阻问题.为便于与实验数据比较,取来流参数:

u ∞=24m /s ,ν=1.46×10-5m 2/s ,δ0=0.

017m ,平片厚度不计.装置参数见图1,h 是从主平面

算起的平片高度,C 是平片弦长,S 为双平片间的

距离,δ0为平片前缘处的边界层厚度.30 弹　道　学　报 第13卷