高 山 史東華,2) 郭永新
*(北京理工大學數學與統計學院,北京 100081)
?(遼寧大學物理學院,沈陽 110036)
柔性梁的動力學廣泛存在于力學系統當中,其研究發展不斷與固體力學、流體力學、非線性動力學、生物力學等其他學科進行交叉[1-4].幾何精確梁模型作為容許大范圍運動、大變形的柔性梁的理想模型,自提出后在分析、數值、應用等方向吸引了大量關注[5-11].在柔性機器人的運動控制[12-13]、DNA動力學的研究[14]中有廣泛的應用.
在幾何精確梁的數值算法方面,傳統數值格式直接離散動力學方程,可以實現連續系統的計算仿真,但是很難保持系統內在的對稱性、守恒性等物理特性,如能量、動量等.Marsden 等[15]提出的變分積分子是一種從Lagrange 函數出發,通過離散變分原理直接得到保結構數值格式的方法.Demoures 等[16]提出的幾何精確梁的李群積分子和李代數積分子可以自然保持能量動量,但沒有充分利用系統對稱性,對于幾何精確梁等對稱性很強的系統離散格式復雜程度很高.在此基礎上,Ball 等[17]對有限維系統提出了Hamel 變分積分子,主要適用于帶對稱性的約束力學系統,特別是帶對稱性的非完整約束力學系統.Shi 等[18-20]進一步應用無窮維幾何中的活動標架法結合變分原理,得到約束無窮維力學系統和經典場論中刻畫運動的標準型—–Hamel 形式,進而得到能用于經典場論動力學計算的保結構格式—–Hamel 場變分積分子.Hamel 場變分積分子可以無需額外的Lagrange 乘子來求解無窮維約束力學系統,特別是非完整力學系統.其是一種分布式算法,每一次計算只需利用待計算點周圍的信息,具有計算效率高的特點.王亮等[7]將Hamel 場變分積分子應用于幾何精確梁,并通過數值算法驗證了Hamel 場變分積分子能長時間保結構的特征.
保結構性質對非線性力學系統的長時間定性行為的了解和定量分析是非常必要的.這里結構包括辛結構、酉結構、體積結構、接觸結構、泊松結構等幾何結構,也包括能量、動量、系統特征值等首次積分結構[21].在天體力學[22-23]、量子力學[24]、電磁學[25]等學科中許多模型數值算法的構造中,盡可能多地保持原系統的內在對稱性、守恒性等物理特征的算法往往表現出很強的數值預測能力和跟蹤能力[21].例如,電磁學中受到廣泛應用的Yee 格式已被證明保持多辛結構和動量映射[26].在保結構算法的理論研究上,鐘萬勰等[27]建立了一套Hamilton 動力系統的辛幾何方法,張素英等[28]在此基礎上,對Poisson 流形上的廣義Hamilton 系統的數值算法進行了研究,提出了廣義Hamilton 系統顯式的保持真解典則性的數值方法.劉世興等[29]研究了特定的非完整系統的辛算法,并將其與Runge-Kutta 法進行對比,說明了辛算法具有高的精確度.滿淑敏等[30]提出了一類求解非完整系統的保結構算法,并證明了其對稱性,說明了其收斂性和長時間保結構的特性.
對Hamel 變分積分子的現有研究結果在數值上充分表現出了其對非線性系統,包括有限維系統、無窮維系統、完整系統與非完整系統,長時間動力學行為的良好預測[7,17,20].本文以幾何精確梁為例,利用離散Noether 定理,研究了Hamel 場變分積分子的保動量性質,以期為進一步研究其保結構性質提供參考依據.
考慮一個在空間中運動的質量均勻的彈性梁.設{E1,E2,E3} 為物質標架的一組基,初始時刻梁沿E2軸水平自然放置(無彈性形變),一端位于原點(見圖1).設梁長為l,密度為ρ,截面面積為A的正方形.

圖1 幾何精確梁初始位形Fig.1 Initial configuration of a geometrically exact beam
在幾何精確梁的模型中,梁的截面做剛體運動,其位形可由

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

和對流應變

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

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

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


其中[·,·]為se(3)上的李括號.
根據變分公式和Hamilton 原理,計算得到自由幾何精確梁的動力學方程

和自由邊界條件

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

相容性條件描述了對流速度和對流應變生成的分布的可積性.在場論框架下,位形g是以時空為底流形,纖維為的主叢S E(3)的截面,相容性條件在幾何上可解釋為主叢的局部平坦性條件.對一般纖維叢情況下的相容性條件的討論參見文獻[19].對于給定初始條件和邊界條件的幾何精確梁,其運動由動力學方程(3)和相容性條件(5)共同確定,將式(3)和式(5) 構成的方程組稱為幾何精確梁的Hamel 場方程組[18-19].
通過離散Hamilton 原理得到Hamel 場變分積分子.設離散梁的空間節點數為K,空間步長為Δh,時間步數為N,時間步長為Δt.梁的位形由序列n=0,1,2,···,N-1,j=1,2,···,K描述,表示在n時刻第j個截面的位形.分別定義離散對流速度和離散對流應變如下

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


