CN106508028B - 一种确定复杂外形飞行器超声速、高超声速有迎角颤振安全边界的方法 - Google Patents

Abstract

本发明提出一种确定复杂外形飞行器超声速、高超声速有迎角颤振安全边界的方法。该方法针对任意三维物体，假设当地气流沿物面切面上任一方向，在切面上建立垂直于物面振动引起的下洗表达式，并针对任意外形翼面和旋成体机身，建立飞行器典型部件的非定常气动力表达式，以及相应的用模态坐标表示的全套计算公式，而复杂外形的当地流通过CFD数值求解，大攻角、复杂外形干扰等都可在当地流中得到反映。本发明用于任意复杂外形飞行器，包括气流三维效应严重的大后掠翼面、旋转体外形的机身、升力体和翼身融合体等复杂外形的飞行器，计算精度与CFD/CSD耦合的数值仿真算法相当，但不要求很大的计算机硬、软件资源，计算速度提高数十倍到100多倍。

Description

一种确定复杂外形飞行器超声速、高超声速有迎角颤振安全边界的方法

技术领域

本发明涉及飞行器总体技术领域，具体涉及一种确定复杂外形飞行器超声速、高超声速有迎角颤振安全边界的方法。

背景技术

目前国内外与本发明相关的现有技术如下：

(1)确定超声速、高超声速飞行器颤振安全边界的主要关键技术之一就是如何确定引起颤振的非定常气动力。目前按不同非定常气动力算法分类的颤振计算方法有两类，一类是基于计算流体力学(CFD)与计算结构动力学(CSD)耦合的数值仿真算法。CFD/CSD耦合的数值仿真方法虽然精度高，可适合任意复杂外形飞行器，但存在需要巨大的计算机硬、软件资源，计算量大，计算效率低等问题，难以在工程设计中推广使用。另一类是传统的工程方法，主要有活塞理论、非定常细长体理论、统一升力面理论和牛顿冲击流理论。这些方法的共同特点是只能用于翼面、舵面或旋转体机身等单独部件，不能用于翼-身组合体、翼-身-舵组合体或翼-身-尾组合体等复杂外形飞行器考虑部件之间的干扰、阻滞和下洗等问题，也不能用于考虑小展弦比翼面等的三元效应，除牛顿冲击流理论外，一般只能用于零迎角或小迎角情况。牛顿冲击流理论可用于较大迎角情况，但通常只适用于马赫数Ma＞7的近似计算。因此，传统工程方法的特点是计算速度快，适合于工程设计应用，但计算精度低，使用范围受限制。

(2)基于上述两类方法存在的问题，寻找一种既能适合于较大攻角复杂外形飞行器颤振非定常气动力计算，又能快速得到计算结果的方法，是业界多年来努力的方向。自上世纪90年代开始发展起来的当地流活塞理论及其与CFD方法的结合，是较好解决上述问题的方法之一。

1990年9月，陈劲松、曹军在“空气动力学学报”第8卷第3期发表的论文“超声速和高超声速翼型非定常气动力的一种近似计算方法”中，对20世纪50年代美国人发明的“活塞理论”方法作了重要改进，提出用平均迎角下翼型表面的当地流速度代替活塞理论中的无穷远来流速度作为参考速度，用翼型表面由微振动产生的附加速度代替下洗速度，仿照活塞理论给出了微振动引起的非定常压力公式，从而使得原来只能用于超声速零迎角薄翼的活塞理论可扩展用于较厚翼型和较大迎角情况。作者将它命名为“当地流活塞理论”，但所提供的公式主要用于二维刚体翼型的非定常气动力计算，未开发用于颤振计算的方法和公式。

1995年6月，杨炳渊、宋伟力在“振动与冲击”第14卷第2期发表的论文“用当地流活塞理论计算大攻角翼面超音速颤振”中，基于陈劲松等提出的当地流活塞理论，针对作弹性振动的大迎角三维翼面，给出了三维振动翼上、下表面的非定常气动力计算公式，并且给出了便于颤振计算的压力分布表达式和用模态坐标表示的非定常广义气动力表达式，从而使当地流活塞理论推广应用于超声速大迎角三维翼面的颤振非定常气动力计算。1999年9月，杨炳渊、宋伟力在上述发明的基础上，在“上海力学”第20卷第3期发表的论文“应用当地流活塞理论的大攻角升力面颤振气动力表达式”中公开了应用当地流活塞理论建立颤振数学模型中非定常广义气动力矩阵的技术方案，即通过片条模型利用样条函数弦向积分、展向分段求和的数值积分方法和公式。

