基于李群局部標架的多柔體系統動力學建模與計算1)

2021-03-24 06:11:32劉鋮胡海巖

力學學報 2021年1期

劉鋮胡海巖

(北京理工大學宇航學院力學系,北京 100081)

引言

多體系統動力學的研究對象是由多個具有運動學約束、存在大范圍相對運動的剛性或柔性部件組成的復雜系統,主要研究內容是這類系統的動力學建模、計算和控制.該學科與多個相鄰學科具有交叉.例如,計算結構力學領域的學者研究柔性結構的大變形動態行為,這相當于一類特殊的多柔體系統動力學問題.

在歷史上,多體系統動力學與計算結構力學沿著各自的途徑發展.早期,多體系統動力學研究針對多剛體系統,關注如何描述剛體大范圍運動以及剛體間的連接關系,即如何表征系統慣性力以及約束力等相互作用力;計算結構力學則研究無剛體運動的復雜結構,關注如何近似描述結構變形及表征結構的內力.由此,兩個學科的研究思路以及研究方法存在明顯區別,導致了一定的學科壁壘.

隨著學科的拓展,上述兩個學科均涉足如何描述具有大變形的柔體動力學問題,在結構力學領域,將這類問題稱為幾何非線性問題.從嚴格意義上看,計算結構力學中的幾何非線性是由于柔體大變形引起的轉動所致.而多柔體系統動力學中的幾何非線性不僅包括柔體大變形引起的轉動,更包含由柔體大范圍剛體運動引起的轉動,其非線性程度遠高于計算結構力學中的幾何非線性問題.因此,從幾何非線性角度來看,多柔體系統動力學建模和計算的難度顯著高于計算結構力學.

梁、板/殼結構是多柔體系統的基本部件.多柔體系統動力學建模的核心問題是:如何精確描述梁、板/殼結構的大范圍剛體運動與彈性變形的耦合.長期以來,為建立梁、板/殼結構發生大轉動時的精確動力學模型,計算結構力學領域與多柔體動力學領域的學者分別提出了多種有限元建模方法,例如幾何精確法、絕對節點坐標法、共旋坐標法等.其中,結構力學領域提出的幾何精確方法(geometrically exact formulation,GEF)[1]與多柔體動力學領域提出的絕對節點坐標法(absolute nodal coordinate formulation,ANCF)[2]應用較為廣泛.從本質上看,這兩類建模方法均屬于非線性有限元方法.前者中早期幾類單元屬于含轉動參數的有限元法,后者屬于無轉動參數有限元方法.關于幾何精確梁建模方法的系統性綜述可參見文獻[3];關于ANCF 法的系統性綜述可參見文獻[4].文獻[5-7]表明,這兩類建模方法在計算精度、計算效率以及實現難易程度上有明顯差異.

計算結構力學中的幾何精確法源自Reissner[8]對大變形梁的幾何精確描述.隨后,經Simo 等[1]完善,形成了R3×SO(3)群中的幾何精確梁模型,也被稱為一維Cosserat 桿模型.該方法用中線位置矢量與截面轉動偽矢量作為節點參數來描述三維梁的大轉動和大變形,中線位置矢量屬于R3線性空間,截面轉動旋轉矩陣屬于SO(3)李群,可視為9 維線性空間中的三維流形.基于截面轉動,在梁中線上建立一個局部標架,在該局部標架下描述梁單元的應變與角速度.由于該方法在局部坐標系中完成轉動的更新,從而可有效規避轉動參數的奇異性問題.但由于轉動參數不屬于線性空間,轉動參數的時空離散、線性化以及更新等算法與傳統有限元存在很大差異,導致該算法比較復雜.這在一定程度上制約了該方法的推廣.考慮到在R3×SO(3)群中建模的復雜性,已有學者構造了R3×R3空間的幾何精確梁單元[9-10].即認為轉動參數屬于R3空間,采用線性空間運算進行轉動變量的線性化及更新,降低了幾何精確梁單元的復雜性.需要指出的是,上述建模所采用的構型空間并不影響單元收斂速度與精度,單元精度僅與應變及其一次變分的計算方法相關.然而,R3×R3空間的幾何精確梁單元無法避免轉動參數的奇異性問題,也難以通過局部標架方法在計算中來消除剛體運動導致的幾何非線性.

最初,Simo 等[1]提出的幾何精確梁模型在全局坐標系中求解位移和轉動偽矢量的增量,梁的力平衡方程包含剛體運動帶來的幾何非線性.Cardona 等[11]將力平衡方程的線性化過程從當前構型切平面轉換到前一收斂步的切平面,并在前一收斂步的構型空間的局部坐標系中求解轉動偽矢量增量.該方法不僅能夠得到對稱的切向剛度矩陣,而且能同時降低剛體運動導致的幾何非線性程度.在數值計算中,這表現為剛度矩陣中的轉動部分可近似滿足剛體轉動的不變性.Sonneville 等[12]進一步采用局部標架描述幾何精確梁的平動,在SE(3)群內統一描述梁的平動與轉動,包括統一插值、線性化、更新等.上述研究指出,基于SE(3)群局部標架的幾何精確梁單元能夠在計算中消除單元剛體運動帶來的幾何非線性,其切線剛度矩陣滿足剛體運動的不變性.相比于R3×SO(3)群中的幾何精確梁方法,在SE(3)群中的建模方法在局部標架中計算線速度及其增量.由于在不同時刻線速度并不屬于同一坐標系,其對應的位移增量與截面轉動偽矢量增量需要在SE(3) 群內統一采用乘法更新.劉鋮等[13]的研究表明,雖然在SE(3)群中統一插值能消除梁單元剪切閉鎖,但對平動與轉動獨立插值以及采用縮減積分技術可進一步提高梁單元的收斂性,并能夠兼容等幾何分析[14]思想.此外,局部標架建模方法還可用于消除部分非完整約束.例如,對于經典力學中的雪橇運動,其質心速度在連體坐標系下一坐標軸方向速度投影為零.若選擇合適的局部標架,可直接消除該自由度,無需施加非完整約束.上述性質也可用于構造可精確描述大轉動大變形的五自由度板殼單元,消除了一自轉(drilling)自由度.總體來看,基于局部標架的幾何精確梁研究尚處于初步階段,例如,不同運動參數化方法對單元計算精度與效率的影響、如何構造局部標架Kirchhoff幾何精確梁單元以及如何構造高頻耗散的時間積分方法等方面尚無相關研究.

幾何精確板殼單元的發展略滯后于幾何精確梁單元,但發展歷程類似.Simo 等[15-17]發表系列論文,在R3×S2群內基于Reissner-Mindlin 中厚板理論提出了考慮剪切的五自由度幾何精確板殼模型,包含3個平動自由度與2 個旋轉自由度.該單元采用中面位置矢量與沿厚度方向的單位矢量作為廣義坐標描述板殼運動,其中位置矢量屬于線性空間R3,而單位矢量屬于二維球面流形S2.為了定義厚度方向單位矢量唯一的旋轉矩陣,將單元厚度方向的單位矢量約束在S2中并沿測地線進行旋轉,約束節點角速度矢量不含沿厚度方向分量,得到相應的旋轉偽矢量.在局部標架中,可自然消去自轉自由度.并可保證沿厚度方向的矢量是單位矢量,不需要增加額外的約束方程.但五自由度板殼單元難以直接處理含約束的多體系統,例如,板殼結構與梁結構的連接問題.為此,Simo 等[18-19]進一步通過單元中面變形梯度張量的極分解,建立了自轉自由度與中面運動之間的關系,基于約束變分原理提出了含自轉自由度的六自由度幾何精確板殼單元.同時,引入自轉自由度后,還可用于改進單元的收斂性,克服對非規則網格的敏感性問題[20].本文研究表明,在SE(3)群中構造的六自由度局部標架幾何精確板殼單元可自然地消除幾何非線性.然而,由于五自由度板殼單元局部標架無法描述中面自轉運動,導致只能消除部分幾何非線性.如何進一步改進五自由度板殼單元,使其同樣能夠完全消除幾何非線性值得進一步研究.

