CN110929336B - 基于多体系统传递矩阵法求解三维机翼线性颤振速度的方法 - Google Patents

Abstract

本发明公开了一种基于多体系统传递矩阵法求解三维机翼线性颤振速度的方法，包括以下过程：基于多体系统传递矩阵法推导机翼的弯扭耦合梁传递矩阵；建立机翼的总传递方程，并求解得到机翼的固有频率和振型；根据Theodorsen非定常气动力理论建立机翼振动位移与气动力的关系；建立机翼的体动力学方程，并将其转换到频域获得机翼的颤振频域模型；对机翼颤振频域模型进行频域求解，得到机翼的颤振速度。本发明可实现机翼颤振速度的快速求解。

Description

基于多体系统传递矩阵法求解三维机翼线性颤振速度的方法

技术领域

本发明属于气动弹性仿真技术领域，具体涉及一种基于多体系统传递矩阵法求解三维机翼线性颤振速度的方法。

背景技术

随着航空工业的不断发展，飞机飞行速度的不断提高，气动弹性问题越发受到人们重视。在飞机设计技术发展的早期，多次出现由于气动弹性而导致的严重事故。经过几十年的发展，气动弹性力学日趋成熟，相应的分析方法在航空航天领域起到了重要的作用。气动弹性问题研究的每一步进展都能对航空技术的发展产生极大的促进，对气动弹性问题的认识也不断深入，解决气动弹性问题的能力也在不断增强。即便如此，现代的飞行器依然没有摆脱气动弹性力学问题的困扰，特别是颤振事故频繁出现。颤振无疑是所有气动弹性问题中最为突出的一个。在航空工程和风工程领域中，飞机机翼、风力机叶片等结构在高流速下，结构在气动力的激励下常常会产生颤振现象。颤振的产生将使结构产生疲劳损耗，甚至在短时间内导致结构破坏。对颤振问题的研究对风力机与飞机的设计制造具有一定实际意义。并且随着飞机柔性增大，操作面以及尾翼更容易引起颤振。这些操作面、尾翼系统本质上是多体系统，对这些系统结构的气动弹性求解，是现在面临的主要问题。

发明内容

本发明的目的在于克服现有技术中的不足，提出了一种基于多体系统传递矩阵法求解三维机翼线性颤振速度的方法，解决了对机翼振动特性求解效率低、速度慢的问题，实现机翼颤振速度的快速求解。

为解决上述技术问题，本发明提供了一种基于多体系统传递矩阵法求解三维机翼线性颤振速度的方法，其特征是，包括以下过程：

基于多体系统传递矩阵法推导机翼的弯扭耦合梁传递矩阵；

建立机翼的总传递方程，并求解得到机翼的固有频率和振型；

根据Theodorsen非定常气动力理论建立机翼振动位移与气动力的关系；

建立机翼的体动力学方程，并将其转换到频域获得机翼的颤振频域模型；

对机翼颤振频域模型进行频域求解，得到机翼的颤振速度。

进一步的，所述基于多体系统传递矩阵法推导机翼的弯扭耦合梁传递矩阵包括：

建立弯扭耦合梁的振动微分方程：

式中，m为线质量，I_α为单位长度转动惯量，EI为弯曲刚度，GJ为扭转刚度，y为机翼弯曲方向位移，θ_x为扭转角度，t为时间，b为机翼弦长的一半，x_α为质心到弹性轴的距离，且当质心在Z正方向，x_α为正；

将物理坐标转化为模态坐标，令

y(x,t)＝Y(x)sinωt,θ_x(x,t)＝Θ_x(x)sinωt (2)

式中，ω为圆频率，Y(x)、Θ_x(x)为y(x,t)、θ_x(x,t)通过模态坐标转换得到；

将(2)带入(1)中，(1)式变为：

消去式(3)中的Y(x)或者Θ_x(x)，得到

式中，

W＝Y或者Θ (5)引入无量纲长度：

ξ＝x/L (6)

式(6)中L为划分后单元弯扭耦合梁长度；

式(4)可改写成无量纲形式：

(D⁶+aD⁴-bD²-abc)W＝0 (7)

其中，

六阶微分方程(7)的通解可以表示为：

W(ξ)＝C₁ coshαξ+C₂ sinhαξ+C₃ cosβξ+C₄ sinβξ+C₅ cosγξ+C₆ sinγξ (9)

式中C₁,......,C₆是常数，并且

