CN113128085B - 一种基于状态空间形式涡格法的阵风减缓分析方法 - Google Patents

Abstract

本发明公开了一种基于状态空间形式涡格法的阵风减缓分析方法，依据给定模型建立气动网格，根据状态空间形式涡格法建立模型附着涡涡强、尾涡涡强与机翼表面诱导速度的关系，插值得到有限元模态及舵面模态广义坐标与机翼表面诱导速度的关系后形成气动弹性耦合方程，气动弹性耦合方程同样为状态空间形式，给定阵风模型作为外激励，可进行气动弹性阵风分析；在气动弹性耦合方程的基础上以机翼某点加速度作为系统输出，经过控制器后形成舵面偏转输入信号，搭建闭环控制系统，进行气动弹性阵风减缓仿真分析。本发明的方法能够提高阵风分析精度，避免频域时域转换，提高分析效率。

Description

一种基于状态空间形式涡格法的阵风减缓分析方法

技术领域

本发明属于结构动力学和气动弹性力学分析技术领域，尤其涉及一种基于状态空间形式涡格法的阵风减缓分析方法。

背景技术

气动弹性力学作为应用力学的分支，主要研究气动力、弹性力与惯性力之间耦合的问题。在气动力作用下弹性结构发生振动与变形，结构弹性运动反过来又会导致气动力的大小及分布发生变化，这种相互作用会带来多种多样的气动弹性现象，包括气动弹性变形、颤振、阵风响应等。飞行器在飞行中遭遇阵风扰动后机体结构会产生颠簸，干扰驾驶员操作，影响乘坐舒适性，增加气动弹性载荷，缩短结构疲劳寿命，严重威胁飞行安全。随着现代飞行器设计要求的逐渐提高，复合材料使用率及结构柔性的增大使得飞行器对阵风扰动的敏感性进一步提高，阵风扰动对于飞行器的影响更加复杂。阵风减缓分析研究对于提高现代飞行器的飞行性能，保障飞行安全具有重要意义。

利用势流理论面元法进行气动力计算是飞行器气动载荷分析的重要方法。面元法在求解精度上相较二维气动力建模方法提高很多且计算效率高，满足理论研究及工程应用需求。偶极子格网法和涡格法是两类常用的气动力计算方法，能够分别在频域和时域上给出气动力计算结果，在气动弹性分析特别是阵风问题分析上应用广泛。

基于偶极子格网法的阵风分析方法是飞行器气动弹性阵风分析的传统方法吗，由于偶极子格网法本身仅能提供频域气动力计算结果，为得到便于系统仿真及控制器设计的状态空间形式动力学方程，需要在偶极子格网法求得的广义非定常气动力影响系数矩阵的基础上进行有理函数拟合，得到气动力时域表达形式。即该方法需要先建立频域气动力模型，再进行频域到时域的拟合，计算较繁琐；同时，有理函数拟合的准确性取决于气动力滞后根的选取，选取依赖经验及反复试验，容易出现较大误差。

涡格法相较于偶极子格网法具有优势。涡格法是时域三维气动力分析方法，不需要简谐振动假设，可以直接计算任意运动下的非定常气动力。传统基于涡格法进行弹性飞行器阵风响应分析主要基于松耦合迭代算法进行结构与气动模型的耦合求解。在该分析框架下针对飞行器全机或部件进行阵风响应分析，计算结果精度令人满意。但随着分析模型复杂性的增加，基于涡格法的松耦合迭代求解方法计算量大，分析效率大幅度降低。同时，松耦合迭代算法也使得分析方法难以应用于控制器设计，无法进行阵风减缓研究。综上，目前气动弹性阵风分析中主流的基于(1)偶极子格网法- 有理函数拟合方法及(2)涡格法-松耦合迭代方法的两类方法均有不足。

状态空间形式的涡格法可以将涡格法气动力模型写为状态空间形式同时不损失精度，Breuker等人在其已发表论文中给出了状态空间形式涡格法的表达形式及基于简单梁模型与状态空间涡格法的气动弹性方程，并未将其应用于基于有限元方法的结构模型，使用范围受到很大限制。同时，已发表文献及专利中未有涉及状态空间形式涡格法能够将气动模型写为状态空间形式，在控制仿真及设计上具有很大优势。

