一种基于气动力不确定降阶的机翼结构气动弹性稳定性分析方法

技术领域

本发明涉及机翼结构气动弹性鲁棒稳定性分析领域，特别涉及一种基于气动力不确定降阶的机翼结构气动弹性稳定性分析方法。

背景技术

气动弹性力学主要研究弹性结构在气动力、弹性力和惯性力等耦合作用下的响应和稳定性问题，它与现代飞行器技术的发展密切相关。动稳定性问题即通常所说的颤振问题是气动弹性力学领域备受关注的一个分支，也是现代飞行器设计中需要首先考虑的问题之一。从振动的观点看，颤振是弹性结构在非定常气动力作用下的一种自激振动，非定常气动力在这一过程中起了十分关键的作用。因此，建立准确高效的非定常气动力模型是开展颤振分析的重要基础。20世纪发展的基于线化理论的各种非定常气动力模型因其建模简便、计算量小而被广泛应用于工程结构的气动弹性分析中，但是这类模型并不适用于跨音速流动、大攻角飞行、气流分离等情况下的非线性气动弹性问题。

随着计算机硬件的快速发展，以跨音速小扰动方程、Euler方程、N-S方程为基础的计算流体力学(CFD)技术因其在非线性气动力预测，特别是跨音速流动模拟方面所展现的优越性而在非定常气动力计算和气动弹性分析中得到了更多的重视。基于非定常CFD技术的时域气动弹性模拟显著提升了气动弹性分析的精度和应用范围，已成为气动弹性领域的研究热点。然而，伴随着CFD技术在时空维上对流动描述越来越精细，更逼近真实物理特性，使得非定常气动力模型的维数越来越高。一般情况下，基于CFD技术的流场求解器的阶数可达10⁴～10⁷，这意味着利用CFD技术开展气动弹性研究需要耗费庞大的计算量和分析时间，一定程度上阻碍了其在系统参数设计、气动弹性优化和颤振主动抑制等方面的进一步应用。

为了克服基于CFD技术的气动弹性分析在计算效率和易设计性方面的局限性，近年来，人们一直致力于寻求高效高精度的低阶非定常气动力模型。对于通常仅涉及微幅振动的气动弹性稳定性分析而言，尽管背景流场在空间维上是非线性的，但非定常气动力关于小幅结构振动在时间维上表现为线性。依据上述动态线化假设发展起来的基于CFD技术的非定常气动力降阶模型(ROM)因其形式简单、计及流动的非线性特征、兼顾计算精度和效率而成为代 替CFD流场求解器的理想选择。根据建模思路的不同，基于CFD技术的非定常气动力ROM方法主要分为两类：基于本征正交分解(POD)技术的非定常流场降阶方法和基于系统辨识技术的模态气动力建模方法。在非定常气动力降阶方面，两类方法精度相当，且均依赖于结构的模态信息，本发明在构建非定常气动力ROM时采用了基于系统辨识技术的模态气动力建模方法。相较于POD方法，系统辨识方法着眼于气动弹性系统的输入输出关系，思路直观且应用方便。实际工程中的气动弹性系统通常是多输入多输出(MIMO)形式的，即含有多阶结构模态和多阶广义气动力，本发明充分利用自回归滑动平均(ARMA)模型在MIMO系统辨识方面的优势，将非定常CFD流场求解器视为待辨识的动力系统，以结构模态位移为输入，广义气动力为输出，构建基于ARMA模型的状态空间形式的非定常气动力ROM，并直接耦合结构运动状态空间方程实现高效高精度的气动弹性分析。