此外,在幾何精確板殼單元中,可選擇不同的轉動參數進行離散.例如,選擇轉動偽矢量增量、全局坐標系下沿厚度方向的單位矢量增量、或是局部標架下沿厚度方向的單位矢量增量等.當然,不同的轉動參數化方法導致算法復雜性及收斂性不同.Sonneville[21]對SE(3)幾何精確板單元建模方法進行了初步研究,但所構造的板單元在收斂性方面存在一定問題,同時單元適用性方面尚不完善.目前,關于SE3 群中的幾何精確板殼單元研究也處于初步階段.

在多柔體系統動力學領域,Shabana[2]于1996 年提出了描述柔體大范圍運動與大變形耦合的絕對節點坐標方法.該方法在處理柔體轉動時,摒棄復雜的轉動參數,采用節點位置矢量以及沿物質標架的斜率矢量來描述柔體運動.由于上述矢量均屬于線性空間,絕對節點坐標法通俗易懂,在多柔體系統動力學領域受到廣泛關注.然而,該方法在計算效率和收斂性方面存在一定的缺陷,在非線性有限元領域得到的關注較少.在該方法中,根據對單元變形的假設,可將絕對節點坐標單元分為全參數(fully parameterized)單元[22]與縮減(slope deficiency)單元[23]兩類.例如,全參數梁單元的每節點含12 個廣義坐標,通過單元中線/面斜率矢量以及截面/厚度方向斜率矢量描述單元的運動與變形,采用三階Hermite 插值對中線/面位置矢量進行插值,同時采用一階Lagrange 插值對截面/厚度斜率矢量進行插值,能夠描述梁的截面變形,可認為屬于一類實體梁單元;縮減梁、板殼單元則分別采用Euler-Bernoulli 梁和Kirchhoff板假設,描述細長梁與薄板殼.此類單元忽略橫向剪切變形,采用中線/面位置矢量與斜率矢量描述單元運動.與上述思路類似的剛體動力學建模方法為自然坐標方法[24],它直接采用正交矩陣元素與剛體任意點的位置矢量作為參數描述剛體運動.由于在同一慣性坐標系下描述多體系統運動,此類建模方法較為簡便.自然坐標方法與絕對節點坐標方法統也因此被統稱為絕對坐標方法.然而,多柔體系統中柔體大范圍剛體運動所帶來的幾何非線性將完全保留于系統動力學方程,會給動力學求解帶來不便.

基于局部標架思想,可在絕對節點坐標方法單元的任意物質點上構造局部標架,從而用李群方法將現有全局坐標系下的單元改造為局部標架下的絕對節點坐標單元.事實上,基于局部標架思想可以改造幾乎所有考慮幾何非線性的梁、板殼以及實體單元,使其在計算中消除剛體運動導致的幾何非線性,同時繼承該單元的原有特征,例如單元收斂性、閉鎖處理方法等.

本文以梁、板殼結構為例,闡述如何發展一套基于李群局部標架(local frame of lie group)的多柔體系統動力學建模和計算方法體系.該體系不同于以往在慣性坐標系下(如絕對坐標系)或在柔體的某個特殊連體坐標系下(如浮動坐標系) 描述柔體運動,而是在柔體任意點的李群局部標架中表征該點的運動和變形,通過李群運算完成柔體物理量在局部標架與全局標架之間的轉換,進而建立和求解多柔體系統動力學方程.

上述方法體系的特點是,對于柔體的大范圍運動和大變形耦合問題,可在計算中消除物體整體運動中所包含的剛體運動,從根本上規避剛體運動導致的幾何非線性.該方法同樣適用于大變形結構力學問題,可消除由大變形所帶來的大轉動剛體運動.消除剛體運動后,多柔體系統的有限元模型與大變形結構力學有限元模型僅含以彎曲應變為主導的幾何非線性,具有類似的數值特性.例如,單元廣義慣性力、彈性力及其雅可比矩陣均滿足對任意剛體運動的不變性,單元幾何剛度矩陣的貢獻可弱化.因此,具有大范圍剛體運動與大變形耦合特征的多柔體系統動力學與僅考慮大變形的計算結構力學可趨于統一.這有望推動多柔體系統動力學與計算結構力學的相互融合,在此基礎上形成新一代的多柔體系統動力學建模和計算軟件.

1 基于SE(3)群局部標架的幾類梁單元

本節以含轉動參數的幾何精確梁與幾種無轉動參數的梁單元為例,闡述基于SE(3)群局部標架的梁單元建模方法.

1.1 幾何精確Timoshenko 梁單元

現以圓截面梁單元為例,給出基于SE(3)群局部標架的幾何精確梁構造方式.該方法可方便地推廣至構造任意截面曲梁單元或變截面梁單元.在幾何精確梁模型中,采用Timoshenko 梁模型假設,即梁具有剛性橫截面,梁的構型由梁的中線位置與剛性橫截面轉動唯一確定.

首先,在三維歐氏空間R3中建立正交的慣性坐標系(O-e1-e2-e3),描述圖1 所示空間梁由初始構型到當前構型的運動,包括梁的剛體運動和變形之耦合.

圖1 三維幾何精確梁的初始構型與當前構型Fig.1 Initial and current configurations of a geometrically exact beam element of three dimensions

不失一般性,設梁在初始構型時的中線沿著慣性坐標系的X軸方向,采用R3中的矢量φ0描述梁的中線位置.以矢量φ0的端點為原點,建立梁的局部連體坐標系,其方向與慣性坐標系一致,記其正交基矢量為此時,局部連體坐標系與慣性坐標系之間的旋轉矩陣可記為R0=I3.

當梁抵達當前構型時,描述梁中線的位置矢量φ0變為矢量φ,上述局部連體坐標系隨著梁的剛性橫截面轉動,其基矢量變為正交基矢量(i1-i2-i3),局部連體坐標系與慣性坐標系間的旋轉矩陣為R.本文簡稱上述局部連體坐標系為局部標架.

基于上述局部標架和SE(3)群理論,梁單元的中線位置矢量及梁的橫截面旋轉矩陣合并為運動張量H,可表示為

其中,ξ1∈[0,L0]表示梁中線上P點在初始構型中的弧長坐標,L0為梁單元長度,u為該點在初始構型局部標架中的位移矢量,?R為初始構型與當前構型之間的旋轉變換矩陣.H矩陣稱為歐幾里得變換,它包含了梁截面作剛體平動和轉動的完整信息.需要指出的是,SE(3)群中的歐幾里得變換矩陣除了具有式(1)這樣的形式,還可表示為六階方陣形式[25].

對于式(1)中在局部標架描述的運動增量,可通過SE(3)群的指數映射,轉換為初始構型下的位移矢量與旋轉矩陣,其表達式為

其中,Θ為描述剛體轉動的偽矢量,符號帽子映射代表對應·的反對稱矩陣,為在局部標架中度量的位移矢量.在SE(3)群中,上述指數映射可表示為

進一步,在初始構型與當前構型中,梁單元上任意點的位置矢量可分別表示為

其中,ξ2與ξ3為梁截面上P點在局部標架中的坐標分量.

對于動力學問題,梁單元在任意時刻的中線速度和截面轉動角速度,可通過運動學方程,也被稱為Poisson 方程,表示為

其中,頂標表示矢量函數對時間t的導數,U與ω 為梁單元中線和橫截面的線速度與角速度在局部標架下的投影.此時,慣性力所做的虛功可表示為

其中,v=δd和δ η 是梁單元中線和橫截面在局部標架下的虛位移及虛轉角,ρA表示單位長度梁的材料密度,J為單位長度梁的慣性矩陣.上述轉動項與傳統方法一致,但由于采用局部標架描述平動,平動慣性力的形式與轉動慣性力的形式一致,存在轉動與平動的耦合項,與絕對導數與相對導數運算相關.

根據連續介質力學,梁單元內任意點的變形梯度張量可表示為

其中,()′表示向量函數對ξ1的偏導數.對于小應變問題,可根據Timoshenko 梁的變形假設,將梁單元內的應變表示為如下矢量形式

其中,axial()表示反對稱矩陣對應的軸矢量,Γ表示梁中線的拉伸應變和橫向剪切應變,κ 的分量κ1表示梁的橫截面扭轉應變,分量κ2與κ3分別表示橫截面繞i2與i3軸的彎曲應變.對于線彈性材料,通過平面應力假設以及截面積分,可將梁單元內力所做的虛功表示為

