带内埋质量超声速悬臂板的颤振特性

Flutter characteristics of supersonic cantilever plate with embedded mass

摘要

为提高具有内埋质量升力面的颤振速度，将其简化为带集中质量的矩形悬臂板。基于假设模态法和一阶活塞理论得到颤振模型，进而考察其颤振特性。首先，基于Kirchhoff薄板理论建立矩形悬臂板的结构动力学模型，并利用假设模态法得到带内埋质量时的模态属性。其次，引入一阶活塞理论近似求解在超声速条件下的模态气动力，代入结构动力学模型得到颤振方程的时域表达。最后，通过引入指数形式的模态位移将颤振方程变换至频域并通过 p - k法求解颤振速度和颤振频率。通过改变内埋件的位置、质量和数目等参数，得到其对超声速悬臂板颤振特性的影响规律。研究结果表明，在悬臂板前缘靠近翼梢区域内埋小质量物体可有效提高升力面的颤振速度，而在翼根附近，颤振速度对内埋质量的变化不敏感。

关键词

颤振 / 内埋质量 / 活塞理论 / 假设模态法 / 超声速悬臂板

Key words

flutter / embedded mass / piston theory / the assumed mode method / supersonic cantilever plate

超声速飞行器的飞行速度包线受到来自结构、气动和控制等方面的共同制约。随着飞行速度的不断增大，飞行器结构会从气流中不断汲取能量而发生振动。当飞行速度达到某一数值时，升力面将维持一种等幅振荡状态，此时的飞行速度被认为是颤振边界。颤振是一种危险的行为，飞行器结构可能会在数秒内发生解体。另一方面，设计先进飞行器时，通常希望其工作于多任务剖面，需要携带的质量附着物也是多样的，例如弹体、主（副）油箱、自毁装置等。同时，为了不影响升力面的气动外形，这些装置通常以内埋形式置于机翼和机身内部。所以，对带内埋质量机翼颤振特性的研究是迫切的。

2003年，Chiba等^［1］采用Rayleigh-Ritz法对附加有弹簧-质量系统的悬臂薄板进行了耦合自由振动分析，阐明了附加质量、弹簧常数及附加系统位置对系统耦合振动特性的影响。之后，赵永辉^［2-3］详细研究了超声速悬臂板及超声速翼面颤振的分析方法。肖艳平等^［4］研究了复合材料机翼的颤振速度随材料耦合刚度的变化及在不同挂载刚度、挂载位置情况下颤振速度随材料耦合刚度的变化。陈晓明等^［5］结合理论分析、数值模拟及实验测量计算了含有集中质量悬臂薄板模型的固有频率，经过对比发现三者有较好的一致性。Aksencer等^［6］使用Ritz法研究了带有附加点质量的旋转层压复合梁的自由振动。Rahmane等^［7］研究了在夹紧/自由边界条件的弯曲振动下，附加质量对复合材料层合板动态特性的影响，并得到了质量大小对固有频率的影响。Pacheco等^［8］采用非线性能量阱来抑制面板颤振并降低极限循环振荡的强度，探究了非线性能量阱的参数如何影响气动弹性响应的本质。Molina等^［9］提出了一种带附加质量的矩阵计算方法，在颤振计算中可有效地提升准确性和计算效率。Bera等^［10］使用可控小翼和旋转质量阻尼器的组合进行了桥面的颤振抑制，桥面在颤振时的垂直响应可通过旋转质量阻尼器进行有效衰减。Zhou等^［11-12］尝试使用非线性能量阱（nonlinear energy sink，NES）来抑制面板颤振，利用Galerkin方法建立了超声速流动中具有NES的二维平板的非线性气动弹性模型。Tian等^［13］设计了一个嵌入局部非线性谐振器的周期性阵列的超材料板，与非超材料板进行比较，证明其有更高的颤振边界和更低的振幅。Bahrami等^［14］提出了一种新型流固耦合有限元模型来研究超声速流动下矩形板的动态稳定性，对具有不同几何特性和边界条件的矩形板进行了计算。祁武超等^［15］基于假设模态法提出一种带集中质量弯扭组合梁模态分析手段，结合片条理论考察了外挂物不同质量及位置分布对机翼颤振的影响。Qi等^［16］用假设模态法求得附着集中质量壁板的模态函数，使用一阶活塞理论和p - k法，对不同质量附着的情况进行颤振分析，预测有质量附着时壁板的颤振速度。Tian等^［17］提出了一种罚函数的刚度矩阵约束关系处理方法，结合一阶活塞理论计算了不同弹性支撑边界条件的薄壁板的颤振特性，并通过增加阻尼来对颤振速度进行增强和恢复。Muñoz等^［18］提出了一种计算亚声速、跨声速和超声速势流中机翼截面在拉普拉斯域下的非定常气动力求解方法，对三自由度机翼截面的主动颤振抑制进行研究，得到了合适的控制方法。Uymaz^［19］ 使用三维Ritz解分析了质量和温度对任意位置承载点质量的功能梯度板的自由振动的影响。Kalosha等^［20］考虑了一个带有附加弹簧质量系统的简支欧拉-伯努利梁的数学模型，该模型由分布式压电致动器和集中力控制。Dadoulis等^［21］研究了在外部基础运动下具有大质量附件梁的解析解。Mahmoud^［22］介绍了使用集中质量传递矩阵对承载任意数量质量点的非均匀和阶梯式轴向功能梯度梁进行自由振动分析的方法。Simiriotis等^［23］介绍了直接预测法、在频域内求解的p - k法和有理插值法等3种方法求解颤振方程，所预测的颤振速度几乎相同。

现有文献主要研究了附加质量对板自由振动特性的影响，而在讨论附加质量对板颤振特性影响方面仍处于起步阶段，尤其是具有多质量附着时板的颤振特性方面未见相关研究成果。另外，在初始设计过程中，需要不断考察板的尺寸效应及内埋件质量属性的变化（大小、位置、数目等因素）所带来的效应，重复建立有限元模型，并利用Nastran进行分析，从而产生极大的计算代价和时间成本。从解析的角度构建气动弹性模型则具有明显优势，一方面可以有效地提高计算效率，另一方面则更容易观察质量参数变化对系统本质属性的影响。因此，本文基于假设模态法得到带内埋质量矩形板的模态特性，并结合一阶活塞理论求解悬臂板的颤振分析模型，通过改变内埋件位置、质量大小和数目等参数，分析其对超声速悬臂板颤振特性的影响规律。

