張騰 劉 鋮 ?,, 張志娟 劉紹奎

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

? (中物院高性能數(shù)值模擬軟件中心,北京 100088)** (北京應(yīng)用物理與計算數(shù)學(xué)研究所,北京 100088)

?? (北京空間飛行器總體設(shè)計部,北京 100094)

引言

多柔體系統(tǒng)的動力學(xué)過程通常包含柔性部件的大范圍剛體運動與彈性變形.柔性部件作大范圍剛體運動時引起的轉(zhuǎn)動與彈性變形引起的轉(zhuǎn)動是導(dǎo)致多柔體系統(tǒng)動力學(xué)過程呈現(xiàn)較高非線性的主要原因.如何有效降低多柔體系統(tǒng)動力學(xué)過程的非線性程度,特別是幾何非線性,是多柔體系統(tǒng)動力學(xué)領(lǐng)域面臨的挑戰(zhàn)之一.

近來,李群局部標(biāo)架(local frame of Lie group,LFLG)[1]思想與幾何精確方法(geometrically exact formulation,GEF)[2]的進(jìn)一步融合,有效規(guī)避了剛體轉(zhuǎn)動導(dǎo)致的幾何非線性.如Sonneville 等[3]、劉鋮和胡海巖[1]基于SE(3) 群(the special Euclidean group)局部標(biāo)架開發(fā)了幾何精確梁單元,數(shù)值結(jié)果表明任意剛體運動下單元廣義彈性力、廣義慣性力及他們的雅可比矩陣保持不變,從而消除了剛體運動帶來的幾何非線性.而且,在處理幾何非線性問題時,牛頓迭代過程中系統(tǒng)的雅可比矩陣不需要更新或者較少更新,這使動力學(xué)過程的計算效率得到明顯改善.最近,Zhang 等[4]基于SE(3)群局部標(biāo)架開發(fā)了含drilling 自由度的幾何精確板/殼單元,數(shù)值結(jié)果表明該單元繼承了李群局部標(biāo)架下幾何精確梁單元中的優(yōu)良特性,如幾何非線性程度降低、計算效率顯著改善等.事實上,李群局部標(biāo)架思想幾乎可以和所有考慮幾何非線性的梁、板/殼以及實體單元相結(jié)合,使其在計算中規(guī)避剛體運動導(dǎo)致的幾何非線性,如李群局部標(biāo)架下基于絕對節(jié)點坐標(biāo)方法(absolute nodal coordinate formulation,ANCF)[5]的三維全參數(shù)梁單元、李群局部標(biāo)架下基于絕對節(jié)點坐標(biāo)方法的三維全參數(shù)板/殼單元等,詳細(xì)內(nèi)容可參見作者之前的工作[1,6].

基于SE(3)群局部標(biāo)架的幾何精確單元雖擁有上述良好特性,但單元的實現(xiàn)過程相對困難,尤其是廣義彈性力及其雅可比矩陣的計算.其難點可以歸納為兩個方面,一是如何保證空間離散過程中插值算法的客觀性.Crisfiled 與Jeleni?[7]曾指出對轉(zhuǎn)動偽矢量直接進(jìn)行有限元插值將違背應(yīng)變的客觀性與材料路徑的無關(guān)性.為此,文獻(xiàn)[1,3]均在SE(3)群上空間離散時采用了相對插值,保證了離散應(yīng)變的客觀性.然而,采用相對插值進(jìn)一步導(dǎo)致廣義彈性力的計算及線性化過程變得復(fù)雜.如Sonneville[8]在其博士論文中通過Newton-Raphson 迭代獲取單元內(nèi)部相對插值的參考點,導(dǎo)致彈性力及其線性化過程更加復(fù)雜.Zhang 等[4]通過數(shù)值算例驗證了選取不同節(jié)點作為參考點對收斂結(jié)果的影響可忽略不計.因此,通常選取單元第一個節(jié)點作為相對插值的參考點,但這仍沒有規(guī)避相對插值帶來的復(fù)雜性.二是如何處理有限元離散與變分操作的先后順序.眾所周知,線性空間中不區(qū)分先離散再線性化(稱為離散一致線性化)與先線性化再離散(稱為連續(xù)一致線性化),但在SE(3)群中不然.張騰等曾按照這兩種方式嚴(yán)格推導(dǎo)了對應(yīng)版本的SE(3)群局部標(biāo)架下幾何精確板/殼單元,發(fā)現(xiàn)連續(xù)一致線性化版本的幾何精確板/殼單元在動力學(xué)仿真過程中遭遇收斂性問題.先離散再線性化版本的板/殼單元不存在收斂性問題,但該版本廣義彈性力的計算及線性化過程比較復(fù)雜.

基于連續(xù)體(continuum-based) 的殼[9](簡稱CB 殼)理論有望為規(guī)避上述滿足客觀性的插值及離散與變分操作順序帶來的復(fù)雜性提供新的解決思路.CB 殼單元的有限變形是建立在三維實體單元上,其廣義彈性力及雅可比矩陣是通過實體單元的廣義彈性力及雅可比矩陣以及實體單元上下表面與參考中面之間的節(jié)點轉(zhuǎn)換矩陣得到.顯然,實體單元的廣義彈性力及雅可比矩陣在計算時只包含平動自由度,無轉(zhuǎn)動自由度.平動屬于線性空間,在線性化過程中不存在有限元離散與變分操作先后順序的問題.而且,在線性空間中插值算法的客觀性可自然滿足,無需使用相對插值即可保證離散應(yīng)變的客觀性.因此,李群局部標(biāo)架思想與CB 殼理論的結(jié)合可以規(guī)避SE(3)群局部標(biāo)架下幾何精確板/殼單元中由插值和離散與變分操作所導(dǎo)致單元實現(xiàn)過程過于復(fù)雜的問題.