其中,N=D1Γ 表示局部標架下的梁截面內力,M=D2κ表示局部標架下的梁截面內力矩.D1與D2為彈性矩陣,可分別表示為

其中,E和G分別為材料的拉伸模量和剪切模量,A,I1,I2分別表示梁截面面積和截面慣性矩,Js表示圣維南扭轉常數,ks1和ks2為相應的剪切修正系數.相應的,外力所做的虛功表示為

其中,P表示在局部標架下梁中線上的分布力,T表示在局部標架下梁截面上的分布力矩,表示加載于梁中線的集中力,表示加載于梁截面的集中力矩.

由式(6)、式(9)以及式(11),可構造無約束梁單元的動力學方程.在此基礎上,通過對梁單元的動力學方程進行空間域有限元離散以及時間域差分離散,可得到一個非線性代數方程組,進而進行數值求解.

在對上述動力學方程進行空間離散時,如何保證插值算法的客觀性至關重要.例如,Jeleni?與Crisfiled[26]指出,若直接用有限元方法對轉動偽矢量插值,將不滿足客觀性條件;即單元剛體轉動將產生虛假應變能,應變更新算法則會對超彈性材料帶來路徑相關性問題.隨后,他們借鑒共旋坐標系概念,提出一種相對插值方法,解決了幾何精確建模方法插值的客觀性問題.Ibrahimbegovic 等[27]指出,若采用更新算法計算單元應變,直接對轉動增量進行插值,同樣能夠滿足離散應變的客觀性.近年來,Bauchau等[28]系統地討論了幾類常用轉動/運動插值算法是否滿足客觀性、矢量性、本征性以及它們的計算效率.其中,插值算法的本征性概念是指插值算法是否依賴某類特殊的轉動參數化方法.

此外,在構造離散系統動力學方程時,需要對應變進行一次變分.在線性空間內,有限元離散與變分操作可交換順序.但在SE3 群中,兩者交換順序并不等價.對于先離散后線性化,一般稱為離散一致線性化;而對于先線性化后離散,則稱為連續一致線性化.為簡單起見,下面給出連續一致線性化過程.通過SE(3)群內的變分運算規則,將幾何精確梁單元的應變矢量在局部標架下進行變分,可表示為

其中,I6為6×6 的單位矩陣.經有限元離散后,內力所做虛功的離散形式可表示為

其中,Ψ為有限元離散形函數矩陣.本文選用一階Lagrange 插值進行離散,并采用選擇性縮減積分處理單元的剪切閉鎖問題.為了得到切線剛度矩陣,再進一步作二次變分(即線性化)可得

其中,第一項包含離散形式的材料剛度矩陣、第二項包含離散形式的幾何剛度矩陣,

從材料剛度矩陣與幾何剛度矩陣的表達式可見,其僅與局部標架下的應變相關,而與剛體運動無關.因此,單元剛度矩陣自然滿足剛體運動的不變性,也就是消除了剛體運動導致的幾何非線性.值得指出,由于SE(3)群中的運算不具有交換性,此處的單元切線剛度矩陣是非對稱的.當然,也可通過轉換至前一收斂步切平面的方式構造對稱形式的剛度矩陣,其具體推導可參見文獻[11].由于上述結果消除了剛體運動導致的幾何非線性,而在多數動力學計算問題中,幾何剛度矩陣可忽略不計.此時,可獲得一個近似的對稱剛度矩陣.值得指出,相比于R3×R3線性空間的幾何精確梁單元,基于SE(3)群局部標架的空間梁單元彈性力及其剛度矩陣表達式非常簡潔,可顯著提高計算效率.

消除剛體運動產生的幾何非線性后,梁的幾何非線性主要由其彎曲應變引起.而隨著單元數目的增加,積分形式的彎曲應變逐漸減少,應變一次變分B矩陣趨于常數矩陣,材料剛度矩陣也將趨于一常數矩陣,梁的幾何非線性程度進一步減弱.因此,對于足夠稠密的有限元網格,系統剛度矩陣可近似為一常數矩陣,與線性有限元具有類似的數值特性.然而,在工程實際問題中,在保證空間離散收斂前提下,為平衡剛度矩陣的更新次數與線性方程組的求解規模,應選擇合適的有限元網格尺寸.上述性質對于任意的局部標架幾何非線性有限單元都是適用的.

在時間離散方面,幾何精確方法最為常用的是李群a 類算法與李群保能量?(角)動量算法.為簡單起見,此處給出對應SE(3)群描述的廣義a 方法的基本離散格式.李群廣義a 方法對前后兩個相鄰時刻運動增量進行差分離散,并通過SE(3)指數映射遞推運動張量H,其具體離散格式為

其中,dn表示系統從n時刻到n+1 時刻在局部標架下的運動增量,αm,αf,β,γ 為廣義a 方法算法的參數,?t為離散的時間步長,Hn=由于通過前后兩個時間步增量進行系統構型更新,運動增量dn的范數通常為小量,不會出現轉動參數的奇異性問題,也可避免轉動參數在周期臨界位置的特殊處理.系統廣義質量矩陣的離散形式可表示

其中,廣義質量矩陣的質量項Mm、陀螺力項Mgyr和離心力項Mcent分別為

其中,v=由上式可知,在動力學計算過程中,時間離散后廣義質量矩陣中包含O(h?2)Mm+O(h?1)Mm+Mcent,其中主要貢獻項O(h?2)Mm為常數矩陣,陀螺項與離心力項相對貢獻較小.對于大多工程問題,陀螺力項與離心力項可幾乎保持不變.若對于極端繞三軸高轉速旋轉動力學問題,需保留陀螺項與離心力項.此外,廣義質量矩陣表達式僅與系統的慣性特性、速度以及加速度相關,與轉動參數無關.這表明,數值計算過程中不會出現奇異性問題.

值得指出,上述基于SE(3)群局部標架的幾何精確梁建模方法可退化為基于SE(3) 群局部標架的剛體動力學建模方法,且具有避免轉動參數奇異性以及降低剛體轉動非線性的優點.

1.2 無轉動參數的梁單元

對于無轉動參數的梁單元,也可建立基于SE(3)群局部標架,進而消除梁單元剛體運動導致的幾何非線性.

首先,考察基于ANCF 描述的全參數梁單元.采用中線位置矢量以及沿中線方向的斜率矢量作為廣義坐標,則梁單元中線上任意一點的位置矢量rmid可通過三階Hermite 插值函數表示為

上式中的形函數矩陣可表示為

其中,S1=1?3ξ2+2ξ3,S2=L0(ξ?2ξ2+ξ3),S3=3ξ2?2ξ3,S4=L0(?ξ2+ξ3),ξ=X/L0為單元局部坐標,L0為梁單元的初始長度,I3為三階單位矩陣.此外,通過一階線性Lagrange 插值得到梁截面的斜率矢量r,Y與r,Z.因此,該梁單元上任意點的位置矢量可表示為

其中,S=[S1I3S2I3S5I3S7I3S3I3S4I3S6I3S8I3]T,S5=Y(1 ?ξ),S6=Yξ,S7=Z(1 ?ξ),S8=Zξ,X=[X Y Z]T.該單元在軸向與橫向的插值階數不同,會產生閉鎖問題,處理方法可詳見文獻[29].

該單元上任意點的局部標架可通過該點的3 個斜率矢量確定.雖然局部標架的構造方法不唯一,但這些方法的數值特性類似.例如,可采用如下局部標架

由于該局部標架完全由梁的平動位移場確定,因此局部標架所對應的旋轉矢量變分可表示為

類似地,可計算出局部標架的角速度以及角加速度,進而可采用基于SE(3) 群描述的廣義a 方法求解系統動力學方程.考慮到篇幅限制,此處對該單元的廣義質量矩陣、切向剛度矩陣等不加以推導,部分細節及切向剛度矩陣高效計算方法可見文獻[30].

除上述全參數ANCF 梁單元外,其余幾類梁單元均采用Euler-Bernoulli 梁的變形假設.若采用位移場的二階導數計算連續形式的彎曲應變,則單元構造需要采納至少C1連續的插值函數.對于ANCF 縮減梁單元,為滿足上述要求,其節點坐標包含節點位置以及節點沿軸向的偏導數矢量,采用三階Hermite插值離散梁中線位移場,如式(18)所示.任意點的局部標架可通過該點在參考構型中的斜率矢量以及在當前構型中的斜率矢量定義為

