一種幾何非線性復(fù)合材料梁瞬態(tài)動(dòng)響應(yīng)計(jì)算方法*

張勇 霍佳波 宋小輝 張凌云

(桂林航天工業(yè)學(xué)院 機(jī)電工程學(xué)院，廣西 桂林 541004)

先進(jìn)的復(fù)合材料具有高比強(qiáng)度、高比模量、優(yōu)越的抗疲勞和耐環(huán)境等性能，在航空、航天、船舶、車輛等領(lǐng)域取得了廣泛的應(yīng)用[1]。高長(zhǎng)細(xì)比復(fù)合材料梁在航空航天領(lǐng)域、汽車領(lǐng)域有著廣泛的用途，比如直升機(jī)旋翼葉片、飛機(jī)機(jī)翼大梁、汽車減振板簧等。細(xì)長(zhǎng)梁結(jié)構(gòu)在受載狀態(tài)下，往往表現(xiàn)出幾何非線性變形的特征。復(fù)合材料因?yàn)榇嬖诟飨虍愋裕谧冃芜^程中，還存在拉伸、彎曲、扭轉(zhuǎn)多自由度之間的耦合以及與截面翹曲變形的復(fù)雜耦合效應(yīng)。另外，作為結(jié)構(gòu)的主承力梁，往往面臨突加載荷的作用。因此對(duì)復(fù)合材料中等變形梁進(jìn)行瞬態(tài)響應(yīng)研究，可充分發(fā)揮復(fù)合材料結(jié)構(gòu)的可設(shè)計(jì)性，對(duì)復(fù)合材料梁結(jié)構(gòu)設(shè)計(jì)有著重要的意義。

幾何非線性中等變形梁模型主要由Hodges等人完成的。1974年，Hodges和Dowell提出了中等變形梁模型，推導(dǎo)過程中使用階次準(zhǔn)則忽略那些由于中等變形引起的高階項(xiàng)，最終給出適應(yīng)于均勻的、各向同性的揮-擺-扭耦合的非線性槳葉運(yùn)動(dòng)方程[2]。隨后，Tong、Kaza等人進(jìn)一步發(fā)展和擴(kuò)充了中等變形梁理論，考慮了結(jié)構(gòu)揮-擺-扭耦合以及預(yù)扭、預(yù)錐等。1985年，Hodges舍棄了中等變形假設(shè)，建立了另一套考慮了剪切變形、彈性軸初始彎曲和預(yù)扭的彎扭耦合旋翼槳葉運(yùn)動(dòng)方程，即小應(yīng)變有限轉(zhuǎn)角旋轉(zhuǎn)梁理論。隨后為了使建立的模型適用于復(fù)合材料槳葉[3]，1990年，Hodges沒有再限制彈性轉(zhuǎn)角的大小，并且采用非傳統(tǒng)的截面三維翹曲的概念，進(jìn)一步發(fā)展了有限轉(zhuǎn)角旋轉(zhuǎn)梁理論[4]。2008年，王益峰[5]等人建立了一種復(fù)合材料柔性梁靜態(tài)響應(yīng)計(jì)算方法，采用中等變形梁理論，對(duì)梁結(jié)構(gòu)進(jìn)行了有限元分析。動(dòng)響應(yīng)常用的求解方法是模態(tài)疊加法，但是模態(tài)疊加法無(wú)法準(zhǔn)確計(jì)算高階模態(tài)。近年來(lái)，國(guó)內(nèi)外學(xué)者提出了通過直接的非線性數(shù)值積分求解旋翼結(jié)構(gòu)動(dòng)響應(yīng)的方法，主要方法有線性加速度法、Newmark法、Wilson-θ法、HHT(Hilber-Hughes-Taylor)法等。2018年，電子科技大學(xué)的石勝兵[6]將時(shí)域微分方程數(shù)值積分的經(jīng)典時(shí)域步進(jìn)算法——Newmark 算法引入到時(shí)域有限差分法中，利用中心差分離散空間偏微分，得到了一種基于Newmark算法的無(wú)條件穩(wěn)定的時(shí)域有限差分法，并將算法引入電磁領(lǐng)域，證明了該算法有效性和高效性。文獻(xiàn)調(diào)研表明，國(guó)內(nèi)外學(xué)者對(duì)幾何非線性中等變形梁理論及其靜態(tài)響應(yīng)進(jìn)行了一定的研究，但是對(duì)瞬態(tài)動(dòng)響應(yīng)的研究較少。Newmark算法具有高穩(wěn)定性和高精度的特點(diǎn)，對(duì)于求解中等變形梁的動(dòng)響應(yīng)問題具有較好的應(yīng)用價(jià)值。

