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

摘要

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

主权项

一种基于传递函数的细长体颤振速度确定方法，其特征在于所述方法利用梁的横向弯曲振动微分方程和细长体准定常气动力模型，得出细长体颤振微分方程，然后对所述细长体颤振微分方程进行Fourier变换，再运用传递函数法求解细长体颤振速度；一、所述方法需要依据的计算公式：a.利用梁的横向弯曲振动微分方程和细长体准定常气动力模型建立细长体颤振微分方程：本发明利用梁的横向弯曲振动微分方程式(1)来描述细长体的弯曲振动：$\frac{{\partial }^2\left(EI\frac{{\partial }^{2}h}{\partial y^2}\right)}{\partial y^2}+m\frac{{\partial }^{2}h}{\partial t^2}+\Delta p=0---\left(1\right)$式中，h为细长体弯曲振动位移，单位为米，EI为细长体抗弯刚度，单位牛顿·米²，m为细长体单位长度质量，单位为千克，△p为细长体单位长度的气动力，单位为牛顿/米，y为细体长轴向坐标，单位为米，t为时间，单位为秒，为细长体弯曲振动位移h对机轴向坐标y的二阶偏导数，为细长体弯曲振动位移h对时间t的二阶偏导数；根据细长体准定常气动力模型，得到细长体单位长度的气动力△p为下式(2)：$\Delta p=-\{\rho A(\frac{{\partial }^{2}h}{\partial t^2}+V\frac{{\partial }^{2}h}{\partial t\partial y})+\rho V\frac{\partial A}{\partial y}(\frac{\partial h}{\partial t}+V\frac{\partial h}{\partial y})+\rho AV(\frac{{\partial }^{2}h}{\partial t\partial y}+V\frac{{\partial }^{2}h}{\partial y^2})\}---\left(2\right)$式中，V为空速，单位为米/秒，A为细长体横截面积，单位为平方米，ρ为空气密度，单位为千克/立方米；将(2)式代入(1)式得到细长体颤振微分方程式(3)：$\frac{{\partial }^2\left(EI\frac{{\partial }^{2}h}{{\mathrm{ay}}^2}\right)}{{\mathrm{ay}}^2}+m\frac{{\partial }^{2}h}{\partial t^2}-\rho A(\frac{{\partial }^{2}h}{\partial t^2}+V\frac{{\partial }^{2}h}{\partial t\partial y})-\rho V\frac{\partial A}{ay}(\frac{\partial h}{\partial t}+V\frac{\partial h}{ay})-\rho AV(\frac{{\partial }^{2}h}{\partial t\partial y}+V\frac{{\partial }^{2}h}{{\mathrm{ay}}^2})=0---\left(3\right)$b.沿细长体轴向划分n个细长体单元，建立细长体单元的颤振微分方程：一般地，细长体的物理参数EI、m、A等沿其轴向是变化的，为了方便求解，按照有限元法，沿细长体轴向划分n个单元，为了使细长体单元的形式统一颤振微分方程，令：$\xi =\frac{y-y_j}{l_j}---\left(4\right)$式中，y_j为第j个细长体单元前节点距离细长体前端点的距离，l_j为第j个细长体单元的长度，ξ为细长体单元内轴向无量纲坐标，由(4)式可知ξ∈(0,1)；当划分单元足够多时，假设EI、m，A在每个细长体单元内近似线性变化是可理的，则在第j个细长体单元内有如下线性表达式：$\left\{\begin{array}{c}EI\left(\xi \right)=EI\left(0\right)+\frac{EI\left(1\right)-EI\left(0\right)}{l_j}\xi \\ m\left(\xi \right)=m\left(0\right)+\frac{m\left(1\right)-m\left(0\right)}{l_j}\xi \\ A\left(\xi \right)=A\left(0\right)+\frac{A\left(1\right)-A\left(0\right)}{l_j}\xi \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(\xi \right)}{l_j^4}\frac{{\partial }^{4}h}{\partial {\xi }^4}+\frac{2\left(EI\right(1)-EI(0\left)\right)}{l_j^4}\frac{{\partial }^{3}h}{\partial {\xi }^3}-\frac{\rho A\left(\xi \right)V^2}{l_j^2}\frac{{\partial }^{2}h}{\partial {\xi }^2}\\ -\frac{2\rho A\left(\xi \right)V}{l_j}\frac{{\partial }^{2}h}{\partial t\partial \xi }-{\mathrm{\rho V}}^2\frac{A\left(1\right)-A\left(0\right)}{l_j^2}\frac{\partial h}{\partial \xi }+[m\left(\xi \right)-\rho A\left(\xi \right)]\frac{{\partial }^{2}h}{\partial {\xi }^2}-\rho V\frac{A\left(1\right)-A\left(0\right)}{l_j}\frac{\partial h}{\partial t}=0\end{array}---\left(6\right)$c.利用传递函数法求解细长体颤振速度①对(6)式中有关时间t的变量作Fourier变换，可得到下式(7)：$\begin{array}{c}\frac{EI\left(\xi \right)}{l_j^4}\frac{{\partial }^{4}h}{\partial {\xi }^4}+\frac{2\left(EI\right(1)-EI(0\left)\right)}{l_j^4}\frac{{\partial }^{3}h}{\partial {\xi }^3}-\frac{\rho A\left(\xi \right)V^2}{l_j^2}\frac{{\partial }^{2}h}{\partial {\xi }^2}\\ -\left[i\omega \frac{2\rho A\left(\xi \right)V}{l_j}+V^2\frac{\sigma A\left(1\right)-\sigma A\left(0\right)}{l_j^2}\right]\frac{\partial h}{\partial \xi }+\left[{\left(i\omega \right)}^2\right[m\left(\xi \right)-\rho A\left(\xi \right)]-\left(i\omega \right)V\frac{\rho A\left(1\right)-\rho A\left(0\right)}{l_j}]h=0\end{array}---\left(7\right)$式中，虚数ω为角频率，单位为弧度/秒；②将细长体单元的颤振微分方程改写为状态空间方程形式；为了便于应用传递函数理论，定义一个状态变量的向量η_e(ξ,ω)如下式(8)：${\eta }_e(\xi ,\omega )={\left[\begin{array}{cccc}h& \frac{\partial h}{\partial y}& \frac{{\partial }^{2}h}{\partial y^2}& \frac{{\partial }^{3}h}{\partial y^3}\end{array}\right]}^T---\left(8\right)$式中，T表示向量转置；基于状态变量向量η_e(ξ,ω)，可将(7)式写成状态方程的形式如下式(9)：$\frac{\partial {\eta }_e(\xi ,\omega )}{\partial \xi }=F_e(\xi ,\omega ,V){\eta }_e(\xi ,\omega )+g_e(\xi ,\omega )---\left(9\right)$式中，F_e(ξ,ω,V)、g_e(ξ,ω)的表达式可根据(7)式得出，具体如下$\left\{\begin{array}{c}{F}_e(\xi ,\omega ,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,\omega )=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(\xi \right)}\\ A_2=\frac{\rho A\left(\xi \right)V^2}{EI\left(\xi \right)}{l}_j^2\\ A_3=i\omega \frac{2\rho A\left(\xi \right)V}{EI\left(\xi \right)}{l}_j^3+V^2\frac{\rho A\left(1\right)-\rho A\left(0\right)}{EI\left(\xi \right)}{l}_j^2\\ A_4=-\frac{{\left(i\omega \right)}^2[m\left(\xi \right)-\rho A\left(\xi \right)]}{EI\left(\xi \right)}{l}_j^4+\left(i\omega \right)V\frac{\rho A\left(1\right)-\rho A\left(0\right)}{EI\left(\xi \right)}{l}_j^3\end{array}\right.---\left(11\right)$根据传递函数理论，方程(9)式的标准解为：η_e(ξ,ω)＝H_e(ξ,ω,V)γ_e(ω)+d_e(ξ,ω)    (12)式中，γ_e(ω)为由边界条件给定的位移或力组成的向量，可根据传递函数方法理论得到表达式如下：${\gamma }_e\left(\omega \right)=\left\{\begin{array}{c}h(0,\omega )\\ \frac{\partial h(0,\omega )}{\partial \xi }\\ h(1,\omega )\\ \frac{\partial h(1,\omega )}{\partial \xi }\end{array}\right\}---\left(13\right)$式中，H_e(ξ,ω,V)、f_e(ξ,ω)的表达式也可根据传递函数方法理论得到，如下：$\left\{\begin{array}{c}{H}_e(\xi ,\omega ,V)={\Phi }_F(\xi ,0,\omega ,V){[M_b{\Phi }_F(0,0,\omega ,V)+N_b{\Phi }_F(1,0,\omega ,V)]}^{-1}\\ d_e(\xi ,\omega )={\int }_0^1(\xi ,\zeta ,\omega ,V)g_e(\zeta ,\omega )d\zeta \\ G(\xi ,\zeta ,\omega ,V)=\left\{\begin{array}{cc}{H}_e(\xi ,\omega ,V)M_b{\Phi }_F(0,\zeta ,\omega ,V)& \zeta <\xi \\ -H_e(\xi ,\omega ,V)N_b{\Phi }_F(1,\zeta ,\omega ,V)& \zeta >\xi \end{array}\right.\\ {\Phi }_F(\xi ,\zeta ,\omega ,V)=e^{F_e(\xi ,\omega ,V)\zeta }\end{array}\right.---\left(14\right)$式中，变量ζ∈(0,1)，M_b为细长体单元前端边界条件选择矩阵，N_b为细长体单元后端边界条件选择矩阵，可根据传递函数方法理论得到，如下：$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)式中给出g_e(y,ω)＝0，由(14)式可知d_e(ξ,ω)＝0，从而(12)式可简化为：η_e(ξ,ω)＝H_e(ξ,ω,V)γ_e(ω)     (16)③根据细长体内力关系建立单元平衡方程；由材料力学梁的弯矩和剪力公式可知，细长体截面上弯矩和剪力的表达式可写为：$\left\{\begin{array}{c}M(y,\omega )=EI\frac{{\partial }^{2}h}{\partial y^2}\\ F(y,\omega )=\frac{\partial \left(EI\frac{{\partial }^{2}h}{\partial y^2}\right)}{\partial y}\end{array}\right.---\left(17\right)$在细长体单元内，并利用EI在单元内的线性变化假设，(17)式可写为：${\sigma }_e(\xi ,\omega )=\left\{\begin{array}{c}M(\xi ,\omega )\\ F(\xi ,\omega )\end{array}\right\}=\left\{\begin{array}{c}\frac{EI\left(\xi \right)}{l_i^2}\frac{{\partial }^{2}h}{\partial {\xi }^2}\\ \frac{EI\left(\xi \right)}{l_i^3}\frac{{\partial }^{3}h}{\partial {\xi }^3}+\frac{EI\left(1\right)-EI\left(0\right)}{l_i^3}\frac{{\partial }^{2}h}{\partial {\xi }^2}\end{array}\right\}---\left(18\right)$引入状态变量向量(8)式，从而(18)可写成如下矩阵形式：σ_e(ξ,ω)＝E_eη(ξ)η_e(ξ,ω)    (19)式中，E_eη(ξ)可根据(18)式得到其表达式如下：$E_{e\eta }\left(\xi \right)=\left[\begin{array}{cccc}0& 0& \frac{EI\left(\xi \right)}{l_i^2}& 0\\ 0& 0& \frac{EI\left(1\right)-EI\left(0\right)}{l_i^3}& \frac{EI\left(\xi \right)}{l_i^3}\end{array}\right]---\left(20\right)$将(16)式代入(19)式，可得到σ_e(ξ,ω)＝E_eη(ξ)H_e(ξ,ω,V)γ_e(ω)    (21)式中，γ_e(ω)可视为细长体单元两节点处的广义位移，如果ξ取0和1，则可得到细长体单元两节点的广义内力与广义位移之间的关系，从而建立单元平衡方程，即：令ξ＝0，ξ＝1，则有${\sigma }_e\left(\omega \right)=\left\{\begin{array}{c}{\sigma }_e(0,\omega )\\ {\sigma }_e(1,\omega )\end{array}\right\}=\left[\begin{array}{c}{E}_{e\eta }\left(0\right)H_e(0,\omega ,V)\\ E_{e\eta }\left(1\right)H_e(1,\omega ,V)\end{array}\right]{\gamma }_e\left(\omega \right)---\left(22\right)$上式即为表征细长体单元节点广义内力与广义位移关系的单元平衡方程，从而可视为单元的广义刚度矩阵，可记为：$K_e(\omega ,V)=\left[\begin{array}{c}{E}_{e\eta }\left(0\right)H_e(0,\omega ,V)\\ E_{e\eta }\left(1\right)H_e(1,\omega ,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}\left\{\det \right[K(\omega ,V)\left]\right\}=0\\ \mathrm{Im}\left\{\det \right[K(\omega ,V)\left]\right\}=0\end{array}\right.---\left(27\right)$矩阵K(ω,V)中有空速V和圆频率ω两个变量，而(27)式恰好有两个方程，可以定解。求解上述方程组时，可能会得到满足方程组多个解，即存在多组(V,ω)能满足方程组(27)式。根据细长体颤振时在某一空速时由稳定转变为不稳定，因而，空速V最小的一组解(V,ω)应为细长体的颤振速度和相应的颤振频率。将得到满足方程组的空速V和圆频率ω，分别记为V_cz和ω_cz。其中，V_cz即为细长体的颤振速度，ω_cz即为细长体的颤振圆频率。二、所述方法的具体步骤如下：步骤(一)：将细长体划分为n个单元，测量细长体各单元节点的下列物理参数：细长体各单元轴向长度[l₁ l₂ … l_j … l_n]，单位为米；细长体单元节点处的横截面积[A₁ A₂ … A_j … A_n+1]，单位为米²细长体单位节点处单位长度质量[m₁ m₂ … m_j … m_n+1]，单位为千克/米；细长体单位节点处抗弯刚度[EI₁ EI₂ … EI_j … EI_n+1]，单位为牛顿·米²；空气密度ρ，单位为千克/米³；步骤(二)：确定细长体飞行的空速V和圆频率ω的大致范围，假设为步骤(三)：在范围内划分合适步长△V和△ω，并进行离散，空速V和圆频率ω的取值为：步骤(四)：取空速V＝V₀，圆频率ω依次取ω₀+j△ω(j＝0,1,2,3…)，将空速V和圆频率ω的取值与步骤(一)中的细长体物理参数代入(11)式，计算各单元系数A₁、A₂、A₃、A₄的值；步骤(五)：将系数A₁、A₂、A₃、A₄的值代入(10)式，得到细长体各单元矩阵F_e(ξ,ω,V)的值；步骤(六)：将矩阵F_e(ξ,ω,V)的值代入(14)式，结合(15)式，计算H_e(ξ,ω,V)的值；步骤(七)：将H_e(ξ,ω,V)的值代入(23)式，结合(20)式，计算单元刚度矩阵K_e(ω,V)的值；步骤(八)：利用K_e(ω,V)组装整体刚度矩阵K(ω,V)；步骤(九)：计算Re{det[K(ω,V)]}和Im{det[K(ω,V)]}的值；步骤(十)：再依次取空速V＝V₀+j△V(j＝1,2,3…)，重复步骤(四)至步骤(十)，计算Re{det[K(ω,V)]}和Im{det[K(ω,V)]}的值；步骤(十一)：确定满足(27)式的空速V的值，即可确定细长体的颤振速度V_cz。