基于传递函数的机翼颤振速度确定方法

摘要

本发明涉及基于传递函数的机翼颤振速度确定方法，所述方法利用机翼弯扭振动微分方程和Theodrosen非定常气动力模型，得出长直机翼颤振微分方程，然后对所述机翼颤振微分方程进行Fourier变换，再运用传递函数法求解机翼颤振速度。本发明的有益效果如下：（1）由于本发明利用机翼弯扭振动微分方程来准确描述机翼振动，而没有采用机翼振动的低阶模态来近似描述机翼振动，所以计算的机翼颤振速度更加准确；（2）本发明求解机翼颤振速度的方法与现有技术相比更加简捷。

主权项

一种基于传递函数的机翼颤振速度确定方法，其特征在于所述方法利用机翼弯扭振动微分方程和Theodrosen非定常气动力模型，得出长直机翼颤振微分方程，然后对所述机翼颤振微分方程进行Fourier变换，再运用传递函数法求解机翼颤振速度；一、所述方法需要依据的计算公式：a.利用机翼弯扭振动微分方程建立机翼颤振微分方程：本发明利用机翼弯扭振动微分方程式(1)来描述机翼弯扭振动：$\left\{\begin{array}{c}EI\frac{{\partial }^{4}h}{\partial y^4}+m\frac{{\partial }^{2}h}{\partial t^2}-{\mathrm{mx}}_{\alpha }\frac{{\partial }^2\alpha }{\partial t^2}-L_h=0\\ GJ\frac{{\partial }^2\alpha }{\partial y^2}-I_{\alpha }\frac{{\partial }^2\alpha }{\partial t^2}+{\mathrm{mx}}_{\alpha }\frac{{\partial }^{2}h}{\partial t^2}+T_{\alpha }=0\end{array}\right.---\left(1\right)$式中，h为机翼弯曲振动位移，单位为米，α为机翼扭转振动转角，单位为弧度，EI为机翼抗弯刚度，单位牛顿·米²，GJ为机翼抗扭刚度，单位为牛顿·米²，m为机翼单位长度质量，单位为千克，I_α为单位长度机翼绕弹性轴的转动惯量，单位为千克·米²，x_α为机翼弹性轴到机翼截面重心的距离，单位为米，L_h为机翼单位长度的升力，单位为牛顿/米，T_α机翼单位长度的扭矩，单位为牛顿，y为机翼展向坐标，单位为米，t为时间，单位为秒，为机翼弯曲振动位移h对机翼展向坐标y的四阶偏导数，为机翼扭转振动转角α对机翼展向坐标y的二阶偏导数，为机翼弯曲振动位移h对时间t的二阶偏导数，为机翼扭转振动转角α对时间t的二阶偏导数；b.根据Theodrosen非定常气动力模型，得到机翼单位长度的升力L_h和机翼单位长度的扭矩T_α为下式(2)：$\left\{\begin{array}{c}{L}_h={\mathrm{\pi \rho b}}^2(-\frac{{\partial }^{2}h}{\partial t^2}+V\frac{\partial \alpha }{\partial t}-b\bar{a}\frac{{\partial }^2\alpha }{\partial t^2})+2\pi \rho VbC\left(k\right)(V\alpha -\frac{\partial h}{\partial t}+b(\frac{1}{2}-\bar{a})\frac{\partial \alpha }{\partial t})\\ T_{\alpha }={\mathrm{\pi \rho b}}^2(-b\bar{a}\frac{{\partial }^{2}h}{\partial t^2}-Vb(\frac{1}{2}-\bar{a})\frac{\partial \alpha }{\partial t}-b^2(\frac{1}{8}+{\bar{a}}^2)\frac{{\partial }^2\alpha }{\partial t^2})\\ +2{\mathrm{\pi \rho Vb}}^2\left(\frac{1}{2}+\bar{a}\right)C\left(k\right)(V\alpha -\frac{\partial h}{\partial t}+b(\frac{1}{2}-\bar{a})\frac{\partial \alpha }{\partial t})\end{array}\right.---\left(2\right)$式中，V为空速，单位为米/秒，ρ为空气密度，单位为千克/立方米，b为机翼的半弦长，单位为米，C(k)为Theodrosen函数，为减缩频率，ω为圆频率，单位为弧度/秒，为机翼弹性轴到机翼弦长中点的距离占半弦长的百分比；由于减缩频率k是圆频率ω和空速V的函数，将C(k)写为C(ω,V)，将(2)式代入(1)式得到机翼颤振微分方程式(3)：$\left\{\begin{array}{c}EI\frac{{\partial }^{4}h}{\partial y^4}+m\frac{{\partial }^{2}h}{\partial t^2}-{\mathrm{mx}}_{\alpha }\frac{{\partial }^2\alpha }{\partial t^2}-{\mathrm{\pi \rho b}}^2\left(-\frac{{\partial }^{2}h}{\partial t^2}+V\frac{\partial \alpha }{\partial t}-b\bar{a}\frac{{\partial }^2\alpha }{\partial t^2}\right)\\ -2\pi \rho VbC(\omega ,V)\left(V\alpha -\frac{\partial h}{\partial t}+b\left(\frac{1}{2}-\bar{a}\right)\frac{\partial \alpha }{\partial t}\right)=0\\ GJ\frac{{\partial }^{2}h}{\partial y^2}-I_{\alpha }\frac{{\partial }^2\alpha }{\partial t^2}+{\mathrm{mx}}_{\alpha }\frac{{\partial }^{2}h}{\partial t^2}+{\mathrm{\pi \rho b}}^2\left(-b\bar{a}\frac{{\partial }^{2}h}{\partial t^2}-Vb\left(\frac{1}{2}-\bar{a}\right)\frac{\partial \alpha }{\partial t}-b^2\left(\frac{1}{8}+{\bar{a}}^2\right)\frac{{\partial }^2\alpha }{\partial t^2}\right)\\ +2{\mathrm{\pi \rho Vb}}^2\left(\frac{1}{2}+\bar{a}\right)C(\omega ,V)\left(V\alpha -\frac{\partial h}{\partial t}+b\left(\frac{1}{2}-\bar{a}\right)\frac{\partial \alpha }{\partial t}\right)=0\end{array}\right.---\left(3\right)$c.利用传递函数法求解机翼颤振速度①对(3)式作Fourier变换，并经过整理可得到下式(4)：$\left\{\begin{array}{c}\frac{{\partial }^4\alpha }{\partial y^4}=A_1(\omega ,V)h+B_1(\omega ,V)\alpha \\ \frac{{\partial }^2\alpha }{\partial y^2}=A_2(\omega ,V)h+B_2(\omega ,V)\alpha \end{array}\right.---\left(4\right)$式中，A₁(ω,V)、A₂(ω,V)、B₁(ω,V)、B₂(ω,V)的具体表达式如下：$\{\begin{array}{c}{A}_1(\omega ,V)=\frac{{\omega }^2[(m-{\mathrm{\pi \rho b}}^2)+2\pi \rho b\frac{V}{\omega }C(\omega ,V)i]}{EI}\\ B_1(\omega ,V)=\frac{{\omega }^2[{\mathrm{\pi \rho b}}^2(b\bar{a}+{\mathrm{mx}}_{\alpha })+({\mathrm{\pi \rho b}}^2\frac{V}{\omega }+2{\mathrm{\pi \rho b}}^2\frac{V}{\omega }C(\omega ,V\left)\right(0.5-\bar{a}\left)\right)i+2\pi \rho b\frac{V^2}{{\omega }^2}C(\omega ,V)]}{EI}\\ A_2(\omega ,V)=\frac{{\omega }^2[({\mathrm{mx}}_{\alpha }-{\mathrm{\pi \rho b}}^3\bar{a})+2\pi \rho \frac{V}{\omega }{b}^2(0.5+\bar{a}\left)C\right(\omega ,V\left)i\right]}{GJ}\\ B_2(\omega ,V)=\frac{{\omega }^2[(I_{\alpha }+{\mathrm{\pi \rho b}}^4(0.125+{\bar{a}}^2\left)\right)+[2\pi \rho \frac{V}{\omega }{b}^3\left(0.5+\bar{a}\right)C\left(\omega ,V\right)\left(0.5-\bar{a}\right)-{\mathrm{\pi \rho b}}^3\frac{V}{\omega }\left(0.5-\bar{a}\right)]i+2{\mathrm{\pi \rho b}}^2\frac{V^2}{{\omega }^2}(0.5+\bar{a}\left)C(\omega ,V)\right]}{GJ}\end{array}---\left(5\right)$其中，虚数$i=\sqrt{-1};$②采用传递函数法求解上述(4)式，确定机翼的颤振速度；为了便于应用传递函数理论，定义状态变量向量如下式(6)：$\eta (y,\omega )={\left[\begin{array}{cccccc}h& \frac{\partial h}{\partial y}& \frac{{\partial }^{2}h}{\partial y^2}& \frac{{\partial }^{3}h}{\partial y^3}& \alpha & \frac{\partial \alpha }{\partial y}\end{array}\right]}^T---\left(6\right)$式中，T表示向量转置；从而，将(4)式写成状态方程的形式如下式(7)：$\frac{\partial \eta (y,\omega )}{\partial y}=F(\omega ,V)\eta (y,\omega )+g(y,\omega )---\left(7\right)$式中，$F(\omega ,V)=\left[\begin{array}{cccccc}0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ A_1(\omega ,V)& 0& 0& 0& B_1(\omega ,V)& 0\\ 0& 0& 0& 0& 0& 1\\ A_2(\omega ,V)& 0& 0& 0& B_2(\omega ,V)& 0\end{array}\right]---\left(8\right)$g(y,ω)＝0边界条件为：M_bη(y＝0,ω)+N_bη(y＝l,ω)＝γ(ω)        (9)式中，η(0,ω)为η(y,ω)在机翼根端y＝0处的值，η(l,ω)为η(y,ω)在机翼梢端y＝l处的值，l为机翼的半展长，单位为米，M_b为机翼根端边界条件选择矩阵，N_b为机翼梢端边界条件选择矩阵，γ(ω)为由边界条件给定的位移或力组成的向量，其表达式分别为：$M_b=\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right],N_b=\left[\begin{array}{cccccc}0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right],\gamma \left(\omega \right)=\left\{\begin{array}{c}h(y=0,\omega )\\ \frac{\partial h}{\partial x}(y=0,\omega )\\ \alpha (y=0,\omega )\\ \frac{{\partial }^{2}h}{\partial x^2}(y=l,\omega )\\ \frac{{\partial }^{3}h}{\partial x^3}(y=l,\omega )\\ \frac{\partial \alpha }{\partial x}(y=l,\omega )\end{array}\right\}$根据传递函数理论，方程(7)式的解为：$\eta (y,\omega )=H(y,\omega ,V)\gamma \left(\omega \right)+{\int }_0^{l}G(y,\xi ,\omega ,V)g(\xi ,\omega )d\xi ---\left(10\right)$式中，G(y,ξ,ω,V)为状态空间方程的域内传递函数，H(y,ω,V)为状态空间方程的边界传递函数，其表达式分别为：$\left\{\begin{array}{c}G(y,\xi ,\omega ,V)=\left\{\begin{array}{c}H(y,\omega ,V)M_b{e}^{-F(\omega ,V)\xi },\xi \le y\\ -H(y,\omega ,V)N_b{e}^{-F(\omega ,V)(l-\xi )},\xi >y\end{array}\right.\\ H(y,\omega ,V)=e^{F(\omega ,V)y}{\left[M_b+N_b{e}^{F(\omega ,V)l}\right]}^{-1}\end{array}\right.---\left(11\right)$式中，变量ξ∈(0,l)，为机翼展向的坐标，单位为米；机翼颤振时弯曲振动位移h和扭转振动转角α的振幅为非零常数，即(10)式有非零解，根据传递函数理论可知，(10)式有非零解的充分必要条件为：det[M_b+N_be^F(ω,V)l]^‑1＝0           (12)令A＝[M_b+N_be^F(ω,V)L]^‑1，由于A为复矩阵，其行列式值等于零的必要条件为矩阵行列式值的实部与虚部均为零，即$\left\{\begin{array}{c}\mathrm{Re}\left[\det A\right]=0\\ \mathrm{Im}\left[\det A\right]=0\end{array}\right.---\left(13\right)$求解上述(13)式，可以得到满足方程组的空速V和圆频率ω，分别记为V_cz和ω_cz；其中，V_cz即为机翼的颤振速度，ω_cz即为机翼的颤振圆频率；二、所述方法的具体步骤如下：步骤(一)：测量长直机翼的下列物理参数：半弦长b，单位为米；半展长l，单位为米；单位长度机翼质量m，单位为千克/米；机翼弹性轴(2)到机翼弦长中点的距离占半弦长的百分比机翼弹性轴(2)到机翼截面重心的距离x_α，单位为米；单位长度机翼绕弹性轴的转动惯量I_α，单位为千克·米²；机翼的抗弯刚度EI，单位为牛顿·米²；机翼的抗扭刚度GJ，单位为牛顿·米²；空气密度ρ，单位为千克/立方米；步骤(二)：确定飞机机翼的空速V和圆频率ω的大致范围，假设为$\left\{\begin{array}{c}V\in (V_0,V_1)\\ \omega \in ({\omega }_0,{\omega }_1)\end{array}\right.;$步骤(三)：在$\left\{\begin{array}{c}V\in (V_0,V_1)\\ \omega \in ({\omega }_0,{\omega }_1)\end{array}\right.$范围内划分合适步长ΔV和Δω，并进行离散，空速V和圆频率ω的取值为：$\left\{\begin{array}{cc}V=V_0+i\Delta V& (i=0,1,2,3...)\\ \omega ={\omega }_0+j\Delta \omega & (j=0,1,2,3...)\end{array}\right.;$步骤(四)：取空速V＝V₀，圆频率ω依次取ω₀+jΔω(j＝0,1,2,3…)，将空速V和圆频率ω的取值与步骤(一)中的机翼各物理参数代入(5)式，得到系数A₁(ω,V)、A₂(ω,V)、B₁(ω,V)、B₂(ω,V)的值；步骤(五)：将系数A₁(ω,V)、A₂(ω,V)、B₁(ω,V)、B₂(ω,V)的值代入(8)式，得到矩阵F(ω,V)的值；步骤(六)：将矩阵F(ω,V)的值代入(13)式，计算Re[det A]和Im[det A]的值；步骤(七)：再依次取空速V＝V₀+jΔV(j＝1,2,3…)，重复步骤(四)至步骤(六)，计算Re[det A]和Im[det A]的值；步骤(八)：确定同时满足(13)式的空速V的值，即为机翼的颤振速度V_cz。