q＝b+a²/3

式(9)中的W(ξ)代表弯曲位移Y和扭转角度Θ_x在不同常数下的解；因此，

Y(ξ)＝A₁ coshαξ+A₂ sinhαξ+A₃ cosβξ+A₄ sinβξ+A₅ cosγξ+A₆ sinγξ (11)

Θ_x(ξ)＝B₁ coshαξ+B₂ sinhαξ+B₃ cosβξ+B₄ sinβξ+B₅ cosγξ+B₆sinγξ(12)

式中A₁,......,A₆和B₁,......,B₆是两组不同的常数；

将式(11)和(12)代入式(3)可以确定常数有如下规律：

式中，

从式(11)、(12)可以得到弯曲角度Θ_z(ξ)，弯矩M_z(ξ)，剪切力Q_y(ξ)和扭矩M_x(ξ)的表达式：

为了推导弯扭耦合梁的传递矩阵，由式(11)、(12)、(15)、(16)、(17)、(18)整理可得

即Z(ξ)＝B(ξ)·a,a＝[A₁,A₂,A₃,A₄,A₅,A₆]^T，其中Z(ξ)为状态矢量，即方程(19)等式左边项，B(ξ)代表方程(19)等式右边的第一个矩阵；令ξ＝0，且根据方程(19)和多体系统传递矩阵法可将方程(19)写成Z_I＝B(0)·a形式；令ξ＝1得到Z_O＝B(1)·a＝B(1)B(0)^-1Z_I＝U_iZ_I；其中Z_O代表模态坐标下元件输出端的状态矢量，Z_I代表模态坐标下元件输入端的状态矢量，U_i为第i段的弯扭耦合梁的传递矩阵；

所以第i段的弯扭耦合梁的传递矩阵U_i为：

U_i＝B(1)·B^-1(0) (20)

其中，

进一步的，所述建立机翼的总传递方程包括：

根据各段弯扭耦合梁传递矩阵，建立机翼总传递方程为：

Z_o＝U_allZ_i (21)

式中U_all＝U_n···U₃U₂U₁，U₁,U₂,U₃···U_n分别代表机翼各段弯扭耦合梁的传递矩阵，Z_i,Z_o分别代表系统模态坐标下输入端、输出端的状态矢量；输入端Z_i＝[0,0,M_z,Q_y,0,M_x]^T，输出端Z_o＝[Y,Θ_z,0,0,Θ_x,0]^T。

进一步的，所述根据Theodorsen非定常气动力理论建立机翼振动位移与气动力的关系包括：

根据机翼坐标系，使用Theodorsen公式计算机翼气动力：

式(22)、(23)中，L表示升力，T_α表示俯仰力矩，U为来流速度，ρ为空气密度，b为半弦长，a为刚心距离弦线中点的距离与半弦长的比值，并在中点后为正，C(K)为Theodorsen函数，其值与折合频率K＝ωb/U有关，ω为圆频率，y表示沉浮位移，

为速度，

为加速度；θ_x表示俯仰位移，

表示俯仰速度，

表示俯仰加速度；

假设Goland机翼各方向的运动为简谐运动，即：

y＝y_se^iωt；θ＝θ_se^iωt (24)

式中，y、θ分别表示物理坐标系下的沉浮和俯仰位移；y_s、θ_s分别表示模态坐标下的沉浮和俯仰位移，ω表示圆频率，t表示时间；

将该式代入到式(22)、(23)中，可以整理得：

令f＝[L₁,T_a1,L₂,T_a2,…L_n,T_an]^T，n表示将机翼分为n段，L_n、T_αn表示机翼第n段上面的升力和俯仰力矩；令v＝v_se^iωt，利用式(25)、(26)可以建立以下位移与气动力的关系：

f＝Bv_se^iωt (27)

式(27)中，

B＝[B₁,B₂,…B_n]^T，v_s＝[v_s1,v_s2,…v_sn]^T (28)

其中，

k_a＝-1+i·(2Π/π)·C(K)/K

k_b＝-0.5+i·(1+2Π/π·C(K))/K+(2Π/π)·C(K)/K²

m_a＝0.5

b_n、l_n、a_n分别表示机翼第n段的半弦长、长度、刚心距离弦线中点的距离与半弦长的比值。B_n为机翼第n段的气动力矩阵；v_sn表示机翼第n段模态坐标下的沉浮和俯仰位移。

