CN103310060A - 一种跨音速极限环颤振分析方法 - Google Patents

Abstract

本发明提出了一种跨音速极限环颤振分析方法，采用描述函数等效线化处理跨音速气动力的非线性特性，将等效线化后的频域气动力系数代入颤振方程，采用频域颤振求解方法获得颤振速度和颤振频率。对于不同的极限环幅值，计算得到不同的颤振速度和颤振频率构成跨音速非线性颤振的极限环特性。本发明通过气动力描述函数将跨音速气动力的非线性特征进行等效线化，在频域内求解颤振方程，能够准确预测出给定极限环幅值下的颤振速度和颤振频率。本发明计算量适中，且易于掌握，并具有很好的鲁棒性。

Description

一种跨音速极限环颤振分析方法

技术领域

本发明涉及飞行器跨音速颤振分析技术领域，具体为一种跨音速极限环颤振分析方法。

背景技术

气动弹性力学是一门研究惯性力、弹性力和气动力相互作用的交叉学科。颤振是弹性结构体在气流中发生动力不稳定的现象，是气动弹性力学的主要问题之一。由于跨音速激波的存在，跨音速流动一般都是非线性的，也就是说跨音速颤振是一种非线性颤振问题。当机翼做微幅振动时，激波移动随机翼运动表现为一种线性关系，称这种气动力为动态线性气动力模型（Dynamical Linear Aerodynamics）；当机翼结构的运动幅值较大时，激波移动随结构运动将表现出非线性关系，称这种气动力为非线性气动力模型（Nonlinear Aerodynamics）。因此，对于机翼跨音速颤振问题，如果采用动态线性气动力模型，则只能计算（微幅运动假设下的）机翼线性颤振边界，而采用非线性气动力模型，则可以分析（依赖于运动幅值的）极限环（LCO）颤振特性。

计算流体力学（CFD）/计算结构动力学（CSD）耦合的时域方法是一种常用的跨音速颤振分析方法，由于没有微幅振动假设，因此它既能计算线性颤振边界又能分析非线性极限环颤振特性。但是CFD/CSD耦合时域方法在进行颤振分析时要计算一系列跨音速速度下的气动弹性响应，根据响应的衰减和发散情况来判定颤振与否，计算量很大，耗时很长。而且它不适合对机翼颤振特性做变结构参数分析。对于非线性颤振问题，对结构和流场的计算需要经历较长计算时间才能达到稳定的极限环，计算量更大。

基于CFD技术的气动力降阶方法（ROM）能够显著提高跨音速时域方法的计算效率，但是多数ROM都是采用动态线性气动力模型，只能用于分析线性颤振问题。基于CFD计算进行跨音速非线性气动力建模往往很困难，利用神经网络可以建立用于跨音速颤振分析的非线性气动力模型，但这种方法需要进行多次训练，对构建模型所需训练信号要求很高。

发明内容

要解决的技术问题

现有的CFD/CSD耦合的跨音速极限环颤振分析方法计算量很大，同时不适合结构的变参数分析，ROM方法不能进行跨音速极限环颤振分析，而非线性气动力建模方法由于对构建模型所需的训练信号较高，不便于掌握。为解决现有技术存在的问题，本发明提出了一种跨音速极限环颤振分析方法。

技术方案

本发明采用描述函数等效线化处理跨音速气动力的非线性特性，将等效线化后的频域气动力系数代入颤振方程，采用频域颤振求解方法获得颤振速度和颤振频率。对于不同的极限环幅值，计算得到不同的颤振速度和颤振频率构成跨音速非线性颤振的极限环特性。

本发明的技术方案为：

所述一种跨音速极限环颤振分析方法，其特征在于：采用以下步骤：

步骤1：应用拉格朗日方程建立二元机翼颤振运动微分方程：