发明内容

为了克服传统气动弹性阵风分析中(1)偶极子格网法-有理函数拟合方法需要进行有理函数拟合导致精度下降、步骤繁琐及(2)涡格法-松耦合迭代方法计算效率低，不便于控制系统设计的不足。已有的基于状态空间形式涡格法的阵风分析不能适用于有限元复杂模型同时未有扩展至阵风减缓控制问题，本发明基于状态空间形式的涡格法建立气动弹性阵风响应计算及气动弹性阵风减缓分析框架，提出一种高效准确的气动弹性阵风响应计算方法，建立气动弹性阵风减缓分析方法。本发明的具体技术方案如下：

一种基于状态空间形式涡格法的阵风减缓分析方法，包括以下步骤：

S1：计算初始化；

建立气动模型，在机翼中弧面沿弦向和展向划分四边形气动网格，气动网格包括在机翼表面的附着涡网格及沿来流方向拖出的尾涡网格；

S2：对步骤S1得到的气动网格建立状态空间形式的气动力方程；

S2-1：在气动网格每个1/4弦线处布置涡线段，四段等强度涡线首尾相连构成一个涡环；选择1/4弦线中点为力的作用点，选择3/4弦线中点为控制点；

S2-2：将气动模型中的涡分为三部分：机翼表面的附着涡、机翼后缘的第一排尾涡以及其他尾涡，设机翼表面附着涡强度列向量为Γ_b，机翼后缘第一排尾涡强度为Γ_w0，其他尾涡强度为Γ_wl，气动控制方程为：

K_bΓ_b+K_w0Γ_w0+K_wlΓ_wl＝-w (1)

其中，w＝(V_∞+V_g-V_b)·n为机翼表面诱导速度，V_∞为来流速度矢量，V_g为来流阵风扰动速度矢量及V_b为壁面运动速度矢量，n为控制点处的法向量列阵；K_b,K_w0,K_wl分别为附着涡、第一排尾涡及其他尾涡的诱导系数矩阵；

S2-3：机翼后缘的第一排尾涡在脱出的过程中保持涡强守恒，关系表达为：

其中，

为机翼后缘第一排尾涡强度Γ_w0关于时间的导数，Δt为时间步长，C₁为保证机翼后缘第一排尾涡与机翼表面的附着涡对应关系正确的系数矩阵，只包含0和1 两种元素；

S2-4：规定其他尾涡脱出后保持涡强不变，关系表达为：

其中，

为其他尾涡强度Γ_wl关于时间的导数，C₂,C₃分别为表征其他尾涡与机翼后缘第一排尾涡位置对应关系的常数提取矩阵，只包含0和1两种元素；

S2-5：综合式(1)(2)(3)得到涡格法状态空间方程形式：

其中，Γ_w＝[Γ_w0 Γ_wl]^T，A_a,B_a为状态空间系数矩阵，仅与气动面几何形状及机翼表面的附着涡、机翼后缘的第一排尾涡以及其他尾涡三者的划分有关：

其中，I为单位阵，O为只包含0元素的零矩阵；

S2-6：给出网格上的气动力表达，气动网格上的压强差表达为：

其中，下标i，j表示沿展向第i个，沿弦向第j个网格中对应物理量矢量或矩阵的分量，Δp_ij为气动网格上压强差分量，V_l,ij为气动网格上当地来流速度分量，Γ_b,ij为气动网格上表面附着涡的涡强分量，ρ_∞为来流大气密度，τ₁,τ₂分别为气动网格中控制点处沿当地速向及弦向方向的切向量；

假设气动面扰动速度远小于来流速度，即认为V_l＝V_∞，方程(7)中涡强的变化量由下述差分方程求解：

作用在气动网格上的气动力表示为：

F_ij＝Δp_ijS_ijn (9)

