基于CFD数值仿真技术的飞行器动导数计算

米百刚，詹浩，朱军

（西北工业大学 翼型叶栅空气动力国防科技重点实验室，陕西 西安 710072）

基于CFD数值仿真技术的飞行器动导数计算

米百刚，詹浩，朱军

（西北工业大学 翼型叶栅空气动力国防科技重点实验室，陕西 西安 710072）

为了研究表征飞行器动态稳定特性的动导数，结合计算流体力学方法，基于滑移网格技术，采用小幅度强迫振动法和差分法两种非定常方法，推导了俯仰组合动导数以及滚转阻尼导数的计算表达式，并结合多参考系模型采用准定常方法计算滚转阻尼导数。建立了基于计算流体力学技术的飞行器动导数计算方法。对国际动导数标模Einner导弹进行了计算验证。两种非定常方法对俯仰组合动导数的计算误差分别为0．65%、6．13%，对滚转阻尼导数的计算误差分别为2．5%、0．06%。求解滚转阻尼的准定常方法对于滚转阻尼导数的计算误差为2．67%。三种方法的计算结果与风洞试验结果吻合的很好。非定常方法计算精度高，准定常方法可以快捷迅速得到结果。本文的数值方法不仅适用于超声速，同样可以用于亚、跨声速时的动导数计算。

滑移网格；多参考系模型；动导数；计算流体力学技术

0 引 言

本文结合CED计算技术发展了飞行器超声速动导数计算方法。以国外动导数标模［17］Einner导弹为例，基于滑移网格技术，对动态气动设计特性中的滚转阻尼导数，以及俯仰组合动导数进行了模拟，并采用多参考系方法，对滚转阻尼导数进行了数值计算。计算结果表明，无论是基于滑移网格的非定常方法，还是结合多参考系模型的准定常方法，都具有比较高的计算精度。

1 计算方法

1.1 控制方程

文中采用的控制方程为三维积分形式的雷诺平均N-S方程。

其中V为任意控制体，W是守恒变量，F为无粘（对流）通矢量项，Fv为粘性通量，∂V为控制体的边界，n为控制体边界单位外法向矢量，Re为计算的雷诺数。

空间离散采用二阶迎风格式——通量差分分裂（Roe-EDS）格式，采用双时间步推进并用隐式LUSGS格式迭代求解，应用当地时间步长、残值光顺、预处理和多重网格加速收敛。湍流模型采用对逆压梯度流动模拟精度较高的k-ωSST湍流模型。

1.2 滑移网格

滑移网格是一种特殊的动网格，它将问题的求解域划分成多个计算域，即动区域和静止区域。各计算域之间采用定义好的分界面进行关联，动区域沿着滑移面进行运动。在滑移网格问题中，动区域运动是相对于静止参考系进行跟踪的，因此，没有运动参考系附加在计算域上，简化了穿过分界面的通量传递。与一般动网格相比，网格的运动仅沿着滑移面进行整体滑移，没有网格变形重构的过程，相对节省了产生新网格所需的计算量，而且运动时各区域网格质量不发生变化，可以提高计算精度。

1.3 多参考系模型

多参考系模型（MRE）是不同旋转或移动速度的每个区域的稳态近似。同样它将计算域分成多个小的子域，每个子域可以有自己的运动方式，或静止或旋转或平移。此方法适用于边界上流动区域几乎均匀混合的情况。以同时存在静止和旋转域情况为例，在两个子域上分别建立静止坐标系和旋转坐标系，控制方程在每个子域分别求解，在交接面上通过将速度换算成绝对速度形式进行流场信息交换。由于进行的是稳态近似，多参考系模型是一种很高效的方法。

1.4 动导数计算方法

1.4.1 强迫振动法

强迫飞行器做小幅度俯仰振荡时，依据傅立叶公式展开，刚体飞行器所受非定常气动力矩可以写为表达式：

其中Mz是瞬态非定常气动力矩值；Mz0是平衡位置处的气动力矩是基频谐波分量的气动力矩幅值；ω是强迫飞行器振动频率；λ是刚体飞行器强迫振动时位移和气动力矩之间的相位差；u（t）为高次谐波分量。

依据气动力导数的概念，上述小幅度强迫俯仰振动的刚体飞行器所受非定常气动力矩还可以表示为：

当刚体飞行器以低频ω做小幅度振荡时，其模型运动方程可以简化为：

展开式（2）并略去高次谐波分量，同时将式（4）代入略去高阶分量，整理同类项可得：

