大變形復(fù)合材料薄板多體系統(tǒng)動(dòng)力學(xué)建模

張煒華, 劉錦陽(上海交通大學(xué) 船舶海洋與建筑工程學(xué)院，上海　200240)

?

大變形復(fù)合材料薄板多體系統(tǒng)動(dòng)力學(xué)建模

張煒華, 劉錦陽(上海交通大學(xué) 船舶海洋與建筑工程學(xué)院，上海200240)

摘要：對(duì)大變形復(fù)合材料薄板的多體系統(tǒng)動(dòng)力學(xué)建模方法進(jìn)行研究。基于Kirchhoff假設(shè)，法線與中面保持垂直，從格林應(yīng)變的表達(dá)式出發(fā)，建立了面內(nèi)應(yīng)變和曲率與絕對(duì)位置坐標(biāo)和斜率的關(guān)系，在此基礎(chǔ)上推導(dǎo)了廣義彈性力陣和彈性力陣對(duì)廣義坐標(biāo)的導(dǎo)數(shù)陣，用絕對(duì)節(jié)點(diǎn)坐標(biāo)方法建立了大變形復(fù)合材料薄板多體系統(tǒng)的動(dòng)力學(xué)方程，用廣義α法和和牛頓迭代法求解微分-代數(shù)混合方程。對(duì)外載荷作用下的復(fù)合材料薄板進(jìn)行數(shù)值仿真，通過與ANSYS的仿真結(jié)果進(jìn)行對(duì)比，驗(yàn)證了該建模方法的準(zhǔn)確性和快速收斂性。最后，將建模方法應(yīng)用于復(fù)合材料太陽帆板展開機(jī)構(gòu)的數(shù)值仿真，分析了不同鋪層情況下驅(qū)動(dòng)力和約束力的振動(dòng)特性。

關(guān)鍵詞：大變形；復(fù)合材料薄板；動(dòng)力學(xué)；空間展開機(jī)構(gòu)

復(fù)合材料具有輕質(zhì)高強(qiáng)度和耐高溫等特點(diǎn)。在航空航天領(lǐng)域，大多數(shù)機(jī)械系統(tǒng)是由多個(gè)柔性復(fù)合材料部件組成的多體系統(tǒng)，如航天器太陽帆板、螺旋槳槳葉和航天器的旋轉(zhuǎn)機(jī)翼都可以看成復(fù)合材料板。每個(gè)柔性部件在外載荷的作用下都可能會(huì)產(chǎn)生較大的彈性變形，因此，考慮幾何非線性的大變形復(fù)合材料柔性板的建模具有挑戰(zhàn)性，結(jié)構(gòu)大變形這一幾何非線性問題已成為力學(xué)、航空航天和機(jī)械等領(lǐng)域研究的難點(diǎn)問題之一。

混合坐標(biāo)法[1]采用剛體大范圍運(yùn)動(dòng)坐標(biāo)和彈性變形坐標(biāo)為廣義坐標(biāo)，建立剛-柔耦合的動(dòng)力學(xué)模型。混合坐標(biāo)法的優(yōu)點(diǎn)是可以在小變形的情況下通過模態(tài)截?cái)喾s小廣義坐標(biāo)的規(guī)模，但是在大變形情況下，模態(tài)截?cái)喾ú辉龠m用，質(zhì)量陣和剛度陣都含有與廣義坐標(biāo)相關(guān)的高次項(xiàng)，導(dǎo)致計(jì)算效率較低。