其中，Δc_ij为气动网格沿弦向几何长度，Δb_ij为气动网格沿展向几何长度，n为控制点处的法向量列阵，F_ij为气动网格上的气动力向量，S_ij为气动网格面积；综合方程 (1)(7)(8)(9)得到气动力矢量F表达式：

其中，B₁,B₂,B₃为系数矩阵；

S2-7：涡格法状态空间方程完整表达为：

状态变量x_aero＝Γ，输入量

输出量y_aero＝F；A_aero＝A_a， B_aero＝[B_aO]，C_aero＝B₃，D_aero＝[B₁ B₂]；

状态空间形式涡格法气动力模型即方程(11)所示，给定输入量后即能够求解得到气动力数值，输入量中包含阵风扰动量，因此适用于计算阵风问题；

S3：建立结构有限元模型；

利用有限元方法建立结构模型，同时给出模态坐标下的结构动力学方程：

其中，M_s,B_s,K_s分别为广义质量阵、阻尼阵及刚度阵，q_s为广义模态坐标，F_s为广义力；模态坐标与物理坐标的关系为：

u_s＝Φq_s (13)

其中，u_s为物理空间坐标，Φ为模态矩阵；

S4：利用三维样条插值理论得到结构-气动耦合关系；

步骤S2得到气动力方程，步骤S3得到结构动力学方程，气动弹性计算需要将气动与结构计算耦合在一起，利用三维样条插值理论耦合有限元模态广义自由度与气动面控制点法向量变化量，表征气动力方程与结构动力学方程之间的耦合关系；

气动面位移W_A与结构位移W_S关系为：

W_A＝EW_S (14)

其中，E为位移插值矩阵，同样的，气动面速度V_A插值能够由结构速度V_S与映射关系求解：

V_A＝EV_S (15)

由某一结构位移u_S导致的气动面上任意点切向量τ的变化Δτ同样能够由映射关系求解：

Δτ＝E_T0u_S (16)

其中，E_T0为物理位移到切向量变化的插值矩阵，考虑物理位移与模态空间广义自由度关系(13)，则有：

Δτ＝E_Tq_S (17)

其中，E_T＝E_T0Φu_S为模态空间广义自由度到切向量变化的插值矩阵；

考虑结构变形后空间法向量与两个切向量的叉乘关系：

其中，n为控制点处的法向量列阵，n₀,τ₁,τ₂分别为未变形状态下控制点处的法向量、气动网格中控制点沿当地速向的切向量及气动网格中控制点沿当地弦向方向的切向量，忽略二阶小量Δτ₁×Δτ₂,结合方程(17)得到变形后法向量与模态空间广义自由度关系：

n＝n₀+Δn＝n₀+E_Nq_S (19)

其中，E_N为模态空间广义自由度即结构模态坐标到法向量变化量的插值矩阵；

考虑结构变形后，诱导速度表达为：

忽略扰动速度与法向量增量相乘形成的高阶项，诱导速度表达为：

w＝V_∞n₀+V_gn₀-V_bn₀+V_∞Δn (21)

壁面运动速度V_b完全由弹性振动导致：

综合方程(19)(21)(22)得到耦合关系：

S5：给出气动弹性系统状态空间方程；

通过步骤S4的结构-气动耦合关系耦合步骤S2中的气动力方程及步骤S3中的结构动力学方程组装得到气动弹性系统状态空间方程；

改写方程(23)为：

其中，

R_1s＝EΦ,R_2s＝V_∞E_N，R₀为来流速度在初始构型下导致的诱导速度量，R_g为单位阵风速度导致的诱导速度系数矩阵，R_1s为模态坐标变化单位量导致的诱导速度系数矩阵，R_2s为模态坐标时间导数变化单位量导致的诱导速度系数矩阵，v_g为空间流场中的阵风速度向量；

为便于引入控制面自由度，引入舵面偏转模态Φ_C，舵面偏转模态与弹性振动模态处理方式一致，即将舵面偏转带来的气动力变化通过对诱导速度处法向量变化Δn 及运动速度V_b的影响体现，则考虑舵面偏转δ的诱导速度为：

