基于李群局部標(biāo)架的多柔體系統(tǒng)動力學(xué)建模與計算1)

劉鋮 胡海巖

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

引言

多體系統(tǒng)動力學(xué)的研究對象是由多個具有運動學(xué)約束、存在大范圍相對運動的剛性或柔性部件組成的復(fù)雜系統(tǒng),主要研究內(nèi)容是這類系統(tǒng)的動力學(xué)建模、計算和控制.該學(xué)科與多個相鄰學(xué)科具有交叉.例如,計算結(jié)構(gòu)力學(xué)領(lǐng)域的學(xué)者研究柔性結(jié)構(gòu)的大變形動態(tài)行為,這相當(dāng)于一類特殊的多柔體系統(tǒng)動力學(xué)問題.

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

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

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

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

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

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

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

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

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

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

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

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

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

1.1 幾何精確Timoshenko 梁單元

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

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

圖1 三維幾何精確梁的初始構(gòu)型與當(dāng)前構(gòu)型Fig.1 Initial and current configurations of a geometrically exact beam element of three dimensions

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

當(dāng)梁抵達(dá)當(dāng)前構(gòu)型時,描述梁中線的位置矢量φ0變?yōu)槭噶喀?上述局部連體坐標(biāo)系隨著梁的剛性橫截面轉(zhuǎn)動,其基矢量變?yōu)檎换噶?i1-i2-i3),局部連體坐標(biāo)系與慣性坐標(biāo)系間的旋轉(zhuǎn)矩陣為R.本文簡稱上述局部連體坐標(biāo)系為局部標(biāo)架.

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

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

對于式(1)中在局部標(biāo)架描述的運動增量,可通過SE(3)群的指數(shù)映射,轉(zhuǎn)換為初始構(gòu)型下的位移矢量與旋轉(zhuǎn)矩陣,其表達(dá)式為

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

進(jìn)一步,在初始構(gòu)型與當(dāng)前構(gòu)型中,梁單元上任意點的位置矢量可分別表示為

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

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

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

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

根據(jù)連續(xù)介質(zhì)力學(xué),梁單元內(nèi)任意點的變形梯度張量可表示為

其中,()′表示向量函數(shù)對ξ1的偏導(dǎo)數(shù).對于小應(yīng)變問題,可根據(jù)Timoshenko 梁的變形假設(shè),將梁單元內(nèi)的應(yīng)變表示為如下矢量形式

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

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

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

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

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

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

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

其中,I6為6×6 的單位矩陣.經(jīng)有限元離散后,內(nèi)力所做虛功的離散形式可表示為

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

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

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

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

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

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

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

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

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

1.2 無轉(zhuǎn)動參數(shù)的梁單元

對于無轉(zhuǎn)動參數(shù)的梁單元,也可建立基于SE(3)群局部標(biāo)架,進(jìn)而消除梁單元剛體運動導(dǎo)致的幾何非線性.

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

上式中的形函數(shù)矩陣可表示為

其中,S1=1?3ξ2+2ξ3,S2=L0(ξ?2ξ2+ξ3),S3=3ξ2?2ξ3,S4=L0(?ξ2+ξ3),ξ=X/L0為單元局部坐標(biāo),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.該單元在軸向與橫向的插值階數(shù)不同,會產(chǎn)生閉鎖問題,處理方法可詳見文獻(xiàn)[29].

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

由于該局部標(biāo)架完全由梁的平動位移場確定,因此局部標(biāo)架所對應(yīng)的旋轉(zhuǎn)矢量變分可表示為

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

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

其中,E=?為并矢符號.關(guān)于ANCF 縮減梁單元構(gòu)造可見文獻(xiàn)[31].若進(jìn)一步需要考慮梁扭轉(zhuǎn)變形,或矩形截面梁的彎曲變形,細(xì)節(jié)推導(dǎo)可見文獻(xiàn)[32],其中局部標(biāo)架則需在式(23)基礎(chǔ)上考慮扭轉(zhuǎn)運動.

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

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

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

本節(jié)基于SE(3)群局部標(biāo)架,介紹含轉(zhuǎn)動參數(shù)的幾何精確板殼單元以及幾類無轉(zhuǎn)動參數(shù)的板殼單元.

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

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

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

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

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

根據(jù)連續(xù)介質(zhì)力學(xué),板殼中任意點的應(yīng)變可分解為薄膜應(yīng)變ε、彎曲應(yīng)變κ 以及橫向剪切應(yīng)變γ,并可在局部標(biāo)架中分別表示為

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

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

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

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

2.2 幾類無轉(zhuǎn)動參數(shù)的板殼單元

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

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

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

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