进一步的，建立机翼的体动力学方程，并将其转换到频域获得机翼的颤振频域模型包括：

根据多体系统传递矩阵法，建立机翼的体动力学方程为：

Mv_tt+Cv_t+Kv＝f (29)

其中，M、C、K分别为质量矩阵、阻尼矩阵、刚度矩阵；v_t,v_tt分别为v对时间t的一阶导数和二阶导数；

结合公式27并且将v＝v_se^iωt代入式(29)中，则有：

-ω²Mv_s+iωCv_s+Kv_s＝Bv_s (30)

由于体动力学方程式29右边的外力项非解耦，把振型写为V＝[V¹,V²]，其中V¹为机翼的第1阶振型，V²为机翼的第2阶振型；令式(30)中的v_s＝Vq，q为机翼振动位移的广义坐标，得到：

-ω²MVq+iωCVq+KVq＝BVq (31)

由于增广特征矢量的正交性(左乘V^T)，因此有：

-ω²V^TMVq+iωV^TCVq+V^TKVq＝V^TBVq (32)

式中，

一般计算中忽略结构的阻尼，因此

其中ω₁、ω₂分别为机翼结构的第一阶圆频率和第二阶圆频率。

进一步的，所述对机翼颤振频域模型进行频域求解包括：

采用U-g法对机翼颤振频域模型进行频域求解。

进一步的，所述采用U-g法对机翼颤振频域模型进行频域求解包括：

机翼颤振速度的求解过程被转化为求解该复特征值的问题，其特征值为：

式中，λ为特征值；λ_Re、λ_Im分别为特征值的实部、虚部。

由此可以得到：

式中，g、ω和U分别为人工阻尼、圆频率和速度；

最后，可以通过以下过程来实现颤振分析：

1)给定密度ρ，预设一折合频率K；

2)在指定折合频率K下计算式(36)，并得到该折合频率下的人工阻尼g、圆频率ω和速度U；

3)一定步长减小折合频率K，并重复第2步，计算各折合频率下的人工阻尼g、圆频率ω和速度U；

4)计算完所有的折合频率K下的数值后，将圆频率转换为频率，绘制相应的速度-人工阻尼U-g图和速度-频率U-f图；

5)从速度-人工阻尼U-g图中找到曲线与x轴相交的点，此处对应的速度U即为颤振速度U_F，在频率U-f图中对应于颤振速度U_F的点即为颤振频率f_F。

与现有技术相比，本发明所达到的有益效果：

(1)本发明提出一种基于多体系统传递矩阵法的机翼线性颤振频域模型快速建模思路，相比于计算流体力学方法过长的计算时间，更能满足工程上快速计算的需求；

(2)本发明采用多体系统传递矩阵法建立机翼动力学模型，能快速求解机翼的振动特性。

(3)本发明建立了机翼的气动弹性体动力学模型和仿真方法。该方法中的传递性、多体性表明该方法同样适用于三维非均直机翼、舵面系统等多体系统结构的气动弹性频域颤振计算，为机翼、舵面系统的结构设计提供参考。

附图说明

图1是本发明方法的流程示意图；

图2是Goland机翼模型示意图；

图3是弯扭耦合梁模型示意图；

图4是基于MSTMM计算出的Goland机翼圆频率计算结果图；

图5是基于MSTMM计算出的Goland机翼在y方向上的第一、第二阶振型图；

图6是基于MSTMM计算出的Goland机翼在θ_x方向上的第一、第二阶振型图；

图7是Goland机翼的速度-人工阻尼图；

图8是Goland机翼的速度-频率图。

具体实施方式

下面结合附图对本发明作进一步描述。以下实施例仅用于更加清楚地说明本发明的技术方案，而不能以此来限制本发明的保护范围。

现有技术中，多体系统传递矩阵法(MSTMM)由芮筱亭院士及其团队建立，用于多体动力学分析。该方法能实现一般多体系统动力学研究，其高计算效率方便运用于工程中。运用于本发明中对机翼振动特性的计算，可以实现机翼气动弹性模型的快速建模与快速计算。

本发明采用MSTMM计算出三维机翼振动特性，并根据MSTMM计算出的机翼振型，提出机翼线性颤振频域模型建模思路。同等计算条件下，采用CFD/CSD(计算流体力学/计算结构动力学双向流固耦合)尽管计算精度高，考虑因素全，但是计算时间过长，难以满足工程上快速计算的需求，本方法更方便运用于工程中对机翼临界颤振速度的快速计算，为实际机翼或其它舵面系统等多体系统结构的气动弹性计算提供参考。

