大展弦比板弯扭耦合受迫振动分析*

王卓 毛晓晔 丁虎 陈立群

(上海大学 力学与工程科学学院 上海力学信息学前沿科学中心,上海 200444)

引言

大展弦比板是工程中的重要结构之一,无论是大展弦比机翼、大跨度桥梁、发动机的叶片等都可以简化为大展弦比板结构进行分析研究,因此对大展弦比板结构进行分析具有重要意义.

随着大展弦比机翼的广泛应用,结构非线性带来的振动问题越来越多[1].Dunn等研究了矩形悬臂机翼在不同弯扭刚度耦合下的非线性、失速、气动弹性行为并利用傅里叶分析进行非线性颤振计算[2].Xie等发现结构大变形引起的几何非线会导致机翼弦向弯曲和扭转耦合运动,从而改变频率和振型[3].Minguet等对结构耦合复合材料叶片的动力学特性进行了分析和实验研究,结果表明静态挠度对叶片扭转和纵向模态和频率有显著影响[4].已经有许多学者采用有限元软件模拟和实验分析的方法对大展弦比结构的振动特性进行了分析.

2012年,Shams等基于非线性Euler-Bernoulli梁行对大展弦比柔性机翼的气动弹性响应进行了预测[5].2014年,Mian对大展弦比的线性和非线性弯曲变形及扭转角进行了分析[6].2019年,Trahair等基于简支梁均匀非弹性屈曲研究了悬臂钢梁的非弹性侧向屈曲[7].2019年,Pezeshky等研究了腹板变形对悬臂梁横向扭转去屈曲强度的影响,量化了不同横向支撑方案对弹性横向扭转屈曲强度的影响[8].2020年,刘燕等对可伸缩复合材料梁的时变非线性振动进行理论研究,分析了可伸缩悬臂复合材料层合梁在外伸与收缩变形过程中的非线性动力学特性[9].2014年,杨智春等发现几何非线性会使扭转固有频率明显下降[10].2021年,田少杰等采用模态叠加法分析了气流激励下叶片的振动响应[11].2021年,张立等研究了主梁偏置角度对弯扭耦合叶片力学性能的影响[12].2022年,朱成秀等通过虚功原理建立轴向运动板的有限元方程,给出了板厚与系统复频率的关系,揭示了板厚对轴向运动三阶剪切变形板稳定性的影响[13].2012年,Zhang等研究了层合板在四边简支及不同荷载激励作用下层合板的非线性动力学振动[14].2022年,韩勤锴等开展了变速旋转圆柱薄壳动力稳定性研究,分别探讨仅考虑周期轴向力、仅考虑变转速以及同时考虑两种时变因素时,系统主参数不稳定区和组合不稳定的变化规律[15].2013年,Zhang等考虑了板的面内激励、横向振动及弯矩激励共同作用下板的非线性动力学响应[16].2011年,Yang研究了四边简支层合板的展向振动及稳定性[17].2011年,Tang等研究了有无内共振时平板的横向非线性自由振动[18].2020～2022年郭翔鹰等研究了不同横向激励、不同材料参数对FGM板动力学特性的影响[19-22].2021年,Shams等研究了含各向异性复合梁的大展弦比机翼在不可压缩流体中的气动弹性失稳问题,研究表明弯曲-扭转耦合刚度对减小或增大机翼的非线性失稳速度有显著影响[23].这些工作研究了大展弦比结构弯曲或扭转变形时的振动特性和非线性响应,但是对于大展弦比结构耦合变形仍然缺乏相关的研究.