传统的气动弹性稳定性分析都是基于参数确定的标称系统展开的，而真实的气动弹性系统会受到各种不确定性因素的影响，如在系统建模过程中存在的各种假设和简化、模态截断、因机理不清而未建模等引起的物理模型的不确定性，网格质量差异、收敛精度不同、计算区域大小等导致的数值计算的不确定性以及结构、气动等物理参数的不精确性或分散性造成的系统参数的不确定性。由于这些不确定因素的存在，使得理论模型不足以准确描述真实系统的动力学行为，特别是其稳定性特性。目前，在实际工程中，通过引入颤振安全裕度来避免飞行器在飞行包线内因各种不确定性因素影响而发生颤振等失稳现象。这种对于不确定性因素一体估计的策略缺乏对不确定性的定量认识，有悖于气动弹性系统精细化分析和设计的发展趋势，甚至会因对不确定性的估计不足而导致灾难性的后果。例如，美国高超声速飞行器X-43A在第一次试飞中正是由于气动设计过程中对于不确定性的模拟不足致使控制系统过高估计了设计冗余而造成失控。因此，合理准确的气动弹性系统不确定性建模很大程度上决定了系统的不确定性颤振边界，是开展不确定性颤振分析的关键。

目前，定量考虑不确定性对于气动弹性动稳定性影响的分析方法主要有两类，即概率颤振分析方法和非概率颤振分析方法。概率颤振分析将不确定性量处理成满足某种概率分布的随机变量，目的是获得概率意义上的偏乐观的“软”的稳定边界，在这个边界内不能保证气动弹性系统绝对安全。概率颤振分析的主要弊端在于其过分依赖于不确定性量的先验信息，需要通过大量的样本实验事先获得不确定性量的分布规律。非概率颤振分析仅需知道不确定性量的边界信息，将不确定性量定量化为不确定但有界变量，可实现贫信息、少数据条件下的不确定性影响分析。从数学的观点看，含不确定性的气动弹性系统已由一个单一的确定的系统转化为一个系统的集合，不确定性量的大小决定了集合的边界，进而确定了系统的稳定性边界。相较于概率颤振分析，非概率颤振分析获得的是一个偏保守的“硬”的鲁棒稳定边界， 在这个边界内能保证气动弹性系统绝对安全。本发明采用非概率颤振分析方法，综合考虑气动力辨识环节的不确定性因素影响，通过将其统一定量化为辨识模型中的不确定但有界的区间噪声序列进行气动弹性系统的不确定性建模和鲁棒稳定性分析。

发明内容

本发明要解决的技术问题为：提供了一种基于非定常气动力不确定降阶的气动弹性系统不确定性建模技术及其鲁棒稳定性边界分析方法。该方法以基于CFD技术的非定常气动力模型降阶方法为基础，综合考虑气动力辨识过程中数值计算和气动参数的不确定性，将其统一定量化为辨识模型中的不确定但有界区间噪声序列，借助区间集员辨识算法实现气动力模型的不确定性降阶。所提供的气动弹性系统不确定性建模思路和稳定性边界预测技术兼顾了计算效率、分析精度和系统鲁棒性。

本发明解决上述技术问题采用的技术方案为：一种基于气动力不确定降阶的机翼结构气动弹性稳定性分析方法，包括以下步骤：

(1)建立机翼结构CSD分析模型并进行模态分析，提取机翼结构各有限元结点归一化后的模态位移信息；

(2)将机翼结构表面作为气动结构耦合交界面，建立机翼结构非定常气动力CFD分析模型，根据机翼结构模态位移信息生成时间历程形式的“3211”位移信号，提取机翼气动结构耦合交界面网格结点的变形时间历程，并据此进行CFD非定常气动力求解器数据训练，获得给定马赫数条件下CFD求解过程的输入，即模态位移历程ξ(k)，和输出，即模态气动力系数历程f_a(k)；

(3)综合考虑气动力辨识过程中数值计算和气动参数的不确定性，将其统一定量化为不确定但有界噪声序列e(k)∈e^I(k)＝[-ω(k),ω(k)]，分别将步骤(2)中的模态位移历程和模态气动力系数历程作为输入和输出，建立离散时间形式的含区间噪声的非定常气动力ARMA辨识模型：