当非定常问题计算足够长时，令ωt＝2nπ，就可以抹去初始效应的影响，气动力矩达到一个周期性稳态值，于是式（5）可以写为：

式（6）即为基于傅立叶及泰勒展开方法推导出的俯仰组合动导数的表达式。

1.4.2 差分法

假定飞行器以两种不同角速度ωz1，ωz2等速上仰运动至相同攻角，依据飞行力学小扰动理论对其力矩进行展开并略去高阶量可得：

式（9）即为基于差分法得到的俯仰组合动导数公式。

1.4.3 准定常方法

前面提出的强迫振动法以及差分法都是非定常方法，可以用来求解纵向，横航向动导数。而一般的飞行器，具有顺流方向的面对称性或者线对称性，横向气动力在0°迎角，定速转动状态下一般保持不变，且横向更多的是考虑稳定滚转特性，因而通常使用准定常方法进行滚转动导数的计算。

飞行器进行小角速度等速滚转时，令驱动力矩为Mx，阻尼力矩为Mkx，根据小扰动理论：

文中考虑零滚转角下的动稳定性导数，θ＝0，因此绕质心转动动力学方程简化为：

注意到方程左端向为0，即可推导出滚转动导数表达式：

其中Ix为飞行器对ox轴的转动惯量。

1.5 边界条件

本文使用的基于滑移网格以及多参考系模型的求解方法要求将流场划分为动域以及静域，两个域之间通过交接面进行数据传递。动域中的网格刚化后沿着交接面进行滑移运动，以此来实现物面的运动。本文涉及到的边界条件有远场边界，物面边界以及交接面边界。其中远场边界使用基于Riemann不变量的无反射边界条件，物面边界无滑移，绝热，交接面使用基于通量守恒的交接面边界。

2 动导数计算公式无量纲化

一般来讲，计算动导数的公式都是无量纲的，因此，需要对本文的几种方法进行无量纲化处理。

2.1 强迫振动法

根据国家标准，减缩频率取为k＝ωl／2V*，其中l为参考长度，V*为远场自由来流速度，代入可得：

2.2 差分法

引入无量纲角速度，

带入原式，

可得无量纲差分法公式：

2.3 准定常方法

3 动导数算例及分析

3.1 计算模型及网格

计算模型为国际标准动导数模型Einner导弹，图1为计算用的Einner导弹模型及网格，其中d＝1m，采用结构网格，根据本文的方法将流场分为动域和静域，两部分通过interface交接面连接，这个交接面的网格不要求完全一致，给复杂模型网格划分带来方便。最终划分的网格总量为390万，其中，动域网格量260万，静域130万。

图1 计算模型及网格Fig.1 Computational model and grid

3.2 算例验证

3.2.1 强迫振动法求解俯仰组合动导数

使用非定常方法求解动导数的关键在于准确地预测瞬时力矩。本文选取M∞＝1．58，α0＝0°为计算状态，减缩频率为k＝0．0158226，振幅为αm＝1°，也即振荡规律可写为

图2为计算得到的瞬时力矩随迎角的变化曲线。可以明显地看出，本文的非定常强迫振动方法能够很好地预测迎角与力矩系数之间的迟滞效应，计算的俯仰力矩系数迟滞环与文献［1］中的非常一致。

图2 力矩系数随迎角变化规律图Fig.2 Time history of pitching moment coefficient with angle of attack in pitching cycle

根据前面得到的动导数无量纲计算公式，可求解俯仰组合动导数，计算结果见表1。

表1 俯仰组合动导数Table 1 Pitchingdynamic stabilityderivative by using small oscillation method

文献［1］中的俯仰组合动导数计算值为-506，误差为3．62%。由此可见本文的计算方法对俯仰组合动导数的计算是比较准确的。

一般使用强迫振动方法计算动导数多使用积分法进行结果辨识，本文则使用稳定后的平衡位置处单点气动数据计算动导数。表2是两种方法的对比，可以看出，两者均能比较准确的得到动导数结果，相比下，本文的单点数据计算量小，能够方便快速得到结果。

表2 辨识方法对比Table 2 Identification methods of dynamic derivatives

为了进一步就本文的小幅度强迫振动方法进行验证，分别计算在马赫数为1．58，1．74，1．88，2．1，2．55五个状态下的俯仰组合动导数。