其中,E=?為并矢符號.關于ANCF 縮減梁單元構造可見文獻[31].若進一步需要考慮梁扭轉變形,或矩形截面梁的彎曲變形,細節推導可見文獻[32],其中局部標架則需在式(23)基礎上考慮扭轉運動.

此外,為減弱對梁中線位移場高階連續性要求,減少單元廣義坐標個數,有學者提出了離散彈性桿(discrete elastic rods,DER)模型[33].該模型通過離散方式描述桿的軸向、彎曲以及扭轉應變,具有計算效率高的優點,在計算幾何領域被廣泛應用.值得指出,該模型的彎曲變形描述與Euler-Bernoulli 梁一致,只是中文直譯為離散彈性桿單元.當不考慮扭轉變形且網格均勻劃分時,桿單元節點處的局部標架也表示為式(23)構造,式中的單位矢量E與t定義如下

近年來,隨著CAE 與CAD 的融合,等幾何分析[14]受到廣泛關注.等幾何分析將構造CAD 幾何造型的樣條函數(如Bézier、B-spline、非均勻有理B樣條NURBS 和T-spline 等)及其控制點與權重直接作為有限元建模的輸入信息.考慮到三階NURBS 樣條在特殊情況下可退化為三階Hermite 插值函數,針對一類特殊等幾何梁單元(三階NURBS 插值,除首尾節點重復度為四外,其余節點矢量重復度為二),可類似對上述ANCF 縮減梁單元的討論來構造形如式(23)的局部標架,對于由一般NURBS 樣條函數離散的梁單元,如何構造局部標架是個難題,尚值得進一步研究.

2 基于SE(3)群局部標架的幾類板殼單元

本節基于SE(3)群局部標架,介紹含轉動參數的幾何精確板殼單元以及幾類無轉動參數的板殼單元.

2.1 三維幾何精確Reissner-Mindlin 板殼單元

Simo 等[15]首先提出了R3×S2幾何精確板殼單元,隨后,并將該單元推廣至可考慮變厚度板殼單元[16]、以及考慮彈塑性本構的板殼單元[17].本文進一步將上述單元推廣至SE(3)群局部標架中,消除剛體運動導致的幾何非線性.

圖2 幾何精確板殼單元中面的初始和當前構型Fig.2 Initial and current configurations of a geometrically exact plate/shell element of three dimensions

在幾何精確板殼理論中,假設板的構型為不可伸展的單方向Cosserat 曲面,即板的厚度不發生變化,沿厚度方向的單位矢量屬于S2空間.在初始構型與當前構型中,板單元上任意點的位置矢量可分別表示為

其中,局部標架對應的旋轉矩陣R可通過厚度方向的單位矢量確定,如式(23).E=t0為初始構型中沿厚度方向的單位矢量,t為當前構型中沿厚度方向的單位矢量.由式(23)可見,局部標架在時間域的演化并非自由運動,其轉動不含沿厚度方向的分量.換言之,沿著該方向的轉動增量與角速度分量均為零.在局部標架中,可根據該約束直接消除一個自轉自由度,即單元的節點自由度為5.

根據連續介質力學,板殼中任意點的應變可分解為薄膜應變ε、彎曲應變κ 以及橫向剪切應變γ,并可在局部標架中分別表示為

此時,式(26)退化為文獻[21]中的應變表達式,彎曲應變與薄膜應變解耦,適用于描述以彎曲變形為主導的板殼結構.若采用式(27) 的彎曲應變,則在板殼單元純彎曲測試中單個單元即可達到解析解精度.進一步,通過與1.2 節類似的時空離散方法、一次變分以及兩次變分等操作,可構造出基于SE(3)群局部標架的五自由度幾何精確板殼單元.

若進一步通過中面變形梯度張量的極分解,建立自轉自由度與中面運動之間的關系,可推導出SE(3)群局部標架的六自由度幾何精確板殼單元.此時,局部標架坐標旋轉矩陣Q可表示為

其中,θ 表示為自轉角度.對于上述六自由度板殼單元,由于所建立的局部標架包含該點的完整的轉動信息,該單元能夠完全消除剛體運動帶來的幾何非線性.而對于五自由度板殼單元,若直接采用式(25)中旋轉矩陣R所定義局部標架則僅能夠消除部分的幾何非線性.然而,五自由度單元的缺陷也可通過構造近似的自轉轉動進行彌補.考慮到篇幅限制,此處不對板殼單元給出具體推導,將在作者后期論文中進一步展示.

對于低階板殼單元,單元存在嚴重的薄膜與剪切閉鎖.本文分別采用Hellinger-Reissner 變分原理和Hu-Washizu 變分原理處理薄膜與剪切閉鎖.相比于文獻[21]研究工作,本文所構造的板殼單元具有更好的數值穩定性,位移場時空離散方法也可進一步提高單元收斂性.

2.2 幾類無轉動參數的板殼單元

現基于SE(3)群局部標架,簡要介紹幾類無轉動參數的板單元建模方法,包括ANCF 全參數板殼單元、ANCF 縮減薄板殼單元、以及等幾何薄板單元.

基于局部標架的ANCF 全參數板單元構造方法與ANCF 全參數梁單元幾乎一樣,此處不再贅述.關于全參數ANCF 板單元建模與高效計算方法可見文獻[34].對于ANCF 縮減薄板殼單元,作者曾基于Kirchhoff假設,引入一個局部笛卡爾標架來描述板殼單元的應變[35].事實上,該局部笛卡爾標架可直接作為局部標架,進而描述板殼單元的速度以及角速度.該局部標架的具體表示為

上述方法可推廣至等幾何薄板殼單元中.與處理等幾何細長梁類似,現考察一類特殊的三階雙變量NURBS 離散方法,其中首尾節點四次重復,其余節點二次重復.該單元可等價于48 自由度的ANCF 縮減板殼單元,其中有限元節點自由度包含2 個一次偏導數以及1 個二次混合偏導數.由此可通過式(29)構造一個局部標架.可證明,在該局部標架下建立的動力學方程能夠規避剛體運動導致的幾何非線性.

事實上,基于SE(3)群局部標架的描述適用于所有描述幾何非線性問題的有限單元.只要將原有有限單元拓展到SE(3)群中,即可實現多柔體系統動力學與現有非線性有限元法的融合.

3 基于局部標架的多柔體系統動力學核心算法

多體系統動力學是典型的交叉學科,其研究內容涉及多剛體系統動力學、結構靜/動力學、分析力學、計算數學等學科.歷史上,多柔體系統動力學與結構靜/動力學沿著各自的途徑發展.因此,在大型CAE 軟件中,多柔體系統動力學模塊與結構靜力/動力學模塊之間并不融合.例如,在MSC 軟件體系中,解決大轉動、大變形、剛柔耦合動力學問題時需采用多體動力學軟件MSC.ADAMS 與有限元軟件MSC.NASTRAN 進行聯合仿真.其原因是,在多體系統動力學軟件開發之初,僅有多剛體系統動力學作為理論基礎,并未將其與有限元理論相融合.現有的其他幾類多體系統動力學軟件也存在類似的問題,這給多柔體系統的動力學計算帶來極大不便.

基于上述SE(3)群局部標架,融合計算幾何力學與非線性有限元理論,有望發展一套新的多柔體系統動力學建模與計算方法體系;結合等幾何分析,可開發一類消除剛體運動非線性的剛體、梁、板、殼以及實體單元,實現多柔體系統動力學與結構大變形動力學建模與計算方法的統一.此外,這樣的多柔體系統動力學軟件需具備用于長歷程動態仿真且高頻耗散可控的通用時間積分算法,能模擬多柔體系統的復雜碰撞問題,同時實現大規模多柔體系統動力學計算的通用并行異步仿真,以滿足復雜系統的動力學仿真需求.上述方法體系的發展涉及如下幾個關鍵科學問題:SE(3)群中三維運動參數化及基本運算規則問題;不同構型空間建模方法之間的轉換關系;含碰撞的多柔體系統長時間計算精度問題;復雜多柔體系統負載平衡圖分解及迭代并行算法高效率預條件算子構造問題.

具體來說,上述多柔體系統動力學建模和計算方法可分解為如下6 個方面.

