王卓 毛曉曄 丁虎 陳立群

(上海大學(xué) 力學(xué)與工程科學(xué)學(xué)院 上海力學(xué)信息學(xué)前沿科學(xué)中心,上海 200444)

引言

大展弦比板是工程中的重要結(jié)構(gòu)之一,無論是大展弦比機(jī)翼、大跨度橋梁、發(fā)動機(jī)的葉片等都可以簡化為大展弦比板結(jié)構(gòu)進(jìn)行分析研究,因此對大展弦比板結(jié)構(gòu)進(jìn)行分析具有重要意義.

隨著大展弦比機(jī)翼的廣泛應(yīng)用,結(jié)構(gòu)非線性帶來的振動問題越來越多[1].Dunn等研究了矩形懸臂機(jī)翼在不同彎扭剛度耦合下的非線性、失速、氣動彈性行為并利用傅里葉分析進(jìn)行非線性顫振計算[2].Xie等發(fā)現(xiàn)結(jié)構(gòu)大變形引起的幾何非線會導(dǎo)致機(jī)翼弦向彎曲和扭轉(zhuǎn)耦合運動,從而改變頻率和振型[3].Minguet等對結(jié)構(gòu)耦合復(fù)合材料葉片的動力學(xué)特性進(jìn)行了分析和實驗研究,結(jié)果表明靜態(tài)撓度對葉片扭轉(zhuǎn)和縱向模態(tài)和頻率有顯著影響[4].已經(jīng)有許多學(xué)者采用有限元軟件模擬和實驗分析的方法對大展弦比結(jié)構(gòu)的振動特性進(jìn)行了分析.

2012年,Shams等基于非線性Euler-Bernoulli梁行對大展弦比柔性機(jī)翼的氣動彈性響應(yīng)進(jìn)行了預(yù)測[5].2014年,Mian對大展弦比的線性和非線性彎曲變形及扭轉(zhuǎn)角進(jìn)行了分析[6].2019年,Trahair等基于簡支梁均勻非彈性屈曲研究了懸臂鋼梁的非彈性側(cè)向屈曲[7].2019年,Pezeshky等研究了腹板變形對懸臂梁橫向扭轉(zhuǎn)去屈曲強(qiáng)度的影響,量化了不同橫向支撐方案對彈性橫向扭轉(zhuǎn)屈曲強(qiáng)度的影響[8].2020年,劉燕等對可伸縮復(fù)合材料梁的時變非線性振動進(jìn)行理論研究,分析了可伸縮懸臂復(fù)合材料層合梁在外伸與收縮變形過程中的非線性動力學(xué)特性[9].2014年,楊智春等發(fā)現(xiàn)幾何非線性會使扭轉(zhuǎn)固有頻率明顯下降[10].2021年,田少杰等采用模態(tài)疊加法分析了氣流激勵下葉片的振動響應(yīng)[11].2021年,張立等研究了主梁偏置角度對彎扭耦合葉片力學(xué)性能的影響[12].2022年,朱成秀等通過虛功原理建立軸向運動板的有限元方程,給出了板厚與系統(tǒng)復(fù)頻率的關(guān)系,揭示了板厚對軸向運動三階剪切變形板穩(wěn)定性的影響[13].2012年,Zhang等研究了層合板在四邊簡支及不同荷載激勵作用下層合板的非線性動力學(xué)振動[14].2022年,韓勤鍇等開展了變速旋轉(zhuǎn)圓柱薄殼動力穩(wěn)定性研究,分別探討僅考慮周期軸向力、僅考慮變轉(zhuǎn)速以及同時考慮兩種時變因素時,系統(tǒng)主參數(shù)不穩(wěn)定區(qū)和組合不穩(wěn)定的變化規(guī)律[15].2013年,Zhang等考慮了板的面內(nèi)激勵、橫向振動及彎矩激勵共同作用下板的非線性動力學(xué)響應(yīng)[16].2011年,Yang研究了四邊簡支層合板的展向振動及穩(wěn)定性[17].2011年,Tang等研究了有無內(nèi)共振時平板的橫向非線性自由振動[18].2020～2022年郭翔鷹等研究了不同橫向激勵、不同材料參數(shù)對FGM板動力學(xué)特性的影響[19-22].2021年,Shams等研究了含各向異性復(fù)合梁的大展弦比機(jī)翼在不可壓縮流體中的氣動彈性失穩(wěn)問題,研究表明彎曲-扭轉(zhuǎn)耦合剛度對減小或增大機(jī)翼的非線性失穩(wěn)速度有顯著影響[23].這些工作研究了大展弦比結(jié)構(gòu)彎曲或扭轉(zhuǎn)變形時的振動特性和非線性響應(yīng),但是對于大展弦比結(jié)構(gòu)耦合變形仍然缺乏相關(guān)的研究.