經(jīng)典CB 殼單元最初由Ahmad 等[10]于1970 年提出.Parisch[11]、Hughes和Liu[12-13]先后將CB 殼模型推廣至可用于非線性分析的格式.不同的是,Parisch 使用節(jié)點局部坐標(biāo)系與全局坐標(biāo)系之間的旋轉(zhuǎn)角描述截面轉(zhuǎn)動,Hughes和Liu 則使用沿厚度方向的矢量變化描述截面轉(zhuǎn)動.目前,CB 殼模型已經(jīng)應(yīng)用于多個領(lǐng)域.如莊茁和成斌斌[14]發(fā)展了基于CB 殼單元的擴展有限元,用于模擬三維任意擴展裂紋問題.Bergman和Yang[15]基于完全的拉格朗日方法采用CB 殼單元對形狀記憶聚合物建立有限元模型,成功應(yīng)用于航空航天領(lǐng)域.Momenan[16]將CB 殼單元用于軟生物組織的動力學(xué)分析.Han 等[17]采用CB 殼單元模擬流固耦合問題中經(jīng)歷大變形的柔性薄殼結(jié)構(gòu).

有限元方法中低階梁、板/殼和實體單元普遍存在薄膜(面內(nèi))、剪切和厚度閉鎖[18]問題,CB 殼單元也不例外.針對本文所提出基于SE(3)群局部標(biāo)架的CB 殼單元,分別采用Hellinger-Reissner 兩場變分原理[19-20]和假設(shè)自然應(yīng)變(assumed natural strain,ANS)方法[21]緩解面內(nèi)閉鎖和橫向剪切閉鎖.

為包含drilling 自由度,本文借鑒Fox和Simo[22]在幾何精確殼中引入運動約束的方法,并將該方法應(yīng)用于5 Dofs CB 殼單元,提出了含有drilling 自由度的CB 殼單元-6 Dofs CB 殼單元.本文引入drilling 自由度的目的有兩個,一是保證CB 殼單元可直接處理含約束的多體系統(tǒng)(如殼與梁、殼與殼連接等結(jié)構(gòu));二是6 Dofs CB 殼單元可完全消除剛體運動帶來的幾何非線性.Zhang 等[4]在最近的工作中通過數(shù)值算例驗證了基于SE(3)群局部標(biāo)架5 Dofs 幾何精確殼單元只能部分消除剛體運動帶來的幾何非線性,而6 Dofs 幾何精確殼單元可完全消除幾何非線性.類似地,本文所提出的6 Dofs CB 殼也可以完全消除剛體運動帶來的幾何非線性.

為解決上述SE(3)群局部標(biāo)架下幾何精確殼單元存在的兩個問題,本文融合李群局部標(biāo)架思想與經(jīng)典CB 殼理論推導(dǎo)基于SE(3) 群局部標(biāo)架的5 DoFs CB 殼單元.該單元可有效規(guī)避幾何精確殼單元中插值帶來的復(fù)雜性,單元離散應(yīng)變張量的客觀性可自然得到滿足.同時,該單元在線性化過程中不存在有限元離散與變分操作先后順序的問題,這進(jìn)一步簡化了單元廣義彈性力及其雅可比矩陣的計算.為提高單元收斂性,本文分別采用Hellinger-Reissner 兩場變分原理和ANS 方法緩解面內(nèi)和剪切閉鎖.為完全消除剛體運動導(dǎo)致的幾何非線性以及方便處理組合體結(jié)構(gòu),在5 Dofs CB 殼單元基礎(chǔ)上引入drilling 自由度,提出了6 Dofs CB 殼單元.最后,通過5 個靜力學(xué)算例驗證了所提出CB 殼單元的收斂精度,一個多體系統(tǒng)動力學(xué)算例驗證了CB 殼單元在計算過程中消除剛體運動帶來的幾何非線性,從而使系統(tǒng)的廣義彈性力、廣義慣性力及其雅可比矩陣滿足剛體運動的不變性.

1 基于SE(3)群局部標(biāo)架的CB 殼單元運動描述及本構(gòu)關(guān)系

1.1 運動描述

在經(jīng)典的CB 殼理論中,通常采用如下的假設(shè):(1)纖維保持直線(修正的Mindlin-Reissner 假設(shè));(2)平面應(yīng)力條件;(3)動量源于纖維的伸長,并且沿纖維方向忽略動量平衡[9].在上述假設(shè)基礎(chǔ)上,本文還假定板/殼單元厚度不變.

圖1 所示為8 節(jié)點CB 殼單元,白色和黑色節(jié)點分別表示從屬節(jié)點(slave nodes)和主控節(jié)點(master nodes).每個從屬節(jié)點有3 個平動自由度,每個主控節(jié)點有5 個獨立的自由度,包括3 個平動和2 個轉(zhuǎn)動.ξ=[ξ η ζ]T為母單元坐標(biāo)(自然坐標(biāo)),ζ對應(yīng)的每一個面稱為層,通常將ζ=0 看作單元參考中面.沿著ζ軸的線稱為纖維,沿著纖維方向的單位矢量稱為方向矢量.e1,e2和e3為慣性坐標(biāo)系下的單位正交基矢量.