Shabana提出的絕對(duì)節(jié)點(diǎn)坐標(biāo)法已廣泛應(yīng)用于大變形柔性多體系統(tǒng)動(dòng)力學(xué)分析，該方法的特點(diǎn)是選取柔性梁或柔性板各單元節(jié)點(diǎn)相對(duì)慣性基的位移和斜率作為廣義坐標(biāo)，建立動(dòng)力學(xué)方程，這樣就不需要引入大范圍運(yùn)動(dòng)變量和變量位移坐標(biāo)，廣義坐標(biāo)全部是全局坐標(biāo)，得到的質(zhì)量陣為常值陣，廣義力的表達(dá)式也比較簡(jiǎn)單。Berzeri等[2]基于絕對(duì)節(jié)點(diǎn)坐標(biāo)法，提出了平面梁的三種幾何非線性動(dòng)力學(xué)模型，研究了各種模型對(duì)大變形問題的適用性，通過大量的數(shù)值計(jì)算發(fā)現(xiàn)，縱向應(yīng)變和曲率的精確計(jì)算是保證計(jì)算精度的關(guān)鍵。在此基礎(chǔ)上，Berzeri等[3]用絕對(duì)節(jié)點(diǎn)坐標(biāo)法研究了動(dòng)力剛化問題，通過數(shù)值對(duì)比驗(yàn)證了該方法的正確性，Omar等[4]進(jìn)一步考慮剪切效應(yīng)，建立了二維梁的有限元模型。在平面梁的基礎(chǔ)上，Dombrowski[5]將絕對(duì)節(jié)點(diǎn)坐標(biāo)法推廣到大變形的空間梁。近幾年來，Dombrowski等[6]用絕對(duì)節(jié)點(diǎn)坐標(biāo)法對(duì)大變形薄板的動(dòng)力學(xué)建模理論進(jìn)行研究，建立了矩形薄板的二維有限元模型，鄒凡等[7]將大變形薄板的動(dòng)力學(xué)建模理論研究推廣到多體系統(tǒng)。Mikkola等[8]進(jìn)一步考慮剪切變形，建立了三維矩形薄板的有限元模型。為了驗(yàn)證絕對(duì)節(jié)點(diǎn)坐標(biāo)法的正確性，Yoo等[9-11]首先用快速相機(jī)對(duì)大變形薄板進(jìn)行實(shí)驗(yàn)研究，并將測(cè)量技術(shù)推廣到大變形復(fù)雜結(jié)構(gòu)多體系統(tǒng)。劉錦陽等[12]對(duì)氣浮臺(tái)-大變形薄板多體系統(tǒng)開展了實(shí)驗(yàn)研究。為了提高絕對(duì)節(jié)點(diǎn)坐標(biāo)法的計(jì)算效率，劉鋮等[13]采用第一類Piola-Kirchhoff應(yīng)力張量的方法推導(dǎo)了彈性力和Jacobi矩陣的解析表達(dá)式，并驗(yàn)證了理論模型的準(zhǔn)確性。趙將等[14]針對(duì)簡(jiǎn)化的“IKAROS” 自旋展開太陽帆系統(tǒng)，采用結(jié)合自然坐標(biāo)方法與絕對(duì)節(jié)點(diǎn)坐標(biāo)方法的絕對(duì)坐標(biāo)方法對(duì)其進(jìn)行建模，研究了黏彈性太陽帆薄膜自旋展開過程的動(dòng)力學(xué)特性。以上工作的研究對(duì)象均為各向同性材料的梁和板。在此基礎(chǔ)上，Liu等[15]用絕對(duì)節(jié)點(diǎn)坐標(biāo)法建立了復(fù)合材料三維矩形板的動(dòng)力學(xué)模型，但由于直接采用連續(xù)介質(zhì)力學(xué)方法描述應(yīng)變能，容易產(chǎn)生泊松鎖定現(xiàn)象。

本文研究大變形復(fù)合材料薄板的多體系統(tǒng)動(dòng)力學(xué)建模方法。不計(jì)復(fù)合材料薄板的面外剪切變形，基于Kirchhoff假設(shè)，從格林應(yīng)變和曲率與絕對(duì)位移的關(guān)系式出發(fā)，利用復(fù)合材料板的本構(gòu)關(guān)系，推導(dǎo)了廣義彈性力陣，在此基礎(chǔ)上用絕對(duì)節(jié)點(diǎn)坐標(biāo)法建立了大變形復(fù)合材料薄板的有限元離散的動(dòng)力學(xué)方程。為了提高數(shù)值仿真效率，用廣義-α法和牛頓迭代法求解微分-代數(shù)混合方程，推導(dǎo)了彈性力陣對(duì)廣義坐標(biāo)的導(dǎo)數(shù)陣。分別對(duì)懸臂板和鉸支板進(jìn)行數(shù)值仿真，將本文方法的仿真結(jié)果與ANSYS非線性有限元的計(jì)算結(jié)果進(jìn)行對(duì)比，驗(yàn)證了本文建模方法的準(zhǔn)確性和單元收斂性，最后將建模方法推廣到大變形柔性多體系統(tǒng)，數(shù)值模擬了復(fù)合材料太陽帆板的展開過程，得到了對(duì)工程實(shí)際具有參考價(jià)值的結(jié)論。

