基于T樣條的變網(wǎng)格等幾何薄板動(dòng)力學(xué)分析1)

王 悅 崔雅琦 於祖慶 蘭 朋,2) 陸念力

* (哈爾濱工業(yè)大學(xué)機(jī)電工程學(xué)院,哈爾濱 150001)

? (河海大學(xué)機(jī)電工程學(xué)院,江蘇常州 213002)

引言

隨著柔性輕質(zhì)結(jié)構(gòu)在航天、汽車(chē)以及船舶工程等領(lǐng)域日益廣泛的應(yīng)用,柔性多體系統(tǒng)的動(dòng)力學(xué)分析也愈加重要.然而,以線彈性小變形假設(shè)為基礎(chǔ)的,將節(jié)點(diǎn)位移和轉(zhuǎn)動(dòng)作為坐標(biāo)的傳統(tǒng)有限單元不適合解決存在大位移、大變形的柔性多體系統(tǒng)的動(dòng)力學(xué)問(wèn)題.此外,對(duì)于傳統(tǒng)結(jié)構(gòu)有限元,將幾何數(shù)據(jù)轉(zhuǎn)換為有限元網(wǎng)格數(shù)據(jù)是困難且耗時(shí)的[1].據(jù)估計(jì),在航空航天、船舶制造和汽車(chē)工業(yè)中,大約80%的分析時(shí)間用于網(wǎng)格生成[2].在計(jì)算機(jī)輔助工程(CAE)領(lǐng)域,有限單元只是計(jì)算機(jī)輔助設(shè)計(jì)(CAD)幾何模型的近似表示而非精確描述.這種CAE 網(wǎng)格與CAD 幾何之間的區(qū)別將會(huì)導(dǎo)致精確性問(wèn)題,例如在殼分析中,隨著幾何缺陷的增加,屈曲荷載有相當(dāng)大的降低[2].

為了準(zhǔn)確刻畫(huà)柔性體大位移、大變形的動(dòng)力學(xué)行為,研究者將彈性力學(xué)與有限元方法和多體動(dòng)力學(xué)理論相結(jié)合,形成了絕對(duì)節(jié)點(diǎn)坐標(biāo)方法(ANCF).ANCF 利用梯度代替轉(zhuǎn)角對(duì)柔性體位移場(chǎng)進(jìn)行插值并與多體系統(tǒng)動(dòng)力學(xué)理論相結(jié)合,使得其具備描述柔性體大變形、大轉(zhuǎn)動(dòng)的能力,從而成為多柔體系統(tǒng)動(dòng)力學(xué)領(lǐng)域的研究熱點(diǎn)[3].此外,研究者又發(fā)現(xiàn)其與計(jì)算機(jī)輔助分析軟件中的非均勻有理B 樣條(NURBS)幾何體家族之間存在線性轉(zhuǎn)化關(guān)系,從而可以實(shí)現(xiàn)從幾何模型到ANCF 單元網(wǎng)格之間的快速轉(zhuǎn)換,避免了網(wǎng)格劃分的步驟.相關(guān)的研究始于纜索單元與B 樣條曲線之間的轉(zhuǎn)化關(guān)系[4],并延伸至有理的ANCF 纜索單元與NURBS 曲線間關(guān)系[5].相關(guān)研究深刻揭示了二者在幾何體形狀描述、連續(xù)性控制等方面的高度相似性[6-7].此后,學(xué)者們將研究拓展至ANCF 雙參數(shù)單元[8-9].研究結(jié)果表明,當(dāng)單元維數(shù)升高,NURBS 幾何體張量積形式的參數(shù)方程會(huì)引入冗余自由度,從而需要對(duì)ANCF 單元進(jìn)行改造從而使其與之匹配.例如增加混合梯度向量,或在板單元中間添加額外節(jié)點(diǎn)等.這一問(wèn)題在三維單元中更加嚴(yán)重.Yu 和Cui[10]提出了一種48 自由度的八節(jié)點(diǎn)實(shí)體梁?jiǎn)卧?可以精確離散使用相同的基函數(shù)階次Bézier 體所代表的幾何模型.Ma 等[11]研究了一種基于三次有理Bézier 體積的三維有理絕對(duì)節(jié)點(diǎn)坐標(biāo)公式(RANCF)流體單元,準(zhǔn)確描述初始曲線形狀的液柱.