(1) 基于SE(3) 群局部標架的建模方法.研究SE(3) 群中運動的參數化描述方法、運動參數插值方法、動力學方程線性化運算規則;基于SE(3)群局部標架,構造含轉動參數的梁板殼單元,例如1.1 節提出的幾何精確梁單元、2.1 節提出的幾何精確板殼單元、Kirchhoff幾何精確梁/板殼單元、離散彈性桿單元;基于SE(3)群局部標架,研究考慮截面變形的任意截面的空間梁建模方法,提出截面運動與中線運動解耦的維數降階動力學建模方法;基于SE(3)群局部標架,構造無轉動參數的梁、板殼、實體單元,例如Kirchhoff板殼單元,實體單元等;分析不同運動參數化對計算精度與效率的影響;借鑒有限元方法中對閉鎖問題的處理手段(如混合變分原理),提高上述單元的收斂性.

(2) 長時間歷程積分器.保能量?(角) 動量時間積分算法的構造依賴于單元運動的參數化描述方法、時空離散格式以及力平衡方程、幾何方程、本構方程的形式.針對上述幾類基于SE(3) 群局部標架的單元,分別構造其保能量?(角)動量時間積分算法,并進一步提出高頻耗散能量衰減動量守恒算法;基于Hamilton 原理,構造一類時空統一離散的變分積分子,包括運動參數李群變分積分子、Hamel 變分積分子,探索變分積分子處理繩索與薄膜系統動力學響應的優勢.

(3)多柔體系統的碰撞模擬算法.借鑒有限元方法研究成果,研究SE(3)群局部標架下的碰撞問題罰函數方法與Lagrange 乘子法;對于碰撞問題的隱式算法,構建保能量?(角) 動量的大步長數值計算方法;對碰撞問題的顯式算法,開發異步變分積分算法,并能夠兼容顯隱式混合時間積分方法.

(4)多柔體系統的分布式通用并行異步算法.多柔體系統模型通常存在多種不同類型的單元以及運動副,提出多柔體系統加權圖的表征方法、基于圖論的負載平衡區域分解方法;提出基于區域分解的多柔體系統并行迭代算法,研究界面粗網格自適應構造方法;當界面材料存在明顯差異時,研究界面問題高效的預條件處理算子構造方法,分析區域分解對界面問題條件數的影響;針對上述區域分解算法,探索多體系統并行異步算法,不同子區域可采用獨立的時間積分步長.

(5) 基于時空后驗誤差估計的自適應算法.針對基于SE(3)群描述的α 類時間積分算法以及保能量?(角)動量時間積分算法,提出時間離散后驗誤差估計方法,并實現變步長時間積分算法;借鑒已有的有限元方法,采用超收斂修復類或殘差類方法進行空間離散后驗誤差估計,借鑒等幾何分析中分層基函數概念實現空間網格自適應;針對時空統一離散的變分積分子,借鑒時空有限元方法發展成果,進行時空后驗誤差估計及自適應算法.

(6)多柔體系統動力學降階算法.研究基于數據驅動的參數化多柔體系統動力學降階方法[36],探索局部標架建模方法在構造降階基矢量以及降階模型動力學仿真的優勢.例如,消除剛體運動的幾何非線性,多柔體系統的降階基矢量僅需描述柔體彈性變形部分,這將有望減少降階模型的基矢量個數.同時,多柔體系統的切線剛度矩陣更新次數降低,能有效提高降階模型的動力學仿真計算效率.此外,通過對比局部標架運動規律,研究單元局部標架自適應合并與分離方法,有望實現柔體的自適應降階算法.當柔性體合并至唯一局部標架時,與傳統浮動坐標降階方法實現兼容.

4 數值算例

本節通過若干數值算例,闡述上述基于SE(3)群局部標架的多柔體系統動力學建模和計算方法的可行性.

4.1 多柔體系統保能量(角)動量算法

對于線性系統,當時間積分算法譜半徑ρ ≤1 時,該算法無條件穩定.但上述性質無法直接推廣至非線性系統,評價總能量是否守恒是評價非線性系統時間積分算法穩定性的重要指標之一.而若采用傳統高頻耗散的a 類算法,例如HHT 或廣義a 方法,可能會引起非線性系統能量增加,引起時間積分算法的發散.因此,對于柔性多體系統動力學仿真,開發保能量?(角)動量算法具有極其重要意義.

Simo 等[37]首次針對R3× SO(3)群幾何精確梁模型,提出了保能量?(角)動量時間積分方法.盡管隨后研究表明[38]該算法并不能夠精確保持系統能量與角動量,但上述方法為非線性系統,尤其是基于SO(3)群的梁板殼有限單元,構造保能量?(角)動量方法提供了基本思路.本小節以SE(3)群局部標架的幾何精確梁單元為例,給出一種保能量?(角)動量時間積分算法.并通過含線性/非線性約束、小/大變形的空間雙擺動力學模型,分析其處理多體系統動力學問題的數值特性.

(1) 考察含線性約束的雙擺模型的小變形動力學問題.如圖3(a) 所示,雙擺的兩根桿分別長0.3 m和0.5 m,其寬度與高度均為0.005 m,材料參數為:ρ=2700 kg/m3,E=7.0 × 1010Pa.桿I 與地面以及兩桿間均通過球鉸連接,系統約束方程均為線性約束.在初始時刻,桿I 靜止放置于X軸上,桿II 繞B點存在初始角速度[3 0 ?4]Trad/s.系統不受外力作用,在初始速度作用下開始運動.通過基于SE(3)群局部標架的幾何精確梁單元對雙擺系統進行離散,采用保能量?(角)動量時間積分算法求解系統動力學方程,仿真時間為100 s,時間步長為h=1.0 × 10?4s,總仿真步數為106.

圖3 柔性雙擺系統的初始構型:(a)線性約束;(b)非線性約束Fig.3 Initial configuration of a flexible double pendulum:(a)Linear constrains;(b)nonlinear constrains

圖4 含線性約束雙擺模型在小變形運動中末端C 點位置Fig.4 Displacement of end point C of the double pendulum with linear constrains subject to small deformation and overall motion

圖4 給出了雙擺模型在前30 s 內末端C點的位置隨時間變化規律,其余70 s 內位置變化趨勢大體類似.由圖可見,小變形雙擺模型長時間歷程動力學響應非常復雜,運動非線性程度較高.由于強非線性系統對初值異常敏感,難以評價算法精度,而系統守恒量能夠作為評價算法精度的一個重要指標.由于雙擺系統僅在與地面連接處受到約束力作用,該系統對慣性坐標系原點的角動量以及系統總能量守恒.圖5 與圖6 分別給出了系統能量、系統沿三個方向角動量的演化規律,其中曲線上下方的數字表示該物理量波動的最大值與最小值.本質上來說,含線性約束多體系統在計算過程中可等效為一個無約束系統,半離散平衡方程為一常微分方程組.對于此類情況,時間積分算法可精確地保持系統能量及角動量,系統能量的相對誤差僅為5.0 × 10?12,系統角動量的相對誤差為7.0 × 10?10.

圖5 含線性約束雙擺模型在小變形運動中的總能量變化Fig.5 Total energy variation of the double pendulum with linear constrains subject to small deformation and overall motion

圖6 含線性約束雙擺模型在小變形運動中的角動量變化Fig.6 Angular momentum variation of the double pendulum with linear constrains subject to small deformation and overall motion

圖7 含線性約束雙擺模型在大變形運動中的總能量變化Fig.7 Total energy variation of the double pendulum with linear constrains subject to large deformation and overall motion

圖8 含線性約束雙擺模型在大變形運動中的角動量變化Fig.8 Angular momentum variation of the double pendulum with linear constrains subject to large deformation and overall motion