1薄板應(yīng)變和曲率的推導(dǎo)

為了表述簡(jiǎn)單，對(duì)于任意列向量a=[a1a2a3]T，本節(jié)采用如下表述：

(1)

大變形復(fù)合材料矩形薄板如圖1所示。

圖1　變形前后的復(fù)合材料薄板Fig.1 Composite thin plate before and after deformation

設(shè)板的長度為a，寬度為b，厚度為h，體密度為ρ。建立固結(jié)在地面的慣性坐標(biāo)系O-XYZ。設(shè)法線n0(x,y,t)為變形后沿法線的單位矢量，基于Kirchhoff假設(shè)，法線n0(x,y,t)與中面保持垂直。用有限元方法對(duì)矩形薄板進(jìn)行離散。將矩形薄板等分為n=n1×n2個(gè)單元，對(duì)第e個(gè)單元建立單元坐標(biāo)系Oe-XeYeZe，則板單元的長度和寬度分別為ae=a/n1，be=b/n2。設(shè)(x,y,0)為薄板中面上任意一點(diǎn)在該單元坐標(biāo)系下的坐標(biāo)，薄板非中面上任意一點(diǎn)k的絕對(duì)位置矢量r在慣性坐標(biāo)系下的坐標(biāo)陣為：

r(x,y,z,t)=r0(x,y,t)+zn0(x,y,t)

(2)

式中，r0為中面上對(duì)應(yīng)點(diǎn)k0的絕對(duì)位置矢量在慣性坐標(biāo)系下的坐標(biāo)陣。

格林應(yīng)變列陣與絕對(duì)位移的關(guān)系為：

(3)

將式(2)代入式(3)，不計(jì)與z2有關(guān)的項(xiàng)，得到：

(4)

考慮到法線n0與中面保持垂直，得到如下關(guān)系式：

r0xTn0=0，r0yTn0=0

(5)

式(5)關(guān)于x求偏導(dǎo)：

r0xTn0x+r0xxTn0=0，r0yTn0x+r0yxTn0=0

(6)

式(5)關(guān)于y求偏導(dǎo)：

r0xTn0y+r0xyTn0=0，r0yTn0y+r0yyTn0=0

(7)

將式(6)和式(7)代入式(4)，并利用關(guān)系式r0xy=r0yx，格林應(yīng)變的表達(dá)式可寫成：

ε=ε0-zκ

(8)

式中，ε0為中面的面內(nèi)應(yīng)變列陣，可表示為：

(9)

κ為中面的曲率列陣，可表示為：

(10)

2質(zhì)量陣和廣義外力陣的推導(dǎo)

基于絕對(duì)節(jié)點(diǎn)坐標(biāo)法，中面上對(duì)應(yīng)點(diǎn)k0的絕對(duì)位置矢量在慣性坐標(biāo)系下的坐標(biāo)陣為：

(11)

式中，S(x,y)為形函數(shù)矩陣，由文獻(xiàn)[6]給出。qe(t)為單元節(jié)點(diǎn)的絕對(duì)位置坐標(biāo)，位置坐標(biāo)分別對(duì)x,y的一階導(dǎo)數(shù)，以及位置坐標(biāo)對(duì)x和y的兩階導(dǎo)數(shù)在慣性基下的坐標(biāo)陣，qe(t)表達(dá)式分別為

qe(t)=[q1e(t)Tq2e(t)Tq3e(t)Tq4e(t)T]T

(12)

式中，qke(t),k=1,…,4為單元4個(gè)節(jié)點(diǎn)的絕對(duì)節(jié)點(diǎn)坐標(biāo)列陣，具體表述方式如下：

qke(t)=[rk(t)Trkx(t)Trky(t)Trkxy(t)T]T,

k=1,…,4

(13)