但由于之前的当地流活塞理论是基于二维翼型推导的，上述改进和发展仍然将当地速度局限在翼型平面内，尽管由于利用片条模型达到向三维翼面推广的目的，但本质上仍然属于二维理论，因此只能用于气流速度方向在翼型平面内的二维翼面，或者可近似认为气流速度方向在翼型平面内的三维翼面，对于复杂组合体外形的翼面将会产生较大的误差。不能用于气流三维效应严重的大后掠翼面，更不能用于旋转体外形的导弹弹身和其它复杂外形的飞行器。

(3)2005年9月，张伟伟、叶正寅等在“力学学报”第37卷第5期发表的论文“基于当地流活塞理论的气动弹性计算方法研究”中，通过动量定理重新推导了当地流活塞理论，得到与现有技术(2)完全等价的非定常压力计算公式，但推导过程不需要台劳级数展开，没有级数收敛条件的限制和忽略高阶项的近似，也不需要等熵假设，从而奠定了当地流活塞理论不仅可用于大攻角大厚度的二维翼面，也可用于各种三维复杂外形的理论基础，而且精度比忽略台劳级数高阶项的原始活塞理论高。据此，重新定义了当地下洗速度

式中n₀为物面变形前的外法线单位矢量，n为物面变形后的外法线单位矢量，V_b为物面振动速度。2007年2月，张陈安、张伟伟、叶正寅等在“工程力学”第24卷第2期发表的论文“基于当地流活塞理论的气动弹性稳定性分析方法研究”中基于上述当地下洗速度公式，给出了用模态坐标表示的非定常广义气动力表达式

其中

公式中当地流速度V₁允许取任意方向，振动速度V_b也允许任意方向，因此是真正的三维理论，可适用于任意外形的三维翼面。但公式中下洗速度的计算必须基于变形后的外法线方向n，要求当地流计算和非定常气动力计算前必须提供翼面在气动载荷作用下的变形数据，振型数据也要求以矢量zx_m＝a_mi+b_mj+c_mk的形式提供，存在许多困难和工作量，而且非定常广义气动力表达式中如何计算广义气动力矩阵A和B积分的方法和公式也未公开。

(4)2007年1月，张陈安、张伟伟、叶正寅等在“宇航学报”第28卷第1期发表的论文“一种基于当地流活塞理论的超音速导弹气动伺服弹性分析方法”中，报导了利用现有技术(3)提供的方法完成全弹气动弹性(开环)和全弹气动伺服弹性(闭环)的临界速度计算，但所提供的公式仍然是与(3)中已公开的公式相同的关于翼面的广义气动力计算公式，未公开任何关于机(弹)身广义气动力的计算公式。

针对上述现有技术的不足，本发明对当地流活塞理论进行改进和扩展，使之适应于任意外形的飞行物体，同时提供配套的方法和公式。在此基础上，提供了一种确定复杂外形飞行器超声速、高超声速有迎角颤振安全边界的新方法。

发明内容

针对本发明所要解决的确定复杂外形飞行器超声速、高超声速有迎角颤振安全边界的关键技术，本发明基于对现有技术当地流活塞理论的改进和扩展，提出了一种计算超声速、高超声速飞行器颤振非定常气动力的方法，所采用的技术方案通过如下步骤实现：

步骤1

针对任意三维物体，假设当地气流沿物面切面上的任一方向代替背景技术(2)中假设当地气流沿二维翼型切线方向，在切面上建立垂直于物面振动引起的下洗表达式

式中ζ＝ζ(ξ，η)表示垂直于物面的振动；坐标oξηζ为切面P上的局部坐标系，原点o为切面与曲面的相切点，ξ轴方向与当地流V_L方向重合，ζ轴与过o点的外法线方向重合(如图1)。然后利用实际结构振动方向与物面外法线方向之间的几何关系，建立当地流活塞理论计算任意外形非定常气动力的表达式。而复杂外形的当地流通过CFD数值求解Euler方程或N-S方程方法获得，大攻角、复杂外形干扰流场、下洗等都可在当地流中得到反映。

