CN103400035A - 一种高可信度快速预测飞行器滚转动导数的方法 - Google Patents

Abstract

本发明公开了一种高可信度快速预测飞行器滚转动导数的方法，包括如下步骤：模型表面网格的生成和空间网格划分；旋转坐标系下气动参数的计算方法：（1）将惯性系下的Navier-Stokes方程变换到旋转坐标系下；（2）对变换后的方程进行数值求解，获得每个状态的流场；（3）对物面压力和粘性应力在表面上积分得到气动力，通过表面力向质心取矩并积分得到作用在质心的气动力矩；滚转动导数的差分计算：分别计算两次不同旋转速度下模型的气动力和气动力矩，然后通过差分方法计算得到模型的滚转动导数。本发明方法计算结果可信度高，同时计算量较需要时间精确求解的非定常强迫振动方法小的多，可以高可信度快速预测飞行器滚转动导数。

Description

一种高可信度快速预测飞行器滚转动导数的方法

技术领域

本发明涉及流体力学领域，具体涉及一种高可信度快速预测飞行器滚转动导数的方法。

背景技术

动导数为飞行器运动中单位速度或角速度变化引起的气动力或力矩的变化，以导数形式表达。飞行器的滚转动导数是设计飞行器导航系统和控制系统以及动态品质分析所不可缺少的原始气动参数，滚转动导数的准确评估对飞机设计和飞行都有重要的意义。随着现在先进战斗机和导弹飞行包线的扩展，以及大迎角和带侧滑情况下的机动飞行等，人们对复杂流场中的飞行器动稳定性研究更加重视。近年来随着国际形势紧张及各国之间高科技竞争的日趋激烈，国家对我国航空航天事业持续加大了关注和投入力度，各类新型飞行器正在加速研发，这些飞行器的飞行品质和操稳特性的分析以及控制律的设计是很重要的研究内容。针对复杂外形以及一些非传统外形飞行器，在气动外形设计的初期发展滚转动导数的预测方法以减少风洞试验的数量，缩短初期设计周期是十分必要的。因此，如何获得滚转动导数是设计者面临的一项重要工作。

飞行试验，风洞实验以及理论计算是目前获取滚转动导数的主要方法。飞行试验的难度大、周期长、费用高。国家发明专利200810137422（旋转流场下的滚转动导数测量系统）公开了一种用于风洞动态综合试验中旋转流场下的滚转动导数测量系统，针对在非定常流场下的气动数据测量，需要研制飞机非定常动态综合实验系统，特别是需要研制旋转流场下的滚转动导数测量系统的问题，提供旋转流场下的滚转动导数测量系统，该发明技术方案由振荡轴控制机构和旋转轴控制机构构成，对于旋转轴控制采用速度控制，对振荡轴控制采用位置随动控制。采用控制工控机、ＣＡＮ总线控制器、旋转电机驱动器、振荡同步信号发生器、驱动旋转电机编码和减速器等装置，它与尾旋测试系统、强迫振荡测试系统等各动态实验系统一起，可以提供飞机大攻角机动性、操纵性、稳定性、失速尾旋预测和分析等全套数据，特别对研究振荡尾旋有重要意义。其主要用于旋转流场下的滚转动导数测量，但是同样难度较大、成本高、同时洞壁干扰，支撑干扰，模型震动等因素也影响测量精度。

随着计算机技术及数值计算方法的发展，通过数值计算获得滚转动导数是一种最理想的方式。经典的滚转动导数计算方法主要采用基于线化理论的气动力模型，如线化升力面方法和细长体理论，这些方法对相对简单的外形如旋成体等构型的稳定性导数预测精度较高、速度也较快，但是其不能考虑到各种非线性效应。当飞行器进入大迎角状态或姿态变化较大时，流场中非定常分离和旋涡运动、非定常运动激波及其与附面层的干扰占主导地位，其气动力系数变化呈现强烈的非线性效应。此时，传统上基于线性小扰动理论的研究方法如摄动法等都不再适于表征飞行器的动稳定性，即在大迎角情况下，线性模型不再适用。