1 带内埋质量悬臂板的颤振方程

考察如图1所示的带内埋质量悬臂板的几何模型，展向和弦向的长度分别为a和b，板厚为h，这里h远小于a和b，即所考虑的悬臂板为薄板。在所建立的全局坐标系中，坐标原点o位于靠近固定端的前缘处，x轴正向为顺气流方向，y轴正向为悬臂板展向方向，z轴正向为悬臂板的法线方向（右手坐标系）。悬臂板在超声速气流中运动，相对来流速度记为 U_{\mathrm{\infty }}。板面内附着有N_p 个内埋件，可将其等效为集中质量，其中第p个内埋件的质量记为m_p，挂载点记为（x_p，y_p，z_p ）。由于附加质量的内埋属性和悬臂板的薄板特征，这里可忽略内埋件z向坐标的变化所带来的影响。

由于内埋质量件的存在，悬臂板的质量属性发生了比较大的变化，相应地其固有特性也会发生改变。本文引入假设模态法以得到带内埋质量悬臂板的模态属性。基于Kirchhoff薄板理论，悬臂板的弯曲势能可写为

\begin{array}{l}{U}_{\mathrm{B}}=\frac{D_b}{\mathrm{2}}{\int }_{\mathrm{0}}^a{\int }_{\mathrm{0}}^b\left[{\left(\frac{{\partial }^{\mathrm{2}}w}{\partial x^{\mathrm{2}}}\right)}^{\mathrm{2}}+{\left(\frac{{\partial }^{\mathrm{2}}w}{\partial y^{\mathrm{2}}}\right)}^{\mathrm{2}}\right]\mathrm{d}x\mathrm{d}y+\\ D_b{\int }_{\mathrm{0}}^a{\int }_{\mathrm{0}}^b\left[\nu \frac{{\partial }^{\mathrm{2}}w}{\partial x^{\mathrm{2}}}\frac{{\partial }^{\mathrm{2}}w}{\partial y^{\mathrm{2}}}+\left(\mathrm{1}-\nu \right){\left(\frac{{\partial }^{\mathrm{2}}w}{\partial x\partial y}\right)}^{\mathrm{2}}\right]\mathrm{d}x\mathrm{d}y\end{array} （1）

式中：D_b 为板的抗弯刚度；w为板的弯曲变形； \nu为泊松比。

T _B和T _A分别为悬臂板发生振动时自身动能附加质量件的动能，则系统的总动能为

\begin{array}{l}T=T_{\mathrm{B}}+T_{\mathrm{A}}\\ =\frac{\mathrm{1}}{\mathrm{2}}{\int }_{\mathrm{0}}^a{\int }_{\mathrm{0}}^b\rho h{\left(\frac{\partial w}{\partial t}\right)}^{\mathrm{2}}\mathrm{d}x\mathrm{d}y\\ +\frac{\mathrm{1}}{\mathrm{2}}{{\displaystyle \sum _{p=\mathrm{1}}^{N_p}}{m}_p\left(\frac{\partial w\left(x_p,y_p,t\right)}{\partial t}\right)}^{\mathrm{2}}\end{array} （2）

式中： \rho为悬臂板的质量密度。

取带内埋质量悬臂板的假设模态为

\begin{array}{l}{\tilde{\psi }}_i\left(x,y\right)={\varphi }_m\left(x\right){\phi }_n\left(y\right),\\ m=\mathrm{1,2},\cdots ,N_m;n=\mathrm{1,2},\cdots ,N_n\end{array} （3）

式中： {\varphi }_m \left(x\right)和 {\phi }_n \left(y\right)分别为两端自由梁和悬臂梁的振型函数；变量i根据 {\varphi }_m和 {\phi }_n组合得出的悬臂板固有频率的大小来排序，具体表达式为^［4］

{\varphi }_m \left(x\right) = \left\{\begin{array}{ll}\mathrm{2},& m=\mathrm{1}\\ \mathrm{2}\left(\mathrm{1}-\mathrm{2}x/b\right),& m=\mathrm{2}\\ R+S,& m\ge \mathrm{3}\end{array}\right. （4）

其中

\begin{array}{l}R=\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{h}\left({\beta }_{m}x/b\right)+\mathrm{c}\mathrm{o}\mathrm{s}\left({\beta }_{m}x/b\right)\\ S={\sigma }_m\left(\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}\left({\beta }_{m}x/b\right)+\mathrm{s}\mathrm{i}\mathrm{n}\left({\beta }_{m}x/b\right)\right)\end{array} （5）

其中

{\beta }_m = \left\{\begin{array}{ll}\mathrm{4.73},& m=\mathrm{3}\\ \mathrm{7.853},& m=\mathrm{4}\\ \left(m-\mathrm{3}/\mathrm{2}\right)\mathrm{\pi },& m\ge \mathrm{5}\end{array}\right. （6）

{\sigma }_m = - \frac{\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{h}{\beta }_m-\mathrm{c}\mathrm{o}\mathrm{s}{\beta }_m}{\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}{\beta }_m-\mathrm{s}\mathrm{i}\mathrm{n}{\beta }_m} （7）

\begin{array}{l}{\phi }_n\left(y\right)=\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{h}\left({\beta }_{n}y/a\right)-\mathrm{c}\mathrm{o}\mathrm{s}\left({\beta }_{n}y/a\right)+\\ {\sigma }_n\left(\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}\left({\beta }_{n}y/a\right)-\mathrm{s}\mathrm{i}\mathrm{n}\left({\beta }_{n}y/a\right)\right)\end{array} （8）

其中

{\beta }_n = \left\{\begin{array}{ll}\mathrm{1.875},& n=\mathrm{1}\\ \mathrm{4.694},& n=\mathrm{2}\\ \left(n-\mathrm{1}/\mathrm{2}\right)\mathrm{\pi },& n\ge \mathrm{3}\end{array}\right. （9）

{\sigma }_n = - \frac{\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{h}{\beta }_n+\mathrm{c}\mathrm{o}\mathrm{s}{\beta }_n}{\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}{\beta }_n+\mathrm{s}\mathrm{i}\mathrm{n}{\beta }_n} （10）