2000年,Hashemi等提出了彎扭耦合梁自由振動分析的動力有限元公式[24].2022年,孔嘉翔等基于變分原理建立了柔性帆板的剛?cè)狁詈夏Ｐ蚚25].2021年,張婷等研究了具有實際背景的彎曲和扭轉(zhuǎn)聯(lián)合作用下梁方程組在非線性邊界條件下的吸引子[26].2020年,楊興等利用Lagrange方法推導(dǎo)了FGM厚板的剛?cè)狁詈蟿恿W(xué)方程,研究了FGM厚板的橫向變形、速度響應(yīng)頻率和固有頻率[27].2011年,Hashemi等將動力有限元技術(shù)應(yīng)用于組合梁的拉伸-扭轉(zhuǎn)自由振動分析,評估了系統(tǒng)的固有頻率和模態(tài)[28].2016年,Carr運用能量法求出了兩端固定和一邊固定一邊簡支梁的近似頻率方程,對開截面均勻薄壁梁的純扭轉(zhuǎn)振動進(jìn)行了分析[29].2017年,郭兵等給出了能夠模擬復(fù)雜變形并滿足邊界條件的水平彎曲和扭轉(zhuǎn)變形函數(shù)并給出了任意荷載作用下梁的彎矩表達(dá)式[30].2009年,Zhang等發(fā)現(xiàn)幾何非線性會加劇弦向彎曲和扭轉(zhuǎn)剛度的耦合[31].2014年,任智毅給出了水平彎曲頻率和扭轉(zhuǎn)頻率發(fā)生模態(tài)交換的存在條件[32].2014年,劉占科等研究了不同類型的荷載組合作用下簡支梁的彈性屈曲-扭轉(zhuǎn)屈曲[33].2015年,Eken等研究了薄壁組合梁的彎扭耦合振動,發(fā)現(xiàn)固有頻率會隨著展弦比和厚度比的增大而增大[34].2022年,Wu等基于Timoshenko梁理論對旋轉(zhuǎn)葉片的軸向彎曲耦合機(jī)理進(jìn)行研究,結(jié)果表明,葉片軸向響應(yīng)對裂紋引起的非線性比彎曲響應(yīng)更加敏感[35].然而,以上大部分研究都基于有限元仿真或者實驗研究結(jié)果,分析了幾何形狀和邊界條件對大展弦比結(jié)構(gòu)固有特性的影響,而且通常采用實驗的手段得到系統(tǒng)的非線性響應(yīng),所以目前仍然缺乏相關(guān)的理論分析.

基于以上調(diào)研,本文將重點研究大展弦比板彎扭耦合受迫振動.第一節(jié)利用Hamilton原理建立兩端簡支的大展弦比模型.第二節(jié)運用Galerkin法對控制方程進(jìn)行截斷,得到一系列常微分方程,介紹了用于數(shù)值求解的DQEM和DQM法.第三節(jié)分析了諧波平衡法的截斷收斂性和諧波收斂性,之后將數(shù)值與解析結(jié)果進(jìn)行比較,驗證了所求響應(yīng)的正確性,最后得到了不同外激勵下系統(tǒng)的幅頻特性曲線,分析了外激勵對系統(tǒng)受迫振動響應(yīng)的影響.最后第四節(jié)得到結(jié)論.

1 力學(xué)模型與控制方程

