变体机翼粘弹性材料连接蒙皮板结构的振动特性

Vibration Characteristics of a Morphing Wings Viscoelastic Material Connected Skin Plate Structure

摘要: 本文针对新能源变体飞机机翼后缘蒙皮拼接结构为研究对象进行了系统的研究，建立柔性蒙皮单元模型及其动力学方程。采用粘弹性材料与复合材料板连接模拟蒙皮柔性部分，将粘弹性材料部分等效为各向同性材料板。基于瑞利里兹法建立了柔性蒙皮单元不同铺层角度下的动力学方程。在自由边界处引入弹性弹簧模拟约束条件，得到了系统的耦合行为与动能势能关系，并探究了耦合角度变化对结构振动特性的影响。结果表明：本文所提出的方法与有限元仿真高度吻合，正确性得到了验证。铺层角度的改变，对固有频率以及振型具有显著影响，选取适当的铺层方向可以满足特定的工况。在一定范围内增加耦合角度，可以使结构整体固有频率明显上升。本文的研究为新能源飞机的机翼后缘变形提供了技术支持。

Abstract: In this paper, a systematic study is carried out for the new energy morphing aircraft wing trailing edge skin connecting structure as the research object, and the flexible skin unit model and its dynamic equations are established. The viscoelastic material is connected with the composite material plate to simulate the flexible part of the skin, and the viscoelastic material part is equivalent to the isotropic material plate. Based on the Rayleigh-Ritz method, the kinetic equations of the flexible skin unit under different layup angles are established. Elastic spring simulation constraints are introduced at the free boundary, and the coupling behavior and kinetic potential energy relationship of the system are obtained, and the effect of the change of coupling angle on the vibration characteristics of the structure is investigated. The results show that the method proposed in this paper is highly consistent with the finite element simulation, and the correctness of the method is verified. The change of the layup angle has a significant effect on the intrinsic frequency as well as the vibration pattern, and the appropriate layup direction can be selected to meet the specific working conditions. Increasing the coupling angle within a certain range can make the overall intrinsic frequency of the structure rise significantly. The research in this paper provides technical support for the deformation of the trailing edge of the wing of a new energy airplane.

1. 引言

随着新能源飞行器的快速发展，传统刚性机翼无法满足航时航程的需求，为实现气动性能优化、减少阻力、提高升力、降低电力消耗从而达到提高航程等目标，通过动态调整机翼的气动外形的变体机翼技术是最有效的技术途径。

目前，关于变体机构的研究已经开展，然而如何在变形过程中保证表面的平滑性、完整性和在飞行过程中的强度问题成为了当前的首要任务。为此，众多学者对变体机翼的表面蒙皮结构展开了分析。周春华等[1]提出了一种半圆波纹型复合材料蒙皮结构，将橡胶材料与复合材料基体叠加使用，虽然实现了机翼表面的光滑变形，但复合材料在高频率的拉压过程中会不可避免地出现疲劳损伤，影响结构的安全性。Francois [2]以减少蒙皮质量为目标，提出一种复合材料蒙皮，通过风洞测试和疲劳循环测试验证结构耐久性，以预估保持结构完整性的极限载荷。Senthil等[3]为满足蒙皮面内弯度低和面外刚度高这一要求，将直线纤维弯曲，通过调整柔度比例来控制面内外刚度，结合数值模拟验证了设计的可行性。Guo [4]通过分割算法引入分段式结构，对蒙皮进行水平无干涉分割实现了蒙皮水平面内无干涉的刚性变形，实现了硬质蒙皮的滑动变形但机翼整体的连续性无法得到保证。

机翼蒙皮作为承担气动载荷的主要部件，常与发动机或气流激励发生不可避免的共振，这种周期性的影响，会导致裂纹的出现，从而对安全性造成影响，所以对机翼蒙皮进行动力学研究是十分重要的。目前大部分文章将机翼蒙皮简化为板与板拼接，板与壳拼接或梁模型进行振动特性分析。为探究耦合结构，首要的是先对单板部分分析；滕兆春等[5]使用广义微分求积法，得到多种边界条件下的变厚度矩形板的振动特性。曾军才等[6]用Fourier级数方法，构造了正交各向异性矩形板在任意边界下的自由振动模型。传统理论计算受限于自由、简支和固支等经典边界，为了解决这一问题，柴玉阳[7]提出在矩形板的边界分布弹簧，以弹簧刚度的改变来模拟不同的边界条件的计算方法。随着上述研究的日趋完善，耦合结构的振动特性研究逐渐被关注。Du等[8]基于瑞利–里兹法得到矩形板的振动方程，通过在两板的接触边放置弹簧来达到相互连接的作用，计算结果与有限元模型相吻合。Shi [9]提出了一种串联求解方法，实现了任意角度的T形板结构的自由振动分析。