基于式（3），悬臂板的z向位移函数 w \left(x,y,t\right)可写为

w \left(x,y,t\right) = {\displaystyle \sum _{i=\mathrm{1}}^r} {\tilde{\psi }}_i\left(x,y\right)q_i\left(t\right) （11）

式中： q_i \left(t\right)为系统第i阶主坐标；r为模态截断阶数。将式（11）代入式（1），可将系统势能改写为

U = \frac{\mathrm{1}}{\mathrm{2}} {\displaystyle \sum _{i=\mathrm{1}}^r} {\displaystyle \sum _{j=\mathrm{1}}^r}{\tilde{k}}_{ij} q_i \left(t\right) q_j \left(t\right) （12）

其中

{\tilde{k}}_{ij} = D_b {\int }_{\mathrm{0}}^a {\int }_{\mathrm{0}}^b\left[\frac{{\partial }^{\mathrm{2}}{\tilde{\psi }}_i}{\partial x^{\mathrm{2}}}\frac{{\partial }^{\mathrm{2}}{\tilde{\psi }}_j}{\partial x^{\mathrm{2}}}+\frac{{\partial }^{\mathrm{2}}{\tilde{\psi }}_i}{\partial y^{\mathrm{2}}}\frac{{\partial }^{\mathrm{2}}{\tilde{\psi }}_j}{\partial y^{\mathrm{2}}}\right] \mathrm{d} x \mathrm{d} y + \mathrm{2} D_b {\int }_{\mathrm{0}}^a {\int }_{\mathrm{0}}^b\left[\nu \frac{{\partial }^{\mathrm{2}}{\tilde{\psi }}_i}{\partial x^{\mathrm{2}}}\frac{{\partial }^{\mathrm{2}}{\tilde{\psi }}_j}{\partial y^{\mathrm{2}}}+\left(\mathrm{1}-\nu \right)\frac{{\partial }^{\mathrm{2}}{\tilde{\psi }}_i}{\partial x\partial y}\frac{{\partial }^{\mathrm{2}}{\tilde{\psi }}_j}{\partial x\partial y}\right] \mathrm{d} x \mathrm{d} y （13）

同时，将式（11）代入式（2），可将系统动能改写为

T = \frac{\mathrm{1}}{\mathrm{2}} {\displaystyle \sum _{i=\mathrm{1}}^r} {\displaystyle \sum _{j=\mathrm{1}}^r}{\tilde{m}}_{ij} {\dot{q}}_i \left(t\right) {\dot{q}}_j \left(t\right) （14）

其中

\begin{array}{l}{\tilde{m}}_{ij}=\rho h{\int }_{\mathrm{0}}^a{\int }_{\mathrm{0}}^b{\tilde{\psi }}_i\left(x,y\right){\tilde{\psi }}_j\left(x,y\right)\mathrm{d}x\mathrm{d}y+\\ {\displaystyle \sum _{p=\mathrm{1}}^{N_p}}{m}_p{\tilde{\psi }}_i\left(x_p,y_p\right){\tilde{\psi }}_j\left(x_p,y_p\right)\end{array} （15）

由式（13）和式（15）得出带内埋质量悬臂板的广义刚度矩阵 \tilde{\mathit{K}} = {\left[{\tilde{k}}_{ij}\right]}_{r\times r}和广义质量矩阵 \tilde{\mathit{M}} = {\left[{\tilde{m}}_{ij}\right]}_{r\times r}，将其代入本征方程，有

\left|\tilde{\mathit{K}}-{\omega }^{\mathrm{2}}\tilde{\mathit{M}}\right| = \mathrm{0} （16）

可得到系统的前r阶固有频率 {\omega }_{\mathrm{1}} , {\omega }_{\mathrm{2}} , \cdots , {\omega }_r和特征向量 {\mathit{a}}_{\mathrm{1}} , {\mathit{a}}_{\mathrm{2}} , \cdots , {\mathit{a}}_r，写成矩阵形式为

\mathit{\Phi } = \left[{\mathit{a}}_{\mathrm{1}},{\mathit{a}}_{\mathrm{2}},\cdots ,{\mathit{a}}_r\right] （17）

记为

\mathit{\psi } = \left[{\psi }_{\mathrm{1}},\cdots ,{\psi }_r\right] （18）

则 \mathit{\psi }与假设模态 \tilde{\mathit{\psi }} = \left[{\tilde{\psi }}_{\mathrm{1}},\cdots ,{\tilde{\psi }}_r\right]之间存在着变换关系，可写为

\mathit{\psi } = \tilde{\mathit{\psi }} \mathit{\Phi } （19）

因此，带内埋质量悬臂板的第i阶模态可写为

{\psi }_i \left(x,y\right) = {\displaystyle \sum _{j=\mathrm{1}}^r} {\tilde{\psi }}_j\left(x,y\right)a_{ji} （20）

利用式（20）可以得到带内埋质量悬臂板系统的位移函数为

w \left(x,y,t\right) = {\displaystyle \sum _{i=\mathrm{1}}^r} {\displaystyle \sum _{j=\mathrm{1}}^r}{\tilde{\psi }}_j\left(x,y\right)a_{ji}{q}_i\left(t\right) （21）

引入主质量参数为

\begin{array}{l}{M}_{ii}={\int }_{\mathrm{0}}^a{\int }_{\mathrm{0}}^b\rho h{\left[{\psi }_i\left(x,y\right)\right]}^{\mathrm{2}}\mathrm{d}x\mathrm{d}y\\ ={\int }_{\mathrm{0}}^a{\int }_{\mathrm{0}}^b\rho h{\left[{\displaystyle \sum _{j=\mathrm{1}}^r}{\tilde{\psi }}_j\left(x,y\right)a_{ji}\right]}^{\mathrm{2}}\mathrm{d}x\mathrm{d}y\end{array} （22）

主刚度参数为