步骤2

将步骤1应用于图2所示的任意对称翼型的三维翼面，得到翼面上某点非定常压力分布：

或

式中上翼面取正号，下翼面取负号；

ρ_L——当地流密度；c_L——当地声速；q_L——当地流动压；M_L——当地流马赫数；μ_w——翼面外法线对z轴的方向余弦；V_xL，V_yL——当地流速度V_L在

x轴、y轴上的分量；M_xL，M_yL——M_L在x轴、y轴上的分量；z＝z(x，y，t)——翼面在z方向的振动扰度。

以来流为参考，便于颤振计算搜索来流临界参数的压力分布表达式为：

式中

B_w(x，y)＝(c_p(x，y)·c_ρ(x，y))^1/2

C_wX(x，y)＝c_Mx(x，y)·c_p(x，y)

C_wY(x，y)＝c_My(x，y)·c_p(x，y)

c_p(x，y)＝p_L/p_∞ c_ρ(x，y)＝ρ_L/ρ_∞

c_Mx(x，y)＝M_xL/M_∞c_My(x，y)＝M_yL/M_∞

步骤3

将步骤1应用于图3所示的轴对称旋成体机身，得到机身表面的非定常压力分布：

式中μ_b＝μ_b(x)是旋成体机身在xoz平面上的母线方程r＝R(x)的外法线对z轴的方向余弦，z＝z(x，t)是机身在z方向的振动，其余符号同技术解决方案(1)。

以来流为参考，便于颤振计算搜索来流临界参数的压力分布表达式为：

式中

B_b(x，θ)＝(c_p(x，θ)·c_ρ(x，θ))^1/2

C_bX(x，θ)＝c_Mx(x，θ)·c_p(x，θ)

c_p(x，θ)＝p_L/p_∞，c_ρ(x，θ)＝ρ_L/ρ_∞，c_Mx(x，θ)＝M_xL/M_∞

步骤4

为进行结构动力学仿真，利用模态叠加法或分枝模态法对运动方程解耦，得到解耦或准解耦的运动微分方程

式中[M]、[C]和[K]分别为结构的模态质量、模态阻尼和模态刚度矩阵；[Ф]为结构模态或分枝模态行阵；{ξ}为对应的模态坐标列阵；{Q_A}为广义气动力列阵

步骤5

利用步骤2和步骤4结合，得到用模态坐标表示的任意对称翼型三维翼面的广义气动力列阵

式中

步骤6

利用步骤3和步骤4结合，得到用模态坐标表示的任意旋转体外形机身的广义气动力列阵

式中

步骤7

针对任意对称翼型三维翼面(含升力体)，如图4所示沿展向将机翼等分成n个片条，各片条中心弦上等分成m小段，机翼分成m×n个小段，利用三次样条函数按各个片条弦向积分、展向分段求和，得到步骤5中的气动力系数矩阵[B^(w)]、[C_x ^(w)]、[C_y ^(w)]中各元素的计算公式

式中A_p、A_q是离散振型沿片条中心弦分段无因次坐标拟合的样条函数三次多项式系数，B_q是离散振型沿展向线(见图5)分段无因次坐标拟合的样条函数三次多项式系数，下标i、j表示对应于按分枝顺序编号的第i阶、第j阶模态，上标(k_w，l_w)中k_w为展向片条的编号，l_w为每片条中心弦上分段的编号。和为各分段内μ_w(x)和B_w(x)、C_wX(x)、C_wY(x)的均值，和为第k_w片条宽度和片条中心弦分段的长度，是各分段线的掠角(见图5)。

步骤8

针对任意旋转体外形机身，沿弹身轴向x等分为m-1段(共m个结点)，圆周方向右半周沿圆周顺弹身子午线划分为n-1区。利用三次样条函数对轴线上的离散振型拟合分段积分后求和，得到步骤6中的气动力系数矩阵[B^(b)]、[C^(b)]中各元素的计算公式