在考虑到机翼在变形过程中，连接处会不可避免地出现应力集中。由于粘弹性材料作为柔性材料具有良好的拉伸性能可以很好地解决这一问题。目前作为夹层被广泛应用于板和梁的减振，并提出了很多分析方法。Kiasat等[10]采用玻尔兹曼叠加积分模型准确地描述了粘弹性材料在夹层中的动态力学行为。Huang [11]使用Biot模型参数法，得到了具有粘弹性材料夹层的板结构的振动特性。Jafari等人[12]提出在计算粘弹性材料固有频率时，将其等效为具有相同力学特性的各向同性材料足以满足计算精度要求。

综上所述，同种材料的相互连接和粘弹性材料的理论研究已经趋于成熟；然而，随着变体机翼后缘设计要求的提升，单一材料体系难以兼顾蒙皮的面外刚度与变形柔性需求。由于复合材料具有质量轻，强度高等优点，粘弹性材料具有较高的柔度特性，本文提出一种新的蒙皮结构设计思路；构建复合材料层合板与粘弹性材料板横向拼接的变体机翼复合材料蒙皮组合结构。此结构在保证了机翼表面高强度的同时，又实现了后缘大曲率的光滑变形。本文从理论计算与有限元仿真相结合，对复合材料板与粘弹性材料拼接的变体机翼后缘蒙皮组合结构开展研究。采用薄板理论推导出能量表达，基于瑞利里兹法得到结构的动力学方程，在自由边和连接边加入可变刚度弹簧，用于模拟结构边界条件以及复合材料与粘弹性材料复杂的耦合行为。揭示复合材料板纤维铺层角度，板间耦合角度对固有频率以及振动特性的影响规律。

2. 变体机翼复合材料蒙皮组合结构动力学建模

2.1. 结构设计

图1为蒙皮拼接单元结构示意图。板1在x₁ = a₁处与板2在x₂ = 0处连接组成，θ为两板之间的耦合夹角，复合材料板1的长度a₁ = 0.1 m，粘弹性材料板2的长度a₂ = 0.05 m，两块板具有相同的宽度b = 0.1 m和厚度h = 0.002 m。边界条件与连接边界由均布弹簧模拟，弹性边界的平动弹簧和转动弹簧的刚度分别由K_bu、K_bv、K_bw和K_b表示，连接边界的平动弹簧和转动弹簧的刚度分别由K_cu、K_cv、K_cw和K_c表示。本文将弹簧刚度的取值设置为：K_bu = K_bv = K_bw = 10¹² N/m²，K_b = 10¹² N/rad以模拟边界的固定支撑，K_cu = K_cv = K_cw = 10¹² N/m²，K_c = 10¹² N/rad以模拟复材层合板与粘弹性材料耦合的刚性连接。

Figure 1. Structure of the skinned unit

图1. 蒙皮单元结构

2.2. 能量表达

耦合板系统的拉格朗日函数可以表示为：

$L=T-V$ (1)

式中：T代表系统的总动能，V代表系统的总势能；在薄板理论的假设下，V与T可以表示为：

