箔片转动数学建模及仿真分析

付俊博,金忠庆

(1.吉林交通职业技术学院,吉林 长春 130000；2.空军航空大学,吉林 长春 130000)

1 引 言

箔片型面源红外诱饵由上千个箔片压缩而成,诱饵发射后箔片迅速扩散形成箔片云。箔片在运动过程中,只受重力和气动力的影响。由于箔片的气动中心与重心不重合,因此存在气动力矩,使得箔片在运动过程中还伴随着转动。箔片的转动使得其所受气动力呈周期性变化,影响箔片的运动过程,并决定着箔片云的运动扩散性能。但由于箔片的转动过程中存在着障碍转动运动进行的气动阻尼,使得箔片的转动过程较为复杂。

目前,国内外的学者对箔片的运动及转动进行了相应的建模研究。邹涛将箔片运动过程中的转动角速度认为是定值,并建立了箔片的运动模型[1]。黄蓓应用非结构动网格方法研究了箔片非定常分离过程,并设计了实验对数值计算结果进行验证,其仿真计算结果表明箔片在降落过程中存在滞空和转动运动[2]。黄蓓还应用高速摄像机设计了一个箔片运动试验,依靠火药燃烧产生的推力将箔片推出发射筒,试验结果表明箔片在高速运动中还存在着高速转动[3]。

由于箔片的转动过程较为复杂,涉及多学科、交叉理论。其内部既包含了数值积分算法、数学建模方法又涵盖了空气动力学、飞行力学、流体力学、理论力学、非定常流动等相关领域的知识[4-5]。同时箔片转动过程还需要考虑气动阻尼对箔片转动力矩的影响。因此,国内外的学者对箔片的运动模型研究较为充分,而对箔片的转动过程及转动模型的研究还不够深入,相关模型建立的较为简单,不能够真实反应箔片转动过程中其角速度的变化规律,进而不能够得到箔片云的真实扩散运动过程。

基于此,本文对箔片的转动现象进行深入研究,建立了相对精确的箔片转动模型,并将仿真结果与实验对比,以验证模型的准确性。然后,应用建立的模型对箔片的转动特性进行分析,并仿真分析了箔片云的扩散运动过程。

2 箔片运动建模

由于箔片的运动需要在地轴系内计算,而箔片的动力学方程需要在航迹坐标系内计算,因此首先需要建立箔片的地轴系Oxgygzg、航迹坐标系Oxhyhzh以及速度坐标系Oxayaza,各坐标系的定义见文献[6]～[7],坐标系之间的关系如图1所示。

图1 箔片各坐标系之间的关系

定义Oxh轴与水平面之间的夹角为航迹俯仰角θ,运动方向向上为正；Oya轴与Oyh轴之间的夹角为速度滚转角γs,速度坐标系绕Oxh轴向右旋转时,γs为正；Oxh轴水平面的投影与Oxgyg之间的夹角为航向角ψs,Oxh轴左偏为正。定义垂直于箔片平面且过圆心的单位向量为箔片轴心线向量nd,Oya轴与Oxa轴所组成的平面为气流对称平面。则向量nd及箔片运动方向都在气流对称平面内。

由于箔片在整个运动过程中都是无动力的,因此在箔片运动过程中只受气动力和重力的影响。由于箔片存在特殊的对称性,且其表面较为平滑,因此气动侧力予以忽略。则其运动过程中只受到气动阻力X、气动升力Y的作用,且其受到的气动力都在气流对称平面内,如图2所示。

图2 箔片受力分析

箔片的动力学方程可以简化为[8]:

(1)

式中,X、Y为气动力;cx、cy为相应的气动力系数;V为箔片的运动速度;γs为箔片的速度滚转角;ψs为箔片的航向角;θ为箔片航迹俯仰角;ρ为当前环境的大气密度。由于箔片存在特殊的对称性,因此其向各个方向运动时,cx、cy都只与箔片的迎角和运动速度有关。

箔片的运动学方程为[9-10]:

(2)

式中,(xg,yg,zg)为箔片在地轴系中的坐标。则:

(3)

式中,ndx、ndy、ndz为nd在地轴系中的三轴分量;nV为运动速度方向的单位向量;nVx、nVy、nVz为其在地轴系中的三轴分量。则:

(4)

式中,ψ为偏航角;ϑ为俯仰角。

定义nY为气动升力的单位向量,nYx、nYy、nYz为其在地轴系上的三轴分量。则:

(5)

因此可以得到速度滚转角γs、迎角α的表达式为:

(6)

(7)

3 箔片转动建模

箔片的运动过程中,由于其压力中心与箔片的几何中心不重合,因此在其运动过程中受到的合外力矩不为零。作用在其上的合外力矩使得箔片高速转动,箔片的轴心线向量nd、迎角α等参数周期性地进行变换,最终导致箔片所受的气动力不断变化。