式中A_p、A_q是离散振型沿机身轴线分段无因次坐标拟合的样条函数三次多项式系数，L_b为机身轴线总长，下标i、j表示对应于机身的第i阶、第j阶模态，上标(k，l)中k为圆周向分区的编号，l为弹身轴线上分段的编号。R^(l)、分别为各分段内R(x)、μ_b(x)、B_b(x，θ)、C_bX(x，θ)的平均值。

与现有技术比较，本发明的技术效果如下：

(1)由于本发明假设当地气流沿物面切面上的任一方向代替背景技术(2)中假设当地气流沿二维翼型切线方向，在切面上建立垂直于物面振动引起的下洗表达式，而当地流参数又通过求解适合于任意外形的Euler方程的CFD数值方法得到，因此本发明改变了背景技术(2)只能用于气流速度方向在翼型平面内的二维翼面，或者可近似认为气流速度方向在翼型平面内的三维翼面的状态，可适合于任意复杂外形飞行器，包括气流三维效应严重的大后掠翼面、旋转体外形的导弹弹身、升力体和翼身融合体等复杂外形的飞行器。

(2)同样基于上述假设，本发明下洗速度计算建立在物面切面上的局部坐标系oξηζ上，不必考虑物面在气动载荷作用下的变形，当地流和下洗速度的计算简单方便；不象背景技术(3)那样基于变形后的外法线方向，需要计算物面在气动载荷作用下的变形。

(3)本发明提供了适合于任意外形翼面或升力体、任意旋转体外形机身包括非定常气动力压力分布、广义气动力和广义气动力系数矩阵计算的成套技术方案和便于编程的计算公式，可直接用于工程计算和软件编制。

(4)本发明与背景技术(1)中基于CFD/CSD耦合的数值仿真算法比较，同样可适合任意复杂外形飞行器，但不需要很大的计算机硬、软件资源，不存在计算量大，计算效率低等问题，易于在工程设计中推广使用。经具体实施方式实施案例二的计算结果比较(图6)，本发明的计算结果与CFD/CSD耦合的数值仿真算法接近，具有相同的计算精度，但计算速度提高数十倍到100多倍。

(5)本发明与背景技术(3)的当地流活塞理论方法比较，具有相同的计算精度(图6)，但避免了当地流计算和非定常气动力计算必须基于气动载荷作用下的物面变形所造成的困难和工作量，而且公开了全套可直接用于工程设计的计算方案和公式。

本发明针对现有技术的不足，对当地流活塞理论进行改进和扩展，使之适应于任意外形的飞行物体工程计算。

附图说明

图1任意曲面上的局部坐标

图2翼面上的法向下洗

图3机身的柱坐标系统

图4机翼单元划分

图5单个分段示意图

图6A实施方案二对称分枝颤振安全边界计算结果比较图

图6B实施方案二反对称分枝颤振安全边界计算结果比较图

图7轴对称布局导弹

图8时域仿真计算流程

图9面对称布局飞行器

具体实施方式

本发明基于对现有技术当地流活塞理论的改进和扩展，提出了一种计算超声速、高超声速飞行器颤振非定常气动力的方法，所采用的技术方案如图1所示ζ＝ζ(ξ，η)表示垂直于物面的振动；坐标oξηζ为切面P上的局部坐标系，原点o为切面与曲面的相切点，ξ轴方向与当地流V_L方向重合，ζ轴与过o点的外法线方向重合；如图2所示为任意对称翼型的三维翼面，得到翼面上某点非定常压力分布；通过图3所示的轴对称旋成体机身，得到机身表面的非定常压力分布：针对任意对称翼型三维翼面如图4所示沿展向将机翼等分成n个片条，各片条中心弦上等分成m小段，机翼分成m×n个小段，利用三次样条函数按各个片条弦向积分、展向分段求和，得到如图5所示的分段无因次坐标拟合的样条函数三次多项式系数，通过实施方案二对称分枝颤振安全边界计算得到图6A、图6B结果比较图，其实现步骤如下：

步骤1

针对任意三维物体，假设当地气流沿物面切面上的任一方向代替背景技术(2)中假设当地气流沿二维翼型切线方向，在切面上建立垂直于物面振动引起的下洗表达式