圖1 8節(jié)點CB 殼單元Fig.1 Configuration of an eight-node CB shell element

從屬節(jié)點和主控節(jié)點之間的關(guān)系可表示為

式中,h0為殼單元厚度,xI+和xI-∈ R3為單元頂面和底面上從屬節(jié)點的位置矢量,xI∈ R3表示參考中面上節(jié)點的位置,I +和I-分別對應(yīng)殼單元上表面和下表面從屬節(jié)點編號,I為主控節(jié)點編號,其中,I +=1,2,···,4,I-=5,6,···,8,I=1,2,···,4.pI∈S2為主控節(jié)點I沿纖維方向的單位矢量,滿足pI=RIe3,S2為R3空間中的單位球[23],RI為當(dāng)前構(gòu)形中主控節(jié)點I的方向矢量相對于e3的旋轉(zhuǎn).以主控節(jié)點I為原點,通過RI可建立節(jié)點I局部隨體坐標(biāo)系與慣性坐標(biāo)系之間的關(guān)系,本文稱此局部隨體坐標(biāo)系為局部標(biāo)架[1].事實上,殼單元中面上任意一點處均存在局部標(biāo)架.在當(dāng)前構(gòu)形中,CB 殼單元上任意一點的位置矢量可表示為

式中,φ為當(dāng)前構(gòu)形中單元的廣義坐標(biāo),為標(biāo)準(zhǔn)等參三維線性形函數(shù)[24],I3為3 × 3 單位矩陣,NI*的表達(dá)式如下

將式(1) 作簡單變分操作可建立局部標(biāo)架下CB 殼單元從屬節(jié)點平動自由度與主控節(jié)點的平動和轉(zhuǎn)動自由度之間的變換關(guān)系,如下

SE(3)群局部標(biāo)架下CB 殼單元的插值本質(zhì)上是在線性空間中對實體單元的從屬節(jié)點進(jìn)行插值,獲得實體單元的廣義彈性力及其雅可比矩陣,之后通過單元層面的轉(zhuǎn)換矩陣TCB得到CB 殼單元的廣義彈性力及其雅可比矩陣.注意,TCB為4 個主控節(jié)點和8 個從屬節(jié)點之間的轉(zhuǎn)換矩陣,具體表達(dá)如下

1.2 本構(gòu)關(guān)系

由連續(xù)介質(zhì)力學(xué)[26]可得Green-Lagrange 應(yīng)變張量的表達(dá)式

F為變形梯度張量,定義如下

式中,J=?x/?ξ,J0=?X/?ξ,X表示初始構(gòu)形中任意一點的位置矢量.

將式(8)代入式(7),重寫Green-Lagrange 應(yīng)變張量

式中,為協(xié)變應(yīng)變張量(the covariant strain tensor),定義如下

兼顧對稱性和下文閉鎖處理的便捷性,引入E和對應(yīng)的Voigt 標(biāo)記

二者之間的轉(zhuǎn)換關(guān)系可表示為

式中,變換矩陣T的具體表示為

分別定義當(dāng)前和初始構(gòu)形中協(xié)變基矢量gi和Gi如下

將式(15)代入式(10),則協(xié)變應(yīng)變的具體表示為

式中,記 (·),i=?(·)/?ξi,i,j取1,2,3.

采用各向同性材料的線彈性本構(gòu)模型,本構(gòu)關(guān)系可表示為

式中,E為楊氏模量,v為泊松比.考慮到有限元建模對象為殼結(jié)構(gòu),引入CB 殼理論中第2 個假設(shè)-平面應(yīng)力條件,即 σzz=0,利用該關(guān)系式有

最終,單元本構(gòu)關(guān)系可表示為

1.3 SE(3)群統(tǒng)一插值

與上文介紹的CB 殼單元的插值過程不同,對于SE(3)群局部標(biāo)架下的幾何精確殼單元,本質(zhì)上是在SE(3) 群上對節(jié)點的相對運動張量H進(jìn)行插值[4],其中H包含了殼單元中面節(jié)點的位置矢量和旋轉(zhuǎn)矩陣,具體插值公式如下

式中,為單元任一點k和參考節(jié)點r之間的相對運動矢量,為單元節(jié)點數(shù),為一階Lagrange 插值形函數(shù),為形函數(shù)的元素,expSE(3)() 為SE(3)群中的指數(shù)映射,具體表達(dá)式如下

式中,Θ為旋轉(zhuǎn)偽矢量,u為在局部標(biāo)架中度量的位移矢量,和TSO(3)(Θ) 分別為特殊正交群SO(3)中的指數(shù)映射和切映射.關(guān)于SE(3)和SO(3)群的基礎(chǔ)知識參見文獻(xiàn)[8].關(guān)于SE(3)群中幾何精確殼更詳細(xì)的插值過程可參見作者之前的工作[4].

注意,SE(3)群局部標(biāo)架下幾何精確殼單元的插值過程是在流形上完成,需要插值相對運動張量才能保證應(yīng)變張量的客觀性,這使得廣義彈性力及其雅可比矩陣的計算比較復(fù)雜.而SE(3)群局部標(biāo)架下CB 殼單元的插值過程是在 R3空間中完成,插值過程比較簡單,大大簡化了廣義彈性力及其雅可比矩陣的計算,同時應(yīng)變張量的客觀性自然滿足.