對(duì)計(jì)算機(jī)輔助工程和計(jì)算機(jī)輔助設(shè)計(jì)進(jìn)行整合的另外一個(gè)方向就是由Hughes 等[2]于2005年提出的等幾何分析(IGA),其概念是將CAD 和有限元分析(FEA)兩個(gè)領(lǐng)域統(tǒng)一起來(lái),使用相同的基礎(chǔ)進(jìn)行幾何描述和分析[12].IGA 方法因?yàn)榫哂袑?shí)現(xiàn)無(wú)縫整合設(shè)計(jì)和分析的潛力[13]而受到了力學(xué)領(lǐng)域的關(guān)注[2].學(xué)者們進(jìn)行了很多針對(duì)殼體和板的等幾何分析[14-17],這些研究多是基于當(dāng)前CAD 系統(tǒng)的行業(yè)標(biāo)準(zhǔn)NURBS.然而,由于NURBS 控制點(diǎn)呈矩形網(wǎng)格分布,從而導(dǎo)致其拓?fù)浣Y(jié)構(gòu)中存在大量冗余控制點(diǎn)[18],進(jìn)而導(dǎo)致分析效率降低.此外,修剪后的NURBS 表面不可避免地存在間隙和重疊是影響CAD,CAM和CAA 系統(tǒng)互操作性的最嚴(yán)重的障礙之一[19].為了克服NURBS 的這些缺點(diǎn),Sederberg 等[20]提出了T 樣條.與NURBS 曲面相比,T 樣條曲面的一個(gè)優(yōu)點(diǎn)是允許局部細(xì)化.由于T 交叉點(diǎn)(T-junction)的存在使得T 樣條允許控制點(diǎn)局部插入到控制網(wǎng)格中,而非整行或整列地增加控制點(diǎn)[20],從而大大減少了多余控制點(diǎn)的數(shù)量[18].T 樣條的另一個(gè)優(yōu)點(diǎn)是T 樣條模型是水密的.多個(gè)NURBS 補(bǔ)丁可以合并成一個(gè)單一的水密T 樣條,任何修剪過(guò)的NURBS 對(duì)象都可以轉(zhuǎn)換為未修剪的T 樣條[21].因此,T 樣條被認(rèn)為是未來(lái)CAD 行業(yè)的新標(biāo)準(zhǔn),將在CAD 與CAE 的集成中發(fā)揮重要作用.Casquero 等[22]利用適合分析的任意度T 樣條曲面進(jìn)行了完全非線性薄殼的結(jié)構(gòu)分析.Casquero 等[23]還利用適合分析的T 樣條來(lái)解決Kirchhoff?Love 殼問(wèn)題.Dimitri 等[24]提出了一種基于T 樣條的等幾何分析方法,應(yīng)用于大變形條件下變形體間的無(wú)摩擦接觸問(wèn)題.

在實(shí)際應(yīng)用中,由于ANCF 單元使用梯度作為節(jié)點(diǎn)向量,其建模過(guò)程較為復(fù)雜.另外,如果將一個(gè)薄板離散為幾個(gè)ANCF 薄板單元,其單元邊界上的梯度不連續(xù),因此格林?拉格朗日應(yīng)變也是不連續(xù)的.而使用IGA 方法可以直接使用CAD 軟件對(duì)分析對(duì)象進(jìn)行幾何建模,而且?guī)缀文Ｐ偷倪B續(xù)性條件可以自然得通過(guò)結(jié)點(diǎn)重復(fù)性保證.為此,本文將在等幾何分析的框架下,開(kāi)展基于T 樣條曲面單元的基爾霍夫薄板動(dòng)力學(xué)分析方法研究.本文的主要貢獻(xiàn)在于:

(1)基于可局部細(xì)化的T 樣條曲面,提出局部細(xì)化策略解決動(dòng)力學(xué)分析中局部應(yīng)變變化較大的問(wèn)題,并根據(jù)T 樣條曲面幾何模型的局部細(xì)化算法創(chuàng)建變自由度系統(tǒng)動(dòng)力學(xué)方程的求解算法;

(2)通過(guò)變網(wǎng)格柔性薄板與剛性球的碰撞問(wèn)題,展示所提出局部細(xì)化T 樣條曲面單元在接觸碰撞問(wèn)題中的應(yīng)用價(jià)值.