2000年,Hashemi等提出了弯扭耦合梁自由振动分析的动力有限元公式[24].2022年,孔嘉翔等基于变分原理建立了柔性帆板的刚柔耦合模型[25].2021年,张婷等研究了具有实际背景的弯曲和扭转联合作用下梁方程组在非线性边界条件下的吸引子[26].2020年,杨兴等利用Lagrange方法推导了FGM厚板的刚柔耦合动力学方程,研究了FGM厚板的横向变形、速度响应频率和固有频率[27].2011年,Hashemi等将动力有限元技术应用于组合梁的拉伸-扭转自由振动分析,评估了系统的固有频率和模态[28].2016年,Carr运用能量法求出了两端固定和一边固定一边简支梁的近似频率方程,对开截面均匀薄壁梁的纯扭转振动进行了分析[29].2017年,郭兵等给出了能够模拟复杂变形并满足边界条件的水平弯曲和扭转变形函数并给出了任意荷载作用下梁的弯矩表达式[30].2009年,Zhang等发现几何非线性会加剧弦向弯曲和扭转刚度的耦合[31].2014年,任智毅给出了水平弯曲频率和扭转频率发生模态交换的存在条件[32].2014年,刘占科等研究了不同类型的荷载组合作用下简支梁的弹性屈曲-扭转屈曲[33].2015年,Eken等研究了薄壁组合梁的弯扭耦合振动,发现固有频率会随着展弦比和厚度比的增大而增大[34].2022年,Wu等基于Timoshenko梁理论对旋转叶片的轴向弯曲耦合机理进行研究,结果表明,叶片轴向响应对裂纹引起的非线性比弯曲响应更加敏感[35].然而,以上大部分研究都基于有限元仿真或者实验研究结果,分析了几何形状和边界条件对大展弦比结构固有特性的影响,而且通常采用实验的手段得到系统的非线性响应,所以目前仍然缺乏相关的理论分析.

基于以上调研,本文将重点研究大展弦比板弯扭耦合受迫振动.第一节利用Hamilton原理建立两端简支的大展弦比模型.第二节运用Galerkin法对控制方程进行截断,得到一系列常微分方程,介绍了用于数值求解的DQEM和DQM法.第三节分析了谐波平衡法的截断收敛性和谐波收敛性,之后将数值与解析结果进行比较,验证了所求响应的正确性,最后得到了不同外激励下系统的幅频特性曲线,分析了外激励对系统受迫振动响应的影响.最后第四节得到结论.

1 力学模型与控制方程

如图1所示,大展弦比板的长度为L,宽度为b,两端为简支端约束,其展向位移为w(x,t),绕中轴线的转为φ(x,t).ρ为大展弦比板的质量密度,A为横截面面积,h为截面高度,E为弹性模量,Ib是梁关于中心轴的横截面惯性矩,It是大展弦比板关于中心轴的横截面极惯性矩,υ为泊松比,G为剪切模量,外激励线性分布在大展弦比板的一侧边缘,其幅值大小为F.

图1 两端简支大展弦比板模型Fig.1 The high aspect ratio plate model with simple support at both ends

图2 大展弦比板弯扭耦合变形图Fig.2 The flexural and torsional coupling deformation diagram of high aspect ratio plate

取大展弦比板微元dV为研究对象,其任意时刻具有的动能为:

(1)

该微元任意时刻具有的势能将由两部分组成:轴向拉压势能以及绕轴向的剪切势能.

首先计算轴向拉压势能,微元在未受到扰动时,其单位弧长为dx,变形之后,其单位弧长变为:

(2)

在受到弯曲扰动后,微元在轴向上的应变为:

(3)

在受到弦向扭转扰动后,微元的应变为:

εφ=(ysinφ+zcosφ)w,xx

(4)

至此可以得到微元在轴向上的应变为:

(5)

因此该大展弦比板的拉压势能为

(6)

因此系统的总势能为:

(7)

据Hamilton变分原理:

(8)

式中Q为外力功,变分后,得到:

(9)

采用Kelvin黏弹性本构关系模型,Λ为黏性阻尼.

(10)

将式(10)代入式(9)后并对其进行分部积分后,可以得到大展弦比板的控制方程为:

(11)

(12)

相应的边界条件为:

(13)

2 求解方法介绍

2.1 谐波平衡法过程

采用n阶Galerkin截断,为了满足两端简支的边界条件,取正弦函数为势函数,仅保留前n阶,将梁的横向位移和扭转位移用广义坐标的形式表示为:

(14)

(15)

其中qi(t),Qi(t)为大展弦比板的广义位移,模态函数sin(iπx/L)为特征函数.权函数取特征函数本身,将广义坐标形式的位移函数(14)和(15)代入两端简支大展弦比板的偏微分方程后,得到含有响应qi(t),Qi(t)的方程H(x,t)和G(x,t),对方程H(x,t)和G(x,t)采用n阶积分截断,将方程的两边同乘以权函数,并在区间[0,L]上积分后得到n个常微分方程,其方程满足:

(16)

(17)

由此得到了关于qi(t),Qi(t)的n个常微分方程组:

j=1,2…,n

(18)

j=1,2…,n

(19)

其中l1-8分别为对应系数,经过简化之后可得.

对控制方程(11)和(12)的解式(14)和(15)中时间函数qi(t)和Qi(t)作出如下假设:

(20)

其中,w和φ为模态阶数,m为谐波阶数.

将式(20)及其导函数代入经过伽辽金截断预处理后得到的n个常微分方程组(18)及方程组(19),可以得到以含t的各阶谐波为未知变量的n个代数方程组.由于t的任意性,代数方程中各谐波项的系数和常数项均为0,所以需要提取代数方程组中(2m+1)×n项谐波系数,并使它们都等于0,此时可以得到(2m+1)×n个非线性齐次方程组,记为:

F(aw,0,aw,1,bw,1,…aw,m,bw,m,Ω)=0

(21)

F(aφ,0,aφ,1,bφ,1,…aφ,m,bφ,m,Ω)=0

(22)

通过方程组(21)和(22),可以得到式(20)中各系数与激励频率Ω的关系,由于方程组(21)和(22)很难用解析方法求解,而只用牛顿迭代数值方法求解会遇到转折点的奇异性问题.故此处采用伪弧长延伸法中的预报-修正进行数值求解,即用解曲线c(s)的弧长来避免遇到转折点的奇异性问题.具体原理如下:

记y=(aw,0,aw,1,bw,1,…,aw,m,bw,m,Ω),H=(2m+1) ×n.方程组的雅各比矩阵Y可以表示为

(23)

设J={J1,J2,…,JH+1}T为方程组的解曲线c(s)在该H+1维空间中的切向量,J中元素的表达式为:

(24)

结合雅各比矩阵Y和J的表达式,易知

Y·J=0

(25)

式(25)说明J确实是方程组(21)的解曲线c(s)的切向量.定义该H+1维空间中解曲线的弧长微段满足几何关系,即

(26)

解曲线的单位切向量可以表示为:

(27)

则伪弧长延伸法的预测过程可表示为:

dy=τdsy(0)=y0

(28)

对式(28),采用经典的欧拉法进行数值求解,有

yp=y0+τ(y0)Δs

(29)

其中yp为预测解,y0为初值,Δs为给定的弧长增量.修正过程采用牛顿迭代法,其计算格式为:

yc,0=yp

i=1,2,…

(30)

在用伪弧长延伸法求解时可能会遇到以下情况:如果‖JT(yc,0)‖→0,则预测过程的欧拉法失效,修正过程的牛顿迭代法也可能会产生奇异点.如果,‖JT(yc,0)‖→∞,则预测过程的欧拉法中yp=y0,导致修正过程的牛顿迭代法又会回到前一时刻的解,这样会一直在该处循环,程序不能停止,也不能得到全频域上的解,故初值的选取及判定每一步迭代对应的‖J‖是否为很小值或很大值显得尤为重要.基于上述伪弧长法,解出谐波系数,进而得到平板弯扭耦合振动的位移响应.

2.2 直接数值仿真方法过程

由于微分求积法(DQM)可以严格满足两个边界条件,因此它适用于具有两个边值问题的连续体仿真,而微分求积单元法(DQEM)适用于具有四个边值问题的连续体仿真.由于式(11)是四阶微分方程,而式(12)是二阶微分方程,应该采用DQEM和DQM一起进行计算.

按照微分求积法,函数f(x,t)在节点xi处对变量x的偏导数可以表示为:

(31)

式中,算子A可以由以下显示表达式计算.

(i,j=1,2...,N;j≠i)

(n=2,3...;j=1,2...,N;j≠i)

(n=2,3…;i=1,2…,N)

(32)

已有文献表明当节点采用非均匀分布时,计算结果精度及收敛性更高,因此本文采用Chebyshev-Gauss-Lobatto非均匀分布,即

(33)

DQM过程将产生N-2个内点动力学方程,联合两个边界条件,方程数目与节点数目匹配,系统可以求解.

微分方程w(x,t)的解函数可以用DQEM表示为

ψ1(x)w(1)(x1)+ψN(x)w(1)(xN)

(34)