(2)為了考察含線性約束的雙擺模型在大變形運動中的角動量及能量守恒特性,將(1)模型中兩根桿的截面寬度與高度均改為0.000 2 m,仿真時間為10 s,時間步長為h=1.0 × 10?5s,總仿真步仍取為106步.此時,計算得到的最大變形率達到15%.圖7 與圖8分別給出系統能量與角動量的時間歷程.由圖可見,該算法同樣能精確地保持系統的守恒量,其中能量相對誤差為1.7 × 10?5,角動量相對誤差為1.0 × 10?5.從理論上來看,該系統的能量嚴格守恒.此處的誤差主要來自求解非線性代數方程,且該誤差隨著非線性方程收斂標準的降低而減少.此外,從能量的時間歷程可見,系統的低頻振動激發起高頻振動響應,而有限元模型無法精確描述高頻振動模態.同時,系統保留的大量高頻響應將導致計算過程中需要非常小的時間步長,導致計算效率低下.因此,對于多柔體系統動力學仿真,有必要研究高頻能量耗散可控的時間積分算法.

(3)本小節考察含非線性約束的保能量?(角)動量算法的數值特性.此時,半離散系統動力學方程為一指標三的微分代數方程組,非線性約束的引入將對多體系統動力學數值特性帶來顯著的改變.

如圖3(b)所示,兩桿間由旋轉鉸相連接,為一典型的含非線性約束多體系統.雙擺模型兩桿的幾何及材料屬性與模型(1) 一致,桿II 繞B點存在初始角速度[3 0 0]Trad/s.仿真時間為100 s,時間步長為h=1.0 × 10?4s,總仿真步仍取為106步.圖9 與圖10分別給出系統總能量與角動量的時間歷程.由圖可見,系統能量相對誤差為4.9×10?9,角動量相對誤差為9.1 × 10?3.非線性約束對系統總能量的守恒特性幾乎無影響.然而,這對系統角動量守恒精度帶來了較大的影響.這是由于直接求解指標三的含非線性約束微分代數方程組時,需要引入h2量級的縮放系數,導致平衡方程中廣義慣性力與廣義約束力數值精度差h2量級.而系統角動量誤差與廣義慣性力、廣義彈性力、廣義約束力的數值精度以及當前構型位置相關.隨著時間積分步長的累計,角動量計算過程中相對誤差僅為9.1 × 10?3,但數值精度仍在可接受范圍內.理論上,可通過降微分代數方程指標技術,提高角動量守恒精度.但在實際工程問題中,需統籌考慮計算效率與計算精度.此外,對于其他類型的非線性約束,都能夠構造相應的廣義約束力,使得系統能量守恒.

圖9 含非線性約束雙擺模型在大變形運動中的總能量變化Fig.9 Total energy variation of the double pendulum with nonlinear constrains subject to large deformation and overall motion

圖10 含非線性雙擺模型在大變形運動中的角動量變化Fig.10 Angular momentum variation of the double pendulum with nonlinear constrains subject to large deformation and overall motion

由于保能量?(角) 動量算法的構造依賴于運動參數化方法、系統平衡、幾何、本構以及約束方程等一些細節因素,針對其他類型局部標架單元的保能量?(角) 動量算法值得進一步研究.需要指出的是,對于非保守系統,可將耗散力所做的功作為系統一部分能量,從數值角度亦可其視為一個守恒系統,系統能量依然是一個守恒量.例如,4.4 小節所涉及的摩擦碰撞動力學問題,保能量?(角)動量算法能夠精確表征法向碰撞勢能與切向摩檫力耗散能量,整體考慮總能量依然守恒.

此外,在計算數學領域,Marsden 等[39]提出了變分積分概念,該算法被證明能在長時間動力學仿真過程中保持系統的辛幾何結構.對于保守系統,離散系統的動量與角動量守恒,離散系統的能量在守恒量附近波動.在歐氏空間內,變分積分算法與上述保能量?(角)動量算法之間的對比可參見文獻[40].從工程應用角度來看,李群SE(3)的變分積分算法仍需要進一步發展.

4.2 廣義a 方法后驗誤差估計方法

本小節通過一個自由剛體運動問題和剛性雙擺運動問題,驗證基于SE(3)群描述的廣義a 方法誤差估計方法的有效性.

本算例將每一時間步長細化10 倍的數值仿真結果視為精確解,由此得到時間離散截斷后驗誤差的參考值.將參考誤差與后驗誤差結果進行對比,以驗證基于SE(3) 群描述的廣義a 方法誤差估計方法的正確性.圖11 給出了剛體位移與轉動的時間離散誤差估計值與參考值的二范數.由于本算例的動力學方程具有高度光滑性,誤差估計值與參考解極其吻合.

圖11 基于SE(3)群描述的廣義a 方法計算單剛體運動的后驗誤差估計Fig.11 Posterior error estimation of the generalized alpha method of SE(3)group for a single rigid body

其次,考察剛性雙擺系統,其幾何參數與材料屬性與4.1 算例(3)中模型完全一致.設該系統在沿著Z軸負方向的重力作用下開始運動.采用基于SE(3)群描述的廣義a 方法進行仿真計算,總仿真時間為1 s,時間步長為1.0 × 10?3s,算法參數譜半徑為0.8.圖12 與圖13 分別給出了桿OB 的估計誤差與參考誤差的二范數.由于后驗誤差計算方法是基于常微分方程推導,而多柔體系統約束的引入,導致估計誤差與參考誤差之間存在一定的差異,但兩者的誤差量級與整體趨勢一致,通過估計誤差進行變步長時間積分算法設計是可行的.

圖12 基于SE(3)群的廣義a 方法對剛性雙擺模型平動位移的后驗誤差估計Fig.12 Posterior error estimation of the generalized alpha method of SE(3)group for the translational displacement of a double pendulum

圖13 基于SE(3)群的廣義a 方法對剛性雙擺模型轉動位移的后驗誤差估計Fig.13 Posterior error estimation of the generalized alpha Lie group SE(3)method for the rotational displacement of a double pendulum

上述基于SE(3) 群描述的廣義a 方法的誤差估計方法可直接推廣至多柔體系統的動力學仿真中.但隨著系統剛度矩陣與阻尼矩陣的引入,時間離散誤差估計值與參考誤差將進一步擴大,但兩者的變化趨勢依然是一致的.

4.3 局部標架殼單元靜力學分析

本小節采用文獻[41]中圓柱殼和球殼靜力學模型,驗證本文所提出的SE(3)群局部標架板殼單元的正確性.

圖14 受拉載荷的圓柱殼初始構形Fig.14 Initial configuration of a pinched cylinder

如圖14 所示,一個無約束的圓柱殼在A,E兩點受到一對大小相等、方向相反的集中力F=40 kN的作用,其中A點為圓柱體上母線的中點,E點為下母線的中點.圓柱殼的長度L為10.35 mm,內徑R與厚度分別為4.953 mm 與0.094 mm,材料楊氏模量E=10.5 × 106N/mm2,泊松比υ=0.312 5.考慮到圓柱殼模型的對稱性,只需要對該模型上半部分的四分之一(ABCD)進行建模.現將在圓柱殼上施加外力F 的過程平均分成100 個增量步,并設靜力學迭代容許誤差為1.0 × 10?6.

本小節分別采用局部標架GESF 殼單元與局部標架縮減ANCF 殼單元離散ABCD圓柱殼,并分析了24×16 與36×24 的兩種網格數目下圓柱殼的靜力學變形.圖15 中給出了殼體上A,B和C三點的位移隨外力增加的變化曲線.當載荷接近20 kN 時,圓柱殼發生了屈曲變形,屈曲前圓柱殼以彎曲變形為主,屈曲后圓柱殼以拉伸變形為主.因此,本算例中GEF 殼單元不能對忽略薄膜應變對彎曲變形的影響.數值結果表明,上述兩類殼單元數值結果與參考文獻均吻合較好,均能夠精確地描述圓柱殼的大變形特性.關于一般性彈性體的屈曲、后屈曲及分叉模擬算法可參見作者前期工作[42].

圖15 圓柱殼上A,B,C 點的位移大小Fig.15 Magnitudes of displacements at nodes A,B and C of a pinched cylinder

其次,本節對圖16 所示的頂部存在18?開口薄球殼進行靜力學計算,以說明SE(3) 局部標架幾何精確殼單元同樣能夠適用于描述薄殼結構的大變形.該球殼的半徑為10.0 mm,楊氏模量為E=6.825 × 107N/mm2,泊松比為υ=0.3,厚度為t=0.01 mm(R/t=1000).該球殼在底部的A,B,C,D四點分別受到集中力F=2λF0的作用,其中F0=1.0 N,λ 取為2.5.考慮到球殼的對稱性,現對其四分之一的模型進行建模.本小節分別采用局部標架GESF 殼單元與局部標架縮減ANCF 殼單元離散四分之一開口球殼,并分析了24×24 與32×32 的兩種網格數目下球殼的靜力學變形.圖17 給出了不同單元數目下球殼上點A與點B的位移隨外載荷變化曲線.數值結果表明,上述兩類殼單元數值結果與參考文獻均吻合較好,本文提出的兩類單元能夠精確描述薄殼結構的大變形特性.