其中，

右上角标C表示该矩阵或矢量对应舵面偏转模态即

为舵面偏转单位量导致的诱导速度系数矩阵，

为舵偏偏转时间导数单位量导致的诱导速度系数矩阵，

为为舵面偏转模态空间广义自由度到法向量变化量的插值矩阵；

结合方程(10)，结构广义外力为：

其中，E_F为力插值矩阵，综合方程(25)(26)：

其中，K₀,K₁,K₂,K₃,K₄,K₅,K₆,K₇,K₈,K₉为对应的系数矩阵，结合涡格法状态空间方程(4)，气动力控制方程为：

综合方程(12)(27)(28)得到气动弹性系统状态空间方程：

其中，A₀，A₁，A₂，A₃，A₄，A₅分别为对应初始状态，结构模态广义坐标，结构模态广义坐标时间导数，阵风速度向量，舵面偏转，舵面偏转时间导数的系数矩阵，A,B 为状态空间模型系数矩阵，

考虑结构信号输出，常用为传感器安置点加速度信号，则：

其中，Φ_r为传感器位置的模态矩阵；

指定阵风模型及对应阵风速度量v_g，

即能够通过(29)(30)仿真得到系统阵风响应；

S6：引入控制系统；

选择结构某位置处传感器测量的信号作为反馈信号，回传控制器，通过控制增益环节计算后转化为控制信号传给舵机，驱动控制面，控制系统方程为：

其中K_con为控制增益矩阵；

方程(29)(30)(31)构成气动弹性阵风减缓控制系统，能够在Matlab软件Simulink模块中搭建实现。

本发明的有益效果在于：

1.传统气动弹性阵风减缓分析方法存在诸多不足。考虑到状态空间形式涡格法具有计算准确、状态空间形式便于控制器设计等优势，本发明对飞行器实际外形建立状态空间形式涡格法控制方程，采用有限元方法建立结构模型，使得本发明具有更好的模型适用性。

2.本发明给出了结构动力学方程与气动控制方程的耦合关系，将状态空间形式涡格法的尾涡涡强变量与有限元模型模态空间自由度耦合在一起，为气动弹性系统状态空间方程推导奠定了基础。

3.本发明给出了气动弹性系统状态空间方程，引入控制环节综合气动弹性系统状态空间方程得到气动弹性阵风减缓分析方程，将状态空间涡格法应用于气动弹性阵风减缓分析中，并适用于有限元方法建立的模型。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，通过参考附图会更加清楚的理解本发明的特征和优点，附图是示意性的而不应理解为对本发明进行任何限制，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，可以根据这些附图获得其他的附图。其中：

图1是本发明分析模型建立流程图；

图2是本发明控制思路；

图3是实施例的飞翼全机有限元模型示意图；

图4是实施例的拆分后机翼及机身有限元模型示意图；

图5是实施例的阵风减缓分析结果。

具体实施方式

为了能够更清楚地理解本发明的上述目的、特征和优点，下面结合附图和具体实施方式对本发明进行进一步的详细描述。需要说明的是，在不冲突的情况下，本发明的实施例及实施例中的特征可以相互组合。

在下面的描述中阐述了很多具体细节以便于充分理解本发明，但是，本发明还可以采用其他不同于在此描述的其他方式来实施，因此，本发明的保护范围并不受下面公开的具体实施例的限制。

如图1-2所示，本发明依据给定模型建立气动网格，根据状态空间形式涡格法建立模型附着涡涡强、尾涡涡强与机翼表面诱导速度的关系，插值得到有限元模态及舵面模态广义坐标与机翼表面诱导速度的关系后形成气动弹性耦合方程，气动弹性耦合方程同样为状态空间形式，给定阵风模型作为外激励，可进行气动弹性阵风分析；在气动弹性耦合方程的基础上以机翼某点加速度作为系统输出，经过控制器后形成舵面偏转输入信号，搭建闭环控制系统，进行气动弹性阵风减缓仿真分析。

为了方便理解本发明的上述技术方案，以下通过具体实施例对本发明的上述技术方案进行详细说明。