其中，f_a(k)是系统输出量的第k次观测值，为p维列向量；ξ(k)是系统输入量的第k次观测值，为q维列向量；e(k)为p维区间噪声序列的第k次观测值；A_i和B_j为待辨识的系统参数矩阵；na和nb分别为输出和输入的延迟级数，θ^T＝[A₁…A_naB₀…B_nb-1]为>a^T(k-1),…,f_a^T(k-na),ξ^T(k),…,ξ^T(k-nb+1)]^T为na·p+nb·q维回归向量；

(4)利用区间数学和集员辨识思想，寻求与训练数据序列{f_a(k),x(k)}和噪声序列{e(k)}相容的系统参数的最小超长方体，给出辨识参数的区间估计和获得离散时间形式的非定常气动力不确定降阶模型，即

其中，

$x_a\left(k\right)={[{f_a}^T(k-1),...,{f_a}^T(k-na),{\xi }^T(k-1),...,{\xi }^T(k-nb+1)]}^T$

$D_a^I=\left[B_0^I\right]$

(5)用步骤(4)中的不确定降阶模型代替步骤(2)中的CFD求解器，并耦合由步骤(1)中机翼结构CSD分析模型提取的结构运动状态方程，建立离散时间形式的不确定气动弹性系统状态空间模型，即

其中，

式(18)中，q为来流动压，为结构状态变量，A_s、B_s、C_s和D_s均为离散时间域内结构运动状态空间方程的系数矩阵。

(6)调整来流动压q，计算该动压条件下气动弹性系统区间状态矩阵特征值实部和虚部的上、下界，即：

其中，和分别为区间状态矩阵取名义值时第i阶特征值的实部和虚部，和分别为区间状态矩阵取名义值时与第i阶特征值所对应的右特征向量的实部和虚部，ΔA_as为的区间半径矩阵，并据此在复平面内绘制不确定气动弹性系统随来流动压变化的根轨迹图；

(7)判断是否完成不确定气动弹性系统根轨迹分析，若未完成，转到步骤(6)，若完成，则由根轨迹穿越复平面单位圆的临界点预测不确定气动弹性系统颤振速度因子的上下界，在获得给定来流动压q下区间状态矩阵特征值实部和虚部的范围后，便可通过优化方法确定区间状态矩阵谱半径的上、下界，即：

当时，不确定气动弹性系统完全稳定；当而时，不确定气动弹性系统不完全稳定；当ρ(A_as)>1时，不确定气动弹性系统完全不稳定。使和ρ(A_as)＝1的来流动压分别为不确定气动弹性系统由完全稳定变为不完全稳定的临界动压q和由不完全稳定变为完全不稳定的临界动压其分别对应不确定气动弹性系统颤振速度因子的下界和上界

(8)判断是否完成全马赫数条件下不确定气动弹性系统的颤振速度边界估计，若未完成，调整计算马赫数，重复步骤(2)～(7)，若完成，给出不确定气动弹性系统颤振速度因子上、 下界随马赫数的变化情况，由此识别不确定气动弹性系统的完全稳定域、不完全稳定域和完全不稳定域，预测不确定气动弹性系统的颤振速度边界，完成不确定气动弹性系统的稳定性分析；

(9)由步骤(6)中已知来流动压q条件下的气动弹性系统区间状态矩阵，还可直接建立不确定气动弹性系统的鲁棒稳定性快速判据，即：

$\underset{i,j=1}{\overset{n}{\Sigma }}{\mathrm{\Delta a}}_{asij}{\sigma }_{max}\left(P_{ij}\right)<1---\left(21\right)$

其中，σ_max(B)表示矩阵B的最大奇异值；Δa_asij为不确定气动弹性系统区间状态矩阵第i行、第j列元素的区间半径；(P为正定对称矩阵，是方程的解；E_ij表示第i行、第j列的元素为1，其他元素为0的[na×p+(nb+1)×q]×[na×p+(nb+1)×q]维矩阵)；