設(shè)q(t)為總體絕對(duì)節(jié)點(diǎn)坐標(biāo)列陣，qe(t)與q(t)的關(guān)系為qe(t)=Beq(t)，其中，Be為單元e的布爾陣。

針對(duì)薄板問題，在計(jì)算慣性力的虛功時(shí)可以忽略絕對(duì)位移沿單元坐標(biāo)系Ze方向的變化，板單元e的慣性力做的虛功為：

(14)

將式(11)代入式(14)：

(15)

式中，單元質(zhì)量陣Me為：

(16)

矩形薄板的慣性力做的虛功為：

(17)

在體積力只考慮重力的情況下，設(shè)重力沿Z的負(fù)向，單位體積力在慣性基下的坐標(biāo)列陣為：

f=[00-ρg]T

(18)

利用關(guān)系式(11)，板單元的體力做的虛功為：

δWge=∫VeδrTfdVe=δqeTQge

(19)

式中，單元體力列陣Qge為：

(20)

體力做的虛功為：

(21)

設(shè)集中力作用于中面上的P點(diǎn)，該集中力在慣性坐標(biāo)系下的坐標(biāo)列陣為F，P點(diǎn)的絕對(duì)位置矢量在慣性坐標(biāo)系下的坐標(biāo)列陣與絕對(duì)節(jié)點(diǎn)坐標(biāo)列陣的關(guān)系為r0P=SPq，集中力做的虛功為：

δWF=δr0PTF=δqTQF

(22)

式中，集中力對(duì)應(yīng)的廣義力陣為：QF=SPTF。

3廣義彈性力陣的推導(dǎo)

對(duì)中面的面內(nèi)應(yīng)變求變分，并利用關(guān)系式(11)，得到：

δε0=Hδqe

(23)

其中：

(24)

基于Kirchhoff假設(shè)，中面的單位法向量在慣性基下的坐標(biāo)陣n0可表示為：

(25)

式中，n為rx與ry的叉積，n為矢量n的模，可表示為

(26)

(27)

對(duì)式(26)求變分，利用式(11)，得到：

(28)

其中：

(29)

對(duì)式(27)求變分：

(30)

將式(25)代入式(10)，得到κ的表達(dá)式為：

(31)

對(duì)式(31)求變分，并利用關(guān)系式(11)，得到：

δκ=Gδqe

(32)

其中：

G=G1+G2+G3

(33)

本文研究的層合板由N層復(fù)合材料板粘合在一起組成的整體結(jié)構(gòu)板，第k層板的變形如圖2所示，其中，Γk為纖維角。

第k層板的應(yīng)力σ=[σ11σ22σ12]T與應(yīng)變的本構(gòu)關(guān)系可表示為[16]：

σ=Dkε

(34)

式中：

Dk=Tk-1DTk-T

(35)

(36)

(37)

式中，E1和E2是沿主軸方向的彈性模量，G12為面內(nèi)剪切模量，v21,v12為泊松比，滿足關(guān)系式v21E1=v12E2。

圖2　復(fù)合材料板面內(nèi)變形Fig.2 In-plane deformation of the composite plate

板單元彈性力做的虛功為：

(38)

其中：

將式(8)和(34)代入式(39)可得第k層層合板彈性力做的虛功的表達(dá)式為：

(40)

將式(40)代入式(38)并化簡(jiǎn)得到N層復(fù)合材料板的單元彈性力虛功表達(dá)式：

δWfe=δWf1+δWf2+δWf3+δWf4

(41)

其中，At,Bt,Ct的具體表達(dá)式為：

(42)

將式(23)和(32)代入式(41)可得：

δWfe=δqeTQfe

(43)

其中，單元廣義彈性力陣的具體表達(dá)式為：

Qfe=Qfe1+Qfe+Qfe3+Qfe4

(44)

彈性力做的虛功為：

(45)

4復(fù)合材料薄板動(dòng)力學(xué)方程

根據(jù)虛功原理，動(dòng)力學(xué)變分方程為：

δWi+δWg+δWF+δWf=0

(46)

將式(17)，(21)，(22)和(45)代入式(46)，復(fù)合材料層合薄板的動(dòng)力學(xué)變分方程可寫為：

(47)