1 基于T 樣條的等幾何基爾霍夫薄板

1.1 T 樣條曲面

T 樣條可視為控制網(wǎng)格上帶有T 交叉點(diǎn)的NURBS 曲面[25],可以通過(guò)控制點(diǎn)與T 樣條基函數(shù)Ti(u,v)相乘得到.T 樣條曲面的定義式如下[26]

式中,ST為T(mén) 樣條曲面上物質(zhì)坐標(biāo)為 (u,v)的點(diǎn)的全局坐標(biāo),Pi=(xi,yi,zi) 為第i個(gè)控制點(diǎn)的全局坐標(biāo),n為控制點(diǎn)的數(shù)量,Ti為第i個(gè)控制點(diǎn)所對(duì)應(yīng)的有理T 樣條基函數(shù),其表達(dá)式為

式 中,混合函數(shù)Bi為B 樣條基函數(shù)Ni(u)和Ri(v)的乘積. ωi是第i個(gè)控制點(diǎn)所對(duì)應(yīng)的權(quán)重,W為權(quán)值與基函數(shù)的乘積之和.

與NURBS 不同的是,每一個(gè)T 樣條基函數(shù)都是基于局部結(jié)點(diǎn)矢量而不是全局結(jié)點(diǎn)矢量來(lái)計(jì)算的.局部結(jié)點(diǎn)矢量的建立依賴于T 網(wǎng)格和錨點(diǎn)的建立,根據(jù)這兩者可以得到T 樣條曲面上每一個(gè)控制點(diǎn)所對(duì)應(yīng)的局部結(jié)點(diǎn)矢量和T 網(wǎng)格和錨點(diǎn)的具體定義參見(jiàn)文獻(xiàn)[26-27],并可依據(jù)參考文獻(xiàn)[28]得到三次T 樣條基函數(shù)以及其導(dǎo)數(shù)的表達(dá)式.

1.2 T 樣條薄板的運(yùn)動(dòng)學(xué)描述

由于基爾霍夫薄板忽略橫向剪切變形,法向量始終垂直于中性層,其運(yùn)動(dòng)學(xué)描述可以簡(jiǎn)化為其中性層的運(yùn)動(dòng)學(xué)描述.為此,首先需要采用T 樣條單元對(duì)薄板中性層進(jìn)行運(yùn)動(dòng)學(xué)描述.

T 樣條單元是由T 樣條基函數(shù)的簡(jiǎn)化連續(xù)線構(gòu)成的物理區(qū)域[29].在等幾何分析中,利用IEN 數(shù)組可以確定每個(gè)樣條單元的連接關(guān)系.IEN 中的I 表示數(shù)組是整數(shù)值,EN 是“element nodes”的首字母縮寫(xiě),表示“單元節(jié)點(diǎn)”,該數(shù)組可以將局部基函數(shù)號(hào)和單元號(hào)映射到相應(yīng)的全局控制點(diǎn)號(hào)[29].對(duì)于第i個(gè)T 樣條單元而言,其形函數(shù)以及單元節(jié)點(diǎn)向量ei可以表示為

式中,n是一個(gè)T 樣條單元中包含的控制點(diǎn)的個(gè)數(shù),IEN(j) 為該T 樣條單元的IEN 數(shù)列中的第j項(xiàng),PIEN(j)為第 I EN(j)個(gè)控制點(diǎn)的坐標(biāo),I為 3 ×3的單位矩陣,TIEN(j)(u,v) 為第IEN(j)個(gè)控制點(diǎn)所對(duì)應(yīng)的有理T 樣條基函數(shù).

綜上所述,基爾霍夫薄板的運(yùn)動(dòng)學(xué)描述為

式中, (u,v)為當(dāng)前構(gòu)型下物質(zhì)點(diǎn)r的參數(shù)坐標(biāo),ei為參數(shù) (u,v)所屬的第i個(gè)T 樣條單元的單元節(jié)點(diǎn)向量,Si(u,v)是第i個(gè)T 樣條單元的形函數(shù).

1.3 T 樣條薄板的彈性模型

基于T 樣條曲面的基爾霍夫薄板的總彈性能可表示為[30]