圖16 頂部18?開口球殼初始構形Fig.16 Initial configuration of a hemispherical shell with 18?hole

圖17 開口球殼上A、B 點的位移大小Fig.17 Magnitudes of displacements at nodes A and B of a hemispherical shell

考慮到板殼單元慣性力與梁單元慣性力計算方法類似,本節中不再加以分析板殼結構的動力學問題.由于局部標架的引入,板殼結構的單元質量矩陣與切線剛度矩陣滿足剛體運動的不變性.在數值計算中,六自由度的板殼單元質量矩陣與切線剛度矩陣的更新次數將急劇減少,具體數值性質可參考4.5小節中的大型桁架結構展開動力學模擬.

4.4 多柔體系統碰撞動力學仿真

本小節通過空間飛網抓捕收口動力學算例以及剛性桿盒碰撞動力學算例,展示基于SE(3)群局部標架描述的多柔體系統與多剛體系統碰撞動力學算法.

首先,將基于SE(3)群局部標架的柔性梁建模方法與罰函數算法相結合,處理碰撞問題.

圖18 為正六邊形的空間飛網,其在完全展開狀態下以一定的速度抓捕正前方的固定星體.當抓捕至適當位置時,通過收口裝置與收口繩進行縮緊網口.現通過基于SE(3)群局部標架的DER 單元對飛網系統進行有限元離散,對飛網與剛性星體、收口繩與剛性環間的碰撞,采用點?線接觸策略與保能量罰函數算法進行模擬,利用4.1 節提出的保能量?(角)動量時間積分算法求解系統動力學方程.

圖19 中給出了飛網抓捕星體以及收口階段的示意圖.在整個抓捕與收口階段,該系統總能量、動量以及角動量均守恒.計算結果表明,罰函數方法能夠較好地處理大變形柔體間的碰撞問題.若進一步考慮線?線接觸、線?面接觸問題,可采用Mortar 方法細化接觸區域,從而得到高精度的碰撞模型,其具體算法細節可見文獻[43].需要指出的是,對于此類強非線性問題,若采用傳統算法(如Newmark-β、廣義a 方法等),即使譜半徑為1,系統能量也無法守恒,甚至會出現總能量增大.

圖18 大型空間飛網展開狀態及被捕獲物體示意圖Fig.18 Deployed configuration of the large tethered-net and the target

圖19 大型空間飛網抓捕與收口動力學模擬Fig.19 Capturing and closing dynamics of the large tethered-net

其次,考察基于SE(3)群局部標架的建模方法與Lagrange 乘子法相結合,解決多體系統接碰撞動力學問題.此時,通過將碰撞時的單邊約束表示為一類互補條件,可通過互補算法求解接觸沖量.以下通過剛性桿與空心盒間的兩點碰撞算例,說明Lagrange 乘子法可直接與上述基于SE(3) 群局部標架的方法相兼容,并與罰函數法進行比較.

圖20 剛性桿與空心盒碰撞模型Fig.20 Contact model between a rigid rod and a hollow box

如圖20 所示,正方形截面桿的邊長為0.2 m.剛性空心盒的外表面長、寬、高分別為 1.8 m,1.6 m,1.4 m,厚度為0.1 m.桿與空心盒的密度均為7850 kg/m3.在初始時刻,桿兩端點在慣性坐標系下的位置坐標分別為[0 ?0.2 1]Tm 與[1 0.2 0]Tm.空心盒不受重力,桿在沿Z 軸負方向的重力作用下開始運動.現采用SE(3)群剛體單元描述剛性桿與空心盒,并通過非光滑廣義α 錐互補方法對摩擦系數為0,0.2,0.4 以及0.6 的4 種情況進行動力學仿真,算法參數譜半徑設為0.8,其中非光滑廣義a 錐互補方法可參見文獻[44].為了驗證該算法的正確性,同時采用罰函數方法對上述幾種情況進行數值仿真.

圖21 分別給出了摩擦系數為0,0.2,0.4 以及0.6時,罰函數方法與非光滑廣義α 錐互補方法得到的桿上端點z方向位移變化曲線.計算結果表明,上述兩類算法在無摩擦時的位移完全一致.當存在摩擦時,兩中方法存在微小差異.這是由于錐互補條件推導過程中引入了放松條件,導致該方法處理滑動摩擦時存在小的誤差.上述研究表明,在SE(3) 群內構造的兩種算法均能處理多柔體系統的多點碰撞問題.但兩者相比較,罰函數方法較為簡單,能方便處理點?面接觸或面?面接觸問題;而非光滑互補算法相對復雜,應用于大變形柔性間接觸問題的難度較大.然而,在處理剛性碰撞問題時,由于罰參數的引入將在系統中引起高頻振動,而采用互補條件描述的非光滑動力學能較為精確地處理剛性碰撞問題.

在本算例中,上述兩種算法均與隱式時間積分算法相結合.罰函數方法與錐互補算法同樣能應用于顯式時間積分算法,將其與異步變分積分算法結合[45],有望能處理復雜多柔體系統的多尺度碰撞問題,但此類算法的計算精度與計算效率仍需進一步研究.

圖21 不同摩擦系數時錐互補方法與罰函數方法對比Fig.21 Comparison between the CCP method and penalty method for the different friction coefficients

4.5 大型桁架結構展開動力學模擬

本小節通過一種空間桁架結構展開動力學模擬,展示基于SE(3) 群局部標架的建模方法在計算效率方面的優勢.圖22 所示為一單胞元可展開空間桁架結構,可用于構建大型空間結構系統.

圖22 單胞元空間桁架結構展開示意圖[46]Fig.22 The deployed configuration of a unit of the space truss structure[46]

該胞元結構由28 根桿通過旋轉鉸與滑移鉸連接而成,其具體描述見文獻[46].胞元的截面為邊長2 m 的正方形,橫桿長度為4 m,桿為空心圓截面,材料的拉伸模量為2.1 × 1011Pa,泊松比為0.3,密度為7850 kg/m3.設桁架結構在端部與地面固連,各胞元通過控制角度β 同時驅動展開,展開時間設為20 s,展開角的驅動為余弦控制函數.在初始時刻,展開角度與展開鎖定時刻角度分別為β0=5?與βf=90?.對每個桁架胞元,采用10 個基于SE(3)群局部標架的幾何精確梁建模,用基于SE(3)群描述的廣義a 方法進行動力學仿真,算法譜半徑為0.8.

首先,進行兩個桁架胞元的展開動力學仿真,分析展開動力學模擬中系統迭代矩陣的復用情況,以說明基于SE(3) 群局部標架的建模方法能有效消除剛體運動的非線性,極大地降低整個系統的幾何非線性.當梁截面內外直徑分別選為0.055 m 與0.08 m時,桁架胞元在展開過程中幾乎不發生彈性變形.此時,系統廣義質量矩陣與廣義剛度矩陣均無需更新,僅根據初始構型下的初值就能完成整個展開過程的動力學仿真.這意味著,此時無須計算梁單元的幾何剛度矩陣,即系統迭代矩陣為對稱矩陣.

為了考察具有大變形的多柔體系統動力學仿真過程具有上述數值特性,將梁截面內外直徑設為0.023 m 與0.024 m,進行同樣工況的動力學仿真.此時桁架胞元的最大變形率達到3.3%,在其展開過程中能夠觀察到明顯的彈性振動.

表1 中給出了兩種不同時間步長所對應的迭代矩陣K與Φq的更新次數以及總的牛頓迭代次數.其中,矩陣K與系統廣義質量矩陣與切線剛度矩陣相關,Φq為約束方程的雅可比矩陣,“NNI≥i”表示一個時間步內,單個時間步內牛頓迭代次數大于等于i步時,迭代矩陣強制更新;“Frozen”表示迭代矩陣永不更新.當牛頓迭代步數超過預設步數時,優先更新約束方程的雅可比矩陣.

