刘时杰, 徐浩军, 薛源, 张久星
(空军工程大学 航空航天工程学院, 陕西 西安 710038)
微下击暴流属于下冲气流的一种,是以垂直风切变为主要特征的复杂风切变区域[1-2]。研究微下击暴流对于飞机起飞、着陆时的飞行品质分析与飞行安全保障具有重要的意义。微下击暴流模型的研究目前有多种方法,例如利用基于风洞实验与CFD数值计算的方式研究其速度场分布[3-4],以及基于融合重叠的方法[5]和基于涡环原理的方法[6-8]。其中Woodfield和Woods[6]最早提出的涡环原理模型是目前应用最广泛的一种微下击暴流模型。
目前的方法所构建的微下击暴流模型多为沿中心轴对称的规则形状,但实际的微下击暴流风场非常复杂,其下沉气流并不一定垂直于地面,其外流范围也不是轴对称的,且一个区域内可能伴随有多个风暴核,甚至还有局部上升气流。并且在其风速的分布范围内亦可能伴随有其他大气流动状况,如不规则突风与飞机尾流产生的湍流。故本文拟在构建不规则微下击暴流模型的基础上考虑突风与湍流的影响,为起飞、着陆与低空飞行时的安全性仿真提供更加准确的风场输入条件。
首先利用涡环原理构造微下击暴流风场模型。依据美式坐标系,建立一个如图1所示的涡环模型。在地面处,风速垂直方向的分量为0。通过在地面上方布置一个强度为Γ的主涡环,同时在对称下方布置一个强度为-Γ的镜像涡环,即可满足地面垂直速度为0的边界条件。图1中,主涡环中心在P(xP,yP,zP)处,镜像涡环中心在(xP,yP,-zP)处。,R为涡环半径,A为空间中某点,坐标为(xA,yA,zA)。
图1 单涡环示意图Fig.1 Single vortex ring
首先讨论单涡环影响下的微下击暴流场构建方法,当涡环面与地面平行时,涡环的流线方程为:
ψ=-(Γ/2π)(r1+r2)[F1(k)-F2(k)]
(1)
式中,Γ为涡环强度;r1为A点距涡环距离最近的距离;r2为A点距涡环最远点的距离;k=|(r2-r1)/(r2+r1)|;F1(k)与F2(k)为圆环的积分函数。若0≤k2≤1,则:
(2)
将式(2)代入式(1),获得主涡环的流线方程为:
(3)
镜像涡环的流线表达式为:
(4)
最终将主涡环流线表达式与镜像涡环流线表达式线性叠加,得到空间点A的流线方程为:
(5)
而后得到点A的诱导速度为:
(6)
(7)
wz=(-1/rA)∂ψ/∂rA
(8)
式中,rA为点A距涡环对称轴的距离。
在涡环的中心轴处rA=0,用式(6)~式(8)计算诱导速度会产生奇点。因此,通过引入涡环的位函数,利用式(10)计算垂直方向的速度wz;由于涡环沿中心轴的对称性,水平速度wx和wy为0。
wz(z)=(Γ/2R)[1/(1+(z/R)2)3/2]
(9)
式中,z为中心轴上各点离涡环中心的距离。
根据式(6)~式(8),在涡环的中心涡丝处,流速趋于无穷大,这不符合微下击暴流风场的实际情况。由于流体粘性的影响,实际涡环有一个涡核的存在,且涡核内速率逐渐减小,至涡核的中心涡丝处减为0。本文运用复合涡原理,将涡核看作半径为r的圆柱,涡核内部的涡量均匀分布,保证涡丝处的流速为0,而涡核外部仍然服从流线方程。从涡核中心到涡核半径处,流速呈线性分布,如图2所示。
图2 涡核内部诱导速度示意图Fig.2 Velocity induced inside the vortex core
设涡核半径为r,首先应该判断点A(xA,yA,zA)是否在涡核内,若满足下式,则点A在涡核内。
(rA-R)2+(zP-zA)2≤r2
(10)
根据点A的坐标、涡环中心的坐标求出与点A、涡环中心点共面的涡核丝坐标点O(xO,yO,zO)。再求出与点A、点O共线,与点A、涡核中心点共面的涡核外壁坐标N,如图2所示。根据式(6)~式(8)求出坐标N的诱导速度,再根据点A、点N距离点O的距离按比例求出A点的诱导速度。
根据美式坐标系确定旋转轴心为涡环的中心点处(见图1),引入俯仰变换矩阵Lθ、滚转变换矩阵Lφ、偏航变换矩阵Lψ,对坐标系及速度进行旋转变换,从而构建不规则速度场。对于主涡环坐标系的变换公式为式(14),对于速度矢量的变换公式为式(15)。镜像涡环的俯仰变换角度和滚转变换角度与主涡环正负相反,偏航变换角度与主涡环相同。
(11)
(12)
(13)
(14)
(15)
经过变换后的坐标排列亦是不规则的,因为飞行实时仿真需要根据微下击暴流模型动态插值每个坐标点的速度矢量,故这种混乱不规则的坐标排列在飞行实时仿真中是无法使用的。故需将变换后坐标点对应的速度矢量拟合到标准的升序排列三维坐标点上,如式(16)。文中采用的拟合方法为线性插值法,由于坐标点较多,达到了100万以上,故拟合计算的过程相对于其他步骤来说耗时较长。
(16)
将标准坐标系内的主涡环与镜像涡环的诱导速度场相叠加,得到涡环1总的诱导速度场[wx1,wy1,wz1]。
由于单个涡环诱导的微下击暴流场垂直风速相对较小,为了更加逼真,通常使用多个涡环诱导微下击暴流场。多个倾斜涡环叠加而成的微下击暴流速度场为:
(17)
式中,[wxn,wyn,wzn]为第n个涡环的诱导速度场。
湍流风场采用Dryden模型[9]得到,其单侧一维频谱函数为:
(18)
式中,Lwx,Lwy,Lwz为湍流尺度;σwx,σwy,σwz为湍流强度。文中Lwx=2Lwy=2Lwz=380 m,σwx=σwy=σwz=3.6 m/s。
在飞行实时仿真中,关键问题在于如何产生高保真度的大气湍流。产生的湍流速度不仅要具有随机性,还必须符合其频谱特性。本文利用白噪声通过一定形式的滤波器来产生湍流速度。白噪声X(t)是一种具有特殊性质的随机过程,它的频谱函数等于常数,其相关函数为脉冲函数,如下式:
ΦXX(ω)=N0,RXX(τ)=N0δ(τ)
让白噪声X(t)通过一个传递函数为G(s)的环节,设其输出为Y(t)。输出Y(t)的频谱函数为:
ΦYY(ω) =G*(iω)ΦXX(ω)G(iω)
=N0G*(iω)G(iω)
由上式可见,若要求输出Y(t)具有特定的频谱ΦYY(ω),应将ΦYY(ω)作因式分解为G*(iω)和G(iω)的乘积,从而找到滤波器应具有的传递函数G(s)。
以Dryden模型(式(18))的第一式为例进行因式分解。将其按空间频率和时间频率的转换关系Φ(ω)=(1/v*)Φ(ω/v*)转化为:
对其进行因式分解:
Φwxwx(ω)=
所以其对应的传递函数为:
同理可得Gwy(s)和Gwz(s)。
根据以上论述,三个方向的大气湍流速度的产生如图3所示。
图3 湍流速度的产生Fig.3 Generation of turbulence velocity
将微下击暴流速度场与尾流引起的湍流速度场在坐标系内进行线性叠加,从而得到受湍流效应影响的微下击暴流速度场。
选择三个涡环来构成不规则微下击暴流风场,坐标范围为x:0~6 km,y:0~6 km,z:0~0.8 km。第1个涡环的中心点坐标为(3.0,3.0,0.6)km,R=0.9 km,r=0.455 km,Γ=18 000 m2/s,俯仰偏角为15°,滚转偏角为5°;第2个涡环的中心点坐标为(2.5,2.5,0.5) km,R=0.7 km,r=0.3 km,Γ=-10 000 m2/s,俯仰偏角为10°,滚转偏角为-5°;第3个涡环的中心点坐标为(3.0,4.0,0.6) km,R=0.8 km,r=0.4 km,Γ=11 000 m2/s,俯仰偏角为25°,滚转偏角为15°。所得结果如图4~图9所示。
从图4~图8中可以看出,相比于规则风场,文中对不规则微下击暴流的构建是成功有效的,构建的风场不再沿轴线对称,比较符合实际情况;并且由于加入了湍流的影响,速度场存在明显的不规则湍动,较好地反映了低空大气湍流对风场随机性的影响。
图4 在剖面y=3 km上的风速矢量图Fig.4 Wind speed vector on the profile y=3 km
图5 在剖面z=0.4 km上的风速矢量图Fig.5 Wind speed vector on the profile z=0.4 km
图6 在剖面y=3 km上的wz云图Fig.6 Cloud map of wz on the profile y=3 km
图7 在剖面z=0.4 km上的wx云图Fig.7 Cloud map of wx on the profile z=0.4 km
图8 在剖面x=3.5 km上的wz云图Fig.8 Cloud map of wz on the profile x=3.5 km
图9 三维空间内wz=-5 m/s的拓扑结构图Fig.9 Topology structure in the 3D space with wz=-5 m/s
本文综合考虑了微下击暴流的随机性与不确定性,通过多涡环建模法构建了具有旋转倾斜效应的不规则微下击暴流速度场,并将其与低空大气湍流相互叠加以反映风场在湍流作用下的随机性。所得风场能较好地描述多因素耦合复杂情况下的大气流动状况,相比于以往对微下击暴流速度场的研究更加深入、系统,也更接近复杂风场的实际情况。对于起飞、着陆阶段飞行安全性仿真分析、复杂风场中的风险规避以及控制技术研究具有积极的意义,亦可以进一步补充完善大气流动数据库。
参考文献:
[1] Fujita T T.Tornadoes and downbursts in the context of generalized planetary scales [J].Journal of Atmospheric Science,1981,38(8):1511-1534.
[2] 朱上翔.微下击暴流的流体动力学模型[J].飞行力学,1984,2(2):59-72.
[3] Qu W,Ji B.Numerical study on formation and diffusion wind fields for thunderstorm microburst[C]//Proc.of the International Conference on Mechanic Automation and Control Engineering.U.S.:IEEE,2010:1389-1392.
[4] 瞿伟廉,王锦文.下击暴流风荷载的数值模拟[J].武汉理工大学学报,2008,30(2):70-74.
[5] Nguyen H,Manuel L,Veers P.Simulation of inflow velocity fields and wind turbine loads during thunderstorm downbursts[R].AIAA-2010-2651,2010.
[6] Woodfield A A,Woods J F.Worldwide experice of wind shear during 1981-1982[R].ADP002705,1983.
[7] DoD.MIL-HDBK-1797 Military Handbook[S].U.S.:DoD,1997.
[8] DoD.MIL-F-8785C Military Specification[S].U.S.:DoD,1980.
[9] Yeager J.Implementation and testing of turbulence models for the F18-HARV simulation[R].NASA CR-1998-206937,1998.