其中N为节点数.φj(x)和w(xj)分别代表插值函数和j点的位移解.w(1)(x1)和w(1)(xN)表示两个新引入的变量,表示边界处的旋转角度.节点xi处的k阶导数可以如下表示

i=1,2…,N

(35)

j=1,N,i=1,2,…,N

(36)

(37)

lj(x)为拉格朗日插值函数,其表达式为

(38)

(39)

(40)

为了提高仿真精度及收敛性,几点分部依然选取Chebyshev-Gauss-Lobatto非均匀分布,即式(33).

3 受迫振动响应算例分析

采用表1的参数进行计算,系统的初始激励幅值为5N,外激励频率为20Hz,数值计算时常取300个周期,每个周期内取50个点,初始位移和初始速度均取为0.

3.1 截断收敛性和谐波收敛性分析

用于数值计算的参数同表1,图3(a)对比了展弦比为10时,二阶、三阶、四阶截断得到的大展弦比中点处弯曲响应的解析计算结果.对比之后发现,二阶截断与三阶和四阶截断得到的结果差异较大,而三阶截断和四阶截断的结果吻合很好.图3(b)对比了展弦比为10时,二阶、三阶截断得到的大展弦比中点处扭转响应的解析计算结果.对比之后发现,二阶截断与三阶得到的结果相同.为了方便后续的计算,故本文后续研究大展均采用三阶截断,即取式(14)和式(15)中的n=3.

图4(a)和(b)对比了展现比为10时三阶谐波次数和五阶谐波次数下大展弦比中点处的弯曲响应和扭转响应计算结果.图中的对比可以发现,三阶谐波系数和五阶谐波系数的计算结果吻合较好,故本文后续研究均采用三阶谐波系数,即取式(20)中的N=3.经过分析,确定了截断收敛性和谐波收敛性,在后文的研究过程中,只需要采用相对应的截断阶数和谐波阶数即可满足计算精度.

(a) 弯曲响应谐波收敛性分析

3.2 数值验证与分析

为了验证3.1节所计算数据的正确性,基于3.1节所确定的截断收敛阶数和谐波收敛阶数,采用数值仿真的方法对解析计算的结果进行验证.同样采用表1的数值进行计算,对于DQEM/DQM,给定所有的初值都设置为零.图5给出了F=5N,ω=20Hz时大展弦比模型弯曲响应的近似解析解和数值解,图6给出了相同条件下大展弦比模型扭转响应的近似解析解和数值解,可以看出此时系统已经进入稳态响应,仿真时长足够,取响应进入稳态最后一个周期的最大位移作为幅值,验证大展弦比模型稳态响应的近似解析解.图中红色点划线表示利用DQM/DQEM计算的得到数值解,蓝色虚线表示运用谐波平衡法计算得到的近似解析结果.经过对比之后可以发现,当系统进入稳态响应后,解析结果与数值结果两者相同,说明近似解析解具有令人满意的计算精度.

图5 大展弦比板中点处稳态弯曲响应时程图Fig.5 Time history of steady-state response

图6 大展弦比板中点处稳态扭转响应时程图Fig.6 Time history of steady-state response

为了分析外激励幅值的影响,图7给出了第一阶模态响应幅值在不同激励频率下随外激励幅值的连续变化曲线,由于控制方程中的三次方项,从图中可以看出随着外激励的增大,幅频特性曲线向右偏转,幅频特性曲线呈现硬特性,并且随着外激励的增大,硬特性增强,同时响应幅值随着外激励的增大而增大,结构的非线性增强.

图7 大展弦比板中点处幅频响应曲线Fig.7 Amplitude-frequency response curve at mid-point

4 结论

本文研究了大展弦比模型的非线性弯扭耦合受迫振动,建立了两端简支的大展弦比模型并利用Hamilton原理推导得到了对应的控制方程,运用谐波平衡法得到了系统受迫振动的稳态响应,无论是对于弯曲振动还是扭转振动来说,三阶截断和三阶谐波系数就能达到很好的收敛性.之后运用DQEM和DQM相结合的方法得到了数值解,对比解析解后发现两者具有较好的精度,验证了解析结果的正确性.最后分析了不同外激励幅值下结构的幅频特性曲线,随着外激励的增大,幅频特性曲线向右弯曲,振动响应幅值增大,结构呈现硬特性.