式中,Umid為只與中性層應(yīng)變 εmid有關(guān)的薄板中性層的能量,Uκ為只與彎曲應(yīng)變 κ 有關(guān)的彎曲應(yīng)變能.E為各向同性線彈性材料的廣義胡克定律彈性系數(shù)矩陣.根據(jù)文獻(xiàn)[31], εmid和 κ 可以表示為

式中,r和r0分別表示當(dāng)前構(gòu)型和初始構(gòu)型中,薄板中性層上任意點(diǎn)P(u,v)的坐標(biāo).表示中性層的單位法向量.T0表示中性層的局部笛卡爾坐標(biāo)系(e0)1?(e0)2?(e0)3與曲 面坐標(biāo)系(g0)1?(g0)2?n0的變換矩陣,其表達(dá)式為

彈性力Qs是彈性能U的相對(duì)于單元節(jié)點(diǎn)向量e的導(dǎo)數(shù).對(duì)于一個(gè)基于T 樣條曲面單元的薄板系統(tǒng),其彈性力表達(dá)式為

切線剛度矩陣JQs是彈性力Qs對(duì)e的導(dǎo)數(shù),其表達(dá)式為

圖1 給出了等幾何基爾霍夫薄板的轉(zhuǎn)換示意圖.首先,體積積分簡(jiǎn)化為厚度h與中性層面積分的乘積.然后,當(dāng)前構(gòu)型下物理空間中的面微元 dA被映射到參數(shù)空間中的面微元 dA0,其定義為

圖1 等幾何基爾霍夫薄板積分過(guò)程Fig.1 The process of isogeometric Kirchhoff thin plate integration

最后,對(duì)薄板的積分表達(dá)式可以寫(xiě)為

式中, ωi和 ωj是積分點(diǎn)u和v所對(duì)應(yīng)的積分系數(shù);J1是參數(shù)空間到父單元的轉(zhuǎn)換矩陣的行列式值,J1的表達(dá)式為

2 網(wǎng)格局部細(xì)化

與NURBS 單元網(wǎng)格相比,T 樣條允許在局部細(xì)化單元網(wǎng)格.在處理碰撞接觸等局部發(fā)生應(yīng)變的劇烈變化的問(wèn)題時(shí)具有優(yōu)勢(shì).本節(jié)將討論T 樣條薄板的網(wǎng)格局部細(xì)分及新控制點(diǎn)的計(jì)算方法.

如圖2 所示,一個(gè)邊界為 [ui,ui+1]×[vi,vi+1]的T 樣條單元被均勻分為4 個(gè)單元,將會(huì)產(chǎn)生兩個(gè)新的結(jié)點(diǎn),umid=(ui+ui+1)/2 和vmid=(vi+vi+1)/2.細(xì)化后,初始T 樣條單元的區(qū)域?qū)⒈? 個(gè)新的子單元所替換,其單元邊界分別為以及[umid,ui+1]每一個(gè)T 樣條單元可以被循環(huán)細(xì)化直到達(dá)到停止細(xì)化的指標(biāo).

圖2 局部細(xì)化策略的示意圖Fig.2 Schematic diagram of local refinement strategy

曲面H初始由一個(gè)雙三次T 樣條S1表示,它的控制點(diǎn)數(shù)量為m,笛卡爾空間中的控制點(diǎn)矩陣P的表達(dá)式在式(15)中給出.由于T 網(wǎng)格中一個(gè)交點(diǎn)對(duì)應(yīng)一個(gè)控制點(diǎn),細(xì)化后的幾何模型將生成新的控制點(diǎn),所以細(xì)化后的曲面S2中控制點(diǎn)個(gè)數(shù)增至n,其控制點(diǎn)坐標(biāo)為P?.初始T 樣條曲面S1被稱為是細(xì)化后的T 樣條曲面S2的子空間(表示為S1?S2)

根據(jù)參考文獻(xiàn)[18],可以得到T 樣條混合函數(shù)的局部細(xì)化算法:初始T 樣條曲面上任意控制點(diǎn)Pi的局部結(jié)點(diǎn)矢量為則由ui和vi得到的B 樣條基函數(shù)Ni(u) 和Ri(v)分別為所以其混合函數(shù)Bi可表示為Bi=Ni(u)·Ri(v).如圖2所示,初始T 樣條曲面被細(xì)化后,兩個(gè)新的結(jié)點(diǎn)=umid和=vmid被插入到原始單元的邊界上.控制點(diǎn)Pi所對(duì)應(yīng)的基函數(shù)Ni(u) 和Ri(v)將使用插入u和v的新的局部結(jié)點(diǎn)矢量表示.因此,基函數(shù)Ni(u)可以寫(xiě) 作之和,Ri(v) 可以寫(xiě)作與之和.一個(gè)初始控制點(diǎn)Pi的混合函數(shù)Bi可以由表示

