Kirchhoff動(dòng)力學(xué)比擬對彈性薄殼的推廣1)

薛紜 陳立群

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

?(哈爾濱工業(yè)大學(xué)(深圳)理學(xué)院力學(xué)系,深圳 518055)

引言

彈性細(xì)桿的平衡與剛體的定點(diǎn)轉(zhuǎn)動(dòng)是兩類完全不同的經(jīng)典力學(xué)問題.1859 年Kirchhoff等[1-2]在平面截面假定下建立的彈性細(xì)桿靜力學(xué)理論,依據(jù)其平衡的Kirchhoff方程與剛體動(dòng)力學(xué)的Euler 方程在數(shù)學(xué)形式上的相似性,提出了Kirchhoff動(dòng)力學(xué)比擬方法.1927 年,Love[3]在其彈性力學(xué)著作中詳細(xì)敘述Kirchhoff理論,明確了兩者在物理概念上的對應(yīng)和在數(shù)學(xué)表達(dá)上的等同關(guān)系.對此,劉延柱[4]作了更明確的敘述,如中心線弧長坐標(biāo)?時(shí)間;截面的姿態(tài)?定點(diǎn)轉(zhuǎn)動(dòng)剛體的姿態(tài);截面的彎扭度?剛體的角速度;截面的抗彎抗扭剛度?剛體對主軸的轉(zhuǎn)動(dòng)慣量等.

促使彈性細(xì)桿力學(xué)近代發(fā)展的是分子生物學(xué).作為脫氧核糖核酸(DNA) 雙螺旋結(jié)構(gòu)的力學(xué)模型,彈性細(xì)桿的Kirchhoff動(dòng)力學(xué)比擬方法成為描述DNA雙螺旋復(fù)雜位形的有力工具[5].Coleman 和Swigon[6]于2004 年由Kirchhoff方程得到有結(jié)和無結(jié)質(zhì)粒超螺旋的精確解.從力學(xué)的角度看,動(dòng)力學(xué)比擬也為連續(xù)彈性細(xì)桿提供新的離散化方法,從而成為3 自由度系統(tǒng),彈性細(xì)桿的位形是截面沿中心線的弧坐標(biāo)歷程.位形、曲率和撓率等幾何概念都可用彎扭度和Euler 角表達(dá)[4],開辟了用剛體動(dòng)力學(xué)的理論和方法研究彈性細(xì)桿靜力學(xué)的新思路,成為求解彈性細(xì)桿小應(yīng)變大位移問題的有效方法.同時(shí)還極大地豐富了彈性細(xì)桿力學(xué)理論.分析力學(xué)方法[7]、非完整約束[8-9]、對稱性和守恒量[10]、Lyapunov 穩(wěn)定性[11],甚至混沌[12]等動(dòng)力學(xué)概念和理論都可移植或應(yīng)用到彈性細(xì)桿靜力學(xué),并賦予新的含義.同時(shí)Kirchhoff動(dòng)力學(xué)比擬也為剛體動(dòng)力學(xué)的運(yùn)動(dòng)歷程提供幾何直觀圖像.

彈性細(xì)桿中心線的二維擴(kuò)展形成彈性薄殼的中面.自然希望將彈性細(xì)桿的Kirchhoff動(dòng)力學(xué)比擬方法推廣到彈性薄殼.

隨著工程技術(shù)的發(fā)展,大型柔性航天器動(dòng)力學(xué)和多柔體系統(tǒng)動(dòng)力學(xué)應(yīng)用日益廣泛[13].文獻(xiàn)[14]簡要回顧了多柔體系統(tǒng)動(dòng)力學(xué)建模的3 類方法:浮動(dòng)坐標(biāo)方法、幾何精確方法和絕對節(jié)點(diǎn)坐標(biāo)方法.對于大變形柔性梁,文獻(xiàn)[15]在浮動(dòng)坐標(biāo)系下使用柔性梁參數(shù)方程形式的精確曲率公式來計(jì)算柔性梁的彎曲變形能,建立了基于浮動(dòng)坐標(biāo)系的考慮曲率縱向變形效應(yīng)的剛耦合動(dòng)力學(xué)模型.