式中ζ＝ζ(ξ，η)表示垂直于物面的振动；坐标oξηζ为切面P上的局部坐标系，原点o为切面与曲面的相切点，ξ轴方向与当地流V_L方向重合，ζ轴与过o点的外法线方向重合。然后利用实际结构振动方向与物面外法线方向之间的几何关系，建立当地流活塞理论计算任意外形非定常气动力的表达式。而复杂外形的当地流通过CFD数值求解Euler方程或N-S方程方法获得，大攻角、复杂外形干扰流场、下洗等都可在当地流中得到反映。

步骤2

将步骤1应用于图2所示的任意对称翼型的三维翼面，得到翼面上某点非定常压力分布：

或

式中 上翼面取正号，下翼面取负号；

ρ_L——当地流密度；c_L——当地声速；q_L——当地流动压；M_L——当地流马赫数；μ_w——翼面外法线对z轴的方向余弦；V_xL，V_yL——当地流速度V_L在

x轴、y轴上的分量；M_xL，M_yL——M_L在x轴、y轴上的分量；z＝z(x，y，t)——翼面在z方向的振动扰度。

以来流为参考，便于颤振计算搜索来流临界参数的压力分布表达式为：

式中

B_w(x，y)＝(c_p(x，y)·c_ρ(x，y))^1/2

C_wX(x，y)＝c_Mx(x，y)·c_p(x，y)

C_wY(x，y)＝c_My(x，y)·c_p(x，y)

c_p(x，y)＝p_L/p_∞ c_ρ(x，y)＝ρ_L/ρ_∞

c_Mx(x，y)＝M_xL/M_∞ c_My(x，y)＝M_yL/M_∞

步骤3

将步骤1应用于图3所示的轴对称旋成体机身，得到机身表面的非定常压力分布：

式中μ_b＝μ_b(x)是旋成体机身在xoz平面上的母线方程r＝R(x)的外法线对z轴的方向余弦，z＝z(x，t)是机身在z方向的振动，其余符号同技术解决方案(1)。

以来流为参考，便于颤振计算搜索来流临界参数的压力分布表达式为：

式中

B_b(x，θ)＝(c_p(x，θ)·c_ρ(x，θ))^1/2

C_bX(x，θ)＝c_Mx(x，θ)·c_p(x，θ)

c_p(x，θ)＝p_L/p_∞，c_ρ(x，θ)＝ρ_L/ρ_∞，c_Mx(x，θ)＝M_xL/M_∞

步骤4

为进行结构动力学仿真，利用模态叠加法或分枝模态法对运动方程解耦，得到解耦或准解耦的运动微分方程

式中[M]、[C]和[K]分别为结构的模态质量、模态阻尼和模态刚度矩阵；[Ф]为结构模态或分枝模态行阵；{ξ}为对应的模态坐标列阵；{Q_A}为广义气动力列阵

步骤5

利用步骤2和步骤4结合，得到用模态坐标表示的任意对称翼型三维翼面的广义气动力列阵

式中

步骤6

利用步骤3和步骤4结合，得到用模态坐标表示的任意旋转体外形机身的广义气动力列阵

式中

步骤7

针对任意对称翼型三维翼面(含升力体)，如图4所示沿展向将机翼等分成n个片条，各片条中心弦上等分成m小段，机翼分成m×n个小段，利用三次样条函数按各个片条弦向积分、展向分段求和，得到步骤5中的气动力系数矩阵[B^(w)]、[C_x ^(w)]、[C_y ^(w)]中各元素的计算公式

式中A_p、A_q是离散振型沿片条中心弦分段无因次坐标拟合的样条函数三次多项式系数，B_q是离散振型沿展向线(见图5)分段无因次坐标拟合的样条函数三次多项式系数，下标i、j表示对应于按分枝顺序编号的第i阶、第j阶模态，上标(k_w，l_w)中k_w为展向片条的编号，l_w为每片条中心弦上分段的编号。和为各分段内μ_w(x)和B_w(x)、C_wX(x)、C_wY(x)的均值，和为第k_w片条宽度和片条中心弦分段的长度，是各分段线的掠角(见图5)。

步骤8

针对任意旋转体外形机身，沿弹身轴向x等分为m-1段(共m个结点)，圆周方向右半周沿圆周顺弹身子午线划分为n-1区。利用三次样条函数对轴线上的离散振型拟合分段积分后求和，得到步骤6中的气动力系数矩阵[B^(b)]、[C^(b)]中各元素的计算公式

