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

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

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

引言

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

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

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

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

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

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

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

1.1 T 樣條曲面

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

式中,ST為T 樣條曲面上物質(zhì)坐標為 (u,v)的點的全局坐標,Pi=(xi,yi,zi) 為第i個控制點的全局坐標,n為控制點的數(shù)量,Ti為第i個控制點所對應的有理T 樣條基函數(shù),其表達式為

式 中,混合函數(shù)Bi為B 樣條基函數(shù)Ni(u)和Ri(v)的乘積. ωi是第i個控制點所對應的權重,W為權值與基函數(shù)的乘積之和.

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

1.2 T 樣條薄板的運動學描述

由于基爾霍夫薄板忽略橫向剪切變形,法向量始終垂直于中性層,其運動學描述可以簡化為其中性層的運動學描述.為此,首先需要采用T 樣條單元對薄板中性層進行運動學描述.

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

式中,n是一個T 樣條單元中包含的控制點的個數(shù),IEN(j) 為該T 樣條單元的IEN 數(shù)列中的第j項,PIEN(j)為第 I EN(j)個控制點的坐標,I為 3 ×3的單位矩陣,TIEN(j)(u,v) 為第IEN(j)個控制點所對應的有理T 樣條基函數(shù).

綜上所述,基爾霍夫薄板的運動學描述為

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

1.3 T 樣條薄板的彈性模型

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

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

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

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

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

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

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

最后,對薄板的積分表達式可以寫為

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

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

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

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

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

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

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

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

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

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

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

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

綜上所述,根據(jù)本節(jié)給出的局部細化策略和算法,可以根據(jù)原T 樣條曲面和插入的結點得到細化后的T 樣條曲面,實現(xiàn)幾何模型拓撲結構的變化.

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

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

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

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

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

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

4 數(shù)值算例

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

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

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

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

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

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

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

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

圖7 局部細化后構件米塞斯應力云圖Fig.7 von-Mises stress nephogram after local refinement

因為該算例屬于平面問題,所以軸向應力 σz=0.根據(jù)拉美方程可以得到圓環(huán)在任意半徑r處米塞斯應力 σs的理論解

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

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

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

4.3 柔性單擺

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

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

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

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

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

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

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

4.4 剛性球與柔性板的接觸

本算例將T 樣條局部細化算法應用于柔體動力學分析中.如圖12 所示,一個剛性球從四邊固定的薄板中心正上方自由下落.接觸發(fā)生后,對薄板中心受沖擊區(qū)域進行局部細化.剛性球與板的參數(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

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

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

一旦小球與薄板發(fā)生接觸,開始計算薄板所受的接觸力.計算切向接觸力時,采用文獻[33]中計及臨界滑動速度v0的平滑化庫倫摩擦模型.接觸探測點i的法向接觸力Fin和切向接觸力Fit表達式分別為

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

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

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

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

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

圖17 計算消耗時間Fig.17 Time consumption for three groups

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

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

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

5 結論

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