本发明的一种基于多体系统传递矩阵法的三维机翼线性颤振频域模型建模方法，参见图1所示。具体实施步骤如下：

步骤一：基于多体系统传递矩阵法推导机翼的弯扭耦合梁传递矩阵；

本发明实施例中以Goland机翼为例，建立三维机翼线性颤振模型，并进行频域分析。Goland机翼模型如图2所示。Goland机翼虽然为均直机翼，但本发明方法通过对机翼分段建模的方法同样适用于非均直机翼或舵面系统等多体系统结构。

Goland机翼模型具体参数如下：机翼弦长c＝1.829m，半展长s＝6.096m，线质量m＝539.6kg/m，前缘到弹性轴的距离0.6096m，弹性轴到中心距离的无量纲量a＝-1/3，前缘到质心的距离0.7925m。建立空间直角坐标系如图2所示，以弹性轴为x轴，以机翼端面弦线方向为z轴，x轴与z轴交点为坐标原点，y轴过坐标原点垂直于xz平面，对z轴的单位展长的转动惯量I_z＝1.111kg·m²/m，对x轴的单位展长的转动惯量I_x＝112.2kg·m²/m，弯曲刚度EI＝9.76×10⁶N·m²，扭转刚度GJ＝9.88×10⁵N·m²。

根据机翼的变形特性，将机翼简化为n段弯扭耦合梁。如图3所示为机翼弯扭耦合梁模型。

根据图3所示的弯扭耦合梁模型中，以弯曲和扭转振动为主，忽略其剪切效应，建立弯扭耦合梁的振动微分方程：

式中，m为线质量，I_α为单位长度转动惯量，EI为弯曲刚度，GJ为扭转刚度，y为机翼弯曲方向位移，θ_x为扭转角度，t为时间，b为机翼弦长的一半，x_α为质心到弹性轴的距离，且当质心在Z正方向，x_α为正。

将物理坐标转化为模态坐标，令

y(x,t)＝Y(x)sinωt,θ_x(x,t)＝Θ_x(x)sinωt (2)

式中，ω为圆频率，Y(x)、Θ_x(x)为y(x,t)、θ_x(x,t)通过模态坐标转换得到。

将(2)带入(1)中，(1)式变为：

消去式(3)中的Y(x)或者Θ_x(x)，得到

式中，

W＝Y或者Θ (5)

引入无量纲长度：

ξ＝x/L (6)

式(6)中L为划分后单元弯扭耦合梁长度。

式(4)可改写成无量纲形式：

(D⁶+aD⁴-bD²-abc)W＝0 (7)

其中，

六阶微分方程(7)的通解可以表示为：

W(ξ)＝C₁ coshαξ+C₂ sinhαξ+C₃ cosβξ+C₄ sinβξ+C₅ cosγξ+C₆ sinγξ (9)

式中C₁,......,C₆是常数，并且

q＝b+a²/3

式(9)中的W(ξ)代表弯曲位移Y和扭转角度Θ_x在不同常数下的解。因此，

Y(ξ)＝A₁ coshαξ+A₂ sinhαξ+A₃ cosβξ+A₄ sinβξ+A₅ cosγξ+A₆ sinγξ (11)

Θ_x(ξ)＝B₁ coshαξ+B₂ sinhαξ+B₃ cosβξ+B₄ sinβξ+B₅ cosγξ+B₆ sinγξ(12)

式中A₁,......,A₆和B₁,......,B₆是两组不同的常数。

将式(11)和(12)代入式(3)可以确定常数有如下规律：

式中，

从式(11)、(12)可以得到弯曲角度Θ_z(ξ)，弯矩M_z(ξ)，剪切力Q_y(ξ)和扭矩M_x(ξ)的表达式：

为了推导弯扭耦合梁的传递矩阵，由式(11)、(12)、(15)、(16)、(17)、(18)整理可得