(10)判断气动弹性系统区间状态矩阵是否满足步骤(9)中的鲁棒稳定性快速判据条件，若不满足，则不确定气动弹性系统完全不稳定或不完全稳定，若满足，则不确定气动弹性系统稳定。

本发明与现有技术相比的优点在于：本发明提供了一种气动弹性系统不确定性建模的新思路，在构建高效高精度的基于CFD技术的非定常气动力降阶模型的过程中综合考虑存在于数值计算和气动参数的不确定性，将其统一定量化为辨识模型中的不确定但有界区间噪声序列，并借助区间集员辨识算法建立不确定性气动力降阶模型，兼顾了气动力模型的精度、计算效率和鲁棒性。同时，将不确定性气动力降阶模型与结构模型相耦合，构建了状态空间形式的不确定性气动弹性系统的数学模型，提供了一种从区间状态矩阵特征值角度出发的预测系统鲁棒稳定性边界的高效方法。本发明所提出的基于CFD技术的不确定性气动力降阶模型构建技术和含区间参数的气动弹性系统鲁棒稳定性分析方法都具有工程实用价值。

附图说明

图1为区间集员辨识算法示意图；

图2为二元Isogai机翼气动弹性模型示意图；

图3为“3211”模态位移训练输入信号图；

图4为“3211”信号输入下模态气动力系数的CFD训练输出和降阶模型输出对比图；

图5为由辨识参数上下界构成的超长方体体积收敛历程图；

图6为不确定气动弹性系统根轨迹图，其中，图6(a)为Ma＝0.75时不确定气动弹性系统的根轨迹图，图6(b)为图6(a)中框选区域的根轨迹放大图；

图7为本发明的方法实现流程图。

具体实施方式

本实例以图2所示的二元Isogai机翼为对象，利用本发明提出的一种基于气动力不确定降阶的机翼结构气动弹性稳定性分析方法对其进行稳定性分析，如图7所示，包括以下步骤：

(1)建立Isogai机翼结构CSD分析模型，该机翼为一个后掠三维机翼的外伸截面，采用NACA 64A010翼型，具有沉浮h(向下为正)和俯仰α(抬头为正)两个自由度，具体结构参数为：b＝0.5m，x_α＝1.8，a＝-2，ω_h/ω_α＝1，μ＝60，其中，b为半弦长，x_α、a分别为弹性轴与翼弦中点(弹性轴位于中点后方时为正)、弹性轴与质心间的无量纲距离，r_α为机翼对弹性轴的无量纲回转半径，ω_h、ω_α分别为沉浮和俯仰模态的解耦固有频率，为质量比，对机翼结构进行模态分析并提取各有限元结点归一化后的模态位移信息；

(2)将机翼结构表面作为气动结构耦合交界面，建立机翼结构非定常气动力CFD分析模型，生成CFD求解器的训练输入信号，即机翼结构“3211”模态位移ξ(k)历程(如图3所示)，本实例中对应为沉浮和俯仰模态位移，并以此为输入进行CFD求解器数据训练，获得给定马赫数(本实例中选取的来流马赫数为0.75)条件下的训练输出，即模态气动力系数f_a(k)历程(如图4所示)，本实例中对应为广义升力系数和广义力矩系数；

(3)综合考虑气动力辨识过程中CFD数值计算和气动参数的不确定性，将其统一定量化为不确定但有界噪声序列e(k)∈e^I(k)＝[-ω(k),ω(k)]，本实例中噪声序列的区间半径ω(k)取为气动力系数绝对值的1％，分别将步骤(2)中的模态位移历程和模态气动力系数历程作为输入和输出，建立离散时间形式的含区间噪声的非定常气动力ARMA辨识模型，即

$f_a\left(k\right)=\underset{i=1}{\overset{na}{\Sigma }}{A}_i{f}_a(k-i)+\underset{j=0}{\overset{nb-1}{\Sigma }}{B}_j\xi (k-j)+e\left(k\right)={\theta }^{T}x\left(k\right)+e\left(k\right)---\left(22\right)$

