Hamel 框架下幾何精確梁的離散動(dòng)量守恒律1)

高 山 史東華,2) 郭永新

*(北京理工大學(xué)數(shù)學(xué)與統(tǒng)計(jì)學(xué)院,北京 100081)

?(遼寧大學(xué)物理學(xué)院,沈陽(yáng) 110036)

引言

柔性梁的動(dòng)力學(xué)廣泛存在于力學(xué)系統(tǒng)當(dāng)中,其研究發(fā)展不斷與固體力學(xué)、流體力學(xué)、非線性動(dòng)力學(xué)、生物力學(xué)等其他學(xué)科進(jìn)行交叉[1-4].幾何精確梁模型作為容許大范圍運(yùn)動(dòng)、大變形的柔性梁的理想模型,自提出后在分析、數(shù)值、應(yīng)用等方向吸引了大量關(guān)注[5-11].在柔性機(jī)器人的運(yùn)動(dòng)控制[12-13]、DNA動(dòng)力學(xué)的研究[14]中有廣泛的應(yīng)用.

在幾何精確梁的數(shù)值算法方面,傳統(tǒng)數(shù)值格式直接離散動(dòng)力學(xué)方程,可以實(shí)現(xiàn)連續(xù)系統(tǒng)的計(jì)算仿真,但是很難保持系統(tǒng)內(nèi)在的對(duì)稱性、守恒性等物理特性,如能量、動(dòng)量等.Marsden 等[15]提出的變分積分子是一種從Lagrange 函數(shù)出發(fā),通過離散變分原理直接得到保結(jié)構(gòu)數(shù)值格式的方法.Demoures 等[16]提出的幾何精確梁的李群積分子和李代數(shù)積分子可以自然保持能量動(dòng)量,但沒有充分利用系統(tǒng)對(duì)稱性,對(duì)于幾何精確梁等對(duì)稱性很強(qiáng)的系統(tǒng)離散格式復(fù)雜程度很高.在此基礎(chǔ)上,Ball 等[17]對(duì)有限維系統(tǒng)提出了Hamel 變分積分子,主要適用于帶對(duì)稱性的約束力學(xué)系統(tǒng),特別是帶對(duì)稱性的非完整約束力學(xué)系統(tǒng).Shi 等[18-20]進(jìn)一步應(yīng)用無(wú)窮維幾何中的活動(dòng)標(biāo)架法結(jié)合變分原理,得到約束無(wú)窮維力學(xué)系統(tǒng)和經(jīng)典場(chǎng)論中刻畫運(yùn)動(dòng)的標(biāo)準(zhǔn)型—–Hamel 形式,進(jìn)而得到能用于經(jīng)典場(chǎng)論動(dòng)力學(xué)計(jì)算的保結(jié)構(gòu)格式—–Hamel 場(chǎng)變分積分子.Hamel 場(chǎng)變分積分子可以無(wú)需額外的Lagrange 乘子來(lái)求解無(wú)窮維約束力學(xué)系統(tǒng),特別是非完整力學(xué)系統(tǒng).其是一種分布式算法,每一次計(jì)算只需利用待計(jì)算點(diǎn)周圍的信息,具有計(jì)算效率高的特點(diǎn).王亮等[7]將Hamel 場(chǎng)變分積分子應(yīng)用于幾何精確梁,并通過數(shù)值算法驗(yàn)證了Hamel 場(chǎng)變分積分子能長(zhǎng)時(shí)間保結(jié)構(gòu)的特征.

保結(jié)構(gòu)性質(zhì)對(duì)非線性力學(xué)系統(tǒng)的長(zhǎng)時(shí)間定性行為的了解和定量分析是非常必要的.這里結(jié)構(gòu)包括辛結(jié)構(gòu)、酉結(jié)構(gòu)、體積結(jié)構(gòu)、接觸結(jié)構(gòu)、泊松結(jié)構(gòu)等幾何結(jié)構(gòu),也包括能量、動(dòng)量、系統(tǒng)特征值等首次積分結(jié)構(gòu)[21].在天體力學(xué)[22-23]、量子力學(xué)[24]、電磁學(xué)[25]等學(xué)科中許多模型數(shù)值算法的構(gòu)造中,盡可能多地保持原系統(tǒng)的內(nèi)在對(duì)稱性、守恒性等物理特征的算法往往表現(xiàn)出很強(qiáng)的數(shù)值預(yù)測(cè)能力和跟蹤能力[21].例如,電磁學(xué)中受到廣泛應(yīng)用的Yee 格式已被證明保持多辛結(jié)構(gòu)和動(dòng)量映射[26].在保結(jié)構(gòu)算法的理論研究上,鐘萬(wàn)勰等[27]建立了一套Hamilton 動(dòng)力系統(tǒng)的辛幾何方法,張素英等[28]在此基礎(chǔ)上,對(duì)Poisson 流形上的廣義Hamilton 系統(tǒng)的數(shù)值算法進(jìn)行了研究,提出了廣義Hamilton 系統(tǒng)顯式的保持真解典則性的數(shù)值方法.劉世興等[29]研究了特定的非完整系統(tǒng)的辛算法,并將其與Runge-Kutta 法進(jìn)行對(duì)比,說(shuō)明了辛算法具有高的精確度.滿淑敏等[30]提出了一類求解非完整系統(tǒng)的保結(jié)構(gòu)算法,并證明了其對(duì)稱性,說(shuō)明了其收斂性和長(zhǎng)時(shí)間保結(jié)構(gòu)的特性.