式中A_p、A_q是离散振型沿机身轴线分段无因次坐标拟合的样条函数三次多项式系数，L_b为机身轴线总长，下标i、j表示对应于机身的第i阶、第j阶模态，上标(k，l)中k为圆周向分区的编号，l为弹身轴线上分段的编号。R^(l)、分别为各分段内R(x)、μ_b(x)、B_b(x，θ)、C_bX(x，θ)的平均值。

实施案例一：

图7所示轴对称布局导弹，考虑在XOZ平面内的弯曲振动。按照分枝模态法，将弹体结构分成弹身、弹翼和尾翼三个分枝。弹体结构的运动可表示为三个分枝模态振型的叠加

Z(x，y，t)＝[Φ]{ξ}

式中

[Ф]＝[Ф_b(x，y)|Ф_w(x，y)|Ф_t(x，y)]

是由弹身分枝n_b个模态、弹翼分枝n_w个模态、尾翼分枝n_t个模态振型组成的行阵，{ξ}为对应的模态坐标列阵。由拉格朗日方程可导出结构的运动微分方程为

式中[M]、[C]和[K]分别为结构的模态质量、模态阻尼和模态刚度矩阵；{Q_A}为弹性振动产生的广义气动力列阵，{Q_C}为舵偏产生的广义控制力列阵。

利用本发明提供的方法分别计算弹身、弹翼和尾翼(舵面)用模态坐标表示的广义气动力

将舵偏角α(t)等效为垂直于舵面的位移Z_r(x，y，t)＝α(t)(x-x_轴)计算舵面的控制力{Q_C}。利用状态空间法将上述运动微分方程和控制系统方程转换为状态方程，即可通过频域的复特征值方法或时域仿真方法求解颤振和气动伺服弹性的稳定性边界参数。时域仿真方法的流程如图8所示。

实施案例二：

图9所示为面对称布局飞行器，利用振型对称性和反对称性的特点，按照分枝模态法将飞行器分为机身机翼对称振型分枝、机身机翼反对称振型分枝和垂尾分枝。根据三个分枝的划分，飞行器的运动表示为三个分枝模态振型的叠加

式中

[Ф_s]、[Ф_a]、[Ф_t]分别是由对称分枝n_s个模态、反对称分枝n_a个模态、垂尾分枝n_t个模态振型函数组成的矩阵；为对应的模态坐标矢量。忽略结构阻尼，由拉格朗日方程可导出结构的运动微分方程

式中，[M]、[C]、[K]分别为结构的模态质量、模态阻尼和模态刚度矩阵，{Q_A}为广义气动力矢量，由机身、机翼、垂尾各部分的广义气动力组成

{Q_C}为副翼舵偏转产生的广义控制力列阵。利用结构对称、振型对称反对称特性和正交性特性，可将运动微分方程解耦为如下两个独立方程组：

利用本发明提供的方法分别计算机身、机翼和垂尾翼用模态坐标表示的广义气动力和副翼舵的广义控制力，考虑机身和机翼的几何关系进行矢量叠加，可分别得到两个独立方程组右边的广义力矢量；定义状态向量对应于对称振型分枝的第一个方程组，取对应于反对称分枝的第二个方程组，取则两个独立方程组)可分别转化为如下形式的状态方程

S＝(2q_∞/M_∞)，R＝S/(M_∞c_∞)

注意式中[M]、[B]、[C]分别采用与两个独立方程相对应的下标。同样可通过频域的复特征值方法或时域仿真方法求解颤振和气动伺服弹性的稳定性边界参数。

Claims (2)

1.一种确定复杂外形飞行器超声速、高超声速有迎角颤振安全边界的方法，其特征在于：

步骤1

针对任意三维物体，假设当地气流沿物面切面上的任一方向，在切面上建立垂直于物面振动引起的下洗表达式

$W_m=\frac{\∂}{\∂t}\mathrm{\ζ}(\mathrm{\ξ},\mathrm{\η})+V_L\frac{\∂}{\∂\mathrm{\ξ}}\mathrm{\ζ}(\mathrm{\ξ},\mathrm{\η})$

式中ζ＝ζ(ξ，η)表示垂直于物面的振动；坐标O ξ η ζ为切面P上的局部坐标系，原点O为切面与曲面的相切点，ξ轴方向与当地流V_L方向重合，ζ轴与过O点的外法线方向重合；