2 基于SE(3)群局部標(biāo)架的CB 殼單元動力學(xué)平衡方程及其線性化

引入動力學(xué)平衡方程之前,先給出應(yīng)變的一次和二次變分的表達(dá)式,如下

CB 殼單元動力學(xué)平衡方程的弱形式為

式中,δWint,δWine,δWLa和δWext分別為廣義彈性力、慣性力、拉格朗日乘子項和外力所做的虛功.

首先,彈性力所做的虛功δWint由殼單元面內(nèi)部分和剪切部分組成,經(jīng)有限元離散后可表示為

εip=和為面內(nèi)應(yīng)變和橫向剪切應(yīng)變,和為對應(yīng)的協(xié)變應(yīng)變,他們之間的關(guān)系可通過轉(zhuǎn)換矩陣和得到[27],和來自于式(13)和式(14),具體表示為

由式(26)～ 式(28)可得到廣義彈性力的表達(dá)式

其中

對式(26)進(jìn)一步作二次變分(即線性化)可得

式(31)中DTCB與TCB類似,具體表示為

式中,和分別為彈性力Fint在節(jié)點I處的平動和轉(zhuǎn)動部分.

其次,慣性力所做虛功以及線性化過程與作者之前的工作[4]一致,此處不再贅述.

然后,拉格朗日乘子項所做虛功δWLa的具體表示為

式中L*=(G1)TR0RTg2-(G2)TR0RTg1,R0和R分別為初始和當(dāng)前構(gòu)形中任意一點的方向矢量相對于e3的旋轉(zhuǎn)矩陣.拉格朗日乘子系數(shù)γ*的選取參見文獻(xiàn)[4,28].注意,此處引入拉格朗日乘子項的目的是建立中面運動和drilling 自由度之間的關(guān)系,從而將5 Dofs CB 殼單元變成6 Dofs CB 殼單元.而且,其對應(yīng)的彈性力和雅可比矩陣是在主控節(jié)點上計算,不需要使用TCB進(jìn)行轉(zhuǎn)換.此處拉格朗日乘子項對應(yīng)彈性力FLa的離散形式及其線性化與作者之前的工作[4]一致,不再贅述.

最后,廣義外力所做虛功δWext的離散表達(dá)式為

式中,b為外界作用的體力,n為外力作用面的法向矢量,p為作用面力大小,ρ為材料密度,Γ0為作用面積.嚴(yán)格地說,牛頓迭代過程中外力的雅可比也有貢獻(xiàn).在動力學(xué)過程中,與廣義彈性力和廣義慣性力的切線剛度相比,廣義外力的雅可比矩陣通?？梢院雎圆挥?故本文不再給出具體形式.

3 基于SE(3)群局部標(biāo)架的CB 殼單元時間離散方法

通過上述空間離散,動力學(xué)平衡方程可表示為如下半離散形式的非線性微分代數(shù)方程組

其中,Φ和Φq分別為約束方程和約束方程的雅可比矩陣,λ為約束方程的拉格朗日乘子,t為時間,v和為系統(tǒng)的非物質(zhì)速度和非物質(zhì)加速度.

對于時間域離散,與幾何精確殼類似,采用李群廣義α 時間積分算法[1,29]求解上述系統(tǒng)動力學(xué)方程.李群廣義α 方法對前后兩個相鄰時刻的運動增量進(jìn)行差分離散,其具體離散格式可表示為

式中,αm，αf，β，γ和矢量a為李群廣義α 算法參數(shù),Δt為系統(tǒng)離散的時間步長,Ψn表示系統(tǒng)從n時刻到n+1 時刻在局部標(biāo)架下的運動增量,運動張量Hn∈SE(3)包含了n時刻殼單元中面作平動和轉(zhuǎn)動的全部信息.聯(lián)立動力學(xué)平衡方程(37)與式(38)并采用牛頓迭代可獲得下一時刻的速度vn+1和加速度,通過SE(3)群中指數(shù)映射可將運動張量Hn遞推到n+1 時刻Hn+1.

4 閉鎖緩解技術(shù)

4.1 基于Hellinger-Reissner 兩場變分原理的閉鎖緩解方法

本文采用Hellinger-Reissner 兩場變分原理處理殼單元面內(nèi)的閉鎖,兩場代表假設(shè)的位移場和假設(shè)的應(yīng)力場,且兩類場變量的插值相互獨立.假設(shè)殼單元面內(nèi)應(yīng)力的分布為

式中

殼單元面內(nèi)部分對應(yīng)的泛函可表示為

式中,為彈性系數(shù)矩陣Cip的逆,對單元面內(nèi)應(yīng)變部分的泛函求變分得

進(jìn)一步,對 δΠip作二次變分(即線性化)得

進(jìn)一步,離散方程組可以表示為

式中,Res為殘差,,Hip和Gip的表達(dá)式如下

式中βip=,靜力縮聚化為單場的形式

4.2 基于ANS 的剪切閉鎖緩解方法

本文使用假設(shè)自然應(yīng)變 (ANS) 法緩解CB 殼單元存在的剪切閉鎖.

下面構(gòu)造假設(shè)應(yīng)變場.如圖2 所示,ANS 的核心是沿ξ和η方向在等參單元四邊中點處布置采樣點,通過對采樣點處假設(shè)自然剪切應(yīng)變線性插值得到當(dāng)前構(gòu)型中假設(shè)自然剪切應(yīng)變,即

圖2 ANS 插值采樣點Fig.2 Interpolation sampling points for ANS

