基于传递函数的细长体颤振速度确定方法.pdf

本发明涉及一种基于传递函数的细长体颤振速度确定方法，属于飞行器气动弹性技术领域。所述方法利用梁的横向弯曲振动微分方程和细长体准定常气动力模型，得出细长体颤振微分方程，然后对所述细长体离散为细长体单元，在单元内对细长体单元颤振微分方程进行Fourier变换，再运用传递函数法求解细长体颤振速度；本发明利用梁的横向弯曲振动微分方程来准确描述细长体的横向振动，而没有采用梁横向振动的低阶模态来近似描述细长体振动，所以计算的细长体颤振速度更加准确，本发明求解细长体颤振速度的方法与现有技术相比更加简捷。

1.一种基于传递函数的细长体颤振速度确定方法，其特征在于所述方法利用梁的横向
弯曲振动微分方程和细长体准定常气动力模型，得出细长体颤振微分方程，然后对所述细
长体颤振微分方程进行Fourier变换，再运用传递函数法求解细长体颤振速度；
一、所述方法需要依据的计算公式：
a.利用梁的横向弯曲振动微分方程和细长体准定常气动力模型建立细长体颤振微分
方程：
本发明利用梁的横向弯曲振动微分方程式(1)来描述细长体的弯曲振动：
$\frac{{\∂}^2\left(EI\frac{{\∂}^{2}h}{\∂y^2}\right)}{\∂y^2}+m\frac{{\∂}^{2}h}{\∂t^2}+\mathrm{\Δ}p=0---\left(1\right)$
式中，h为细长体弯曲振动位移，单位为米，EI为细长体抗弯刚度，单位牛顿·米2，m为细
长体单位长度质量，单位为千克，△p为细长体单位长度的气动力，单位为牛顿/米，y为细体
长轴向坐标，单位为米，t为时间，单位为秒，为细长体弯曲振动位移h对机轴向坐标y的
二阶偏导数，为细长体弯曲振动位移h对时间t的二阶偏导数；
根据细长体准定常气动力模型，得到细长体单位长度的气动力△p为下式(2)：
$\mathrm{\Δ}p=-\{\mathrm{\ρ}A(\frac{{\∂}^{2}h}{\∂t^2}+V\frac{{\∂}^{2}h}{\∂t\∂y})+\mathrm{\ρ}V\frac{\∂A}{\∂y}(\frac{\∂h}{\∂t}+V\frac{\∂h}{\∂y})+\mathrm{\ρ}AV(\frac{{\∂}^{2}h}{\∂t\∂y}+V\frac{{\∂}^{2}h}{\∂y^2})\}---\left(2\right)$
式中，V为空速，单位为米/秒，A为细长体横截面积，单位为平方米，ρ为空气密度，单位
为千克/立方米；
将(2)式代入(1)式得到细长体颤振微分方程式(3)：
$\frac{{\∂}^2\left(EI\frac{{\∂}^{2}h}{{\mathrm{ay}}^2}\right)}{{\mathrm{ay}}^2}+m\frac{{\∂}^{2}h}{\∂t^2}-\mathrm{\ρ}A(\frac{{\∂}^{2}h}{\∂t^2}+V\frac{{\∂}^{2}h}{\∂t\∂y})-\mathrm{\ρ}V\frac{\∂A}{ay}(\frac{\∂h}{\∂t}+V\frac{\∂h}{ay})-\mathrm{\ρ}AV(\frac{{\∂}^{2}h}{\∂t\∂y}+V\frac{{\∂}^{2}h}{{\mathrm{ay}}^2})=0---\left(3\right)$
b.沿细长体轴向划分n个细长体单元，建立细长体单元的颤振微分方程：
一般地，细长体的物理参数EI、m、A等沿其轴向是变化的，为了方便求解，按照有限元
法，沿细长体轴向划分n个单元，为了使细长体单元的形式统一颤振微分方程，令：
$\mathrm{\ξ}=\frac{y-y_j}{l_j}---\left(4\right)$
式中，yj为第j个细长体单元前节点距离细长体前端点的距离，lj为第j个细长体单元的
长度，ξ为细长体单元内轴向无量纲坐标，由(4)式可知ξ∈(0,1)；
当划分单元足够多时，假设EI、m，A在每个细长体单元内近似线性变化是可理的，则在
第j个细长体单元内有如下线性表达式：
$\left\{\begin{array}{c}EI\left(\mathrm{\ξ}\right)=EI\left(0\right)+\frac{EI\left(1\right)-EI\left(0\right)}{l_j}\mathrm{\ξ}\\ m\left(\mathrm{\ξ}\right)=m\left(0\right)+\frac{m\left(1\right)-m\left(0\right)}{l_j}\mathrm{\ξ}\\ A\left(\mathrm{\ξ}\right)=A\left(0\right)+\frac{A\left(1\right)-A\left(0\right)}{l_j}\mathrm{\ξ}\end{array}\right.---\left(5\right)$
式中，EI(0)、m(0)、A(0)为EI、m、A在第j个细长体单元前节点处的值，EI(1)、m(1)、A(1)
为EI、m、A在第j个细长体单元后节点处的值，EI(ξ)、m(ξ)、A(ξ)为EI、m、A在第j个细长体单
元内ξ处的值；
将(4)式与(5)式代入(3)式，可得到细长体单元的颤振微分方程：
$\begin{array}{c}\frac{EI\left(\mathrm{\ξ}\right)}{l_j^4}\frac{{\∂}^{4}h}{\∂{\mathrm{\ξ}}^4}+\frac{2\left(EI\right(1)-EI(0\left)\right)}{l_j^4}\frac{{\∂}^{3}h}{\∂{\mathrm{\ξ}}^3}-\frac{\mathrm{\ρ}A\left(\mathrm{\ξ}\right)V^2}{l_j^2}\frac{{\∂}^{2}h}{\∂{\mathrm{\ξ}}^2}\\ -\frac{2\mathrm{\ρ}A\left(\mathrm{\ξ}\right)V}{l_j}\frac{{\∂}^{2}h}{\∂t\∂\mathrm{\ξ}}-{\mathrm{\ρV}}^2\frac{A\left(1\right)-A\left(0\right)}{l_j^2}\frac{\∂h}{\∂\mathrm{\ξ}}+\[m\left(\mathrm{\ξ}\right)-\mathrm{\ρ}A\left(\mathrm{\ξ}\right)\]\frac{{\∂}^{2}h}{\∂{\mathrm{\ξ}}^2}-\mathrm{\ρ}V\frac{A\left(1\right)-A\left(0\right)}{l_j}\frac{\∂h}{\∂t}=0\end{array}---\left(6\right)$
c.利用传递函数法求解细长体颤振速度
①对(6)式中有关时间t的变量作Fourier变换，可得到下式(7)：
$\begin{array}{c}\frac{EI\left(\mathrm{\ξ}\right)}{l_j^4}\frac{{\∂}^{4}h}{\∂{\mathrm{\ξ}}^4}+\frac{2\left(EI\right(1)-EI(0\left)\right)}{l_j^4}\frac{{\∂}^{3}h}{\∂{\mathrm{\ξ}}^3}-\frac{\mathrm{\ρ}A\left(\mathrm{\ξ}\right)V^2}{l_j^2}\frac{{\∂}^{2}h}{\∂{\mathrm{\ξ}}^2}\\ -\[i\mathrm{\ω}\frac{2\mathrm{\ρ}A\left(\mathrm{\ξ}\right)V}{l_j}+V^2\frac{\mathrm{\σ}A\left(1\right)-\mathrm{\σ}A\left(0\right)}{l_j^2}\]\frac{\∂h}{\∂\mathrm{\ξ}}+\[{\left(i\mathrm{\ω}\right)}^2\[m\left(\mathrm{\ξ}\right)-\mathrm{\ρ}A\left(\mathrm{\ξ}\right)\]-\left(i\mathrm{\ω}\right)V\frac{\mathrm{\ρ}A\left(1\right)-\mathrm{\ρ}A\left(0\right)}{l_j}\]h=0\end{array}---\left(7\right)$
式中，虚数ω为角频率，单位为弧度/秒；
②将细长体单元的颤振微分方程改写为状态空间方程形式；
为了便于应用传递函数理论，定义一个状态变量的向量ηe(ξ,ω)如下式(8)：
${\mathrm{\η}}_e(\mathrm{\ξ},\mathrm{\ω})={\left[\begin{array}{cccc}h& \frac{\∂h}{\∂y}& \frac{{\∂}^{2}h}{\∂y^2}& \frac{{\∂}^{3}h}{\∂y^3}\end{array}\right]}^T---\left(8\right)$
式中，T表示向量转置；
基于状态变量向量ηe(ξ,ω)，可将(7)式写成状态方程的形式如下式(9)：
$\frac{\∂{\mathrm{\η}}_e(\mathrm{\ξ},\mathrm{\ω})}{\∂\mathrm{\ξ}}=F_e(\mathrm{\ξ},\mathrm{\ω},V){\mathrm{\η}}_e(\mathrm{\ξ},\mathrm{\ω})+g_e(\mathrm{\ξ},\mathrm{\ω})---\left(9\right)$
式中，Fe(ξ,ω,V)、ge(ξ,ω)的表达式可根据(7)式得出，具体如下
$\left\{\begin{array}{c}{F}_e(\mathrm{\ξ},\mathrm{\ω},V)=\left[\begin{array}{cccc}0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\\ A_4& A_3& A_2& A_1\end{array}\right]\\ g_e(y,\mathrm{\ω})=0\end{array}\right.---\left(10\right)$
式中，
$\left\{\begin{array}{c}{A}_1=-2\frac{EI\left(1\right)-EI\left(0\right)}{EI\left(\mathrm{\ξ}\right)}\\ A_2=\frac{\mathrm{\ρ}A\left(\mathrm{\ξ}\right)V^2}{EI\left(\mathrm{\ξ}\right)}{l}_j^2\\ A_3=i\mathrm{\ω}\frac{2\mathrm{\ρ}A\left(\mathrm{\ξ}\right)V}{EI\left(\mathrm{\ξ}\right)}{l}_j^3+V^2\frac{\mathrm{\ρ}A\left(1\right)-\mathrm{\ρ}A\left(0\right)}{EI\left(\mathrm{\ξ}\right)}{l}_j^2\\ A_4=-\frac{{\left(i\mathrm{\ω}\right)}^2\[m\left(\mathrm{\ξ}\right)-\mathrm{\ρ}A\left(\mathrm{\ξ}\right)\]}{EI\left(\mathrm{\ξ}\right)}{l}_j^4+\left(i\mathrm{\ω}\right)V\frac{\mathrm{\ρ}A\left(1\right)-\mathrm{\ρ}A\left(0\right)}{EI\left(\mathrm{\ξ}\right)}{l}_j^3\end{array}\right.---\left(11\right)$
根据传递函数理论，方程(9)式的标准解为：
ηe(ξ,ω)＝He(ξ,ω,V)γe(ω)+de(ξ,ω) (12)
式中，γe(ω)为由边界条件给定的位移或力组成的向量，可根据传递函数方法理论得
到表达式如下：
${\mathrm{\γ}}_e\left(\mathrm{\ω}\right)=\left\{\begin{array}{c}h(0,\mathrm{\ω})\\ \frac{\∂h(0,\mathrm{\ω})}{\∂\mathrm{\ξ}}\\ h(1,\mathrm{\ω})\\ \frac{\∂h(1,\mathrm{\ω})}{\∂\mathrm{\ξ}}\end{array}\right\}---\left(13\right)$
式中，He(ξ,ω,V)、fe(ξ,ω)的表达式也可根据传递函数方法理论得到，如下：
$\left\{\begin{array}{c}{H}_e(\mathrm{\ξ},\mathrm{\ω},V)={\mathrm{\Φ}}_F(\mathrm{\ξ},0,\mathrm{\ω},V){\[M_b{\mathrm{\Φ}}_F(0,0,\mathrm{\ω},V)+N_b{\mathrm{\Φ}}_F(1,0,\mathrm{\ω},V)\]}^{-1}\\ d_e(\mathrm{\ξ},\mathrm{\ω})={\∫}_0^1(\mathrm{\ξ},\mathrm{\ζ},\mathrm{\ω},V)g_e(\mathrm{\ζ},\mathrm{\ω})d\mathrm{\ζ}\\ G(\mathrm{\ξ},\mathrm{\ζ},\mathrm{\ω},V)=\left\{\begin{array}{cc}{H}_e(\mathrm{\ξ},\mathrm{\ω},V)M_b{\mathrm{\Φ}}_F(0,\mathrm{\ζ},\mathrm{\ω},V)& \mathrm{\ζ}\<\mathrm{\ξ}\\ -H_e(\mathrm{\ξ},\mathrm{\ω},V)N_b{\mathrm{\Φ}}_F(1,\mathrm{\ζ},\mathrm{\ω},V)& \mathrm{\ζ}\>\mathrm{\ξ}\end{array}\right.\\ {\mathrm{\Φ}}_F(\mathrm{\ξ},\mathrm{\ζ},\mathrm{\ω},V)=e^{F_e(\mathrm{\ξ},\mathrm{\ω},V)\mathrm{\ζ}}\end{array}\right.---\left(14\right)$
式中，变量ζ∈(0,1)，Mb为细长体单元前端边界条件选择矩阵，Nb为细长体单元后端边
界条件选择矩阵，可根据传递函数方法理论得到，如下：
$M_b=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right],N_b=\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 1& 0& 0& 0\\ 0& 1& 0& 0\end{array}\right]---\left(15\right)$
由于(10)式中给出ge(y,ω)＝0，由(14)式可知de(ξ,ω)＝0，从而(12)式可简化为：
ηe(ξ,ω)＝He(ξ,ω,V)γe(ω) (16)
③根据细长体内力关系建立单元平衡方程；
由材料力学梁的弯矩和剪力公式可知，细长体截面上弯矩和剪力的表达式可写为：
$\left\{\begin{array}{c}M(y,\mathrm{\ω})=EI\frac{{\∂}^{2}h}{\∂y^2}\\ F(y,\mathrm{\ω})=\frac{\∂\left(EI\frac{{\∂}^{2}h}{\∂y^2}\right)}{\∂y}\end{array}\right.---\left(17\right)$
在细长体单元内，并利用EI在单元内的线性变化假设，(17)式可写为：
${\mathrm{\σ}}_e(\mathrm{\ξ},\mathrm{\ω})=\left\{\begin{array}{c}M(\mathrm{\ξ},\mathrm{\ω})\\ F(\mathrm{\ξ},\mathrm{\ω})\end{array}\right\}=\left\{\begin{array}{c}\frac{EI\left(\mathrm{\ξ}\right)}{l_i^2}\frac{{\∂}^{2}h}{\∂{\mathrm{\ξ}}^2}\\ \frac{EI\left(\mathrm{\ξ}\right)}{l_i^3}\frac{{\∂}^{3}h}{\∂{\mathrm{\ξ}}^3}+\frac{EI\left(1\right)-EI\left(0\right)}{l_i^3}\frac{{\∂}^{2}h}{\∂{\mathrm{\ξ}}^2}\end{array}\right\}---\left(18\right)$
引入状态变量向量(8)式，从而(18)可写成如下矩阵形式：
σe(ξ,ω)＝Eeη(ξ)ηe(ξ,ω) (19)
式中，Eeη(ξ)可根据(18)式得到其表达式如下：
$E_{e\mathrm{\η}}\left(\mathrm{\ξ}\right)=\left[\begin{array}{cccc}0& 0& \frac{EI\left(\mathrm{\ξ}\right)}{l_i^2}& 0\\ 0& 0& \frac{EI\left(1\right)-EI\left(0\right)}{l_i^3}& \frac{EI\left(\mathrm{\ξ}\right)}{l_i^3}\end{array}\right]---\left(20\right)$
将(16)式代入(19)式，可得到
σe(ξ,ω)＝Eeη(ξ)He(ξ,ω,V)γe(ω) (21)
式中，γe(ω)可视为细长体单元两节点处的广义位移，如果ξ取0和1，则可得到细长体
单元两节点的广义内力与广义位移之间的关系，从而建立单元平衡方程，即：
令ξ＝0，ξ＝1，则有
${\mathrm{\σ}}_e\left(\mathrm{\ω}\right)=\left\{\begin{array}{c}{\mathrm{\σ}}_e(0,\mathrm{\ω})\\ {\mathrm{\σ}}_e(1,\mathrm{\ω})\end{array}\right\}=\left[\begin{array}{c}{E}_{e\mathrm{\η}}\left(0\right)H_e(0,\mathrm{\ω},V)\\ E_{e\mathrm{\η}}\left(1\right)H_e(1,\mathrm{\ω},V)\end{array}\right]{\mathrm{\γ}}_e\left(\mathrm{\ω}\right)---\left(22\right)$
上式即为表征细长体单元节点广义内力与广义位移关系的单元平衡方程，从而
可视为单元的广义刚度矩阵，可记为：
$K_e(\mathrm{\ω},V)=\left[\begin{array}{c}{E}_{e\mathrm{\η}}\left(0\right)H_e(0,\mathrm{\ω},V)\\ E_{e\mathrm{\η}}\left(1\right)H_e(1,\mathrm{\ω},V)\end{array}\right]---\left(23\right)$
④将细长体单元平衡方程组装形成细长体整体平衡方程；
按照有限元法进行单元组装，可得到细长体整体平衡方程为
K(ω,V)γ(ω)＝F(ω) (24)
式中，K(ω,V)可视为整体刚度矩阵，γ(ω)可视为整体节点位移向量，F(ω)为各单元
节点内力拼装成的向量。
⑤确定颤振速度；
对于细长体，颤振时不需虑重力，本专利中将细长体所受气动力的影响囊括在刚度矩
阵K(ω,V)，除此以外细长体不受其它力的作用，从而外载荷为零，根据单元节点内力与外
载荷平衡，可得出
F(ω)＝0 (25)
当细长体颤振时，γ(ω)应有非零解，此时整体平衡方程必须满足条件为
det[K(ω,V)]＝0 (26)
由于K(ω,V)为复矩阵，其行列式值等于零的必要条件为矩阵行列式值的实部与虚部
均为零，即
$\left\{\begin{array}{c}\mathrm{Re}\{\det \[K(\mathrm{\ω},V)\]\}=0\\ \mathrm{Im}\{\det \[K(\mathrm{\ω},V)\]\}=0\end{array}\right.---\left(27\right)$
矩阵K(ω,V)中有空速V和圆频率ω两个变量，而(27)式恰好有两个方程，可以定解。求
解上述方程组时，可能会得到满足方程组多个解，即存在多组(V,ω)能满足方程组(27)式。
根据细长体颤振时在某一空速时由稳定转变为不稳定，因而，空速V最小的一组解(V,ω)应
为细长体的颤振速度和相应的颤振频率。将得到满足方程组的空速V和圆频率ω，分别记为
Vcz和ωcz。其中，Vcz即为细长体的颤振速度，ωcz即为细长体的颤振圆频率。
二、所述方法的具体步骤如下：
步骤(一)：将细长体划分为n个单元，测量细长体各单元节点的下列物理参数：
细长体各单元轴向长度[l1 l2 … lj … ln]，单位为米；
细长体单元节点处的横截面积[A1 A2 … Aj … An+1]，单位为米2
细长体单位节点处单位长度质量[m1 m2 … mj … mn+1]，单位为千克/米；
细长体单位节点处抗弯刚度[EI1 EI2 … EIj … EIn+1]，单位为牛顿·米2；
空气密度ρ，单位为千克/米3；
步骤(二)：确定细长体飞行的空速V和圆频率ω的大致范围，假设为
步骤(三)：在范围内划分合适步长△V和△ω，并进行离散，空速V和圆频率
ω的取值为：
步骤(四)：取空速V＝V0，圆频率ω依次取ω0+j△ω(j＝0,1,2,3…)，将空速V和圆频率
ω的取值与步骤(一)中的细长体物理参数代入(11)式，计算各单元系数A1、A2、A3、A4的值；
步骤(五)：将系数A1、A2、A3、A4的值代入(10)式，得到细长体各单元矩阵Fe(ξ,ω,V)的
值；
步骤(六)：将矩阵Fe(ξ,ω,V)的值代入(14)式，结合(15)式，计算He(ξ,ω,V)的值；
步骤(七)：将He(ξ,ω,V)的值代入(23)式，结合(20)式，计算单元刚度矩阵Ke(ω,V)的
值；
步骤(八)：利用Ke(ω,V)组装整体刚度矩阵K(ω,V)；
步骤(九)：计算Re{det[K(ω,V)]}和Im{det[K(ω,V)]}的值；
步骤(十)：再依次取空速V＝V0+j△V(j＝1,2,3…)，重复步骤(四)至步骤(十)，计算Re
{det[K(ω,V)]}和Im{det[K(ω,V)]}的值；
步骤(十一)：确定满足(27)式的空速V的值，即可确定细长体的颤振速度Vcz。

基于传递函数的细长体颤振速度确定方法

技术领域

本发明涉及一种基于传递函数的细长体颤振速度确定方法，属于飞行器气动弹性
技术领域。

背景技术

对于飞机机身、导弹弹体、机翼外挂火箭等，当长细比较大时(一般大于10左右)，
称为细长体。细长体颤振是指发生在飞行器飞行中的动不稳定性，此时的飞行速度称为颤
振速度。在飞行达到颤振速度时所发生的自激振动，大多数都会造成灾难性的后果。颤振分
析就是求出颤振发生的条件，亦即求出颤振速度，并寻求在飞行器飞行速度范围内避免颤
振发生的措施，进而寻找提高颤振临界速度的方法。

进行细长体颤振分析和计算时，一方面需要知道细长体结构的动力学特性，另一
方面需要知道细长体结构表面非定常流动的空气动力特性。目前，工程中为了避免分析计
算过于复杂，在细长体结构动力学计算时尽量减小细长体振动自由度，通常忽略高阶模态，
仅利用细长体振动的低阶模态来近似表示细长体的弯曲振动位移。然后，再结合细长体准
定常气动理论，求解细长体颤振速度。

实际上，细长体是无限多自由度的弹性体，在颤振分析中也没有对自由度的限制。
因此，上述计算方法的缺点是：为了避免过于复杂的计算而减少了细长体振动的自由度，对
细长体实际振动作出了近似，使计算结果的准确性降低了。

发明内容

本发明所要解决的技术问题是提供一种计算方法简捷、计算结果较准确的基于传
递函数的细长体颤振速度确定方法。

本发明解决其技术问题所需采用的技术方案：

一种基于传递函数法的细长体颤振速度确定方法，所述方法利用梁的横向弯曲振
动微分方程和细长体准定常气动力模型，得出细长体颤振微分方程，然后对所述细长体离
散为细长体单元，在单元内对细长体单元颤振微分方程进行Fourier变换，再运用传递函数
法求解细长体颤振速度；

一、所述方法需要依据的计算公式：

a.利用梁的横向弯曲振动微分方程和细长体准定常气动力模型建立细长体颤振
微分方程：

本发明利用梁的横向弯曲振动微分方程式(1)来描述细长体的弯曲振动：

$\frac{{\∂}^2\left(EI\frac{{\∂}^{2}h}{\∂y^2}\right)}{\∂y^2}+m\frac{{\∂}^{2}h}{\∂t^2}+\mathrm{\Δ}p=0---\left(1\right)$

式中，h为细长体弯曲振动位移，单位为米，EI为细长体抗弯刚度，单位牛顿·米2，
m为细长体单位长度质量，单位为千克，△p为细长体单位长度的气动力，单位为牛顿/米，y
为细体长轴向坐标，单位为米，t为时间，单位为秒，为细长体弯曲振动位移h对机轴向坐
标y的二阶偏导数，为细长体弯曲振动位移h对时间t的二阶偏导数；

根据细长体准定常气动力模型，得到细长体单位长度的气动力△p为下式(2)：

$\mathrm{\Δ}p=-\{\mathrm{\ρ}A(\frac{{\∂}^{2}h}{\∂t^2}+V\frac{{\∂}^{2}h}{\∂t\∂y})+\mathrm{\ρ}V\frac{\∂A}{\∂y}(\frac{\∂h}{\∂t}+V\frac{\∂h}{\∂y})+\mathrm{\ρ}AV(\frac{{\∂}^{2}h}{\∂t\∂y}+V\frac{{\∂}^{2}h}{\∂y^2})\}---\left(2\right)$

式中，V为空速，单位为米/秒，A为细长体横截面积，单位为平方米，ρ为空气密度，
单位为千克/立方米；

将(2)式代入(1)式得到细长体颤振微分方程式(3)：

$\frac{{\∂}^2\left(EI\frac{{\∂}^{2}h}{{\mathrm{ay}}^2}\right)}{\∂y^2}+m\frac{{\∂}^{2}h}{\∂t^2}-\mathrm{\ρ}A(\frac{{\∂}^{2}h}{\∂t^2}+V\frac{{\∂}^{2}h}{\∂t\∂y})-\mathrm{\ρ}V\frac{\∂A}{\∂y}(\frac{\∂h}{\∂t}+V\frac{\∂h}{\∂y})-\mathrm{\ρ}AV(\frac{{\∂}^{2}h}{\∂t\∂y}+V\frac{{\∂}^{2}h}{\∂y^2})=0---\left(3\right)$

b.沿细长体轴向划分n个细长体单元，建立细长体单元的颤振微分方程：

一般地，细长体的物理参数EI、m、A等沿其轴向是变化的，为了方便求解，按照有限
元法，沿细长体轴向划分n个单元，为了使细长体单元的形式统一颤振微分方程，令：

$\mathrm{\ξ}=\frac{y-y_j}{l_j}---\left(4\right)$

式中，yj为第j个细长体单元前节点距离细长体前端点的距离，lj为第j个细长体单
元的长度，ξ为细长体单元内轴向无量纲坐标，由(4)式可知ξ∈(0,1)；

当划分单元足够多时，假设EI、m，A在每个细长体单元内近似线性变化是可理的，
则在第j个细长体单元内有如下线性表达式：

$\left\{\begin{array}{c}EI\left(\mathrm{\ξ}\right)=EI\left(0\right)+\frac{EI\left(1\right)-EI\left(0\right)}{l_j}\mathrm{\ξ}\\ m\left(\mathrm{\ξ}\right)=m\left(0\right)+\frac{m\left(1\right)-m\left(0\right)}{l_j}\mathrm{\ξ}\\ A\left(\mathrm{\ξ}\right)=A\left(0\right)+\frac{A\left(1\right)-A\left(0\right)}{l_j}\mathrm{\ξ}\end{array}\right.---\left(5\right)$

式中，EI(0)、m(0)、A(0)为EI、m、A在第j个细长体单元前节点处的值，EI(1)、m(1)、
A(1)为EI、m、A在第j个细长体单元后节点处的值，EI(ξ)、m(ξ)、A(ξ)为EI、m、A在第j个细长
体单元内ξ处的值；

将(4)式与(5)式代入(3)式，可得到细长体单元的颤振微分方程：

$\begin{array}{c}\frac{EI\left(\mathrm{\ξ}\right)}{l_j^4}\frac{{\∂}^{4}h}{\∂{\mathrm{\ξ}}^4}+\frac{2\left(EI\right(1)-EI(0\left)\right)}{l_j^4}\frac{{\∂}^{3}h}{\∂{\mathrm{\ξ}}^3}-\frac{\mathrm{\ρ}A\left(\mathrm{\ξ}\right)V^2}{l_j^2}\frac{{\∂}^{2}h}{\∂{\mathrm{\ξ}}^2}\\ -\frac{2\mathrm{\ρ}A\left(\mathrm{\ξ}\right)V}{l_j}\frac{{\∂}^{2}h}{\∂t\∂\mathrm{\ξ}}-{\mathrm{\ρV}}^2\frac{A\left(1\right)-A\left(0\right)}{l_j^2}\frac{\∂h}{\∂\mathrm{\ξ}}+\[m\left(\mathrm{\ξ}\right)-\mathrm{\ρ}A\left(\mathrm{\ξ}\right)\]\frac{{\∂}^{2}h}{\∂t^2}-\mathrm{\ρ}V\frac{A\left(1\right)-A\left(0\right)}{l_j}\frac{\∂h}{\∂t}=0\end{array}---\left(6\right)$

c.利用传递函数法求解细长体颤振速度

①对(6)式中有关时间t的变量作Fourier变换，可得到下式(7)：

$\begin{array}{c}\frac{EI\left(\mathrm{\ξ}\right)}{l_j^4}\frac{{\∂}^{4}h}{\∂{\mathrm{\ξ}}^4}+\frac{2\left(EI\right(1)-EI(0\left)\right)}{l_j^4}\frac{{\∂}^{2}h}{\∂{\mathrm{\ξ}}^3}-\frac{\mathrm{\ρ}A\left(\mathrm{\ξ}\right)V^2}{l_j^2}\frac{{\∂}^{2}h}{\∂{\mathrm{\ξ}}^2}\\ -\[i\mathrm{\ω}\frac{2\mathrm{\ρ}A\left(\mathrm{\ξ}\right)V}{l_j}+V^2\frac{\mathrm{\ρ}A\left(1\right)-\mathrm{\ρ}A\left(0\right)}{l_j^2}\]\frac{\∂h}{\∂\mathrm{\ξ}}+\[{\left(i\mathrm{\ω}\right)}^2\[m\left(\mathrm{\ξ}\right)-\mathrm{\ρ}A\left(\mathrm{\ξ}\right)\]-\left(i\mathrm{\ω}\right)V\frac{\mathrm{\ρ}A\left(1\right)-\mathrm{\ρ}A\left(0\right)}{l_j}\]h=0\end{array}---\left(7\right)$

式中，虚数ω为角频率，单位为弧度/秒；

②将细长体单元的颤振微分方程改写为状态空间方程形式；

为了便于应用传递函数理论，定义一个状态变量的向量ηe(ξ,ω)如下式(8)：

${\mathrm{\η}}_e(\mathrm{\ξ},\mathrm{\ω})={\left[\begin{array}{cccc}h& \frac{\∂h}{\∂y}& \frac{{\∂}^{2}h}{\∂y^2}& \frac{{\∂}^{3}h}{\∂y^3}\end{array}\right]}^T---\left(8\right)$

式中，T表示向量转置；

基于状态变量向量ηe(ξ,ω)，可将(7)式写成状态方程的形式如下式(9)：

$\frac{\∂{\mathrm{\η}}_e(\mathrm{\ξ},\mathrm{\ω})}{\∂\mathrm{\ξ}}=F_e(\mathrm{\ξ},\mathrm{\ω},V){\mathrm{\η}}_e(\mathrm{\ξ},\mathrm{\ω})+g_e(\mathrm{\ξ},\mathrm{\ω})---\left(9\right)$

式中，Fe(ξ,ω,V)、ge(ξ,ω)的表达式可根据(7)式得出，具体如下

$\left\{\begin{array}{c}{F}_e(\mathrm{\ξ},\mathrm{\ω},V)=\left[\begin{array}{cccc}0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\\ A_4& A_3& A_2& A_1\end{array}\right]\\ g_e(y,\mathrm{\ω})=0\end{array}\right.---\left(10\right)$

式中，

$\left\{\begin{array}{c}{A}_1=-2\frac{EI\left(1\right)-EI\left(0\right)}{EI\left(\mathrm{\ξ}\right)}\\ A_2=\frac{\mathrm{\ρ}A\left(\mathrm{\ξ}\right)V^2}{EI\left(\mathrm{\ξ}\right)}{l}_j^2\\ A_3=i\mathrm{\ω}\frac{2\mathrm{\ρ}A\left(\mathrm{\ξ}\right)V}{EI\left(\mathrm{\ξ}\right)}{l}_j^2+V^2\frac{\mathrm{\ρ}A\left(1\right)-\mathrm{\ρ}A\left(0\right)}{EI\left(\mathrm{\ξ}\right)}{l}_j^2\\ A_4=-\frac{{\left(i\mathrm{\ω}\right)}^2\[m\left(\mathrm{\ξ}\right)-\mathrm{\ρ}A\left(\mathrm{\ξ}\right)\]}{EI\left(\mathrm{\ξ}\right)}{l}_j^4+\left(i\mathrm{\ω}\right)V\frac{\mathrm{\ρ}A\left(1\right)-\mathrm{\ρ}A\left(0\right)}{EI\left(\mathrm{\ξ}\right)}{l}_j^3\end{array}\right.---\left(11\right)$

根据传递函数理论，方程(9)式的标准解为：

ηe(ξ,ω)＝He(ξ,ω,V)γe(ω)+de(ξ,ω) (12)

式中，γe(ω)为由边界条件给定的位移或力组成的向量，可根据传递函数方法理
论得到表达式如下：

${\mathrm{\γ}}_e\left(\mathrm{\ω}\right)=\left\{\begin{array}{c}h(0,\mathrm{\ω})\\ \frac{\∂h(0,\mathrm{\ω})}{\∂\mathrm{\ξ}}\\ h(1,\mathrm{\ω})\\ \frac{\∂h(1,\mathrm{\ω})}{\∂\mathrm{\ξ}}\end{array}\right\}---\left(13\right)$

式中，He(ξ,ω,V)、fe(ξ,ω)的表达式也可根据传递函数方法理论得到，如下：

$\left\{\begin{array}{c}{H}_e(\mathrm{\ξ},\mathrm{\ω},V)={\mathrm{\Φ}}_F(\mathrm{\ξ},0,\mathrm{\ω},V){\[M_b{\mathrm{\Φ}}_F(0,0,\mathrm{\ω},V)+N_b{\mathrm{\Φ}}_F(1,0,\mathrm{\ω},V)\]}^{-1}\\ d_e(\mathrm{\ξ},\mathrm{\ω})={\∫}_0^{1}G(\mathrm{\ξ},\mathrm{\ζ},\mathrm{\ω},V)g_e(\mathrm{\ζ},\mathrm{\ω})d\mathrm{\ζ}\\ \mathrm{G(}\mathrm{\ξ},\mathrm{\ζ},\mathrm{\ω},V\mathrm{)=}\left\{\begin{array}{cc}{H}_e(\mathrm{\ξ},\mathrm{\ω},V)M_b{\mathrm{\Φ}}_F(0,\mathrm{\ζ},\mathrm{\ω},V)& \mathrm{\ζ}\<\mathrm{\ξ}\\ -H_e(\mathrm{\ξ},\mathrm{\ω},V)N_b{\mathrm{\Φ}}_F(1,\mathrm{\ζ},\mathrm{\ω},V)& \mathrm{\ζ}\>\mathrm{\ξ}\end{array}\right.\\ {\mathrm{\Φ}}_F(\mathrm{\ξ},\mathrm{\ζ},\mathrm{\ω},V)=e^{F_e(\mathrm{\ξ},\mathrm{\ω},V)\mathrm{\ζ}}\end{array}\right.---\left(14\right)$

式中，变量ζ∈(0,1)，Mb为细长体单元前端边界条件选择矩阵，Nb为细长体单元后
端边界条件选择矩阵，可根据传递函数方法理论得到，如下：

$M_b=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right],N_b=\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 1& 0& 0& 0\\ 0& 1& 0& 0\end{array}\right]---\left(15\right)$

由于(10)式中给出ge(y,ω)＝0，由(14)式可知de(ξ,ω)＝0，从而(12)式可简化
为：

ηe(ξ,ω)＝He(ξ,ω,V)γe(ω) (16)

③根据细长体内力关系建立单元平衡方程；

由材料力学梁的弯矩和剪力公式可知，细长体截面上弯矩和剪力的表达式可写
为：

$\left\{\begin{array}{c}M(y,\mathrm{\ω})=EI\frac{{\∂}^{2}h}{\∂y^2}\\ F(y,\mathrm{\ω})=\frac{\∂\left(EI\frac{{\∂}^{2}h}{\∂y^2}\right)}{\∂y}\end{array}\right.---\left(17\right)$

在细长体单元内，并利用EI在单元内的线性变化假设，(17)式可写为：

${\mathrm{\σ}}_e(\mathrm{\ξ},\mathrm{\ω})=\left\{\begin{array}{c}M(\mathrm{\ξ},\mathrm{\ω})\\ F(\mathrm{\ξ},\mathrm{\ω})\end{array}\right\}=\left\{\begin{array}{c}\frac{EI\left(\mathrm{\ξ}\right)}{l_i^2}\frac{{\∂}^{2}h}{\∂{\mathrm{\ξ}}^2}\\ \frac{EI\left(\mathrm{\ξ}\right)}{l_i^3}\frac{{\∂}^{3}h}{\∂{\mathrm{\ξ}}^3}+\frac{EI\left(1\right)-EI\left(0\right)}{l_i^3}\frac{{\∂}^{2}h}{\∂{\mathrm{\ξ}}^2}\end{array}\right\}---\left(18\right)$

引入状态变量向量(8)式，从而(18)可写成如下矩阵形式：

σe(ξ,ω)＝Eeη(ξ)ηe(ξ,ω) (19)

式中，Eeη(ξ)可根据(18)式得到其表达式如下：

$E_{e\mathrm{\η}}\left(\mathrm{\ξ}\right)=\left[\begin{array}{cccc}0& 0& \frac{EI\left(\mathrm{\ξ}\right)}{l_i^2}& 0\\ 0& 0& \frac{EI\left(1\right)-EI\left(0\right)}{l_i^3}& \frac{EI\left(\mathrm{\ξ}\right)}{l_i^3}\end{array}\right]---\left(20\right)$

将(16)式代入(19)式，可得到

σe(ξ,ω)＝Eeη(ξ)He(ξ,ω,V)γe(ω) (21)

式中，γe(ω)可视为细长体单元两节点处的广义位移，如果ξ取0和1，则可得到细
长体单元两节点的广义内力与广义位移之间的关系，从而建立单元平衡方程，即：

令ξ＝0，ξ＝1，则有

${\mathrm{\σ}}_e\left(\mathrm{\ω}\right)=\left\{\begin{array}{c}{\mathrm{\σ}}_e(0,\mathrm{\ω})\\ {\mathrm{\σ}}_e(1,\mathrm{\ω})\end{array}\right\}=\left[\begin{array}{c}{E}_{e\mathrm{\η}}\left(0\right)H_e(0,\mathrm{\ω},V)\\ E_{e\mathrm{\η}}\left(1\right)H_e(1,\mathrm{\ω},V)\end{array}\right]{\mathrm{\γ}}_e\left(\mathrm{\ω}\right)---\left(22\right)$

上式即为表征细长体单元节点广义内力与广义位移关系的单元平衡方程，从而
可视为单元的广义刚度矩阵，可记为：

$K_e(\mathrm{\ω},V)=\left[\begin{array}{c}{E}_{e\mathrm{\η}}\left(0\right)H_e(0,\mathrm{\ω},V)\\ E_{e\mathrm{\η}}\left(1\right)H_e(1,\mathrm{\ω},V)\end{array}\right]---\left(23\right)$

④将细长体单元平衡方程组装形成细长体整体平衡方程；

按照有限元法进行单元组装，可得到细长体整体平衡方程为

K(ω,V)γ(ω)＝F(ω) (24)

式中，K(ω,V)可视为整体刚度矩阵，γ(ω)可视为整体节点位移向量，F(ω)为各
单元节点内力拼装成的向量。

⑤确定颤振速度；

对于细长体，颤振时不需虑重力，本专利中将细长体所受气动力的影响囊括在刚
度矩阵K(ω,V)，

除此以外细长体不受其它力的作用，从而外载荷为零，根据单元节点内力与外载
荷平衡，可得出

F(ω)＝0 (25)

当细长体颤振时，γ(ω)应有非零解，此时整体平衡方程必须满足条件为

det[K(ω,V)]＝0 (26)

由于K(ω,V)为复矩阵，其行列式值等于零的必要条件为矩阵行列式值的实部与
虚部均为零，即

$\left\{\begin{array}{c}\mathrm{Re}\{\det \[K(\mathrm{\ω},V)\]\}=0\\ \mathrm{Im}\{\det \[K(\mathrm{\ω},V)\]\}=0\end{array}\right.---\left(27\right)$

矩阵K(ω,V)中有空速V和圆频率ω两个变量，而(27)式恰好有两个方程，可以定
解。求解上述方程组时，可能会得到满足方程组多个解，即存在多组(V,ω)能满足方程组
(27)式。根据细长体颤振时在某一空速时由稳定转变为不稳定，因而，空速V最小的一组解
(V,ω)应为细长体的颤振速度和相应的颤振频率。将得到满足方程组的空速V和圆频率ω，
分别记为Vcz和ωcz。其中，Vcz即为细长体的颤振速度，ωcz即为细长体的颤振圆频率。

二、所述方法的具体步骤如下：

步骤(一)：将细长体划分为n个单元，测量细长体各单元节点的下列物理参数：

细长体各单元轴向长度[l1 l2 … lj … ln]，单位为米；

细长体单元节点处的横截面积[A1 A2 … Aj … An+1]，单位为米2

细长体单位节点处单位长度质量[m1 m2 … mj … mn+1]，单位为千克/米；

细长体单位节点处抗弯刚度[EI1 EI2 … EIj … EIn+1]，单位为牛顿·米2；

空气密度ρ，单位为千克/米3；

步骤(二)：确定细长体飞行的空速V和圆频率ω的大致范围，假设为

步骤(三)：在范围内划分合适步长△V和△ω，并进行离散，空速V和圆
频率ω的取值为：

步骤(四)：取空速V＝V0，圆频率ω依次取ω0+j△ω(j＝0,1,2,3…)，将空速V和圆
频率ω的取值与步骤(一)中的细长体物理参数代入(11)式，计算各单元系数A1、A2、A3、A4的
值；

步骤(五)：将系数A1、A2、A3、A4的值代入(10)式，得到细长体各单元矩阵Fe(ξ,ω,V)
的值；

步骤(六)：将矩阵Fe(ξ,ω,V)的值代入(14)式，结合(15)式，计算He(ξ,ω,V)的值；

步骤(七)：将He(ξ,ω,V)的值代入(23)式，结合(20)式，计算单元刚度矩阵Ke(ω,
V)的值；

步骤(八)：利用Ke(ω,V)组装整体刚度矩阵K(ω,V)；

步骤(九)：计算Re{det[K(ω,V)]}和Im{det[K(ω,V)]}的值；

步骤(十)：再依次取空速V＝V0+j△V(j＝1,2,3…)，重复步骤(四)至步骤(十)，计
算Re{det[K(ω,V)]}和Im{det[K(ω,V)]}的值；

步骤(十一)：确定满足(27)式的空速V的值，即可确定细长体的颤振速度Vcz。

本发明的有益效果如下：

(1)由于本发明利用梁的横向弯曲振动微分方程来准确描述细长体的横向振动，
而没有采用梁横向振动的低阶模态来近似描述细长体振动，所以计算的细长体颤振速度更
加准确。

(2)本发明求解细长体颤振速度的方法与现有技术相比更加简捷。

附图说明

图1为细长体正视图；

图2为细长体截面图；

图3为细长体结构示意图；

图4为细长体单元剖分示意图；

图5为实施例1的Re[detA]和Im[detA]的等值线图(V∈(0,50),ω∈(0,200π))、

在附图中：1细长体、2第i个细长体、3前节点、4后节点。

具体实施方式

如附图1-5所示，本发明的实施例1如下：

某细长体长度为1.8米，重量22千克，最大截面直径为0.025平方米。

本实施例1的具体计算步骤如下：

步骤(一)：将细长体划分为10个单元，测量细长体各单元节点的下列物理参数：

各单元轴向长度：l1＝l2＝…＝l10＝0.18，单位为米；

各单元节点处的横截面积：A1＝0,A2＝A3＝…＝A10＝0.025,A11＝0.02，单位为米2

各单元节点处的单位长度质量m1＝10,m2＝m3…＝m10＝20，单位为千克/米；

各单元节点处的抗弯刚度EI1＝0.5×103,EI2＝…＝EI11＝1.1×103，单位为牛
顿·米2；

空气密度ρ＝1.225，单位为千克/米3；

步骤(二)：确定飞机机翼的空速V和圆频率ω的大致范围为划分步长
依据发明内容部分中步骤(三)至步骤(十一)，计算Re[detA]和Im[detA]的值。图
3给出了在范围内Re[detA]和Im[detA]的等值线图。从图3中可以发现，在给定
范围内存在满足(26)式的点A，其坐标为(14.7,3020)，它表示细长体颤振速度
Vcz＝3020m/s，颤振频率为14.7Hz(颤振圆频率ωcz＝92.32rad/s)。