$\left\{\begin{array}{c}m\stackrel{\·\·}{h}+S_{\mathrm{\α}}\stackrel{\·\·}{\mathrm{\α}}+K_{h}h=-\frac{1}{2}\mathrm{\ρ}{V}^2\left(2b\right)c_1\\ S_{\mathrm{\α}}\stackrel{\·\·}{h}+I_{\mathrm{\α}}\stackrel{\·\·}{\mathrm{\α}}+K_{\mathrm{\α}}\mathrm{\α}=\frac{1}{2}\mathrm{\ρ}{V}^2{\left(2b\right)}^2{c}_m\end{array}\right.$

其中b为半弦长；m为质量；ρ为空气密度；V为速度；c_l为升力系数；c_m为力矩系数；S_α为机翼对刚心的质量静矩，S_α＝mx_αb；I_α为机翼对刚心的质量惯性矩，

K_h为线弹簧刚度，

ω_h为沉浮自由度的“部分”频率；K_α为扭转弹簧刚度，

ω_α为俯仰自由度的“部分”频率；二元机翼刚心位于弦线中点后ab处，a为系数；

步骤2：引入无量纲质量

将步骤1中的二元机翼颤振运动微分方程转化矩阵形式：

$\left[M\right]\left\{\stackrel{\·\·}{\mathrm{\ξ}}\right\}+{\mathrm{\ω}}_{\mathrm{\α}}^2\left[K\right]\left\{\mathrm{\ξ}\right\}=\frac{V^2}{\mathrm{\π}{\mathrm{\μb}}^2}\left\{f\right\}$

其中 $\left\{\mathrm{\ξ}\right\}=\left\{\begin{array}{c}\frac{h}{b}\\ \mathrm{\α}\end{array}\right\}$ 为广义位移，

为沉浮位移，α为俯仰位移； $\left\{f\right\}=\left\{\begin{array}{c}{-c}_l\\ {2c}_m\end{array}\right\}$ 为广义气动力； $\left[M\right]=\left[\begin{array}{cc}1& x_{\mathrm{\α}}\\ x_{\mathrm{\α}}& r_{\mathrm{\α}}^2\end{array}\right]$ 为质量阵； $\left[K\right]=\left[\begin{array}{cc}{\left(\frac{{\mathrm{\ω}}_h}{{\mathrm{\ω}}_{\mathrm{\α}}}\right)}^2& 0\\ 0& r_{\mathrm{\α}}^2\end{array}\right]$ 为刚度阵；

步骤3：计算线性颤振边界：

取沉浮位移和俯仰位移为

$\left\{\begin{array}{c}h/b=0\\ \mathrm{\α}={\mathrm{\α}}_0\sin \left(\mathrm{\ωt}\right)\end{array}\right.$

其中α₀不大于0.5°；通过最小二乘法拟合得到由俯仰位移引起的线性广义气动力的形式为

为待定系数；则由俯仰位移引起的跨音速动态线性频域气动力系数为

取沉浮位移和俯仰位移为

$\left\{\begin{array}{c}h/b={(h/b)}_0\sin \left(\mathrm{\ωt}\right)\\ \mathrm{\α}=0\end{array}\right.$

其中(h/b)₀不大于0.05；通过最小二乘法拟合得到由沉浮位移引起的线性广义气动力的形式为

为待定系数；则由沉浮位移引起的跨音速动态线性频域气动力系数为

机翼的跨音速线性广义气动力由俯仰位移和沉浮位移所引起的线性广义气动力叠加得到

$\begin{array}{c}{\left\{\begin{array}{c}\left({-c}_l\right)\\ \left({2c}_m\right)\end{array}\right\}}_L=\left\{\begin{array}{c}{(-c_l)}^{L,h}\\ {\left(2c_m\right)}^{L,h}\end{array}\right\}+\left\{\begin{array}{c}{\left({-c}_l\right)}^{L,\mathrm{\α}}\\ {\left(2c_m\right)}^{L,\mathrm{\α}}\end{array}\right\}=\left\{\begin{array}{c}{Q}_l^{L,h}\\ Q_m^{L,h}\end{array}\right\}{(h/b)}_0+\left\{\begin{array}{c}{Q}_l^{L,\mathrm{\α}}\\ Q_m^{L,\mathrm{\α}}\end{array}\right\}{\mathrm{\α}}_0\\ =\left[\begin{array}{cc}{Q}_l^{L,h}& Q_l^{L,\mathrm{\α}}\\ Q_m^{L,h}& Q_m^{L,\mathrm{\α}}\end{array}\right]\left\{\begin{array}{c}{(h/b)}_0\\ {\mathrm{\α}}_0\end{array}\right\}=\left[Q_0\right]\left\{\begin{array}{c}{(h/b)}_0\\ {\mathrm{\α}}_0\end{array}\right\}\end{array}$

其中频域气动力系数矩阵[Q₀]为马赫数及减缩频率的函数[Q₀(Ma,k)]；将机翼的跨音速线性广义气动力代入步骤2中二元机翼颤振运动微分方程的矩阵形式内，并将矩阵形式转化到频域内，得到频域颤振方程

${-\mathrm{\ω}}^2\left[M\right]\left\{{\mathrm{\ξ}}_0\right\}+{\mathrm{\ω}}_{\mathrm{\α}}^2\left[K\right]\left\{{\mathrm{\ξ}}_0\right\}=\frac{V^2}{{\mathrm{\π\μb}}^2}\left[Q_0(\mathrm{Ma},k)\right]\left\{{\mathrm{\ξ}}_0\right\}$

对频域颤振方程通过频域颤振分析方法直接求解得到无量纲化的线性颤振速度

无量纲化的线性颤振频率及线性颤振模态 $\left\{u_L\right\}={\left\{\begin{array}{cc}{\left(\frac{h}{\mathrm{b\α}}\right)}_f^L& 1\end{array}\right\}}^T;$

步骤4：计算跨音速极限环特性：

步骤4.1：根据步骤3得到的线性颤振模态和给定的俯仰极限环幅值α_0,NL，得到相应的沉浮极限环幅值

步骤4.2：取沉浮位移和俯仰位移为

$\left\{\begin{array}{c}h/b=0\\ \mathrm{\α}={\mathrm{\α}}_{0,\mathrm{NL}}\sin \left(\mathrm{\ωt}\right)\end{array}\right.$

通过最小二乘法拟合得到由俯仰位移引起的非线性广义气动力的形式为

为待定系数；相应的由俯仰位移引起的跨音速气动力描述函数为

取沉浮位移和俯仰位移为

$\left\{\begin{array}{c}h/b={(h/b)}_{0,\mathrm{NL}}\sin \left(\mathrm{\ωt}\right)\\ \mathrm{\α}=0\end{array}\right.$

通过最小二乘法拟合得到由沉浮位移引起的非线性广义气动力的形式为

为待定系数；相应的由沉浮位移引起的跨音速气动力描述函数为

机翼的跨音速非线性广义气动力由俯仰位移和沉浮位移所引起的非线性广义气动力叠加得到

$\begin{array}{c}\left\{\begin{array}{c}\left({-c}_l\right)\\ \left({2c}_m\right)\end{array}\right\}=\left\{\begin{array}{c}{\left({-c}_l\right)}^{\mathrm{NL},h}\\ {\left({2c}_m\right)}^{\mathrm{NL},h}\end{array}\right\}+\left\{\begin{array}{c}{\left({-c}_l\right)}^{\mathrm{NL},\mathrm{\α}}\\ {\left({2c}_m\right)}^{\mathrm{NL},\mathrm{\α}}\end{array}\right\}=\left\{\begin{array}{c}{Q}_l^{\mathrm{NL},h}\\ Q_m^{\mathrm{NL},h}\end{array}\right\}{(h/b)}_{0,\mathrm{NL}}+\left\{\begin{array}{c}{Q}_l^{\mathrm{NL},\mathrm{\α}}\\ Q_m^{\mathrm{NL},\mathrm{\α}}\end{array}\right\}{\mathrm{\α}}_{0,\mathrm{NL}}\\ =\left[\begin{array}{cc}{Q}_l^{\mathrm{NL},h}& Q_l^{\mathrm{NL},\mathrm{\α}}\\ Q_m^{\mathrm{NL},h}& Q_m^{\mathrm{NL},\mathrm{\α}}\end{array}\right]\left\{\begin{array}{c}{(h/b)}_{0,\mathrm{NL}}\\ {\mathrm{\α}}_{0,\mathrm{NL}}\end{array}\right\}=\left[Q_D\right]\left\{\begin{array}{c}{(h/b)}_{0,\mathrm{NL}}\\ {\mathrm{\α}}_{0,\mathrm{NL}}\end{array}\right\}\end{array}$

其中频域气动力系数矩阵[Q_D]为马赫数、减缩频率及给定俯仰位移幅值的函数[Q_D(Ma,k,α_0,NL)]；将机翼的跨音速非线性广义气动力代入步骤2中二元机翼颤振运动微分方程的矩阵形式内，并将矩阵形式转化到频域内，得到跨音速气动力非线性颤振方程

${-\mathrm{\ω}}^2\left[M\right]\left\{{\mathrm{\ξ}}_{0,\mathrm{NL}}\right\}+{\mathrm{\ω}}_{\mathrm{\α}}^2\left[K\right]\left\{{\mathrm{\ξ}}_{0,\mathrm{NL}}\right\}=\frac{V^2}{{\mathrm{\π\μb}}^2}\left[Q_D(\mathrm{Ma},k,{\mathrm{\α}}_{0,\mathrm{NL}})\right]\left\{{\mathrm{\ξ}}_{0,\mathrm{NL}}\right\}$

对跨音速气动力非线性颤振方程通过频域颤振分析方法直接求解得到给定极限环幅值 $\left\{\begin{array}{c}{\mathrm{\α}}_{0,\mathrm{NL}}\\ {(h/b)}_{0,\mathrm{NL}}\end{array}\right.$ 下的无量纲化的颤振速度

无量纲化的颤振频率

步骤5：重复步骤4，得到不同给定极限环幅值下的无量纲化颤振速度和无量纲化颤振频率，得到考虑气动力非线性的跨音速极限环颤振特性。

有益效果

本发明通过气动力描述函数将跨音速气动力的非线性特征进行等效线化，在频域内求解颤振方程，能够准确预测出给定极限环幅值下的颤振速度和颤振频率。本发明计算量适中，且易于掌握，并具有很好的鲁棒性。

附图说明

图1二元机翼颤振系统示意图。

图2气动力系统示意图。

图3本发明所计算的极限环的速度特性（a）和频率特性（b）与前方法的对比图。

具体实施方式

下面结合具体实施例描述本发明：

参照附图1，本实施例的二元机翼模型采用NACA64A010翼型，计算马赫数为0.8，其他的结构参数如下：

$a=-0.3,x_{\mathrm{\α}}=0.25,r_{\mathrm{\α}}^2=0.75,{\mathrm{\ω}}_h/{\mathrm{\ω}}_a=0.5,\mathrm{\μ}=75$

步骤1：应用拉格朗日方程建立二元机翼颤振运动微分方程：

$\left\{\begin{array}{c}m\stackrel{\·\·}{h}+S_{\mathrm{\α}}\stackrel{\·\·}{\mathrm{\α}}+K_{h}h=-\frac{1}{2}{\mathrm{\ρV}}^2\left(2b\right)c_1\\ S_{\mathrm{\α}}\stackrel{\·\·}{h}+I_{\mathrm{\α}}\stackrel{\·\·}{\mathrm{\α}}+K_{\mathrm{\α}}\mathrm{\α}=\frac{1}{2}{\mathrm{\ρV}}^2{\left(2b\right)}^2{c}_m\end{array}\right.$

其中b为半弦长；m为质量；ρ为空气密度；V为速度；c_l为升力系数；c_m为力矩系数（抬头为正）；S_α为机翼对刚心的质量静矩，S_α＝mx_αb；I_α为机翼对刚心的质量惯性矩，

K_h为线弹簧刚度，ω_h为沉浮自由度的“部分”频率；K_α为扭转弹簧刚度，ω_α为俯仰自由度的“部分”频率；二元机翼刚心位于弦线中点后ab处，a为系数；

步骤2：引入无量纲质量

将步骤1中的二元机翼颤振运动微分方程转化矩阵形式：

$\left[M\right]\left\{\stackrel{\·\·}{\mathrm{\ξ}}\right\}+{\mathrm{\ω}}_{\mathrm{\α}}^2\left[K\right]\left\{\mathrm{\ξ}\right\}=\frac{V^2}{\mathrm{\π}{\mathrm{\μb}}^2}\left\{f\right\}$

其中 $\left\{\mathrm{\ξ}\right\}=\left\{\begin{array}{c}\frac{h}{b}\\ \mathrm{\α}\end{array}\right\}$ 为广义位移，

为沉浮位移，α为俯仰位移； $\left\{f\right\}=\left\{\begin{array}{c}{-c}_l\\ {2c}_m\end{array}\right\}$ 为广义气动力； $\left[M\right]=\left[\begin{array}{cc}1& x_{\mathrm{\α}}\\ x_{\mathrm{\α}}& r_{\mathrm{\α}}^2\end{array}\right]$ 为质量阵； $\left[K\right]=\left[\begin{array}{cc}{\left(\frac{{\mathrm{\ω}}_h}{{\mathrm{\ω}}_{\mathrm{\α}}}\right)}^2& 0\\ 0& r_{\mathrm{\α}}^2\end{array}\right]$ 为刚度阵；

本实施例中的质量阵为 $\left[\begin{array}{cc}1& 0.25\\ 0.25& 0.75\end{array}\right],$ 刚度阵为 $\left[\begin{array}{cc}0.25& 0\\ 0& 0.75\end{array}\right].$

步骤3：计算线性颤振边界：

取沉浮位移和俯仰位移为

$\left\{\begin{array}{c}h/b=0\\ \mathrm{\α}={\mathrm{\α}}_0\sin \left(\mathrm{\ωt}\right)\end{array}\right.$

其中α₀不大于0.5°，本实施例中α₀＝0.2°；通过最小二乘法拟合得到由俯仰位移引起的线性广义气动力的形式为

为待定系数；则由俯仰位移引起的跨音速动态线性频域气动力系数为

取沉浮位移和俯仰位移为

$\left\{\begin{array}{c}h/b={(h/b)}_0\sin \left(\mathrm{\ωt}\right)\\ \mathrm{\α}=0\end{array}\right.$

其中(h/b)₀不大于0.05，本实施例中(h/b)₀＝0.02；通过最小二乘法拟合得到由沉浮位移引起的线性广义气动力的形式为

为待定系数；则由沉浮位移引起的跨音速动态线性频域气动力系数为

机翼的跨音速线性广义气动力由俯仰位移和沉浮位移所引起的线性广义气动力叠加得到

$\begin{array}{c}{\left\{\begin{array}{c}\left({-c}_l\right)\\ \left({2c}_m\right)\end{array}\right\}}_L=\left\{\begin{array}{c}{(-c_l)}^{L,h}\\ {\left(2c_m\right)}^{L,h}\end{array}\right\}+\left\{\begin{array}{c}{\left({-c}_l\right)}^{L,\mathrm{\α}}\\ {\left(2c_m\right)}^{L,\mathrm{\α}}\end{array}\right\}=\left\{\begin{array}{c}{Q}_l^{L,h}\\ Q_m^{L,h}\end{array}\right\}{(h/b)}_0+\left\{\begin{array}{c}{Q}_l^{L,\mathrm{\α}}\\ Q_m^{L,\mathrm{\α}}\end{array}\right\}{\mathrm{\α}}_0\\ =\left[\begin{array}{cc}{Q}_l^{L,h}& Q_l^{L,\mathrm{\α}}\\ Q_m^{L,h}& Q_m^{L,\mathrm{\α}}\end{array}\right]\left\{\begin{array}{c}{(h/b)}_0\\ {\mathrm{\α}}_0\end{array}\right\}=\left[Q_0\right]\left\{\begin{array}{c}{(h/b)}_0\\ {\mathrm{\α}}_0\end{array}\right\}\end{array}$

在进行跨音速下的线性颤振分析时，气动力非线性效应可以忽略，频域气动力系数矩阵[Q₀]为马赫数及减缩频率的函数[Q₀(Ma,k)]；k为减缩频率，与气流速度（马赫数）的关系为 $k=\frac{\mathrm{\ωb}}{V}.$

本实施例中在马赫数0.8下计算出减缩频率依次为0.08、0.10和0.15的动态线性频域气动力系数，如表1所示：

表1马赫数0.8，动态线性的频域气动力系数

将机翼的跨音速线性广义气动力代入步骤2中二元机翼颤振运动微分方程的矩阵形式内，并将矩阵形式转化到频域内，得到频域颤振方程

${-\mathrm{\ω}}^2\left[M\right]\left\{{\mathrm{\ξ}}_0\right\}+{\mathrm{\ω}}_{\mathrm{\α}}^2\left[K\right]\left\{{\mathrm{\ξ}}_0\right\}=\frac{V^2}{{\mathrm{\π\μb}}^2}\left[Q_0(\mathrm{Ma},k)\right]\left\{{\mathrm{\ξ}}_0\right\}$

频域气动力系数[Q₀(Ma,k)]为给定马赫数下的一系列减缩频率所对应的离散数值,对频域颤振方程通过频域颤振分析方法直接求解得到无量纲化的线性颤振速度

无量纲化的线性颤振频率

及线性颤振模态 $\left\{u_L\right\}={\left\{\begin{array}{cc}{\left(\frac{h}{\mathrm{b\α}}\right)}_f^L& 1\end{array}\right\}}^T;$

本实施例中利用频域颤振求解方法求解得到：

和 ${\left(\frac{h}{\mathrm{b\α}}\right)}_f^L=6.1865+1.6642i.$

步骤4：计算跨音速极限环特性：

步骤4.1：根据步骤3得到的线性颤振模态和给定的俯仰极限环幅值α_0,NL，得到相应的沉浮极限环幅值

本实施例中，给定的俯仰极限环幅值α_0,NL＝2°,根据步骤3得到的线性颤振模态计算得到沉浮极限环幅值

步骤4.2：取沉浮位移和俯仰位移为

$\left\{\begin{array}{c}h/b=0\\ \mathrm{\α}={\mathrm{\α}}_{0,\mathrm{NL}}\sin \left(\mathrm{\ωt}\right)\end{array}\right.$

通过最小二乘法拟合得到由俯仰位移引起的非线性广义气动力的形式为

为待定系数；相应的由俯仰位移引起的跨音速气动力描述函数为

取沉浮位移和俯仰位移为

$\left\{\begin{array}{c}h/b={(h/b)}_{0,\mathrm{NL}}\sin \left(\mathrm{\ωt}\right)\\ \mathrm{\α}=0\end{array}\right.$

通过最小二乘法拟合得到由沉浮位移引起的非线性广义气动力的形式为

为待定系数；相应的由沉浮位移引起的跨音速气动力描述函数为

机翼的跨音速非线性广义气动力由俯仰位移和沉浮位移所引起的非线性广义气动力叠加得到

$\begin{array}{c}\left\{\begin{array}{c}\left({-c}_l\right)\\ \left({2c}_m\right)\end{array}\right\}=\left\{\begin{array}{c}{\left({-c}_l\right)}^{\mathrm{NL},h}\\ {\left({2c}_m\right)}^{\mathrm{NL},h}\end{array}\right\}+\left\{\begin{array}{c}{\left({-c}_l\right)}^{\mathrm{NL},\mathrm{\α}}\\ {\left({2c}_m\right)}^{\mathrm{NL},\mathrm{\α}}\end{array}\right\}=\left\{\begin{array}{c}{Q}_l^{\mathrm{NL},h}\\ Q_m^{\mathrm{NL},h}\end{array}\right\}{(h/b)}_{0,\mathrm{NL}}+\left\{\begin{array}{c}{Q}_l^{\mathrm{NL},\mathrm{\α}}\\ Q_m^{\mathrm{NL},\mathrm{\α}}\end{array}\right\}{\mathrm{\α}}_{0,\mathrm{NL}}\\ =\left[\begin{array}{cc}{Q}_l^{\mathrm{NL},h}& Q_l^{\mathrm{NL},\mathrm{\α}}\\ Q_m^{\mathrm{NL},h}& Q_m^{\mathrm{NL},\mathrm{\α}}\end{array}\right]\left\{\begin{array}{c}{(h/b)}_{0,\mathrm{NL}}\\ {\mathrm{\α}}_{0,\mathrm{NL}}\end{array}\right\}=\left[Q_D\right]\left\{\begin{array}{c}{(h/b)}_{0,\mathrm{NL}}\\ {\mathrm{\α}}_{0,\mathrm{NL}}\end{array}\right\}\end{array}$

其中频域气动力系数矩阵[Q_D]为马赫数、减缩频率及给定俯仰位移幅值的函数[Q_D(Ma,k,α_0,NL)]；

本实施例中在马赫数0.8下计算出减缩频率依次为0.08、0.10和0.15的气动力描述函数，如表2所示：

表2马赫数0.8，俯仰极限环幅值α₀＝2°气动力描述函数

将机翼的跨音速非线性广义气动力代入步骤2中二元机翼颤振运动微分方程的矩阵形式内，并将矩阵形式转化到频域内，得到跨音速气动力非线性颤振方程

${-\mathrm{\ω}}^2\left[M\right]\left\{{\mathrm{\ξ}}_{0,\mathrm{NL}}\right\}+{\mathrm{\ω}}_{\mathrm{\α}}^2\left[K\right]\left\{{\mathrm{\ξ}}_{0,\mathrm{NL}}\right\}=\frac{V^2}{{\mathrm{\π\μb}}^2}\left[Q_D(\mathrm{Ma},k,{\mathrm{\α}}_{0,\mathrm{NL}})\right]\left\{{\mathrm{\ξ}}_{0,\mathrm{NL}}\right\}$

对跨音速气动力非线性颤振方程通过频域颤振分析方法直接求解得到给定极限环幅值 $\left\{\begin{array}{c}{\mathrm{\α}}_{0,\mathrm{NL}}\\ {(h/b)}_{0,\mathrm{NL}}\end{array}\right.$ 下的无量纲化的颤振速度无量纲化的颤振频率

本实施例中利用频域颤振求解方法求解得到：

步骤5：重复步骤4，得到不同给定极限环幅值下的无量纲化颤振速度和无量纲化颤振频率，得到考虑气动力非线性的跨音速极限环颤振特性。

本实施例中计算出其他极限环幅值下的颤振速度和颤振频率，如表3所示，并与CFD/CSD耦合方法的结果比较如图3所示，本方法与前方法结果吻合很好。这说明采用气动力描述函数方法进行跨音速非线性颤振分析能够很好的反映出跨音速极限环颤振特性，是一种有效的跨音速气动力非线性颤振分析途径。

表3本方法计算得到的极限环特性

Claims (1)

1.一种跨音速极限环颤振分析方法，其特征在于：采用以下步骤：

步骤1：应用拉格朗日方程建立二元机翼颤振运动微分方程：

$\left\{\begin{array}{c}m\stackrel{\·\·}{h}+S_{\mathrm{\α}}\stackrel{\·\·}{\mathrm{\α}}+K_{h}h=-\frac{1}{2}\mathrm{\ρ}{V}^2\left(2b\right)c_1\\ S_{\mathrm{\α}}\stackrel{\·\·}{h}+I_{\mathrm{\α}}\stackrel{\·\·}{\mathrm{\α}}+K_{\mathrm{\α}}\mathrm{\α}=\frac{1}{2}\mathrm{\ρ}{V}^2{\left(2b\right)}^2{c}_m\end{array}\right.$

其中b为半弦长；m为质量；ρ为空气密度；V为速度；c_l为升力系数；c_m为力矩系数；S_α为机翼对刚心的质量静矩，S_α＝mx_αb；I_α为机翼对刚心的质量惯性矩，

K_h为线弹簧刚度，

ω_h为沉浮自由度的“部分”频率；K_α为扭转弹簧刚度，

ω_α为俯仰自由度的“部分”频率；二元机翼刚心位于弦线中点后ab处，a为系数；

步骤2：引入无量纲质量将步骤1中的二元机翼颤振运动微分方程转化矩阵形式：

$\left[M\right]\left\{\stackrel{\·\·}{\mathrm{\ξ}}\right\}+{\mathrm{\ω}}_{\mathrm{\α}}^2\left[K\right]\left\{\mathrm{\ξ}\right\}=\frac{V^2}{\mathrm{\π}{\mathrm{\μb}}^2}\left\{f\right\}$

其中 $\left\{\mathrm{\ξ}\right\}=\left\{\begin{array}{c}\frac{h}{b}\\ \mathrm{\α}\end{array}\right\}$ 为广义位移，

为沉浮位移，α为俯仰位移； $\left\{f\right\}=\left\{\begin{array}{c}{-c}_l\\ {2c}_m\end{array}\right\}$ 为广义气动力； $\left[M\right]=\left[\begin{array}{cc}1& x_{\mathrm{\α}}\\ x_{\mathrm{\α}}& r_{\mathrm{\α}}^2\end{array}\right]$ 为质量阵； $\left[K\right]=\left[\begin{array}{cc}{\left(\frac{{\mathrm{\ω}}_h}{{\mathrm{\ω}}_{\mathrm{\α}}}\right)}^2& 0\\ 0& r_{\mathrm{\α}}^2\end{array}\right]$ 为刚度阵；

步骤3：计算线性颤振边界：

取沉浮位移和俯仰位移为

$\left\{\begin{array}{c}h/b=0\\ \mathrm{\α}={\mathrm{\α}}_0\sin \left(\mathrm{\ωt}\right)\end{array}\right.$

其中α₀不大于0.5°；通过最小二乘法拟合得到由俯仰位移引起的线性广义气动力的形式为

为待定系数；则由俯仰位移引起的跨音速动态线性频域气动力系数为

取沉浮位移和俯仰位移为

$\left\{\begin{array}{c}h/b={(h/b)}_0\sin \left(\mathrm{\ωt}\right)\\ \mathrm{\α}=0\end{array}\right.$

其中(h/b)₀不大于0.05；通过最小二乘法拟合得到由沉浮位移引起的线性广义气动力的形式为

为待定系数；则由沉浮位移引起的跨音速动态线性频域气动力系数为

机翼的跨音速线性广义气动力由俯仰位移和沉浮位移所引起的线性广义气动力叠加得到

$\begin{array}{c}{\left\{\begin{array}{c}\left({-c}_l\right)\\ \left({2c}_m\right)\end{array}\right\}}_L=\left\{\begin{array}{c}{(-c_l)}^{L,h}\\ {\left(2c_m\right)}^{L,h}\end{array}\right\}+\left\{\begin{array}{c}{\left({-c}_l\right)}^{L,\mathrm{\α}}\\ {\left(2c_m\right)}^{L,\mathrm{\α}}\end{array}\right\}=\left\{\begin{array}{c}{Q}_l^{L,h}\\ Q_m^{L,h}\end{array}\right\}{(h/b)}_0+\left\{\begin{array}{c}{Q}_l^{L,\mathrm{\α}}\\ Q_m^{L,\mathrm{\α}}\end{array}\right\}{\mathrm{\α}}_0\\ =\left[\begin{array}{cc}{Q}_l^{L,h}& Q_l^{L,\mathrm{\α}}\\ Q_m^{L,h}& Q_m^{L,\mathrm{\α}}\end{array}\right]\left\{\begin{array}{c}{(h/b)}_0\\ {\mathrm{\α}}_0\end{array}\right\}=\left[Q_0\right]\left\{\begin{array}{c}{(h/b)}_0\\ {\mathrm{\α}}_0\end{array}\right\}\end{array}$

其中频域气动力系数矩阵[Q₀]为马赫数及减缩频率的函数[Q₀(Ma,k)]；将机翼的跨音速线性广义气动力代入步骤2中二元机翼颤振运动微分方程的矩阵形式内，并将矩阵形式转化到频域内，得到频域颤振方程

${-\mathrm{\ω}}^2\left[M\right]\left\{{\mathrm{\ξ}}_0\right\}+{\mathrm{\ω}}_{\mathrm{\α}}^2\left[K\right]\left\{{\mathrm{\ξ}}_0\right\}=\frac{V^2}{{\mathrm{\π\μb}}^2}\left[Q_0(\mathrm{Ma},k)\right]\left\{{\mathrm{\ξ}}_0\right\}$

对频域颤振方程通过频域颤振分析方法直接求解得到无量纲化的线性颤振速度

无量纲化的线性颤振频率

及线性颤振模态 $\left\{u_L\right\}={\left\{\begin{array}{cc}{\left(\frac{h}{\mathrm{b\α}}\right)}_f^L& 1\end{array}\right\}}^T;$

步骤4：计算跨音速极限环特性：

步骤4.1：根据步骤3得到的线性颤振模态和给定的俯仰极限环幅值α_0,NL，得到相应的沉浮极限环幅值

步骤4.2：取沉浮位移和俯仰位移为

$\left\{\begin{array}{c}h/b=0\\ \mathrm{\α}={\mathrm{\α}}_{0,\mathrm{NL}}\sin \left(\mathrm{\ωt}\right)\end{array}\right.$

通过最小二乘法拟合得到由俯仰位移引起的非线性广义气动力的形式为

为待定系数；相应的由俯仰位移引起的跨音速气动力描述函数为

取沉浮位移和俯仰位移为

$\left\{\begin{array}{c}h/b={(h/b)}_{0,\mathrm{NL}}\sin \left(\mathrm{\ωt}\right)\\ \mathrm{\α}=0\end{array}\right.$

通过最小二乘法拟合得到由沉浮位移引起的非线性广义气动力的形式为

为待定系数；相应的由沉浮位移引起的跨音速气动力描述函数为

机翼的跨音速非线性广义气动力由俯仰位移和沉浮位移所引起的非线性广义气动力叠加得到

$\begin{array}{c}\left\{\begin{array}{c}\left({-c}_l\right)\\ \left({2c}_m\right)\end{array}\right\}=\left\{\begin{array}{c}{\left({-c}_l\right)}^{\mathrm{NL},h}\\ {\left({2c}_m\right)}^{\mathrm{NL},h}\end{array}\right\}+\left\{\begin{array}{c}{\left({-c}_l\right)}^{\mathrm{NL},\mathrm{\α}}\\ {\left({2c}_m\right)}^{\mathrm{NL},\mathrm{\α}}\end{array}\right\}=\left\{\begin{array}{c}{Q}_l^{\mathrm{NL},h}\\ Q_m^{\mathrm{NL},h}\end{array}\right\}{(h/b)}_{0,\mathrm{NL}}+\left\{\begin{array}{c}{Q}_l^{\mathrm{NL},\mathrm{\α}}\\ Q_m^{\mathrm{NL},\mathrm{\α}}\end{array}\right\}{\mathrm{\α}}_{0,\mathrm{NL}}\\ =\left[\begin{array}{cc}{Q}_l^{\mathrm{NL},h}& Q_l^{\mathrm{NL},\mathrm{\α}}\\ Q_m^{\mathrm{NL},h}& Q_m^{\mathrm{NL},\mathrm{\α}}\end{array}\right]\left\{\begin{array}{c}{(h/b)}_{0,\mathrm{NL}}\\ {\mathrm{\α}}_{0,\mathrm{NL}}\end{array}\right\}=\left[Q_D\right]\left\{\begin{array}{c}{(h/b)}_{0,\mathrm{NL}}\\ {\mathrm{\α}}_{0,\mathrm{NL}}\end{array}\right\}\end{array}$

其中频域气动力系数矩阵[Q_D]为马赫数、减缩频率及给定俯仰位移幅值的函数[Q_D(Ma,k,α_0,NL)]；将机翼的跨音速非线性广义气动力代入步骤2中二元机翼颤振运动微分方程的矩阵形式内，并将矩阵形式转化到频域内，得到跨音速气动力非线性颤振方程

${-\mathrm{\ω}}^2\left[M\right]\left\{{\mathrm{\ξ}}_{0,\mathrm{NL}}\right\}+{\mathrm{\ω}}_{\mathrm{\α}}^2\left[K\right]\left\{{\mathrm{\ξ}}_{0,\mathrm{NL}}\right\}=\frac{V^2}{{\mathrm{\π\μb}}^2}\left[Q_D(\mathrm{Ma},k,{\mathrm{\α}}_{0,\mathrm{NL}})\right]\left\{{\mathrm{\ξ}}_{0,\mathrm{NL}}\right\}$

对跨音速气动力非线性颤振方程通过频域颤振分析方法直接求解得到给定极限环幅值 $\left\{\begin{array}{c}{\mathrm{\α}}_{0,\mathrm{NL}}\\ {(h/b)}_{0,\mathrm{NL}}\end{array}\right.$ 下的无量纲化的颤振速度

无量纲化的颤振频率

步骤5：重复步骤4，得到不同给定极限环幅值下的无量纲化颤振速度和无量纲化颤振频率，得到考虑气动力非线性的跨音速极限环颤振特性。

Images