精確Cosserat彈性桿動力學(xué)的分析力學(xué)方法*

薛紜 翁德瑋 陳立群

1)(上海應(yīng)用技術(shù)學(xué)院機械工程學(xué)院,上海 201418)

2)(上海大學(xué),上海市應(yīng)用數(shù)學(xué)和力學(xué)研究所,上海 200072)

(2012年7月11日收到;2012年8月29日收到修改稿)

1 引言

彈性桿力學(xué)具有悠久的歷史和廣泛的應(yīng)用背景.20世紀(jì)70年代,彈性桿被作為脫氧核糖核酸(DNA)等生物大分子的力學(xué)模型研究其位形的拓撲性質(zhì)、平衡及其穩(wěn)定性,進一步推動了彈性桿力學(xué)的發(fā)展[1?4].鑒于應(yīng)用背景的復(fù)雜和多樣性,彈性桿力學(xué)面臨新的問題:其一是極端細長導(dǎo)致小應(yīng)變累積成超大位移,形成異常復(fù)雜的幾何形態(tài);其二是全新的約束形式,如自接觸等,對組蛋白的纏繞包括了經(jīng)典約束具有的單面、非定常、非完整等特征[2,5].因分析力學(xué)的建模方法在約束處理和數(shù)值計算上的優(yōu)勢[6,7],文獻[8—12]建立了Kirchhoff彈性桿靜力學(xué)和動力學(xué)的分析力學(xué)理論框架,文獻[13,14]將近代分析力學(xué)中的對稱性理論引入彈性桿靜力學(xué),得到了由對稱性導(dǎo)致的守恒量.

Kirchhoff彈性桿靜力學(xué)和動力學(xué)的分析力學(xué)方法已有所研究,由于忽略了彈性桿截面因變形導(dǎo)致隨弧坐標(biāo)的平移,即拉壓和剪切變形,使得變分原理中與主矢有關(guān)的項以及關(guān)于主矢的3個動力學(xué)方程不具有分析力學(xué)的形式,并且因主矢與主矩的耦合使得主矩方程的表達也不盡如人意[8?12],給應(yīng)用帶來了不便.在彎扭基礎(chǔ)上再計及軸向拉壓和截面剪切變形的彈性桿模型是Kirchhoff模型的改進[15,16],這被稱之為Cosserat彈性桿或精確模型,Kirchhoff彈性桿是其特例.

本文引入了關(guān)于主矢的本構(gòu)方程后形成的精確Cosserat彈性桿的分析力學(xué)方法,其變分原理和動力學(xué)方程完全具有分析力學(xué)的形式,這既有助于充分利用分析力學(xué)的現(xiàn)有成果為彈性桿力學(xué)提供新方法,又有利于拓展分析力學(xué)的應(yīng)用領(lǐng)域.

精確Cosserat彈性桿運動和變形的幾何關(guān)系見文獻[16,17].本文基于位形空間的虛位移概念,建立彈性桿精確模型的d’Alembert-Lagrange原理和Hamilton原理,在線性本構(gòu)關(guān)系下導(dǎo)出完全具有分析力學(xué)形式的各種運動微分方程,且以弧坐標(biāo)和時間雙自變量為特征.為利用分析力學(xué)的現(xiàn)有成果,如對稱性和守恒量等以及在數(shù)值仿真時實施辛算法鋪平道路.

2 彈性桿精確模型的位形描述

依據(jù)平面截面假定,以彈性桿的截面為對象,因變形導(dǎo)致截面隨弧坐標(biāo)發(fā)生轉(zhuǎn)動和斜向平移,如圖1所示.建立慣性坐標(biāo)系(O-ξηζ)和與截面固結(jié)的形心主軸坐標(biāo)系(P-xyz),沿坐標(biāo)軸的單位基矢量列陣分別為和ep=(e1(t)e2(,t)e3(,t))T,其中=(s,t)為變形后的弧坐標(biāo),s為未變形時的弧坐標(biāo),e3為截面的外法矢,指向弧坐標(biāo)的正向.兩組基有關(guān)系ep=Q eI,Q為單位正交陣.(P-xyz)的位置和姿態(tài)用截面形心相對慣性系的矢徑r=O P的坐標(biāo)陣和Euler角列陣描述,r,qE的坐標(biāo)依次用 qα,(α=1,···,6)表示.截面的運動方程為