對(duì)于由NB個(gè)復(fù)合材料薄板組成的多體系統(tǒng)，設(shè)系統(tǒng)廣義坐標(biāo)陣為q，第i個(gè)柔性體的廣義坐標(biāo)陣為qi，多體系統(tǒng)的動(dòng)力學(xué)變分方程為：

(48)

或?qū)懗桑?/p>

(49)

其中，M=diag(M1,…，MNB)，Qg=[Qg1T，…，QgNBT]T，QF=[QF1T，…，QFNBT]T，Qf=[Qf1T，…，QfNBT]T。

設(shè)柔性多體系統(tǒng)的約束方程為：

Φ(q,t)=0

(50)

含拉格朗日乘子的第一類拉格朗日動(dòng)力學(xué)方程為：

(51)

式中，Qc=-ΦqTλ為廣義約束力，Φq為約束方程的雅可比矩陣，λ為由拉格朗日乘子構(gòu)成的列陣。

將式(51)和式(50)聯(lián)立構(gòu)成微分代數(shù)混合方程：

(52)

其中：

(53)

本文采用廣義-α法結(jié)合牛頓-拉斐遜方法進(jìn)行數(shù)值仿真，以提高計(jì)算效率。數(shù)值方法的詳細(xì)計(jì)算步驟可參考文獻(xiàn)[14]。

5廣義彈性力陣對(duì)廣義坐標(biāo)的導(dǎo)數(shù)陣的推導(dǎo)

下面詳細(xì)介紹導(dǎo)數(shù)陣?Qf/?q的計(jì)算方法。結(jié)合式(44)，可將Qfe分成四項(xiàng)分別進(jìn)行推導(dǎo)。

(54)

其中，Qfe的第一項(xiàng)可表示為：

(55)

其中：

J1=Atε0=[J11J12J13]T

(56)

Qfe的第二項(xiàng)可表示為：

(57)

其中：

J2=Btκ=[J21J22J23]T

(58)

Qfe的第三項(xiàng)可表示為：

(59)

其中：

J3=Btε0=[J21J22J23]T

(60)

Vc中每一項(xiàng)的表達(dá)式為：

(61)

(62)

(63)

Qfe的第四項(xiàng)可表示為：

(64)

其中：

J4=Ctκ=[J41J42J43]T

(65)

Vd中每一項(xiàng)的表達(dá)式為：

(66)

(67)

(68)

結(jié)合qe=Beq和式(54)可得到：

(69)

6仿真計(jì)算

6.1外載荷作用下矩形薄板的數(shù)值仿真

本文驗(yàn)證所用的大變形復(fù)合材料矩形薄板的各項(xiàng)參數(shù)如下：長度為l=1.5 m，寬度b=2.0 m，厚度h=1.5×10-3m，體密度為ρ=1 560 kg/m3，彈性模量為E1=181 GPa，E2=103 GPa，G12=7.7 GPa，泊松比為v12=0.28，v21=v12E2/E1，復(fù)合材料板纖維角分布為Γ=[0°90°0°]T。

分別對(duì)下面兩種情況進(jìn)行動(dòng)力學(xué)分析。第一種情況為受載荷作用的懸臂-自由板，如圖3所示。以O(shè)為原點(diǎn)建立慣性坐標(biāo)系XYZ，矩形薄板一條邊固支，在另一條邊的中點(diǎn)C處施加一個(gè)沿Z軸負(fù)向，大小分別為10 N和20 N的恒力F。將本文模型分別與ANSYS的大變形模型和小變形模型進(jìn)行比較，以驗(yàn)證本文模型的準(zhǔn)確性。ANSYS驗(yàn)證單元取SHELL281單元。

圖3　大變形復(fù)合材料懸臂板Fig.3 Cantilevered composite plate with large deformation

圖4(a)和(b)為當(dāng)F的大小分別等于10 N和20 N時(shí)，用本文模型、ANSYS大變形模型和小變形模型仿真得到的C點(diǎn)的撓度曲線。從圖中不難看出，當(dāng)變形超過柔性部件本身尺寸的10%時(shí)，小變形理論已經(jīng)不再適用，而這時(shí)，本文模型和ANSYS大變形模型得到的結(jié)果一致，驗(yàn)證了本文模型的準(zhǔn)確性。而通過比較兩種不同大小的力的作用，可以發(fā)現(xiàn)懸臂板的振動(dòng)頻率并沒有變化。這也和實(shí)際情況相符。