\begin{array}{l}{K}_{ii}=D_b{\int }_{\mathrm{0}}^a{\int }_{\mathrm{0}}^b\left[\frac{{\partial }^{\mathrm{2}}{\psi }_i}{\partial x^{\mathrm{2}}}\frac{{\partial }^{\mathrm{2}}{\psi }_i}{\partial x^{\mathrm{2}}}+\frac{{\partial }^{\mathrm{2}}{\psi }_i}{\partial y^{\mathrm{2}}}\frac{{\partial }^{\mathrm{2}}{\psi }_i}{\partial y^{\mathrm{2}}}\right]\mathrm{d}x\mathrm{d}y+\\ \mathrm{2}{D}_b{\int }_{\mathrm{0}}^a{\int }_{\mathrm{0}}^b\left[\nu \frac{{\partial }^{\mathrm{2}}{\psi }_i}{\partial x^{\mathrm{2}}}\frac{{\partial }^{\mathrm{2}}{\psi }_i}{\partial y^{\mathrm{2}}}+\left(\mathrm{1}-\nu \right)\frac{{\partial }^{\mathrm{2}}{\psi }_i}{\partial x\partial y}\frac{{\partial }^{\mathrm{2}}{\psi }_i}{\partial x\partial y}\right]\mathrm{d}x\mathrm{d}y\end{array} （23）

由于悬臂板的刚度属性并没有因为附加的内埋质量而发生改变，所以，主刚度参数 K_{ii}也可以被改写为

\begin{array}{l}{K}_{ii}=D_b{\int }_{\mathrm{0}}^a{\int }_{\mathrm{0}}^b\left[\frac{{\partial }^{\mathrm{2}}{\tilde{\psi }}_i}{\partial x^{\mathrm{2}}}\frac{{\partial }^{\mathrm{2}}{\tilde{\psi }}_i}{\partial x^{\mathrm{2}}}+\frac{{\partial }^{\mathrm{2}}{\tilde{\psi }}_i}{\partial y^{\mathrm{2}}}\frac{{\partial }^{\mathrm{2}}{\tilde{\psi }}_i}{\partial y^{\mathrm{2}}}\right]\mathrm{d}x\mathrm{d}y+\\ \mathrm{2}{D}_b{\int }_{\mathrm{0}}^a{\int }_{\mathrm{0}}^b\left[\nu \frac{{\partial }^{\mathrm{2}}{\tilde{\psi }}_i}{\partial x^{\mathrm{2}}}\frac{{\partial }^{\mathrm{2}}{\tilde{\psi }}_i}{\partial y^{\mathrm{2}}}+\left(\mathrm{1}-\nu \right)\frac{{\partial }^{\mathrm{2}}{\tilde{\psi }}_i}{\partial x\partial y}\frac{{\partial }^{\mathrm{2}}{\tilde{\psi }}_i}{\partial x\partial y}\right]\mathrm{d}x\mathrm{d}y\end{array} （24）

则系统的第i阶固有频率为

{\omega }_i = \sqrt[]{\frac{K_{ii}}{M_{ii}}} , i = \mathrm{1,2} , \cdots , r （25）

忽略系统的结构阻尼，可以得到用模态坐标表示的系统运动方程为

M_{ii} {\ddot{q}}_i \left(t\right) + K_{ii} q_i \left(t\right) = Q_i \left(t\right) , i = \mathrm{1,2} , \cdots , r （26）

其中 Q_i为第i阶模态气动力，需要通过作用在悬臂板上的气动条件得到，式（26）的矩阵形式为

\mathit{M} \ddot{\mathit{q}} ( t ) + \mathit{K} \dot{\mathit{q}} ( t ) = \mathit{Q} ( t ) （27）

基于一阶活塞理论求解超声速条件下悬臂板上的模态气动力。由于所考虑的悬臂板为等厚度的，板的厚度函数关于顺气流方向的梯度为0，所以可以忽略板厚对上下表面压差的影响。

对于悬臂板上表面内的点，其 z向速度为

\begin{array}{l}{v}_{zU}=\left(U_{\mathrm{\infty }}\frac{\partial }{\partial x}+\frac{\partial }{\partial t}\right)w\left(x,y,t\right)\\ =\left(U_{\mathrm{\infty }}\frac{\partial }{\partial x}+\frac{\partial }{\partial t}\right)\left({\displaystyle \sum _{i=\mathrm{1}}^r}{\displaystyle \sum _{j=\mathrm{1}}^r}{\tilde{\psi }}_j\left(x,y\right)a_{ji}{q}_i(t)\right)\end{array} （28）

而对于下表面，有

\begin{array}{l}{v}_{zL}=-\left(U_{\mathrm{\infty }}\frac{\partial }{\partial x}+\frac{\partial }{\partial t}\right)w\left(x,y,t\right)\\ =-\left(U_{\mathrm{\infty }}\frac{\partial }{\partial x}+\frac{\partial }{\partial t}\right)\left({\displaystyle \sum _{i=\mathrm{1}}^r}{\displaystyle \sum _{j=\mathrm{1}}^r}{\tilde{\psi }}_j\left(x,y\right)a_{ji}{q}_i(t)\right)\end{array} （29）

利用式（28）和（29），悬臂板下表面与上表面的压差系数 \mathrm{\Delta } p \left(x,y,t\right)为

\begin{array}{l}\mathrm{\Delta }p(x,y,t)=-\frac{\mathrm{4}{q}_d}{M{a}_{\mathrm{\infty }}}\left({\displaystyle \sum _i^r}{\displaystyle \sum _j^r}\frac{\partial {\tilde{\psi }}_j\left(x,y\right)}{\partial x}{a}_{ji}{q}_i(t)\frac{\partial w}{\partial x}\right.\\ \left.+\frac{\mathrm{1}}{U_{\mathrm{\infty }}}{\displaystyle \sum _i^r}{\displaystyle \sum _j^r}{\tilde{\psi }}_j\left(x,y\right)a_{ji}{\dot{q}}_i(t)\right)\end{array} （30）

式中： q_d为动压， q_d = {\rho }_{\mathrm{\infty }} U_{\mathrm{\infty }}^{\mathrm{2}} / \mathrm{2}； M a_{\mathrm{\infty }}为未扰动空气的马赫数。将系统的位移函数带入式（30）中，得到第i阶模态气动力