對(duì)Hamel 變分積分子的現(xiàn)有研究結(jié)果在數(shù)值上充分表現(xiàn)出了其對(duì)非線性系統(tǒng),包括有限維系統(tǒng)、無(wú)窮維系統(tǒng)、完整系統(tǒng)與非完整系統(tǒng),長(zhǎng)時(shí)間動(dòng)力學(xué)行為的良好預(yù)測(cè)[7,17,20].本文以幾何精確梁為例,利用離散Noether 定理,研究了Hamel 場(chǎng)變分積分子的保動(dòng)量性質(zhì),以期為進(jìn)一步研究其保結(jié)構(gòu)性質(zhì)提供參考依據(jù).

1 幾何精確梁的動(dòng)力學(xué)模型

考慮一個(gè)在空間中運(yùn)動(dòng)的質(zhì)量均勻的彈性梁.設(shè){E1,E2,E3} 為物質(zhì)標(biāo)架的一組基,初始時(shí)刻梁沿E2軸水平自然放置(無(wú)彈性形變),一端位于原點(diǎn)(見圖1).設(shè)梁長(zhǎng)為l,密度為ρ,截面面積為A的正方形.

圖1 幾何精確梁初始位形Fig.1 Initial configuration of a geometrically exact beam

在幾何精確梁的模型中,梁的截面做剛體運(yùn)動(dòng),其位形可由

描述,其中φ:R×[0,l] → R3為中線位置函數(shù),Λ:R×[0,l] →S O(3) 為截面旋轉(zhuǎn)矩陣.引入對(duì)流速度

和對(duì)流應(yīng)變

其中頂標(biāo)·和上標(biāo)′分別表示對(duì)時(shí)間變量t和空間變量s求導(dǎo),e2為E2方向的單位向量.此處利用了李代數(shù)se(3) 與R6的等同.對(duì)流速度ζ包含了截面運(yùn)動(dòng)的角速度ω和平移速度π的信息,對(duì)流應(yīng)變包含了彈性拉伸、彎曲和剪切信息.從而,幾何精確梁的Lagrange 密度可寫作

其中〈·,·〉為se*(3)與se(3)配對(duì),K,P 分別為常對(duì)角陣.從而作用泛函為

令η=g-1δg為獨(dú)立變分,易知η與ζ,γ滿足如下變分公式[20]

其中[·,·]為se(3)上的李括號(hào).

根據(jù)變分公式和Hamilton 原理,計(jì)算得到自由幾何精確梁的動(dòng)力學(xué)方程

和自由邊界條件

其中[·,·]*為se(3)李括號(hào)[·,·]的對(duì)偶括號(hào).同時(shí),由式(1)和式(2)可知,ζ和γ滿足相容性條件

相容性條件描述了對(duì)流速度和對(duì)流應(yīng)變生成的分布的可積性.在場(chǎng)論框架下,位形g是以時(shí)空為底流形,纖維為的主叢S E(3)的截面,相容性條件在幾何上可解釋為主叢的局部平坦性條件.對(duì)一般纖維叢情況下的相容性條件的討論參見文獻(xiàn)[19].對(duì)于給定初始條件和邊界條件的幾何精確梁,其運(yùn)動(dòng)由動(dòng)力學(xué)方程(3)和相容性條件(5)共同確定,將式(3)和式(5) 構(gòu)成的方程組稱為幾何精確梁的Hamel 場(chǎng)方程組[18-19].

2 幾何精確梁Hamel 場(chǎng)變分積分子

通過離散Hamilton 原理得到Hamel 場(chǎng)變分積分子.設(shè)離散梁的空間節(jié)點(diǎn)數(shù)為K,空間步長(zhǎng)為Δh,時(shí)間步數(shù)為N,時(shí)間步長(zhǎng)為Δt.梁的位形由序列n=0,1,2,···,N-1,j=1,2,···,K描述,表示在n時(shí)刻第j個(gè)截面的位形.分別定義離散對(duì)流速度和離散對(duì)流應(yīng)變?nèi)缦?/p>