即Z(ξ)＝B(ξ)·a,a＝[A₁,A₂,A₃,A₄,A₅,A₆]^T，其中Z(ξ)为状态矢量，即方程(19)等式左边项，B(ξ)代表方程(19)等式右边的第一个矩阵。令ξ＝0，且根据方程(19)和多体系统传递矩阵法可将方程(19)写成Z_I＝B(0)·a形式。令ξ＝1得到Z_O＝B(1)·a＝B(1)B(0)^-1Z_I＝U_iZ_I。其中Z_O代表模态坐标下元件输出端的状态矢量，Z_I代表模态坐标下元件输入端的状态矢量，U_i为第i段的弯扭耦合梁的传递矩阵。

所以第i段的弯扭耦合梁的传递矩阵U_i为：

U_i＝B(1)·B^-1(0) (20)

其中，

U的下标i为上文中对机翼的分段，i表示具体的某一段，U的矩阵中各基本参数都表示第i段的参数，所以不标出i了，此处仅为推导传递矩阵形式，各段传递矩阵形式都相同，但参数根据不同形状的弯扭耦合梁都不同。

步骤二：建立机翼的总传递方程，并求解机翼的固有频率和振型；

状态矢量写为[Y,Θ_z,M_z,Q_y,Θ_x,M_x]^T，其中Y为沿坐标轴y位移对应的模态坐标列阵，Θ_z为该点处相对于平衡位置在z轴的角位移对应的模态坐标列阵，M_z为沿坐标轴z内力矩对应的模态坐标列阵，Q_y为沿坐标轴y内力对应的模态坐标列阵，Θ_x为该点处相对于平衡位置在x轴的角位移对应的模态坐标列阵，M_x为沿坐标轴x内力矩对应的模态坐标列阵。

根据各段弯扭耦合梁传递矩阵，建立机翼总传递方程为：

Z_o＝U_allZ_i (21)

式中U_all＝U_n···U₃U₂U₁，U₁,U₂,U₃···U_n分别代表机翼各段弯扭耦合梁的传递矩阵，Z_i,Z_o分别代表系统模态坐标下输入端、输出端的状态矢量；根据系统两端边界点，固定端位移为0，自由端力和力矩为0，因此可以根据边界条件确定：输入端Z_i＝[0,0,M_z,Q_y,0,M_x]^T，输出端Z_o＝[Y,Θ_z,0,0,Θ_x,0]^T。

求解机翼总传递方程，得到机翼的振动特性。其中计算得到机翼圆频率值如图4所示。图中横坐标为圆频率(圆频率即为频率*2π，频率和圆频率为固有频率的两种表达形式)，纵坐标表示|Δ|值的大小，当|Δ|值接近于0时即可以求出圆频率，每条竖线下对应的就是该阶模态的圆频率。如图5所示为机翼在y方向上的第一、第二阶振型，如图6所示为机翼在θ_x方向上的第一、第二阶振型。从图中可看出，第一阶振型在θ_x方向上几乎没有振动，说明第一阶振型以弯曲为主，第二阶振型以扭转为主。计算得到的固有频率和振型V，用于下一步三维机翼颤振模型建模中。

步骤三：使用Theodorsen非定常流理论模拟机翼气动力，建立机翼位移与气动力的关系；

根据机翼坐标系，使用Theodorsen公式计算机翼气动力(包括升力和俯仰力矩)：

式(22)、(23)中，L表示升力，T_α表示俯仰力矩，U为来流速度，ρ为空气密度，b为半弦长，a为刚心距离弦线中点的距离与半弦长的比值，并在中点后为正，C(K)为Theodorsen函数，其值与折合频率K＝ωb/U有关，ω为圆频率，y表示沉浮位移，

为速度，

为加速度；θ_x表示俯仰位移，

表示俯仰速度，

表示俯仰加速度。

假设Goland机翼各方向的运动为简谐运动，即：

y＝y_se^iωt；θ＝θ_se^iωt (24)

式中，y、θ分别表示物理坐标系下的沉浮和俯仰位移；y_s、θ_s分别表示模态坐标下的沉浮和俯仰位移，ω表示圆频率，t表示时间。

将该式代入到式(22)、(23)中，可以整理得：

令f＝[L₁,T_a1,L₂,T_a2,…L_n,T_an]^T，n表示将机翼分为n段，L_n、T_αn表示机翼第n段上面的升力和俯仰力矩；令v＝v_se^iωt，利用式(25)、(26)可以建立以下位移与气动力的关系：

f＝Bv_se^iωt (27)

式(27)中，

B＝[B₁,B₂,…B_n]^T，v_s＝[v_s1,v_s2,…v_sn]^T (28)

其中，

k_a＝-1+i·(2Π/π)·C(K)/K