20世纪90年代以来，随着计算流体力学方法的发展，多重网格法、局部步长法、隐式残差光顺法、隐式迭代等许多加速流场收敛的技术得以成功的采用同时随着计算机技术的发展以及并行计算等技术的应用使采用求解Navier-Stokes方程加湍流模型的方法来预测滚转动导数的方法成为可能。目前CFD(Computational Fluid Dynamics)技术已经逐步用于滚转动导数的计算，主要在惯性系下采用数值模拟强迫振动或者自由衰减运动两种方法求解流体方程，再通过数值方法辨识出各类导数。这种基于CFD技术的数值方法可以考虑到流场的非线性特性，适用更大的迎角范围等条件，可用于开展复杂外形的气动力计算，但是这类方法需要开展时间精确求解，计算量大，当模拟的外形较复杂时对计算资源的需求更加庞大，不能满足型号设计的要求。

寻求一种高可信度快速预测飞行器滚转动导数方法对加快飞行器设计是十分必要的。

发明内容

为了克服现有技术的上述缺点，本发明提供了一种通过在旋转坐标系下求解定常流体力学控制方程来获取高可信度飞行器滚转动导数的方法。

本发明解决其技术问题所采用的技术方案是：一种高可信度快速预测飞行器滚转动导数的方法，包括如下步骤：

第一步、模型表面网格的生成和空间网格划分：

第二步、旋转坐标系下气动参数的计算方法：

（1）将惯性系下的Navier-Stokes方程变换到旋转坐标系下；

（2）对变换后的方程进行数值求解，获得每个状态的流场；

（3）对物面压力和粘性应力在表面上积分得到气动力，通过表面力向质心取矩并积分得到作用在质心的气动力矩；

第三步、滚转动导数的差分计算：分别计算两次不同旋转速度下模型的气动力和气动力矩，然后通过差分方法计算得到模型的滚转动导数。

本发明主要针对滚转动导数的计算。传统的计算方法中由于旋转速度的存在，需要把模型考虑成运动问题，通过求解惯性系下的Navier-Stokes方程来获得，这样涉及到动网格生成、非定常的流动计算等复杂的数值模拟过程。为此，本发明将着眼点放在模型坐标系上，将原本在惯性系中旋转的模型看做是在模型坐标系下静止不动的，而是远场来流在旋转，如果模型的旋转速度保持不变的话，那么远场的旋转来流是一种不随时间变化的稳态旋流。基于这种原理，本发明将原先在惯性坐标系下的Navier-Stokes方程投影到模型坐标系下，得到了旋转坐标系下的流动控制方程，这样将原本惯性系下的非定常流动动网格模拟问题转化为非惯性系下定常的流动问题了。所以要计算滚转动导数，就是要计算各个运动状态下的气动参数，例如在旋转角速度为ω1时的气动力F1和气动力矩M1。

通过计算不同运动状态下的气动参数，即计算不同旋转速度（例如：旋转角速度分别为ω1、ω2）下的气动参数（例如：气动力F、气动力矩M等）后，通过简单差分方法就可以直接计算得到相应的滚转动导数。

与现有技术相比，本发明的积极效果是：本发明方法计算结果可信度高，同时计算量较需要时间精确求解的非定常强迫振动方法小的多，可以高可信度快速预测飞行器滚转动导数。

附图说明

本发明将通过例子并参照附图的方式说明，其中：

图1为惯性坐标系与非惯性坐标系示意图；

图2为流体粒子在惯性系及非惯性系中运动关系示意图；

图3为旋转坐标系示意图；

图4为典型外形物面及空间网格示意图；

图5为气动力收敛过程示意图；

图6为不同旋转速度下的流场结构；

图7本方法计算结果与强迫振动方法计算结果对比。

具体实施方式

一种高可信度快速预测飞行器滚转动导数的方法，包括以下步骤：

第一步、模型表面网格的生成和空间网格划分：

首先，对需要获得滚转动导数的模型进行表面网格和空间网格划分。由于目前所需求解滚转动导数的模型日益复杂，本发明采用多块对接网格技术保证网格质量（指正交性、Y⁺、增长率等），物面附近划分边界层网格，以精确捕捉边界层内的信息。带控制舵导弹标模的物面及空间网格如图4所示。

第二步、旋转坐标系下气动参数的计算方法：

根据本发明的思想，将着眼点放在和模型固结在一起的旋转坐标系下，将原先在惯性坐标系下的Navier-Stokes方程投影到模型旋转坐标系，这样将原本惯性系下的非定常流动动网格模拟问题转化为非惯性系下定常的流动问题。通过定常计算的方法获取各个旋转速度下的气动力和气动力矩。