如圖1所示,大展弦比板的長度為L,寬度為b,兩端為簡支端約束,其展向位移為w(x,t),繞中軸線的轉(zhuǎn)為φ(x,t).ρ為大展弦比板的質(zhì)量密度,A為橫截面面積,h為截面高度,E為彈性模量,Ib是梁關(guān)于中心軸的橫截面慣性矩,It是大展弦比板關(guān)于中心軸的橫截面極慣性矩,υ為泊松比,G為剪切模量,外激勵線性分布在大展弦比板的一側(cè)邊緣,其幅值大小為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,變形之后,其單位弧長變?yōu)?

(2)

在受到彎曲擾動后,微元在軸向上的應(yīng)變?yōu)?

(3)

在受到弦向扭轉(zhuǎn)擾動后,微元的應(yīng)變?yōu)?

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

(4)

至此可以得到微元在軸向上的應(yīng)變?yōu)?

(5)

因此該大展弦比板的拉壓勢能為

(6)

因此系統(tǒng)的總勢能為:

(7)

據(jù)Hamilton變分原理:

(8)

式中Q為外力功,變分后,得到:

(9)

采用Kelvin黏彈性本構(gòu)關(guān)系模型,Λ為黏性阻尼.

(10)

將式(10)代入式(9)后并對其進(jìn)行分部積分后,可以得到大展弦比板的控制方程為:

(11)

(12)

相應(yīng)的邊界條件為:

(13)

2 求解方法介紹

2.1 諧波平衡法過程

采用n階Galerkin截斷,為了滿足兩端簡支的邊界條件,取正弦函數(shù)為勢函數(shù),僅保留前n階,將梁的橫向位移和扭轉(zhuǎn)位移用廣義坐標(biāo)的形式表示為:

(14)

(15)

其中qi(t),Qi(t)為大展弦比板的廣義位移,模態(tài)函數(shù)sin(iπx/L)為特征函數(shù).權(quán)函數(shù)取特征函數(shù)本身,將廣義坐標(biāo)形式的位移函數(shù)(14)和(15)代入兩端簡支大展弦比板的偏微分方程后,得到含有響應(yīng)qi(t),Qi(t)的方程H(x,t)和G(x,t),對方程H(x,t)和G(x,t)采用n階積分截斷,將方程的兩邊同乘以權(quán)函數(shù),并在區(qū)間[0,L]上積分后得到n個常微分方程,其方程滿足:

(16)

(17)

由此得到了關(guān)于qi(t),Qi(t)的n個常微分方程組:

j=1,2…,n

(18)

j=1,2…,n

(19)

其中l(wèi)1-8分別為對應(yīng)系數(shù),經(jīng)過簡化之后可得.

對控制方程(11)和(12)的解式(14)和(15)中時間函數(shù)qi(t)和Qi(t)作出如下假設(shè):

(20)

其中,w和φ為模態(tài)階數(shù),m為諧波階數(shù).

將式(20)及其導(dǎo)函數(shù)代入經(jīng)過伽遼金截斷預(yù)處理后得到的n個常微分方程組(18)及方程組(19),可以得到以含t的各階諧波為未知變量的n個代數(shù)方程組.由于t的任意性,代數(shù)方程中各諧波項的系數(shù)和常數(shù)項均為0,所以需要提取代數(shù)方程組中(2m+1)×n項諧波系數(shù),并使它們都等于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)中各系數(shù)與激勵頻率Ω的關(guān)系,由于方程組(21)和(22)很難用解析方法求解,而只用牛頓迭代數(shù)值方法求解會遇到轉(zhuǎn)折點的奇異性問題.故此處采用偽弧長延伸法中的預(yù)報-修正進(jìn)行數(shù)值求解,即用解曲線c(s)的弧長來避免遇到轉(zhuǎn)折點的奇異性問題.具體原理如下:

記y=(aw,0,aw,1,bw,1,…,aw,m,bw,m,Ω),H=(2m+1) ×n.方程組的雅各比矩陣Y可以表示為

(23)

設(shè)J={J1,J2,…,JH+1}T為方程組的解曲線c(s)在該H+1維空間中的切向量,J中元素的表達(dá)式為:

(24)

結(jié)合雅各比矩陣Y和J的表達(dá)式,易知

Y·J=0

(25)

式(25)說明J確實是方程組(21)的解曲線c(s)的切向量.定義該H+1維空間中解曲線的弧長微段滿足幾何關(guān)系,即

(26)

解曲線的單位切向量可以表示為:

(27)

則偽弧長延伸法的預(yù)測過程可表示為:

dy=τdsy(0)=y0

(28)

對式(28),采用經(jīng)典的歐拉法進(jìn)行數(shù)值求解,有

yp=y0+τ(y0)Δs

(29)

其中yp為預(yù)測解,y0為初值,Δs為給定的弧長增量.修正過程采用牛頓迭代法,其計算格式為:

yc,0=yp

i=1,2,…

(30)

在用偽弧長延伸法求解時可能會遇到以下情況:如果‖JT(yc,0)‖→0,則預(yù)測過程的歐拉法失效,修正過程的牛頓迭代法也可能會產(chǎn)生奇異點.如果,‖JT(yc,0)‖→∞,則預(yù)測過程的歐拉法中yp=y0,導(dǎo)致修正過程的牛頓迭代法又會回到前一時刻的解,這樣會一直在該處循環(huán),程序不能停止,也不能得到全頻域上的解,故初值的選取及判定每一步迭代對應(yīng)的‖J‖是否為很小值或很大值顯得尤為重要.基于上述偽弧長法,解出諧波系數(shù),進(jìn)而得到平板彎扭耦合振動的位移響應(yīng).

2.2 直接數(shù)值仿真方法過程

由于微分求積法(DQM)可以嚴(yán)格滿足兩個邊界條件,因此它適用于具有兩個邊值問題的連續(xù)體仿真,而微分求積單元法(DQEM)適用于具有四個邊值問題的連續(xù)體仿真.由于式(11)是四階微分方程,而式(12)是二階微分方程,應(yīng)該采用DQEM和DQM一起進(jìn)行計算.

按照微分求積法,函數(shù)f(x,t)在節(jié)點xi處對變量x的偏導(dǎo)數(shù)可以表示為:

(31)

式中,算子A可以由以下顯示表達(dá)式計算.

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

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

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

(32)

已有文獻(xiàn)表明當(dāng)節(jié)點采用非均勻分布時,計算結(jié)果精度及收斂性更高,因此本文采用Chebyshev-Gauss-Lobatto非均勻分布,即

(33)

DQM過程將產(chǎn)生N-2個內(nèi)點動力學(xué)方程,聯(lián)合兩個邊界條件,方程數(shù)目與節(jié)點數(shù)目匹配,系統(tǒng)可以求解.

微分方程w(x,t)的解函數(shù)可以用DQEM表示為

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

(34)

其中N為節(jié)點數(shù).φj(x)和w(xj)分別代表插值函數(shù)和j點的位移解.w(1)(x1)和w(1)(xN)表示兩個新引入的變量,表示邊界處的旋轉(zhuǎn)角度.節(jié)點xi處的k階導(dǎo)數(shù)可以如下表示

i=1,2…,N

(35)

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

(36)

(37)

lj(x)為拉格朗日插值函數(shù),其表達(dá)式為

(38)

(39)

(40)

為了提高仿真精度及收斂性,幾點分部依然選取Chebyshev-Gauss-Lobatto非均勻分布,即式(33).

3 受迫振動響應(yīng)算例分析

采用表1的參數(shù)進(jìn)行計算,系統(tǒng)的初始激勵幅值為5N,外激勵頻率為20Hz,數(shù)值計算時常取300個周期,每個周期內(nèi)取50個點,初始位移和初始速度均取為0.

3.1 截斷收斂性和諧波收斂性分析

用于數(shù)值計算的參數(shù)同表1,圖3(a)對比了展弦比為10時,二階、三階、四階截斷得到的大展弦比中點處彎曲響應(yīng)的解析計算結(jié)果.對比之后發(fā)現(xiàn),二階截斷與三階和四階截斷得到的結(jié)果差異較大,而三階截斷和四階截斷的結(jié)果吻合很好.圖3(b)對比了展弦比為10時,二階、三階截斷得到的大展弦比中點處扭轉(zhuǎn)響應(yīng)的解析計算結(jié)果.對比之后發(fā)現(xiàn),二階截斷與三階得到的結(jié)果相同.為了方便后續(xù)的計算,故本文后續(xù)研究大展均采用三階截斷,即取式(14)和式(15)中的n=3.