圖4　C點(diǎn)撓度Fig.4 Deflection of point C

第二種情況為受載荷作用的鉸支板，如圖5所示。以O(shè)為原點(diǎn)建立慣性坐標(biāo)系XYZ，矩形薄板四個(gè)角點(diǎn)均通過球鉸和固定支座相連接，并在板中心C處施加一個(gè)沿Z軸負(fù)向，大小分別為2 N和10 N的恒力F。將本文模型和ANSYS大變形理論進(jìn)行了收斂速度的比較。ANSYS驗(yàn)證單元取SHELL281單元。

圖5　大變形復(fù)合材料鉸支板Fig.5 Simply-supported composite plate with large deformation

圖6(a)～(c)為鉸支板在2 N力的作用下的撓度曲線，圖7為10 N力作用下的撓度曲線。通過比較圖6(a)～(c)不難看出，用本文模型，單塊板僅需36個(gè)單元就收斂，而收斂到同樣精度ANSYS大變形模型需要64個(gè)單元。通過比較圖6和圖7不難發(fā)現(xiàn)，和懸臂板不同的是鉸支板的振動(dòng)頻率會(huì)隨著作用力增大而增高，這是因?yàn)樗狞c(diǎn)鉸支，增大作用力會(huì)增加面內(nèi)應(yīng)力，從而引起剛化現(xiàn)象。從圖中可以看到，用本文方法收斂到相同精度所需要的單元數(shù)少于ANSYS大變形模型，體現(xiàn)了本文方法的高效性。

6.2復(fù)合材料太陽帆板展開數(shù)值仿真

本算例所用的太陽帆板為四塊具有相同尺寸和相同物理參數(shù)的復(fù)合材料薄板，板的參數(shù)如下：長度為l=1.5 m，寬度為b=2.0 m，厚度為h=1.5×10-3m，體密度為ρ=1 560 kg/m3，彈性模量為E1=181 GPa，E2=103 GPa，G12=7.7 GPa，泊松比為v12=0.28，v21=v12E2/E1。每塊板劃分為6×6個(gè)單元。

圖6　2 N力作用下C點(diǎn)撓度Fig.6 Deflection of point C under the force of 2 N

圖7　10 N力作用下C點(diǎn)撓度Fig7 Deflection of point C under the force of 10 N

下面對(duì)以下三種情況進(jìn)行數(shù)值仿真。第一種情況是直接將薄板簡(jiǎn)化為剛體模型進(jìn)行仿真，第二種情況復(fù)合材料板纖維角為對(duì)稱分布：

Γ=[-45°45°-45°45°90°0°0°90°45°-45°45°-45°]T

第三種情況復(fù)合材料纖維角為反對(duì)稱分布：

Γ=[-45°45°-45°45°90°0°0°-90°-45°45°-45°45°]T

圖8　復(fù)合材料太陽帆板Fig.8 Composite solar plate

帆板展開過程中的1，2號(hào)節(jié)點(diǎn)作用于帆板的Y方向約束力，以及Z方向的約束力和驅(qū)動(dòng)力如圖9～圖11所示：

圖9　Y方向約束力Fig.9 Constrained force along Y axis

圖10　Z方向約束力Fig.10 Constrained force along Z axis

圖11　驅(qū)動(dòng)力Fig.11 Driving force

從上圖不難看出，在展開過程中太陽帆板由于彈性變形導(dǎo)致約束力和驅(qū)動(dòng)力產(chǎn)生明顯的震蕩，但是平均值和剛體模型仿真計(jì)算結(jié)果一致。可以看出，纖維角為對(duì)稱分布和反對(duì)稱分布的復(fù)合材料板沿Y和Z方向的約束力基本重合，而沿X方向的驅(qū)動(dòng)力則在減速時(shí)有較大的區(qū)別，與對(duì)稱分布的復(fù)合材料板相比，反對(duì)稱復(fù)合材料板在受壓時(shí)沿X方向的驅(qū)動(dòng)力會(huì)有更為明顯的震蕩。此外可以看出，復(fù)合材料板的節(jié)點(diǎn)1和節(jié)點(diǎn)2沿Y方向的約束力大小相等，方向相反，震蕩幅值較大，不容忽視。