5 數(shù)值算例

本章節(jié)討論5 個常見的靜力學(xué)驗證算例和一個多體系統(tǒng)動力學(xué)算例.為考察SE(3)群中CB 殼單元的收斂精度,下文算例中將其他殼單元的收斂結(jié)果作為參考解.其中,GESF 表示作者最近的工作-SE(3)群局部標(biāo)架下經(jīng)過閉鎖處理的6 Dofs 一階幾何精確殼單元[4],ANCF 表示作者早期工作-基于絕對節(jié)點坐標(biāo)方法的縮減(slope deficiency)板/殼單元[30],Shell181 表示有限元軟件ANSYS 中的四邊形殼單元[31].本部分算例的數(shù)值結(jié)果均由基于SE(3)群局部標(biāo)架的6 Dofs CB 殼單元計算得到.值得一提的是,作者之前的工作[4]已經(jīng)對5 Dofs 與 6 Dofs 幾何精確殼單元在收斂性方面進(jìn)行測試,二者收斂精度類似.區(qū)別是,SE(3)群中6 Dofs 殼單元能夠完全消除剛體運動帶來的幾何非線性,而5 Dofs 只能部分消除.為避免重復(fù),本節(jié)不再提供5 Dofs CB 殼的數(shù)值結(jié)果.

5.1 懸臂板末端受均布彎矩

本節(jié)考察CB 殼在大變形大轉(zhuǎn)動方面的能力.

如圖3 所示的懸臂板,該模型經(jīng)常作為幾何非線性分析的測試算例[32],長、寬和厚分別為L=12 m，B=1 m和h=0.1 m.楊氏模量為E=6.825 ×107N/mm2,泊松比為v=0.0.在懸臂板最右端施加繞Y方向的均布彎矩Mmax=2πEI0/L,其中截面慣性矩I0=Bh3/12.由圖4可知,本文所提出SE(3)群局部標(biāo)架下CB 殼單元計算結(jié)果與解析解吻合較好.圖5給出了不同載荷下網(wǎng)格為20 × 2 的變形構(gòu)形圖.

圖3 受均布彎矩的懸臂板Fig.3 Cantilever plate under the distributed moments

圖4 懸臂板末端載荷-位移曲線Fig.4 Load-displacement curves at the end of the cantilever plate

圖5 不同載荷下變形構(gòu)形Fig.5 Deformed configurations under different loads

5.218 °開口球殼受集中載荷拉壓

本節(jié)考察CB 殼單元對不可展薄殼彎曲的模擬和緩解薄膜閉鎖的能力.

圖6 所示為帶有18°開口的球殼,該模型已成為殼單元測試的標(biāo)準(zhǔn)算例[33-34].球殼半徑R=10 mm,材料的楊氏模量E=6.825 × 107N/mm2,泊松比v=0.3,厚度分別為h=0.04 (R/h=250) mm和h=0.01(R/h=1000) mm,底部受大小相等且相互垂直的一對拉力和一對壓力作用,不同厚度對應(yīng)的載荷大小為F=200 N (h=0.04 mm)與F=5 N (h=0.01 mm).鑒于幾何對稱性,只需對其1/4 結(jié)構(gòu)建立有限元模型.AM和BN 為對稱平面,MN和AB 為自由邊.

圖618 °開口球殼Fig.6 Spherical shell with an 18° hole

由圖7和圖8 中點A和B的位移曲線可知,當(dāng)h=0.04 mm 時,12 × 12 個CB 殼單元得到收斂結(jié)果,當(dāng)h=0.01 mm 時,24 × 24 個CB 殼單元得到收斂結(jié)果,驗證了CB 殼單元結(jié)果的正確性.圖9 給出了在最大載荷下網(wǎng)格數(shù)目為16 × 16和24 × 24 的兩種厚度薄殼的變形構(gòu)形.

圖7 點A和B 位移曲線(h=0.04 mm)Fig.7 Displacement curves of points A and B (h=0.04 mm)

圖8 點A和B 位移曲線(h=0.01 mm)Fig.8 Displacement curves of points A and B (h=0.01 mm)

圖9 兩種厚度的殼對應(yīng)的變形構(gòu)形Fig.9 Deformed configurations corresponding to two thicknesses of shell

將ABAQUS 中厚度為0.04 mm、網(wǎng)格為24 ×24 的S4 R 殼單元點A和B的計算結(jié)果PA=3.41和PB=5.88 mm 作為參考解,統(tǒng)計CB 殼和GESF殼不同網(wǎng)格的計算結(jié)果如表1 中所示,(·%)表示與參考解的誤差.由此可得,與GESF 殼單元相比,CB 殼單元具有類似的收斂精度.同時,也驗證了閉鎖緩解技術(shù)有效改善了CB 殼單元的收斂精度.

表1 不同網(wǎng)格下兩種殼單元在點A和B 的位移結(jié)果對比(h=0.04 mm)Table 1 Comparison of displacement results of two shell elements at points A and B under different meshes(h=0.04 mm)