k_b＝-0.5+i·(1+2Π/π·C(K))/K+(2Π/π)·C(K)/K²

m_a＝0.5

b_n、l_n、a_n分别表示机翼第n段的半弦长、长度、刚心距离弦线中点的距离与半弦长的比值。B_n为机翼第n段的气动力矩阵。v_sn表示机翼第n段模态坐标下的沉浮和俯仰位移。

步骤四：基于多体系统传递矩阵法和Theodorsen非定常气动力理论建立机翼的体动力学模型，即三维机翼的颤振频域模型；

根据多体系统传递矩阵法，建立机翼的体动力学方程为：

Mv_tt+Cv_t+Kv＝f (29)

其中，M、C、K分别为质量矩阵、阻尼矩阵、刚度矩阵。v_t,v_tt分别为v对时间t的一阶导数和二阶导数。

结合公式27并且将v＝v_se^iωt代入式(29)中，则有：

-ω²Mv_s+iωCv_s+Kv_s＝Bv_s (30)

由于体动力学方程式29右边的外力项非解耦，所以需要公式两边同时乘以所有振型，因此把振型写为V＝[V¹,V²]，其中V¹为机翼的第1阶振型，V²为机翼的第2阶振型。令式(30)中的v_s＝Vq，q为机翼振动位移的广义坐标，因为将物理坐标变换到广义坐标下，才能进一步进行频域求解，得到：

-ω²MVq+iωCVq+KVq＝BVq (31)

由于增广特征矢量的正交性(左乘V^T)，因此有：

-ω²V^TMVq+iωV^TCVq+V^TKVq＝V^TBVq (32)

式中，

一般计算中忽略结构的阻尼，因此

其中ω₁、ω₂分别为机翼结构的第一阶圆频率和第二阶圆频率。

步骤五：基于U-g法对三维机翼线性颤振模型进行频域求解；

对于线性颤振问题，可以利用U-g法对其进行求解。在该方法中需要引入人工结构阻尼g，在式(29)的右边添加人工阻尼力项，D＝-igKv_se^iωt，且假设机翼结构本身阻尼C＝0，因此式(29)最终可以转化为：

因此该机翼颤振速度的求解过程被转化为求解该复特征值的问题，其特征值为：

式中，λ为特征值；λ_Re、λ_Im分别为特征值的实部、虚部。

由此可以得到：

式中，g、ω和U分别为人工阻尼、圆频率和速度。

最后，可以通过以下过程来实现基于MSTMM的速度-人工阻尼U-g法颤振分析：

1)给定密度ρ，预设一折合频率K；

2)在指定折合频率K下计算式(36)，并得到该折合频率下的人工阻尼g、圆频率ω和速度U；

3)一定步长减小折合频率K，并重复第2步，计算各折合频率下的人工阻尼g、圆频率ω和速度U；

4)计算完所有的折合频率K下的数值后，将圆频率转换为频率，即可绘制相应的速度-人工阻尼U-g图和速度-频率U-f图；

5)从速度-人工阻尼U-g图中找到曲线与x轴相交的点，即g为0的点，此处对应的速度U即为颤振速度U_F，在频率U-f图中对应于颤振速度U_F的点即为颤振频率f_F。

图7、8为采用本文方法计算出的Goland机翼的速度-人工阻尼图和速度-频率图。动力学模型是考虑弯曲和扭转的双自由度模型所以有两条曲线，弯曲是指弯曲运动，扭转是指扭转运动。如图所示，随着无量纲速度的增加，扭转分支对应的人工阻尼g先变为0，所以Goland机翼的颤振速度U_F＝55.65m/s＝182ft/s，对应图8中颤振频率f_F＝3.68Hz。

同样，本方法对其他非均直机翼、舵面系统建立颤振模型，只需修改简化模型和模型中具体参数即可。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和变型，这些改进和变型也应视为本发明的保护范围。

Claims (4)

1.一种基于多体系统传递矩阵法求解三维机翼线性颤振速度的方法，其特征是，包括以下过程：

基于多体系统传递矩阵法推导机翼的弯扭耦合梁传递矩阵；

建立机翼的总传递方程，并求解得到机翼的固有频率和振型；

根据Theodorsen非定常气动力理论建立机翼振动位移与气动力的关系；

建立机翼的体动力学方程，并将其转换到频域获得机翼的颤振频域模型；