對于幾何非線性問題的建模與數(shù)值計(jì)算的主要方法有初應(yīng)力法[16]、幾何大變形約束法[17]、大變形耦合方法等[18].在大變形情況下,文獻(xiàn)[19]提出了用絕對節(jié)點(diǎn)坐標(biāo)法建模,解決了用混合坐標(biāo)法建模計(jì)算效率較低的問題,文獻(xiàn)[20]進(jìn)一步用大變形實(shí)驗(yàn)驗(yàn)證了絕對節(jié)點(diǎn)坐標(biāo)法的正確性.文獻(xiàn)[21]不使用微分幾何中的矩陣及其退化概念,在彈性中面的原始位形和當(dāng)前位形之外,定義了平面參考位形,在這3 個(gè)位形上定義正交坐標(biāo)系,通過映射表達(dá)他們間的關(guān)系.用Cartesian 坐標(biāo)表達(dá)應(yīng)變.

文獻(xiàn)[22]指出,“顯然,研究薄殼大變形大位移比研究剛體的一般空間運(yùn)動(dòng)更復(fù)雜,后者的自由度為6,而對可變形物體卻有無窮多個(gè)自由度.”將Kirchhoff動(dòng)力學(xué)比擬方法應(yīng)用到彈性薄殼,有望將彈性薄殼化作離散力學(xué)系統(tǒng),從而經(jīng)典動(dòng)力學(xué)理論和方法可以用到彈性薄殼力學(xué),為彈性薄殼靜力學(xué)注入新的概念和方法.

經(jīng)典的彈性細(xì)桿和彈性薄殼的現(xiàn)代研究在對象,以及研究內(nèi)容和研究方法上都賦予了全新的內(nèi)涵[23].將彈性細(xì)桿Kirchhoff動(dòng)力學(xué)比擬的思想,推廣到彈性薄殼是很有意義的工作.在彈性中面上建立正交軸系并定義姿態(tài)坐標(biāo)和彎扭度.用Lamé系數(shù)、彎扭度和姿態(tài)坐標(biāo)表達(dá)中面偏微分方程、一點(diǎn)的應(yīng)變和應(yīng)力狀態(tài)、截面的內(nèi)力,以及平衡方程等.此外,彎扭度和Lamé系數(shù)還可表達(dá)曲面微分幾何的基本概念和公式,這將另文討論.

本文在直法線假設(shè)下[24-26],基于Kirchhoff動(dòng)力學(xué)比擬方法和剛體動(dòng)力學(xué)的觀念,用Euler 角、彎扭度和Lamé系數(shù)討論彈性薄殼的靜力學(xué)建模.

文中指標(biāo)取值i,j=1,2;k=1,2,3.

1 彈性薄殼中面的軸系

設(shè)O?ξηζ 為慣性參照系,沿坐標(biāo)軸方向的單位基矢量為eξ,eη,eζ.變形前彈性中面上的任意點(diǎn)p變形后為P,中面方程分別為

其中q1,q2為廣義弧坐標(biāo).設(shè)r和R具有所需要的各階連續(xù)偏導(dǎo)數(shù)

為位移函數(shù),如圖1 所示.

圖1 變形前和變形后中面上點(diǎn)的矢徑和位移矢量Fig.1 The position vector and displacement vector of the point on the mid-surface before and after deformation

變形后中面上沿坐標(biāo)線的切矢量、單位切矢量和弧長微分依次為

其中Lamé系數(shù)Ti=|Ti|.在中面上一點(diǎn)P建立兩個(gè)軸系其中

式中,N0為變形后中面的單位法矢量,sin φ=物理上要求如圖2 所示.文獻(xiàn)上稱為點(diǎn)P的軸系,或Cosserat 方向子或活動(dòng)標(biāo)架.可以用線性表為

圖2 彈性中面上兩軸系的相對位置Fig.2 The position relation between two frames on the elastic mid-surface

對中面的任意曲線qi=qi(q)的單位切矢量

和其q坐標(biāo)線的弧長微分dS,將式(5)第1 式代入式(6),化作

得到其弧長微分關(guān)系為

其中dS3是方向的弧長微分.

同理,上述概念可以在變形前的彈性中面上定義,并用下標(biāo)0 表示

式中dsi為變形前中面坐標(biāo)線的弧長微分.

2 彈性薄殼的彎扭度及其相容方程