$\{\begin{array}{l}V=V_{1-in}+V_{2-in}+V_{1B}+V_{2B}+V_C\\ T=T_1+T_2\end{array}$ (2)

其中$V_{i-in}\left(i=1,2\right)$ 代表了第i块板的总势能，$V_{iB}\left(i=1,2\right)$ 代表了第i块板弹性边界的总势能，V_C代表了耦合边界的总势能；$T_i\left(i=1,2\right)$ 代表了第i块板的总动能。

以板1为例，总势能与总动能的表达式为：

$\{\begin{array}{l}{T}_1=\frac{1}{2}{\displaystyle {\int }_V\rho \left[{\left(\frac{\partial u_0}{\partial t}\right)}^2+{\left(\frac{\partial v_0}{\partial t}\right)}^2+{\left(\frac{\partial w_0}{\partial t}\right)}^2\right]\text{d}V}\\ V_{1-in}=\frac{1}{2}{\displaystyle {\int }_V{\epsilon }^{\text{T}}\left[S\right]\epsilon \text{d}V}\end{array}$ (3)

弹性边界势能表达：

$V_{1B}=V_{1Bx}+V_{1By}$ (4)

$\begin{array}{c}{V}_{1Bx}=\frac{1}{2}{{\displaystyle {\int }_0^b\left[k_{ux0}{u}_0^2+k_{vx0}{v}_0^2+k_{wx0}{w}_0^2+k_{\theta x0}{\left(\frac{\text{d}w}{\text{d}x}\right)}^2\right]}|}_{x=0}\text{d}y\\ +\frac{1}{2}{{\displaystyle {\int }_0^b\left[k_{uxa}{u}_0^2+k_{vxa}{v}_0^2+k_{wxa}{w}_0^2+k_{\theta xa}{\left(\frac{\text{d}w}{\text{d}x}\right)}^2\right]}|}_{x=a}\text{d}y\end{array}$ (5)

$\begin{array}{c}{V}_{1By}=\frac{1}{2}{{\displaystyle {\int }_0^a\left[k_{uy0}{u}_0^2+k_{vy0}{v}_0^2+k_{wy0}{w}_0^2+k_{\theta y0}{\left(\frac{\text{d}w}{\text{d}y}\right)}^2\right]}|}_{y=0}\text{d}x\\ +\frac{1}{2}{{\displaystyle {\int }_0^a\left[k_{uyb}{u}_0^2+k_{vyb}{v}_0^2+k_{wyb}{w}_0^2+k_{\theta yb}{\left(\frac{\text{d}w}{\text{d}y}\right)}^2\right]}|}_{y=b}\text{d}x\end{array}$ (6)

板1与板2间的耦合势能：

$\begin{array}{c}{V}_C=\frac{1}{2}{\displaystyle {\int }_0^b[K_{cu\_12}{\left({{u_1|}_{x1=a_1}+{u_2|}_{x_2=0}\cos \theta -w_2|}_{x_2=0}\sin \theta \right)}^2}\\ +K_{cw\_12}{\left({w_1|}_{x_1=a_1}+{w_2|}_{x_2=0}\cos \theta +{u_2|}_{x_2=0}\sin \theta \right)}^2\\ +K_{cv\_12}{\left({v_1|}_{x_1=a_1}-{v_2|}_{x_2=0}\right)}^2+K_{c\_12}{\left({\partial w_1/\partial x_1|}_{x_1=a_1}-{\partial w_2/\partial x_2|}_{x_2=0}\right)}^2]\text{d}y\end{array}$ (7)

式中ρ是材料密度，w为z方向的横向位移，u和v分别为x和y方向的平面内位移。

薄板的位移场可表示为：

$u=u_0-z\frac{\partial w_0}{\partial x},v=v_0-z\frac{\partial w_0}{\partial y},w=w_0$ (8)

由此可以得到应变–位移表达式：

${\epsilon }_x=\frac{\partial u}{\partial x},{\epsilon }_y=\frac{\partial v}{\partial y},{\gamma }_{xy}=\frac{\partial v}{\partial x}+\frac{\partial u}{\partial y}$ (9)

式中ε_x和ε_y为线应变；γ_xy为切应变，然后根据第k层主方向上的本构方程以及坐标转换矩阵$T^{-1}$ 得到与材料主方向成任意角度θ的应力–应变关系：

$T^{-1}=\left[\begin{array}{ccc}{\cos }^2\theta & {\sin }^2\theta & -2\sin \theta \cos \theta \\ {\sin }^2\theta & {\cos }^2\theta & 2\sin \theta \cos \theta \\ \sin \theta \cos \theta & -\sin \theta \cos \theta & {\cos }^2\theta -{\sin }^2\theta \end{array}\right]$ (10)

${\left[\begin{array}{l}{\sigma }_x\hfill \\ {\sigma }_y\hfill \\ {\tau }_{xy}\hfill \end{array}\right]}_k={\left[\begin{array}{ccc}{\overline{Q}}_{11}& {\overline{Q}}_{12}& {\overline{Q}}_{16}\\ {\overline{Q}}_{12}& {\overline{Q}}_{22}& {\overline{Q}}_{26}\\ {\overline{Q}}_{16}& {\overline{Q}}_{26}& {\overline{Q}}_{66}\end{array}\right]}_k{\left[\begin{array}{c}{\epsilon }_x\\ {\epsilon }_y\\ {\gamma }_{xy}\end{array}\right]}_k$ (11)

式中$\overline{Q}=T^{-1}Q{\left(T^{-1}\right)}^{\text{T}}$，再通过对层数积分得到复合材料层合板的刚度矩阵S。

2.3. 位移函数

在瑞利–里兹法中，构造合适的位移容许函数起着至关重要的作用，结构位移的未知函数u₀，v₀，w₀可以写成：

$\begin{array}{l}{u}_0\left(x,y,t\right)=\xi \left(x,y\right){\text{e}}^{-i\omega t}\\ v_0\left(x,y,t\right)=\zeta \left(x,y\right){\text{e}}^{-i\omega t}\\ w_0\left(x,y,t\right)=\varsigma \left(x,y\right){\text{e}}^{-i\omega t}\end{array}$ (12)

式中ω为结构固有圆频率；$\xi \left(x,y\right)$ ，$\zeta \left(x,y\right)$ 和$\varsigma \left(x,y\right)$ 为模态形函数，可以写成如下多项式形式：

$\begin{array}{l}\xi \left(x,y\right)={\displaystyle \sum _{m=1}^{M_c}{\displaystyle \sum _{n=1}^{N_c}{a}_{mn}}}{\phi }_m^u\left(x\right){\varphi }_n^u\left(y\right)\\ \zeta \left(x,y\right)={\displaystyle \sum _{m=1}^{M_c}{\displaystyle \sum _{n=1}^{N_c}{b}_{mn}}}{\phi }_m^v\left(x\right){\varphi }_n^v\left(y\right)\\ \varsigma \left(x,y\right)={\displaystyle \sum _{m=1}^{M_c}{\displaystyle \sum _{n=1}^{N_c}{c}_{mn}}}{\phi }_m^w\left(x\right){\varphi }_n^w\left(y\right)\end{array}$ (13)

式中M_c和N_c表示截断特征多项式项数；a_mn，b_mn和c_mn代表未知的多项式系数；${\phi }_m^u\left(x\right)$ ，${\varphi }_n^u\left(y\right)$ ，${\phi }_m^v\left(x\right)$ ，${\varphi }_n^v\left(y\right)$ ，${\phi }_m^w\left(x\right)$ 和${\varphi }_n^w\left(y\right)$ 为特征多项式。然后利用GramSchmidt方法构造的多项式，从而计算弹性约束矩形板的模态形函数。然后需要确定多项式系数a_mn，b_mn和c_mn使能量函数L达到最小值：

$\frac{\partial L}{\partial a_{mn}}=0,\frac{\partial L}{\partial b_{mn}}=0,\frac{\partial L}{\partial c_{mn}}=0$ (14)

式中$m=1,2,\cdots ,M_c$ ；$n=1,2,\cdots ,N_c$ 。

将以上各式整合，得到一组关于未知系数的线性方程组，进一步矩阵化可得：

$\left(K-{\omega }^{2}M\right)E=0$ (15)

其中$E={\left[X_1,X_2\right]}^{\text{T}}$ ，X₁为板1未知系数组成的系数列向量，X₂为板2未知系数组成的系数列向量，K为刚度矩阵，M为质量矩阵；通过求解特征方程(15)，可以获得矩形板结构的固有频率及模态振型函数。

3. 变体机翼粘弹性材料连接蒙皮动力学特性

本节通过有限元仿真验证了理论计算方法的有效性，并进一步探究了铺层角度与耦合角度对振动特性的影响。复合材料层合板每层由相同的单向复合材料层(T300)对称铺设，材料参数分别为：密度ρ₁ = 1512 kg/m³，x方向杨氏模量E₁ = 149 × 10⁹ N/m²，y方向E₂ = 8 × 10⁹ N/m²，泊松比μ₁ = 0.32。粘弹性材料使用各向同性材料(Dyad609)，材料参数为：密度ρ₂ = 1250 kg/m³，弹性模量E_V = 8.27 × 10⁸ N/m²，泊松比μ₂ = 0.499。位移函数的截断阶数为M = N = 7。表1给出了耦合板在多种铺层角度情况下，理论计算的前5阶固有频率和ABAQUS仿真的计算结果的对比，可以明显地看出结果高度一致，最大误差仅仅为1.12%。

Table 1. First 5 orders of natural frequency (θ = 0˚)

表1. 前5阶固有频率(θ = 0˚)

铺层方式	阶次	理论/Hz	ABAQUS/Hz	误差/%
[0/90]5s	1	242.06	241.87	0.07
	2	369.74	367.76	0.53
	3	795.84	792.02	0.48
	4	846.58	841.44	0.60
	5	850.59	842.39	0.96
[−45/45]5s	1	208.04	207.07	0.12
	2	489.30	488.06	0.91
	3	715.04	711.73	0.51
	4	856.09	848.60	1.12
	5	861.13	851.72	0.19
[30/60]5s	1	203.25	202.88	0.18
	2	347.78	344.56	0.92
	3	701.60	699.12	0.35
	4	748.79	751.42	0.35
	5	830.91	825.24	0.68

为进一步验证方法的正确性，以及研究铺层角度对振动特性的影响，特别地在图2中分别展示了耦合角度为0˚下[0/90]_5s、[−45/45]_5s和[30/60]_5s三种铺层方式的前三阶模态振型图。图2中(a)、(c)和(e)各图为理论结算所得模态振型，(b)、(d)和(f)为有限元仿真所得模态振型。从图2中可以明显地看到，此方法与仿真所计算结果具有极高的吻合度。

(a) 理论计算[0/90]_5s铺层前三阶模态振型

(b) 有限元仿真[0/90]_5s铺层前三阶模态振型

(c) 理论计算[−45/45]_5s铺层前三阶模态振型

(d) 有限元仿真[−45/45]_5s铺层前三阶模态振型

(e) 理论计算[30/60]_5s铺层前三阶模态振型

(f) 有限元仿真[30/60]_5s铺层前三阶模态振型

Figure 2. Comparison of the first 3 orders of modal shapes of the coupling plate structure

图2. 耦合板结构前3阶模态振型对比

为研究耦合角度对结构固有频率的影响，在铺层方式为[0/90]_5s的基础上改变耦合角度，通过改变板1与板2之间的耦合角度θ，表2得到θ取值为0˚、10˚和20˚的前5阶固有频率，与有限元仿真吻合良好，最大误差仅为0.96%，进一步验证了所构建模型的准确性，可以简便地对任意角度耦合的板进行计算。为方便观察固有频率随耦合角度变化的趋势为进一步探究耦合角度变化对固有频率的影响，将0~90˚的前两阶固有频率变化趋势进行分析，如图3所示。

Table 2. First 5 orders of natural frequency (Lamination mode: [0/90]_5s)

表2. 结构的前5阶固有频率(铺层方式：[0/90]_5s)

耦合角度	阶次	理论/Hz	ABAQUS/Hz	误差/%
0˚	1	242.11	241.87	0.10
	2	369.76	367.76	0.54
	3	795.84	792.02	0.48
	4	846.56	841.44	0.60
	5	850.58	842.39	0.96
10˚	1	486.51	484.77	0.35
	2	549.86	549.08	0.14
	3	804.07	806.92	0.35
	4	861.15	859.28	0.21
	5	864.22	864.69	0.05
20˚	1	687.39	690.29	0.004
	2	718.71	717.91	0.001
	3	869.91	866.25	0.004
	4	900.65	899.9	0.0008
	5	968.34	969.01	0.0006

Figure 3. Variation of the intrinsic frequency of the structure of the skinned splicing unit with the angle of coupling

图3. 蒙皮拼接单元结构固有频率随耦合角度变化关系

从图3可以看出耦合系统的第1阶和第2阶固有频率在耦合角0˚~20˚内，有明显的上升趋势，通过对比，1阶频率变化速率高于2阶。20˚~30˚内，频率的增量逐渐减小。为仔细观察变化趋势，特别地将耦合角度为40˚~90˚部分做了放大处理，可以看出二者频率趋于稳定，总体呈非线性变化。

4. 结论

本文提出了一种新的针对新能源变体飞机机翼后缘的蒙皮结构，构建的模型可以适应任意耦合角度；并计算了固有频率和振型与有限元方法结果进行了对比，验证了方法的正确性和可靠性。得出以下主要结论：

(1) 不同铺层角度下，复合材料蒙皮组合结构的模态和频率变化明显且有限元与理论分析的结果误差较小。这可以为新能源变体飞机机翼后缘的蒙皮结构的设计提供理论支撑。

(2) 复合材料蒙皮组合结构前两阶固有频率随耦合角度变化趋势都是先快速增长后趋于稳定。在耦合角处于0˚~20˚的范围内变化，第一阶固有频率的增长速度低于第二阶固有频率。