式中,初始T 樣條曲面S1中的第i個(gè)混合函數(shù)Bi被分解為四個(gè)新混合函數(shù)與其對(duì)應(yīng)系數(shù)乘積之和.如果等于細(xì)化后T 樣條曲面第j個(gè)混合函數(shù)那么S1中的混合函數(shù)Bi可以被寫(xiě)作S2中混合函數(shù)的線性組合

因?yàn)镾1?S2,由參數(shù) (u,v)定義的S1上的點(diǎn)r1和S2上的點(diǎn)r2是相等的,也就是r1(u,v)≡r2(u,v).根據(jù)式(1)和式(2),r1和r2可表示為

根據(jù)Bi與B?j之間的轉(zhuǎn)換關(guān)系,可得到之間的轉(zhuǎn)換關(guān)系,轉(zhuǎn)換方程如式(19)和式(20)所示

式中,P4和P?4是分別由控制點(diǎn)組成的矩陣.對(duì)于轉(zhuǎn)換矩陣Mm,其第j行第i列為式(17)中的系數(shù).

根據(jù)式(20),可以得到細(xì)化后T 樣條的齊次坐標(biāo)矩陣然而,本文中T 樣條單元的節(jié)點(diǎn)向量是由控制點(diǎn)笛卡爾坐標(biāo)構(gòu)成的而非齊次坐標(biāo).因此,控制點(diǎn)的笛卡爾坐標(biāo)和權(quán)重的計(jì)算表達(dá)示為

綜上所述,根據(jù)本節(jié)給出的局部細(xì)化策略和算法,可以根據(jù)原T 樣條曲面和插入的結(jié)點(diǎn)得到細(xì)化后的T 樣條曲面,實(shí)現(xiàn)幾何模型拓?fù)浣Y(jié)構(gòu)的變化.

3 變網(wǎng)格T 樣條薄板的求解算法

本文中基于T 樣條的柔性等幾何薄板系統(tǒng)的動(dòng)力學(xué)方程可表示為

圖3 變網(wǎng)格系統(tǒng)動(dòng)力學(xué)分析流程圖Fig.3 Dynamic analysis flow chart of variable mesh system

然后根據(jù)轉(zhuǎn)換矩陣Mm得到細(xì)化后新系統(tǒng)的節(jié)點(diǎn)坐標(biāo),節(jié)點(diǎn)速度和節(jié)點(diǎn)加速度,如下式(24)所示

更新幾何模型后,約束方程也會(huì)發(fā)生變化,所以將會(huì)得到新的布爾矩陣.根據(jù)可以得到新系統(tǒng)的廣義坐標(biāo)廣 義 速 度和廣義加速度其計(jì)算公式如下量將會(huì)得到更新.通過(guò)與的乘積,將會(huì)得到新系統(tǒng)的廣義質(zhì)量陣,廣義外力,廣義彈性力和廣

獲取了上述變量之后,系統(tǒng)動(dòng)力學(xué)方程中的變義接觸力.將更新后的動(dòng)力學(xué)方程輸入到廣義α 算法中進(jìn)行求解,便可得到網(wǎng)格局部細(xì)化、自由度增加后的新系統(tǒng)下一時(shí)間步的廣義變量實(shí)現(xiàn)了對(duì)變自由度系統(tǒng)動(dòng)力學(xué)方程的求解.

4 數(shù)值算例

4.1 受端彎矩作用的懸臂薄板

圖4 為一末端受彎矩M的懸臂薄板.該薄板長(zhǎng)度L=12 m, 寬度b=1 m , 厚度h=0.1 m, 彈性模量E=1.2 MPa,泊松比 ν =0,慣性矩圖5 給出了不同端彎矩下的薄板構(gòu)型圖.根據(jù)文獻(xiàn)[32] 可知,取最大彎矩為當(dāng)末端所受彎矩M=M0,此時(shí)薄板將彎曲成一個(gè)半徑為的圓.當(dāng)末端所受彎矩為M=0.75M0,0.5M0和 0 .25M0時(shí),薄板的端截面轉(zhuǎn)角分別為 1 .5π , π和 0 .5π.