步骤2

将步骤1应用于任意对称翼型的三维翼面，得到翼面上某点非定常压力分布：

$\mathrm{\Δ}p(x,y,t)=\±{\mathrm{\ρ}}_L{c}_L{\mathrm{\μ}}_w(\frac{\∂}{\∂t}+V_{xL}\frac{\∂}{\∂x}+V_{yL}\frac{\∂}{\∂y})Z(x,y,t)$

或

$\mathrm{\Δ}p(x,y,t)=\±\frac{2q_L{\mathrm{\μ}}_w}{M_L^2{c}_L}(\frac{\∂}{\∂t}+M_{xL}{c}_L\frac{\∂}{\∂x}+M_{yL}{c}_L\frac{\∂}{\∂y})Z(x,y,t)$

式中上翼面取正号，下翼面取负号；

ρ_L——当地流密度；c_L——当地声速；q_L——当地流动压；M_L——当地流马赫数；

μ_w——翼面外法线对Z轴的方向余弦；V_xL，V_yL——当地流速度V_L在x轴、y轴上的分量；M_xL，M_yL——M_L在x轴、y轴上的分量；z＝z(x，y，t)——翼面在Z方向的振动扰度，

步骤3

将步骤1应用于轴对称旋成体机身，得到机身表面的非定常压力分布：

$\mathrm{\Δ}p(x,\mathrm{\θ},t)=\frac{2q_L}{M_L^2{c}_L}{\mathrm{\μ}}_b\left(x\right)sin\mathrm{\θ}(\frac{\∂}{\∂t}+M_{xL}{C}_L\frac{\∂}{\∂x})Z(x,t)$

式中μ_b＝μ_b(x)是旋成体机身在XOZ平面上的母线方程r＝R(x)的外法线对Z轴的方向余弦，z＝z(x，t)是机身在z方向的振动；

步骤4

利用模态叠加法或分枝模态法对运动方程解耦，得到解耦或准解耦的运动微分方程

$\[M\]\left\{\stackrel{\·\·}{\mathrm{\ξ}}\right\}+\[C\]\left\{\stackrel{\·}{\mathrm{\ζ}}\right\}+\[K\]\left\{\mathrm{\ξ}\right\}=\left\{Q_A\right\}$

式中[M]、[C]和[K]分别为结构的模态质量、模态阻尼和模态刚度矩阵；[Φ]为结构模态或分枝模态行阵；{ξ}为对应的模态坐标列阵；{Q_A}为广义气动力列阵

$\left\{Q_A\right\}=\underset{\mathrm{\σ}}{\∫\∫}{\[\mathrm{\Φ}\]}^T\mathrm{\Δ}p(x,y,t)\cos (\mathrm{\Δ}p,^z)d\mathrm{\σ}$

步骤5

利用步骤2和步骤4结合，得到用模态坐标表示的任意对称翼型三维翼面的广义气动力列阵

$\left\{Q_A^{\left(w\right)}\right\}=-\frac{2q_{\mathrm{\∞}}}{M_{\mathrm{\∞}}}\[\frac{1}{M_{\mathrm{\∞}}{c}_{\mathrm{\∞}}}\[B^{\left(w\right)}\]\left\{\stackrel{\·}{\mathrm{\ξ}}\right\}+(\[{C_x}^{\left(w\right)}\]+\[{C_y}^{\left(w\right)}\])\left\{\mathrm{\ξ}\right\}\]$

式中

步骤6

利用步骤3和步骤4结合，得到用模态坐标表示的任意旋转体外形机身的广义气动力列阵

$\left\{Q_A^{\left(b\right)}\right\}=-\frac{2q_{\mathrm{\∞}}}{M_{\mathrm{\∞}}}\[\frac{1}{M_{\mathrm{\∞}}{c}_{\mathrm{\∞}}}\[B^{\left(b\right)}\]\left\{\stackrel{\·}{\mathrm{\ξ}}\right\}+\[C^{\left(b\right)}\]\left\{\mathrm{\ξ}\right\}\]$

式中

