CN102866637B - 一种基于二次降阶的带操纵面机翼非定常气动力模拟方法 - Google Patents

Abstract

本发明提出了一种基于二次降阶的带操纵面机翼非定常气动力模拟方法，首先利用面样条插值方法确定主翼面、操纵面上激振点和拾振点的位置，进而得到相应的频域降阶气动力模型，采用最小状态法进行频域降阶气动力模型从频域到时域的转换，获得初步气动力降阶模型；以操纵面气动力、铰链力矩等效为前提，将操纵面上多处激振点的力等效到两个激振点上，以二次气动力降阶前后对应结构的动态气动弹性响应差异最小为目标函数，采用优化算法进行了二次气动力降阶激振点位置优化。本发明通过对带操纵面机翼的非定常气动力进行二次降阶，在满足试验精度的前提下，减少气动力的阶数，进而减少地面颤振试验中激振器的个数，从而降低激振力控制系统的设计难度。

Description

一种基于二次降阶的带操纵面机翼非定常气动力模拟方法

技术领域

本发明涉及飞行器气动弹性试验技术领域，具体为一种基于二次降阶的带操纵面机翼非定常气动力模拟方法，主要用于地面颤振试验，通过面样条插值降阶和等效降阶方法对带操纵面机翼的非定常气动力模型进行两次降阶，将机翼上分布的非定常气动力通过气动力降阶，等效到有限个加载点处的激振器上。

背景技术

对于每一种新型号的飞机，其颤振速度的获取是必不可少的。地面颤振试验作为一种新兴的颤振研究手段，是将机翼上分布的非定常气动力首先按照气动力降阶的思想，等效为有限个加载点处的非定常气动力，再用激振器在这些加载点处对机翼进行激励，激振器的输出力就是等效后的气动力，而这些集中的非定常气动力是根据机翼上若干个拾振点处的响应（位移、速度、加速度）按照气动力降阶方法计算后得到的。该试验可直接对真实结构进行颤振试验，并通过激振器来模拟给定速度下的气动力，在试验中就可以观察到真实结构的颤振，进而得到颤振速度。

从公开的文献可看到，地面颤振试验最新的研究进展是：在AIAA2011-1942编号文献“GVT-based ground flutter test without wind tunnel”中Zeng等公开了一种对矩形平板和长直机翼的地面颤振试验方法；在SCIENCE CHINA Technological Sciences 2012,55(9):2482-2488的文献“Studies on Aeroservoelasticity Semi-physical Simulation Test forMissiles”中，吴志刚等在考虑了飞行控制系统基础上，针对导弹模型公开了一种气动伺服弹性半物理仿真试验方法。

这些公开的地面颤振试验技术中，主要都是对不带操纵面的机翼气动力降阶建模的研究，而对带操纵面的三元机翼的非定常气动力降阶建模和地面颤振试验尚无相关研究。

发明内容

要解决的技术问题

在地面颤振试验中，通过激振器输出的激振力来模拟分布式的气动力。理论上讲，无论是为了更精确地描述气动面上分布式的气动力，还是多部件试验（如全机试验）的需要，我们都希望使用更多的激振器。但为了保证激振力的准确加载，需进行激振力控制系统的设计，如果激振器数目过多，激振力控制系统的设计难度将大大增加，尤其是在操纵面上，激振器的数目往往有很大的限制。因此在实际颤振试验中，希望在满足试验精度的前提下，尽量减少气动力的阶数（这样，对应的激振器的数目就会相应减少）。

为此，本发明提出了一种基于二次降阶的带操纵面机翼非定常气动力模拟方法，对带操纵面机翼的非定常气动力进行二次降阶，在得到合理的气动力模型的同时，降低激振力控制系统设计的难度。

技术方案

本发明提出的基于二次降阶的带操纵面机翼非定常气动力模拟方法，首先通过面样条插值方法（AIAA 2011-1942编号文献“GVT-based ground futter test without windtunnel”中公开），分别确定主翼面、操纵面上激振点和拾振点的位置，进而得到相应的频域降阶气动力模型，采用最小状态法（Journal of Aircraft,Vol.19,March 1982,pp.221-227中的文献“Design for Active Flutter Suppression and Gust Alleviation UsingState-Space Aeroelastic Modeling”公开）进行频域降阶气动力模型从频域到时域的转换，获得初步气动力降阶模型。为进一步减少操纵面上所需激振器的个数，以操纵面气动力、铰链力矩等效为前提，将操纵面上多处激振点的力等效到两个激振点上。同时，为了得到较优的等效结果，在优选方案中，以二次气动力降阶前后对应结构的动态气动弹性响应差异最小为目标函数，采用优化算法进行了二次气动力降阶激振点位置优化。

本发明的技术方案为：

所述一种基于二次降阶的带操纵面机翼非定常气动力模拟方法，其特征在于：包括以下步骤：

步骤1：建立带操纵面机翼的网格模型，采用面样条插值方法分别确定带操纵面机翼中的主机翼和操纵面上的激振点和拾振点的位置，并得到带操纵面机翼的频域降阶气动力模型；采用最小状态法将频域降阶气动力模型转换到时域空间，得到带操纵面机翼的初步时域气动力降阶模型；所述带操纵面机翼的初步时域气动力降阶模型分为主机翼的初步时域气动力降阶模型和操纵面的初步时域气动力降阶模型；

步骤2：在带操纵面机翼的操纵面上选取两个点作为二次降阶后的激振点，所述两个点中一个在操纵面的铰链轴上，另一个在铰链轴外；