由于指數(shù)映射exp :se(3) →S E(3)是局部微分同胚,在exp:se(3)→S E(3)與相距較近,即在exp 的值域中時(shí),上式唯一確定了離散對(duì)流速度,離散對(duì)流應(yīng)變的情形類似.從而,離散Lagrange 密度為相應(yīng)的作用和為

注:上式中第一項(xiàng)的求和范圍應(yīng)為n=0,1,2,···,N-1,j=1,2,···,K,第二項(xiàng)的求和范圍應(yīng)為n=1,2,···,N-1,j=0,1,2,···,K.為了后續(xù)計(jì)算表達(dá)式的簡(jiǎn)潔,本文將上式求和范圍統(tǒng)一寫為n=0,1,2,···,N-1,j=0,1,2,···,K,將指標(biāo)溢出的變量統(tǒng)一記為0.

及離散相容性條件

3 Hamel 場(chǎng)變分積分子保動(dòng)量性質(zhì)

本節(jié)首先證明自由的(不受外力、外力矩) 的幾何精確梁系統(tǒng)具有動(dòng)量守恒性質(zhì),隨后通過離散Noether 定理[31]證明Hamel 場(chǎng)變分積分子可以保持精確離散動(dòng)量.

根據(jù)經(jīng)典力學(xué)理論,一個(gè)不受外力(矩) 的完整力學(xué)系統(tǒng)總動(dòng)量保持不變,故自由運(yùn)動(dòng)的幾何精確梁系統(tǒng)具有動(dòng)量守恒性質(zhì).動(dòng)量守恒往往是在空間坐標(biāo)下進(jìn)行驗(yàn)證的,由于Hamel 場(chǎng)變分積分子是在活動(dòng)標(biāo)架下給出的,為建立其離散動(dòng)量守恒,用活動(dòng)標(biāo)架表述的梁的動(dòng)量.

根據(jù)Noether 定理,令φ:[a,b]×[0,l]×(-?0,?0)→R,?0＞0 為任意一個(gè)以時(shí)空和變分變量為自變量的C1光滑的函數(shù),并且滿足φ(a,s,?)=φ(b,s,?)=φ(t,0,?)=φ(t,l,?)=φ(t,s,0)=0,?t∈[a,b],?s∈[0,l].取定李代數(shù)中任意一個(gè)元素ξ∈se(3),構(gòu)造變分曲線

其中g(shù)(t,s)為Hamel 場(chǎng)方程的任意一個(gè)精確解.下文中均用表示變量依賴于變分參數(shù)?.于是的偏導(dǎo)數(shù)和無(wú)窮小變分δg可分別寫作

從而,變分曲線對(duì)應(yīng)的對(duì)流速度、對(duì)流應(yīng)變和體坐標(biāo)系下的變分分別為

將變分曲線代入作用泛函并對(duì)其變分,得到

將式(3)代入,得到

其中Ad*表示Ad 算子的伴隨算子.由于δφ和ξ的任意性可知

此式為場(chǎng)論框架下的動(dòng)量律.進(jìn)一步,對(duì)式(11) 在[0,l]上積分,利用空間邊界條件(4)可以得到即系統(tǒng)總動(dòng)量守恒.將此表達(dá)式改寫到空間坐標(biāo)系中,即可得到角動(dòng)量守恒和平動(dòng)動(dòng)量守恒的經(jīng)典表達(dá)式.下面,在上述連續(xù)情形動(dòng)量守恒證明的基礎(chǔ)上,利用離散Noether 定理進(jìn)一步證明Hamel 變分積分子也可以自然保持系統(tǒng)離散動(dòng)量.定義變分序列

體坐標(biāo)下的變分為

根據(jù)式(13)和式(14)和BCH 公式[32]可知

其中算子I:se(3)→se(3)*為

同理可得

根據(jù)離散Noether 定理,將變分序列代入離散Lagrange 密度并對(duì)其作用和變分,得到

其中最后一個(gè)等式利用了變分的定義,可知小量o(?)的變分為零.將離散Hamel 方程代入上式可得

對(duì)任意給定的ξ,由于的任意性可知

上式為式(11)的離散對(duì)應(yīng)形式.

上式對(duì)j在1,2,···,K上求和并用離散邊界條件(9)可得

其中I*為I 的伴隨算子.

綜上,可以得到如下定理定理3.1.