圖1 彈性桿截面的平移

對于給定的姿態(tài) qα=qα(s,t),(α=4,5,6),截面隨弧坐標(biāo)的平移存在如下關(guān)系:?r=e3?s+?w(見圖 1),r(q1,q2,q3) 和 w(q1,···,q6)分別為中心線上一點P的矢徑和位移矢量.于是得到關(guān)于中心線的微分方程

其中γ=?sw為P點的應(yīng)變矢量,其主軸分量γi=γ·ei,(i=1,2,3)中的前兩個為截面的剪應(yīng)變,第三個為拉壓應(yīng)變,存在關(guān)系

Kirchhoff彈性桿是取γ=0[8?12].

截面的彎扭度ω和角速度?為截面姿態(tài)角對弧坐標(biāo)s和時間t全偏導(dǎo)數(shù)的線性組合

式中波浪號表示相對主軸坐標(biāo)系(P-xyz)求導(dǎo),ν=?tr,ω和? 的Euler角表達將使(6)式成為恒等式.(2),(3),(6)和(7)式組成彈性桿精確模型的運動學(xué)方程.

值得指出的是,應(yīng)變矢γ和彎扭度ω都是反映中心線彎曲的特征量.

3 彈性桿截面的虛位移

對于變形后的彈性桿,截面隨弧坐標(biāo)繞形心轉(zhuǎn)動的同時還存在平移.定義如下彈性桿截面的虛位移[8].

定義(截面的虛位移) 約束所允許的、與弧坐標(biāo)s和時間t變化無關(guān)的、假想的截面無限小位移定義為彈性桿截面的虛位移,它可分解為隨形心的虛平移和相對形心的虛角位移,分別記為δr和δΦ,它對應(yīng)于如下的變分定義

其中與虛位移δr和δΦ對應(yīng)的廣義坐標(biāo)變分依次記為δpqα,(α=1,2,3)和δrqα,(α=4,5,6).對僅為廣義坐標(biāo)qα,(α=1,2,3)的函數(shù)的變分,有δ=δp,同理,對僅為廣義坐標(biāo)qα,(α=4,5,6)的函數(shù)的變分,有δ=δr.

定義微分和變分δ,δp服從交換關(guān)系

和運算的普遍性.

由(2)和(5)式知,截面的虛位移存在如下關(guān)系:

并有

與Kirchhoff模型類似,可以證明存在以下運動學(xué)關(guān)系[13]

或

以及

或

其中波浪號表示相對主軸坐標(biāo)系求導(dǎo)或變分.這些關(guān)系也可以用Euler角驗證.

若彈性桿被限制在慣性空間(O-ξηζ)中的固定光滑曲面上,約束方程為

在約束的基本假設(shè)下[8,9,11,12],(14)式可以化作對截面位形的約束方程

式中 (ξ,η,?)為截面的形心坐標(biāo),b(s,μ(s,t))為截面與約束曲面接觸點相對截面形心的矢徑,bξ,bη,bζ是其投影.(15)式的變分可化為

(16)式可簡寫為

式中Aα為(16)式中坐標(biāo)變分前的系數(shù).方程(16)等價于理想約束條件

其中fC為曲面對桿的分布約束力集度.

4 彈性桿精確模型動力學(xué)的微分變分原理

考察原長為?s的微段桿(圖1),作用于其上的內(nèi)力主矢F,主矩M,慣性力以及約束力的“虛功率”與Kirchhoff模型的區(qū)別僅在于(2)式和?sˉs?=1.不計因變形產(chǎn)生的慣性力,并利用理想約束條件,可以得到如下的“虛功率”表達式