步骤3：将步骤1初步时域气动力降阶模型中操纵面上的N个激振点的激振力等效到步骤2中选定的两个激振点上，并保持拾振点与步骤1中操纵面初步时域气动力降阶模型上的拾振点相同，得到操纵面的二次时域气动力降阶模型；将操纵面的二次时域气动力降阶模型与主机翼的初步时域气动力降阶模型合并得到带操纵面机翼的二次时域气动力降阶模型；所述激振力等效的要求为：

$F_h+F_s=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{F}_i$

$F_s\×d_s=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{F}_i\×d_i$

其中F_h、F_s分别表示步骤2中选定的两个激振点上激振力，F_i表示步骤1初步时域气动力降阶模型中操纵面上的第i个激振点的激振力；d_s表示步骤2中操纵面铰链轴外激振点的激振力相对于铰链轴的力臂值，d_i表示步骤1初步时域气动力降阶模型中操纵面上第i个激振点的激振力相对于铰链轴的力臂值。

所述一种基于二次降阶的带操纵面机翼非定常气动力模拟方法，其特征在于：

步骤4：在模态坐标下建立带操纵面机翼结构的气动弹性运动方程：

$\left[M\right]\left\{\stackrel{\·\·}{\mathrm{\η}}\right\}+\left[B\right]\left\{\stackrel{\·}{\mathrm{\η}}\right\}+\left[K\right]\left\{\mathrm{\η}\right\}=\left\{F\right\}$

其中[M]、[B]、[K]∈R^n×n分别是带操纵面机翼结构的广义质量、广义阻尼和广义刚度矩阵，{n}∈Rⁿ是广义坐标向量，{F}∈Rⁿ是广义气动力向量，n是模态阶数；将带操纵面机翼结构的气动弹性运动方程写为状态空间形式，得到：

$\left\{\begin{array}{c}\stackrel{\·}{\mathrm{\η}}\\ \stackrel{\·\·}{\mathrm{\η}}\end{array}\right\}=\left[\begin{array}{cc}0& O\\ -M^{-1}K& -M^{-1}B\end{array}\right]\left\{\begin{array}{c}\mathrm{\η}\\ \stackrel{\·}{\mathrm{\η}}\end{array}\right\}\left\{\begin{array}{c}0\\ M^{-1}F\end{array}\right\}$

步骤5：将步骤1中带操纵面机翼的初步时域气动力降阶模型代入步骤4中带操纵面机翼结构的气动弹性运动方程的状态空间形式中，建立带操纵面机翼的一次降阶颤振仿真模型；将步骤3中带操纵面机翼的二次时域气动力降阶模型代入步骤4中带操纵面机翼结构的气动弹性运动方程的状态空间形式中，建立带操纵面机翼的二次降阶颤振仿真模型；在带操纵面机翼的一次降阶颤振仿真模型和二次降阶颤振仿真模型中给带操纵面机翼结构施加相同的飞行速度参数和相同的脉冲激励，得到带操纵面机翼结构二次降阶前的时域气弹响应和二次降阶后的时域气弹响应；

步骤6：重复步骤2至步骤6，以

${\mathrm{obj}}_2=\underset{t=0}{\overset{t_n}{\mathrm{\Σ}}}{(r_e\left(t\right)-r\left(t\right))}^2$

取最小值为目标函数，对操纵面上两个二次降阶后的激振点位置采用优化算法寻优，其中t_n是步骤5中颤振仿真模型总的仿真时长，r(t)、r_e(t)分别为t时刻，带操纵面机翼结构上一个参考点的二次降阶前的时域气弹响应值和二次降阶后的时域气弹响应值；以寻优结果建立操纵面的二次时域气动力降阶优化模型，将操纵面的二次时域气动力降阶优化模型与主机翼的初步时域气动力降阶模型合并得到带操纵面机翼的二次时域气动力降阶优化模型。

有益效果

本发明通过对带操纵面机翼的非定常气动力进行二次降阶，在满足试验精度的前提下，减少气动力的阶数，进而减少地面颤振试验中激振器的个数，从而降低激振力控制系统的设计难度。

附图说明

图1：地面颤振试验系统原理图；

图2：带操纵面的机翼模型；（图中Thickness=8.5表示模型厚度为8.5mm）

图3：带操纵面机翼气动网格图；

图4：8激振点/8拾振点布置图；

图5：采用Matlab软件建立的带操纵面机翼颤振仿真模型框图；

图6：6激振点/8拾振点布置图；

图7：二次降阶前后，带操纵面机翼的时域气弹响应对比；

具体实施方式

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

参照附图1，飞行器机翼结构的地面颤振试验原理为：试验中先预设一个飞行条件，由传感器和信号采集系统获取结构的运动信息，通过非定常气动力实时生成软件计算出相应的分布气动载荷，经过降阶后再由激振器系统施加到机翼结构上，如此反复，从而达到“结构—气动力(由激振器模拟)—结构”相互作用的动态气动弹性耦合现象。在不同飞行条件下，只需对机翼结构施加一个初始激励，根据机翼结构的振动响应是否收敛，即可判断出颤振现象是否发生，并得到颤振速度。激振力控制系统是保证各激振器实际输出力与气动力生成系统的输出力信号误差最小。

而本实施例就是通过对带操纵面机翼的非定常气动力进行二次降阶，在满足试验精度的前提下，减少气动力的阶数，进而减少地面颤振试验中激振器的个数，从而降低激振力控制系统的设计难度。

本实例中的气弹分析均是基于Nastran的SOL 145进行的。对图2所示的带操纵面后掠翼进行气动力降阶处理，采用根部固支的约束方式，材料为铝，E＝70Gpa，v＝0.3，ρ＝2750kg/m³。对机翼划分结构网格如图4所示。