其中，x(k)＝[f_a^T(k-1),…,f_a^T(k-na),ξ^T(k),…,ξ^T(k-nb+1)]^T，本实例中，na和nb均取3；

(4)利用区间数学和集员辨识思想，寻求与训练数据序列{f_a(k),x(k)}和噪声序列{e(k)}相容的系统参数的最小超长方体，给出如表1所示的辨识参数的区间估计和由辨识参数上下界构成的超长方体体积的收敛历程如图5所示，根据参数辨识结果建立离散时间形式的非定常气动力不确定降阶模型，即：

$\begin{array}{c}{x}_a(k+1)=A_a^I{x}_a\left(k\right)+B_a^I\xi \left(k\right)\\ f_a\left(k\right)=C_a^I{x}_a\left(k\right)+D_a^I\xi \left(k\right)\end{array}---\left(23\right)$

表1 不确定性气动力降阶模型参数的区间估计

(5)用步骤(4)中的不确定降阶模型代替步骤(2)中的CFD求解器，并耦合由步骤(1)中机翼结构CSD分析模型提取的结构运动状态方程，建立离散时间形式的不确定气动弹性系统状态空间模型，即：

$\begin{array}{c}{x}_{as}(k+1)=A_{as}^I{x}_{as}\left(k\right)\\ \xi \left(k\right)=C_{as}{x}_{as}\left(k\right)\end{array}---\left(24\right)$

(6)调整来流动压q，计算该动压条件下气动弹性系统区间状态矩阵特征值实部和虚部的上、下界，即：

$\begin{array}{c}{\underset{‾}{\lambda }}_{ir}={\lambda }_{ir}^c-\left|v_{ir}^c{|}^T{\mathrm{\Delta A}}_{as}\right|v_{ir}^c|-|v_{im}^c{|}^T{\mathrm{\Delta A}}_{as}\left|v_{im}^c\right|\\ {\bar{\lambda }}_{ir}={\lambda }_{ir}^c+\left|v_{ir}^c{|}^T{\mathrm{\Delta A}}_{as}\right|v_{ir}^c|+|v_{im}^c{|}^T{\mathrm{\Delta A}}_{as}\left|v_{im}^c\right|\\ {\underset{‾}{\lambda }}_{im}={\lambda }_{im}^c-\left|v_{ir}^c{|}^T{\mathrm{\Delta A}}_{as}\right|v_{im}^c|-|v_{im}^c{|}^T{\mathrm{\Delta A}}_{as}\left|v_{ir}^c\right|\\ {\bar{\lambda }}_{im}={\lambda }_{im}^c+\left|v_{ir}^c{|}^T{\mathrm{\Delta A}}_{as}\right|v_{im}^c|+|v_{im}^c{|}^T{\mathrm{\Delta A}}_{as}\left|v_{ir}^c\right|\end{array}---\left(25\right)$

并据此在复平面内绘制给定马赫数下(本例中选取的来流马赫数为0.75)不确定气动弹性系统随来流动压变化的根轨迹图，如(b)

图6所示，不确定气动弹性系统在给定来流动压的特征值为由其上下界围成的一个矩形区域， 而由此生成的根轨迹呈带状；

(7)判断是否完成不确定气动弹性系统根轨迹分析，若未完成，转到步骤(6)，若完成，则由根轨迹穿越复平面单位圆的临界点预测不确定气动弹性系统颤振速度因子的上下界，在获得给定来流动压q下区间状态矩阵特征值实部和虚部的范围后，便可通过优化方法确定区间状态矩阵谱半径的上、下界，即：

