劉鋮 胡海巖

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

引言

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

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

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

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

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

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

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

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

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

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

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

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

1 基于SE(3)群局部標(biāo)架的幾類(lèi)梁?jiǎn)卧?/h2>

本節(jié)以含轉(zhuǎn)動(dòng)參數(shù)的幾何精確梁與幾種無(wú)轉(zhuǎn)動(dòng)參數(shù)的梁?jiǎn)卧獮槔?闡述基于SE(3)群局部標(biāo)架的梁?jiǎn)卧７椒?

1.1 幾何精確Timoshenko 梁?jiǎn)卧?/h3>

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

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

圖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)型時(shí)的中線沿著慣性坐標(biāo)系的X軸方向,采用R3中的矢量φ0描述梁的中線位置.以矢量φ0的端點(diǎn)為原點(diǎn),建立梁的局部連體坐標(biāo)系,其方向與慣性坐標(biāo)系一致,記其正交基矢量為此時(shí),局部連體坐標(biāo)系與慣性坐標(biāo)系之間的旋轉(zhuǎn)矩陣可記為R0=I3.

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

基于上述局部標(biāo)架和SE(3)群理論,梁?jiǎn)卧闹芯€位置矢量及梁的橫截面旋轉(zhuǎn)矩陣合并為運(yùn)動(dòng)張量H,可表示為

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

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

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

進(jìn)一步,在初始構(gòu)型與當(dāng)前構(gòu)型中,梁?jiǎn)卧先我恻c(diǎn)的位置矢量可分別表示為

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

對(duì)于動(dòng)力學(xué)問(wèn)題,梁?jiǎn)卧谌我鈺r(shí)刻的中線速度和截面轉(zhuǎn)動(dòng)角速度,可通過(guò)運(yùn)動(dòng)學(xué)方程,也被稱(chēng)為Poisson 方程,表示為

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

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

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

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

其中,axial()表示反對(duì)稱(chēng)矩陣對(duì)應(yīng)的軸矢量,Γ表示梁中線的拉伸應(yīng)變和橫向剪切應(yīng)變,κ 的分量κ1表示梁的橫截面扭轉(zhuǎn)應(yīng)變,分量κ2與κ3分別表示橫截面繞i2與i3軸的彎曲應(yīng)變.對(duì)于線彈性材料,通過(guò)平面應(yīng)力假設(shè)以及截面積分,可將梁?jiǎn)卧獌?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)造無(wú)約束梁?jiǎn)卧膭?dòng)力學(xué)方程.在此基礎(chǔ)上,通過(guò)對(duì)梁?jiǎn)卧膭?dòng)力學(xué)方程進(jìn)行空間域有限元離散以及時(shí)間域差分離散,可得到一個(gè)非線性代數(shù)方程組,進(jìn)而進(jìn)行數(shù)值求解.

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

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

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

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

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

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

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

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

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

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

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

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

1.2 無(wú)轉(zhuǎn)動(dòng)參數(shù)的梁?jiǎn)卧?/h3>

對(duì)于無(wú)轉(zhuǎn)動(dòng)參數(shù)的梁?jiǎn)卧?也可建立基于SE(3)群局部標(biāo)架,進(jìn)而消除梁?jiǎn)卧獎(jiǎng)傮w運(yùn)動(dòng)導(dǎo)致的幾何非線性.

首先,考察基于ANCF 描述的全參數(shù)梁?jiǎn)卧?采用中線位置矢量以及沿中線方向的斜率矢量作為廣義坐標(biāo),則梁?jiǎn)卧芯€上任意一點(diǎn)的位置矢量rmid可通過(guò)三階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為梁?jiǎn)卧某跏奸L(zhǎng)度,I3為三階單位矩陣.此外,通過(guò)一階線性Lagrange 插值得到梁截面的斜率矢量r,Y與r,Z.因此,該梁?jiǎn)卧先我恻c(diǎn)的位置矢量可表示為

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

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

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

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

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

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

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

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

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

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

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

Simo 等[15]首先提出了R3×S2幾何精確板殼單元,隨后,并將該單元推廣至可考慮變厚度板殼單元[16]、以及考慮彈塑性本構(gòu)的板殼單元[17].本文進(jìn)一步將上述單元推廣至SE(3)群局部標(biāo)架中,消除剛體運(yùn)動(dòng)導(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)型中,板單元上任意點(diǎn)的位置矢量可分別表示為

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

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

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

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

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

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

2.2 幾類(lèi)無(wú)轉(zhuǎn)動(dòng)參數(shù)的板殼單元

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

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

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

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

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

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

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

具體來(lái)說(shuō),上述多柔體系統(tǒng)動(dòng)力學(xué)建模和計(jì)算方法可分解為如下6 個(gè)方面.

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

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

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

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

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

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

4 數(shù)值算例

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

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

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

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