圖10 給出了采用16 × 16 網(wǎng)格時CB 殼單元和GESF 殼單元在不同載荷步(load step,LS)下牛頓迭代次數(shù)的分布.表2 中還給出了在配置處理器Intel Core i7-11700 (3.60 GHz)及32 GB 內(nèi)存的PC 機Matlab 編譯環(huán)境中,CB 殼與GESF 殼在收斂網(wǎng)格16 × 16 下彈性力及其Jacobian 矩陣的計算耗時和總牛頓迭代次數(shù)(number of iterations,NOI).盡管上述計算程序并未嚴(yán)格進(jìn)行計算效率優(yōu)化,但仿真耗時表明,由于CB 殼單元插值的簡潔性,其彈性力及雅可比矩陣的計算效率遠(yuǎn)高于同階數(shù)的GESF殼.對于0.01 mm 厚度的薄殼,可得出與0.04 mm 厚度的殼在收斂精度、牛頓迭代次數(shù)和計算時間等方面類似的結(jié)論.

圖10 不同載荷步下迭代次數(shù)(h=0.04 mm,16 × 16)Fig.10 Number of iterations at different load steps(h=0.04 mm,16 × 16)

表2 不同載荷步下彈性力及其雅可比矩陣計算時間和迭代次數(shù)對比(h=0.04 mm,16 × 16)Table 2 Comparison of computational time of elastic forces and their Jacobian matrices and number of iterations under different load steps (h=0.04 mm,16 × 16)

5.3 圓柱殼中心集中受拉

本節(jié)考察CB 殼單元模擬殼結(jié)構(gòu)發(fā)生以大薄膜應(yīng)變主導(dǎo)彈性變形的能力.

如圖11 所示,圓柱殼長、內(nèi)徑和厚度分別為L=10.35 mm,R=4.953 mm和h=0.094 mm.楊氏模量E=1.05 × 107N/mm2,泊松比v=0.3125.圓柱殼中心受集中載荷大小F=40 kN.鑒于幾何對稱性,只需對其1/8 結(jié)構(gòu)建立有限元模型.由于該圓柱殼在拉伸過程中發(fā)生了屈曲,本文采用Lam和Morley 提出的弧長算法[35]求解平衡方程,實現(xiàn)平衡路徑的自動追蹤.由圖12可得,使用24 × 16 單元網(wǎng)格時,CB 殼單元計算結(jié)果已經(jīng)收斂,且收斂結(jié)果與GESF(24 × 16) 結(jié)果一致.可知,同等網(wǎng)格下CB 殼單元收斂精度與GESF 類似.圖13 繪制了不同視角下圓柱殼最后的變形構(gòu)形.值得一提的是,針對此算例,Zhang 等[4]曾使用SE(3)群中幾何精確殼單元對此進(jìn)行計算.若忽略薄膜應(yīng)變中的二階項,則直接導(dǎo)致點B,C的位移曲線與收斂結(jié)果之間存在較大偏差,保留此項后方可得出正確結(jié)果.因此,對于發(fā)生以大薄膜應(yīng)變主導(dǎo)彈性變形的殼結(jié)構(gòu),使用本文提出的CB 殼單元仍然可以得出正確的數(shù)值結(jié)果.

圖11 圓柱殼初始構(gòu)形Fig.11 Initial configuration of the pinched cylinder

圖12 點A,B和C 載荷位移曲線Fig.12 Load-displacement curves of points A,B and C

圖13 三個不同視角下圓柱殼的變形構(gòu)形Fig.13 Deformed configurations for cylindrical shell from three different perspectives

5.4 中心受擠壓圓柱形屋頂?shù)暮笄治?/h3>

本節(jié)考慮圖14 所示中心受擠壓的圓柱形屋頂,該算例被廣泛用于檢驗單元模擬殼結(jié)構(gòu)發(fā)生屈曲變形的正確性[19,36-38].

圖14 中心受擠壓的鉸接圓柱形屋頂Fig.14 Hinged cylindrical roof subjected to a central pinching force

圓柱形屋頂邊界條件為直邊鉸接、曲邊自由,中心點C在負(fù)Z方向受擠壓載荷Fmax=3000 N.楊氏模量E=3102.75 N/mm2,泊松比為v=0.3.該模型長L=254 mm,半徑R=2540 mm,對應(yīng)的圓心角θ=0.1 rad.鑒于幾何性對稱性,只需對結(jié)構(gòu)的1/4 建立有限元模型.考慮兩種厚度h=12.7 mm和h=6.35 mm 對結(jié)構(gòu)發(fā)生屈曲變形的影響.隨著載荷的增加,圓柱形屋頂在變形過程中出現(xiàn)突然的跳躍(snap-through)和曲率的較大變化.采用弧長法[35]計算得到點C在Z方向的位移變化如圖15 所示,曲線Ⅰ和Ⅱ分別對應(yīng)厚度為12.7 mm和6.35 mm 時C點的位移變化.與曲線Ⅰ相比,曲線Ⅱ中屈曲載荷相對較小,故厚度變化對結(jié)構(gòu)屈曲變形有較大影響.CB 殼單元的計算結(jié)果與Sze 等[36]中使用ABAQUS中S4 R 殼單元計算結(jié)果非常一致,這說明CB 殼單元在處理發(fā)生后屈曲結(jié)構(gòu)的問題時,同樣表現(xiàn)良好.

圖15 圓柱形屋頂中心的載荷位移曲線Fig.15 Load-displacement curves at the center of the cylindrical roof

5.5 直角懸臂板受均布力/力矩

本節(jié)考察CB 殼單元在處理組合體結(jié)構(gòu)的數(shù)值表現(xiàn),Simo 等[39]和Li 等[40]曾使用該模型檢驗殼單元的精度.

