基于计算流体力学与刚体动力学耦合的高速旋转(7)
[27] WEINACHT P. Projectile performance, stability, and free-flight motion prediction using computational fluid dynamics[J].Journal of Spacecraft and Rockets,2004, 45(2): 257-263.
[28] 徐辉雯,陈少松.二维修正弹修正组件反旋与不旋气动特性对比[J].弹道学报,2018,30(4):32-37.
XU H W, CHEN S of aerodynamic characteristics of two-dimension trajectory correction projectile with reverse-rotation and non-rotation PGK[J].Journal of Ballistics, 2018, 30(4):32-37.(in Chinese)
[29] 钟杨威.旋转稳定二维弹道修正弹弹道特性分析与制导方法研究[D].南京:南京理工大学,2015.
ZHONG Y W. Study on trajectory characteristics and guidance methods for spin stabilized two-dimensional trajectory correction projectiles[D]. Nanjing: Nanjing University of Science and Technology, 2015. (in Chinese)
0 引言先进的武器装备是国防建设的必要条件,射程、精度、威力等是武器系统的主要战术技术指标,武器系统的研制常从弹道设计[1]开始,因此弹道计算在整个武器系统研发过程中有着非常重要的地位。传统的弹道计算方法是利用基于气动力和力矩系数的弹道模型,通过龙格- 库塔(Runge-Kutta)法或阿当姆斯预报- 校正法进行数值求解。传统弹道计算方法需要建立弹箭气动力模型,对于旋成体弹丸,气动模型相对准确;而对于气动外形复杂甚至非对称的弹箭,建立准确的气动模型非常困难。因此,对于外形复杂的弹丸,如何准确获取其动态变化过程和弹道飞行特性,是弹道学所面临的一个难题。近十多年发展起来的计算流体力学(CFD)和刚体动力学(RBD)耦合计算方法,为弹道仿真提供了一种新的技术途径。相对于传统方法,该方法能够实时模拟弹箭在大气中真实的飞行情况,即使再复杂的弹箭,也能从实时变化的流场参数中获取气动力和力矩,在无需提供气动模型情况下完成弹道计算。另外,该方法能够实时提供弹箭飞行过程中弹体表面和周围空气的压力、密度、温度等全息流场参数,为研究人员提供丰富的数据资料,有利于新型弹箭的研制。文献[2]通过欧拉角建立弹体坐标系与地面坐标系的关系,在弹体坐标系下采用任意拉格朗日- 欧拉(ALE)控制方程对流场进行求解,首次将CFD/RBD耦合技术成功应用到弹道计算中。进而以某尾翼稳定弹为对象,采用阿当姆斯预报- 校正法对流场和6自由度刚体动力学方程组进行耦合求解,并通过试验结果验证了方法的有效性。该方法被进一步应用于超音速、跨音速和亚音速弹丸模拟中,通过弹道仿真结果和试验结果对比再一次验证了方法的正确性[3-5]。文献[6]在耦合CFD/RBD技术中引入入口边界条件,研究了喷流对弹道的影响规律,验证了在耦合计算中加入进口边界来模拟脉冲修正弹弹道的可行性。文献[7-11]通过弹丸飞行动力学模型推导出气动参数估计模型,根据耦合CFD/RBD计算结果估计出弹丸的气动系数。文献[12-16]在耦合CFD/RBD技术中加入飞行控制系统(FCS),建立耦合CFD/RBD/FCS方法并模拟了鸭舵受控尾翼弹的姿态运动,通过与实验数据对比验证了方法的正确性。文献[17]采用CFD/RBD和非结构动网格技术计算均匀流场下不同密度平板的运动轨迹,分析了碎片对飞机造成的二次伤害。文献[18]基于自主开发软件HUNS3D,结合非结构嵌套网格技术,使用阿当姆斯预估- 校正法耦合求解6自由度刚体运动方程,研究了子母弹分离过程。文献[19]采用刚性动网格和CFD技术,在规定运动下模拟了带鸭舵弹丸的运动过程,研究了锥形运动下弹丸非定常空气动力学特性。综上所述,文献[19]采用的是规定弹丸运动的方法,不能直接用于研究弹丸的自由运动。文献[17-18]这类方法的计算域通常包含物体运动轨迹,并不适用于研究射程达几十千米的弹丸。相比较而言,文献[2-16]基于弹体系建立流动关系更为合理。然而,对于高速旋转弹丸,基于弹体系建立流动关系也存在不足之处,高速旋转弹丸会使远场产生很大的附加速度,造成远场计算不准确[20]。现有的CFD/RBD耦合方法一般采用阿当姆斯预报- 校正法,该方法较容易实现,但不能自启动[1]。在流场初始耦合时只能采用当前值来代替过去值,这必将在初始阶段损失时间精度,且不易变步长。本文针对已有研究的不足,基于ALE形式的流动模型,建立控制体表面运动和弹轴坐标系运动的耦合模型。在此基础上提出一种基于4阶Runge-Kutta法的CFD/RBD耦合计算方法,为研究高速旋转弹丸的真实弹道和流场提供参考。1 CFD数值计算方法1.1 基于弹轴运动的ALE形式控制方程弹丸飞行时姿态不断变化,与弹丸固连的计算网格具有附加速度,ALE形式的Navier-Stokes方程(简称N-S方程)能够有效地求解具有网格速度的流场。由于弹丸高速旋转,在弹体系下建立流动关系时,远场流动速度相对于壁面附近非常大,从而导致计算不准确[20]。而弹轴坐标系有效去除了弹体高速自转,克服了这一问题。基于弹轴运动建立具有ALE形式的N-S方程并对流动进行求解,其积分形式为(1)式中:W为守恒变量;Ω为控制体体积;S为控制体表面的外法向矢量;源项Q=ρ(ωa×v),ρ为空气密度,ωa为弹轴系角速度,v为流动速度矢量;为黏性通量;为对流通量,(2)Fc为静态网格下的对流通量,v?Ω为附加在控制体表面的反变速度。当弹轴系的平动速度矢量和转动角速度矢量分别为vp和ωa时,附加在控制体表面的反变速度为v?Ω=n{vp+ωa×(r-rg)},(3)式中:n=[nx,ny,nz]T为控制体表面的单位外法向矢量,nx、ny、nz分别为单位外法向矢量在x轴、y轴、z轴上的分量;r为控制体表面微元的矢径;rg为弹丸质心矢?修正简单低耗散迎风矢通量分裂格式(1)式中对流通量的计算方法是CFD的研究热点之一。简单低耗散迎风矢通量分裂格式(SLAU)是著名迎风型矢通量分裂格式(AUSM)系列的一个发展方向,具有低耗散、格式简洁、低马赫数下无需调整参数等优点,其改进型SLAU2[21]在压力通量中加入了与来流马赫数相关的耗散项,继承了原格式优点,同时实现了比其他全速格式更简洁的形式。引入弹丸运动后,基于SLAU2的计算格式为(4)式中:F1/2为控制面上的数值计算通量;为质量通量;ψ+、ψ-分别为采用控制体表面左侧和右侧流动参数计算出的辅助向量ψ,ψ=[1,vx,vy,vz,H]T,vx、vy、vz分别为流动速度矢量v在x轴、y轴、z轴上的分量,H为单位质量气体总焓;为压力通量;N=[0,nx,ny,nz,v?Ω]T。(4)式中参数和的计算公式如下:(5)式中:ρl和ρr分别表示控制体表面左边和右边的空气密度;vn,l和vn,r分别为控制体表面左边和右边的总反变速度,vn,l=nvl-v?Ω,vn,r=nvr-v?Ω,vl、vr分别为控制面左右两侧流动速度矢量;和χ均为辅助变量;Δp为压力差,Δp=pr-pl,pl和pr分别为控制面左右两侧压力;cij为第i个控制体的第j个面的音速,该参数的常用数值计算方法有代数平均法、几何平均法、最大临界音速法和Roe平均法等[21]。文献[22]通过数值实验发现,采用代数平均法进行三维弹丸亚音速、跨音速和超音速外部绕流模拟时,计算结果具有较好的稳定性和收敛性,cij的代数平均法表达式为cij=(cl+cr)/2,(6)式中:cl和cr分别为控制面左右两侧音速。压力通量为(7)式中:和分别为控制面左右两侧压力相关系数;χm为辅助变量,它考虑了弹轴系的运动,(8)的表达式分别为(9)M为法向流动速度与音速比值,M=(v·n)/cij.1.3 边界条件、时间推进及湍流模型边界条件包含壁面边界条件和远场边界条件。其中,壁面采用绝热和无滑移条件,远场采用黎曼边界条件[20]。双时间步[17]可以有效地增大非定常流场计算时间步长即外层真实时间步长、提高效率,在工程上被广泛采用。本文内层伪时间采用具有总变差减小(TVD)性质的3阶Runge-Kutta法,结合当地时间步和残值光顺等加速收敛方法进行推进。由于网格附加了弹丸的运动,计算当地时间时无黏谱半径应考虑附加网格速度,当地时间表达式为(10)式中:σ为Corant-Friedrichs-Lewy(CFL)数;为系数,为控制体表面上的黏性谱半径。引入弹丸质心速度后,表达式为(11)式中:下标i表示控制体第i个面;nf为面的个数;vi为控制体第i个面上的流动速度矢量;ni为控制体第i个面的单位外法向矢量;ci和ΔSi分别为离散后控制体第i个面的音速和面积。计算(1)式中黏性通量的湍流黏性时采用Spalart-Allmaras(S-A)湍流模型,该模型具有计算量小、稳定性好、计算精度满足工程要求等优点,在气动力计算和航空航天等领域得到广泛应用。S-A模型发展了较多版本,本文通过某弹丸三维外部绕流数值实验发现,标准可压S-A湍流模型[23]对于高雷诺数和低雷诺数的适应能力较强,且在尾部网格不对称情况下能够收敛到较好的结果,因此本文采用该湍流模型计算湍流黏性系数。2 基于CFD的弹丸飞行动力学模型CFD坐标系通常以弹头为原点,Oxc轴与弹轴重合指向弹尾为正,Oyc轴垂直Oxc轴向上为正,Ozc轴满足右手法则。其他坐标系无特殊说明均参考文献[1],弹轴坐标系无特殊说明均指第1弹轴坐标系。旋转稳定弹丸质心运动的动力学方程组通常采用弹道坐标系Otxtytzt下的力系形式[1]。为了使CFD计算的空气动力和空气动力矩能够耦合到弹道模型中,需建立转换关系。首先将CFD坐标系下的力系转换到第1弹轴坐标系下,再将第1弹轴坐标系下力系转到第2弹轴坐标系下,最后将第2弹轴坐标系下的力系转到弹道坐标系下。CFD坐标系到第1弹轴坐标系的转换矩阵为Lac=diag[-1 1 -1].(12)第1弹轴坐标系到弹道坐标系的转换矩阵为Lva=(13)式中:δa和δd分别为高低攻角和方向攻角,总攻角为第2弹轴坐标系到第1弹轴坐标系转过的角度。由(12)式和(13)式将弹道坐标系空气动力转到CFD坐标系下,得到弹道坐标系下弹丸质心运动的动力学方程组为(14)式中:vt为弹丸质心运动速度;θa为速度高低角;ψd为速度方向角;φd为弹轴方位角;m为弹丸质量;g为重力加速度;Fcξ、Fcη、Fcζ分别为空气动力在CFD坐标系Oxc轴、Oyc轴、Ozc轴上的分量;Lc为系数矩阵,Lc=diag[1/m1/mvtcosψd1/mvt].(15)弹丸质心运动的运动学方程组[1]为(16)式中:vpx、vpy、vpz分别为弹丸质心速度矢量vp在地面坐标系Oexeyeze中Oexe轴、Oeye轴、Oeze轴上的分量。由(12)式和基准坐标系到弹轴坐标系的转换矩阵[1],得到控制体表面附加平动速度分量为(17)式中:vcξ、vcη、vcζ分别为控制体表面附加平动速度矢量;φa为弹轴高低角。根据(12)式,得到弹丸总空气动力矩矢量M在弹轴坐标系下的3个分量为(18)式中:Mcξ、Mcη、Mcζ分别为CFD数值计算获得的弹丸总空气动力矩矢量在CFD坐标系3轴上的投影。将(18)式代入弹丸在弹轴坐标系下绕质心转动的动力学方程组[1],可得(19)式中:C为极转动惯量;A为赤道转动惯量;ωξ、ωη和ωζ分别为弹丸角速度矢量在弹轴坐标系上3个轴上的投影分量。由(12)式得到壁面角速度为(ωcwξ,ωcwη,ωcwζ)=(-ωξ,ωη,-ωζ).(20)弹轴坐标系角速度矢量ωa在弹轴坐标系3轴分量为(ωaξ,ωaη,ωaζ)=(ωζtanφd,ωη,ωζ)[1],通过转换矩阵(12)式得到控制体表面附加角速度为(ωcξ,ωcη,ωcζ)=(-ωζtanφd,ωη,-ωζ),(21)对于高速旋转弹丸,由于|ωζtanφd|远小于|ωξ|,由(3)式可知,在弹轴坐标系下控制体表面附加速度远小于弹体坐标系,特别是远场。弹道方程组中绕质心转动的运动学方程组和δa、δd、β3个角关系方程参见文献[1],不再赘述。3 CFD/RBD耦合计算方法本文弹道和流场耦合计算框架如图1所示,本节主要研究图1中的耦合模块。图1 CFD/RBD耦合计算框架Fig.1 Framework of CFD/RBD coupling calculation(14)式、(16)式、(19)式和文献[1]中的绕心转动运动学方程组以及3个角关系方程,共同组成了由CFD数值计算的空气动力和空气动力矩驱动弹道方程。(17)式、(20)式和(21)式将弹道方程解算出的状态量反馈给CFD求解器,从而形成了耦合机制。为叙述方便,(17)式、(20)式和(21)式统称为接口模型。耦合计算采用经典的4阶Runge-Kutta法。已知n时间层面值(tn,y1,n,y2,n,…,ym,n),tn表示数值计算到第n步的时间,ym,n表示第n步下第m个参数,步长为h. 则求解n+1时间层面各函数值的公式为(22)式中:(23)fi表示第i个微分方程的右端表达式。由(23)式可知,每一个Runge-Kutta子步都需要气动数据的支持。有2种较简单的气动力和力矩耦合方法:第1种为定值法,借鉴系数冻结法的思想,将n时间层面CFD计算出的瞬态气动力和力矩冻结,即在Runge-Kutta 4个子步中保持不变;第2种为插值法,利用n-1和n时间层面的气动力和力矩,在4个子步中进行线性外插。本文提出一种紧耦合实现方法,使每一个子步中的气动力和力矩都为CFD计算出的瞬态结果。具体实现步骤如下:步骤1保存n时间层面的流场参数,并积分出CFD坐标系下弹丸受到的空气动力和力矩,通过弹道方程右端子式计算出ki1.步骤2将流场参数设置为步骤1保存的n时间层面;根据ki2右端括号中的参数形式,利用步骤1计算出的ki1和接口模型计算出控制体表面附加运动速度。将时间步长设为0.5h,采用双时间步将流场推进到tn+0.5h时刻;积分出弹丸受到的空气动力和力矩。利用弹道方程右端表达式计算出ki2.步骤3将流动情况返回到n时间层面;根据ki3右端参数形式,利用步骤2计算出的ki2和接口模型计算出控制体表面附加运动速度;双时间步步长设置为0.5h将流场推进到tn+0.5h时刻;积分出作用在弹丸上的空气动力和力矩;利用弹道方程右端表达式计算出ki3.步骤4将流动数设置为tn时刻的值,通过步骤3计算出的ki3来获得ki4右端参数值,由接口模型获取控制体表面附加运动情况;以h为步长,采用双时间步将流场推进到tn+h时刻并获取弹丸上受到的气动力和力矩;由弹道方程右端表达式计算出ki4.步骤5由(22)式计算出n+1时间层面的弹道参数,并通过接口模型计算出控制体表面附加运动;将流场返回至tn时刻,以h为步长,利用双时间步将流场推进至tn+h时刻;积分出当前的空气动力和力矩。图2 流场和弹道紧耦合计算过程Fig.2 Flow field and ballistic tight coupling calculation process通过以上5步可完成流场和弹道方程从n时间层面到n+1时间层面完整的4阶Runge-Kutta紧耦合计算过程,计算流程示意图如图2所示,图2中j=1,2,3,4表示Runge-Kutta4个步骤。需要注意的是,步骤2~步骤4中涉及到的tn+0.5h和tn+h时刻只是中间过程,步骤5中采用双时间步推进到的tn+h时刻为n+1时间层面。4 仿真验证与分析4.1 M549弹丸定态运动和规定运动模拟计算模型采用美国M549旋成体弹丸,具体结构尺寸参见文献[24]。利用Solidworks商业软件进行几何建模,如图3所示。图3 M549弹丸模型Fig.3 M549 projectile model将几何模型导入CFD前处理软件ICEM中进行网格划分,如图4所示。文献[25]对旋成体弹丸的网格划分进行了详细研究,并通过与试验数据对比验证了网格划分的合理性。本文在网格划分过程中部分技术参数参考此文献,如远场位置、近壁第1层网格厚度和黏性层网格延展率等。为使弹丸头部和尾部网格具有良好的正交性,进行O型剖分,网格总数约为186万。图4 M549弹丸计算域网格Fig.4 Computational domain mesh of M549 projectile通过使弹丸做定态飞行和规定运动获得弹丸的气动特性,将结果与Spinner数据[24]进行对比,以检验本文CFD数值方法的可靠性。首先对弹丸定态飞行进行模拟,仿真条件为:马赫数Ma=2.0;高低攻角δa分别为0°和3°;方向攻角δd=0°;弹丸转速为1 360 rad/s. 得到除俯仰阻尼力矩系数导数以外的气动系数,其中零升阻力系数CD0、法向力系数导数Cnα和静力矩系数导数Cmα等静态气动系数误差相对较小,均在6%以内;滚转阻尼力矩系数导数Clp和马格努斯力矩系数导数Cnpα等动导数误差相对较大,分别为9.91%和7.36%,结果如表1所示。为了得到令弹丸做小幅俯仰运动。规定方向攻角δd=0°;一个振荡周期200个计算点;高低攻角δa随时间变化规律为(24)式中:δa0为高低攻角平衡位置,为振幅,为角速度,ω=2vtk/D,D为弹径,k为减缩频率,k=0.03.仿真得到的俯仰力矩系数Cm迟滞回线如图5所示。通过对一个振荡周期进行积分处理[26]得到仿真结果误差为8.04%,如表1所示。图5 俯仰力矩系数迟滞回线Fig.5 Hysteresis loop of pitching moment coefficient表1 马赫数为2时M549弹丸气动系数Tab.1 Aerodynamic coefficients of M549 projectile forMa=2气动系数本文数据Spinner数据[24]误差/%零升阻力系数Cd00...26法向力系数导数Cnα2...65静力矩系数导数Cmα3...70滚转尼力矩系数导数俯仰阻尼力矩系数导数Cmq+Cmα·24...04马格努斯力矩系数导数Cnpα0...36总之,在马赫数为2.0时,静态气动系数和动态气动系数的仿真结果相对Spinner数据的误差分别小于6%和10%,数值仿真的气动系数较可靠 M549弹丸自由飞行仿真验证与分析文献[27]基于传统弹道计算方法,采用旋成体弹丸气动模型和CFD仿真气动系数计算某旋成体弹丸弹道,并与靶道试验结果对比,验证了CFD数值方法结合传统弹道计算方法预测旋成体弹丸实际弹道的可行性。本文借鉴文献[27]的方法计算M549旋成体弹丸弹道,并以该弹道结果为基准,研究第3节中3种力和力矩耦合方法以及不同时间步长对耦合计算方法模拟弹道能力的影响。弹道起始条件为:射角45°;射向0°;高低攻角和飞行马赫数分别沿用4.1节中的3°和2.0;方向攻角0°;弹轴无摆动。流场仿真条件为:双时间步子迭代最大计算步数为200;收敛条件为10-2;流场初始化采用4.1节中对应条件下的计算收敛结果。对于传统6自由度弹道计算,时间步长通常取5 ms左右[1]。为了研究3种力和力矩耦合方法求解弹道的能力,取步长h=5 ms进行仿真,得到攻角关系曲线如图6所示。通过观察图6中的攻角变化情况发现:采用定值法效果最差,仿真攻角结果发散,计算弹道失败;插值法和紧耦合法均模拟出一般旋转稳定弹丸的双圆运动,其中紧耦合法计算结果与基准结果更加接近,即紧耦合精度最好。图6h=5 ms时3种耦合方法对比Fig.6 Comparison of 3 coupling methods forh=5 ms对于非定常流场,为了准确反映流动细节,时间步长要小得多。流场与弹道耦合计算的主要矛盾是非定常气动参数的收敛程度,对某155 mm口径弹丸进行非定常气动力计算时发现,时间步长取0.5 ms和0.05 ms时滚转力矩系数差别仅为4.67%[28],取时间步长为0.5 ms时兼具精度和效率。图7 弹道参数仿真结果Fig.7 Simulated results of trajectory parameters采用紧耦合方法,在时间步长分别为5 ms和0.5 ms条件下进行仿真,仿真结果与基准结果如图7所示。由图7(a)可见:总攻角δt初始为3°,在陀螺稳定力矩作用下逐渐衰减;步长5 ms下总攻角计算结果偏大,衰减明显不足;步长0.5 ms下衰减明显加快,且衰减速度能够跟随基准结果;从二者均值上看,初值均在4°左右,当时间为1.805 s时,前者衰减到1.714 8°,后者衰减到1.170 4°,小步长比大步长衰减效率提高了约24%. 从图7(b)和图7(c)中弹轴摆动曲线以及图7(e)中速度方向角曲线可以看出:当步长为5 ms时仿真结果在第1个快圆运动周期内与基准结果基本一致,但之后不能保持,虽然曲线形状基本一致,但幅值偏大;当步长为0.5 ms时,仿真结果与基准结果始终保持较高的符合度。小步长下速度高低角与基准结果基本重合,大步长下在波谷处略低,如图7(d)所示。弹道侧偏ze反映了地面坐标系Oexeyeze中弹丸偏离射击面Oexeye[1]的程度,如图7(f)所示,仿真结果与基准结果符合较好,在t=1.80 s时,0.5 ms步长和5 ms步长下的弹道侧偏量与基准结果偏差分别为2.18%和4.92%.4.3 修正组件不转时二维弹道修正弹模拟与分析旋转稳定二维弹道修正弹通过后体高速旋转来保持飞行稳定。控制时,修正组件滚转角固定在某个位置,对弹道进行修正。此时亦可采用本文方法进行仿真。计算模型和网格划分如图8所示,网格划分技术参数与4.1节相同,网格数约260万。仿真时,修正组件滚转角固定在图8(a)位置即滚转角为0°;初始高低攻角为0°;后体初始转速1 500 rad/s;其余条件与4.2节相同。图8中:F1和F3表示一对差动舵受到的空气动力;F2和F4表示一对控制舵受到的空气动力。图8 二维弹道修正弹计算模型和网格Fig.8 Calculation model and meshes of two-dimensional trajectory correction projectile图9给出了步长0.5 ms下紧耦合法计算的攻角曲线及特征点P1~P6,攻角变化规律为双圆运动,并与解析方法[29]计算的结果对比,两种方法在频率和幅值方面大致相似,但不完全相同。从攻角收敛趋势来看,控制力向上舵偏角产生的攻角向下,与鸭舵布局的尾翼稳定弹相反,这一结论与文献[29]一致。图9 攻角仿真结果Fig.9 Simulated result of angle of attack图10 弹丸表面压力分布Fig.10 Pressure distribution on the surface of projectile为了分析弹丸飞行过程中表面压力变化情况,观察特征点处的流场,如图10所示,视图方向垂直于弹道坐标系中的Otxtyt平面,特征点参数如表2所示。由起始阶段图10(a)可以看出,与超音速来流直接接触的头部和鸭舵稍部压力较高,最大处超过0.17 MPa;前后体连接处压力分布在0.112~0.156 MPa之间,其中鸭舵正后方处在鸭舵激波内部,压力比周围高;图10(b)中弹轴摆动到了左上方,弹头部压力出现不对称,迎风区压力增大;随着弹轴从左上方摆到左下方,弹头迎风区变成右上方,且由于低头使得水平舵有效迎角减小,压力降低,如图10(c)所示;当弹轴从左下方位置摆到图10(d)中下方略偏右位置时,上方差动舵有效迎角变大同时处于迎风区,下方有效迎角减小同时处于背风区,使得前者压力增大后者压力减小;当弹轴从图10(d)位置略微上摆至图10(e)位置时,水平舵有效迎角增大,压力也增大,直至摆动到图10(f)时,弹丸近似完成一个慢圆周期。紧耦合计算中,随着弹轴摆动,弹丸表面迎风区域和背风区域不断变化,形成时变且复杂的压力分布,压力分布进而又影响弹轴摆动,这与弹丸实际飞行过程是类似的。综上所述可知,与普通弹道计算方法相比,紧耦合法计算不仅能向设计人员提供弹道数据,还能提供流场参数,展现更加实际的弹丸飞行过程。表2 二维弹道修正弹攻角仿真结果Tab.2 Simulated results of angle of attack of two-dimensional trajectory correction projectile特征点t/sδa/(°)δd/(°)××××.×××10-15 结论耦合计算流体力学和刚体动力学方法是外弹道学一个新的发展方向。本文考虑高速旋转弹丸的特点,利用壁面速度表征旋成体表面的高速旋转,结合高速旋转弹丸弹道模型,将弹轴坐标系平动速度和转动角速度附加于控制体表面,建立了基于ALE形式新的流动模型。在此基础上提出一种适合于高速旋转弹丸的基于4阶Runge-Kutta法流场和弹道耦合数值计算方法。最后通过对比传统方法的仿真结果和解析结果验证了紧耦合计算方法的有效性。得到主要结论如下:1)在进行耦合计算时,气动力和力矩耦合模型对计算结果影响较大,采用线性插值和提供实时瞬态值即紧耦合时计算结果收敛性好,精度高。2)时间步长的选取很重要,在合理时间步长下进行紧耦合计算时,可以有效模拟弹丸飞行过程。3)紧耦合方法模拟具有复杂外形的旋转稳定弹飞行时,能够反映出弹丸的运动特性,流场变化规律符合客观实际。本文提出的耦合计算方法可为新型旋转稳定弹,如加装了精确制导组件的二维弹道修正弹设计和计算提供一种新的理论和方法。今后可进一步考虑控制部件的运动,为新型旋转稳定弹控制组件的设计提供基础。参考文献(References)[1] 韩子鹏. 弹箭外弹道学[M].北京: 北京理工大学出版社, 2008: Z P. Exterior ballistics of projectiles[M]. Beijing: Beijing Institute of Technology Press, 2008: 120-143. (in Chinese)[2] SAHU J. Time-accurate numerical prediction of free-flight aerodynamics of a finned projectile[J].Journal of Spacecraft and Rockets, 2008, 45(5): 946-954.[3] SAHU J. Computations of unsteady aerodynamics of a spinning body at transonic speeds[C]∥Proceedings of the 27th AIAA Applied Aerodynamics Conference.San Antonio, TX, US: AIAA, 2009.[4] SAHU J. Virtual fly-out simulations of a spinning projectile from subsonic to supersonic speeds[C]∥Proceedings of the 29th AIAA Applied Aerodynamics Conference. Honolulu, HI, US: AIAA, 2011.[5] SAHU J. Unsteady free-flight aerodynamics of a spinning projectile at a high transonic speed[C]∥Proceedings of the AIAA Atmospheric Flight Mechanics Conference and Exhibit.Honolulu, HI, US: AIAA, 2008.[6] SAHU J, FRESCONI F, HEAVEY K R. Unsteady aerodynamic simulations of a finned projectile at a supersonic speed with jet interaction[C]∥ Proceedings of the 32nd AIAA Applied Aerodynamics Conference.Atlanta, GA, US: AIAA, 2014.[7] SAHU J. Time-accurate numerical prediction of free flight aerodynamics of projectiles[J]. Journal of Spacecraft and Rockets, 2006, 45(5): 946-954.[8] JENNA S, COSTELLO M, SAHU J. Projectile aerodynamic coefficient estimation using integrated CFD/RBD and flight control system modeling[C]∥ Proceedings of the AIAA Atmospheric Flight Mechanics Conference.Chicago, IL, US: AIAA, 2009.[9] STAHL J, COSTELLO M. Estimation of projectile aerodynamic coefficients using coupled CFD/RBD simulation results [C]∥ Proceedings of the AIAA Atmospheric Flight Mechanics Conferenc.Toronto, ON, CAN: AIAA, 2010.[10] SAHU J, GATTO S, COSTELLO M. Using computational fluid dynamics-rigid body dynamic (CFD-RBD) results to generate aerodynamic models for projectile flight simulation: ARL-TR-4270[R]. Harford, MD, US: Army Research Laboratory, Aberdeen Proving Ground, 2007.[11] COSTELLO M, GATTO S, SAHU J. Using CFD/RBD results to generate aerodynamic models for projectile flight simulation[C]∥Proceedings of the AIAA Atmospheric Flight Mechanics Conference and Exhibit.Hilton Head, SC, US:AIAA, 2007.[12] SAHU J. Unsteady aerodynamic simulations of a canard-controlled projectile at low transonic speeds[C]∥ Proceedings of the AIAA Guidance, Navigation, and Control Conference.Portland, OR, US: AIAA, 2011.[13] SAHU J, FRESCONI F. Aeromechanics and control of projectile roll using coupled simulation techniques[J]. Journal of Spacecraft and Rockets, 2015, 52(3): 944-957.[14] SAHU J. Unsteady aerodynamic simulations of a canard-controlled projectile at low transonic speeds[C]∥ Proceedings of the AIAA Guidance, Navigation, and Control Conference.Portland, OR, US: AIAA, 2011.[15] WEINACHT P. Coupled CFD/GN&C modeling for a smart material canard actuator: ARL-TR-4265[R]. Harford, MD, US: Army Research Lab Aberdeen Proving Ground, 2007.[16] SAHU J, FRECONI F, HEAVEY K R. Control performance, aerodynamic modeling, and validation of coupled simulation techniques for guided projectile roll dynamics: ARL-TR-6997[R]. Harford, MD, US: Army Research Laboratory, Aberdeen Proving Ground, 2014.[17] SATHYANARAYANA N D, HOFFMANN K. Trajectory prediction using coupled CFD-RBD with dynamic meshing[C]∥Proceedings of the AIAA Scitech 2019 Forum.San Diego, CA, US: AIAA, 2019.[18] 靳晨晖,王刚,王泽汉.子母弹多体分离过程的非定常CFD/RBD数值仿真[J].气体物理,2018, 3(4): C H, WANG G, WANG Z H. Numerical simulation of unsteady CFD/RBD in multibody separation process of cluster munitions[J]. Physics of Gases, 2018, 3(4): 47-63.(in Chinese)[19] 彭程,郭洋.细长体自旋与锥动耦合运动作用下的流动结构研究[J].兵工学报,2018,39(3): C, GUO Y. Research on flow field structure of slender spinning missile under coupling of conical and spinning motions[J].Acta Armamentarii, 2018, 39(3): 519-527. (in Chinese)[20] 肖中云.旋翼流场数值模拟方法研究[D].绵阳:中国空气动力研究与发展中心, Z Y. Investigation of computational modeling techniques for rotor flowfields[D].Mianyang: China Aerodynamics Research and Development Center, 2007. (in Chinese)[21] KITAMURA K, SHIMA E. A new pressure flux for AUSM-family schemes for hypersonic heating computations[C]∥Proceedings of the 20th AIAA Computational Fluid Dynamics Conference.Honolulu, HI, US: AIAA, 2011.[22] KITAMURA K, SHIMA E. Towards shock-stable and accurate hypersonic heating computations: A new pressure flux for AUSM-family schemes[J]. Journal of Computational Physics, 2013, 245: 62-83.[23] RAHMAN M M, AGARWA R K, LAMPINEN M J, et al. Wall-distance-free version of Spalart-Allmaras turbulence model[J]. AIAA Journal, 2015, 53(10): 3016-3028.[24] STUREK W B. Application of CFD to the aerodynamics of spinning shell[C]∥Proceedings of the 22nd Aerospace Sciences Meeting. Reno, NV, US: AIAA, 1984.[25] DESPIRITO J, HEAVEY K. CFD computation of Magnus moment and roll damping moment of a spinning projectile[C]∥Proceedings of the AIAA Atmospheric Flight Mechanics Conference and Exhibit.Providence, RI, US: AIAA, 2004.[26] VISHAL A, BHAGWANDIN J S. Numerical prediction of pitch damping stability derivatives for finned projectiles[J]. Journal of Spacecraft and Rockets,2014,51(5): 1603-1618.[27] WEINACHT P. Projectile performance, stability, and free-flight motion prediction using computational fluid dynamics[J].Journal of Spacecraft and Rockets,2004, 45(2): 257-263.[28] 徐辉雯,陈少松.二维修正弹修正组件反旋与不旋气动特性对比[J].弹道学报,2018,30(4):32-37.XU H W, CHEN S of aerodynamic characteristics of two-dimension trajectory correction projectile with reverse-rotation and non-rotation PGK[J].Journal of Ballistics, 2018, 30(4):32-37.(in Chinese)[29] 钟杨威.旋转稳定二维弹道修正弹弹道特性分析与制导方法研究[D].南京:南京理工大学, Y W. Study on trajectory characteristics and guidance methods for spin stabilized two-dimensional trajectory correction projectiles[D]. Nanjing: Nanjing University of Science and Technology, 2015. (in Chinese)
文章来源:《弹道学报》 网址: http://www.tdxbzz.cn/qikandaodu/2021/0205/335.html