式中ρ為彈性桿沿中心線的線密度,J為截面的慣量并矢,在主軸坐標(biāo)系下的坐標(biāo)陣為J=其中Ji=ρIi/A,A為截面積,I1,I2為截面對主軸x,y的慣性矩,I3為截面對主軸z的極慣性矩,且有I3=I1+I2;f(s,t)和m(s,t)為沿中心線作用的連續(xù)分布力和分布力偶.

彈性桿精確模型動力學(xué)的 d’Alembert-Lagrange原理表述為:受有理想雙面約束的彈性桿在任意時刻的真實運動不同于運動學(xué)上的可能運動僅在于真實運動對于任意的虛位移,有

(19)和(20)式也可以直接從矢量形式的動力學(xué)方程導(dǎo)出.

設(shè)彈性桿服從線性本構(gòu)關(guān)系,用主軸分量表示為

下面將用矢量表達的“虛功率”(19)式化作分析力學(xué)形式.

根據(jù)(8)—(13)式和分析力學(xué)的經(jīng)典推導(dǎo),對(α=1,2,3)可以導(dǎo)出如下關(guān)系:

式中撇號和點號分別表示?s和?t,Sp=為截面的平移彈性應(yīng)變能和平移動能.

對(α=4,5,6)有以下關(guān)系:

以及[9]

或

可以導(dǎo)出

(25)和(26)式推導(dǎo)從略[9].(22)和(27)式的推導(dǎo)如下:

及

其中用到了關(guān)系

于是,(20)式化作Euler-Lagrange形式

式中,Λp=Sp?Tp,Λr=Sr?Tr,fα=f·?r/?qα,(α=1,2,3),fα=m·Ξα,(α=4,5,6).注意到有關(guān)系

(28)式可寫成

其中

為彈性桿動力學(xué)的Lagrange函數(shù).(29)式可化作Nielsen形式

特殊地,對于除端點外不受力作用的彈性桿平衡時,動能為零,即Tp=Tr=0,由矢量形式的平衡微分方程知截面主矢F為常矢量,不失一般性,設(shè)

此守恒量不能預(yù)先嵌入上面的彈性應(yīng)變能函數(shù)Sr和Sp,但可以證明存在標(biāo)量函數(shù)V,使得

其中

顯然,當(dāng)Ki→∞時,(34)式與Kirchhoff模型一致[1]于是,(29)式化為

其中U=Sr?V.(35)式是利用守恒量(32)式的彈性桿靜力學(xué)的d’Alembert-Lagrange原理.

5 彈性桿精確模型的動力學(xué)方程及其平衡問題的首次積分

對于除端點外不受約束的自由彈性桿,廣義坐標(biāo)變分均為獨立,由原理(20),(29)和(31)式分別導(dǎo)出矢量形式的動力學(xué)方程

和Lagrange方程

以及Nielsen方程

(37)和(38)式統(tǒng)一表達了彈性桿動力學(xué)的全部方程.對于受有形如(14)式的曲面約束,導(dǎo)出帶乘子的Lagrange方程

從原理(35)式導(dǎo)出平衡微分方程

討論(40)式的首次積分:

1)若取q1=ξ,q2=η,q3=ζ,則因?Sp/?qα=0(α=1,2,3),從式(40a)導(dǎo)出循環(huán)積分

此積分的力學(xué)意義是主矢在慣性坐標(biāo)軸的分量為常量,這和平衡的矢量方程結(jié)果一致;

2)若取Euler角為姿態(tài)坐標(biāo),與Kirchhoff模型相同的是,仍有?U/?ψ=0,即進動角仍是循環(huán)坐標(biāo),從式(40b)導(dǎo)出循環(huán)積分

3)因Sp,U都不顯含s,存在廣義能量積分

6 彈性桿精確模型的Hamilton原理和Hamilton正則方程

將(20)式乘d s·d t后再對s和t積分,化作積分變分原理

其中用到了d-δ交換關(guān)系(9)式和端點變分條件

當(dāng)彈性桿服從線性本構(gòu)關(guān)系(21),且不計主動分布力f(s,t)和m(s,t)時,(44)式可進一步化作