步骤1：建立带操纵面机翼的网格模型，设定气弹分析的条件为：空气密度ρ_空=1.226kg/m³，马赫数Ma=0.8，减缩频率k=0.001、0.01、0.05、0.115、0.3、0.6、1.0，带操纵面机翼半弦长b=0.46m。划分的气动力网格如图3所示。

分别在主翼面和操纵面上各布置初始的4个激振点和4个拾振点，通过面样条插值方法以及遗传算法优化分别确定带操纵面机翼中的8个激振点和8个拾振点的位置（如图4所示，▲表示激振点，●表示拾振点，激振点和拾振点重合时用■表示），并得到带操纵面机翼的频域降阶气动力模型。

所述的面样条插值方法以及遗传算法优化过程采用的是AIAA 2011-1942编号文献“GVT-based ground flutter test without wind tunnel”中公开的方法：

地面颤振试验中，假设分析对象是一个机翼，通过下式来表示机翼上各结点位移与气动力之间的关系：

{f_a(ω)}＝q_∞[A(ω)]{z}

其中，q_∞是动压，ω是简谐振动圆频率，A(ω)是对应频域内的气动力影响系数（Aerodynamics Influence Coefficient,AIC）阵，z、f_a分别是各气动网格点处的垂向位移、气动力。

对式{f_a(ω)}＝q_∞[A(ω)]{z}采用面样条插值的初步降阶过程为：

1、建立从传感器位置处到气动网格下洗点处的位移变换关系：

{z}_NP×1＝[G_S]_NP×NS{z_sen}_NS×1

其中[G_s]是从传感器位置处到气动网格下洗点处的位置变换矩阵，{z_sen}_NS×1表示NS个传感器位置处的位移，{z}_NP×1表示NP个气动网格下洗点处的位移；

2、建立从气动网格上升力作用点处到激振器位置处的气动力变换关系：

{f_shak}_NA×1＝[G_f]_NA×NP{f_a}_NP×1

其中[G_f]是从气动网格上升力作用点处到激振器位置处的气动力变换矩阵，{f_shak}_NA×1表示NA个激振器位置处的气动力，{f_a}_NP×1表示NP个气动网格升力作用点处的气动力；

3、将{z}_NP×1＝[G_S]_NP×NS{z_sen}_NS×1和{f_shak}_NA×1＝[G_f]_NA×NP{f_a}_NP×1代入{f_a(ω)}＝q_∞[A(ω)]{z}中，从而直接对频域内的AIC阵进行降阶：

$\left\{f_{\mathrm{shak}}\right\}=q_{\∞}\left[\stackrel{\‾}{A}\left(\mathrm{\ω}\right)\right]\left\{z_{\mathrm{sen}}\right\}$

其中降阶后的AIC阵为 $\left[\stackrel{\‾}{A}\left(\mathrm{\ω}\right)\right]=\left[G_f\right]\left[A\left(\mathrm{\ω}\right)\right]\left[G_S\right],$ {f_shak}即为初步的带操纵面机翼频域一次降阶气动力模型。

根据降阶前后的AIC阵，可分别建立机翼的颤振方程，使用v-g法或p-k法等方法可得到频域颤振结果。由于选择插值点的位置不同时（即不同的激振点/拾振点位置），降阶前后求得的频域颤振结果相差较大，还需要进行激振点/拾振点的位置优化，本实施例中采用遗传算法优化，优化的目标函数为：

${\mathrm{obj}}_1={\left(\frac{V_f-V_0}{V_0}\right)}^2+{\left(\frac{{\mathrm{\ω}}_f-{\mathrm{\ω}}_0}{{\mathrm{\ω}}_0}\right)}^2$

取最小值，其中V₀、ω₀分别是降阶前计算得到的颤振速度、颤振频率，V_f、ω_f分别是降阶后计算得到的颤振速度、颤振频率。通过优化过程，可得到最优的激振点/拾振点位置，进而得到带操纵面机翼的频域降阶气动力模型。

本实施例中，不同减缩频率下，带操纵面机翼的频域降阶AIC阵为：

k=0.001：

[2.5598-0.0007i,-2.6315-0.0005i,0.2970-0.0019i,-0.2253+0.0004i,-0.0833-0.0004i,0.0346+0.0002i,-0.0744-0.0002i,0.1232+0.0004i;

-0.9107+0.0022i,0.8340-0.0022i,0.1257+0.0011i,-0.0490-0.0004i,0.1236+0.0001i,-0.0530-0.0000i,0.1058+0.0000i,-0.1764-0.0001i;

0.4747+0.0008i,-0.8124-0.0012i,0.7892-0.0004i,-0.4515-0.0006i,0.3822-0.0004i,-0.1864+0.0002i,0.2672-0.0003i,-0.4629+0.0003i;

-0.3200+0.0004i,0.2678-0.0004i,0.1014+0.0004i,-0.0491-0.0003i,0.3020-0.0000i,-0.1737-0.0000i,0.1408-0.0000i,-0.2691-0.0000i;

-0.0363+0.0000i,0.0793-0.0000i,-0.0993+0.0001i,0.0563+0.0000i,0.5249+0.0001i,-0.3431-0.0002i,0.1347+0.0001i,-0.3164-0.0002i;

0.0018+0.0000i,-0.0020-0.0000i,0.0006-0.0000i,-0.0004+0.0000i,0.0068+0.0000i,-0.0177-0.0000i,-0.0335-0.0000i,0.0443+0.0000i;

-0.0396+0.0000i,0.0677-0.0001i,-0.0657+0.0002i,0.0376-0.0000i,0.3452+0.0001i,-0.1535-0.0000i,0.2811+0.0000i,-0.4729-0.0002i;

0.0030-0.0000i,-0.0090+0.0000i,0.0139-0.0000i,-0.0078-0.0000i,-0.0804-0.0000i,0.0664+0.0000i,0.0163-0.0000i,-0.0024-0.0000i]