实施例1

采用一大展弦机翼模型，机翼模型主梁为翼尖到翼根线性增大的十字梁。主梁位置位于弦长的40％。主梁在MSC.Nastran中以CBEAM梁单元模拟，翼肋以大刚度梁单元模拟，质量特性以集中质量点模拟。模型部分参数如表1所示。

表1飞翼模型参数

参数名	参数值
		机翼展长/m	1.542
翼根弦长/m	0.261
		翼梢弦长/m	0.069
扭转角/Deg	-2.0
		展弦比	9.3
质量/kg	2.8199
		翼型	超临界翼型

S1：计算初始化；

建立气动模型，在机翼中弧面沿弦向划分5个网格，展向划分40个网格，拖出的尾涡网格4弦向划分50个网格，展向划分40个网格；机翼网格划分如图3所示。

S2：对步骤S1得到的气动网格建立状态空间形式的气动力方程；

涡格法状态空间方程完整表达为：

状态变量x_aero＝Γ_w，输入量

输出量y_aero＝F；A_aero＝A_a， B_aero＝[B_aO]，C_aero＝B₃，D_aero＝[B₁ B₂]。

S3：建立结构有限元模型；

利用有限元方法建立结构模型，有限元模型如图4所示。有限元模型在MSC.Nastran软件中建立，计算中采用前6阶模态给出模态空间结构动力学方程，模态信息如表2所示。

表2模态信息表

S4：利用三维样条插值理论得到结构-气动耦合关系；

综合方程(19)(21)(22)得到耦合关系：

S5：给出气动弹性系统状态空间方程；

为便于引入控制面自由度，引入舵面偏转模态Φ_C，面积为0.0000168m²，平均弦长为0.04m，据翼根距离0.8045m。

本实施例中采用“Sin”型离散阵风模型，数学表达式为：

其中，w_g为阵风速度，w_g0为阵风速度幅值，L_g为阵风尺度，t为仿真时间，x_g为飞行器位置，x₀为飞行器初始位置，一般可置零，V为来流速度。离散阵风指定为一维阵风，在展向保持均匀不变，阵风速度垂直于来流方向，考虑前述状态空间形式表达，阵风加载形式为：

v_g＝[0 0 w_g(l)]^T (36)

综合方程(12)(27)(28)得到气动弹性系统状态空间方程：

其中，A,B为状态空间模型系数矩阵，

考虑结构信号输出，为加速度传感器信号，加速度传感器布置于距翼尖1/3翼展处，则：

其中，Φ_r为加速度传感器位置的模态矩阵。

S6：引入控制系统；

选择加速度传感器信号作为反馈信号，回传控制器，通过控制增益环节计算后转化为控制信号传给舵机，驱动控制面。本实施例采用经典PID控制律进行阵风减缓分析。

本实施例阵风减缓分析在Matlab软件Simulink模块中搭建实现。图5为阵风减缓分析结果，开环状态表示为引入控制系统的机翼阵风响应分析结果，闭环状态表示引入控制系统后的机翼阵风减缓分析结果，由图可知，本发明基于状态空间形式涡格法及结构有限元模型构建气动弹性系统阵风响应分析的方法，能够提高阵风分析精度，避免频域时域转换，提高分析效率。

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (1)

1.一种基于状态空间形式涡格法的阵风减缓分析方法，其特征在于，包括以下步骤：

S1：计算初始化；

建立气动模型，在机翼中弧面沿弦向和展向划分四边形气动网格，气动网格包括在机翼表面的附着涡网格及沿来流方向拖出的尾涡网格；

S2：对步骤S1得到的气动网格建立状态空间形式的气动力方程；

S2-1：在气动网格每个1/4弦线处布置涡线段，四段等强度涡线首尾相连构成一个涡环；选择1/4弦线中点为力的作用点，选择3/4弦线中点为控制点；

S2-2：将气动模型中的涡分为三部分：机翼表面的附着涡、机翼后缘的第一排尾涡以及其他尾涡，设机翼表面附着涡强度列向量为Γ_b，机翼后缘第一排尾涡强度为Γ_w0，其他尾涡强度为Γ_wl，气动控制方程为：