\begin{array}{l}{Q}_i(t)={\int }_{\mathrm{0}}^a{\int }_{\mathrm{0}}^b\mathrm{\Delta }p\left(x,y,t\right){\psi }_i\left(x,y\right)\mathrm{d}x\mathrm{d}y\\ =-\frac{\mathrm{4}{q}_d}{M{a}_{\mathrm{\infty }}}{\displaystyle \sum _{j=\mathrm{1}}^r}\left({\overline{D}}_{ij}{q}_j(t)+\frac{\mathrm{1}}{U_{\mathrm{\infty }}}{\overline{E}}_{ij}{\dot{q}}_j(t)\right)\end{array} （31）

式中

\begin{array}{l}{\overline{D}}_{ij}={\int }_{\mathrm{0}}^a{\int }_{\mathrm{0}}^b\\ \left({\displaystyle \sum _{s=\mathrm{1}}^r}{\tilde{\psi }}_s\left(x,y\right)a_{si}\right)\left({\displaystyle \sum _{s=\mathrm{1}}^r}\frac{\partial {\tilde{\psi }}_s\left(x,y\right)}{\partial x}{a}_{sj}\right)\mathrm{d}x\mathrm{d}y\end{array} （32）

\begin{array}{l}{\overline{E}}_{ij}={\int }_{\mathrm{0}}^a{\int }_{\mathrm{0}}^b\\ \left({\displaystyle \sum _{s=\mathrm{1}}^r}{\tilde{\psi }}_s\left(x,y\right)a_{si}\right)\left({\displaystyle \sum _{s=\mathrm{1}}^r}{\tilde{\psi }}_s\left(x,y\right)a_{sj}\right)\mathrm{d}x\mathrm{d}y\end{array} （33）

从式（31）可以看出，作用在悬臂板上的气动力系数被分成了两部分：一部分为 q_j ( t )前面的系数 {\overline{D}}_{ij}，起到气动刚度的作用；另外一部分为 {\dot{q}}_j ( t )前面的系数 {\overline{E}}_{ij}，起到气动阻尼的作用。

将式（31）代入式（26）中，可得到带内埋质量悬臂板颤振方程的时域表达形式为

\begin{array}{l}{M}_{ii}{\ddot{q}}_i\left(t\right)+K_{ii}{q}_i\left(t\right)\\ =-\frac{\mathrm{4}{q}_d}{M{a}_{\mathrm{\infty }}}{\displaystyle \sum _{j=\mathrm{1}}^r}\left({\overline{D}}_{ij}{q}_j(t)+\frac{\mathrm{1}}{U_{\mathrm{\infty }}}{\overline{E}}_{ij}{\dot{q}}_j(t)\right),\\ i=\mathrm{1,2},\cdots ,r\end{array} （34）

式（34）的矩阵形式为

\begin{array}{l}\mathit{M}\ddot{\mathit{q}}(t)+\mathit{K}\dot{\mathit{q}}(t)\\ =q_d\left(-\frac{\mathrm{4}}{M{a}_{\mathrm{\infty }}}\right)\left(\overline{\mathit{D}}\mathit{q}(t)+\frac{\mathrm{1}}{U_{\mathrm{\infty }}}\mathit{E}\dot{\mathit{q}}(t)\right)\end{array} （35）

2 颤振分析

当颤振发生时，系统作谐振荡运动。p-k法将被引入以在频域内求解式（35）所示的颤振方程。为此，模态位移向量 \mathit{q} \left(t\right)被重写为

\mathit{q} \left(t\right) = \overline{\mathit{q}} e^{st} （36）

式中： \overline{\mathit{q}}是 \mathit{q} \left(t\right)的振幅，与时间无关； s = \left(\gamma +i\right) \omega， \gamma是瞬态衰减率系数，它决定了振动衰减的速度。引入无量纲参数

p = \left(\gamma +\mathrm{1}\right) k （37）

式中：k为系统缩减频率。式（36）可以被重写为

\mathit{q} \left(t\right) = \overline{\mathit{q}} e^{p\frac{U_{\mathrm{\infty }}}{b}t} （38）

将式（38）代入到式（35）中，可以得到

\left[{\left(\frac{\mathrm{2}{U}_{\mathrm{\infty }}}{b}\right)}^{\mathrm{2}}\mathit{M}{p}^{\mathrm{2}}+\mathit{K}-\mathrm{2}{q}_d\mathit{Q}\left(k,M{a}_{\mathrm{\infty }}\right)\right] \overline{\mathit{q}} = \mathrm{0} （39）

式中： \mathit{Q} \left(k,M{a}_{\mathrm{\infty }}\right) = - \frac{\mathrm{4}}{M{a}_{\mathrm{\infty }}} \left(\overline{\mathit{D}}+\frac{\mathrm{2}p}{a}\overline{\mathit{E}}\right)，为了便于计算，把 \mathit{Q} \left(k,M{a}_{\mathrm{\infty }}\right)分为实部和虚部两部分。

\mathit{Q} \left(k,M{a}_{\mathrm{\infty }}\right) = {\mathit{Q}}_{\mathrm{R}\mathrm{e}} \left(k,M{a}_{\mathrm{\infty }}\right) + \mathrm{i} {\mathit{Q}}_{\mathrm{I}\mathrm{m}} \left(k,M{a}_{\mathrm{\infty }}\right) （40）

于是颤振方程可改写为

\begin{array}{l}\left[{\left(\frac{\mathrm{2}{U}_{\mathrm{\infty }}}{a}\right)}^{\mathrm{2}}\mathit{M}{p}^{\mathrm{2}}+\mathit{K}-\mathrm{2}{q}_d{\mathit{Q}}_{\mathrm{R}\mathrm{e}}\left(k,M{a}_{\mathrm{\infty }}\right)\right.\\ \left.-\mathrm{i}\mathrm{2}{q}_d{\mathit{Q}}_{\mathrm{I}\mathrm{m}}\left(k,M{a}_{\mathrm{\infty }}\right)\right]\overline{\mathit{q}}=\mathrm{0}\end{array} （41）