2.1 彎扭度的定義

以變形前或變形后中面坐標(biāo)線弧長微分為基準(zhǔn),即定義式中偏導(dǎo)分別對變形前和變形后的坐標(biāo)線弧長進(jìn)行,彎扭度分別記為三者的關(guān)系為

彎扭度的分量形式記為

若矢量V寫作V=VV0,其中V0為單位矢量.其偏導(dǎo)數(shù)用彎扭度表為

以變形前廣義弧坐標(biāo)微分dq為基準(zhǔn),沿任意方向的彎扭度定義為

其中沿中面上曲線qi=qi(q)的單位切矢量為用沿坐標(biāo)線的基于變形后弧長的彎扭度表為

利用彎扭度的定義式(15),平面的彎扭度特征是

在曲面微分幾何中,Rodrigues 方程指出了曲面主方向和主曲率的幾何特征[28].當(dāng)曲面的坐標(biāo)線為曲率線時(shí),有

式中κni為兩個(gè)主曲率.這也是Rodrigues 方程的彎扭度表達(dá).它指出了主方向的彎扭度特征,進(jìn)一步明確了彎扭度的幾何意義.

2.2 等距面上的Lamé系數(shù)和坐標(biāo)線夾角

板殼中與中面等距離的面稱為等距面.在直法線假設(shè)下,變形前和變形后等距面方程分別為其中h為彈性薄殼的厚度,z為厚度坐標(biāo)(?h/2 ≤z≤h/2),N0為N0的原始值.

計(jì)算變形后等距面坐標(biāo)線的夾角.等距面上沿坐標(biāo)線的切矢量為

略去高階微量,Lamé系數(shù)和弧長微分為

其中用到了對z的一階近似.單位切向量為

等距面上坐標(biāo)線夾角φz由下式給出

在略去高階微量后化作

同理,變形前沿坐標(biāo)線的切矢量和單位切矢量有類似形式,這里不再贅述.

2.3 彎扭度的相容性

定義的彎扭度ωij不是互相獨(dú)立的,存在相容協(xié)調(diào)關(guān)系.設(shè)彈性中面不可撕裂或折疊,這要求彎扭度是關(guān)于廣義弧坐標(biāo)的連續(xù)可微函數(shù),并滿足相容協(xié)調(diào)關(guān)系.

2.3.1 彎扭度ω1j,ω2j 的相容性

其在軸系(e21,e22,e23)上投影的矩陣式為

證明:對軸系(e11,e12,e13)和矩陣Φ,容易驗(yàn)證如下關(guān)系

將式(10)代入式(15),即可導(dǎo)出式(32).再將ωij向軸系(ei1,ei2,ei3)投影,并注意到式(10),即得到式(33).

值得說明的是,式(32) 對任意方向的彎扭度也成立,即

其中ωiα由式(20)定義.

2.3.2 彎扭度ωi1,ωi2 的相容性

彈性中面上同一個(gè)軸系沿不同坐標(biāo)線方向彎扭度ωi1,ωi2的相容協(xié)調(diào)關(guān)系有下面四種等價(jià)形式

其中波浪號表示相對軸系(ei1,ei2,ei3)的偏導(dǎo)數(shù).

證明:設(shè)a為給定軸系(ei1,ei2,ei3)中的固定常矢量,由式(19)依次對q1,q2求混合偏導(dǎo)數(shù).由求偏導(dǎo)次序的可交換性,以及三重矢積恒等式,得到

因a為任意的固定矢量,于是就導(dǎo)出式(36).再由式(19),從式(36)導(dǎo)出式(37).同理可導(dǎo)出式(38)和式(39).

2.3.3 彎扭度相容方程的獨(dú)立性

可以證明式(36)給出的兩個(gè)式子可以互相推導(dǎo),即

將式(32)代入式(41)等號左邊即可推出等號右邊.

2.3.4 中面上彎扭度與切矢量的相容性

彈性中面上彎扭度與一點(diǎn)的切矢量也必須是相容的.切矢量滿足如下關(guān)系

這既是曲線積分

與路徑無關(guān)的條件,保證變形后的彈性中面無裂縫;也是兩沿坐標(biāo)線切矢量的相容協(xié)調(diào)關(guān)系.

注意到Ti=式(42)寫作