K_bΓ_b+K_w0Γ_w0+K_wlΓ_wl＝-w (1)

其中，w＝(V_∞+V_g-V_b)·n为机翼表面诱导速度，V_∞为来流速度矢量，V_g为来流阵风扰动速度矢量及V_b为壁面运动速度矢量，n为控制点处的法向量列阵；K_b,K_w0,K_wl分别为附着涡、第一排尾涡及其他尾涡的诱导系数矩阵；

S2-3：机翼后缘的第一排尾涡在脱出的过程中保持涡强守恒，关系表达为：

其中，

为机翼后缘第一排尾涡强度Γ_w0关于时间的导数，Δt为时间步长，C₁为保证机翼后缘第一排尾涡与机翼表面的附着涡对应关系正确的系数矩阵，只包含0和1两种元素；

S2-4：规定其他尾涡脱出后保持涡强不变，关系表达为：

其中，

为其他尾涡强度Γ_wl关于时间的导数，C₂,C₃分别为表征其他尾涡与机翼后缘第一排尾涡位置对应关系的常数提取矩阵，只包含0和1两种元素；

S2-5：综合式(1)(2)(3)得到涡格法状态空间方程形式：

其中，Γ_w＝[Γ_w0 Γ_wl]^T，A_a,B_a为状态空间系数矩阵，仅与气动面几何形状及机翼表面的附着涡、机翼后缘的第一排尾涡以及其他尾涡三者的划分有关：

其中，I为单位阵，O为只包含0元素的零矩阵；

S2-6：给出网格上的气动力表达，气动网格上的压强差表达为：

其中，下标i，j表示沿展向第i个，沿弦向第j个网格中对应物理量矢量或矩阵的分量，Δp_ij为气动网格上压强差分量，V_l,ij为气动网格上当地来流速度分量，Γ_b,ij为气动网格上表面附着涡的涡强分量，ρ_∞为来流大气密度，τ₁,τ₂分别为气动网格中控制点处沿当地速向及弦向方向的切向量；

假设气动面扰动速度远小于来流速度，即认为V_l＝V_∞，方程(7)中涡强的变化量由下述差分方程求解：

作用在气动网格上的气动力表示为：

F_ij＝Δp_ijS_ijn (9)

其中，Δc_ij为气动网格沿弦向几何长度，Δb_ij为气动网格沿展向几何长度，n为控制点处的法向量列阵，F_ij为气动网格上的气动力向量，S_ij为气动网格面积；综合方程(1)(7)(8)(9)得到气动力矢量F表达式：

其中，B₁,B₂,B₃为系数矩阵；

S2-7：涡格法状态空间方程完整表达为：

状态变量x_aero＝Γ_w，输入量

输出量y_aero＝F；A_aero＝A_a，B_aero＝[B_a O]，C_aero＝B₃，D_aero＝[B₁ B₂]；

状态空间形式涡格法气动力模型即方程(11)所示，给定输入量后即能够求解得到气动力数值，输入量中包含阵风扰动量，因此适用于计算阵风问题；

S3：建立结构有限元模型；

利用有限元方法建立结构模型，同时给出模态坐标下的结构动力学方程：

其中，M_s,B_s,K_s分别为广义质量阵、阻尼阵及刚度阵，q_s为广义模态坐标，F_s为广义力；模态坐标与物理坐标的关系为：

u_s＝Φq_s (13)

其中，u_s为物理空间坐标，Φ为模态矩阵；

S4：利用三维样条插值理论得到结构-气动耦合关系；

步骤S2得到气动力方程，步骤S3得到结构动力学方程，气动弹性计算需要将气动与结构计算耦合在一起，利用三维样条插值理论耦合有限元模态广义自由度与气动面控制点法向量变化量，表征气动力方程与结构动力学方程之间的耦合关系；

气动面位移W_A与结构位移W_S关系为：

W_A＝EW_S (14)