注:上式中第一項的求和范圍應為n=0,1,2,···,N-1,j=1,2,···,K,第二項的求和范圍應為n=1,2,···,N-1,j=0,1,2,···,K.為了后續計算表達式的簡潔,本文將上式求和范圍統一寫為n=0,1,2,···,N-1,j=0,1,2,···,K,將指標溢出的變量統一記為0.

及離散相容性條件

本節首先證明自由的(不受外力、外力矩) 的幾何精確梁系統具有動量守恒性質,隨后通過離散Noether 定理[31]證明Hamel 場變分積分子可以保持精確離散動量.
根據經典力學理論,一個不受外力(矩) 的完整力學系統總動量保持不變,故自由運動的幾何精確梁系統具有動量守恒性質.動量守恒往往是在空間坐標下進行驗證的,由于Hamel 場變分積分子是在活動標架下給出的,為建立其離散動量守恒,用活動標架表述的梁的動量.
根據Noether 定理,令φ:[a,b]×[0,l]×(-?0,?0)→R,?0>0 為任意一個以時空和變分變量為自變量的C1光滑的函數,并且滿足φ(a,s,?)=φ(b,s,?)=φ(t,0,?)=φ(t,l,?)=φ(t,s,0)=0,?t∈[a,b],?s∈[0,l].取定李代數中任意一個元素ξ∈se(3),構造變分曲線

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

從而,變分曲線對應的對流速度、對流應變和體坐標系下的變分分別為

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

將式(3)代入,得到

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


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


體坐標下的變分為

根據式(13)和式(14)和BCH 公式[32]可知
其中算子I:se(3)→se(3)*為

同理可得


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

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

對任意給定的ξ,由于的任意性可知

上式為式(11)的離散對應形式.
上式對j在1,2,···,K上求和并用離散邊界條件(9)可得

其中I*為I 的伴隨算子.
綜上,可以得到如下定理定理3.1.
幾何精確梁的Hamel 場變分積分子(8) 和(10)在邊界條件(9)下滿足如下定義的離散動量守恒

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

這里的首項可作為連續動量(11) 的直接離散,從而Hamel 變分積分子在O(Δt) 誤差范圍內可以保持此離散動量的首項[7].
下面通過兩個例子分別從分析和數值角度驗證上述結論.
例一考慮一個沿E3方向以v0初速度做勻速直線運動的幾何精確梁,不妨設梁長l=1,梁的質量和轉動慣量均為1,那么系統總動量非零分量為Jlin3=v0.
在Hamel 框架下,梁的位形可寫作

其中y,z分別為截面中線在E2,E3方向的坐標.從而對流速度有表達式

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

考慮系統的離散動量,由于在此場景下非平凡的對稱群為沿E3軸做平移運動的加法群,指數映射exp退化為線性映射,于是可得

E3方向上離散動量表達式為

這與空間坐標系下的經典動量表達一致.
例二考慮初始位形如圖1 所示的不受外力的幾何精確梁.其參數[33]如下:梁長l=1,橫截面是邊長為0.1 的正方形,密度ρ=1000,楊氏模量E=107,泊松比ν=0.35.
取時間步長Δt=10-4s,空間節點數K=10,給定梁的初始對流速度為


初始對流應變為

通過Hamel 場變分積分子(8)~(10)迭代計算10 000步,得到幾何精確梁的角動量和線動量如圖2 和圖3 所示,說明了Hamel 場變分積分子可以長時間以O(Δt)精度保持系統角動量和線動量.

圖2 梁的角動量Fig.2 Evolution of angular momentum

圖3 梁的線動量Fig.3 Evolution of linear momentum
對于不受外力、外力矩的力學系統,動量守恒是一個基本的守恒律.對幾何精確梁模型,傳統方法往往將其看作大量剛體的串聯或位形空間為無窮維流形的無窮維力學系統來研究其動量,而本文在協變場論觀點下,從時空角度用活動標架法研究Hamel 場變分積分子的保動量性質,進一步說明了Hamel 場變分積分子保結構的性質,同時也間接說明其具有長時間數值穩定性.本文所給的證明方法也同樣適用于一般經典場論場景下的Hamel 場變分積分子.下一步將探究Hamel 形式及變分積分子的保能量、保辛結構等保結構性質,并將其應用到殼、膜等更復雜力學系統的建模與仿真中.