$\[B^{\left(b\right)}\]=\underset{b}{\∫\∫}{{\mathrm{\μ}}_b}^2\left(x\right)B_b(x,\mathrm{\θ}){\sin }^2\mathrm{\θ}{\[{\mathrm{\Φ}}^{\left(b\right)}\]}^T\[{\mathrm{\Φ}}^{\left(b\right)}\]d\mathrm{\σ}$

$\[C^{\left(b\right)}\]=\underset{b}{\∫\∫}{{\mathrm{\μ}}_b}^2\left(x\right)C_{bX}(x,\mathrm{\θ}){\sin }^2\mathrm{\θ}{\[{\mathrm{\Φ}}^{\left(b\right)}\]}^T\frac{\∂}{\∂x}\[{\mathrm{\Φ}}^{\left(b\right)}\]d\mathrm{\σ}$

步骤7

针对任意对称翼型三维翼面，含升力体，沿展向将机翼等分成n个片条，各片条中心弦上等分成m小段，机翼分成m×n个小段，利用三次样条函数按各个片条弦向积分、展向分段求和，得到步骤5中的气动力系数矩阵[B^(w)]、[C_x ^(w)]、[C_y ^(w)]中各元素的计算公式

式中A_p、A_q是离散振型沿片条中心弦分段无因次坐标拟合的样条函数三次多项式系数，B_q是离散振型沿展向线分段无因次坐标拟合的样条函数三次多项式系数，下标i、j表示对应于按分枝顺序编号的第i阶、第j阶模态，上标(k_w，l_w)中k_w为展向片条的编号，l_w为每片条中心弦上分段的编号；和为各分段内μ_w(x)和B_w(x)、C_wX(x)、C_wY(x)的均值，和为第k_w片条宽度和片条中心弦分段的长度，是各分段线的掠角；

步骤8

针对任意旋转体外形机身，沿弹身轴向x等分为m-1段，共m个结点，圆周方向右半周沿圆周顺弹身子午线划分为n-1区；利用三次样条函数对轴线上的离散振型拟合分段积分后求和，得到步骤6中的气动力系数矩阵[B^(b)]、[C^(b)]中各元素的计算公式

$B_{ij}^{\left(b\right)}=\frac{{\mathrm{\πL}}_b}{(m-1)(n-1)}\underset{k=1}{\overset{n-1}{\Σ}}{\sin }^2\left(\frac{2k-n}{2(n-1)}\mathrm{\π}\right)\underset{l=1}{\overset{m-1}{\Σ}}{R}^{\left(l\right)}{\left({\mathrm{\μ}}_b^{\left(l\right)}\right)}^2{B}_b^{(k,l)}\underset{p=0}{\overset{3}{\Σ}}\underset{q=0}{\overset{3}{\Σ}}{\left(A_p^{\left(l\right)}\right)}_i{\left(A_q^{\left(l\right)}\right)}_j/(p+q+1)$

$C_{ij}^{\left(b\right)}=\frac{\mathrm{\π}}{(n-1)}\underset{k=1}{\overset{n-1}{\Σ}}{\sin }^2\left(\frac{2k-n}{2(n-1)}\mathrm{\π}\right)\underset{l=1}{\overset{m-1}{\Σ}}{R}^{\left(l\right)}{\left({\mathrm{\μ}}_b^{\left(l\right)}\right)}^2{C}_{bX}^{(k,l)}\underset{p=0}{\overset{3}{\Σ}}\underset{q=0}{\overset{3}{\Σ}}q{\left(A_p^{\left(l\right)}\right)}_i{\left(A_q^{\left(l\right)}\right)}_j/(p+q)$

式中A_p、A_q是离散振型沿机身轴线分段无因次坐标拟合的样条函数三次多项式系数，L_b为机身轴线总长，下标i、j表示对应于机身的第i阶、第j阶模态，上标(k，l)中k为圆周向分区的编号，l为弹身轴线上分段的编号；R^(l)、分别为各分段内R(x)、μ_b(x)、B_b(x，θ)、C_bX(x，θ)的平均值。

2.根据权利要求1所述的一种确定复杂外形飞行器超声速、高超声速有迎角颤振安全边界的方法，其特征在于：复杂外形当地流通过CFD数值求解Euler方程或N-S方程方法获得，大攻角、复杂外形干扰流场、下洗都可在当地流中得到反映。