3.1 箔片转动模型

在箔片的运动过程中,其所受的合外力矩主要由升力产生,力矩方程在速度坐标系内可表示为[11]:

(8)

式中,J为箔片转动惯量；R为箔片半径；xF为箔片压力中心与箔片几何中心之间的距离[12]；Md为气动阻尼力矩；ω为转动角速度；Ya为气动升力在速度坐标系中的值。箔片转动过程中所受的气动力和力矩示意图如图3所示。

图3 气动力和力矩分析

(9)

定义nYa为升力向量nY在速度坐标系中的值,则其值可以表示为:

(10)

nYay为nYa在Oya轴上的值,则Ya可以表示为:

Ya=|Y|nYay/|nYay|

(11)

3.2 气动阻尼的求解

箔片的转动过程中,同时还伴随着障碍转动运动进行的气动阻尼力矩Md。如果将箔片的运动过程分解为质心平移运动和绕质心转动运动,则箔片的质心平移运动可以由箔片的运动方程求得,而箔片绕质心的转动运动主要由气动升力产生的力矩M以及气动阻尼力矩Md求得。

设箔片的转动角速度为ω,dS为箔片表面到箔片中心距离为r的面源,如图4所示。

图4 箔片转动力矩分析

则转动过程中,dS所受的气动阻力为:

(12)

(13)

联立以上公式,则可以得到箔片的运动模型和转动模型。

由于箔片的角速度垂直于气流对称平面,气流对称平面又重合于速度坐标系的Oxaya平面,则Δt时刻后,轴心线向量nd在速度坐标系内可以表示为:

nda=[-sin(α+ωΔt),cos(α+ωΔt),0]

(14)

式中,nda为nd在速度坐标系中的值,因此nd可以表示为:

(15)

4 箔片转动实验与模型验证

4.1 箔片气动力系数计算

本节采用CFD对箔片的气动力系数进行数值计算。计算模型的选取如下:控制方程为Navier-Stokes方程,并采用PIM-PLE算法对其进行求解；空间离散方法为空间二阶精度线性插值算法和基于有限体积的空间离散算法；时间离散方法为二阶精度的后向差分算法。初始条件设置如下:箔片表面设置成无滑移壁面,计算域的出口、入口以及其他外壁全部设置成压力远场,湍流计算模型选取SST湍流模型[13]。

采用多面体网格对箔片的气动力系数进行计算,计算域为边长5000 mm的正方体,箔片位于计算域的正中心。箔片直径为50 mm,厚度1 mm,边界层一共设置12层,第一层厚度为10-6m,箔片表面及边界层内的网格如图5所示。计算域内的网格数总数为213万,计算精度能够满足要求。则高度H=0 km、环境温度T=288.15 K、大气压力P=101325 Pa、入口来流速度为100 m/s时,不同迎角下箔片气动力系数CFD的计算结果如表1所示。

图5 多面体网格及箔片边界层网格示意图

表1 箔片气动力系数计算结果

4.2 箔片转动实验

由于箔片运动过程中其姿态、位置变化较为明显,且与其初始状态有关,实验中难以测得其运动过程中的角速度变化规律。因此本节设计了箔片在质心静止状态下的转动实验,箔片在实验中只存在转动过程质心不存在运动,并将实验结果与建立的转动模型计算结果进行对比。

箔片的直径为50 mm,厚度1 mm,尺寸与CFD数值计算的尺寸保持一致。在箔片表面固定有一个穿过其圆形紧贴箔片表面的转轴,转轴固定在一个低速风洞内,如图6所致。则当风洞工作时,箔片只能够绕转轴转动而不能够自由运动。箔片的右侧放置一个高速相机,距其0.5 m。箔片转动时,通过高速相机拍摄到的每一帧画面,计算出箔片的转动角度,进而得到其转动角速度。

图6 箔片转速测试试验示意图

将实验风速设置为20 m/s,为了与实验一致仿真中将箔片质心位置设置为静止不变,只考虑其转动问题不考虑运动问题,即认为公式(1)中的V=20、θ=0、ψs=0不变,公式(2)中的xg、yg、zg始终等于零。则初始迎角α=90°和α=10°时,实验结果和转动模型计算的结果对比图如图7所示。

图7 实验结果与仿真结果对比

由图7中的结果可知,箔片初始迎角α=90°时,箔片基本处于静止状态,几乎没有转动运动,其转动角度速度远低于α=10°。由此可知,箔片的转动角速度与其初始的角度关系很大。由于箔片在转动过程中其气动力矩周期性变化,因此箔片在转动过程中其角速度时刻变化,变化曲线类似于正弦曲线。而箔片在转动过程中气动阻尼力矩与转动方向相反,使得箔片的转动角速度逐渐较小。