3 基于局部標(biāo)架的多柔體系統(tǒng)動力學(xué)核心算法

多體系統(tǒng)動力學(xué)是典型的交叉學(xué)科,其研究內(nèi)容涉及多剛體系統(tǒng)動力學(xué)、結(jié)構(gòu)靜/動力學(xué)、分析力學(xué)、計算數(shù)學(xué)等學(xué)科.歷史上,多柔體系統(tǒng)動力學(xué)與結(jié)構(gòu)靜/動力學(xué)沿著各自的途徑發(fā)展.因此,在大型CAE 軟件中,多柔體系統(tǒng)動力學(xué)模塊與結(jié)構(gòu)靜力/動力學(xué)模塊之間并不融合.例如,在MSC 軟件體系中,解決大轉(zhuǎn)動、大變形、剛?cè)狁詈蟿恿W(xué)問題時需采用多體動力學(xué)軟件MSC.ADAMS 與有限元軟件MSC.NASTRAN 進(jìn)行聯(lián)合仿真.其原因是,在多體系統(tǒng)動力學(xué)軟件開發(fā)之初,僅有多剛體系統(tǒng)動力學(xué)作為理論基礎(chǔ),并未將其與有限元理論相融合.現(xiàn)有的其他幾類多體系統(tǒng)動力學(xué)軟件也存在類似的問題,這給多柔體系統(tǒng)的動力學(xué)計算帶來極大不便.

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

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

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

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

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

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

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

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

4 數(shù)值算例

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

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

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

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

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