直角懸臂板的幾何及材料參數(shù)如圖16 中所示.該板右端固定,左端分別受均布力和均布力矩兩種工況.λ0表示施加的均布力或均布力矩的大小.圖17和圖18 給出了直角板左端節(jié)點的X方向和Z方向位移變化.可知,CB 殼單元收斂結(jié)果與shell181,GESF 殼單元收斂結(jié)果吻合良好.圖16 還給出了均布力和均布力矩兩種工況下直角懸臂板變形構(gòu)形.

圖16 直角懸臂板幾何、材料參數(shù)及變形構(gòu)形Fig.16 Geometry,material parameters and deformed configurations of right-angle cantilever plate

圖17 均布力作用下直角板左端載荷-位移曲線Fig.17 Load-displacement curves of the left end of the right-angle plate subjected to uniform forces

圖18 均布力矩作用下直角板左端載荷-位移曲線Fig.18 Load-displacement curves of the left end of the right-angle plate subjected to uniform moments

5.6 空間薄圓柱殼雙擺多體系統(tǒng)動力學(xué)

本節(jié)通過一個空間薄圓柱殼雙擺模型的多體系統(tǒng)動力學(xué)仿真考察CB 殼單元在動力學(xué)過程中描述大變形和剛體位移以及消除幾何非線性方面的表現(xiàn),Liu 等[30]、Shi 等[41]曾研究該模型在重力作用下的動力學(xué)響應(yīng).

如圖19 所示,曲面ACEFDB 由平面I和II 將兩個半徑均為0.3 m 的圓柱殼裁剪形成,裁剪平面I和II 與平面XOZ之間的夾角滿足關(guān)系tanθ=0.5.圓柱殼Y方向長度為1.2 m,厚度為0.01 m,楊氏模量為E=1.0 × 109Pa,材料密度為7810 kg/m3,泊松比為0.3.圓柱殼I 左端AB 簡支,圓柱殼I和II 通過沿Y軸方向的CD 邊圓柱鉸接.重力加速度為9.81 m/s2,總仿真時間為1.5 s.本算例中,牛頓迭代收斂標(biāo)準(zhǔn)設(shè)為廣義坐標(biāo)修正量‖Δ‖2＜1.0 × 10-6.

圖19 由兩個圓柱殼組成的空間雙擺模型Fig.19 A double pendulum composed of two parts of cylindrical shells

首先,通過與作者早期工作中基于ANCF 的殼單元[30]對比,驗證SE(3)群局部標(biāo)架下CB 殼單元描述動力學(xué)問題的正確性.圖20 為不同單元數(shù)量下點F在Z方向的位移變化.當(dāng)每個殼單元網(wǎng)格數(shù)目為15 × 15 時,該模型計算結(jié)果與ANCF 殼單元數(shù)值結(jié)果吻合良好,驗證了動力學(xué)仿真過程中CB 殼單元的收斂精度.若以圖19 中圓弧AC和CE 的長度變化分別作為兩個圓柱殼發(fā)生彈性變形的度量,變形率即為當(dāng)前和初始構(gòu)形中兩弧長之差與初始構(gòu)形中弧長的比值,變形率為負(fù)數(shù)時表示當(dāng)前構(gòu)形中的弧長小于初始構(gòu)形的弧長.在1.5 s 仿真過程中,圓柱殼I和II 的最大變形率可達(dá)58%和30%.圖21給出了網(wǎng)格數(shù)目為25 × 16 下雙擺圓柱殼在不同時刻的構(gòu)形圖,可以看出該模型發(fā)生了較大的變形和剛體位移.

圖20 不同單元數(shù)目下點F 在Z 方向的位移曲線Fig.20 Displacement curve of point F in Z-direction with different number of elements

圖21 雙擺變形構(gòu)形Fig.21 Deformed configurations of the double pendulum

其次,通過分析系統(tǒng)質(zhì)量矩陣和切剛度矩陣更新情況驗證SE(3)群中CB 殼單元能夠有效降低剛體運動帶來的幾何非線性.若忽略廣義外力的雅可比矩陣,系統(tǒng)迭代矩陣的具體表示為

其中,M和K分別為系統(tǒng)的廣義質(zhì)量矩陣和切剛度矩陣.表3 記錄了總仿真時間的前1 s 雙擺模型的迭代矩陣主要部分Jq=Δt2β(M+K) 的更新次數(shù)(number of update,NOU)和圓柱殼在雙擺過程中牛頓迭代總次數(shù)NOI.“NNI≥i”表示在一個時間步內(nèi),若牛頓迭代的次數(shù)大于等i,則矩陣Jq強制更新.“Frozen”表示系統(tǒng)從初始構(gòu)形計算得到矩陣Jq之后不再更新.

表3 迭代矩陣更新次數(shù)及牛頓迭代總次數(shù)(E=1.0 GPa)Table 3 The number of times that the iteration matrix is updated and the number of total Newton iterations (E=1.0 GPa)