将式（41）中 \mathrm{i} \mathrm{2} q_d {\mathit{Q}}_{\mathrm{I}\mathrm{m}} \left(k,M{a}_{\mathrm{\infty }}\right)改写为

\mathrm{i} \mathrm{2} q_d {\mathit{Q}}_{\mathrm{I}\mathrm{m}} \left(k,M{a}_{\mathrm{\infty }}\right) = \mathrm{2} q_d \frac{{\mathit{Q}}_{\mathrm{I}\mathrm{m}}\left(k,M{a}_{\mathrm{\infty }}\right)}{k} p （42）

显然，式（42）只有 \gamma = \mathrm{0}等号才严格成立

记

\begin{array}{l}\overline{\mathit{A}}={\left(\frac{\mathrm{2}{U}_{\mathrm{\infty }}}{b}\right)}^{\mathrm{2}}\mathit{M}\\ \overline{\mathit{B}}=-q_d\frac{{\mathit{Q}}_{\mathrm{I}\mathrm{m}}\left(k,M{a}_{\mathrm{\infty }}\right)}{k}\\ \overline{\mathit{C}}=\mathit{K}-q_d{\mathit{Q}}_{\mathrm{R}\mathrm{e}}\left(k,M{a}_{\mathrm{\infty }}\right)\end{array} （43）

可以得到带内埋质量悬臂板颤振系统的二次特征值问题为

\left[\overline{\mathit{A}}{p}^{\mathrm{2}}+\overline{\mathit{B}}p+\overline{\mathit{C}}\right] \overline{\mathit{q}} = \mathrm{0} （44）

式中：参数p的计算需要迭代进行。当系统收敛到当前流速 U_{\mathrm{\infty }}时，可以得到人工结构阻尼g和频率f的计算式为

g \approx \mathrm{2} \gamma = \mathrm{2} \frac{\mathrm{R}\mathrm{e}\left(p\right)}{k} （45）

f = \frac{k{U}_{\mathrm{\infty }}}{b\mathrm{\pi }} （46）

3 数值算例

对于超声速流中的悬臂板系统，其物理参数如表1所示。当系统在海平面附近飞行时，进行颤振计算的空气参考密度 {\rho }_{\mathrm{\infty }} = \mathrm{1.226} kg/m³，空气密度比为1.0。在类似机翼的气动弹性系统中，当系统处于颤振临界时，通常表现为气动阻尼占主导成分，而结构阻尼对气动弹性系统的影响较小。所以，为了清晰了解气动弹性系统的本质属性，在计算过程中忽略了结构阻尼，并给出偏于“保守”的设计方案。

图2为 M a_{\mathrm{\infty }} = \mathrm{2.0}时悬臂板颤振速度和颤振频率随模态截断数r的变化情况。由图2可知，为了得到收敛的颤振速度和颤振频率，至少应选取系统的前15阶固有模态参与颤振计算。为了保证计算精度，同时不增加额外的计算量，在后续的颤振计算过程中，均采用前16阶模态参与计算。

取悬臂板的前16阶模态参与计算，忽略结构阻尼，在 M a_{\mathrm{\infty }} = \mathrm{2.0}的情况下利用本文所提方法计算得到非匹配的颤振速度和颤振频率，如图3所示。从图3可知，无内埋质量件的悬臂板系统颤振速度为552.71 m/s，颤振频率为41.25 Hz，颤振发生在第2阶模态。由图3b可知，当颤振发生时，第2阶模态趋向与第1阶模态发生耦合。

为了验证本文所提出方法的准确性，使用MSC.Nastran软件计算了相同的问题，得到颤振速度为555.63 m/s，颤振频率为40.96 Hz。颤振速度和颤振频率分别相差0.53%和0.71%，计算结果误差控制在1%以内，验证了本文方法的正确性。

3.1 单个内埋集中质量的情形

首先考察在悬臂板中心点处布置内埋集中质量的情形。内埋质量件的质量 m_p为50 g（约是悬臂板质量的0.12倍），同样取前16阶模态参与颤振计算。由于有内埋质量的存在，悬臂板-质量系统的固有频率呈下降趋势。

基于本文所提方法计算了带集中质量悬臂板的颤振速度和颤振频率，如图4所示。从图4可知，悬臂板的颤振速度为611.39 m/s，颤振频率为34.97 Hz。相比无内埋质量的情形，颤振速度提升了58.68 m/s，上升比例为10.62%，而颤振频率则下降了6.28 Hz，下降比例为15.22%。颤振依然发生在第2阶模态，当颤振发生时，同样是第2阶模态与第1阶模态发生耦合，并且耦合程度变得更强。

对比无内埋质量和在悬臂板中心点内埋质量件的情形可以发现，系统的颤振特性发生了明显的变化。为进一步考察内埋质量的挂载位置对系统颤振速度的影响规律，分别在悬臂板上均匀选取156个挂载位置进行单个内埋集中质量的挂载，得到挂载位置对颤振速度和颤振频率随挂载位置变化的规律，如图5所示。

从图5a可知，在靠近翼根沿展向a/3范围内，颤振速度对内埋件位置的变化不明显，大约在1.7% 的范围内浮动，颤振也都发生在第2阶模态。在沿展向的同一站位，当内埋质量从前缘向后缘变化时，颤振速度会在某个位置发生突降，最大下降幅度可达到60%。当内埋质量沿展向变化时（中心线前的迎风区域），会出现颤振速度先增加再减小然后再增加的趋势，先后出现两次峰值。第一次峰值出现在沿展向约0.6a位置的前缘处，颤振速度为697.83 m/s，而第二次峰值出现在约0.8a位置的前缘处，颤振速度达到710.47 m/s。相较不带内埋质量的情形，颤振速度提高幅度达到28.54%。相较而言，在靠近后缘处布置内埋质量，会降低系统的颤振速度，尤其是在沿展向站位约0.7a处的后缘附近，颤振速度下降幅值达到最大。

而由图5b可知，除了将内埋质量布置在展向站位0.75 a附近的位置（靠近前缘处），颤振频率的变化对于内埋质量位置的变化整体上不敏感。在这个位置附近，颤振频率有一个突增，这是由于系统参数改变从而导致颤振模态发生了跳转。当悬臂板上的质量属性发生变化时，结构的模态属性会发生相应的变化，继而气动力也将重新分布。当质量参数（包括大小、数目和位置）变化至某个临界值时，所激发的颤振模态将会从原有模态跳转至其他模态。

