[发明专利]一种基于状态空间形式涡格法的阵风减缓分析方法

权利要求书

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模块中搭建实现。