$\begin{array}{c}\underset{‾}{\rho }\left(A_{as}\right)=\underset{i=1,2,\mathrm{...}}{\max }\left\{\underset{{\lambda }_{ir}\in {\lambda }_{ir}^I,{\lambda }_{im}\in {\lambda }_{im}^I}{\min }\right\{\sqrt{{\lambda }_{ir}^2+{\lambda }_{im}^2}\left\}\right\}\\ \bar{\rho }\left(A_{as}\right)=\underset{i=1,2,\mathrm{...}}{\max }\left\{\underset{{\lambda }_{ir}\in {\lambda }_{ir}^I,{\lambda }_{im}\in {\lambda }_{im}^I}{\max }\right\{\sqrt{{\lambda }_{ir}^2+{\lambda }_{im}^2}\left\}\right\}\end{array}---\left(26\right)$

当时，不确定气动弹性系统完全稳定；当而ρ(A_as)≤1时，不确定气动弹性系统不完全稳定；当ρ(A_as)>1时，不确定气动弹性系统完全不稳定。使和ρ(A_as)＝1的来流动压分别为不确定气动弹性系统由完全稳定变为不完全稳定的临界动压q和由不完全稳定变为完全不稳定的临界动压其分别对应不确定气动弹性系统颤振速度因子的下界和上界本实例中，在来流马赫数为0.75的条件下，由上述方法获得的颤振速度因子的下界和上界分别为

(8)判断是否完成全马赫数条件下不确定气动弹性系统的颤振速度边界估计，若未完成，调整计算马赫数，重复步骤(2)～(7)，若完成，给出不确定气动弹性系统颤振速度因子上、下界随马赫数的变化情况，由此识别不确定气动弹性系统的完全稳定域、不完全稳定域和完全不稳定域，其中，由颤振速度因子上、下界包裹的区域为气动弹性系统的不完全稳定域，其含义为当飞行速度处于该区域时气动弹性系统可能稳定也可能不稳定，颤振速度因子上界以上的区域为气动弹性系统的完全不稳定域，其含义为当飞行速度处于该区域时气动弹性系统不稳定，颤振速度因子下界以下的区域为气动弹性系统的完全稳定域，其含义为当飞行速度处于该区域时气动弹性系统稳定，据此获得不确定气动弹性系统的颤振速度边界，完成不确定气动弹性系统的稳定性分析；

(9)由步骤(6)中已知来流动压q条件下的气动弹性系统区间状态矩阵，还可直接建立不确定气动弹性系统的鲁棒稳定性快速判据，即：

$\underset{i,j=1}{\overset{n}{\Sigma }}{\mathrm{\Delta a}}_{asij}{\sigma }_{max}\left(P_{ij}\right)<1---\left(27\right)$

(10)判断气动弹性系统区间状态矩阵是否满足步骤(9)中的鲁棒稳定性快速判据条件，若不满足，则不确定气动弹性系统完全不稳定或不完全稳定，若满足，则不确定气动弹 性系统稳定。

综上所述，本发明提出了一种基于气动力不确定降阶的机翼结构气动弹性稳定性分析方法，该方法在构建高效高精度的基于CFD技术的非定常气动力降阶模型的过程中综合考虑存在于数值计算和气动参数的不确定性，将其统一定量化为辨识模型中的不确定但有界区间噪声序列，所建立的不确定性气动力降阶模型，兼顾了气动力模型的精度、计算效率和鲁棒性。本发明中，不确定性气动力降阶模型是以描述多输入多输出系统的动态线化形式的ARMA模型为基础，通过引入反映不确定性因素的区间噪声序列，利用区间集员辨识算法建立的。本发明将不确定性气动力降阶模型与结构模型相耦合，构建了状态空间形式的不确定性气动弹性系统的数学模型，提出了一种从区间状态矩阵特征值角度出发的预测气动弹性系统鲁棒稳定性边界即颤振速度因子上下界的高效方法。另外，本发明基于所建立的气动弹性系统区间状态矩阵，还提供了一种含区间参数的不确定气动弹性系统鲁棒稳定性的快速判据。

以上仅是本发明的具体步骤，对本发明的保护范围不构成任何限制；其可扩展应用于机翼结构气动弹性鲁棒稳定性分析领域，凡采用等同变换或者等效替换而形成的技术方案，均落在本发明权利保护范围之内。