本文在Hodges中等變形非線性梁理論的基礎(chǔ)上，運(yùn)用VABS算法得出復(fù)合材料梁的剖面特性，然后結(jié)合有限元法，引入Newmark算法，開展中等變形非線性復(fù)合材料梁的瞬態(tài)動(dòng)響應(yīng)計(jì)算研究。以某飛機(jī)機(jī)翼翼梁算例，得到了翼梁在初始垂直速度下的位移、速度和加速度響應(yīng)，并分析了附加集中質(zhì)量位置對(duì)計(jì)算結(jié)果的影響。

1 復(fù)合材料梁剖面特性計(jì)算模型

1.1 局部坐標(biāo)系與鋪層角定義

針對(duì)復(fù)合材料層合結(jié)構(gòu)，建立鋪層局部坐標(biāo)系，如圖1所示，其中x為梁剖面的坐標(biāo)系，e為沿材料主方向建立的坐標(biāo)系，y為單層局部坐標(biāo)系。定義θ1和θ3用于表征單層鋪層角，其中θ3為材料主坐標(biāo)系e與局部坐標(biāo)系y的夾角，θ1是鋪層坐標(biāo)系x與單層局部坐標(biāo)系y的夾角。

圖1 復(fù)合材料鋪層角度定義

1.2 剖面特性計(jì)算

對(duì)于簡(jiǎn)單形狀的各向同性材料梁剖面，可以采用材料力學(xué)方法得到梁剖面的質(zhì)量-剛度特性，而復(fù)合材料梁沿展向的剖面剛度特性和質(zhì)量特性比較復(fù)雜，本文采用有限元的方法進(jìn)行計(jì)算。

VABS[7]程序由Wenbin Yu等人研發(fā)，使用變分漸近數(shù)值方法能夠計(jì)算具有任意幾何和材料特性的梁的一維橫截面剛度常數(shù)，并進(jìn)行橫向剪切和Vlasov細(xì)化，并通過多次試驗(yàn)驗(yàn)證了該程序的準(zhǔn)確性[8]。本文基于Hypermesh軟件和VABS程序，運(yùn)用matlab軟件進(jìn)行文本處理，最終建立了一套剖面特性計(jì)算方法，剖面特性計(jì)算流程如圖2所示。

圖2 基于Hypermesh-VABS的剖面參數(shù)計(jì)算流程圖

2 幾何非線性梁模型

2.1 坐標(biāo)系及幾何非線性變形模型

圖3 梁?jiǎn)卧冃巫鴺?biāo)系

未變形坐標(biāo)系和變形坐標(biāo)系下，單位向量的轉(zhuǎn)換關(guān)系為：

(1)

(2)

(3)

圖4 歐拉角幾何關(guān)系示意圖

根據(jù)歐拉角變換關(guān)系，考慮梁在拉伸、彎曲、扭轉(zhuǎn)方向的耦合變形關(guān)系以及幾何變形非線性，將TDU簡(jiǎn)化到保留二階小量，可得式(4)。

(4)

2.2 漢密爾頓方程

對(duì)于非保守系統(tǒng)，漢密爾頓作用量的表達(dá)式為：

(5)

其中：δU、δT和δW分別表示應(yīng)變能變分、動(dòng)能變分和外力虛功變分，并且各個(gè)變分項(xiàng)由各單元的變分疊加得到，即：

其中：N表示單元總數(shù)；δUi、δTi和δWi分別表示第i個(gè)梁?jiǎn)卧膽?yīng)變能變分、動(dòng)能變分和外力虛功變分。

2.2.1 應(yīng)變能變分

保留二階小量的幾何非線性梁模型的正應(yīng)變和切應(yīng)變表達(dá)式如下：

(6)

進(jìn)而可推導(dǎo)第i個(gè)梁?jiǎn)卧膽?yīng)變能變分：

(7)

其中式(7)中，各變量表達(dá)式如下：

式中：EA表示梁的軸向拉伸剛度；EIy表示梁繞z軸的彎曲剛度；EIz表示梁繞y軸的彎曲剛度；GJ表示梁的扭轉(zhuǎn)剛度。其中包含“︿”的量表示與剖面的翹曲有關(guān)。

2.2.2 動(dòng)能變分