k=0.01：

[2.5599-0.0071i,-2.6316-0.0049i,0.2969-0.0190i,-0.2253+0.0038i,-0.0835-0.0036i,0.0347+0.0019i,-0.0745-0.0021i,0.1233+0.0042i;

-0.9107+0.0224i,0.8341-0.0223i,0.1258+0.0109i,-0.0490-0.0040i,0.1236+0.0008i,-0.0530-0.0004i,0.1059+0.0003i,-0.1765-0.0012i;

0.4747+0.0083i,-0.8124-0.0123i,0.7892-0.0037i,-0.4515-0.0058i,0.3821-0.0039i,-0.1864+0.0018i,0.2672-0.0026i,-0.4629+0.0032i;

-0.3200+0.0042i,0.2678-0.0038i,0.1014+0.0043i,-0.0491-0.0029i,0.3020-0.0001i,-0.1737-0.0003i,0.1408-0.0002i,-0.2692-0.0005i;

-0.0363+0.0000i,0.0793-0.0003i,-0.0993+0.0013i,0.0563+0.0005i,0.5249+0.0009i,-0.3431-0.0015i,0.1347+0.0005i,-0.3164-0.0015i;

0.0018+0.0000i,-0.0020-0.0000i,0.0006-0.0000i,-0.0004+0.0000i,0.0068+0.0000i,-0.0177-0.0002i,-0.0335-0.0001i,0.0443+0.0003i;

-0.0396+0.0002i,0.0677-0.0005i,-0.0656+0.0015i,0.0376-0.0001i,0.3452+0.0008i,-0.1535-0.0005i,0.2811+0.0003i,-0.4729-0.0020i;

0.0030-0.0000i,-0.0090+0.0001i,0.0139-0.0001i,-0.0078-0.0001i,-0.0804-0.0002i,0.0664+0.0004i,0.0163-0.0000i,-0.0024-0.0001i]

k=0.05：

[2.5620-0.0356i,-2.6329-0.0248i,0.2941-0.0948i,-0.2250+0.0189i,-0.0869-0.0175i,0.0364+0.0090i,-0.0765-0.0104i,0.1266+0.0207i;

-0.9109+0.1122i,0.8360-0.1114i,0.1283+0.0545i,-0.0492-0.0200i,0.1249+0.0037i,-0.0537-0.0020i,0.1066+0.0016i,-0.1778-0.0060i;

0.4736+0.0419i,-0.8115-0.0618i,0.7901-0.0184i,-0.4510-0.0291i,0.3808-0.0193i,-0.1858+0.0090i,0.2664-0.0130i,-0.4618+0.0158i;

-0.3206+0.0211i,0.2682-0.0192i,0.1022+0.0216i,-0.0488-0.0143i,0.3024-0.0007i,-0.1739-0.0016i,0.1410-0.0008i,-0.2695-0.0022i;

-0.0363+0.0002i,0.0794-0.0017i,-0.0992+0.0066i,0.0564+0.0023i,0.5249+0.0046i,-0.3431-0.0076i,0.1347+0.0025i,-0.3164-0.0077i;

0.0018+0.0002i,-0.0020-0.0002i,0.0006-0.0002i,-0.0004+0.0002i,0.0068+0.0002i,-0.0176-0.0011i,-0.0335-0.0005i,0.0443+0.0015i;

-0.0397+0.0010i,0.0678-0.0025i,-0.0655+0.0075i,0.0377-0.0003i,0.3452+0.0040i,-0.1535-0.0025i,0.2812+0.0013i,-0.4729-0.0100i;

0.0030-0.0001i,-0.0090+0.0003i,0.0139-0.0006i,-0.0078-0.0006i,-0.0804-0.0008i,0.0664+0.0020i,0.0163-0.0001i,-0.0024-0.0003i]

k=0.115：

[2.5724-0.0814i,-2.6395-0.0575i,0.2822-0.2160i,-0.2238+0.0430i,-0.1001-0.0353i,0.0432+0.0182i,-0.0844-0.0213i,0.1397+0.0432i;

-0.9122+0.2582i,0.8447-0.2566i,0.1390+0.1242i,-0.0503-0.0459i,0.1297+0.0064i,-0.0562-0.0035i,0.1097+0.0025i,-0.1829-0.0117i;

0.4697+0.0982i,-0.8087-0.1436i,0.7944-0.0413i,-0.4490-0.0674i,0.3753-0.0431i,-0.1832+0.0200i,0.2632-0.0292i,-0.4572+0.0352i;

-0.3232+0.0489i,0.2700-0.0446i,0.1060+0.0496i,-0.0474-0.0330i,0.3038-0.0021i,-0.1747-0.0034i,0.1419-0.0022i,-0.2710-0.0046i;

-0.0365+0.0006i,0.0796-0.0040i,-0.0989+0.0153i,0.0569+0.0052i,0.5252+0.0105i,-0.3431-0.0174i,0.1348+0.0058i,-0.3164-0.0177i;

0.0018+0.0003i,-0.0020-0.0005i,0.0005-0.0005i,-0.0004+0.0004i,0.0068+0.0006i,-0.0176-0.0025i,-0.0335-0.0011i,0.0443+0.0035i;

-0.0402+0.0023i,0.0684-0.0059i,-0.0648+0.0173i,0.0381-0.0008i,0.3454+0.0092i,-0.1536-0.0057i,0.2814+0.0030i,-0.4728-0.0231i;

0.0030-0.0002i,-0.0090+0.0007i,0.0139-0.0014i,-0.0079-0.0013i,-0.0804-0.0017i,0.0663+0.0046i,0.0163-0.0002i,-0.0023-0.0006i]