由此可见，当只考虑单个内埋质量时，悬臂板面内存在两个较佳的挂载区域，即沿展向站位0.6a和0.8a附近的位置（靠近前缘处）。但值得注意的是，在这两者之间，有一个颤振速度的极小值点。

根据上述分析，位于悬臂板后缘的内埋质量会降低系统的颤振速度，所以不宜在此区域布置内埋件。基于此，仅考虑在悬臂板的前缘线L ₁、1/4弦线L ₂和悬臂板中心线位置上L ₃布置不同质量内埋件的情形，如图6所示。在不同的线上沿展向均匀选取8个位置，分别记为 P_{\mathrm{1}} ~ P_{\mathrm{8}}（不包含固定端和自由端）。

考虑到悬臂板的承载能力，所考虑内埋件的质量在10~15 g 内变化，分别得到在工况L ₁~L ₃下悬臂板的颤振特性，如图7所示。

由图7a可知，在前缘线上的不同挂点处，颤振速度会随着内埋质量的增大而出现变化。当内埋质量较小时，几乎所有情况下的颤振速度均有所增加。但随着内埋质量的增加， P_{\mathrm{2}}、 P_{\mathrm{3}}、 P_{\mathrm{5}}、 P_{\mathrm{6}}、 P_{\mathrm{7}}均先后出现拐点，而最早出现拐点的是 P_{\mathrm{6}}。不同的是， P_{\mathrm{4}}和 P_{\mathrm{8}}的颤振速度却是呈现单调递增趋势。而由图7b可知，除了 P_{\mathrm{4}}和 P_{\mathrm{8}}的颤振频率随着内埋质量的增加稍有降低外，其他挂点处的颤振频率均发生跃迁，这是由于颤振模态发生跳转所致。

由图7c可知，在悬臂板的1/4弦线上布置内埋质量，随着内埋质量的增大，除了 P_{\mathrm{5}}、 P_{\mathrm{6}}和 P_{\mathrm{7}}的颤振速度曲线有明显的拐点外，其余工况均随着内埋件质量的增大而提高或者维持较高的水平。这一点可以从图7d看出， P_{\mathrm{5}}、 P_{\mathrm{6}}和 P_{\mathrm{7}}的颤振频率曲线在内埋质量达到某个值时发生了阶跃，同样是由于颤振模态发生跳转造成的。

当内埋质量点沿弦向后移至悬臂板中心线时（工况 L_{\mathrm{3}}），图7e显示除了 P_{\mathrm{6}}，其余内埋点的颤振速度曲线均随内埋质量的增加而单调递增，尤其以 P_{\mathrm{5}}的情形最佳。由图7f所示的颤振频率曲线可知， P_{\mathrm{6}}在内埋质量为90 g 附近发生了模态跳转，而其余情形的颤振频率均随着内埋质量的增大而减小。

由以上分析可知，内埋件的挂载位置及质量大小会对系统颤振速度和颤振频率产生重要影响。适宜的挂载位置有两个区域，一个是靠近翼根的区域，在该区域挂载颤振速度和颤振频率对内埋件质量参数不敏感；另一个区域是在靠近翼梢从悬臂板中心线附近到前缘之间的位置，当在此区域挂载时，除了 P_{\mathrm{6}}附近区域，颤振速度均有大幅提升。然而，当在悬臂板展向中心线后缘方向挂载时，颤振速度则通常呈现下降趋势。

3.2 细长体的情形

当内埋件表现为细长体（沿弦向分布）时，通常使用2~3个挂点与悬臂板面连接。一般将细长体的挂点取为3个，并将其质量属性向该3个挂点进行简化。取细长体总质量为100 g，提出了6个质量分配方案，分别为S ₁~S ₆，如表2所示。

根据3.1节的分析可知，为了保证挂载内埋质量件后，颤振速度不低于原有水平，应尽可能选择悬臂板中心线前的迎风区域进行挂载。所以，选取内埋质量的挂载点与图6所示的位置类似，即将细长体的3个挂点分别设置为 P_i , ( i = \mathrm{1,2} , \cdots , \mathrm{8} )情形下的3个点，分别为前缘点、1/4弦线点和中心线点。需要注意的是3.1节中的情形是只有内埋质量点，而这里是沿弦向同时有3个简化的集中质量点。

将细长体所在位置由翼根逐步移至翼梢区域，系统的颤振特性发生了明显的变化，6种质量分配方案下的颤振速度和颤振频率曲线如图8所示。由图8可知，随着细长体挂载位置的展向变化，所有的颤振速度图像均呈现出两次峰值，即在 P_{\mathrm{4}}或 P_{\mathrm{5}}及 P_{\mathrm{7}}的位置，而 P_{\mathrm{4}}则为颤振速度曲线的陡降点，在此区域不宜布置内埋质量。不同质量分配方案下悬臂板的颤振速度和颤振频率曲线均具有类似的变化趋势，颤振速度的最大差异控制在7.0%左右。

3.3 多个内埋集中质量的情形（展向布置）

对于沿展向布置多个内埋集中质量的情形，一般考虑沿悬臂板展向均匀选取8个点进行挂载，如图6中沿 L_i ( i = \mathrm{1,2} , \mathrm{3} )的8个点所示。内埋件总质量控制在100 g以内，均匀分布在8个点上。以悬臂板中心线的情形为例，记为工况ML ₃，绘制了 U_{\mathrm{\infty }} - g图和 U_{\mathrm{\infty }} - f图，如图9a、9b所示。

由图9a、9b可知，悬臂板的颤振速度为609.47 m/s，相应的颤振频率为32.03 Hz，颤振发生在第2阶模态，为第1阶模态和第2阶模态耦合所致。相比未带内埋质量的情形，颤振速度提升了10.27%，而颤振频率则下降了22.35%。