圖4 受端彎矩作用的懸臂薄板Fig.4 Cantilever subjected to end bending moment

圖5 端彎矩作用下,變形后懸臂梁的構(gòu)形圖Fig.5 Configuration of deformed cantilever under external moments

4.2 受內(nèi)壓圓環(huán)

一平面圓環(huán)構(gòu)件,在其中心孔邊緣受沿半徑向外方向的均布載荷P=10 MPa.根據(jù)對(duì)稱性,取圓環(huán)構(gòu)件的四分之一并在兩側(cè)施加滑動(dòng)邊界約束,如圖6(a)所示.圓環(huán)構(gòu)件的內(nèi)與外徑分別為R1=1 m和R2= 2 m,厚度h= 0.1 m.材料密度ρ=7800kg/m3,楊氏模量E=210 GPa,泊松比 ν =0.3.

圖6 四分之一圓環(huán)構(gòu)件及其局部細(xì)化示意圖Fig.6 Diagram of a quarter of circular thin plate and its local refinement diagram

采用T 樣條曲面單元對(duì)構(gòu)件離散,單元網(wǎng)格分布圖如圖6(b) 所示.因?yàn)樽畲竺兹箲?yīng)力出現(xiàn)在AB邊界附近,內(nèi)部區(qū)域被局域細(xì)化.模型共含有28 個(gè)單元和72 個(gè)控制點(diǎn).圖7 給出了局部細(xì)化后構(gòu)件的米塞斯應(yīng)力云圖.

圖7 局部細(xì)化后構(gòu)件米塞斯應(yīng)力云圖Fig.7 von-Mises stress nephogram after local refinement

因?yàn)樵撍憷龑儆谄矫鎲?wèn)題,所以軸向應(yīng)力 σz=0.根據(jù)拉美方程可以得到圓環(huán)在任意半徑r處米塞斯應(yīng)力 σs的理論解

式中, σr為徑向應(yīng)力, σθ為環(huán)向應(yīng)力,K為圓環(huán)外徑與內(nèi)徑之比,K=R2:R1.

在有限元軟件ANSYS 中采用SHELL63 單元運(yùn)行相同算例并對(duì)比構(gòu)件中最大和最小應(yīng)力結(jié)果如表1 所示.可以看出,相較于ANSYS 結(jié)果,IGA 結(jié)果更加接近于理論解而且與理論解的誤差小于1‰.

表1 圓環(huán)構(gòu)件最小和最大米塞斯應(yīng)力Table 1 Minimum and maximum von-Mises stress on circular thin plate

4.3 柔性單擺

該算例測(cè)試了基于T 樣條的等幾何薄板的動(dòng)力學(xué)特性.如圖8 所示,一個(gè)柔性單擺鉸接固定在A 處.柔性薄板單擺在重力加速度g=9.81 m/s2的作用下墜落.薄板的參數(shù)在表2 中給出.

表2 材料參數(shù)Table 2 Material parameters

圖8 一端鉸接的柔性薄板單擺Fig.8 Flexible thin plate pendulum

圖9 給出了含有2 × 2 個(gè)T 樣條單元的薄板在1 s 內(nèi)動(dòng)能Ek、重力勢(shì)能Eg、彈性勢(shì)能Ee及總能量Et的變化情況.可以看出,基于T 樣條的等幾何薄板滿足能量守恒定律.為了檢驗(yàn)本文提出的方法的收斂性,分別使用2 × 2,4 × 4 和8 × 8 的T 樣條曲面單元對(duì)圖8 所示的單擺進(jìn)行離散.圖10 給出了單擺自由端C 的z向位移曲線.可以看出,基于T 樣條的等幾何薄板具有較好的收斂性.圖11 給出了含有2 ×2 個(gè)T 樣條單元的單擺在1 s 內(nèi)的連續(xù)構(gòu)型.

圖9 能量曲線Fig.9 The energy balance curve

圖10 等幾何薄板中T 樣條單元的收斂性Fig.10 Convergence of T-spline element in thin plate