幾何精確梁的Hamel 場(chǎng)變分積分子(8) 和(10)在邊界條件(9)下滿足如下定義的離散動(dòng)量守恒

Jd表達(dá)式中級(jí)數(shù)的存在技術(shù)上是李代數(shù)se(3)非交換性在BCH 公式中的反映,而根本上由于在Hamel 框架下,離散對(duì)流速度處于位形與之間,故位形和離散對(duì)流速度值無(wú)法完全匹配,需離散動(dòng)量表達(dá)式中出現(xiàn)級(jí)數(shù)來(lái)彌補(bǔ).這種不匹配對(duì)于Hamel 場(chǎng)變分積分子的構(gòu)造是必要的,可以使所得變分積分子具有更高精度[20].值得注意的是,上式中的離散動(dòng)量可以展開為

這里的首項(xiàng)可作為連續(xù)動(dòng)量(11) 的直接離散,從而Hamel 變分積分子在O(Δt) 誤差范圍內(nèi)可以保持此離散動(dòng)量的首項(xiàng)[7].

4 算例

下面通過兩個(gè)例子分別從分析和數(shù)值角度驗(yàn)證上述結(jié)論.

例一考慮一個(gè)沿E3方向以v0初速度做勻速直線運(yùn)動(dòng)的幾何精確梁,不妨設(shè)梁長(zhǎng)l=1,梁的質(zhì)量和轉(zhuǎn)動(dòng)慣量均為1,那么系統(tǒng)總動(dòng)量非零分量為Jlin3=v0.

在Hamel 框架下,梁的位形可寫作

其中y,z分別為截面中線在E2,E3方向的坐標(biāo).從而對(duì)流速度有表達(dá)式

在se(3) 與R6的等同下,ζ=(0,0,0,0,0,)T.令a=(0,0,0,0,0,1)T,由式(12)可得系統(tǒng)在E3方向上的動(dòng)量為

考慮系統(tǒng)的離散動(dòng)量,由于在此場(chǎng)景下非平凡的對(duì)稱群為沿E3軸做平移運(yùn)動(dòng)的加法群,指數(shù)映射exp退化為線性映射,于是可得

E3方向上離散動(dòng)量表達(dá)式為

這與空間坐標(biāo)系下的經(jīng)典動(dòng)量表達(dá)一致.

例二考慮初始位形如圖1 所示的不受外力的幾何精確梁.其參數(shù)[33]如下:梁長(zhǎng)l=1,橫截面是邊長(zhǎng)為0.1 的正方形,密度ρ=1000,楊氏模量E=107,泊松比ν=0.35.

取時(shí)間步長(zhǎng)Δt=10-4s,空間節(jié)點(diǎn)數(shù)K=10,給定梁的初始對(duì)流速度為

初始對(duì)流應(yīng)變?yōu)?/p>

通過Hamel 場(chǎng)變分積分子(8)～(10)迭代計(jì)算10 000步,得到幾何精確梁的角動(dòng)量和線動(dòng)量如圖2 和圖3 所示,說(shuō)明了Hamel 場(chǎng)變分積分子可以長(zhǎng)時(shí)間以O(shè)(Δt)精度保持系統(tǒng)角動(dòng)量和線動(dòng)量.

圖2 梁的角動(dòng)量Fig.2 Evolution of angular momentum

圖3 梁的線動(dòng)量Fig.3 Evolution of linear momentum

5 結(jié)論

對(duì)于不受外力、外力矩的力學(xué)系統(tǒng),動(dòng)量守恒是一個(gè)基本的守恒律.對(duì)幾何精確梁模型,傳統(tǒng)方法往往將其看作大量剛體的串聯(lián)或位形空間為無(wú)窮維流形的無(wú)窮維力學(xué)系統(tǒng)來(lái)研究其動(dòng)量,而本文在協(xié)變場(chǎng)論觀點(diǎn)下,從時(shí)空角度用活動(dòng)標(biāo)架法研究Hamel 場(chǎng)變分積分子的保動(dòng)量性質(zhì),進(jìn)一步說(shuō)明了Hamel 場(chǎng)變分積分子保結(jié)構(gòu)的性質(zhì),同時(shí)也間接說(shuō)明其具有長(zhǎng)時(shí)間數(shù)值穩(wěn)定性.本文所給的證明方法也同樣適用于一般經(jīng)典場(chǎng)論場(chǎng)景下的Hamel 場(chǎng)變分積分子.下一步將探究Hamel 形式及變分積分子的保能量、保辛結(jié)構(gòu)等保結(jié)構(gòu)性質(zhì),并將其應(yīng)用到殼、膜等更復(fù)雜力學(xué)系統(tǒng)的建模與仿真中.