在所需要获取滚转动导数的飞行状态（包括马赫数和姿态角）附近针对需要求的导数变量，采用以下数值方法进行不同导数变量的展开求解，以获取不同旋转速度下的流场结构（如图6所示）。按照本发明的思想，首先就是需要将惯性系下的Navier-Stokes方程变换到旋转坐标系（和模型固结在一起的非惯性系）下，其次对变换后的方程建立相应的数值求解方法，最后对数值结果进行积分运算得到模型的气动参数。具体的数值求解方法如下：

（1）Navier-Stokes方程从惯性系到旋转坐标系的变换关系推导

涉及到流体力学控制方程在惯性系及非惯性系之间的转换关系，在本方法中，计算网格（图1中用正方形小格示意）在非惯性坐标系中进行处理。惯性坐标系在空间中位置固定，其三维坐标用B（X，Y，Z）表示，三维方向矢量用（I，J，K）表示。非惯性坐标系在空间中可以旋转，其三维坐标用b（x，y，z）表示，三维方向矢量用（i，j，k）表示，绕三个轴旋转速度分别用（ω_x，ω_y，ω_z）表示，其组成的旋转矢量用ω来表示。惯性坐标系（固定系）与非惯性坐标系（运动系）的示意图如图1所示。

流体粒子在惯性系及非惯性系中运动关系示意图如图2所示，其中B表示流体粒子到惯性坐标系原点的矢径，b表示流体粒子到非惯性坐标系原点的矢径，C表示非惯性坐标系原点到惯性坐标系原点的矢径。