7結(jié)論

目前，大變形復(fù)合材料薄板建模的最大難題是廣義彈性力陣和關(guān)于廣義坐標(biāo)的導(dǎo)數(shù)陣的準(zhǔn)確推導(dǎo)。本文基于Kirchhoff假設(shè)，從格林應(yīng)變與位移的關(guān)系式出發(fā)，建立了面內(nèi)應(yīng)變和曲率與絕對(duì)節(jié)點(diǎn)坐標(biāo)的關(guān)系式，在此基礎(chǔ)上，推導(dǎo)了廣義彈性力陣，用絕對(duì)節(jié)點(diǎn)坐標(biāo)法建立了大變形復(fù)合材料薄板的動(dòng)力學(xué)方程組，結(jié)合隱式廣義-α法和牛頓迭代法求解。為了提高牛頓迭代法的收斂效率，推導(dǎo)了廣義彈性力陣對(duì)廣義坐標(biāo)的導(dǎo)數(shù)矩陣的精確表達(dá)式，從而解決了長期存在的導(dǎo)數(shù)矩陣計(jì)算不準(zhǔn)確引起的數(shù)值發(fā)散或計(jì)算效率低的問題。通過數(shù)值仿真得到了以下結(jié)論：

(1) 本文建立的復(fù)合材料薄板的動(dòng)力學(xué)模型適用于大變形問題。與ANSYS的仿真結(jié)果進(jìn)行對(duì)比表明，用本文方法收斂到相同精度所需要的單元數(shù)少于ANSYS大變形模型，體現(xiàn)了本文方法的高效性。

(2) 將建模方法成功應(yīng)用于復(fù)合材料太陽帆板展開機(jī)構(gòu)的數(shù)值仿真，對(duì)不同鋪層情況下的帆板的數(shù)值仿真表明，在帆板展開過程中，和同性材料板不同的是，復(fù)合材料薄板沿轉(zhuǎn)動(dòng)軸方向的約束力震蕩顯著。對(duì)稱分布板和反對(duì)稱分布板在加速階段驅(qū)動(dòng)力基本吻合，但是在減速階段反對(duì)稱板的驅(qū)動(dòng)力會(huì)有更為顯著的震蕩。

參 考 文 獻(xiàn)

[ 1 ] Likins P W. Finite element appendage equations for hybrid coordinate dynamic analysis[J].Journal of Solids and Structures, 1972, 8：790-831.

[ 2 ] Berzeri M, Shabana A A. Development of simple models for the elastic forces in absolute nodal co-ordinate formulation[J]. Journal of Sound and Vibration, 2000, 235(4)：539-565.

[ 3 ] Berzeri M, Shabana A A. Study of the centrifugal stiffening effect using the finite element absolute nodal coordinate formulation[J].Multibody System Dynamics,2002(7)：357-387.

[ 4 ] Omar M A, Shabana A A. A two-dimensional shear deformation beam for large rotation and deformation problems[J].Journal of Sound and Vibration, 2001, 243(3)： 565-576.

[ 5 ] Dombrowski S V. Analysis of large flexible body deformation in multibody systems using absolute coordinates[J].Multibody System Dynamics,2002(8)： 409-432.

[ 6 ] Dmitrochenko O N,Pogorelov D Yu. Generalization of plate finite elements for absolute nodal coordinate formulation[J].Multibody System Dynamics,2003(10)：17-43.

[ 7 ] 鄒凡，劉錦陽.大變形薄板多體系統(tǒng)的動(dòng)力學(xué)建模[J].應(yīng)用力學(xué)學(xué)報(bào)，2010，27(4)：740-745.

ZOU Fan, LIU Jin-yang. Dynamic modeling for thin plate multibody system with large deformation[J].Chinese Journal of Applied Mechanics, 2010, 27(4)：740-745.

[ 8 ] Mikkola A M, Shabana A A. A non-incremental finite element procedure for the analysis of large deformation of plates and shells in mechanical system applications[J].Multibody System Dynamics,2003(9)：283-309.