圖4(a)和(b)對比了展現(xiàn)比為10時三階諧波次數(shù)和五階諧波次數(shù)下大展弦比中點處的彎曲響應(yīng)和扭轉(zhuǎn)響應(yīng)計算結(jié)果.圖中的對比可以發(fā)現(xiàn),三階諧波系數(shù)和五階諧波系數(shù)的計算結(jié)果吻合較好,故本文后續(xù)研究均采用三階諧波系數(shù),即取式(20)中的N=3.經(jīng)過分析,確定了截斷收斂性和諧波收斂性,在后文的研究過程中,只需要采用相對應(yīng)的截斷階數(shù)和諧波階數(shù)即可滿足計算精度.

(a) 彎曲響應(yīng)諧波收斂性分析

3.2 數(shù)值驗證與分析

為了驗證3.1節(jié)所計算數(shù)據(jù)的正確性,基于3.1節(jié)所確定的截斷收斂階數(shù)和諧波收斂階數(shù),采用數(shù)值仿真的方法對解析計算的結(jié)果進(jìn)行驗證.同樣采用表1的數(shù)值進(jìn)行計算,對于DQEM/DQM,給定所有的初值都設(shè)置為零.圖5給出了F=5N,ω=20Hz時大展弦比模型彎曲響應(yīng)的近似解析解和數(shù)值解,圖6給出了相同條件下大展弦比模型扭轉(zhuǎn)響應(yīng)的近似解析解和數(shù)值解,可以看出此時系統(tǒng)已經(jīng)進(jìn)入穩(wěn)態(tài)響應(yīng),仿真時長足夠,取響應(yīng)進(jìn)入穩(wěn)態(tài)最后一個周期的最大位移作為幅值,驗證大展弦比模型穩(wěn)態(tài)響應(yīng)的近似解析解.圖中紅色點劃線表示利用DQM/DQEM計算的得到數(shù)值解,藍(lán)色虛線表示運用諧波平衡法計算得到的近似解析結(jié)果.經(jīng)過對比之后可以發(fā)現(xiàn),當(dāng)系統(tǒng)進(jìn)入穩(wěn)態(tài)響應(yīng)后,解析結(jié)果與數(shù)值結(jié)果兩者相同,說明近似解析解具有令人滿意的計算精度.

圖5 大展弦比板中點處穩(wěn)態(tài)彎曲響應(yīng)時程圖Fig.5 Time history of steady-state response

圖6 大展弦比板中點處穩(wěn)態(tài)扭轉(zhuǎn)響應(yīng)時程圖Fig.6 Time history of steady-state response

為了分析外激勵幅值的影響,圖7給出了第一階模態(tài)響應(yīng)幅值在不同激勵頻率下隨外激勵幅值的連續(xù)變化曲線,由于控制方程中的三次方項,從圖中可以看出隨著外激勵的增大,幅頻特性曲線向右偏轉(zhuǎn),幅頻特性曲線呈現(xiàn)硬特性,并且隨著外激勵的增大,硬特性增強(qiáng),同時響應(yīng)幅值隨著外激勵的增大而增大,結(jié)構(gòu)的非線性增強(qiáng).

圖7 大展弦比板中點處幅頻響應(yīng)曲線Fig.7 Amplitude-frequency response curve at mid-point

4 結(jié)論

本文研究了大展弦比模型的非線性彎扭耦合受迫振動,建立了兩端簡支的大展弦比模型并利用Hamilton原理推導(dǎo)得到了對應(yīng)的控制方程,運用諧波平衡法得到了系統(tǒng)受迫振動的穩(wěn)態(tài)響應(yīng),無論是對于彎曲振動還是扭轉(zhuǎn)振動來說,三階截斷和三階諧波系數(shù)就能達(dá)到很好的收斂性.之后運用DQEM和DQM相結(jié)合的方法得到了數(shù)值解,對比解析解后發(fā)現(xiàn)兩者具有較好的精度,驗證了解析結(jié)果的正確性.最后分析了不同外激勵幅值下結(jié)構(gòu)的幅頻特性曲線,隨著外激勵的增大,幅頻特性曲線向右彎曲,振動響應(yīng)幅值增大,結(jié)構(gòu)呈現(xiàn)硬特性.