圖11 含有2 × 2 個(gè)T 樣條單元的單擺的連續(xù)構(gòu)型Fig.11 Configuration of the pendulum with 2 × 2 T-spline elements

4.4 剛性球與柔性板的接觸

本算例將T 樣條局部細(xì)化算法應(yīng)用于柔體動(dòng)力學(xué)分析中.如圖12 所示,一個(gè)剛性球從四邊固定的薄板中心正上方自由下落.接觸發(fā)生后,對(duì)薄板中心受沖擊區(qū)域進(jìn)行局部細(xì)化.剛性球與板的參數(shù)表3中給出.

表3 接觸算例的參數(shù)Table 3 Parameters of the contact example

圖12 自由墜落剛性球與柔性薄板的碰撞Fig.12 The collision between a rigid ball and a flexible thin plate

在本算例中,采用罰值法來(lái)完成剛體球與柔性板的接觸實(shí)現(xiàn).在每個(gè)T 樣條單元中,剛體球與均勻檢測(cè)點(diǎn)的距離用式(27)表示

式中,p為嵌入深度,r(u,v)為薄板中性層上的接觸檢測(cè)點(diǎn),rs為任意點(diǎn)檢測(cè)點(diǎn)在薄板上表面所對(duì)應(yīng)的物質(zhì)點(diǎn),rc為球心位置,h為薄板厚度,R為球心半徑,n表示接觸檢測(cè)點(diǎn)的單位法向量.

一旦小球與薄板發(fā)生接觸,開(kāi)始計(jì)算薄板所受的接觸力.計(jì)算切向接觸力時(shí),采用文獻(xiàn)[33]中計(jì)及臨界滑動(dòng)速度v0的平滑化庫(kù)倫摩擦模型.接觸探測(cè)點(diǎn)i的法向接觸力Fin和切向接觸力Fit表達(dá)式分別為

式中,vp為嵌入速度,k和c表示剛度系數(shù)與阻尼系數(shù).在式中, μ 為薄板與球之間的摩擦系數(shù),vt為切向接觸速度,v0為假定的臨界滑動(dòng)速度,vet為單元切向接觸速度.

圖13 介紹了三種網(wǎng)格構(gòu)型圖.網(wǎng)格1 和網(wǎng)格2 都對(duì)表面進(jìn)行了均勻的網(wǎng)格細(xì)化,其單元數(shù)量分別為14 × 14 和10 × 10.網(wǎng)格3 的初始網(wǎng)格較為粗糙,含有10 × 10 個(gè)單元.當(dāng)接觸發(fā)生后,薄板的接觸區(qū)域被局部細(xì)化,單元數(shù)量由100 變?yōu)?12 個(gè).三組薄板的詳細(xì)信息見(jiàn)表4.

圖13 三組不同的網(wǎng)格構(gòu)型Fig.13 The mesh refinement for three groups

表4 三組薄板的詳細(xì)信息Table 4 The detailed information of three groups

圖14 給出了三組小球球心在1 s 內(nèi)的z向位移曲線.圖15 和圖16 則為三組薄板的質(zhì)心z向位移曲線與彈性能曲線及其局部放大圖.在t=0.216 s時(shí),剛性球與柔性板首次接觸,對(duì)薄板進(jìn)行如圖13所示的3 種不同網(wǎng)格細(xì)化策略.三組球的質(zhì)心垂直位移具有較好的一致性,薄板質(zhì)心的垂直位移曲線在1 s 內(nèi)也呈現(xiàn)相似的趨勢(shì).從薄板的質(zhì)心位移與彈性能計(jì)算結(jié)果來(lái)看,允許局部網(wǎng)格細(xì)化的T 樣條等幾何薄板可以達(dá)到與退化為NURBS 曲面的均勻網(wǎng)格板相同的計(jì)算精度,且相較于網(wǎng)格2,局部細(xì)化的網(wǎng)格3 所對(duì)應(yīng)曲線的變化趨勢(shì)更接近網(wǎng)格1.

圖14 球心z 向位移Fig.14 Vertical displacement of the center of the ball

圖15 薄板質(zhì)心的z 向位移Fig.15 Vertical displacement of the plate’s centroid

圖16 薄板彈性能曲線Fig.16 The elastic energy curve of thin plate