k=0.3：

[2.6547-0.2201i,-2.6930-0.1409i,0.2072-0.5405i,-0.2182+0.1098i,-0.1529-0.0395i,0.0702+0.0202i,-0.1190-0.0255i,0.1970+0.0636i;

-0.9193+0.6822i,0.9051-0.6786i,0.2111+0.3057i,-0.0577-0.1174i,0.1436-0.0061i,-0.0633+0.0028i,0.1202-0.0075i,-0.2014-0.0074i;

0.4560+0.2718i,-0.8010-0.3865i,0.8369-0.0996i,-0.4399-0.1805i,0.3450-0.0994i,-0.1690+0.0460i,0.2456-0.0688i,-0.4331+0.0810i;

-0.3392+0.1331i,0.2805-0.1218i,0.1340+0.1282i,-0.0380-0.0870i,0.3088-0.0119i,-0.1776-0.0055i,0.1452-0.0093i,-0.2777-0.0059i;

-0.0379+0.0020i,0.0812-0.0112i,-0.0963+0.0405i,0.0608+0.0133i,0.5271+0.0268i,-0.3426-0.0452i,0.1354+0.0151i,-0.3167-0.0459i;

0.0020+0.0008i,-0.0022-0.0012i,0.0002-0.0014i,-0.0002+0.0010i,0.0067+0.0015i,-0.0173-0.0064i,-0.0334-0.0028i,0.0440+0.0091i;

-0.0432+0.0070i,0.0723-0.0166i,-0.0595+0.0459i,0.0411-0.0025i,0.3467+0.0231i,-0.1542-0.0145i,0.2827+0.0074i,-0.4726-0.0596i;

0.0029-0.0005i,-0.0089+0.0018i,0.0140-0.0036i,-0.0085-0.0034i,-0.0806-0.0046i,0.0660+0.0120i,0.0163-0.0004i,-0.0021-0.0016i]

k=0.6：

[2.9104-0.5318i,-2.8488-0.2037i,-0.0037-0.9538i,-0.1961+0.2090i,-0.1631+0.0195i,0.0741-0.0096i,-0.1371+0.0187i,0.2308+0.0128i;

-0.8965+1.4059i,1.0812-1.4019i,0.4091+0.5070i,-0.0793-0.2161i,0.1204-0.0308i,-0.0508+0.0146i,0.1096-0.0336i,-0.1875+0.0182i;

0.4625+0.5789i,-0.8120-0.7950i,1.0175-0.2155i,-0.4252-0.3671i,0.2710-0.1516i,-0.1340+0.0698i,0.2010-0.1131i,-0.3738+0.1287i;

-0.3746+0.2921i,0.3001-0.2677i,0.2345+0.2362i,-0.0091-0.1768i,0.3044-0.0351i,-0.1760-0.0049i,0.1442-0.0268i,-0.2803+0.0023i;

-0.0416+0.0070i,0.0855-0.0269i,-0.0843+0.0843i,0.0745+0.0249i,0.5333+0.0503i,-0.3408-0.0886i,0.1377+0.0287i,-0.3176-0.0900i;

0.0025+0.0014i,-0.0029-0.0021i,-0.0013-0.0028i,0.0003+0.0023i,0.0064+0.0035i,-0.0161-0.0131i,-0.0333-0.0054i,0.0429+0.0178i;

-0.0512+0.0194i,0.0828-0.0400i,-0.0363+0.0932i,0.0506-0.0080i,0.3482+0.0409i,-0.1552-0.0260i,0.2864+0.0121i,-0.4705-0.1144i;

0.0027-0.0010i,-0.0086+0.0038i,0.0143-0.0078i,-0.0108-0.0068i,-0.0815-0.0091i,0.0648+0.0239i,0.0160-0.0009i,-0.0010-0.0029i]

k=1.0：

[3.2883-1.1021i,-3.0734-0.1312i,-0.2752-1.2790i,-0.1690+0.2971i,-0.1211+0.0315i,0.0525-0.0131i,-0.1087+0.0451i,0.1930-0.0107i;

-0.6956+2.3304i,1.3773-2.4146i,0.7029+0.6300i,-0.0928-0.3134i,0.0938-0.0209i,-0.0379+0.0083i,0.0882-0.0390i,-0.1561+0.0099i;

0.5675+0.9412i,-0.9010-1.3123i,1.3926-0.4896i,-0.4009-0.5897i,0.1820-0.1523i,-0.0924+0.0683i,0.1363-0.1289i,-0.2887+0.1377i;

-0.3908+0.5203i,0.2897-0.4939i,0.4479+0.3113i,0.0551-0.2980i,0.2866-0.0577i,-0.1687-0.0074i,0.1338-0.0466i,-0.2735+0.0105i;

-0.0429+0.0204i,0.0856-0.0601i,-0.0433+0.1442i,0.1060+0.0335i,0.5423+0.0720i,-0.3348-0.1404i,0.1411+0.0421i,-0.3180-0.1422i;

0.0030+0.0015i,-0.0034-0.0025i,-0.0058-0.0036i,0.0020+0.0048i,0.0067+0.0067i,-0.0140-0.0223i,-0.0324-0.0082i,0.0398+0.0282i;

-0.0569+0.0455i,0.0930-0.0849i,0.0291+0.1438i,0.0675-0.0225i,0.3408+0.0570i,-0.1517-0.0368i,0.2893+0.0129i,-0.4577-0.1771i;

0.0022-0.0018i,-0.0076+0.0074i,0.0140-0.0151i,-0.0168-0.0111i,-0.0839-0.0144i,0.0625+0.0393i,0.0150-0.0014i,0.0021-0.0044i]