然而，与此情形不同的是，当内埋件位置向前缘靠近时，颤振模态可能会发生跳转。考虑工况ML ₁情形，绘制了 U_{\mathrm{\infty }} - g图和 U_{\mathrm{\infty }} - f图，如图9c、9d所示。由图9c、9d可知，悬臂板颤振速度为695.17 m/s，颤振频率为19.51 Hz。颤振发生在第1阶弯曲模态，为第2阶模态与第1阶模态耦合所致，但同时第3阶模态和第4阶模态之间出现耦合趋势。

在悬臂板1/4弦向布置内埋件的情形，即工况ML ₂颤振图像与图9类似，得到的颤振速度为673.34 m/s，颤振频率为24.13 Hz，颤振发生在第1阶模态，依然是第2阶模态和第1阶模态耦合所致。

由此可以看出，当沿展向布置的多个内埋件逐渐向前缘移动时，悬臂板的颤振速度会有所提高，但却伴随着颤振频率的下降。期间，颤振模态会发生迁移，从第2阶扭转模态跳转至第1阶弯曲模态。

4 结论

基于假设模态法及一阶活塞理论得到带内埋质量超声速悬臂板在模态坐标下的动力学方程，引入p-k法在不同内埋工况下进行颤振计算，得到以下主要结论：

（1）内埋质量的存在会影响到悬臂板模态属性，并进而影响到悬臂板的颤振特性。随着内埋件质量、位置和数目的变化，悬臂板的颤振速度和颤振频率会在相当大的范围内发生变化，颤振模态也会发生跳转。

（2）在仅考虑单个内埋件时，悬臂板上存在两处较佳的挂载区域。一处是靠近翼根区域，颤振速度和颤振频率对内埋件质量不敏感；另一处是靠近翼梢区域的前缘位置，在此区域挂载可有效地提高悬臂板颤振速度。

（3）对于细长体，不同质量分配方案下的颤振速度和颤振频率曲线均具有类似的变化趋势，颤振速度的最大差异控制在7.0%左右。

（4）对于沿展向布置多个集中质量的情形，内埋件位置前移会导致悬臂板颤振速度的升高和颤振频率的下降，同时，颤振模态也会发生跳转。

参考文献

[1] Chiba M， Sugimoto T.Vibration characteristics of a cantilever plate with attached spring-mass system［J］.Journal of Sound and Vibration，2003，260（2）：237-263.

[2] 赵永辉.气动弹性力学与控制［M］.北京： 科学出版社，2007.

[3] 赵永辉，黄锐.高等气动弹性力学与控制［M］.北京： 科学出版社，2015.

[4] 肖艳平，杨翊仁.复合材料机翼/外挂系统颤振特性分析［J］.应用力学学报，2017，34（3）：502-506，611.

[5] 陈晓明，冯志华，黄萌萌.附加集中质量悬臂板的有限元分析及实验验证［J］.机床与液压，2018，46（7）：14-17.

[6] Aksencer T， Aydogdu M.Vibration of a rotating composite beam with an attached point mass［J］.Composite Structures，2018，190：1-9.

[7] Rahmane A， Benmansour T， Bouakba M，et al.Factors influencing the effect of attached mass on dynamic behavior of laminate composite plates using Taguchi technique［J］.Journal of the Brazilian Society of Mechanical Sciences and Engineering，2018，40：66.

[8] Pacheco D R Q， Marques F D， Ferreira A J M.Panel flutter suppression with nonlinear energy sinks： numerical modeling and analysis［J］.International Journal of Non-Linear Mechanics，2018，106：108-114.

[9] Molina J E， Nunez P G F.Influence of mass matrix models on flutter computations［J］.Journal of Aerospace Engineering，2018，31（3）：04018020.

[10] Bera K K， Chandiramani N K.Aeroelastic flutter control of a bridge using rotating mass dampers and winglets［J］.Journal of Vibration and Control，2020，26（23/24）：2185-2192.

[11] Zhou J， Xu M L， Xia W.Passive suppression of panel flutter using a nonlinear energy sink［J］.International Journal of Aerospace Engineering，2020，2020：8896744.

[12] Zhou J， Xu M L， Zha J，et al.The suppression of nonlinear panel flutter response at elevated temperatures using a nonlinear energy sink［J］.Meccanica，2021，56（1）：41-57.

[13] Tian W， Zhao T， GU Y，et al.Supersonic flutter control and optimization of metamaterial plate［J］.Chinese Journal of Aeronautics，2021，34（11）：15-20.

[14] Bahrami-Torabi H， Kerboua Y， Lakis A A.Finite element model to investigate the dynamic instability of rectangular plates subjected to supersonic airflow［J］.Journal of Fluids and Structures，2021，103：103267.

[15] 祁武超，张贺铭，田素梅.外挂物对大展弦比直机翼颤振特性的影响［J］.计算力学学报，2022，39（5）：598-607.

[16] Qi W C， Wang M， Tian S M.Effects of mass attachments on flutter characteristics of thin-walled panels［J］.Aerospace，2022，9（12）：748-748.

[17] Tian S M， Wang M， Qi W C.Effects of elastically supported boundaries on flutter characteristics of thin-walled panels［J］.Energies，2022，15（19）：7088.

[18] Muñoz Á， GarcíaFogeda P.Active flutter suppre-ssion of a wing section in a compressible flow［J］.Aerospace，2022，9（12）：804.

[19] Uymaz B.Vibration of CFFF functionally graded plates with an attached point mass at an arbitrary point in thermal environment［J］.Journal of Vibration Engineering & Technologies，2023，11（4）：1867-1888.

[20] Kalosha J I， Zuyev A L.Asymptotic stabilization of a flexible beam with attached mass［J］.Ukrainian Mathematical Journal，2022，73（10）：1537-1550.

[21] Dadoulis G I， Manolis G D.A note on analytical solutions for vibrations of beams with an attached large mass［J］.Archive of Applied Mechanics，2022，92（7）：1973-1982.

[22] Mahmoud M A.Free vibrations of tapered and stepped，axially functionally graded beams with any number of attached masses［J］.Engineering Structures，2022，267：114696.

[23] Simiriotis N， Palacios R.A numerical investigation on direct and data-driven flutter prediction me-

[24] thods［J］.Journal of Fluids and Structures，2023，117：103835.

基金资助

国家自然科学基金(12002218)

辽宁省教育厅项目(JYT2020034)