其中，E为位移插值矩阵，同样的，气动面速度V_A插值能够由结构速度V_S与映射关系求解：

V_A＝EV_S (15)

由某一结构位移u_S导致的气动面上任意点切向量τ的变化Δτ同样能够由映射关系求解：

Δτ＝E_T0u_S (16)

其中，E_T0为物理位移到切向量变化的插值矩阵，考虑物理位移与模态空间广义自由度关系(13)，则有：

Δτ＝E_Tq_S (17)

其中，E_T＝E_T0Φu_S为模态空间广义自由度到切向量变化的插值矩阵；

考虑结构变形后空间法向量与两个切向量的叉乘关系：

其中，n为控制点处的法向量列阵，n₀,τ₁,τ₂分别为未变形状态下控制点处的法向量、气动网格中控制点沿当地速向的切向量及气动网格中控制点沿当地弦向方向的切向量，忽略二阶小量Δτ₁×Δτ₂,结合方程(17)得到变形后法向量与模态空间广义自由度关系：

n＝n₀+Δn＝n₀+E_Nq_S (19)

其中，E_N为模态空间广义自由度即结构模态坐标到法向量变化量的插值矩阵；

考虑结构变形后，诱导速度表达为：

忽略扰动速度与法向量增量相乘形成的高阶项，诱导速度表达为：

w＝V_∞n₀+V_gn₀-V_bn₀+V_∞Δn (21)

壁面运动速度V_b完全由弹性振动导致：

综合方程(19)(21)(22)得到耦合关系：

S5：给出气动弹性系统状态空间方程；

通过步骤S4的结构-气动耦合关系耦合步骤S2中的气动力方程及步骤S3中的结构动力学方程组装得到气动弹性系统状态空间方程；

改写方程(23)为：

其中，R₀＝V_∞n₀,

R_1s＝EΦ,R_2s＝V_∞E_N，R₀为来流速度在初始构型下导致的诱导速度量，R_g为单位阵风速度导致的诱导速度系数矩阵，R_1s为模态坐标变化单位量导致的诱导速度系数矩阵，R_2s为模态坐标时间导数变化单位量导致的诱导速度系数矩阵，v_g为空间流场中的阵风速度向量；

为便于引入控制面自由度，引入舵面偏转模态Φ_C，舵面偏转模态与弹性振动模态处理方式一致，即将舵面偏转带来的气动力变化通过对诱导速度处法向量变化Δn及运动速度V_b的影响体现，则考虑舵面偏转δ的诱导速度为：

其中，

右上角标C表示该矩阵或矢量对应舵面偏转模态即

为舵面偏转单位量导致的诱导速度系数矩阵，

为舵面 偏转时间导数单位量导致的诱导速度系数矩阵，

为舵面偏转模态空间广义自由度到法向量变化量的插值矩阵；

结合方程(10)，结构广义外力为：

其中，E_F为力插值矩阵，综合方程(25)(26)：

其中，K₀,K₁,K₂,K₃,K₄,K₅,K₆,K₇,K₈,K₉为对应的系数矩阵，结合涡格法状态空间方程(4)，气动力控制方程为：

综合方程(12)(27)(28)得到气动弹性系统状态空间方程：

其中，A₀，A₁，A₂，A₃，A₄，A₅分别为对应初始状态，结构模态广义坐标，结构模态广义坐标时间导数，阵风速度向量，舵面偏转，舵面偏转时间导数的系数矩阵，A,B为状态空间模型系数矩阵，

考虑结构信号输出，常用为传感器安置点加速度信号，则：

其中，Φ_r为传感器位置的模态矩阵；

指定阵风模型及对应阵风速度量v_g，

即能够通过(29)(30)仿真得到系统阵风响应；

S6：引入控制系统；

选择结构某位置处传感器测量的信号作为反馈信号，回传控制器，通过控制增益环节计算后转化为控制信号传给舵机，驱动控制面，控制系统方程为：

其中K_con为控制增益矩阵；

方程(29)(30)(31)构成气动弹性阵风减缓控制系统，能够在Matlab软件Simulink模块中搭建实现。

Images