由于在地面颤振试验中，为了实现降阶气动力的实时加载，需要将频域降阶气动力模型转换到时域空间，这里采用了Journal of Aircraft,Vol.19,March 1982,pp.221-227.中的文献“Design for Active Flutter Suppression and Gust Alleviation Using State-SpaceAeroelastic Modeling”公开的最小状态法进行降阶气动力模型从频域到时域的转换，得到带操纵面机翼的初步时域气动力降阶模型。

所述最小状态法公式如下式：

$\left[\stackrel{\‾}{A}\left(t\right)\right]=\left[{\stackrel{\‾}{A}}_0\right]+\frac{b}{V}\left[{\stackrel{\‾}{A}}_1\right]s+\frac{b^2}{V^2}\left[{\stackrel{\‾}{A}}_2\right]s^2+\left[\stackrel{\‾}{D}\right]{\left[s\right[I]-\frac{V}{b}[\stackrel{\‾}{R}\left]\right]}^{-1}\left[\stackrel{\‾}{E}\right]s$

其中s是拉普拉斯变量，b是半弦长，V是来流速度，是对进行最小状态法处理得到的系数矩阵，是时域的降阶AIC阵。用该式取代之前式中的得到用于实时气动力生成的气动力计算公式：

$\left\{f_{\mathrm{shak}}\right\}=q_{\∞}\left[\right[{\stackrel{\‾}{A}}_0]+\frac{b}{V}[{\stackrel{\‾}{A}}_1]s+\frac{b^2}{V^2}[{\stackrel{\‾}{A}}_2\left]s^2\right]\left\{z_{\mathrm{sen}}\right\}+q_{\∞}\left[\stackrel{\‾}{D}\right]{\left[s\right[I]-\frac{V}{b}[\stackrel{\‾}{R}\left]\right]}^{-1}\left[\stackrel{\‾}{E}\right]\left\{{\stackrel{\·}{z}}_{\mathrm{sen}}\right\}.$

本实施例中，通过最小状态法，得到实时气动力生成的气动力计算公式对应的矩阵为：

$\left[{\stackrel{\‾}{A}}_0\right]:$

[2.5598,-2.6315,0.2970,-0.2253,-0.0833,0.0346,-0.0744,0.1232;-0.9107,0.8340,0.1257,-0.0490,0.1236,-0.0530,0.1058,-0.1764;0.4747,-0.8124,0.7892,-0.4515,0.3822,-0.1864,0.2672,-0.4629;-0.3200,0.2678,0.1014,-0.0491,0.3020,-0.1737,0.1408,-0.2691;-0.0363,0.0793,-0.0993,0.0563,0.5249,-0.3431,0.1347,-0.3164;0.0018,-0.0020,0.0006,-0.0004,0.0068,-0.0177,-0.0335,0.0443;-0.0396,0.0677,-0.0657,0.0376,0.3452,-0.1535,0.2811,-0.4729;0.0030,-0.0090,0.0139,-0.0078,-0.0804,0.0664,0.0163,-0.0024]

$\left[{\stackrel{\‾}{A}}_1\right]:$

[-0.7873,-0.4342,-1.7428,0.3656,-0.0199,0.0171,-0.0188,0.1191;2.2786,-2.3065,0.9677,-0.3970,-0.1036,0.0450,-0.0733,0.0402;0.9194,-1.2943,-0.3613,-0.6054,-0.3521,0.1630,-0.2416,0.2866；0.4558,-0.4135,0.4060,-0.2903,-0.0670,-0.0059,-0.0475,0.0037;0.0078,-0.0378,0.1316,0.0444,0.0870,-0.1496,0.0486,-0.1506;0.0027,-0.0040,-0.0042,0.0035,0.0054,-0.0216,-0.0092,0.0301;0.0257,-0.0573,0.1471,-0.0082,0.0724,-0.0460,0.0218,-0.1942;-0.0016,0.0060,-0.0121,-0.0112,-0.0152,0.0399,-0.0014,-0.0053]

$\left[{\stackrel{\‾}{A}}_2\right]:$

[-0.9407,0.5470,0.8180,-0.1051,0.2365,-0.1528,0.1806,-0.3813;0.0568,-0.6672,-0.8229,0.1501,0.0762,-0.0069,0.0187,0.0519;0.1025,-0.0675,-0.5866,-0.0953,0.3562,-0.1649,0.2106,-0.2724;0.1631,-0.0945,-0.3473,-0.1360,0.0085,0.0060,0.0030,0.0284;0.0128,-0.0177,-0.0318,-0.0604,-0.0236,-0.0111,-0.0059,-0.0029;-0.0021,0.0026,0.0044,-0.0023,0.0014,-0.0053,-0.0010,0.0054;0.0306,-0.0436,-0.0652,-0.0433,-0.0066,0.0046,-0.0134,-0.0155;0.0011,-0.0013,-0.0014,0.0100,0.0030,0.0060,0.0009,-0.0042]

$\left[\stackrel{\‾}{D}\right]:$

[-7.0920,-1.1837,-1.4817;

6.8868,2.7627,1.8032;

3.1351,0.4691,2.6583;

2.1733,0.3942,0.9388;

0.2288,0.0493,0.1109;

-0.0265,-0.0062,-0.0131;

0.4699,0.1124,0.2265;

0.0063,0.0025,0.0031]

$\left[\stackrel{\‾}{E}\right]:$