[ 9 ] Yoo W S, Lee J H, Park S J, et al. Large deflection analysis of a thin plate: computer simulations and experiments[J].Multibody System Dynamics,2004(11)：185-208.

[10] Yoo W S, Kim M S, Mun S H, et al. Large displacement of beam with base motion：flexible multibody simulations and experiments[J].Computing Methods Applied Mechanical Engineering,2006, 195：7036-7051.

[11] Yoo W S, Kim M S, Mun S H, et al. Developments of multibody system dynamics：computer simulations and experiments[J].Multibody System Dynamics,2007(18)：35-58.

[12] 劉錦陽，鄒凡，余征躍.氣浮臺(tái)-薄板剛-柔耦合系統(tǒng)動(dòng)力學(xué)建模和實(shí)驗(yàn)研究[J].振動(dòng)與沖擊，2012,31(21)：11-21.

LIU Jin-yang, ZOU Fan, YU Zheng-yue.Dynamic modeling and tests for a hub-plate rigid-flexible coupled system[J].Journal of Vibration and Shock, 2012, 31(21)：11-21.

[13] 劉鋮,田強(qiáng),胡海巖. 基于絕對(duì)節(jié)點(diǎn)坐標(biāo)法的多柔體系統(tǒng)動(dòng)力學(xué)高效計(jì)算方法 [J].力學(xué)學(xué)報(bào)，2010，42(2)：1197-1205.

LIU Cheng, TIAN Qiang, HU Hai-yan. Efficient computational method for dynamics of flexible multibody systems based on absolute nodal coordinate[J].Chinese Journal of Theoretical and Applied Mechanics,2010,45(2)：1197-1205.

[14] 趙將,劉鋮,田強(qiáng),等.黏彈性薄膜太陽帆自旋展開動(dòng)力學(xué)分析 [J].力學(xué)學(xué)報(bào)，2013，45(2)：746-754.

ZHAO Jiang, LIU Cheng, TIAN Qiang, et al. Dynamic analysis of spinning deployment of a solar sail composed of viscoelastic membranes[J].Chinese Journal of Theoretical and Applied Mechanics, 2013, 45(2)：746-754.

[15] Liu C, Tian Q，Hu H Y. Dynamics of a large scale rigid-flexible multibody system composed of composite laminated plates[J].Multibody System Dynamics, 2011，26(3)：283-305.

[16] Qatu M S. Vibration of laminated shells and plates [M]. San Diego：Elsevier, 2004： 7-12.

Dynamic modeling of composite thin-plate multibody systems with large deformation

ZHANGWei-hua,LIUJin-yang(School of Naval architecture, Ocean and Civil Engineering, Shanghai Jiao Tong University, Shanghai 200240, China)

Abstract:Dynamic modeling theory of composite thin-plate multibody systems with large deformation was investigated.Based on Kirchhoff’s assumption that the normal vector is always perpendicular to the central surface, the relation among in-plane strains, curvatures, absolute nodal coordinates and absolute gradients were derived according to the definition of Green strain, and the generalized elastic-force vector and its differentiation with respect to the generalized coordinates were derived.Equations of motion of the composite thin-plate multibody systems with large deformation were derived based on absolute nodal-coordinate formulation.Generalized method and the Newton-Raphson method were used to solve the differential-algebraic equations.Simulation of a composite thin plate applied with an external force was conducted.Comparison of the present simulation results with those obtained by ANSYS software verifies the accuracy and effectiveness of the formulation.Finally, the proposed formulation is used for numerical simulation of composite solar-array deployment mechanisms.The vibration characteristics of the driving force and constraint forces are analyzed in case of different panel layers.

Key words:large deformation; composite-laminated plate; dynamics; deployment mechanism

中圖分類號(hào):O313.7

文獻(xiàn)標(biāo)志碼：A

DOI:10.13465/j.cnki.jvs.2016.08.005

通信作者劉錦陽 女，教授，博士生導(dǎo)師，1964年9月生

收稿日期：2015-02-02修改稿收到日期：2015-04-25

基金項(xiàng)目：國家自然科學(xué)基金資助項(xiàng)目(11272203;11132007)

第一作者 張煒華 男，碩士生，1991年8月生

E-mail:liujy@sjtu.edu.cn