从图7中的对比结果可知,转动模型仿真的角速度略大于实验结果,这主要与实验中箔片的转轴与轴承之间还存在摩擦力有关,仿真结果与实验结果总体一致,因此建立的转动模型较为可信,能够满足精度的要求。

5 仿真结果及现象分析

5.1 箔片转动仿真

设载机水平飞行,飞行高度为8 km,飞行马赫数Ma=0.8。箔片1-4从载机上垂直向上抛撒,抛撒速度为20 m/s,箔片直径均为50 mm。箔片1-4初始姿态角完全一致,初始俯仰角ϑ=45°,初始偏航角ψ=-30°。箔片1发射后依据文中建立的模型对其运动和转动过程进行求解,箔片2、箔片3、箔片4转动角速度不变,分别为50 rad/s、100 rad/s以及200 rad/s。以箔片抛撒时刻其位置在水平面上的投影为坐标原点,载机飞行方向为正X轴建立坐标系。比较四个箔片运动和转动过程,如图8、图9所示。

图8 箔片运动轨迹对比图

由图8的仿真结果可知,箔片的转动角速度对箔片的运动过程影响较大,箔片2、箔片3和箔片4相对于箔片1向Z轴运动的距离都较大,且箔片的转动角速度越小其向Z轴运动的距离越长。由于箔片的初始偏航角为负,因此箔片气动力在Z轴正方向具有一定的分量,导致箔片具有沿Z轴正方向运动的趋势。箔片的初始转速越小,其所受气动力周期性变化时间越长,最终导致箔片在Z轴上运动距离越长。

图9 箔片特性对比图

图9(a)的仿真结果可知,箔片抛撒后在气动阻力的作用下速度迅速降低,2.5 s后箔片几乎停止运动。图9(b)的仿真可知,箔片1的初始转动角速度较大,但随着速度的下降转动角速度迅速减小。由于抛撒时刻箔片的运动速度较大,因此其所受的初始转动力矩很大,导致箔片的初始转动速度较大。随着箔片运动速度的下降,箔片受到的转动力矩减小,最终导致箔片转动角速度越来越小。当箔片转动到α<0°时,气动力矩与转动角速度相反,在气动阻尼和气动力矩的作用下,箔片转动角速度迅速降低。因此从图9(b)还可以看出,箔片1的转动角速度呈高频大幅震荡状态,且振幅随着转动角速度的减小而下降。

5.2 箔片云运动扩散仿真

面源红外诱饵发射后由于大气扰动等随机因素的影响,箔片的初始姿态有一定的差别。经过与实测数据的对比分析,近似认为箔片初始俯仰角ϑ0服从U(0,π/2)的均匀分布,初始偏航角ψ0服从U(0,2π)的均匀分布。由于箔片的初始姿态已知,则根据文中建立的箔片运动模型和转动模型可以得到箔片云的整体运动扩散过程。

设载机水平飞行,飞行高度为8 km,飞行马赫数Ma=0.8。面源红外诱饵由一千个箔片压缩而成,箔片直径均为50 mm,诱饵发射初始速度为20 m/s,速度相对于载机垂直向上。则诱饵发射后,箔片云的运动扩散过程如图10、图11所示。

从图10、图11中的仿真结果可以看出,面源诱饵发射后,箔片云大致呈与X轴存在一定夹角的锥形分布。箔片云扩散速度较快,0.5 s时刻,箔片云在X轴扩散长度就达到了59 m,Y轴长度达37 m,Z轴扩散长度达40 m。1 s时刻,箔片云继续扩散,但扩散速度明显减慢。

图10 0.5 s时刻箔片云扩散图

图11 1 s时刻箔片云扩散图

6 结 论

本文研究了箔片在运动过程中的转动问题。通过定义箔片的轴心线向量,确定了箔片迎角和速度滚转角的求解方法,并建立了箔片的运动模型。然后,建立了箔片的转动模型,并研究了在转动过程中气动阻尼的求解方法。通过CFD流场计算方法,计算箔片的气动力系数,并对比箔片转速的实验结果和模型仿真结果,验证了模型的可信性。最终,仿真分析了箔片的转动特性及箔片云的运动扩散特性。所得到的主要结论如下:

(1)箔片的转动角速度与其初始迎角有关,初始迎角接近90°时,其转动角速度最小；

(2)箔片的转动角速度与运动速度有关,且随着运动速度的下降箔片转动角速度越来越小；

(3)箔片在运动过程中其转动角速度呈高频大幅震荡状态,且当来流风速较小时,其角速度变化曲线类似于正弦函数；

(4)箔片云扩散速度较快,诱饵发射1 s后箔片云已基本扩散成型,箔片云大致呈与载机飞行方向存在一定夹角的锥形分布,且箔片云在载机飞行方向上扩散能力最强。