為了比較不同網(wǎng)格構(gòu)型的薄板對(duì)計(jì)算效率的影響,圖17 給出了三組薄板的仿真時(shí)間柱狀圖.

圖17 計(jì)算消耗時(shí)間Fig.17 Time consumption for three groups

三組算例均是在配備Intel Core i5 CPU 和8GB RAM 的筆記本電腦上執(zhí)行的.由圖17 可知,網(wǎng)格2 與網(wǎng)格3 用時(shí)接近,約為網(wǎng)格1 計(jì)算用時(shí)的一半.值得注意的是,由于在接觸區(qū)域附近執(zhí)行了網(wǎng)格細(xì)化,使得求解收斂更快,網(wǎng)格3 雖然自由度數(shù)略多于網(wǎng)格2,用時(shí)反而更少.這也進(jìn)一步體現(xiàn)了T 樣條單元網(wǎng)格局部細(xì)化的優(yōu)勢(shì).

圖18 展示了帶有112 個(gè)T 樣條單元局部細(xì)化薄板在不同時(shí)刻的米塞斯應(yīng)力云圖,分別為t=0.18 s接觸前的薄板發(fā)生最大變形,t=0.218 s發(fā)生接觸后首次使用細(xì)化后模型的時(shí)刻,t=0.548 s接觸后的發(fā)生最大變形的時(shí)刻以及仿真結(jié)束時(shí)刻t=1 s.

圖18 局部細(xì)化薄板的米塞斯應(yīng)力云圖Fig.18 von-Mises stress distribution of the locally refined thin plate

5 結(jié)論

本文提出了一種基于T 樣條曲面的等幾何分析方法,建立了使用T 樣條曲面單元離散的基爾霍夫薄板模型.單元基函數(shù)和節(jié)點(diǎn)坐標(biāo)分別為T(mén) 樣條基函數(shù)和控制點(diǎn)坐標(biāo),無(wú)需網(wǎng)格劃分,既保證模型的幾何精確,又能在沒(méi)有約束方程的情況下保證期望的連續(xù)性條件.給出了基于T 樣條的薄板運(yùn)動(dòng)學(xué)模型、彈性力模型及其雅克比矩陣的計(jì)算方法.為了體現(xiàn)T 樣條局部細(xì)化特性在減少單元數(shù)目、提高計(jì)算效率方面的優(yōu)勢(shì),提出了一種基于T 樣條單元的局部網(wǎng)格細(xì)化算法,包括單元網(wǎng)格拓?fù)涞母淖兒拖鄳?yīng)的新控制點(diǎn)坐標(biāo)計(jì)算方法.將該細(xì)化算法與廣義α 法相結(jié)合,建立了帶有網(wǎng)格局部細(xì)化的變自由度系統(tǒng)動(dòng)力學(xué)方程的求解算法.通過(guò)受端彎矩的懸臂薄板以及環(huán)形構(gòu)件受內(nèi)壓的靜力學(xué)算例證明了T 樣條單元彈性力模型的正確性.柔性擺算例驗(yàn)證了T 樣條單元在動(dòng)力學(xué)問(wèn)題中的收斂性和機(jī)械能守恒特性.最后,為驗(yàn)證T 樣條局部細(xì)化特性在模擬接觸碰撞等柔性體局部發(fā)生應(yīng)變劇烈變化等問(wèn)題中的優(yōu)勢(shì),建立了剛性球落在柔性板上的動(dòng)力學(xué)實(shí)例并進(jìn)行了仿真.首先,采用10 × 10 單元網(wǎng)格對(duì)柔性板進(jìn)行離散.接觸發(fā)生后,對(duì)沖擊區(qū)域進(jìn)行局部細(xì)化,得到112 個(gè)單元的網(wǎng)格.對(duì)10 × 10 和14 × 14 兩種網(wǎng)格進(jìn)行了仿真,并將仿真結(jié)果作為基準(zhǔn).可以看到,112 單元網(wǎng)格的計(jì)算結(jié)果與14 × 14 網(wǎng)格具有良好的一致性,但所消耗的時(shí)間只有其1/2.以上算例證明了基于T 樣條的局部網(wǎng)格細(xì)化等幾何薄板在柔性多體系統(tǒng)動(dòng)力學(xué)分析中的應(yīng)用價(jià)值.