可以得出如下矢量表达式：B＝C+b    (1）

对应的速度表达式为 $\frac{\mathrm{dB}}{\mathrm{dt}}=\frac{\mathrm{dC}}{\mathrm{dt}}+\frac{\mathrm{db}}{\mathrm{dt}}---\left(2\right)$

通过矢量算法当dt趋近于0时（Meirovitch，L.，“Methods ofAnalytical Dynamics，”McGraw-Hill，New York，NY，1970.）得到：

$\stackrel{\·}{B}=\stackrel{\·}{C}+\stackrel{\·}{b}+\mathrm{\ω}\×b---\left(3\right)$

远场边界上由于

将以上带入方程（3）中可以得到非惯性系下的远场边界条件为，如图3所示：

$\stackrel{\·}{b_{\∞}}=u_{\∞}-\mathrm{\ω}\×b_{\∞}---\left(4\right)$

利用类似的运算方式可以得到，加速度关系式为：

$\stackrel{\·\·}{B}=\stackrel{\·\·}{C}+\stackrel{\·\·}{b}+2\mathrm{\ω}\×\stackrel{\·}{b}+\mathrm{\ω}\×(\mathrm{\ω}\×b)---\left(5\right)$

为旋转坐标系相对惯性系的加速度，其表达式为

$\stackrel{\·\·}{C}=\mathrm{\ω}\×\stackrel{\·}{C},\stackrel{\·\·}{C}=\mathrm{\ω}\×(-u_{\∞})---\left(6\right)$

所以可以得到伪加速度表达式为：

$\stackrel{\·\·}{B}-\stackrel{\·\·}{b}=\mathrm{\ω}\×(-u_{\∞})+2\mathrm{\ω}\×\stackrel{\·}{b}+\mathrm{\ω}\×(\mathrm{\ω}\×b)---\left(7\right)$

（2）旋转坐标系下流动方程的数值求解：

惯性坐标系下的流体力学控制方程

$\frac{\∂Q}{\∂t}+(\frac{\∂E}{\∂x}+\frac{\∂F}{\∂y}+\frac{\∂G}{\∂z})-(\frac{{\∂E}_v}{\∂x}+\frac{{\∂G}_v}{\∂y}+\frac{{\∂G}_v}{\∂z})=0---\left(8\right)$

$Q=\left[\begin{array}{c}\mathrm{\ρ}\\ \mathrm{\ρu}\\ \mathrm{\ρv}\\ \mathrm{\ρw}\\ e\end{array}\right]E=\left[\begin{array}{c}\mathrm{\ρu}\\ {\mathrm{\ρu}}^2+p\\ \mathrm{\ρuv}\\ \mathrm{\ρuw}\\ (e+p)u\end{array}\right]F=\left[\begin{array}{c}\mathrm{\ρv}\\ \mathrm{\ρvu}\\ {\mathrm{\ρv}}^2+p\\ \mathrm{\ρvw}\\ (e+p)v\end{array}\right]G=\left[\begin{array}{c}\mathrm{\ρw}\\ \mathrm{\ρwu}\\ \mathrm{\ρwv}\\ {\mathrm{\ρw}}^2+p\\ (e+p)w\end{array}\right]$

其中： $e=\frac{p}{\mathrm{\γ}-1}+\frac{1}{2}\mathrm{\ρ}(u^2+v^2+w^2)---\left(9\right)$

$E_v=\left[\begin{array}{c}0\\ {\mathrm{\τ}}_{\mathrm{xx}}\\ {\mathrm{\τ}}_{\mathrm{xy}}\\ {\mathrm{\τ}}_{\mathrm{xz}}\\ {\mathrm{u\τ}}_{\mathrm{xx}}+{\mathrm{v\τ}}_{\mathrm{xy}}+{\mathrm{w\τ}}_{\mathrm{xz}}-{\stackrel{\·}{q}}_x\end{array}\right]$ $\begin{array}{c}{\mathrm{\τ}}_{\mathrm{xx}}=\frac{2}{3}\mathrm{\μ}(2\frac{\∂u}{\∂x}-\frac{\∂v}{\∂y}-\frac{\∂w}{\∂z})\\ {\mathrm{\τ}}_{\mathrm{yy}}=\frac{2}{3}\mathrm{\μ}(2\frac{\∂v}{\∂y}-\frac{\∂u}{\∂x}-\frac{\∂w}{\∂z})\\ {\mathrm{\τ}}_{\mathrm{zz}}=\frac{2}{3}\mathrm{\μ}(2\frac{\∂v}{\∂y}-\frac{\∂u}{\∂x}-\frac{\∂w}{\∂z})\end{array}$

$F_v=\left[\begin{array}{c}0\\ {\mathrm{\τ}}_{\mathrm{yx}}\\ {\mathrm{\τ}}_{\mathrm{yy}}\\ {\mathrm{\τ}}_{\mathrm{yz}}\\ {\mathrm{u\τ}}_{\mathrm{yx}}+{\mathrm{v\τ}}_{\mathrm{yy}}+{\mathrm{w\τ}}_{\mathrm{yz}}-{\stackrel{\·}{q}}_y\end{array}\right]$ $\begin{array}{c}{\mathrm{\τ}}_{\mathrm{xy}}=\mathrm{\μ}(\frac{\∂u}{\∂y}+\frac{\∂v}{\∂x})={\mathrm{\τ}}_{\mathrm{yx}}\\ {\mathrm{\τ}}_{\mathrm{yz}}=\mathrm{\μ}(\frac{\∂w}{\∂y}+\frac{\∂v}{\∂z})={\mathrm{\τ}}_{\mathrm{zy}}\\ {\mathrm{\τ}}_{\mathrm{xz}}=\mathrm{\μ}(\frac{\∂u}{\∂z}+\frac{\∂w}{\∂x})={\mathrm{\τ}}_{\mathrm{zx}}\end{array}$

$G_v=\left[\begin{array}{c}0\\ {\mathrm{\τ}}_{\mathrm{zx}}\\ {\mathrm{\τ}}_{\mathrm{zy}}\\ {\mathrm{\τ}}_{\mathrm{zz}}\\ {\mathrm{u\τ}}_{\mathrm{zx}}+{\mathrm{v\τ}}_{\mathrm{zy}}+{\mathrm{w\τ}}_{\mathrm{zz}}-{\stackrel{.}{q}}_z\end{array}\right]$ $\begin{array}{c}{\stackrel{\·}{q}}_x=-\mathrm{\κ}\frac{\∂T}{\∂x}\\ {\stackrel{\·}{q}}_y=-\mathrm{\κ}\frac{\∂T}{\∂y}\\ {\stackrel{\·}{q}}_z=-\mathrm{\κ}\frac{\∂T}{\∂z}\end{array}$     (10）

式中ρ,p,T,u,v,w分别为流场的密度、压力、温度、以及x，y，z方向上的速度分量，而μ,κ则分别是流体的动力粘性系数和导热系数。

上式也可简化为

$\frac{\∂Q}{\∂t}=R\left(Q\right)---\left(11\right)$

根据旋转坐标系与惯性坐标系下运动关系推导，旋转坐标系下的流体力学

控制方程需要在上式右端项上添加一个源项：

$\frac{\∂Q}{\∂t}=R\left(Q\right)+S---\left(12\right)$

源项表达式为：

S＝r[0 S_x S_y S_z S_e]^T       （13）

连续性方程与惯性坐标系下一致，S_xS_yS_z为旋转坐标系相对惯性坐标系加速度的三个分量，S_e由当地速度与坐标系相对加速度点乘获得，源项展开表达为：

$S=r\left[\begin{array}{c}0\\ {(\mathrm{\ω}\×\left({-u}_{\∞}\right)+2\mathrm{\ω}\×\stackrel{\·}{b}+\mathrm{\ω}\×(\mathrm{\ω}\×b))}_x\\ {(\mathrm{\ω}\×\left({-u}_{\∞}\right)+2\mathrm{\ω}\×\stackrel{\·}{b}+\mathrm{\ω}\×(\mathrm{\ω}\×b))}_y\\ {(\mathrm{\ω}\×(-u_{\∞})+2\mathrm{\ω}\×\stackrel{\·}{b}+\mathrm{\ω}\×(\mathrm{\ω}\×b))}_z\\ b\·(\mathrm{\ω}\×(-u_{\∞})+2\mathrm{\ω}\×\stackrel{\·}{b}+\mathrm{\ω}\×(\mathrm{\ω}\×b))\end{array}\right]---\left(14\right)$

其中b＝[x y z]， $\stackrel{\·}{b}=\left[\begin{array}{ccc}u& v& w\end{array}\right]$

（3）气动力和气动力矩的积分计算：

在第（2）步求解每个状态的流场后，对物面压力和粘性应力在表面上积分得到气动力F，表面力向质心取矩并积分得到作用在质心的气动力矩M，具体计算公式如下：

                          (15)

式中，p为物面压力，I为单元矩阵，τ为粘性应力张量，其分量形式见式（1）；r为力矩参考点到物面积分点的矢径；ds为物面积分面元矢量。

在不同旋转速度下气动力和气动力矩的收敛曲线如图5所示。

第三步、滚转动导数的差分计算：

上一步建立了求解给定运动状态下的气动力和气动力矩的数值方法，这样分别计算两次不同旋转速度下模型的气动力和气动力矩，然后通过差分方法就可以计算得到模型的滚转动导数信息，相应的差分公式如下：

$F^{\mathrm{\ω}}=\frac{F_{\left(\mathrm{\ω}1\right)}-F_{\left(\mathrm{\ω}2\right)}}{\mathrm{\Δ\ω}}---\left(16\right)$

$M^{\mathrm{\ω}}=\frac{M_{\left(\mathrm{\ω}1\right)}-M_{\left(\mathrm{\ω}2\right)}}{\mathrm{\Δ\ω}}---\left(17\right)$

式中：F^ω表示气动力的滚转动导数，F_(ω1)表示角速度为ω1时的气动力，F_(ω2)表示角速度为ω2时的气动力，Δω表示角速度差，即Δω=ω₁-ω₂，M^ω表示气动力矩的滚转动导数，M_(ω1)表示角速度为ω1时的气动力矩，M_(ω2)表示角速度为ω2时的气动力矩。

滚转力矩对滚转角速度的导数与采用非定常强迫振动计算结果的对比如图7所示，采用本方法获取的滚转动导数（F^ω和M^ω）与非定常方法精度相当，但求解速度快一个量级以上。

Claims (4)

1.一种高可信度快速预测飞行器滚转动导数的方法，其特征在于：包括如下步骤：

第一步、模型表面网格的生成和空间网格划分；

第二步、旋转坐标系下气动参数的计算方法：

（1）将惯性系下的Navier-Stokes方程变换到旋转坐标系下；

（2）对变换后的方程进行数值求解，获得每个状态的流场；

（3）对物面压力和粘性应力在表面上积分得到气动力，通过表面力向质心取矩并积分得到作用在质心的气动力矩；

第三步、滚转动导数的差分计算：分别计算两次不同旋转速度下模型的气动力和气动力矩，然后通过差分方法计算得到模型的滚转动导数。

2.根据权利要求1所述的一种高可信度快速预测飞行器滚转动导数的方法，其特征在于：第一步所述的模型表面网格的生成方法采用多块对接网格技术，以保证网格质量。

3.根据权利要求1所述的一种高可信度快速预测飞行器滚转动导数的方法，其特征在于：第一步所述的空间网格划分方法为在物面附近划分边界层网格，以精确捕捉边界层内的信息。

4.根据权利要求1所述的一种高可信度快速预测飞行器滚转动导数的方法，其特征在于：所述旋转坐标系是指和模型固结在一起的非惯性系。

Images