Z型折叠板内共振下非线性振动特性研究
郭翔鹰 张杨 张伟
摘要：Z型折疊板是一类复杂多体结构，在工程实际中应用广泛。以大型的空间可展开结构为实际工程背景，考虑了Z型折叠板折叠过程中各部件的运动耦合关系，基于Reddy经典板理论和von-Karman几何非线性应变位移关系，利用Hamilton原理建立了Z型折叠板的非线性动力学控制方程。利用有限元分析得到了Z型折叠板结构的模态函数，运用二阶Galerkin截断，得到了Z型折叠板二自由度的非线性常微分方程。考虑系统主参数共振-1∶2内共振的情况，通过摄动分析得到系统四维直角坐标形式的平均方程，最后利用数值模拟方法研究了横向激励对Z型折叠板非线性动力学特性的影响。
关键词： 非线性动力学； Z型折叠板； 内共振； 数值模拟
中图分类号： O322文献标志码： A文章编号： 1004-4523（2018）02-0183-15
DOI：10.16385/j.cnki.issn.1004-4523.2018.02.001
引言
可展开（折叠）结构具有悠久的研究历史和工程应用。折叠板结构是可展结构体系中最常用的一种，近年来在航空航天和建筑领域中得到了广泛的应用。在折叠状态下，结构体积较小，可用于运输或存储；在外力作用或系统内部驱动下，结构逐步展开，最终达到完全展开的工作状态，然后锁定为稳定状态[1-2]。其中，Z型折叠板结构是此类折叠结构中最为经典的一种，在航空领域主要被应用于可变体飞行器机翼结构，而其中运动过程中由于变形引起的结构非线性动态特性问题是结构平稳运行的关键。
在Z型折叠板结构振动特性的研究方面，国内外学者进行了大量的研究，如：Lee等[3-4]分别利用有限元方法和高阶板理论研究了复合材料折叠板的振动特性。Pal等[5-8]分别使用有限元方法，高阶层合板理论研究了复合材料层合板折叠结构的自由振动特性。Topal等[9] 利用一阶板理论建立了具有对称折叠角度的复合材料折叠板的动力学模型，通过MFD方法进行频率优化，数值分析得出层合板的长厚度比、夹角、板的长度以及边界条件会影响到其频率特性。Jian等[10]建立了有限元模型，用一种条状网格划分单元，模拟了折叠板的静态特征和动力学特性分析。Liu等[11]对折叠板多体系统进行了动力学建模，并用数值方法研究了其动态特性。张伟等[12]建立了变角度Z型梁的动力学方程，计算了结构的固有频率和振动模态，并通过数值仿真和实验验证。此外，Zhang等[13-14]研究了复合材料正交铺设层合板结构在外激励作用下的非线性振动问题，分析了系统不同参数下的振动响应。浙江大学的赵孟良和关富玲等[15-16]研究了空间可展结构在外载荷作用下的运动特性。
然而，Z型折叠结构属于多体结构，在外激励作用下将产生大幅的非线性耦合振动响应，这部分的研究成果目前还很有限。本文根据实际工程背景，研究了一类Z型折叠板结构在一定频率的外激励作用下发生1∶2内共振情况下的非线性动力学特性。
1动力学模型和方程〖2〗1.1动力学模型以大型的空间可展开结构为实际工程背景，建立如图1所示的力学模型，所研究的Z型折叠板采用碳纤维复合材料层合板，折叠结构由三部分组成，靠近固定端的部分称为内板，用Ω1表示，完成翻转动作的称为中间板，用Ω2表示，中间板外侧部分为外板，用Ω3表示。图中O1x1y1，O2x2y2，O3x3y3分别为内板、中间板、外板的局部坐标系，全局坐标系Oxy与局部坐标系O1x1y1重合，各个板表面受到的横向简谐激励表示为P1，P2和P3，如图1所示。
模型中，各个部分之间通过刚性铰链相连接，通过内部的作动器和机械结构进行折叠与展开动作，机构在变形时，内板与固定端相连接，中间板在一定的角度范围内转动，外板始终与内板保持平行。
本文研究过程中，满足以下条件：
1）Z型折叠板3块板是等厚度，等宽度的；
2）Z型折叠板3块板在外激励作用下不会产生材料本身的破坏。
图1Z型折叠板的动力学模型
Fig.1Mechanical model of the Z-type folding plates第2期郭翔鹰，等： Z型折叠板内共振下非线性振动特性研究振 动 工 程 学 报第31卷1.2动力学方程的建立
本文选取折叠板材料为正交铺设的碳纤维复合材料层合板，不考虑剪切效应，因此（σz ）z=+ /-h2=0，（τxz）z=+/-h2=0，（τyz）z=+/-h2=0，其中h为层合板结构厚度。
材料的本构关系如下
σxx
σyy
τyz
τzx
τxy=1112000
2122000
004400
000550
000066εxx
εyy
γyz
γzx
γxy（1）
式中ij为转换弹性常数，定义为ij=T-1QijT-1T（2）并且T-1=cos2θsin2θ0
sin2θ〖〗cos2θ0
00cos2θ-sin2θ（3）正交铺设复合材料层合板的弹性模量Qij（i=1，2，4，5，6； j=1，2，4，5，6）可以表示为如下形式：Q11=Q22=E1-ν2，Q12=Q21=νE1-ν2，
Q44=Q55=Q66=E2（1+ν）（4）内力和弯矩的关系可表示为：Nxxi
Nyyi
Nzzi=∑Nk=1∫zk+1zkσxxi
σyyi
σzzidz （5）
Mxxi
Myyi
Mzzi=∑Nk=1∫zk+1zkσxxi
σyyi
σzzizdz（6）本文所研究的Z型可折叠结构在横向外载荷的作用下会产生较大的变形，因此这里用von-Karman大变形理论进行分析，各个板的非线性应变表达式为εx
εy
γxy=ε0+zε1（7）式中ε0=u0ixi+12w0ixi2
v0iyi+12w0iyi2
v0ixi+u0iyi+w0ixiw0iyi（8）
ε1=-ω20ix2i
-ω20iy2i
-2ω20ixiyi （9）式中下标i表示第i个板（i=1，2，3）。
根据正交铺设碳纤维增强复合材料的结构特点，可得到内力与应变关系为N
M=AB
BDε1
ε2 （10）式中剛度矩阵Aij， Bij， Dij， 用如下形式表示Aij，Bij，Dij=∫h2-h2ij1，z，z2dz（11）根据经典的层合板理论[17]，3块板上任意一点在局部坐标系中的位移分别表示为以下形式：ui=u0i（xi，yi，zi，t）-ziw0ixi（12a）
vi=v0i（xi，yi，zi，t）-ziw0iyi（12b）
wi=w0i（xi，yi，zi，t），（i=1，2，3）（12c）式中u0i，v0i，w0i是在板Ωi中性层上任意一点在X， Y， Z方向的位移。
根据本文模型中所建立的坐标转换关系，Z型折叠板结构在全局坐标系下的位移场可描述为：r1=R1u1
v1
w1 （13a）
r2 = r1 L1 cosw1′（L1 ）
v2
w2 + R2 u2
v2
w2 （13b）
r3 = r2 L2 cosw2′（L2 ）
v2
w2 （L2 ） + u3
v3
w3 （13c）式中Li为各个板的长度， Ri 为坐标转换矩阵，形式如下R1=1
1
1 （14a）
R2=cosθ0-sinθ
0〖〗10
sinθ0cosθ （14b）将上述表达式（1）～（14）代入Hamilton原理中，整理展开后得到用广义位移表示的Z型折叠板结构的非线性动力学方程如下。
内板的非线性动力学方程为：
A11（2u01x21+w01x12w01x21）+A12（2u01x21+
w01x12w01x21）+A66（v01x1+u01y1+w01x1w01y1）=
I02u01t2-I13w01t2x1（15a）
A12（2v01y21+w01y12w01y21）+A22（2v01y21+
w01y12w01y21）+A66（v01x1+u01y1+w01x1w01y1）=
I02v01t2-I1（3w01〖〗t2y1）（15b）
A11[2u01x21w01x1+u01x12w01x21+3〖〗2w01x122w01x21]+
A12[2u01x21w01x1+u01x12w01x21+32w01x122w01x21]+
A12[2v01y21w01y1+v01y12w01y21+32w01y122w01y21]+
A22[2v01y21w01y1+v01y12w01y21+32w01y122w01y21]-
（D11+D12）4w01x41-（D22+D12）4w01y41-
8D664ω01x21y21+A66[2v01x12ω01x1y1+2u01y12ω01x1y1+
4w01x1w01y12ω01x1y1+v201x1y1w01x1+2u01x1y1w01x1+ 2v01x21w01y1+2u01x1y1w01y1+2w01x21w01y12+
2w01y21w01x12]+P1=I02w01t2+I1（3u01t2x1+
3v01t2y1）-I2（3w01t2x1+3w01t2y1）（15c）
中间板的非线性动力学方程为：
A11（2u02x22+w02x22w02x22）+A12（2u02x22+
w02x22w02x22）+A66（v02x2+u02y2+w02x2w02y2）=
I0（sinkθ+coskθ）2u02t2+I0（sinkθ+
coskθ）（2w01（L1）tx1+）u02-I1（sinkθ+
coskθ）3w01t2x2-I1（sinkθ+coskθ）A12（2v02〖〗y22+
w02y22w02y22）+A22（2v02y22+w02y22w02y22）+A66（v02x2+
u02y2+w02x2w02y2）=I02v02t2-I1（3w02t2y2）（15d）
A12（2v02y22+w02y22w02y22）+A22（2v02y22+
w02y22w02y22）+A66（v02〖〗x2+u02y2+w02x2w02y2）=
I02v02t2-I1（3w02t2y2） （15e）
A11[2u02x22w02x2+u02x22w02x22+w02x222w02x22]+
A12[2u02x22w02x2+u02〖〗x22w02x22+
32w02x222w02x22]+A12[2v02y22w02y2+
v02y22w02y22+32w02y222w02y22]+A22[2v02y22w02y2+
v02y22w02y22+32w02y222w02y22]-（D11+
D12）4w02x42-（D22+D12）4w02y42-8D664ω02x22y22+
A66[2v02x22ω02x2y2+2u02y22ω02x2y2+
4w02x2w02y22ω02x2y2+v202x2y2w02〖〗x2+
2u02x2y2w02x2+2v02x22w02y2+2u02x2y2w02y2+
2w02x22w02y22+2w02y22w02x22]+P2=
-I02sinkθ2w02t2+I02coskθ（2w01（L1）tx1+
）w02+I1（sinkθ+coskθ）3u02t2x2+I1（sinkθ+
coskθ）（2w01（L1）tx1+）u02x2+I1+3w02t2y2-
I2（sinkθ+coskθ）4w02t2x22+I2（sinkθ+
coskθ）（2w01（L1）tx1+）w202x22-I24w02t2y22 （15f）
外板的非线性动力学方程为：
A11（2u03x23+w03x32w03x23）+A12（2u03x23+
w03〖〗x32w03x23）+A66（v03x3+u03y3+w03〖〗x3w03y3）=
I02u03t2-I13w03t2x3+（-θ）（15g）
A12（2v03y23+w03y32w03〖〗y23）+A22（2v03y23+
w03y32w03y23）+A66（v03x3+u03y3+w03x3w03y3）=
I02v03t2-I1（3w03t2y3）（15h）
A11[2u03x23w03x3+u03x32w03x23+32w03x322w03x23]+
A12[2u03x23w03x3+u03x32w03x23+32w03x322w03x23]+
A12[2v03y23w03y3+v03y32w03y23+32w03y322w03y23]+
A22[2v03y23w03y3+v03y32w03y23+32w03y322w03y23]-
（D11+D12）4w03x43-（D22+D12）4w03y43-D664ω03x23y23+
A66[2v03x32ω03x3y3+2u03y32ω03x3y3+
4w03x3w03y32ω03x3y3+v203x3y3w03x3+2u03x3y3w03x3+
2v03x23w03y3+2u03x3y3w03y3+2w03x23w03y32+
2w03y23w03x32]+P3=I02w03t2+I13u03t2x3+
I13u03t2y3-I24w03t2x23-I24w03t2y23（15i）
結构整体的边界条件满足如下等式：u1（0）=0，u1 x1 = L1 = u2 （0） = L1 （16a）
u2 x2 = L2 = u3 （0） = L1 + L2 cosθ（16b）
u3 x3 = L3 = L1 + L2 cosθ + L3 （16c）
w1（0）=0，w1（L1）=w2（0）=0（16d）
w3 （0） = w2 （L2 ） = L2 sinθ（16e）
w1′（0）=0，w2″（0）=0，w3（0）=0（16f）
w01 （0）y1 x1 = 0 = 0 （16g）
2w2（0）y22x1 = L1 = 2w3（0）y23x2 = L2 = 0（16h）
w01 （b）y1 x1 = 0 = 0（16i）
2w2（b）y22x1 = L1 = 2w3（b）y23x2 = L2 = 0（16j）
2w03 tx3 x3 = L3 = L1 + L2 cosθ + L3 （16k）
（A12+A22）[v0iyi+12（w0iyi）2]y=0，b=0 （16l）
2w02 tx2 x2 = L2 = L1 + L2 cosθ （16m）2有限元模态分析
上述所建立的动力学控制方程是偏微分方程，利用数学方法直接求解极为困难，因此本文将采用Galerkin离散将偏微分方程转换到常微分方程后进行求解。对于单一结构模型，其模态函数通常是可以直接根据经验公式假设的，但Z型折叠板结构为多体结构，很难确定折叠板在外激振力作用下结构的振动模态。因此首先通过ANSYS有限元方法，对Z型折叠板进行模态分析和谐响应分析，得到Z型折叠板结构的固有频率和模态。通过研究结构模态振型，确定系统的模态函数形式。最后，将无量纲形式的Z型折叠板结构的非线性动力学控制方程通过Galerkin方法进行二阶离散，得到可求解的常微分方程组。
2.1有限元模型
航空领域实际应用中，Z型折叠板结构在折叠展开运动过程中的角速度很小，且折叠角度在0～150°的范围内，因此本文将板折叠过程分解，看作是不同折叠角度板结构的准静态慢变过程，选取几个特定折叠角度来建立Z型碳纤维复合材料层合板结构的力学模型。
本章研究中取折叠角度为60°，90°和120°作为典型参数值，进行对比分析和讨论。
有限元模型统一采用四边形板单元，如图2所示；设置弹性模量为5×105 MPa，泊松比为0.3。Z型折叠板的有限元模型尺寸参数如表1所示。
图2模型网格划分示意图
Fig.2Meshing of the Z-type folding plate wing of the three angles
表1Z型折叠板有限元模型的尺寸参数
Tab.1Geometric parameters of the finite element model for the Z-type folding wing
板长度/m厚度/m宽度/m内板20.012.1中间板10.012.1外板40.012.12.2模态分析
通过对建立的有限元模型施加初始条件，进行模态分析，得到Z型折叠板结构不同折叠角度下前5阶固有频率如表2所示。
表2不同折叠角度前5阶固有频率（单位：Hz）
Tab.2The first five order natural frequencies of different folding angles （Unit：Hz）
阶数60°90°120°10.560180.617540.6944322.711402.358002.6580033.225903.319903.6326043.550903.731604.2836057.367306.220905.96140
下面列出Z型折叠板结构在60°，90°，120°的折叠角度下的模态振型图，如图3～5所示。
由有限元模态分析结果可以看出，不同折叠角度下的Z型折叠板结构前5阶模态可表示为： 第1图3折叠角度为60°时Z型折叠板前5阶模态振型图
Fig.3The first five order mode shapes of the Z-type folding angle 60°图4折叠角度为90°时Z型折叠板前5阶模态振型图
Fig.4The first five order mode shapes of the Z-type folding angle 90°圖5折叠角度为120°时Z型折叠板前5阶模态振型图
Fig.5The first five order mode shapes of the Z-type folding angle 120°
阶为弯曲振动，第2阶为扭转振动，第3阶为弯曲振动，第4阶为弯扭耦合振动，第5阶为弯曲振动。通过以上分析可以发现，Z型折叠板结构的前5阶振动模态的形式与悬臂板结构前5阶振动模态的形式相似。
3Galerkin离散
由于结构在共振情况下会发生剧烈的大幅振动，容易产生失稳及结构整体破坏等。因此，对于结构在共振情况下动力学响应的研究是十分必要的。
在有限元分析结果中可以发现，Z型折叠板结构在折叠角度为60°时，第4阶固有频率几乎是第5阶固有频率1/2，存在1∶2内共振的情况，因此，本文将对Z型折叠板结构在主参数共振-1∶2内共振的情况下的动力学特性进行深入分析。
根据上述的分析结果，可选取悬臂板结构的振动模态函数形式作为Z型折叠板在折叠角度为60°时的模态函数，利用Galerkin方法对方程（15）进行二阶截断，得到系统常微分形式的非线性动力学方程。在满足位移边界条件的前提下，选取3个方向的基本函数分别如下：
Ui（x，y，t） = u1 （t）sinx/2πLi cosπyb +
u2 （t）sin3x/2πLi cos2πyb（17a）
Vi （x，y，t） = v1 （t）cosx/2πLi sinπyb +
v2 （t）cos3x/2πLi sin2πyb（17b）
Wi（x，y，t）=w1（t）Xi1（x）Yi1（y）+
w2（t）Xi2（x）Yi2（y） （17c）
根据4阶常微分方程的通解，可将方程（17c）中的Xij（x）和Yi1（y）取为Xij （x） = Ai1 coshki x-Ai2 coski x-
Ai3 βi1 sinhki x + Ai4 sinki x，
Yij （x） = Bi1 coshki x-Bi2 coski x-
Bi3 βi1 sinhki x + Bi4 sinki x（18）式中Xij（x）是沿x轴方向的固支-自由梁函数，Yij（y）为y轴方向的自由-自由梁函数，ki1和ki2为特征方程的根，并有如下关系coski1coshki1+1=0
coski2coshki2-1=0（19）将系统实际参数L1=2 m，L2=1 m，L3=4 m代入式（17）可求得X，Y方向的模态函数，将方程（18），（19）代入边界条件求得模态参数Aij和Bij，得到Xij（x）和Yij（y）的函数表达式，如下：
X11 =-4.9×10 - 5coshk11 x +
2.1×10 - 5cosk11 x - 6.6×10 - 5β11 ·
sinhk11 x-sink11 x （20a）
X12 = 0.25coshk12 x + 27.3cosk12 x -
0.25β12 sinhk12 x - sink12 x（20b）
X21 = -1.4×10 - 6coshk21 x +
1.3×10 - 6cosk21 x - 2.5×10 - 7β21 ·
sinhk21 x-sink21 x（20c）
X22 = 1.2×10 - 4coshk22 x-1.4cosk22 x -
0.12β22 sinhk22 x-sink22 x（20d）
X31 = -1.5×10 - 4coshk31 x +
3.4×10 - 4cosk31 x - 3×10 - 4β31 ·
sinhk31 x-sink31 x（20e）
X32 =0.25coshk32 x + 27.3cosk32 x -
0.25β32 sinhk32 x - sink32 x （20f）
同时，将外激励也进行离散，表示为
Pi=Fi1sin3πxlisinπyb+Fi2sinπxlisin3πyb （21）
用二阶Galerkin方法离散方程（15），并将离散后的面内位移u0，v0用横向位移w0表示，整理可得到Z型折叠板横向振动的常微分运动控制方程：
内板方程：
11+μ111+21w11+α112w12+α113w311+α114w312+ α115w12w211+α116w11w212=f11cos（ξt）（22a）
12+μ212+α121w11+22w12+α123w311+α123w312+α125w12w211+α126w11w212=f12cos（ξt）（22b）
中间板方程：
21+μ321+α21121w21+α21222w21+α21321w22+α21422w22+21w21+α216w22+α217w21u21+
α218w21u22+α219w22u21+α2110w22u22+
（α2111/θ）w321+（α2112/θ）w322+α2113（sinθ+
cosθ）w22w221+α2114（sinθ+cosθ）w21w222=
f21cos（ξt） （22c）
22+μ422+α22121w21+α22222w21+α22321w22+α22422w22+α225w21+22w22+α227w21u21+
α228w21u22+α229w22u21+α2210w22u22+w321+（α2212/θ）w322+α2213（sinθ+cosθ）w22w221+
α2214（sinθ+cosθ）w21w222=f22cos（ξt） （22d）
外板方程：
31+μ531+21w31+（α311+F1cosΩt）ω31+
α312w32+α313w331+α314w332+α315w32w231+
α316w31w232=f31cos（ξt）（22e）
32+μ632+22w32+（α322+F2cosΩt）ω32+
α321w31+α323w331+α324w332+α325w32w231+
α326w31w232=f32cos（ξt）（22f）
式中折叠角θ为3块板的运动方程的连接参数，体现了3块板之间的耦合运动关系。
4摄动分析
对于较复杂的非线性常微分方程，很难求出其精确解，需要用近似解析的方法求其渐近解来替代精确解。摄动分析是近似解析的一种方法，包括直接摄动法、多尺度法以及KBM法等。
因此，本章基于系统主参数共振-1∶2内共振的共振关系，使用多尺度方法进行摄动分析。系统共振关系表示如下：21=14Ω2+εσ1， 22=Ω2+εσ2（23）式中1，2是相应线性系统的第1阶、第2阶固有频率。σ1，σ2为系统的调谐参数，为了方便分析，令Ω=1。
经过计算得到系统的直角坐标下平均方程为：
内板平均方程：
11=-12u1x11-σ11x12-3α113x211x12- 3α113x312-2α116x12x213+x214（24a）
12=-12u1x12+σ11x11+3α113x11x212+x311+ 2α116x11x213+x214（24b）
13=-12u2x13-12σ12x14-32α124x213x14+x314-α125x14x211+x212（24c）
14=-12u2x14+12σ12x13+32α124x214x13+x313+α125x13x211+x212-14f12（24d）
中間板平均方程：
21=-12u3x21-σ11x22+14α212x21x24-x22x23- 3α214x22x221+x222-2α217x21x223+x224（25a）
22=-12u3x22+σ11x21-14α212x21x23-x22x24+ 3α214x21x221+x222+2α217x21x223+x224（25b）
23=-12u4x23-32α225x24x223+x224-
α226x24x221+x222-14σ2x4（25c）
24=-12u4x24+32α225x23x223+x224+
α226x23x221+x222+14σ2x3-14f22 （25d）
外板方程：
31=-12u5x31-σ1x32-3α313x231x32+x332-
2α316x32x233+x234+12f31x32（26a）
32=-12u5x32+σ1x31+3α313x31x232+x331+
2α316x31x233+x234+12f31x32（26b）
33=-12u6x33-12σ2x34-32α324x233x34+x334-α325x34x231+x232 （26c）
34=-12u6x34+12σ2x33+32α324x234x33+x333+α325x33x231+x232-12f32（26d）
5数值模拟
根据数值分析结果发现，Z型折叠板内板的振动幅值很小且多为周期性颤振[18]，考虑到本文篇幅的限制，不再详述，主要讨论Z型折叠机板在橫向激励作用下中间板和外板的非线性振动特性。
5.1幅频响应特性分析
通过数值求解系统的四维平均方程，利用matlab软件绘制3块板的幅频响应曲线，选取外激励幅值和系统阻尼系数为控制参数，研究参数对系统幅频特性的影响。
首先，根据实际参数的取值范围，经过无量纲处理后，选取参数为μ1=0.25，μ3=0.22，μ4=0.14， σ3=-0.015，σ4=-0.014，α214=0.31，α217=0.0002，α225=-0.108，α226=3.5，μ5=0.37，μ6=0.56，σ5=1.77，σ6=1.88，α313=2.49，α316=-3.27，α324=7.16，α325=6.81，α313=2.49。
将2块板的初始条件均设为x10=1.44，x20=1.55，x30=1.35，x40=-1.799。令外激励幅值fi （i=2，3）的值分别为50和100，研究结构幅频响应曲线的变化。图中蓝色曲线和红色曲线分别表示结构第4阶模态和第5阶模态的幅频响应曲线，横坐标为调频参数σi （i=2，3），纵坐标为振动幅值ai （i=1，2）， a1， a2 分别表示第4阶和第5阶的振动幅值。
由图6和7可知，随着外激励幅值的增加，幅频
图6不同外激励幅值下的中间板Ω2幅频响应曲线
Fig.6The frequency-response curves of the middle plate Ω2 to the external excitation amplitude f2
图7不同外激励幅值下的外板Ω3幅频响应曲线
Fig.7The frequency-response curves of the outer plate Ω3 to the external excitation amplitude f3
响应曲线的形态发生了不规律变化。随着调频参数的增加，系统振幅除中间板第5阶模态以外都呈现减小的趋势，并且会出现多值现象和跳跃现象。
为研究系统阻尼系数对幅频响应曲线的影响，分别令阻尼系数μi （i=3，5）的值为0.3和0.8。绘制如下幅频响应曲线。
观察图8和9发现，随着阻尼的增加，系统的振动幅值降低，幅频响应曲线形态发生变化，且对第4阶模态的频响特性影响较大。随着调频参数σ的增加，系统会出现多值和跳跃的现象。
图8不同阻尼下的中间板Ω2幅频响应曲线
Fig.8The frequency-response curves of the middle plate Ω2 to the damping coefficient μ3
图9不同阻尼下的外板Ω3幅频响应曲线
Fig.9The frequency-response curves of the outer plate Ω3 to the damping coefficient μ5
5.2非线性振动响应分析
为了研究Z型折叠板系统在主参数共振-1∶2内共振的情况下的非线性振动响应特性，选取外激励幅值fi为控制参数，研究横向激励幅值对系统产生周期运动和混沌运动的影响。
固定上述结构参数，改变外激励的幅值，运用四阶龙格库塔法对系统运动方程进行数值求解，得到中间板和外板的混沌分叉图，如图10和11所示。
图10中间板Ω2随外激励幅值变化的分叉图
Fig.10The bifurcation diagram of the middle plate Ω2 with the transverse excitation f2
图11外板随外激励幅值变化的分叉图
Fig.11The bifurcation diagram of the outer plate Ω3 with the transverse excitation f3
由图中可以看出，在外激励幅值增大的过程中，结构会出现周期运动-混沌运动-周期运动-混沌运动的变化，说明共振情况下，外激励幅值的变化对于系统的运动稳定性具有重要的影响。
为了更好地描述上述分叉图的特性，首先对中间板随外激励幅值变化的振动特性给出了具体的分析：图12，13和14分别给出了中间板在外激励幅值为3，10和35时的波形图、三维相图和Poincaré截面，此时中间板的振动是从单倍周期进入短暂的混沌运动之后又变为倍周期运动。继续增大外激励幅值到，中间板的振动形式趋于复杂，逐渐变为如图15所示的概周期运动，最后进入图16所示的混沌运动，且不会伴随外激励幅值继续增大而改变运动形态。
图12f2=3时中间板的周期运动
Fig.12The periodic motion of the middle plate when f2=3图13f2=10时中间板的混沌运动
Fig.13The chaotic motion of the middle plate when f2=10
外板随外激励幅值变化的振动特性如图17～20所示。当外激励幅值小于10时，外板呈现不规律的混沌运动，随着外激励幅值的增大，外板的振动形式趋于平稳，呈现周期运动形式；当外激励幅值增大至25～38之间时，再次进入混沌运动，继续增大外激励幅值会发现外板再次变为规律的周期运动。图17～20分别是外激励幅值为3， 15， 35和45时外板的波形图、三维相图和Poincaré截面。图14f2=35时中间板的4倍周期运动
Fig.14The period-4 motion of the middle plate when f2=35
圖15f2=65时中间板的概周期运动
Fig.15The quasi-period motion of the middle plate when f2=65图16f2=80时中间板的混沌运动
Fig.16The chaotic motion of the middle plate when f2=80
图17f3=3时外板的混沌运动
Fig.17The chaotic motion of the outer plate when f3=3图18f3=15时外板的周期运动
Fig.18The periodic motion of the outer plate when f3=15
图19f3=35时外板的混沌运动
Fig.19The chaotic motion of the outer plate when f3=35图20f3=45时外板的周期运动
Fig.20The periodic motion of the outer plate when f3=45此外，中间板第5阶模态的振动幅值比第4阶模态的振动幅值大，外板第4，5阶模态的振动幅值基本相差不大，这是由于系统在与第5阶固有频率对应的外激励作用下，在这两阶模态存在1∶2的关系时发生了耦合的内共振现象。
6结论
本文利用Hamilton原理建立了在外激励作用下Z型折叠板的几何非线性动力学方程，并对系统在主参数共振-1∶2内共振情况下的非线性动力学行为进行了摄动分析，得到系统4自由度的平均方程。利用数值方法分析了系统的幅频响应特性和混沌分叉特性。
数值结果表明，不同的外激励幅值和阻尼系数会对系统的频响特性产生一定的影响，且随着调频参数σ的增大，对应的系统振动幅值会出现多值和跳跃的现象。
选取一定的参数和初始条件，通过数值模拟发现，在主参数共振-1∶2内共振的共振关系下，当外激励的频率与系统第5阶固有频率相同时，只改变外激励幅值时，中间板会出现单倍周期-混沌-概周期-混沌运动，外板会出现混沌-单倍周期-概周期-混沌运动，由此可见外激励的改变会对系统的非线性动力学特性产生显著的影响，且系统的第4阶模态对应的幅值也会产生明显的变化，说明此非线性系统的不同模态振动之间存在复杂的耦合关系。
因此，在研究Z型折叠板这一类结构的非线性动力学行为时，不应该只考虑单一的模态振动，还应考虑多阶模态之间的相互作用，以便更好地利用或控制其运动形式，为实际工程提供重要的理论依据。
参考文献：
[1]韩运龙. 折叠板壳结构的设计与分析[D]. 南京：东南大学， 2011：02.
Han Yunlong. Design and analysis of foldable plates structures[D]. Nanjing： Southeast University， 2011：02.
[2]陈务军， 关富玲， 陈向阳. 可折叠航天结构展开动力学分析[J]. 计算力学学报， 1999，16：4.
Chen Wujun， Guan Fulin， Chen Xingyang. Dynamic analysis for deployment process of foldable aerospace structures[J]. Chinese Journal of Computational Mechanics， 1999， 16：4.
[3]Lee S Y， Wooh S C. Finite element vibration analysis of composite box structures using the high order plate theory[J]. Journal of Sound and Vibration， 2004，277（4-5）： 801—814.
[4]Lee S Y， Wooh S C， Yhim S S. Dynamic behavior of folded composite plates analyzed by the third order plate theory[J]. International Journal of Solids and Structures， 2004，41（7）：1879 —1892.
[5]Pal S， Gu H， Niyogi A. Application of folded plate formulation in analyzing stiffened laminated composite and sandwich folded plate vibration[J]. Journal of Reinforced Plastics and Composites， 2008，27（7）：693 —710.
[6]Haldar S， Sheikh A H. Free vibration analysis of isotropic and composite folded plates using a shear flexible element[J]. Finite Elements in Analysis and Design， 2005，42（3）：208—226.
[7]Hernández E， Hervella-Nieto L. Finite element approximation of free vibration of folded plates[J]. Computer Methods in Applied Mechanics and Engineering， 2009，198（15-16）：1360—1367.
[8]Peng L X， Kitipornchai S， Liew K M. Free vibration analysis of folded plate structures by the FSDT mesh-free method[J]. Computational Mechanics， 2006，39（6）：799—814.
[9]Topal U， Uzman . Frequency optimization of laminated folded composite plates[J]. Materials & Design， 2009，30（3）：494 —501.
[10]Jiang R J， Au F T K. A general finite strip for the static and dynamic analyses of folded plates[J]. Thin-Walled Structures， 2011，49（10）：1288—1294.
[11]Liu C， Tian Q， Hu H. Dynamics of a large scale rigid-flexible multibody system composed of composite laminated plates[J]. Multibody System Dynamics， 2011，26（3）：283—305.
[12]Zhang W， Hu W H， Cao D X， et al. Vibration frequencies and modes of a Z-shaped beam with variable folding angles[J]. Journal of Vibration and Acoustics， 2016，138（4）：041004.
[13]Guo X Y， Zhang W， Yao M H. Nonlinear dynamics of angle-ply composite laminated thin plate with third-order shear deformation[J]. Science China Technological Sciences， 2010，53（3）：612—622.
[14]Zhang W， Guo X Y， Lai S K. Research on periodic and chaotic oscillations of composite laminated plates with one-to-one internal resonance[J]. International Journal of Nonlinear Sciences and Numerical Simulation， 2009，10（11-12）：1567—1583.
[15]趙孟良， 关富玲， 吴开成. 空间可展板壳结构的展开分析[J]. 浙江大学学报（工学版）， 2006，40（11）：1837 —1841.
Zhao Mengliang， Guang Fulin， Wu Kaicheng. Deployment analysis of deployable space panel and shell structure[J]. Journal of Zhejiang University（Engineering Science）， 2006，40（11）：1837—1841.
[16]关富玲， 张惠峰， 韩克良. 二维可展板壳结构展开过程分析[J]. 工程设计学报， 2008，15（5）：351—356.
Guan Fuling， Zhang Huifeng， Han Keliang. Deployment analysis of two-dimensional deployable panel and shell structures[J]. Journal of Engineering Design， 2008，15（5）：351—356.
[17]Reddy J N. Mechanics of Laminated Composite Plates and Shells： Theory and Analysis[M]. CRC Press， 2004.
[18]朴金丽.Z 型折叠机翼的非线性动力学研究[D].北京：北京工业大学，2016：05.
Piao Jinli. Nonlinear vibrations for the Z-type folding wings of the morphing aircraft[D]. Beijing： Beijing University of Technology， 2016：05.
Nonlinear vibration characteristics of Z-type folding plates
with internal resonance
GUO Xiang-ying， ZHANG Yang， ZHANG Wei
（Beijing Key Laboratory of Nonlinear Vibrations and Strength of Mechanical Structures，
Beijing University of Technology， Beijing 100124， China）
Abstract： The Z-type folding plate is a kind of complex multi-body structure， which is widely applied to many engineering fields. Based on the classical laminated plate theory and the von Karman type equation， the nonlinear dynamic equations of the Z-type folding plates are obtained by using the Hamilton′s principle. The mode functions of the Z-type folding plates are analyzed with the ANSYS. Then， the Galerk in procedure is used to obtain the normal differential governing equations of the nonlinear system. The case of primary parametric resonance 1∶2 inner resonance is considered. Based on the averaged equation obtained with the method of multiple scales， the numerical simulation is performed to indicate the nonlinear dynamical characteristics of the system.
Keywords： nonlinear dynamics； Z-type folding plates； inner resonance； numerical simulation
摘要：Z型折疊板是一类复杂多体结构，在工程实际中应用广泛。以大型的空间可展开结构为实际工程背景，考虑了Z型折叠板折叠过程中各部件的运动耦合关系，基于Reddy经典板理论和von-Karman几何非线性应变位移关系，利用Hamilton原理建立了Z型折叠板的非线性动力学控制方程。利用有限元分析得到了Z型折叠板结构的模态函数，运用二阶Galerkin截断，得到了Z型折叠板二自由度的非线性常微分方程。考虑系统主参数共振-1∶2内共振的情况，通过摄动分析得到系统四维直角坐标形式的平均方程，最后利用数值模拟方法研究了横向激励对Z型折叠板非线性动力学特性的影响。
关键词： 非线性动力学； Z型折叠板； 内共振； 数值模拟
中图分类号： O322文献标志码： A文章编号： 1004-4523（2018）02-0183-15
DOI：10.16385/j.cnki.issn.1004-4523.2018.02.001
引言
可展开（折叠）结构具有悠久的研究历史和工程应用。折叠板结构是可展结构体系中最常用的一种，近年来在航空航天和建筑领域中得到了广泛的应用。在折叠状态下，结构体积较小，可用于运输或存储；在外力作用或系统内部驱动下，结构逐步展开，最终达到完全展开的工作状态，然后锁定为稳定状态[1-2]。其中，Z型折叠板结构是此类折叠结构中最为经典的一种，在航空领域主要被应用于可变体飞行器机翼结构，而其中运动过程中由于变形引起的结构非线性动态特性问题是结构平稳运行的关键。
在Z型折叠板结构振动特性的研究方面，国内外学者进行了大量的研究，如：Lee等[3-4]分别利用有限元方法和高阶板理论研究了复合材料折叠板的振动特性。Pal等[5-8]分别使用有限元方法，高阶层合板理论研究了复合材料层合板折叠结构的自由振动特性。Topal等[9] 利用一阶板理论建立了具有对称折叠角度的复合材料折叠板的动力学模型，通过MFD方法进行频率优化，数值分析得出层合板的长厚度比、夹角、板的长度以及边界条件会影响到其频率特性。Jian等[10]建立了有限元模型，用一种条状网格划分单元，模拟了折叠板的静态特征和动力学特性分析。Liu等[11]对折叠板多体系统进行了动力学建模，并用数值方法研究了其动态特性。张伟等[12]建立了变角度Z型梁的动力学方程，计算了结构的固有频率和振动模态，并通过数值仿真和实验验证。此外，Zhang等[13-14]研究了复合材料正交铺设层合板结构在外激励作用下的非线性振动问题，分析了系统不同参数下的振动响应。浙江大学的赵孟良和关富玲等[15-16]研究了空间可展结构在外载荷作用下的运动特性。
然而，Z型折叠结构属于多体结构，在外激励作用下将产生大幅的非线性耦合振动响应，这部分的研究成果目前还很有限。本文根据实际工程背景，研究了一类Z型折叠板结构在一定频率的外激励作用下发生1∶2内共振情况下的非线性动力学特性。
1动力学模型和方程〖2〗1.1动力学模型以大型的空间可展开结构为实际工程背景，建立如图1所示的力学模型，所研究的Z型折叠板采用碳纤维复合材料层合板，折叠结构由三部分组成，靠近固定端的部分称为内板，用Ω1表示，完成翻转动作的称为中间板，用Ω2表示，中间板外侧部分为外板，用Ω3表示。图中O1x1y1，O2x2y2，O3x3y3分别为内板、中间板、外板的局部坐标系，全局坐标系Oxy与局部坐标系O1x1y1重合，各个板表面受到的横向简谐激励表示为P1，P2和P3，如图1所示。
模型中，各个部分之间通过刚性铰链相连接，通过内部的作动器和机械结构进行折叠与展开动作，机构在变形时，内板与固定端相连接，中间板在一定的角度范围内转动，外板始终与内板保持平行。
本文研究过程中，满足以下条件：
1）Z型折叠板3块板是等厚度，等宽度的；
2）Z型折叠板3块板在外激励作用下不会产生材料本身的破坏。
图1Z型折叠板的动力学模型
Fig.1Mechanical model of the Z-type folding plates第2期郭翔鹰，等： Z型折叠板内共振下非线性振动特性研究振 动 工 程 学 报第31卷1.2动力学方程的建立
本文选取折叠板材料为正交铺设的碳纤维复合材料层合板，不考虑剪切效应，因此（σz ）z=+ /-h2=0，（τxz）z=+/-h2=0，（τyz）z=+/-h2=0，其中h为层合板结构厚度。
材料的本构关系如下
σxx
σyy
τyz
τzx
τxy=1112000
2122000
004400
000550
000066εxx
εyy
γyz
γzx
γxy（1）
式中ij为转换弹性常数，定义为ij=T-1QijT-1T（2）并且T-1=cos2θsin2θ0
sin2θ〖〗cos2θ0
00cos2θ-sin2θ（3）正交铺设复合材料层合板的弹性模量Qij（i=1，2，4，5，6； j=1，2，4，5，6）可以表示为如下形式：Q11=Q22=E1-ν2，Q12=Q21=νE1-ν2，
Q44=Q55=Q66=E2（1+ν）（4）内力和弯矩的关系可表示为：Nxxi
Nyyi
Nzzi=∑Nk=1∫zk+1zkσxxi
σyyi
σzzidz （5）
Mxxi
Myyi
Mzzi=∑Nk=1∫zk+1zkσxxi
σyyi
σzzizdz（6）本文所研究的Z型可折叠结构在横向外载荷的作用下会产生较大的变形，因此这里用von-Karman大变形理论进行分析，各个板的非线性应变表达式为εx
εy
γxy=ε0+zε1（7）式中ε0=u0ixi+12w0ixi2
v0iyi+12w0iyi2
v0ixi+u0iyi+w0ixiw0iyi（8）
ε1=-ω20ix2i
-ω20iy2i
-2ω20ixiyi （9）式中下标i表示第i个板（i=1，2，3）。
根据正交铺设碳纤维增强复合材料的结构特点，可得到内力与应变关系为N
M=AB
BDε1
ε2 （10）式中剛度矩阵Aij， Bij， Dij， 用如下形式表示Aij，Bij，Dij=∫h2-h2ij1，z，z2dz（11）根据经典的层合板理论[17]，3块板上任意一点在局部坐标系中的位移分别表示为以下形式：ui=u0i（xi，yi，zi，t）-ziw0ixi（12a）
vi=v0i（xi，yi，zi，t）-ziw0iyi（12b）
wi=w0i（xi，yi，zi，t），（i=1，2，3）（12c）式中u0i，v0i，w0i是在板Ωi中性层上任意一点在X， Y， Z方向的位移。
根据本文模型中所建立的坐标转换关系，Z型折叠板结构在全局坐标系下的位移场可描述为：r1=R1u1
v1
w1 （13a）
r2 = r1 L1 cosw1′（L1 ）
v2
w2 + R2 u2
v2
w2 （13b）
r3 = r2 L2 cosw2′（L2 ）
v2
w2 （L2 ） + u3
v3
w3 （13c）式中Li为各个板的长度， Ri 为坐标转换矩阵，形式如下R1=1
1
1 （14a）
R2=cosθ0-sinθ
0〖〗10
sinθ0cosθ （14b）将上述表达式（1）～（14）代入Hamilton原理中，整理展开后得到用广义位移表示的Z型折叠板结构的非线性动力学方程如下。
内板的非线性动力学方程为：
A11（2u01x21+w01x12w01x21）+A12（2u01x21+
w01x12w01x21）+A66（v01x1+u01y1+w01x1w01y1）=
I02u01t2-I13w01t2x1（15a）
A12（2v01y21+w01y12w01y21）+A22（2v01y21+
w01y12w01y21）+A66（v01x1+u01y1+w01x1w01y1）=
I02v01t2-I1（3w01〖〗t2y1）（15b）
A11[2u01x21w01x1+u01x12w01x21+3〖〗2w01x122w01x21]+
A12[2u01x21w01x1+u01x12w01x21+32w01x122w01x21]+
A12[2v01y21w01y1+v01y12w01y21+32w01y122w01y21]+
A22[2v01y21w01y1+v01y12w01y21+32w01y122w01y21]-
（D11+D12）4w01x41-（D22+D12）4w01y41-
8D664ω01x21y21+A66[2v01x12ω01x1y1+2u01y12ω01x1y1+
4w01x1w01y12ω01x1y1+v201x1y1w01x1+2u01x1y1w01x1+ 2v01x21w01y1+2u01x1y1w01y1+2w01x21w01y12+
2w01y21w01x12]+P1=I02w01t2+I1（3u01t2x1+
3v01t2y1）-I2（3w01t2x1+3w01t2y1）（15c）
中间板的非线性动力学方程为：
A11（2u02x22+w02x22w02x22）+A12（2u02x22+
w02x22w02x22）+A66（v02x2+u02y2+w02x2w02y2）=
I0（sinkθ+coskθ）2u02t2+I0（sinkθ+
coskθ）（2w01（L1）tx1+）u02-I1（sinkθ+
coskθ）3w01t2x2-I1（sinkθ+coskθ）A12（2v02〖〗y22+
w02y22w02y22）+A22（2v02y22+w02y22w02y22）+A66（v02x2+
u02y2+w02x2w02y2）=I02v02t2-I1（3w02t2y2）（15d）
A12（2v02y22+w02y22w02y22）+A22（2v02y22+
w02y22w02y22）+A66（v02〖〗x2+u02y2+w02x2w02y2）=
I02v02t2-I1（3w02t2y2） （15e）
A11[2u02x22w02x2+u02x22w02x22+w02x222w02x22]+
A12[2u02x22w02x2+u02〖〗x22w02x22+
32w02x222w02x22]+A12[2v02y22w02y2+
v02y22w02y22+32w02y222w02y22]+A22[2v02y22w02y2+
v02y22w02y22+32w02y222w02y22]-（D11+
D12）4w02x42-（D22+D12）4w02y42-8D664ω02x22y22+
A66[2v02x22ω02x2y2+2u02y22ω02x2y2+
4w02x2w02y22ω02x2y2+v202x2y2w02〖〗x2+
2u02x2y2w02x2+2v02x22w02y2+2u02x2y2w02y2+
2w02x22w02y22+2w02y22w02x22]+P2=
-I02sinkθ2w02t2+I02coskθ（2w01（L1）tx1+
）w02+I1（sinkθ+coskθ）3u02t2x2+I1（sinkθ+
coskθ）（2w01（L1）tx1+）u02x2+I1+3w02t2y2-
I2（sinkθ+coskθ）4w02t2x22+I2（sinkθ+
coskθ）（2w01（L1）tx1+）w202x22-I24w02t2y22 （15f）
外板的非线性动力学方程为：
A11（2u03x23+w03x32w03x23）+A12（2u03x23+
w03〖〗x32w03x23）+A66（v03x3+u03y3+w03〖〗x3w03y3）=
I02u03t2-I13w03t2x3+（-θ）（15g）
A12（2v03y23+w03y32w03〖〗y23）+A22（2v03y23+
w03y32w03y23）+A66（v03x3+u03y3+w03x3w03y3）=
I02v03t2-I1（3w03t2y3）（15h）
A11[2u03x23w03x3+u03x32w03x23+32w03x322w03x23]+
A12[2u03x23w03x3+u03x32w03x23+32w03x322w03x23]+
A12[2v03y23w03y3+v03y32w03y23+32w03y322w03y23]+
A22[2v03y23w03y3+v03y32w03y23+32w03y322w03y23]-
（D11+D12）4w03x43-（D22+D12）4w03y43-D664ω03x23y23+
A66[2v03x32ω03x3y3+2u03y32ω03x3y3+
4w03x3w03y32ω03x3y3+v203x3y3w03x3+2u03x3y3w03x3+
2v03x23w03y3+2u03x3y3w03y3+2w03x23w03y32+
2w03y23w03x32]+P3=I02w03t2+I13u03t2x3+
I13u03t2y3-I24w03t2x23-I24w03t2y23（15i）
結构整体的边界条件满足如下等式：u1（0）=0，u1 x1 = L1 = u2 （0） = L1 （16a）
u2 x2 = L2 = u3 （0） = L1 + L2 cosθ（16b）
u3 x3 = L3 = L1 + L2 cosθ + L3 （16c）
w1（0）=0，w1（L1）=w2（0）=0（16d）
w3 （0） = w2 （L2 ） = L2 sinθ（16e）
w1′（0）=0，w2″（0）=0，w3（0）=0（16f）
w01 （0）y1 x1 = 0 = 0 （16g）
2w2（0）y22x1 = L1 = 2w3（0）y23x2 = L2 = 0（16h）
w01 （b）y1 x1 = 0 = 0（16i）
2w2（b）y22x1 = L1 = 2w3（b）y23x2 = L2 = 0（16j）
2w03 tx3 x3 = L3 = L1 + L2 cosθ + L3 （16k）
（A12+A22）[v0iyi+12（w0iyi）2]y=0，b=0 （16l）
2w02 tx2 x2 = L2 = L1 + L2 cosθ （16m）2有限元模态分析
上述所建立的动力学控制方程是偏微分方程，利用数学方法直接求解极为困难，因此本文将采用Galerkin离散将偏微分方程转换到常微分方程后进行求解。对于单一结构模型，其模态函数通常是可以直接根据经验公式假设的，但Z型折叠板结构为多体结构，很难确定折叠板在外激振力作用下结构的振动模态。因此首先通过ANSYS有限元方法，对Z型折叠板进行模态分析和谐响应分析，得到Z型折叠板结构的固有频率和模态。通过研究结构模态振型，确定系统的模态函数形式。最后，将无量纲形式的Z型折叠板结构的非线性动力学控制方程通过Galerkin方法进行二阶离散，得到可求解的常微分方程组。
2.1有限元模型
航空领域实际应用中，Z型折叠板结构在折叠展开运动过程中的角速度很小，且折叠角度在0～150°的范围内，因此本文将板折叠过程分解，看作是不同折叠角度板结构的准静态慢变过程，选取几个特定折叠角度来建立Z型碳纤维复合材料层合板结构的力学模型。
本章研究中取折叠角度为60°，90°和120°作为典型参数值，进行对比分析和讨论。
有限元模型统一采用四边形板单元，如图2所示；设置弹性模量为5×105 MPa，泊松比为0.3。Z型折叠板的有限元模型尺寸参数如表1所示。
图2模型网格划分示意图
Fig.2Meshing of the Z-type folding plate wing of the three angles
表1Z型折叠板有限元模型的尺寸参数
Tab.1Geometric parameters of the finite element model for the Z-type folding wing
板长度/m厚度/m宽度/m内板20.012.1中间板10.012.1外板40.012.12.2模态分析
通过对建立的有限元模型施加初始条件，进行模态分析，得到Z型折叠板结构不同折叠角度下前5阶固有频率如表2所示。
表2不同折叠角度前5阶固有频率（单位：Hz）
Tab.2The first five order natural frequencies of different folding angles （Unit：Hz）
阶数60°90°120°10.560180.617540.6944322.711402.358002.6580033.225903.319903.6326043.550903.731604.2836057.367306.220905.96140
下面列出Z型折叠板结构在60°，90°，120°的折叠角度下的模态振型图，如图3～5所示。
由有限元模态分析结果可以看出，不同折叠角度下的Z型折叠板结构前5阶模态可表示为： 第1图3折叠角度为60°时Z型折叠板前5阶模态振型图
Fig.3The first five order mode shapes of the Z-type folding angle 60°图4折叠角度为90°时Z型折叠板前5阶模态振型图
Fig.4The first five order mode shapes of the Z-type folding angle 90°圖5折叠角度为120°时Z型折叠板前5阶模态振型图
Fig.5The first five order mode shapes of the Z-type folding angle 120°
阶为弯曲振动，第2阶为扭转振动，第3阶为弯曲振动，第4阶为弯扭耦合振动，第5阶为弯曲振动。通过以上分析可以发现，Z型折叠板结构的前5阶振动模态的形式与悬臂板结构前5阶振动模态的形式相似。
3Galerkin离散
由于结构在共振情况下会发生剧烈的大幅振动，容易产生失稳及结构整体破坏等。因此，对于结构在共振情况下动力学响应的研究是十分必要的。
在有限元分析结果中可以发现，Z型折叠板结构在折叠角度为60°时，第4阶固有频率几乎是第5阶固有频率1/2，存在1∶2内共振的情况，因此，本文将对Z型折叠板结构在主参数共振-1∶2内共振的情况下的动力学特性进行深入分析。
根据上述的分析结果，可选取悬臂板结构的振动模态函数形式作为Z型折叠板在折叠角度为60°时的模态函数，利用Galerkin方法对方程（15）进行二阶截断，得到系统常微分形式的非线性动力学方程。在满足位移边界条件的前提下，选取3个方向的基本函数分别如下：
Ui（x，y，t） = u1 （t）sinx/2πLi cosπyb +
u2 （t）sin3x/2πLi cos2πyb（17a）
Vi （x，y，t） = v1 （t）cosx/2πLi sinπyb +
v2 （t）cos3x/2πLi sin2πyb（17b）
Wi（x，y，t）=w1（t）Xi1（x）Yi1（y）+
w2（t）Xi2（x）Yi2（y） （17c）
根据4阶常微分方程的通解，可将方程（17c）中的Xij（x）和Yi1（y）取为Xij （x） = Ai1 coshki x-Ai2 coski x-
Ai3 βi1 sinhki x + Ai4 sinki x，
Yij （x） = Bi1 coshki x-Bi2 coski x-
Bi3 βi1 sinhki x + Bi4 sinki x（18）式中Xij（x）是沿x轴方向的固支-自由梁函数，Yij（y）为y轴方向的自由-自由梁函数，ki1和ki2为特征方程的根，并有如下关系coski1coshki1+1=0
coski2coshki2-1=0（19）将系统实际参数L1=2 m，L2=1 m，L3=4 m代入式（17）可求得X，Y方向的模态函数，将方程（18），（19）代入边界条件求得模态参数Aij和Bij，得到Xij（x）和Yij（y）的函数表达式，如下：
X11 =-4.9×10 - 5coshk11 x +
2.1×10 - 5cosk11 x - 6.6×10 - 5β11 ·
sinhk11 x-sink11 x （20a）
X12 = 0.25coshk12 x + 27.3cosk12 x -
0.25β12 sinhk12 x - sink12 x（20b）
X21 = -1.4×10 - 6coshk21 x +
1.3×10 - 6cosk21 x - 2.5×10 - 7β21 ·
sinhk21 x-sink21 x（20c）
X22 = 1.2×10 - 4coshk22 x-1.4cosk22 x -
0.12β22 sinhk22 x-sink22 x（20d）
X31 = -1.5×10 - 4coshk31 x +
3.4×10 - 4cosk31 x - 3×10 - 4β31 ·
sinhk31 x-sink31 x（20e）
X32 =0.25coshk32 x + 27.3cosk32 x -
0.25β32 sinhk32 x - sink32 x （20f）
同时，将外激励也进行离散，表示为
Pi=Fi1sin3πxlisinπyb+Fi2sinπxlisin3πyb （21）
用二阶Galerkin方法离散方程（15），并将离散后的面内位移u0，v0用横向位移w0表示，整理可得到Z型折叠板横向振动的常微分运动控制方程：
内板方程：
11+μ111+21w11+α112w12+α113w311+α114w312+ α115w12w211+α116w11w212=f11cos（ξt）（22a）
12+μ212+α121w11+22w12+α123w311+α123w312+α125w12w211+α126w11w212=f12cos（ξt）（22b）
中间板方程：
21+μ321+α21121w21+α21222w21+α21321w22+α21422w22+21w21+α216w22+α217w21u21+
α218w21u22+α219w22u21+α2110w22u22+
（α2111/θ）w321+（α2112/θ）w322+α2113（sinθ+
cosθ）w22w221+α2114（sinθ+cosθ）w21w222=
f21cos（ξt） （22c）
22+μ422+α22121w21+α22222w21+α22321w22+α22422w22+α225w21+22w22+α227w21u21+
α228w21u22+α229w22u21+α2210w22u22+w321+（α2212/θ）w322+α2213（sinθ+cosθ）w22w221+
α2214（sinθ+cosθ）w21w222=f22cos（ξt） （22d）
外板方程：
31+μ531+21w31+（α311+F1cosΩt）ω31+
α312w32+α313w331+α314w332+α315w32w231+
α316w31w232=f31cos（ξt）（22e）
32+μ632+22w32+（α322+F2cosΩt）ω32+
α321w31+α323w331+α324w332+α325w32w231+
α326w31w232=f32cos（ξt）（22f）
式中折叠角θ为3块板的运动方程的连接参数，体现了3块板之间的耦合运动关系。
4摄动分析
对于较复杂的非线性常微分方程，很难求出其精确解，需要用近似解析的方法求其渐近解来替代精确解。摄动分析是近似解析的一种方法，包括直接摄动法、多尺度法以及KBM法等。
因此，本章基于系统主参数共振-1∶2内共振的共振关系，使用多尺度方法进行摄动分析。系统共振关系表示如下：21=14Ω2+εσ1， 22=Ω2+εσ2（23）式中1，2是相应线性系统的第1阶、第2阶固有频率。σ1，σ2为系统的调谐参数，为了方便分析，令Ω=1。
经过计算得到系统的直角坐标下平均方程为：
内板平均方程：
11=-12u1x11-σ11x12-3α113x211x12- 3α113x312-2α116x12x213+x214（24a）
12=-12u1x12+σ11x11+3α113x11x212+x311+ 2α116x11x213+x214（24b）
13=-12u2x13-12σ12x14-32α124x213x14+x314-α125x14x211+x212（24c）
14=-12u2x14+12σ12x13+32α124x214x13+x313+α125x13x211+x212-14f12（24d）
中間板平均方程：
21=-12u3x21-σ11x22+14α212x21x24-x22x23- 3α214x22x221+x222-2α217x21x223+x224（25a）
22=-12u3x22+σ11x21-14α212x21x23-x22x24+ 3α214x21x221+x222+2α217x21x223+x224（25b）
23=-12u4x23-32α225x24x223+x224-
α226x24x221+x222-14σ2x4（25c）
24=-12u4x24+32α225x23x223+x224+
α226x23x221+x222+14σ2x3-14f22 （25d）
外板方程：
31=-12u5x31-σ1x32-3α313x231x32+x332-
2α316x32x233+x234+12f31x32（26a）
32=-12u5x32+σ1x31+3α313x31x232+x331+
2α316x31x233+x234+12f31x32（26b）
33=-12u6x33-12σ2x34-32α324x233x34+x334-α325x34x231+x232 （26c）
34=-12u6x34+12σ2x33+32α324x234x33+x333+α325x33x231+x232-12f32（26d）
5数值模拟
根据数值分析结果发现，Z型折叠板内板的振动幅值很小且多为周期性颤振[18]，考虑到本文篇幅的限制，不再详述，主要讨论Z型折叠机板在橫向激励作用下中间板和外板的非线性振动特性。
5.1幅频响应特性分析
通过数值求解系统的四维平均方程，利用matlab软件绘制3块板的幅频响应曲线，选取外激励幅值和系统阻尼系数为控制参数，研究参数对系统幅频特性的影响。
首先，根据实际参数的取值范围，经过无量纲处理后，选取参数为μ1=0.25，μ3=0.22，μ4=0.14， σ3=-0.015，σ4=-0.014，α214=0.31，α217=0.0002，α225=-0.108，α226=3.5，μ5=0.37，μ6=0.56，σ5=1.77，σ6=1.88，α313=2.49，α316=-3.27，α324=7.16，α325=6.81，α313=2.49。
将2块板的初始条件均设为x10=1.44，x20=1.55，x30=1.35，x40=-1.799。令外激励幅值fi （i=2，3）的值分别为50和100，研究结构幅频响应曲线的变化。图中蓝色曲线和红色曲线分别表示结构第4阶模态和第5阶模态的幅频响应曲线，横坐标为调频参数σi （i=2，3），纵坐标为振动幅值ai （i=1，2）， a1， a2 分别表示第4阶和第5阶的振动幅值。
由图6和7可知，随着外激励幅值的增加，幅频
图6不同外激励幅值下的中间板Ω2幅频响应曲线
Fig.6The frequency-response curves of the middle plate Ω2 to the external excitation amplitude f2
图7不同外激励幅值下的外板Ω3幅频响应曲线
Fig.7The frequency-response curves of the outer plate Ω3 to the external excitation amplitude f3
响应曲线的形态发生了不规律变化。随着调频参数的增加，系统振幅除中间板第5阶模态以外都呈现减小的趋势，并且会出现多值现象和跳跃现象。
为研究系统阻尼系数对幅频响应曲线的影响，分别令阻尼系数μi （i=3，5）的值为0.3和0.8。绘制如下幅频响应曲线。
观察图8和9发现，随着阻尼的增加，系统的振动幅值降低，幅频响应曲线形态发生变化，且对第4阶模态的频响特性影响较大。随着调频参数σ的增加，系统会出现多值和跳跃的现象。
图8不同阻尼下的中间板Ω2幅频响应曲线
Fig.8The frequency-response curves of the middle plate Ω2 to the damping coefficient μ3
图9不同阻尼下的外板Ω3幅频响应曲线
Fig.9The frequency-response curves of the outer plate Ω3 to the damping coefficient μ5
5.2非线性振动响应分析
为了研究Z型折叠板系统在主参数共振-1∶2内共振的情况下的非线性振动响应特性，选取外激励幅值fi为控制参数，研究横向激励幅值对系统产生周期运动和混沌运动的影响。
固定上述结构参数，改变外激励的幅值，运用四阶龙格库塔法对系统运动方程进行数值求解，得到中间板和外板的混沌分叉图，如图10和11所示。
图10中间板Ω2随外激励幅值变化的分叉图
Fig.10The bifurcation diagram of the middle plate Ω2 with the transverse excitation f2
图11外板随外激励幅值变化的分叉图
Fig.11The bifurcation diagram of the outer plate Ω3 with the transverse excitation f3
由图中可以看出，在外激励幅值增大的过程中，结构会出现周期运动-混沌运动-周期运动-混沌运动的变化，说明共振情况下，外激励幅值的变化对于系统的运动稳定性具有重要的影响。
为了更好地描述上述分叉图的特性，首先对中间板随外激励幅值变化的振动特性给出了具体的分析：图12，13和14分别给出了中间板在外激励幅值为3，10和35时的波形图、三维相图和Poincaré截面，此时中间板的振动是从单倍周期进入短暂的混沌运动之后又变为倍周期运动。继续增大外激励幅值到，中间板的振动形式趋于复杂，逐渐变为如图15所示的概周期运动，最后进入图16所示的混沌运动，且不会伴随外激励幅值继续增大而改变运动形态。
图12f2=3时中间板的周期运动
Fig.12The periodic motion of the middle plate when f2=3图13f2=10时中间板的混沌运动
Fig.13The chaotic motion of the middle plate when f2=10
外板随外激励幅值变化的振动特性如图17～20所示。当外激励幅值小于10时，外板呈现不规律的混沌运动，随着外激励幅值的增大，外板的振动形式趋于平稳，呈现周期运动形式；当外激励幅值增大至25～38之间时，再次进入混沌运动，继续增大外激励幅值会发现外板再次变为规律的周期运动。图17～20分别是外激励幅值为3， 15， 35和45时外板的波形图、三维相图和Poincaré截面。图14f2=35时中间板的4倍周期运动
Fig.14The period-4 motion of the middle plate when f2=35
圖15f2=65时中间板的概周期运动
Fig.15The quasi-period motion of the middle plate when f2=65图16f2=80时中间板的混沌运动
Fig.16The chaotic motion of the middle plate when f2=80
图17f3=3时外板的混沌运动
Fig.17The chaotic motion of the outer plate when f3=3图18f3=15时外板的周期运动
Fig.18The periodic motion of the outer plate when f3=15
图19f3=35时外板的混沌运动
Fig.19The chaotic motion of the outer plate when f3=35图20f3=45时外板的周期运动
Fig.20The periodic motion of the outer plate when f3=45此外，中间板第5阶模态的振动幅值比第4阶模态的振动幅值大，外板第4，5阶模态的振动幅值基本相差不大，这是由于系统在与第5阶固有频率对应的外激励作用下，在这两阶模态存在1∶2的关系时发生了耦合的内共振现象。
6结论
本文利用Hamilton原理建立了在外激励作用下Z型折叠板的几何非线性动力学方程，并对系统在主参数共振-1∶2内共振情况下的非线性动力学行为进行了摄动分析，得到系统4自由度的平均方程。利用数值方法分析了系统的幅频响应特性和混沌分叉特性。
数值结果表明，不同的外激励幅值和阻尼系数会对系统的频响特性产生一定的影响，且随着调频参数σ的增大，对应的系统振动幅值会出现多值和跳跃的现象。
选取一定的参数和初始条件，通过数值模拟发现，在主参数共振-1∶2内共振的共振关系下，当外激励的频率与系统第5阶固有频率相同时，只改变外激励幅值时，中间板会出现单倍周期-混沌-概周期-混沌运动，外板会出现混沌-单倍周期-概周期-混沌运动，由此可见外激励的改变会对系统的非线性动力学特性产生显著的影响，且系统的第4阶模态对应的幅值也会产生明显的变化，说明此非线性系统的不同模态振动之间存在复杂的耦合关系。
因此，在研究Z型折叠板这一类结构的非线性动力学行为时，不应该只考虑单一的模态振动，还应考虑多阶模态之间的相互作用，以便更好地利用或控制其运动形式，为实际工程提供重要的理论依据。
参考文献：
[1]韩运龙. 折叠板壳结构的设计与分析[D]. 南京：东南大学， 2011：02.
Han Yunlong. Design and analysis of foldable plates structures[D]. Nanjing： Southeast University， 2011：02.
[2]陈务军， 关富玲， 陈向阳. 可折叠航天结构展开动力学分析[J]. 计算力学学报， 1999，16：4.
Chen Wujun， Guan Fulin， Chen Xingyang. Dynamic analysis for deployment process of foldable aerospace structures[J]. Chinese Journal of Computational Mechanics， 1999， 16：4.
[3]Lee S Y， Wooh S C. Finite element vibration analysis of composite box structures using the high order plate theory[J]. Journal of Sound and Vibration， 2004，277（4-5）： 801—814.
[4]Lee S Y， Wooh S C， Yhim S S. Dynamic behavior of folded composite plates analyzed by the third order plate theory[J]. International Journal of Solids and Structures， 2004，41（7）：1879 —1892.
[5]Pal S， Gu H， Niyogi A. Application of folded plate formulation in analyzing stiffened laminated composite and sandwich folded plate vibration[J]. Journal of Reinforced Plastics and Composites， 2008，27（7）：693 —710.
[6]Haldar S， Sheikh A H. Free vibration analysis of isotropic and composite folded plates using a shear flexible element[J]. Finite Elements in Analysis and Design， 2005，42（3）：208—226.
[7]Hernández E， Hervella-Nieto L. Finite element approximation of free vibration of folded plates[J]. Computer Methods in Applied Mechanics and Engineering， 2009，198（15-16）：1360—1367.
[8]Peng L X， Kitipornchai S， Liew K M. Free vibration analysis of folded plate structures by the FSDT mesh-free method[J]. Computational Mechanics， 2006，39（6）：799—814.
[9]Topal U， Uzman . Frequency optimization of laminated folded composite plates[J]. Materials & Design， 2009，30（3）：494 —501.
[10]Jiang R J， Au F T K. A general finite strip for the static and dynamic analyses of folded plates[J]. Thin-Walled Structures， 2011，49（10）：1288—1294.
[11]Liu C， Tian Q， Hu H. Dynamics of a large scale rigid-flexible multibody system composed of composite laminated plates[J]. Multibody System Dynamics， 2011，26（3）：283—305.
[12]Zhang W， Hu W H， Cao D X， et al. Vibration frequencies and modes of a Z-shaped beam with variable folding angles[J]. Journal of Vibration and Acoustics， 2016，138（4）：041004.
[13]Guo X Y， Zhang W， Yao M H. Nonlinear dynamics of angle-ply composite laminated thin plate with third-order shear deformation[J]. Science China Technological Sciences， 2010，53（3）：612—622.
[14]Zhang W， Guo X Y， Lai S K. Research on periodic and chaotic oscillations of composite laminated plates with one-to-one internal resonance[J]. International Journal of Nonlinear Sciences and Numerical Simulation， 2009，10（11-12）：1567—1583.
[15]趙孟良， 关富玲， 吴开成. 空间可展板壳结构的展开分析[J]. 浙江大学学报（工学版）， 2006，40（11）：1837 —1841.
Zhao Mengliang， Guang Fulin， Wu Kaicheng. Deployment analysis of deployable space panel and shell structure[J]. Journal of Zhejiang University（Engineering Science）， 2006，40（11）：1837—1841.
[16]关富玲， 张惠峰， 韩克良. 二维可展板壳结构展开过程分析[J]. 工程设计学报， 2008，15（5）：351—356.
Guan Fuling， Zhang Huifeng， Han Keliang. Deployment analysis of two-dimensional deployable panel and shell structures[J]. Journal of Engineering Design， 2008，15（5）：351—356.
[17]Reddy J N. Mechanics of Laminated Composite Plates and Shells： Theory and Analysis[M]. CRC Press， 2004.
[18]朴金丽.Z 型折叠机翼的非线性动力学研究[D].北京：北京工业大学，2016：05.
Piao Jinli. Nonlinear vibrations for the Z-type folding wings of the morphing aircraft[D]. Beijing： Beijing University of Technology， 2016：05.
Nonlinear vibration characteristics of Z-type folding plates
with internal resonance
GUO Xiang-ying， ZHANG Yang， ZHANG Wei
（Beijing Key Laboratory of Nonlinear Vibrations and Strength of Mechanical Structures，
Beijing University of Technology， Beijing 100124， China）
Abstract： The Z-type folding plate is a kind of complex multi-body structure， which is widely applied to many engineering fields. Based on the classical laminated plate theory and the von Karman type equation， the nonlinear dynamic equations of the Z-type folding plates are obtained by using the Hamilton′s principle. The mode functions of the Z-type folding plates are analyzed with the ANSYS. Then， the Galerk in procedure is used to obtain the normal differential governing equations of the nonlinear system. The case of primary parametric resonance 1∶2 inner resonance is considered. Based on the averaged equation obtained with the method of multiple scales， the numerical simulation is performed to indicate the nonlinear dynamical characteristics of the system.
Keywords： nonlinear dynamics； Z-type folding plates； inner resonance； numerical simulation