表1 迭代矩陣與約束方程Jacobian 矩陣更新次數及牛頓迭代總次數Table 1 The number of times that the iteration matrix and Jacobian of constrain equations are updated as well as the number of total Newton iterations

首先考察時間步長為2.0 × 10?3s 時的工況.當NNI≥1 時,即每個時間步中牛頓迭代均更新迭代矩陣,共需進行32 410 次牛頓迭代,平均每個時間步需要3.2 次牛頓迭代;當NNI≥4 時,即若每個時間步中牛頓迭代達到四步時更新一次迭代矩陣,則迭代矩陣僅需更新19 次,約束方程雅可比矩陣需更新3083次,且總迭代次數較之前增加32%.而當迭代矩陣不更新時,總迭代次數比NNI≥4 時略有增大.由此可說明,基于SE(3)群局部標架的建模方法能極大地降低剛體運動導致的幾何非線性,迭代矩陣可視為一個常數矩陣.

再考察時間步長進一步減少的情況.此時,系統迭代矩陣中的廣義質量矩陣占比提升.一般來說,柔體的廣義質量矩陣與彈性變形相關,而切向剛度矩陣與變形梯度相關.因此,相比于切向剛度矩陣,廣義質量矩陣往往變化較小.因此,時間步長為1.0 × 10?3s,NNI≥4 時,迭代矩陣幾乎不用更新.同時,約束方程雅可比矩陣次數也大幅降低.這是由于迭代矩陣趨于常數矩陣,從而減少了對約束方程雅可比矩陣的更新次數.

在大規模的多柔體系統動力學隱式算法仿真中,最耗時的兩部分是計算系統迭代矩陣和求解高維線性代數方程組.基于SE(3)群局部標架的建模方法能最大限度地降低系統迭代矩陣更新次數,由此可將迭代矩陣提前進行計算并進行LU 分解,或用于設計迭代算法的預條件算子.相比于其他建模方法,基于SE(3)群局部標架的建模方法的計算效率顯著提高.

最后,基于上述建模方法,對具有8 個胞元的桁架結構進行展開動力學模擬.圖23 給出了不同時刻該桁架結構的展開構型.

圖23 具有8 個胞元的空間桁架展開動力學模擬Fig.23 Deployment simulation of a truss structure with eight units

4.6 基于區域分解的多柔體系統動力學通用并行算法

本小節通過環形桁架?索網天線與大型空間飛網的展開動力學問題,展示一種基于有限元網格撕裂對接(finite element tearing and interconnecting dualprimal,FETI-DP)的多柔體系統動力學通用并行算法.其中,環形桁架?索網天線展開動力學算例是為了說明并行算法能夠適用于含多種不同類型單元與多種不同運動副的復雜多體系統動力學計算;大型空間飛網展開動力學算例是為了并行算法能夠處理百萬廣義坐標量級的柔性多體系統動力學仿真.在該方法中,采用Dirichlet 預條件算子迭代來并行求解界面問題.該方法已成功應用于求解大規模結構靜/動力學問題,具體算法見文獻[47-48].

圖24 4 m 直徑的環形桁架?索網天線[49]:(a)地面實驗系統;(b)動力學模擬過程Fig.24 Deployable hoop truss antenna of 4 m in diameter[49]:(a)Experiment system;(b)dynamic simulation

圖24 是由中國空間技術研究西安分院與北京理工大學聯合研制的4 m 直徑環形桁架?索網天線模型,其幾何參數和材料參數詳見文獻[49].采用基于SE(3)群局部標架的幾何精確梁單元和DER 單元分別對環形桁架和索網進行建模,將桁架劃分為18 個單元,將索網劃分為20 個單元,共5.3 萬個自由度.在該多柔體系統動力學模型中,存在眾多不同類型的運動副,如球鉸、旋轉鉸、滑移鉸以及齒輪副.如圖25 所示,通過多層圖分解技術,將其有限元網格分解為16 區域,對不同區域通過不同顏色進行表征,實線表示索網系統,虛線表示環形桁架.經區域分解后,每個區約含3800 個自由度,在區域界面上約有1000個自由度.

圖25 環形桁架?索網天線的區域分解示意圖Fig.25 Domain decomposition of the deployable hoop truss antenna

由于在作者的前期研究中已系統分析過該天線的展開動力學特性[50],此處僅給出索網天線在無重力環境下的收回動力學過程,但仿真過程并沒有考慮索網間以及與桁架間的碰撞問題.圖26 給出了不同時刻該天線的構型圖,驗證了上述并行計算方法可適用于含復雜約束的多柔體系統動力學仿真.

最后,通過對大型空間飛網系統進行展開動力學分析,說明該并行算法能夠處理自由度達百萬量級規模的多柔體系統大變形動力學仿真問題.本算例的空間飛網拓撲結構與4.4 節中飛網一致,但是增加了幾何尺寸以及網目個數.對于此類超柔性的多體系統,只有采用較為稠密的網格才能較精確地描述其動力學特性.本研究通過基于SE(3) 群局部標架的DER 單元對該空間飛網系統進行離散,將每一段繩索離散為80 個DER 單元,系統自由度達到202萬.由于此處僅關心并行算法的可行性,故不再詳細介紹繩索的幾何參數、材料參數以及具體展開工況.通過多層圖分解技術,將該系統的有限元網格分解為72 個區域,每個區域具有2.8 萬個自由度,區域界面上具有0.5 萬個自由度.圖27 給出了空間飛網在不同時刻的展開構型.

圖27 大型空間飛網展開動力學的并行模擬Fig.27 Parallel simulation of the deployment of a large tethered-net based on the domain decomposition

需要指出的是,本算法在處理界面問題時的效率仍有較大提升空間.在粗網格構造以及界面問題高效預條件算子等方面,還需要進入深入研究.目前,僅實現了空間梁結構的展開動力學并行計算.對于大型薄膜結構的展開動力學問題,還值得進一步研究.

5 結束語

本文基于SE(3)局部標架,討論了如何發展一套新的多柔體系統動力學建模與計算方法體系.該方法體系在SE(3)局部標架下描述單元應變、應力、運動參數、運動參數變分及增量等物理量,可有效避免大范圍剛體運動導致的幾何非線性.對于數值計算而言,該方法體系具有如下優勢:通過最少轉動參數(三參數)可完全避免轉動奇異性問題;系統廣義質量矩陣的主項為常數,滿足剛體轉動的不變性,而且對大多數問題廣義質量矩陣無須更新;單元切線剛度矩陣能滿足剛體轉動的不變性,可最大限度地減少切線剛度矩陣的更新次數,有效提高計算效率.在計算精度方面,該方法體系能繼承歐氏空間中有限單元的性質,即已有提高收斂性的單元技術均能直接應用于由SE(3)局部標架描述的有限單元;由于在該局部標架下描述柔體運動,可自然消去部分自由度來滿足約束方程,例如構造五自由度的板殼單元.

上述新的多柔體系統動力學建模和計算方法體系包含如下核心方法和算法:基于SE(3)群局部標架的復雜多柔體系統建模方法;多柔體系統動力學的長時間歷程時間積分器;多柔體系統的多尺度碰撞算法;多柔體系統的分布式通用并行異步算法;基于時空后驗誤差估計的自適應算法;多柔體系統的動力學降階算法等.本文通過若干簡單算例和工程案例,簡要介紹了部分算法及其可行性.

目前,對上述方法和算法的研究仍處于初步階段,尚需進行完善并使其相互融合,進而發展成為一套新的多柔體系統動力學建模和計算方法體系,并形成相應的柔性多體系統動力學軟件.

致謝

感謝北京理工大學田強教授在絕對節點坐標全參數梁單元與廣義a 方法方面提供的幫助;感謝北京理工大學史東華副教授在幾何力學與李群理論方面提供的幫助;感謝哈爾濱工業大學任輝教授在時間積分算法誤差估計方面的幫助;感謝清華大學趙治華副教授在幾何精確板殼單元方面的幫助.感謝李培博士、羅凱博士、以及研究生侯云森、張騰、葉子晟、湯惠穎、孫德巍的辛勤工作.最后,還需要感謝北京空間飛行器總體設計部與上海宇航系統工程研究所對本文研究工作的資助.