(1) 考察含線性約束的雙擺模型的小變形動(dòng)力學(xué)問(wèn)題.如圖3(a) 所示,雙擺的兩根桿分別長(zhǎng)0.3 m和0.5 m,其寬度與高度均為0.005 m,材料參數(shù)為:ρ=2700 kg/m3,E=7.0 × 1010Pa.桿I 與地面以及兩桿間均通過(guò)球鉸連接,系統(tǒng)約束方程均為線性約束.在初始時(shí)刻,桿I 靜止放置于X軸上,桿II 繞B點(diǎn)存在初始角速度[3 0 ?4]Trad/s.系統(tǒng)不受外力作用,在初始速度作用下開(kāi)始運(yùn)動(dòng).通過(guò)基于SE(3)群局部標(biāo)架的幾何精確梁?jiǎn)卧獙?duì)雙擺系統(tǒng)進(jìn)行離散,采用保能量?(角)動(dòng)量時(shí)間積分算法求解系統(tǒng)動(dòng)力學(xué)方程,仿真時(shí)間為100 s,時(shí)間步長(zhǎng)為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 含線性約束雙擺模型在小變形運(yùn)動(dòng)中末端C 點(diǎn)位置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點(diǎn)的位置隨時(shí)間變化規(guī)律,其余70 s 內(nèi)位置變化趨勢(shì)大體類(lèi)似.由圖可見(jiàn),小變形雙擺模型長(zhǎng)時(shí)間歷程動(dòng)力學(xué)響應(yīng)非常復(fù)雜,運(yùn)動(dòng)非線性程度較高.由于強(qiáng)非線性系統(tǒng)對(duì)初值異常敏感,難以評(píng)價(jià)算法精度,而系統(tǒng)守恒量能夠作為評(píng)價(jià)算法精度的一個(gè)重要指標(biāo).由于雙擺系統(tǒng)僅在與地面連接處受到約束力作用,該系統(tǒng)對(duì)慣性坐標(biāo)系原點(diǎn)的角動(dòng)量以及系統(tǒng)總能量守恒.圖5 與圖6 分別給出了系統(tǒng)能量、系統(tǒng)沿三個(gè)方向角動(dòng)量的演化規(guī)律,其中曲線上下方的數(shù)字表示該物理量波動(dòng)的最大值與最小值.本質(zhì)上來(lái)說(shuō),含線性約束多體系統(tǒng)在計(jì)算過(guò)程中可等效為一個(gè)無(wú)約束系統(tǒng),半離散平衡方程為一常微分方程組.對(duì)于此類(lèi)情況,時(shí)間積分算法可精確地保持系統(tǒng)能量及角動(dòng)量,系統(tǒng)能量的相對(duì)誤差僅為5.0 × 10?12,系統(tǒng)角動(dòng)量的相對(duì)誤差為7.0 × 10?10.

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

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

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

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

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

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

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

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

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

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

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

4.2 廣義a 方法后驗(yàn)誤差估計(jì)方法

本小節(jié)通過(guò)一個(gè)自由剛體運(yùn)動(dòng)問(wèn)題和剛性雙擺運(yùn)動(dòng)問(wèn)題,驗(yàn)證基于SE(3)群描述的廣義a 方法誤差估計(jì)方法的有效性.

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

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

圖11 基于SE(3)群描述的廣義a 方法計(jì)算單剛體運(yùn)動(dòng)的后驗(yàn)誤差估計(jì)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ù)方向的重力作用下開(kāi)始運(yùn)動(dòng).采用基于SE(3)群描述的廣義a 方法進(jìn)行仿真計(jì)算,總仿真時(shí)間為1 s,時(shí)間步長(zhǎng)為1.0 × 10?3s,算法參數(shù)譜半徑為0.8.圖12 與圖13 分別給出了桿OB 的估計(jì)誤差與參考誤差的二范數(shù).由于后驗(yàn)誤差計(jì)算方法是基于常微分方程推導(dǎo),而多柔體系統(tǒng)約束的引入,導(dǎo)致估計(jì)誤差與參考誤差之間存在一定的差異,但兩者的誤差量級(jí)與整體趨勢(shì)一致,通過(guò)估計(jì)誤差進(jìn)行變步長(zhǎng)時(shí)間積分算法設(shè)計(jì)是可行的.

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

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

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

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

如圖14 所示,一個(gè)無(wú)約束的圓柱殼在A,E兩點(diǎn)受到一對(duì)大小相等、方向相反的集中力F=40 kN的作用,其中A點(diǎn)為圓柱體上母線的中點(diǎn),E點(diǎn)為下母線的中點(diǎn).圓柱殼的長(zhǎng)度L為10.35 mm,內(nèi)徑R與厚度分別為4.953 mm 與0.094 mm,材料楊氏模量E=10.5 × 106N/mm2,泊松比υ=0.312 5.考慮到圓柱殼模型的對(duì)稱(chēng)性,只需要對(duì)該模型上半部分的四分之一(ABCD)進(jìn)行建模.現(xiàn)將在圓柱殼上施加外力F 的過(guò)程平均分成100 個(gè)增量步,并設(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三點(diǎn)的位移隨外力增加的變化曲線.當(dāng)載荷接近20 kN 時(shí),圓柱殼發(fā)生了屈曲變形,屈曲前圓柱殼以彎曲變形為主,屈曲后圓柱殼以拉伸變形為主.因此,本算例中GEF 殼單元不能對(duì)忽略薄膜應(yīng)變對(duì)彎曲變形的影響.數(shù)值結(jié)果表明,上述兩類(lèi)殼單元數(shù)值結(jié)果與參考文獻(xiàn)均吻合較好,均能夠精確地描述圓柱殼的大變形特性.關(guān)于一般性彈性體的屈曲、后屈曲及分叉模擬算法可參見(jiàn)作者前期工作[42].

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

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

圖16 頂部18?開(kāi)口球殼初始構(gòu)形Fig.16 Initial configuration of a hemispherical shell with 18?hole

圖17 開(kāi)口球殼上A、B 點(diǎn)的位移大小Fig.17 Magnitudes of displacements at nodes A and B of a hemispherical shell

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

5 結(jié)束語(yǔ)

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

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

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

致謝

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