其分量形式為

這是關(guān)于彎扭度和Lamé系數(shù)的約束方程.由式(33),式(45)右端的彎扭度可用ω1i表達(dá),化作

式(32) 和式(36) 及其等階形式,以及約束方程(46),在彎扭度用Euler 角及其偏導(dǎo)數(shù)表示時(shí)再作進(jìn)一步討論.

3 軸系姿態(tài)和彎扭度中面偏微分方程

3.1 軸系姿態(tài)的Euler 角

建立慣性參照系O?ξηζ.借助Euler 角將軸系(ei1,ei2,ei3) 用慣性參照系的基矢量(eξ,eη,eζ)表達(dá).為此將O?ξηζ 平移到彈性中面的P點(diǎn),繞軸的如下三次轉(zhuǎn)動(dòng)到達(dá)軸系(ei1,ei2,ei3)的位置P?xyz,如圖3 所示[4,29]

圖3 軸系(ei1, ei2, ei3)姿態(tài)的Euler 角Fig.3 Euler angles of attitude of the frame(ei1, ei2, ei3)

基的變換關(guān)系用矩陣表為

其中Qi為方向余弦矩陣

由圖2 知,ψ2=ψ1,θ2=θ1,φ2=φ1+φ,可略去ψ,θ 的下標(biāo);ψ=ψ(q1,q2),θ=θ(q1,q2),φi=φi(q1,q2)為軸系(ei1,ei2,ei3) 相對于慣性系的Euler 角;字母s,c 分別表示sin,cos.

兩方向余弦矩陣Qi之間有如下關(guān)系

式中Φ由式(10)定義.

當(dāng)給定矩陣Qi,可以通過式(48)從對應(yīng)元素相等解出相應(yīng)的Euler 角[4,29].表達(dá)軸系(ei1,ei2,ei3)在慣性參照系O?ξηζ 中的姿態(tài)還可以用Cardano 角或Euler 參數(shù)等,這里不再贅述.

3.2 彎扭度的Euler 角表示

由剛體動(dòng)力學(xué)知,用Euler 角表達(dá)的彎扭度ωij在基(ei1,ei2,ei3)下的矩陣式為[4,29]

將此彎扭度的Euler 角表達(dá)式代入彎扭度的相容方程(32) 和方程(36)～(39),都成為恒等式.表明彎扭度的相容方程對Euler 角不構(gòu)成約束.彎扭度的Euler 角表達(dá)式將使式(45)和式(46)轉(zhuǎn)化為對Euler角和Lamé系數(shù)及其偏導(dǎo)數(shù)的約束方程.

3.3 中面偏微分方程

在慣性參照系O?ξηζ 中,變形后彈性中面上一點(diǎn)P的直角坐標(biāo)為

式(3)給出中面的偏微分方程的矢量形式

當(dāng)已知Euler 角和Lamé系數(shù)ψ,θ,φ,φ;T1,T2隨廣義弧坐標(biāo)qi的變化規(guī)律,且滿足約束方程(46),式中彎扭度由式(50)給出,即

并且已知邊界條件,如

就可以求得彈性中面的位形.

特殊地,當(dāng)僅有i=1 時(shí),式(53)即為彈性細(xì)桿的軸線微分方程,此時(shí)R為彈性細(xì)桿中心線的矢徑,為其單位切矢量.可見這是彈性細(xì)桿軸線微分方程向彈性薄殼中面的推廣.

4 應(yīng)變和應(yīng)力及其彎扭度表達(dá)

一點(diǎn)沿任意q坐標(biāo)線方向的線應(yīng)變定義為[30]

式中dS是微分弧段ds在變形后的長度,T=dS/dq,T0=ds/dq為相應(yīng)的Lamé系數(shù).可近似表為[30]

4.1 等距面上一點(diǎn)的線應(yīng)變合成定理

等距面上一點(diǎn)的線應(yīng)變合成定理:彈性薄殼變形后等距面上一點(diǎn)Pz沿坐標(biāo)線qi的線應(yīng)變等于中面上點(diǎn)P沿該坐標(biāo)線的線應(yīng)變?chǔ)舉i和點(diǎn)Pz相對于點(diǎn)P的線應(yīng)變?chǔ)舝i,以及原始相對應(yīng)變負(fù)值?εri0的代數(shù)和,表示為