图3的结果表明，重心位置在Xg＝5．0d时，本文所计算的动导数与风洞试验数据较为一致，反映出在超音速范围内，随着马赫数的增加，俯仰组合动导数的逐渐减小，且趋势变缓和。由于目前的技术限制，动导数试验结果的误差带相当大，与试验结果相比，本文采用数值模拟技术较为准确地计算了Einner导弹的俯仰组合动导数，在文中给定的马赫数范围内很好地模拟出动导数的变化趋势。

图3 不同马赫数下的俯仰组合动导数Fig.3 Pitching dynamic stability derivative versus Mach number

3.2.2 差分法求解俯仰组合动导数

本文的计算初始迎角均为0°，导弹以5°／s，10°／s的角速度分别等速上仰至5°。然后根据得到的瞬时力矩即可求出俯仰组合动导数。计算结果见表3。

表3 差分法求解俯仰组合动导数Table 3 Pitching dynamic stability derivative by usingdifferential method

表3的结果表明，使用差分法也能比较准确地计算俯仰组合动导数。

3.2.3 准定常方法求解滚转动导数

给定导弹以5°／s的速度绕轴滚转，使用前面所述的准定常方法，计算得到的滚转阻尼导数为：

同样使用强迫振动法以及差分法求解滚转动导数，并与准定常方法的结果进行比较，可得表4。

由表4可以看出，准定常虽然计算量小，耗时短，但是精度较低。而强迫振动以及差分方法精度比较高。对于动导数的数值计算，关键在于计算非定常气动力，非定常方法基于这一点提出，因而精度也稍高，相对而言，准定常方法是一种近似的方法，对流场特性细节的捕捉能力稍弱，所以精度较低。

表4 滚转动导数求解结果Table 4 Rolling damp dynamic derivative

4 结 论

采用基于滑移网格的非定常方法以及基于多参考系模型的准定常方法对Einner导弹标准模型动导数进行了计算研究和分析。相比以往的使用网格重构进行计算的方法，滑移网格的引入使得计算量减小，并且由于不存在网格变形带来的影响，因此可以提高求解的精度。由本文的算例可以看出所发展的非定常方法以及准定常方法能够较为准确地计算俯仰及滚转动导数。其中小幅度强迫振动法以及差分法两种非定常方法的关键在于准确预测瞬时力矩；准定常方法虽然精度稍差，但是计算量小，快捷。同时，本文的三种方法不仅可以用于超音速动导数辨识，也同样适用于亚音速以及跨音速的动导数计算。

本文采用的动导数辨识方法给出了与风洞试验较为吻合的数值计算结果，以及一致的动导数的变化趋势，表明所采用的动导数计算方法的正确性，具备较好的工程应用前景。

［1］GREEN L L，SPENCE A M，MURPHY P C．Computation methods fordynamic stability and controlderivatives［R］．AIAA 2004-0015，2004．

［2］SHI A M．Numerical analysis of aircraft unsteady flow and aeroelastic problems［D］．Xi’an：Northwestern Polytechnical University，2005．（in Chinese）

史爱明．飞行器非定常流场和气动弹性问题的数值分析研究［D］．西安：西北工业大学，2005．

［3］RONCH A D，VALLESPIN D．Computation ofdynamicderivatives using CED［R］．AIAA 2010-4817，2010．

［4］GREEN L L，SPENCE A M，MURPHY P C．Computation methods fordynamic stability and controlderivatives［R］．AIAA 2004-0015，2004．

［5］RONCH A D，GHOREYSHI M，BADCOCK K J，et al．Linear frequencydomain and harmonic balance predicitions ofdynamicderivatives［R］．AIAA 2010-4699，2010．

［6］PARK M A，GREEN L L．Steady-state computation of constant rotational ratedynamic stabilityderivatives［R］．AIAA 2000-4321，2000．

［7］KRAMER，BRIAN R．Experimental evaluation of superposition techniques applied todynamic aerodynamics［R］．AIAA 2002-0700，2002．

［8］ERDAL O，HSSAN U A．CED predicitions ofdynamicderivatives for missile［R］．AIAA 2002-0276，2002．

［9］SOO H P，YOONSIK K，JANG H K．Prediction ofdynamicdamping coefficients using unsteadydual-time stepping method［R］．AIAA 2002-0715，2002

［10］SCOTT M M，ELORET C．A reduced-frequency approach for calculatingdynamicderivatives［R］．AIAA 2005-840，2005．

［11］YUAN X X，ZHANG H X，XIE Y E．The pitching static／dynamicderivatives computation based on CED methods［J］．ACTA Aerod ynamica Sinica，2005，23（4）：458-463．（in Chinese）