其中Λ已由(30)式定義.(46)式就是彈性桿精確模型動力學(xué)的Hamilton原理.直接計算變分,從原理(46)式可以導(dǎo)出方程(37)和(38).

定義正則變量qα,以及

從而解出

定義彈性桿的Hamilton函數(shù)

直接計算偏導(dǎo)數(shù),并注意到(37)式,導(dǎo)出彈性桿的Hamilton正則方程

3×6 個變量 qα,psα,ptα,(α=1,···,6),(48)式共有3×6方程,方程組封閉.正則方程(48)顯示了彈性桿動力學(xué)的特殊性,其相空間是3×6維.

7 結(jié)論

基于位形空間上的虛位移定義建立的精確Cosserat彈性桿動力學(xué)的d’Alembert-Lagrange原理,可以導(dǎo)出Lagrange方程和Nielsen方程.建立的Hamilton原理和Hamilton正則方程保持了其經(jīng)典形式.從而將存在不足的Kirchhoff彈性桿動力學(xué)的分析力學(xué)方法推廣到計及拉壓和剪切變形情形,使彈性桿的全部動力學(xué)方程都表達成分析力學(xué)的形式,完成了分析力學(xué)方法對彈性桿動力學(xué)的移植.雙自變量特征也為分析力學(xué)的進一步研究提出了新的問題.

[1]Liu Y Z 2006 Nonlinear Mechanics of Thin Elastic Rod—Theoritical Basis of Mechanical Model of DNA(Beijing:Tsinghua University Press,Springer)p1432(in Chinese)[劉延柱2006彈性細桿的非線性力學(xué)—–DNA力學(xué)模型的理論基礎(chǔ)(北京:清華大學(xué)出版社Springer)第1432頁]

[2]Liu Y Z 2003 Mech.Eng.25 1(in Chinese)[劉延柱2003力學(xué)與實踐25 1]

[3]Ouyang Z C 2003 Physics 32 728(in Chinese)[歐陽鐘燦2003物理32 728]

[4]Li M,Ouyang Z C 2003 Science 55 15(in Chinese)[黎明,歐陽鐘燦2003科學(xué)55 15]

[5]Malacinski G M(translated by Wei Q)2005 Essentials of Molecular Biology(4th Ed.)(Bejing:Chemical Industry Press)p59(in Chinese)[喬治M馬拉森斯基著(魏群譯)2005分子生物學(xué)精要(北京:化學(xué)工業(yè)出版社)第59頁]

[6]Westcott T P,Tobias I,Olson WK 1995 J.Phys.Chem.99 17926

[7]Pozo Coronado L M 2000 Physica D 141 248

[8]Xue Y,Liu Y Z,Chen L Q 2005 Chin.J.Theor.Appl.Mech.37 485(in Chinese)[薛紜,劉延柱,陳立群2005力學(xué)學(xué)報37 485]

[9]Xue Y,Liu Y Z 2006 Acta Phys.Sin.55 3845(in Chinese)[薛紜,劉延柱2006物理學(xué)報55 3845]

[10]Xue Y,Liu Y Z 2006 Chinese Quarterly of Mechanics 27 550(in Chinese)[薛紜,劉延柱2006力學(xué)季刊27 550]

[11]Xue Y,Wen D W 2009 Acta Phys.Sin.58 34(in Chinese)[薛紜,翁德瑋2009物理學(xué)報58 34]

[12]Xue Y,Shang H L 2009 Chin.Phys.Lett.26 074501

[13]Zhao WJ,Weng Y Q,Fu JL 2007 Chin.Phys.Lett.24 2773

[14]Wang P,Xue Y,Liu Y L 2012 Chin.Phys.B 21 070203

[15]Cao D Q,Tucker R W 2008 Int.J.Solids Struct.45 460

[16]Xue Y,Weng D W,Chen L Q 2009 Chinese Quarterly of Mechanics 30 116(in Chinese)[薛紜,翁德瑋,陳立群2009力學(xué)季刊30 116]

[17]Xue Y,Weng D W 2011 Mech.Eng.33 65(in Chinese)[薛紜,翁德瑋2011力學(xué)與實踐33 65]