其中ur由式(25) 定義,為其原始值由式(24)定義.

若變形前彈性中面為曲面,?N0/?qi≠0.將式(24)代入式(57),此時(shí)式中ε 為略去二階以上小量便得到式(58).

下面用彎扭度表達(dá)εri,εri0.由式(25),得

代入式(61),取j=i,相對線應(yīng)變?chǔ)舝i用彎扭度表達(dá)為

即

在應(yīng)變的一次近似下,式中可取cos φ=0,sin φ=1,Ti/Ti0=1,寫作

同理,可得

4.2 等距面上一點(diǎn)的剪應(yīng)變

等距面上一點(diǎn)的剪應(yīng)變合成定理:變形后彈性等距面上一點(diǎn)Pz的剪應(yīng)變在一次近似下等于彈性中面上一點(diǎn)P的剪應(yīng)變?chǔ)胑xy和點(diǎn)Pz相對于點(diǎn)P的剪應(yīng)變?chǔ)胷xy,以及變形前點(diǎn)pz的“剪應(yīng)變”負(fù)值的代數(shù)和.表示為

推導(dǎo)如下:變形前彈性薄殼等距面坐標(biāo)線一般不是正交的.因此,按照剪應(yīng)變的定義,在小應(yīng)變下有

這就得到了式(68).

用彎扭度表達(dá)相對剪應(yīng)變.由式(70)得

在應(yīng)變的一次近似下,取T20/T2=T10/T1=1,式(77)和式(79)化作

4.3 等距面上一點(diǎn)沿任意方向的線應(yīng)變

這結(jié)果是顯然的.

可以驗(yàn)證,薄板的經(jīng)典大撓度理論是上述結(jié)果的特例.

4.4 應(yīng)力及其彎扭度表達(dá)

圖4 直法線假設(shè)下計(jì)及的應(yīng)變所對應(yīng)的應(yīng)力Fig.4 The stress corresponding to the strain considered under the straight normal assumption

設(shè)材料是線彈性的.由廣義Hooke 定律知

式中,E,G,μ分別為材料的拉壓彈性模量、剪切彈性模量和Poisson 比,三者滿足關(guān)系G=E/2(1+μ).

注意到式(84),將應(yīng)變式(58) 和式(68) 代入式(86),式中應(yīng)變分量由彎扭度和Lamé系數(shù)表達(dá),得

式中,σex,σey,τexy為薄膜應(yīng)力,σrx,σry為彎曲應(yīng)力,τrxy為扭轉(zhuǎn)剪應(yīng)力,而σex0,σey0,σrx0,σry0,τrxy0不具有應(yīng)力的意義.

5 彈性薄殼截面上的分布內(nèi)力和本構(gòu)方程

5.1 彈性薄殼截面上內(nèi)力的簡化

微分單元截面上應(yīng)力向中線簡化的結(jié)果為一分布主矩和分布主矢,如圖6 所示.

圖5 微分單元及其計(jì)及應(yīng)變的應(yīng)力狀態(tài)Fig.5 Differential unit and its stress state considering strain

圖6 微分單元及其內(nèi)力:分布主矩和分布主矢Fig.6 Differential unit and its internal force:Distributed principalmoment and distributed principal vector

其中,M1x,M3y為分布扭矩,M1y,M3x為分布彎矩,F1x,F3y為法向分布力,F1y,F1z,F3x,F3z為切向分布力.具體表達(dá)如下

以下直接寫為

可以驗(yàn)證,內(nèi)力主矩和主矢表達(dá)式與經(jīng)典的表達(dá)是一致的[24-26].

5.2 本構(gòu)方程

由式(100)～式(105),得到本構(gòu)方程

其中D為抗彎剛度.

根據(jù)直法線假定,F1z,F3z是獨(dú)立變量.

本構(gòu)方程也可以在應(yīng)變的一次近似下表達(dá),這里從略.

6 平衡方程和動(dòng)力學(xué)比擬

6.1 平衡方程

設(shè)沿彈性薄殼中面作用有隨體或固定的分布載荷

列出變形后微分單元的平衡方程并整理得

將式(3)、式(9)和式(97)代入,化作