在未變形坐標(biāo)系中，選取未變形梁?jiǎn)卧先我庖稽c(diǎn)P(x,0,0)，發(fā)生變形后移動(dòng)到了p′(x+u,v,w)，并且包含該點(diǎn)的剖面相對(duì)于變形后彈性軸旋轉(zhuǎn)了θ1。

梁?jiǎn)卧先我庖稽c(diǎn)關(guān)于左側(cè)節(jié)點(diǎn)的速度可表達(dá)如式(8)。

(8)

速度及速度的變分表達(dá)式為：

(9)

(10)

第i個(gè)梁?jiǎn)卧膭?dòng)能表達(dá)式和動(dòng)能變分表達(dá)式分別為：

(11)

(12)

其中：ρs表示梁的線密度。

將式(9)及式(10)代入動(dòng)能變分表達(dá)式(12)中，然后通過分部積分得到動(dòng)能變分表達(dá)式為：

(13)

其中：

(14)

進(jìn)一步，將動(dòng)能變分表達(dá)式寫成關(guān)于六個(gè)自由度的表達(dá)形式，可得：

Tv′δv′+Tw′δw′+Tφδφ+TF)dx

(15)

其中：

2.2.3 外力虛功變分

針對(duì)本文求解的瞬態(tài)動(dòng)響應(yīng)問題，單元外力計(jì)入單元體力和節(jié)點(diǎn)六力素，可得外力虛功變分為：

(16)

2.3 有限元模型

2.3.1 單元離散

如圖5所示，選用2節(jié)點(diǎn)12自由度梁?jiǎn)卧Ｐ蛯?duì)復(fù)合材料梁進(jìn)行單元離散，每個(gè)節(jié)點(diǎn)包含6個(gè)自由度，分別為拉伸位移u，沿z軸的位移v，繞y軸的轉(zhuǎn)角v′，沿y軸的位移w，繞z軸的轉(zhuǎn)角w′，繞x軸的扭轉(zhuǎn)角φ。

圖5 2節(jié)點(diǎn)12自由度梁?jiǎn)卧Ｐ?/p>

單元節(jié)點(diǎn)的位移向量為：

(17)

單元形函數(shù)采用非線性插值，節(jié)點(diǎn)位移到單元位移的換算關(guān)系如式(18)。

(18)

其中：

Hu1=Hφ1=1-s

Hu2=Hφ2=s

Hv1=Hw1=2s3-3s2+1

Hv1′=Hw1′=li(s3-2s2+s)

Hv2=Hw2=-2s3+3s2

Hv2′=Hw2′=li(s3-s2)

s=xj/lj

根據(jù)漢密爾頓原理，將單元形函數(shù)和節(jié)點(diǎn)位移向量帶入，可得單元的漢密爾頓作用量為：

(19)

其中：[M]j表示第j個(gè)單元的質(zhì)量矩陣，[C]j表示第j個(gè)單元的阻尼矩陣，[K]j表示第j個(gè)單元的剛度矩陣。

2.3.2 結(jié)構(gòu)阻尼模型

幾何非線性梁振動(dòng)過程中，起動(dòng)能耗散作用的因素主要是結(jié)構(gòu)阻尼。結(jié)構(gòu)阻尼是指機(jī)械系統(tǒng)振動(dòng)時(shí)，其內(nèi)部產(chǎn)生交變的應(yīng)變，通過材料內(nèi)摩擦而消耗能力，是反映結(jié)構(gòu)體系振動(dòng)過程中能量耗散特征的參數(shù)。結(jié)構(gòu)振動(dòng)時(shí)耗能因素較多，但影響程度有所不同。文獻(xiàn)[10]中提出了四種阻尼模型的建模方法，分別為總體Rayleigh阻尼法，單元Rayleigh阻尼法，整體阻尼比法，單元阻尼比法。對(duì)于復(fù)合材料梁，各單元段鋪層可能存在差異，本文選擇單元Rayleigh阻尼模型。

Rayleigh阻尼用于單一阻尼特性的簡(jiǎn)單阻尼結(jié)構(gòu)體系中，各部分阻尼可用統(tǒng)一的阻尼參數(shù)代表，此種結(jié)構(gòu)中，取阻尼c=αm+βk，其中，α，β通過式(20)確定。

(20)

其中：ξm，ξn為給定兩個(gè)典型模態(tài)振型的阻尼比，工程應(yīng)用和研究中，可以通過模態(tài)試驗(yàn)測(cè)得。

2.3.3 矩陣組裝