表3 的數(shù)值結(jié)果表明:以時間步長為5.0 ×10-4s 工況為例,當(dāng)Jq矩陣在每個牛頓迭代步更新時,需要更新的次數(shù)為7022.當(dāng)NNI≥3 時(即每個時間步長中牛頓迭代次數(shù)等于或超過3 步時更新一次Jq矩陣)總牛頓迭代次數(shù)增加了5.5%,Jq矩陣的更新次數(shù)下降了78.9%.圖22 給出了系統(tǒng)Jq矩陣更新次數(shù)的分布以及圓柱殼I和II 的變形率變化曲線,可知在0.357 s 時Jq發(fā)生了第一次更新,這段時間主要以圓柱殼I 的變形為主,此時圓柱殼I和II 的變形率約為33%和0.4%.0.357 s 至0.408 s 期間Jq不需要更新.直至0.409 s 發(fā)生了第二次更新,此時圓柱殼I和II 的變形率約為52%和2%.從0 s 至0.409 s 期間,雙擺模型變形率的變化基本保持平穩(wěn).0.409 s 之后,Jq在每一個時間步更新的次數(shù)和頻率有所增加,這是由于圓柱殼I和II 在發(fā)生相對較大變形的同時變形率出現(xiàn)了較大的波動所致.事實上,在仿真計算中,每個物體甚至每個單元可以根據(jù)當(dāng)前的變形程度自適應(yīng)的選擇是否更新單元的廣義質(zhì)量矩陣和切剛度矩陣,其魯棒性的判斷標(biāo)準(zhǔn)值得進(jìn)一步研究.

圖22 Jq 更新次數(shù)分布及圓柱殼的變形率曲線Fig.22 Distribution of the number of times that Jq is updated and curves of deformation rate of cylindrical shells

若考慮較小時間步長的工況,此時系統(tǒng)的迭代矩陣中廣義質(zhì)量矩陣的占比提升.通常,柔體的廣義質(zhì)量矩陣與系統(tǒng)的彈性變形相關(guān),而切剛度矩陣與系統(tǒng)彈性變形的梯度相關(guān).與切剛度矩陣相比,廣義質(zhì)量矩陣變化較小.因此,時間步長為1.0 × 10-4s時,總牛頓迭代次數(shù)僅增加了2.4%,Jq矩陣更新次數(shù)下降了98.8%.與總迭代次數(shù)相比,Jq矩陣的更新次數(shù)基本上可以忽略不計.

若將上述模型中楊氏模量設(shè)置為1.0 × 1011Pa,其余參數(shù)保持不變,與表3 類似,表4 記錄了該楊氏模量下Jq的更新次數(shù)和牛頓迭代總次數(shù).考察時間步長為1.0 × 10-3s 時的工況,當(dāng)Jq矩陣在每個牛頓迭代步更新時,需要更新的次數(shù)為2972.當(dāng)系統(tǒng)出現(xiàn)Frozen 現(xiàn)象時,總牛頓迭代次數(shù)增加了8.8%,Jq矩陣不需要更新.

表4 迭代矩陣更新次數(shù)及牛頓迭代總次數(shù)(E=100 GPa)Table 4 The number of times that the iteration matrix is updated and the number of total Newton iterations (E=100 GPa)

由上述數(shù)值結(jié)果可知,基于SE(3)群局部標(biāo)架的CB 殼單元繼承了李群局部標(biāo)架方法體系下非線性有限單元(如基于SE(3)群局部標(biāo)架的幾何精確梁單元、基于SE(3) 群局部標(biāo)架的幾何精確殼單元)的優(yōu)勢[1],即在計算中消除了剛體運動帶來的幾何非線性,使動力學(xué)系統(tǒng)的廣義慣性力、廣義彈性力及其雅可比矩陣滿足剛體運動的不變性.

值得一提的是,在多體系統(tǒng)求解過程中,除減小質(zhì)量矩陣及切線剛度矩陣更新次數(shù)外,如何有效利用上述性質(zhì)提高線性方程組的直接解法或迭代算法的計算效率仍值得進(jìn)一步探索.

6 結(jié)論

本文融合李群局部標(biāo)架方法和經(jīng)典CB 殼理論開發(fā)了基于SE(3)群局部標(biāo)架的CB 殼單元.從單元計算角度看,廣義彈性力及其雅可比矩陣的線性化過程不需要復(fù)雜的插值運算和變分操作,大大簡化了廣義彈性力及其雅可比矩陣的計算過程,而且不需要使用相對插值即可保證插值算法的客觀性和路徑無關(guān)性.從幾何非線性角度看,該單元繼承了李群局部標(biāo)架方法體系中非線性有限單元(如SE(3)群局部標(biāo)架下的幾何精確梁單元、幾何精確板/殼單元) 的優(yōu)勢,即消除了由剛體運動帶來的幾何非線性,廣義慣性力、廣義彈性力及其雅可比矩陣滿足剛體運動的不變性.從收斂性角度看,閉鎖緩解技術(shù)有效改善了SE(3)群局部標(biāo)架下CB 殼單元的收斂精度,使得CB 殼單元與GESF 殼單元在收斂精度方面類似,具有漸進(jìn)二階收斂率.從計算時間角度看,與GESF 殼單元相比,CB 殼單元在計算效率方面具有較大提升.此外,對于組合體結(jié)構(gòu),例如殼與殼連接結(jié)構(gòu),含有drilling 自由度的CB 殼單元可直接進(jìn)行處理.

借鑒上述基于SE(3)群局部標(biāo)架的CB 殼單元的研究思路,后續(xù)可融合李群局部標(biāo)架方法與基于連續(xù)體的梁理論,展開基于SE(3) 群局部標(biāo)架的CB 梁單元的相關(guān)研究工作,實現(xiàn)對SE(3)群中幾何精確梁單元插值過程的簡化.進(jìn)一步,可在上述研究基礎(chǔ)上,結(jié)合共旋坐標(biāo)法(Co-rotational formulation,CRF)[42-45]的建模思想,提出基于SE(3)群的共旋CB 梁單元和CB 殼單元,可進(jìn)一步提高計算效率.