对机翼颤振频域模型进行频域求解，得到机翼的颤振速度；

所述基于多体系统传递矩阵法推导机翼的弯扭耦合梁传递矩阵包括：

将机翼简化为n段弯扭耦合梁，建立弯扭耦合梁的振动微分方程：

式中，m为线质量，I_α为单位长度转动惯量，EI为弯曲刚度，GJ为扭转刚度，y为机翼弯曲方向位移，θ_x为扭转角度，t为时间，b为机翼弦长的一半，x_α为质心到弹性轴的距离，且当质心在Z正方向，x_α为正；

将物理坐标转化为模态坐标，令

y(x,t)＝Y(x)sinωt,θ_x(x,t)＝Θ_x(x)sinωt (2)

式中，Y(x)代表机翼弯曲方向的位移的模态坐标，Θ_x(x)代表机翼弯曲方向的角度的模态坐标，ω为圆频率，Y(x)、Θ_x(x)为y(x,t)、θ_x(x,t)通过模态坐标转换得到；

将(2)带入(1)中，(1)式变为：

消去式(3)中的Y(x)或者Θ_x(x)，得到

式中，

W＝Y或者Θ (5)

引入无量纲长度：

ξ＝x/L (6)

式(6)中L为划分后单元弯扭耦合梁长度；

式(4)可改写成无量纲形式：

(D⁶+aD⁴-bD²-abc)W＝0 (7)

其中，

六阶微分方程(7)的通解可以表示为：

W(ξ)＝C₁coshαξ+C₂sinhαξ+C₃cosβξ+C₄sinβξ+C₅cosγξ+C₆sinγξ (9)

式中C₁,......,C₆是常数，并且

q＝b+a²/3

式(9)中的W(ξ)代表弯曲位移Y和扭转角度Θ_x在不同常数下的解；因此，

Y(ξ)＝A₁coshαξ+A₂sinhαξ+A₃cosβξ+A₄sinβξ+A₅cosγξ+A₆sinγξ (11)

Θ_x(ξ)＝B₁coshαξ+B₂sinhαξ+B₃cosβξ+B₄sinβξ+B₅cosγξ+B₆sinγξ (12)

式中A₁,......,A₆和B₁,......,B₆是两组不同的常数；

将式(11)和(12)代入式(3)可以确定常数有如下规律：

式中，

从式(11)、(12)可以得到弯曲角度Θ_z(ξ)，弯矩M_z(ξ)，剪切力Q_y(ξ)和扭矩M_x(ξ)的表达式：

为了推导弯扭耦合梁的传递矩阵，由式(11)、(12)、(15)、(16)、(17)、(18)整理可得

即Z(ξ)＝B(ξ)·a,a＝[A₁,A₂,A₃,A₄,A₅,A₆]^T，其中Z(ξ)为状态矢量，即方程(19)等式左边项，B(ξ)代表方程(19)等式右边的第一个矩阵；令ξ＝0，且根据方程(19)和多体系统传递矩阵法可将方程(19)写成Z_I＝B(0)·a形式；令ξ＝1得到Z_O＝B(1)·a＝B(1)B(0)^-1Z_I＝U_iZ_I；其中Z_O代表模态坐标下元件输出端的状态矢量，Z_I代表模态坐标下元件输入端的状态矢量，U_i为第i段的弯扭耦合梁的传递矩阵；

所以第i段的弯扭耦合梁的传递矩阵U_i为：

U_i＝B(1)·B^-1(0) (20)

其中，

建立机翼的体动力学方程，并将其转换到频域获得机翼的颤振频域模型包括：

根据多体系统传递矩阵法，建立机翼的体动力学方程为：

Mv_tt+Cv_t+Kv＝f (29)

其中，M、C、K分别为质量矩阵、阻尼矩阵、刚度矩阵；v_t,v_tt分别为v对时间t的一阶导数和二阶导数；

结合公式27并且将v＝v_se^iωt代入式(29)中，则有：

-ω²Mv_s+iωCv_s+Kv_s＝Bv_s (30)

由于体动力学方程式29右边的外力项非解耦，把振型写为V＝[V¹,V²]，其中V¹为机翼的第1阶振型，V²为机翼的第2阶振型；

令式(30)中的v_s＝Vq，q为机翼振动位移的广义坐标，得到：

-ω²MVq+iωCVq+KVq＝BVq (31)

由于增广特征矢量的正交性(左乘V^T)，因此有：

