大展弦比板彎扭耦合受迫振動分析*

王卓 毛曉曄 丁虎 陳立群

(上海大學 力學與工程科學學院 上海力學信息學前沿科學中心,上海 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相結合的方法得到了數值解,對比解析解后發現兩者具有較好的精度,驗證了解析結果的正確性.最后分析了不同外激勵幅值下結構的幅頻特性曲線,隨著外激勵的增大,幅頻特性曲線向右彎曲,振動響應幅值增大,結構呈現硬特性.