袁先旭，张涵信，谢昱飞．基于CED方法的俯仰静，动导数数值计算［J］．空气动力学学报，2005，23（4）：458-463．

［12］SUN T，GAO Z H，HUANG J T．Identify of aircraftdynamicderivatives based on CED technology and analysis of reduce frequency［J］．Elight Dynamics，2011，29（4）：15-18．（in Chinese）

孙涛，高正红，黄江涛．基于CED的动导数计算与减缩频率影响分析［J］．飞行力学，2011，29（4）：15-18．

［13］SHI A M，YANG Y N，YE Z Y．A more accurate method for calculating transonicdynamicderivatives（TDDs）using present state-of-the-art CED［J］．Journal of Northwestern Polytechnical University，2008，26（1）：11-14．（in Chinese）

史爱明，杨永年，叶正寅．结合CED技术的跨音速动导数计算方法研究［J］．西北工业大学学报，2008，26（1）：11-14．

［14］LU X C，YE Z Y，ZHANG W W．A high efficient method for computingdynamicderivatives of supersonic／hypersonic aircraft［J］．Aeronautical Computing Technique，2008，38（3）：28-31．（in Chinese）

卢学成，叶正寅，张伟伟．超音速、高超声速飞行器动导数的高效计算方法［J］．航空计算技术，2008，38（3）：28-31．

［15］SUN Z W，CHENG Z Y，BAI J Q，et al．A high efficient method for computingdynamicderivatives of aircraft based on quasisteady CED method［J］．Elight Dynamics，2010，28（2）：28-30．（in Chinese）

孙智伟，程泽荫，白俊强，等．基于准定常的飞行器动导数的高效计算方法［J］．飞行力学，2010，28（2）：28-30．

［16］JIANG S J，LIU Y Q，DANG M L．A calculation method of aircraft roll-damping moment coefficientderivative based on steady NS equation［J］．Journal of Projectiles，Rockets，Missiles and Guidance，2008，28（1）：180-182．

将胜炬，刘玉琴，党明利．基于定常NS方程的飞行器滚转阻尼力矩系数导数计算方法［J］．弹箭与制导学报，2008，28（1）：180-182．

［17］GREENWELL D L．Erequency effects ondynamic stabilityderivatives obtained from small-amplitude oscillatory testing［J］．Journal of Aircraft，1998，35（5）：776-783．

Calculation ofdynamicderivatives for aircraft based on CFD technique

MI Baigang，ZHAN Hao，ZHU Jun
（National Key Laboratory of Aerodynamic Design and Research，Northwestern Polytechnical University，Xi’an 710072，China）

To study thedynamic stabilityderivatives of aircraft with CED technology，an unsteady method for calculating pitchingdynamic stabilityderivative and rollingdampdynamicderivative isdeveloped by using small amplitude oscillation model anddifferential model，and a quasi-steady method is also built to calculate rollingdampdynamicderivative；the unsteady method is based on sliding mesh technique while the quasi-steady method on both sliding mesh technique and moving reference frame model；compared to thedynamic mesh used before，the sliding mesh technique can improve the computational precision as it avoids the effect of meshdeformation．The finner missile is taken as example to calculate thedynamicderivatives．The errors of two unsteady methods to calculate longitudinal combinedderivatives are 0．65%and 6．13%respectively，while errors to lateral combinedderivatives are 2．5%and 0．06%，and the error of quasi-steady method to obtain lateral combinedderivative is 2．67%．The research shows the calculation result is consistent with the experimentaldata，although the quasi-steady method is less accurate than the unsteady one，nevertheless it can get result more quickly．Moreover，both the methods can be used for supersonic，subsonic and transonicdynamicderivative calculation．

sliding mesh；moving reference frame model；dynamicderivative；CED technology

V211.3

Adoi：10.7638／kqdlxxb-2012.0206

0258-1825（2014）06-0834-06

2012-12-13；

2013-02-25

米百刚（1989-），男，陕西富平人，博士研究生，研究方向：计算流体力学．E-mail:mibaigang＠mail.nwpu.edu.cn

米百刚，詹浩，朱军．基于CED数值仿真技术的飞行器动导数计算［J］．空气动力学学报，2014，32（6）：834-839．

10．7638／kqdlxxb-2012．0206 MI B G，ZH AN H，ZHU J．Calculation ofdynamicderivatives for aircraft based on CED technique［J］．ACTA Aerodynamica Sinica，2014，32（6）：834-839．