-ω²V^TMVq+iωV^TCVq+V^TKVq＝V^TBVq (32)

式中，

计算中忽略结构的阻尼，因此

其中ω₁、ω₂分别为机翼结构的第一阶圆频率和第二阶圆频率。

2.根据权利要求1所述的一种基于多体系统传递矩阵法求解三维机翼线性颤振速度的方法，其特征是，所述建立机翼的总传递方程包括：

根据各段弯扭耦合梁传递矩阵，建立机翼总传递方程为：

Z_o＝U_allZ_i (21)

式中U_all＝U_n···U₃U₂U₁，U₁,U₂,U₃···U_n分别代表机翼各段弯扭耦合梁的传递矩阵，Z_i,Z_o分别代表系统模态坐标下输入端、输出端的状态矢量；输入端Z_i＝[0,0,M_z,Q_y,0,M_x]^T，输出端Z_o＝[Y,Θ_z,0,0,Θ_x,0]^T。

3.根据权利要求1所述的一种基于多体系统传递矩阵法求解三维机翼线性颤振速度的方法，其特征是，所述根据Theodorsen非定常气动力理论建立机翼振动位移与气动力的关系包括：

根据机翼坐标系，使用Theodorsen公式计算机翼气动力：

式(22)、(23)中，L表示升力，T_α表示俯仰力矩，U为来流速度，ρ为空气密度，b为半弦长，a为刚心距离弦线中点的距离与半弦长的比值，并在中点后为正，C(K)为Theodorsen函数，其值与折合频率K＝ωb/U有关，ω为圆频率，y表示沉浮位移，

为速度，

为加速度；θ_x表示俯仰位移，

表示俯仰速度，

表示俯仰速度；

假设Goland机翼各方向的运动为简谐运动，即：

y＝y_se^iωt；θ＝θ_se^iωt (24)

式中，y、θ分别表示物理坐标系下的沉浮和俯仰位移；y_s、θ_s分别表示模态坐标下的沉浮和俯仰位移，ω表示圆频率，t表示时间；

将该式代入到式(22)、(23)中，可以整理得：

令f＝[L₁,T_a1,L₂,T_a2,...L_n,T_an]^T，n表示将机翼分为n段，L_n、T_αn表示机翼第n段上面的升力和俯仰力矩；令v＝v_se^iωt，利用式(25)、(26)可以建立以下位移与气动力的关系：

f＝Bv_se^iωt (27)

式(27)中，B为矩阵，具体数值如下式所示：

B＝[B₁,B₂,...B_n]^T， v_s＝[v_s1,v_s2,...v_sn]^T (28)

其中，

k_a＝-1+i·(2Π/π)·C(K)/K

k_b＝-0.5+i·(1+2Π/π·C(K))/K+(2Π/π)·C(K)/K²

m_a＝0.5

b_n、l_n、a_n分别表示机翼第n段的半弦长、长度、刚心距离弦线中点的距离与半弦长的比值；B_n为机翼第n段的气动力矩阵；v_sn表示机翼第n段模态坐标下的沉浮和俯仰位移。

4.根据权利要求3所述的一种基于多体系统传递矩阵法求解三维机翼线性颤振速度的方法，其特征是，所述对机翼颤振频域模型进行频域求解包括：

采用U-g法对机翼颤振频域模型进行频域求解；

所述采用U-g法对机翼颤振频域模型进行频域求解包括：

机翼颤振速度的求解过程转化为求解该复特征值的问题，其特征值为：

式中，λ为特征值；λ_Re、λ_Im分别为特征值的实部、虚部；

由此可以得到：

式中，g、ω和U分别为人工阻尼、圆频率和来流速度；

最后，可以通过以下过程来实现颤振分析：

1)给定密度ρ，预设一折合频率K；

2)在指定折合频率K下计算式(36)，并得到该折合频率下的人工阻尼g、圆频率ω和来流速度U；

3)步长减小折合频率K，并重复第2步，计算各折合频率下的人工阻尼g、圆频率ω和来流速度U；

4)计算完所有的折合频率K下的数值后，将圆频率转换为频率，绘制相应的速度-人工阻尼U-g图和速度-频率U-f图；

5)从速度-人工阻尼U-g图中找到曲线与x轴相交的点，此处对应的速度U即为颤振速度U_F，在频率U-f图中对应于颤振速度U_F的点即为颤振频率f_F。

Images