圖3 柔性雙擺系統(tǒng)的初始構(gòu)型:(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 內(nèi)末端C點的位置隨時間變化規(guī)律,其余70 s 內(nèi)位置變化趨勢大體類似.由圖可見,小變形雙擺模型長時間歷程動力學(xué)響應(yīng)非常復(fù)雜,運動非線性程度較高.由于強(qiáng)非線性系統(tǒng)對初值異常敏感,難以評價算法精度,而系統(tǒng)守恒量能夠作為評價算法精度的一個重要指標(biāo).由于雙擺系統(tǒng)僅在與地面連接處受到約束力作用,該系統(tǒng)對慣性坐標(biāo)系原點的角動量以及系統(tǒng)總能量守恒.圖5 與圖6 分別給出了系統(tǒng)能量、系統(tǒng)沿三個方向角動量的演化規(guī)律,其中曲線上下方的數(shù)字表示該物理量波動的最大值與最小值.本質(zhì)上來說,含線性約束多體系統(tǒng)在計算過程中可等效為一個無約束系統(tǒng),半離散平衡方程為一常微分方程組.對于此類情況,時間積分算法可精確地保持系統(tǒng)能量及角動量,系統(tǒng)能量的相對誤差僅為5.0 × 10?12,系統(tǒng)角動量的相對誤差為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步.此時,計算得到的最大變形率達(dá)到15%.圖7 與圖8分別給出系統(tǒng)能量與角動量的時間歷程.由圖可見,該算法同樣能精確地保持系統(tǒng)的守恒量,其中能量相對誤差為1.7 × 10?5,角動量相對誤差為1.0 × 10?5.從理論上來看,該系統(tǒng)的能量嚴(yán)格守恒.此處的誤差主要來自求解非線性代數(shù)方程,且該誤差隨著非線性方程收斂標(biāo)準(zhǔn)的降低而減少.此外,從能量的時間歷程可見,系統(tǒng)的低頻振動激發(fā)起高頻振動響應(yīng),而有限元模型無法精確描述高頻振動模態(tài).同時,系統(tǒng)保留的大量高頻響應(yīng)將導(dǎo)致計算過程中需要非常小的時間步長,導(dǎo)致計算效率低下.因此,對于多柔體系統(tǒng)動力學(xué)仿真,有必要研究高頻能量耗散可控的時間積分算法.

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

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

圖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

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

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

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

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

首先,考察自由剛體的運動.設(shè)初始時刻,剛體質(zhì)心位于慣性坐標(biāo)系原點,剛體的連體坐標(biāo)系與慣性坐標(biāo)系重合.剛體在連體坐標(biāo)系下的主轉(zhuǎn)動慣量分別為3783.42 kgm2,3776.16 kgm2,2630.94 kgm2,慣性積分別為352.9 kgm2,484.74 kgm2,469.74kgm2.初始時刻,剛體質(zhì)心速度為[10 10 10]Tm/s,剛體角速度為[10 10 10]Trad/s.此外,剛體在質(zhì)心處受到[50|sin(10t)|52|sin(10t+π/3)|50|sin(10t+π/2)|]TkN 的集中力作用,同時受到[50|sin(10t)| 52|sin (10t+π/3)|50|sin(10t+π/2)|]TkN·m 的集中力矩作用.現(xiàn)采用基于SE(3)群描述的廣義a 方法進(jìn)行仿真計算,總仿真時間為1 s,時間步長為10?3s,算法譜半徑為0.8.

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

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

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

圖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 方法對剛性雙擺模型轉(zhuǎn)動位移的后驗誤差估計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 方法的誤差估計方法可直接推廣至多柔體系統(tǒng)的動力學(xué)仿真中.但隨著系統(tǒng)剛度矩陣與阻尼矩陣的引入,時間離散誤差估計值與參考誤差將進(jìn)一步擴(kuò)大,但兩者的變化趨勢依然是一致的.

4.3 局部標(biāo)架殼單元靜力學(xué)分析

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

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

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

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

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

其次,本節(jié)對圖16 所示的頂部存在18?開口薄球殼進(jìn)行靜力學(xué)計算,以說明SE(3) 局部標(biāo)架幾何精確殼單元同樣能夠適用于描述薄殼結(jié)構(gòu)的大變形.該球殼的半徑為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.考慮到球殼的對稱性,現(xiàn)對其四分之一的模型進(jìn)行建模.本小節(jié)分別采用局部標(biāo)架GESF 殼單元與局部標(biāo)架縮減ANCF 殼單元離散四分之一開口球殼,并分析了24×24 與32×32 的兩種網(wǎng)格數(shù)目下球殼的靜力學(xué)變形.圖17 給出了不同單元數(shù)目下球殼上點A與點B的位移隨外載荷變化曲線.數(shù)值結(jié)果表明,上述兩類殼單元數(shù)值結(jié)果與參考文獻(xiàn)均吻合較好,本文提出的兩類單元能夠精確描述薄殼結(jié)構(gòu)的大變形特性.

圖16 頂部18?開口球殼初始構(gòu)形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

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

4.4 多柔體系統(tǒng)碰撞動力學(xué)仿真

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

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

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

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

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

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

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

圖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.在初始時刻,桿兩端點在慣性坐標(biāo)系下的位置坐標(biāo)分別為[0 ?0.2 1]Tm 與[1 0.2 0]Tm.空心盒不受重力,桿在沿Z 軸負(fù)方向的重力作用下開始運動.現(xiàn)采用SE(3)群剛體單元描述剛性桿與空心盒,并通過非光滑廣義α 錐互補(bǔ)方法對摩擦系數(shù)為0,0.2,0.4 以及0.6 的4 種情況進(jìn)行動力學(xué)仿真,算法參數(shù)譜半徑設(shè)為0.8,其中非光滑廣義a 錐互補(bǔ)方法可參見文獻(xiàn)[44].為了驗證該算法的正確性,同時采用罰函數(shù)方法對上述幾種情況進(jìn)行數(shù)值仿真.

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

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

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

4.5 大型桁架結(jié)構(gòu)展開動力學(xué)模擬

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

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

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

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

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

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

表1 迭代矩陣與約束方程Jacobian 矩陣更新次數(shù)及牛頓迭代總次數(shù)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 時的工況.當(dāng)NNI≥1 時,即每個時間步中牛頓迭代均更新迭代矩陣,共需進(jìn)行32 410 次牛頓迭代,平均每個時間步需要3.2 次牛頓迭代;當(dāng)NNI≥4 時,即若每個時間步中牛頓迭代達(dá)到四步時更新一次迭代矩陣,則迭代矩陣僅需更新19 次,約束方程雅可比矩陣需更新3083次,且總迭代次數(shù)較之前增加32%.而當(dāng)?shù)仃嚥桓聲r,總迭代次數(shù)比NNI≥4 時略有增大.由此可說明,基于SE(3)群局部標(biāo)架的建模方法能極大地降低剛體運動導(dǎo)致的幾何非線性,迭代矩陣可視為一個常數(shù)矩陣.

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

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

最后,基于上述建模方法,對具有8 個胞元的桁架結(jié)構(gòu)進(jìn)行展開動力學(xué)模擬.圖23 給出了不同時刻該桁架結(jié)構(gòu)的展開構(gòu)型.

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

4.6 基于區(qū)域分解的多柔體系統(tǒng)動力學(xué)通用并行算法

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

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

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

圖25 環(huán)形桁架?索網(wǎng)天線的區(qū)域分解示意圖Fig.25 Domain decomposition of the deployable hoop truss antenna

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

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

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

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

5 結(jié)束語

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

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

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

致謝

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