這形式上與Kirchhoff方程和Euler 方程相似,但多了一個(gè)自變量.與經(jīng)典薄殼平衡微分方程也是一致的[24-26].

值得指出式(117) 的N0軸投影式不是獨(dú)立的.不失一般性,當(dāng)坐標(biāo)線qi是曲率線時(shí),φ=π/2,且有式(23)知,=0,可以證明式(117)沿N0的投影式是恒等式.需要注意的是,變形后的中面位形是待求的.

靜力方程式(116)和式(117)是2 自變量的偏微分方程組,與約束方程(46) 聯(lián)立,方程數(shù)位8,所含的未知量是ψ,θ,φ1,φ,F1z,F3z,T1,T2,個(gè)數(shù)相等.因此,方程組封閉.值得一提的是,這里的約束方程是運(yùn)算和變形規(guī)則要求,因而無需對應(yīng)的約束力,文獻(xiàn)上稱此類約束為內(nèi)約束.

由復(fù)合函數(shù)的求導(dǎo)規(guī)則和式(9),有

據(jù)此,可將對q3的偏導(dǎo)數(shù)化作對qi的偏導(dǎo)數(shù).

6.2 彈性薄殼的Kirchhoff動(dòng)力學(xué)比擬

式(115)與定點(diǎn)轉(zhuǎn)動(dòng)剛體的Euler 方程和彈性細(xì)桿Kirchhoff方程有一定相似性.按照Kirchhoff動(dòng)力學(xué)比擬,彈性薄殼靜力學(xué)與剛體動(dòng)力學(xué)有很好的對應(yīng)關(guān)系,主要有:時(shí)間坐標(biāo)t?中面廣義弧坐標(biāo)qi,Euler 角ψ,θ,φ–Euler 角ψ,θ,φi,角速度??彎扭度ωij;剛體的運(yùn)動(dòng)過程?彈性薄殼的位形,動(dòng)量矩?內(nèi)力主矩等.彈性薄殼可以取剛體的姿態(tài)坐標(biāo)和Lamé系數(shù)以及厚度方向的切向分布力為未知函數(shù),按照求解剛體運(yùn)動(dòng)的思路求解彈性薄殼的位形.

7 算例

設(shè)薄板寬度為b,長度為l,厚度為h,不計(jì)自重.在板兩端作用分布力偶

如圖7 所示.變形前彈性中面的矢徑為

Lamé系數(shù)Ti0=1,軸系基矢量為

圖7 變形前后的彈性薄板Fig.7 Thin elastic plate before and after deformation

N0方向的投影式為恒等式.再將式(124)和式(125)代入式(119),解得

考慮到約束方程式(46),解為

積分得

此時(shí)的彎扭度和中面的應(yīng)力為

由式(107)得到

注意到φ=π/2 和式(129),式(109)化作

解得

積分得

其中C為積分常數(shù).設(shè)=0,導(dǎo)出C=?6ml/(Eh3).

由式(53)得到變形后的中面偏微分方程

其中,K=12m/(Eh3).這就是薄板在兩邊分布力偶作用下的中面方程.是一圓弧面,圓心軸為平行于ξ軸的直線段(q1,l/2,1/K),半徑為1/K.這結(jié)論是熟知的,如圖8 所示.

圖8 彈性中面的位形Fig.8 Configuration of elastic mid-splane

8 結(jié)論

用中面軸系的姿態(tài)坐標(biāo)和彎扭度,以及Lamé系數(shù)表達(dá)了彈性薄殼中面的位形、一點(diǎn)的應(yīng)變狀態(tài)和應(yīng)力狀態(tài)、以及截面上的分布主矩和分布主矢,得到的平衡偏微分方程與定點(diǎn)轉(zhuǎn)動(dòng)剛體的Euler 方程和彈性細(xì)桿Kirchhoff方程有一定的相似性.使中面位形成為廣義弧坐標(biāo)歷程.將彈性桿的Kirchhoff動(dòng)力學(xué)比擬方法推廣到彈性薄殼.

在Kirchhoff動(dòng)力學(xué)比擬下,彈性薄殼的位形空間是是一個(gè)6 維空間,因?yàn)殡p自變量,且存在3 個(gè)非完整約束,因此自由度為9.可以引入完整和非完整約束的概念,并形成彈性薄殼的分析力學(xué)方法.