按照有限元法，依據(jù)各單元間的公共節(jié)點(diǎn)關(guān)系組裝各個(gè)單元便可獲得總質(zhì)量矩陣、總阻尼矩陣、總剛度矩陣和總載荷向量。依據(jù)漢密爾頓原理，將總體質(zhì)量、阻尼和剛度矩陣帶入，便可得到整個(gè)復(fù)合材料梁的漢密爾頓作用量表達(dá)式：

(21)

由于實(shí)際的位移向量δqT為任意量，所以式(21)可以轉(zhuǎn)化為如下的有限元形式的復(fù)合材料梁運(yùn)動(dòng)方程：

(22)

其中：M為梁的總體質(zhì)量矩陣，C為梁的總體阻尼矩陣，K為梁的總體剛度矩陣，F(xiàn)為梁的外載荷總向量。

3 瞬態(tài)動(dòng)響應(yīng)求解模型

本文采用Newmark方法對(duì)時(shí)間進(jìn)行離散，得到復(fù)合材料梁的時(shí)域有限元方法。基于Newmark法的位移、速度、加速度滿足以下條件：

(23)

將公式(23)的第三式代入前兩式，可得：

(24)

將公式(24)進(jìn)一步改寫成:

(25)

(26)

其中：

b4=γb1Δt,b5=1+γb2Δt,b6=(1+γb3-γ)Δt

(27)

γ和β為Newmark算法參數(shù)，γ∈[0,1]，β∈[0,0.5]其中當(dāng)γ≥0.5，β≥0.25(0.5+γ2)時(shí)，Newmark法是無(wú)條件穩(wěn)定的，所以本文的Newmark參數(shù)取γ=0.5，β=0.25。

基于Newmark法的復(fù)合材料梁瞬態(tài)動(dòng)響應(yīng)計(jì)算流程如圖6所示。

圖6 Newmark瞬態(tài)動(dòng)響應(yīng)計(jì)算流程

4 算例計(jì)算及討論

以變截面混雜復(fù)合材料機(jī)翼盒型曲梁作為算例，玻璃纖維和碳纖維均為預(yù)浸編織布，翼梁鋪層角為[(0/90)G/(0/90)C/(±45)G/(±45)C/(0/90)G]，翼梁根部剖面幾何尺寸如圖7所示。翼梁長(zhǎng)度為2 m，無(wú)上反，沿拋物線曲線后掠，梁軸線函數(shù)為y=0.1x2，翼梁剖面梢根比為0.5，翼梁幾何模型如圖8所示。復(fù)合材料參數(shù)摘自文獻(xiàn)[11]。

圖7 算例復(fù)合材料曲梁根部剖面示意圖

圖8 算例復(fù)合材料曲梁三維模型

算例工況：初始垂向速度瞬態(tài)響應(yīng)分析，模擬飛機(jī)在應(yīng)急著陸過程中垂直下降率過大的工況。

本文計(jì)算下降率2 m/s工況下翼梁0.6 s內(nèi)的瞬態(tài)動(dòng)響應(yīng)，算得振動(dòng)響應(yīng)時(shí)間為0.493 5 s。翼尖節(jié)點(diǎn)加速度、速度、加速度響應(yīng)如圖9～圖11所示，響應(yīng)過程中翼尖最大位移為0.031 6 m，最大速度為4.441 m/s，最大加速度為1 561.4 m/s2。分析圖像可以看出，0～0.12 s內(nèi)響應(yīng)位移較大，響應(yīng)曲線表現(xiàn)出較強(qiáng)的非線性特征；0.12～0.6 s內(nèi)位移值較少，曲線特征為標(biāo)準(zhǔn)的正弦衰減曲線。

圖9 2 m/s下降率著陸0～0.6 s翼尖位移響應(yīng)曲線

圖10 2 m/s下降率著陸0～0.6 s翼尖速度響應(yīng)曲線

圖11 2 m/s下降率著陸0～0.6 s翼尖加速度響應(yīng)曲線

在飛機(jī)機(jī)翼上，會(huì)有發(fā)動(dòng)機(jī)、副油箱等集中質(zhì)量掛載點(diǎn)。本文在計(jì)算過程中考慮翼梁有附加集中質(zhì)量地對(duì)著陸響應(yīng)的影響，分別在機(jī)翼根部、中部和翼尖位置附加1 kg的集中載荷，初始動(dòng)能相同的情況下，計(jì)算得到機(jī)翼翼尖位移、速度和加速度動(dòng)響應(yīng)如圖12～圖14所示。