[-0.0043,0.0019,0.0060,-0.0008,0.0149,-0.0073,0.0086,-0.0132;

0.0094,0.0020,-0.0015,0.0006,-0.0118,0.0062,-0.0067,0.0109;

-0.0017,0.0008,-0.0049,0.0023,-0.0155,0.0075,-0.0089,0.0137]

$\left[\stackrel{\‾}{R}\right]:$

[-0.2000,0,0;

0,-0.1790,0;

0,0,-0.1580]

自此，通过步骤1得到了带操纵面机翼的初步时域气动力降阶模型，该模型分为主机翼的初步时域气动力降阶模型和操纵面的初步时域气动力降阶模型。下面主要对操纵面进行二次降阶。

步骤2：在带操纵面机翼的操纵面上选取两个点作为二次降阶后的激振点，所述两个点中一个在操纵面的铰链轴上，另一个在铰链轴外。其目的就是要通过两个点来替代操纵面上第一次降阶后得到的多个激振点。

步骤3：将步骤1操纵面的初步时域气动力降阶模型中操纵面上的N个激振点的激振力等效到步骤2中选定的两个激振点上，并保持拾振点与步骤1中操纵面初步时域气动力降阶模型上的拾振点相同，得到操纵面的二次时域气动力降阶模型；将操纵面的二次时域气动力降阶模型与主机翼的初步时域气动力降阶模型合并得到带操纵面机翼的二次时域气动力降阶模型；所述激振力等效的要求为：

$F_h+F_s=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{F}_i$

$F_s\×d_s=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{F}_i\×d_i$

其中F_h、F_s分别表示步骤2中选定的两个激振点上激振力，F_i表示步骤1初步时域气动力降阶模型中操纵面上的第i个激振点的激振力；d_s表示步骤2中操纵面铰链轴外激振点的激振力相对于铰链轴的力臂值，d_i表示步骤1初步时域气动力降阶模型中操纵面上第i个激振点的激振力相对于铰链轴的力臂值。

此时，由于步骤2中的两个激振点位置选取不同，会导致带操纵面机翼的时域气弹响应相差比较大，所以，优选方案中需要对步骤2中的两个激振点位置进行优化。

步骤4：在模态坐标下建立带操纵面机翼结构的气动弹性运动方程：

$\left[M\right]\left\{\stackrel{\·\·}{\mathrm{\η}}\right\}+\left[B\right]\left\{\stackrel{\·}{\mathrm{\η}}\right\}+\left[K\right]\left\{\mathrm{\η}\right\}=\left\{F\right\}$

其中[M]、[B]、[K]∈R^n×n分别是带操纵面机翼结构的广义质量、广义阻尼和广义刚度矩阵，{n}∈Rⁿ是广义坐标向量，{F}∈Rⁿ是广义气动力向量，n是模态阶数；将带操纵面机翼结构的气动弹性运动方程写为状态空间形式，得到：

$\left\{\begin{array}{c}\stackrel{\·}{\mathrm{\η}}\\ \stackrel{\·\·}{\mathrm{\η}}\end{array}\right\}=\left[\begin{array}{cc}0& O\\ -M^{-1}K& -M^{-1}B\end{array}\right]\left\{\begin{array}{c}\mathrm{\η}\\ \stackrel{\·}{\mathrm{\η}}\end{array}\right\}\left\{\begin{array}{c}0\\ M^{-1}F\end{array}\right\}.$

步骤5：将步骤1中带操纵面机翼的初步时域气动力降阶模型代入步骤4中带操纵面机翼结构的气动弹性运动方程的状态空间形式中，建立带操纵面机翼的一次降阶颤振仿真模型；将步骤3中带操纵面机翼的二次时域气动力降阶模型代入步骤4中带操纵面机翼结构的气动弹性运动方程的状态空间形式中，建立带操纵面机翼的二次降阶颤振仿真模型；在带操纵面机翼的一次降阶颤振仿真模型和二次降阶颤振仿真模型中给带操纵面机翼结构施加相同的飞行速度参数和相同的脉冲激励，得到带操纵面机翼结构二次降阶前的时域气弹响应和二次降阶后的时域气弹响应；

需要注意的是，步骤4中的操纵面机翼结构的气动弹性运动方程状态空间形式是在模态坐标下，而带操纵面机翼的初步时域气动力降阶模型和二次时域气动力降阶模型是在物理坐标下，在代入时需要分别对拾振点处的位移和激振点的输出力进行模态坐标变换：

{z_sen}＝[φ_z]{n}

{F}＝[φ_f]{f_shak}

其中[φ_z]是模态位移到物理位移的模态坐标变换矩阵，[φ_f]是物理气动力到模态气动力的模态坐标变换矩阵。

建立的带操纵面机翼的一次降阶颤振仿真模型和二次降阶颤振仿真模型框图如图5所示，该仿真模型是采用Matlab软件建立，其中disturbation模块表示在t=0时刻，对结构施加的一个脉冲激励。

步骤6：重复步骤2至步骤6，以

${\mathrm{obj}}_2=\underset{t=0}{\overset{t_n}{\mathrm{\Σ}}}{(r_e\left(t\right)-r\left(t\right))}^2$

取最小值为目标函数，对操纵面上两个二次降阶后的激振点位置采用遗传算法寻优，其中t_n是步骤5中颤振仿真模型总的仿真时长，r(t)、r_e(t)分别为t时刻，带操纵面机翼结构上一个参考点的二次降阶前的时域气弹响应值和二次降阶后的时域气弹响应值，本实施例中参考点为图4中的传感器M点，通过优化得到的最佳激振器位置如图6所示。以寻优结果建立操纵面的二次时域气动力降阶优化模型，将操纵面的二次时域气动力降阶优化模型与主机翼的初步时域气动力降阶模型合并得到带操纵面机翼的二次时域气动力降阶优化模型。