圖12 附加集中質(zhì)量位置對(duì)位移響應(yīng)的影響

圖13 附加集中質(zhì)量位置對(duì)速度響應(yīng)的影響

圖14 附加集中質(zhì)量位置對(duì)加速度響應(yīng)的影響

計(jì)算得到三種狀態(tài)下的動(dòng)能耗散時(shí)間、翼尖最大垂向位移、最大垂向速度和最大垂向加速度如表1所示。計(jì)算結(jié)果表明，附加質(zhì)量越靠近翼尖， 動(dòng)能耗散時(shí)間越長(zhǎng)，阻尼振動(dòng)響應(yīng)周期越長(zhǎng)，翼尖最大垂向位移響應(yīng)越大，翼尖最大垂向速度越小，翼尖最大垂向加速度越小。

表1 附加集中質(zhì)量位置對(duì)振動(dòng)響應(yīng)的影響結(jié)果匯總

根據(jù)表1中的數(shù)據(jù)，分析附加質(zhì)量位置對(duì)動(dòng)能耗散總時(shí)間和翼尖最大垂向位移兩個(gè)關(guān)鍵參數(shù)的影響，如圖15所示?？梢钥闯?，動(dòng)能耗散總時(shí)間和翼尖最大垂向位移均隨著附加質(zhì)量的展向位置的增加而單調(diào)增加，將集中質(zhì)量布置在翼根對(duì)降低動(dòng)響應(yīng)時(shí)間和位移幅值最有利。對(duì)比翼根和翼中附加質(zhì)量的計(jì)算結(jié)果，翼尖最大垂向位移增加了32.18%，動(dòng)能耗散時(shí)間增加了33.31%；對(duì)比翼中和翼尖附加質(zhì)量的計(jì)算結(jié)果，翼尖最大垂向位移增加了154.98%，動(dòng)能耗散時(shí)間增加了0.24%。

圖15 附加質(zhì)量位置對(duì)最大位移和耗散時(shí)間的影響對(duì)比

5 總結(jié)

本文將高長(zhǎng)細(xì)比復(fù)合材料梁分解成非線性一維梁模型和線性二維剖面特性分析兩部分，首先根據(jù)梁剖面的幾何特性和鋪層特征，運(yùn)用VABS算法得出梁剖面的質(zhì)量矩陣、剛度矩陣等，然后將基于Hodges中等變形梁模型，建立中等變形梁模型，將二維剖面特性參數(shù)帶入梁模型中。選用Newmark算法，完成一套基于VABS程序、Hodges中等變形梁模型和Newmark算法的幾何非線性復(fù)合材料梁瞬態(tài)動(dòng)響應(yīng)高效計(jì)算方法。最后以某飛機(jī)機(jī)翼翼梁為算例，得到了翼梁在初始垂直速度下的位移、速度和加速度響應(yīng)，并分析了附加集中質(zhì)量位置對(duì)計(jì)算結(jié)果的影響。

在工程應(yīng)用中，復(fù)合材料主承力梁的結(jié)構(gòu)設(shè)計(jì)思路往往是從二維剖面設(shè)計(jì)再到三維結(jié)構(gòu)設(shè)計(jì)，且初步設(shè)計(jì)階段并未將動(dòng)響應(yīng)考慮在內(nèi)。本文所研究的算法具有高效性和高精度性，對(duì)考慮動(dòng)響應(yīng)的復(fù)合材料梁的高效優(yōu)化迭代設(shè)計(jì)具有一定的應(yīng)用前景。通過對(duì)飛機(jī)機(jī)翼大梁響應(yīng)進(jìn)行算例分析，得出了以下結(jié)論：

1)中等變形復(fù)合材料梁在初始響應(yīng)階段，由于梁撓度較大，在位移、速度、加速度時(shí)域曲線中，表現(xiàn)出典型的幾何非線性響應(yīng)特性。

2)中等變形復(fù)合材料梁在響應(yīng)穩(wěn)定衰減階段，梁撓度變小，中等變形梁退化為小變形梁動(dòng)響應(yīng)問題，位移、速度、加速度時(shí)域曲線為標(biāo)準(zhǔn)正弦衰減型號(hào)。

3)附加集中質(zhì)量位置對(duì)復(fù)合材料中等變形梁的動(dòng)響應(yīng)指標(biāo)參數(shù)均有不同程度的影響，本文的研究?jī)?nèi)容可用于指導(dǎo)結(jié)構(gòu)方案布局設(shè)計(jì)。