通过时域仿真，对降阶前后的带操纵面机翼结构时域气弹响应的比较，如图4中R点的响应在图7所示，由图7可知，以机翼结构上某一点的气动弹性响应差异最小作为二次降阶的目标函数，通过优化方法可以确定操纵面上2个激振点的最优布置方案。

各次降阶前后的颤振速度比较如表1所示：

表1带操纵面机翼的颤振特性

由表1的结果知，采用此方案，得到的降阶前后的颤振速度误差很小，满足地面颤振试验的要求。因此，采用本方法，可以对带操纵面机翼的地面颤振试验进行气动力降阶处理，在保证试验精度的前提下，有效减少激振点数目，得到正确的颤振模拟试验结果。

Claims (2)

1.一种基于二次降阶的带操纵面机翼非定常气动力模拟方法，其特征在于：包括以下步骤：

步骤1：建立带操纵面机翼的网格模型，采用面样条插值方法分别确定带操纵面机翼中的主机翼和操纵面上的激振点和拾振点的位置，并得到带操纵面机翼的频域降阶气动力模型；采用最小状态法将频域降阶气动力模型转换到时域空间，得到带操纵面机翼的初步时域气动力降阶模型；所述带操纵面机翼的初步时域气动力降阶模型分为主机翼的初步时域气动力降阶模型和操纵面的初步时域气动力降阶模型；

步骤2：在带操纵面机翼的操纵面上选取两个点作为二次降阶后的激振点，所述两个点中一个在操纵面的铰链轴上，另一个在铰链轴外；

步骤3：将步骤1初步时域气动力降阶模型中操纵面上的N个激振点的激振力等效到步骤2中选定的两个激振点上，并保持拾振点与步骤1中操纵面初步时域气动力降阶模型上的拾振点相同，得到操纵面的二次时域气动力降阶模型；将操纵面的二次时域气动力降阶模型与主机翼的初步时域气动力降阶模型合并得到带操纵面机翼的二次时域气动力降阶模型；所述激振力等效的要求为：

$\begin{array}{c}{F}_h+F_S=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{F}_i\\ F_s\×d_s=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{F}_i\×d_i\end{array}$

其中F_h、F_s分别表示步骤2中选定的两个激振点上激振力，F_i表示步骤1初步时域气动力降阶模型中操纵面上的第i个激振点的激振力；d_s表示步骤2中操纵面铰链轴外激振点的激振力相对于铰链轴的力臂值，d_i表示步骤1初步时域气动力降阶模型中操纵面上第i个激振点的激振力相对于铰链轴的力臂值。

2.根据权利要求1所述的一种基于二次降阶的带操纵面机翼非定常气动力模拟方法，其特征在于：

步骤4：在模态坐标下建立带操纵面机翼结构的气动弹性运动方程：

$\left[M\right]\left\{\stackrel{\·\·}{\mathrm{\η}}\right\}+\left[B\right]\left\{\stackrel{\·}{\mathrm{\η}}\right\}+\left[K\right]\left\{\mathrm{\η}\right\}=\left\{F\right\}$

其中[M]、[B]、[K]∈R^n×n分别是带操纵面机翼结构的广义质量、广义阻尼和广义刚度矩阵，{η}∈Rⁿ是广义坐标向量，{F}∈Rⁿ是广义气动力向量，n是模态阶数；将带操纵面机翼结构的气动弹性运动方程写为状态空间形式，得到：

$\left\{\begin{array}{c}\stackrel{\·}{\mathrm{\η}}\\ \stackrel{\·\·}{\mathrm{\η}}\end{array}\right\}=\left[\begin{array}{cc}0& I\\ {-M}^{-1}K& {-M}^{-1}B\end{array}\right]\left\{\begin{array}{c}\mathrm{\η}\\ \stackrel{\·}{\mathrm{\η}}\end{array}\right\}+\left\{\begin{array}{c}0\\ M^{-1}F\end{array}\right\}$

步骤5：将步骤1中带操纵面机翼的初步时域气动力降阶模型代入步骤4中带操纵面机翼结构的气动弹性运动方程的状态空间形式中，建立带操纵面机翼的一次降阶颤振仿真模型；将步骤3中带操纵面机翼的二次时域气动力降阶模型代入步骤4中带操纵面机翼结构的气动弹性运动方程的状态空间形式中，建立带操纵面机翼的二次降阶颤振仿真模型；在带操纵面机翼的一次降阶颤振仿真模型和二次降阶颤振仿真模型中给带操纵面机翼结构施加相同的飞行速度参数和相同的脉冲激励，得到带操纵面机翼结构二次降阶前的时域气弹响应和二次降阶后的时域气弹响应；

步骤6：重复步骤2至步骤5，以

${\mathrm{obj}}_2=\underset{t=0}{\overset{t_n}{\mathrm{\Σ}}}{(r_e\left(t\right)-r\left(t\right))}^2$

取最小值为目标函数，对操纵面上两个二次降阶后的激振点位置采用优化算法寻优，其中t_n是步骤5中颤振仿真模型总的仿真时长，r(t)、r_e(t)分别为t时刻，带操纵面机翼结构上一个参考点的二次降阶前的时域气弹响应值和二次降阶后的时域气弹响应值；以寻优结果建立操纵面的二次时域气动力降阶优化模型，将操纵面的二次时域气动力降阶优化模型与主机翼的初步时域气动力降阶模型合并得到带操纵面机翼的二次时域气动力降阶优化模型。