劉登學(xué)張友良劉高敏,?

?(中國科學(xué)院武漢巖土力學(xué)研究所巖土力學(xué)與工程國家重點(diǎn)實(shí)驗(yàn)室，武漢430071)

?(中國科學(xué)院大學(xué)，北京100049)

基于適合分析T樣條的高階數(shù)值流形方法1)

劉登學(xué)?,?,2)張友良?劉高敏?,?

?(中國科學(xué)院武漢巖土力學(xué)研究所巖土力學(xué)與工程國家重點(diǎn)實(shí)驗(yàn)室，武漢430071)

?(中國科學(xué)院大學(xué)，北京100049)

數(shù)值流形方法是一種非常靈活的數(shù)值計(jì)算方法，連續(xù)體的有限單元方法和塊體系統(tǒng)的非連續(xù)變形分析方法只是這一數(shù)值方法的特例.數(shù)值流形方法中高階位移函數(shù)的構(gòu)造可通過提高權(quán)函數(shù)的階次來實(shí)現(xiàn)，這種方法往往需要沿單元邊界配置適當(dāng)?shù)倪厓?nèi)節(jié)點(diǎn)，這些結(jié)點(diǎn)的出現(xiàn)增加了前處理的復(fù)雜性，特別是對于大型復(fù)雜的空間問題.另一方面，在數(shù)值流形方法中可通過縮小單元尺寸(h加密)來提高求解精度.當(dāng)模擬裂紋擴(kuò)展時(shí)，這種細(xì)化策略可用來克服裂紋尖端的奇異性.一個(gè)傳統(tǒng)的解決方案是細(xì)化整個(gè)網(wǎng)格，但這會導(dǎo)致計(jì)算效率的顯著降低.將適合分析的T樣條(analysis-suitable T-spline，AST)引入數(shù)值流形方法中來建立高階數(shù)值流形方法的分析格式，有效的避免了該問題的出現(xiàn).AST樣條基函數(shù)具有線性無關(guān)，單位分解，局部加密等許多重要性質(zhì)，使得其非常適合用于工程設(shè)計(jì)及分析.在引入AST樣條后，可通過改變數(shù)學(xué)覆蓋的構(gòu)造形式建立不同階次的數(shù)值流形方法分析格式；AST樣條自身的局部加密性質(zhì)也使得數(shù)值流形方法中的數(shù)學(xué)網(wǎng)格局部加密更容易實(shí)現(xiàn).算例結(jié)果表明：隨著AST樣條基函數(shù)階次的提高，數(shù)值流形方法的計(jì)算結(jié)果有了明顯的改善；基于AST樣條基函數(shù)的數(shù)值流形方法在保持計(jì)算精度的前提下降低了自由度的數(shù)量.

高階數(shù)值流形方法，線性相關(guān)，適合分析的T樣條，局部加密

引言

1992年，石根華提出了在巖石工程中具有重要影響的數(shù)值流形法[1-2](numerical manifold method, NMM).該方法非常適合斷續(xù)節(jié)理巖體的模擬，是一個(gè)在統(tǒng)一數(shù)學(xué)框架下的巖石小變形開裂與非連續(xù)大位移分析的一體化計(jì)算方法.數(shù)值流形方法自提出以來，已在巖土工程諸多領(lǐng)域得到成功應(yīng)用，如隧道的開挖模擬[3-5]、巖體的錨固[6-7]、巖體裂紋擴(kuò)展[8-12]以及巖體的飽和滲流分析[13].通?？梢酝ㄟ^以下幾種方式來提高數(shù)值流形方法的求解精度：覆蓋函數(shù)采用高階次的多項(xiàng)式；提高權(quán)函數(shù)的階次；減小流形單元的尺寸，即單元加密.姜清輝等[14]采用一階覆蓋函數(shù)代替常覆蓋函數(shù)，建立了三維高階數(shù)值流形方法的分析格式.蘇海東等[15-16]基于平面三角形數(shù)學(xué)網(wǎng)格和多項(xiàng)式覆蓋函數(shù)，提出高階流形法的兩種初應(yīng)力的處理方法.Wang等[17]提出了一種二階數(shù)值流形方法模型用以解決具有自由表面水流的非線性問題.然而，高階數(shù)值流形方法在提高計(jì)算精度的同時(shí)，也會引起總體剛度矩陣奇異，即線性相關(guān).針對該問題Zheng等[18]提出在數(shù)值流形方法中物理覆蓋上的局部位移函數(shù)采用一階泰勒展開式，建立了線性無關(guān)的高階數(shù)值流形方法，郭朝旭等[19]提出了改進(jìn)的LDLT算法，可快速穩(wěn)定的求得一個(gè)特解.蔡永昌等[20]提出了一種基于獨(dú)立覆蓋的高階流形方法，消除了高階流形方法特有的線性相關(guān)帶來的總體剛度矩陣奇異性的問題.Fan等[21]將位移的導(dǎo)數(shù)作為自由度引入到數(shù)值流形方法中，建立了一種新型的高階數(shù)值流形方法，避免了線性相關(guān)的問題.在提高權(quán)函數(shù)階次方面，通常是通過沿單元邊界配置適當(dāng)?shù)膬?nèi)結(jié)點(diǎn)來實(shí)現(xiàn)，這些結(jié)點(diǎn)的出現(xiàn)增加了前處理的復(fù)雜性，特別是對于大型復(fù)雜的空間問題.

等幾何分析 (isogeometric analysis)是 Hughes等[22]于 2005年提出的一個(gè)將計(jì)算機(jī)鋪助設(shè)計(jì)和計(jì)算機(jī)輔助工程相接合的一種數(shù)值分析方法.它的基本思想是將用來表示幾何造型的基函數(shù)同時(shí)用來定義未知場變量，如位移、溫度等.最常用的基函數(shù)是非均勻有理B樣條(non-uniform rational B spline, NURBS).NURBS有很多優(yōu)點(diǎn)：可以進(jìn)行任意表面的建模；擁有很多優(yōu)良的數(shù)學(xué)特性，如非負(fù)性、單位分解特性、基函數(shù)的線性無關(guān)性、局部支撐性等.隨著研究的深入，NURBS本身的不足也不斷暴露出來[23]，如 NURBS采用張量積形式，使得局部網(wǎng)格細(xì)化效率不高，很多形狀無法用一個(gè)簡單封閉的NURBS面表示；面與面交界處有間隙或重疊等.為了克服這些不足，許多學(xué)者開始將T樣條用于等幾何分析中.Sederberg等[24]最早將T樣條引入計(jì)算機(jī)制圖，并且提出了T樣條的局部加密算法；Bazilevs等[25]探索了基于T樣條基函數(shù)的等幾何分析，用于簡單的流體力學(xué)和結(jié)構(gòu)力學(xué)的計(jì)算，獲得了很好的效果.Zhang等[26]基于等幾何分析的思想，推導(dǎo)了基于二次B樣條的9節(jié)點(diǎn)數(shù)值流形方法分析格式，提出了一種數(shù)值流形T樣條局部加密方法.有些情況下，T樣條的局部加密會出現(xiàn)不適合分析的單元；另外T樣條的混合函數(shù)的線性無關(guān)性在一般T網(wǎng)格上不能保證，可能導(dǎo)致剛度矩陣奇異.為此Li等[27]提出了適合分析的T樣條(analysis-suitable T-spline，AST)來克服上述問題.本文將AST樣條引入到數(shù)值流形方法中，建立了基于AST樣條函數(shù)的高階數(shù)值流形方法分析格式.

1 數(shù)值流形方法簡介

1.1 數(shù)值流形方法基本概念

數(shù)值流形方法是基于三個(gè)基本概念建立起來的：數(shù)學(xué)覆蓋、物理覆蓋、流形單元.數(shù)學(xué)覆蓋系統(tǒng)是許多可重疊小片的并集，覆蓋整個(gè)問題區(qū)域.每一個(gè)小片稱為一個(gè)數(shù)學(xué)覆蓋，記為Mi(i=1,2,··,nM).物理覆蓋系統(tǒng)由數(shù)學(xué)覆蓋和物理網(wǎng)格兩者組成，物理網(wǎng)格包括材料體的邊界、裂縫、塊體邊界、不同材料區(qū)域的交接面.如果裂縫或者塊體邊界把一個(gè)數(shù)學(xué)覆蓋Mi分成兩個(gè)或更多的完全不連續(xù)區(qū)域，這些區(qū)域稱為物理覆蓋，記為因此物理覆蓋是不連續(xù)裂縫對數(shù)學(xué)覆蓋的再剖分.正是由于引入了數(shù)學(xué)覆蓋和物理覆蓋兩套覆蓋系統(tǒng)，數(shù)值流形方法可以統(tǒng)一地解決連續(xù)和非連續(xù)問題.流形單元定義為幾個(gè)物理覆蓋的公共部分.

為了更清楚的說明這3個(gè)概念，給出如圖1所示的一個(gè)例子.對于圖中含一條裂紋的不規(guī)則多邊形板，數(shù)學(xué)覆蓋系統(tǒng)采用規(guī)則矩形網(wǎng)格，共享結(jié)點(diǎn)i的4個(gè)矩形組成一個(gè)數(shù)學(xué)覆蓋，記為Mi.數(shù)學(xué)覆蓋M1與裂紋相交形成兩個(gè)物理覆蓋數(shù)學(xué)覆蓋M2，M3，M4，M5沒有與物理網(wǎng)格相交，各自形成一個(gè)物理覆蓋p2，p3，p4，p5，物理覆蓋與數(shù)學(xué)覆蓋相同.如圖中陰影所示，流形單元E(p2,p3,p4,p5)為物理覆蓋p2，p3，p4，p5的公共區(qū)域，而流形單元E(p6,p7,p8,p9)為物理覆蓋p6，p7，p8，p9的公共區(qū)域.流形單元的形狀可以是任意的,但是其對應(yīng)的物理覆蓋是唯一的.

圖1 數(shù)值流形方法中的覆蓋和單元Fig.1 A simple example to illustrate the basic concepts of NMM

1.2 覆蓋位移函數(shù)

數(shù)值流形方法中數(shù)學(xué)覆蓋獨(dú)立于物理網(wǎng)格，因此一般常用的有限元網(wǎng)格總可以被引入到數(shù)值流形方法中.這些有限元網(wǎng)格必須足夠大來覆蓋整個(gè)問題區(qū)域.為了使權(quán)函數(shù)的建立更加簡單，程序?qū)崿F(xiàn)更加容易，在二維數(shù)值流形方法中通常采用三角形或者是四邊形網(wǎng)格.如圖1所示，數(shù)值流形方法中采用規(guī)則四邊形網(wǎng)格的數(shù)學(xué)覆蓋，雙線性的四邊形有限元形函數(shù)作為數(shù)學(xué)覆蓋的權(quán)函數(shù).

這些覆蓋函數(shù)通過定義在數(shù)學(xué)覆蓋上的權(quán)函數(shù)wi(x,y)連接在一起

在整個(gè)材料體上總體位移函數(shù)可以表示為

1.3 數(shù)值流形方法模擬不連續(xù)問題

數(shù)值流形方法模擬不連續(xù)問題時(shí)裂紋尖端存在奇異性，由于多項(xiàng)式不能準(zhǔn)確描述裂紋尖端的應(yīng)力場，往往會造成求解精度差的問題.這一問題可以通過將附加函數(shù)加入到覆蓋函數(shù)中來解決.如圖2所示，含有裂紋尖端的數(shù)學(xué)覆蓋M1,M2,M3和M4各自形成一個(gè)物理覆蓋，這些物理覆蓋稱之為奇異的物理覆蓋(singular physical cover)[28].這些物理覆蓋的覆蓋函數(shù)通過加入附加函數(shù)來考慮不連續(xù)尖端場的影響

式中Uj為奇異物理覆蓋的附加函數(shù)

式中，cj為附加函數(shù)的未知量矩陣，ns為裂紋尖端所在單元的奇異物理覆蓋的數(shù)量.Φ為奇異矩陣的基函數(shù)矩陣，可以表示為

式中，(r,θ)為裂紋尖端局部坐標(biāo)系中的坐標(biāo)，r是任意一點(diǎn)到裂紋尖端的距離，θ是裂紋擴(kuò)展角.基于這種強(qiáng)化方法，裂紋尖端可落在流形單元的任意位置.即使在相對粗糙的數(shù)學(xué)網(wǎng)格中，我們也可以精確地得到應(yīng)力強(qiáng)度因子.

圖2 裂紋尖端處的奇異物理覆蓋Fig.2 Singular physical covers at crack tip in NMM

2 T樣條

2.1T樣條基函數(shù)

圖3 一個(gè)T網(wǎng)格Fig.3 A T-mesh

T樣條基函數(shù)是由其自身的局部節(jié)點(diǎn)向量定義的，這些節(jié)點(diǎn)向量是從T網(wǎng)格中推導(dǎo)出來的.T樣條基函數(shù)的局部節(jié)點(diǎn)向量按照以下規(guī)則確定：當(dāng)階次p為奇數(shù)時(shí)，先將錨點(diǎn)位置所代表的節(jié)點(diǎn)值存入局部向量中，然后再分別沿錨點(diǎn)的上下左右4個(gè)方向，由錨點(diǎn)位置出發(fā)，記錄遇到的(p+1)/2個(gè)正交邊并將每一條正交邊對應(yīng)的節(jié)點(diǎn)值存入相應(yīng)的局部向量中；當(dāng)p為偶數(shù)時(shí)，分別沿錨點(diǎn)的上下左右4個(gè)方向，由錨點(diǎn)位置出發(fā)，直到遇到(p/2+1)個(gè)正交邊，將每一條正交邊對應(yīng)的節(jié)點(diǎn)值存入相應(yīng)的局部向量中.如圖3中的錨點(diǎn)A，當(dāng)p=1時(shí)，其對應(yīng)的局部向量為ui=[u4,u5,u8]和vi=[v3,v4,v5]，該基函數(shù)的支撐域(support)見圖4.當(dāng)p=2時(shí)，矩形中心點(diǎn)D為一錨點(diǎn)，相應(yīng)基函數(shù)的局部向量為ui=[u3,u4,u5,u8]和vi=[v3,v4,v5,v6]，其支撐域見圖5.

對于一個(gè)給定的T網(wǎng)格和階次p，利用上述過程確定T樣條局部向量ui和vi之后，T樣條基函數(shù)Ti(u,v)可定義為

其中N[ui](u)和N[vi](v)可由Cox-de Boor[29-30]遞推公式導(dǎo)出

圖4 錨點(diǎn)A對應(yīng)的局部節(jié)點(diǎn)向量Fig.4 Inferring local vectors for anchorA

圖5 錨點(diǎn)D對應(yīng)的局部節(jié)點(diǎn)向量Fig.5 Inferring local vectors for anchorD

2.2T結(jié)點(diǎn)擴(kuò)展

一個(gè)T結(jié)點(diǎn)擴(kuò)展是一條線段，通常包括一個(gè)面擴(kuò)展和一個(gè)邊擴(kuò)展.從T結(jié)點(diǎn)自身出發(fā)沿著它缺失的那一個(gè)方向進(jìn)行延長直到于第[(p+1)/2]個(gè)正交邊相交所形成的延長線段即為T結(jié)點(diǎn)的面擴(kuò)展，而T結(jié)點(diǎn)的線擴(kuò)展是從T結(jié)點(diǎn)沿著相反的方向出發(fā)直到于第[p/2]個(gè)正交邊相交所經(jīng)歷的線段.p為T樣條的階次，[(p+1)/2]和[p/2]為分別對(p+1)/2和p/2取整.稱兩個(gè)T節(jié)點(diǎn)擴(kuò)展是相交的，如果它們相交于一個(gè)都在它們內(nèi)部的點(diǎn)；稱兩個(gè)T結(jié)點(diǎn)擴(kuò)展是接觸的，如果它們相交于一點(diǎn)，并且該點(diǎn)為它們至少一個(gè)的端點(diǎn).如圖6所示，點(diǎn)虛線為T網(wǎng)格中T結(jié)點(diǎn)的面擴(kuò)展，線段虛線為T結(jié)點(diǎn)的線擴(kuò)展.可以看出T結(jié)點(diǎn)A和T結(jié)點(diǎn)C的擴(kuò)展是相交的，T結(jié)點(diǎn)A和T結(jié)點(diǎn)B的擴(kuò)展是接觸的.當(dāng)把所有的T結(jié)點(diǎn)擴(kuò)展加入到T網(wǎng)格中后，T網(wǎng)格變?yōu)閿U(kuò)展的T網(wǎng)格，記作Text.

2.3 適合分析的T樣條

適合分析的 T樣條 (analysis-suitable T-spline，AST)定義在一個(gè)限制的T網(wǎng)格上，即T網(wǎng)格中不存在豎直的T結(jié)點(diǎn)擴(kuò)展與水平的T結(jié)點(diǎn)擴(kuò)展相交或接觸的情況.AST樣條繼承了T樣條所有重要的數(shù)學(xué)特性而且具有高效的局部細(xì)分性.AST樣條基函數(shù)具有以下重要性質(zhì)：

圖6 擴(kuò)展T網(wǎng)格Fig.6 The extended T-mesh.Face extensions are represented by dotted arrows and edge extensions are represented by dashed arrows.The T-junctions are denoted by black dots

(1)單位分解性,即

(2)非負(fù)性

(3)線性無關(guān)性

(4)緊支性

圖7(a)中所示的T網(wǎng)格是適合分析的，因?yàn)樵谙鄳?yīng)的Text不存在豎直的T結(jié)點(diǎn)擴(kuò)展與水平的T結(jié)點(diǎn)擴(kuò)展相交或接觸的情況.

3 基于AST樣條基函數(shù)的高階數(shù)值流形方法

高階位移函數(shù)的構(gòu)造可通過提高權(quán)函數(shù)的階次來實(shí)現(xiàn)，這種方法往往需要沿單元邊界配置適當(dāng)?shù)倪厓?nèi)節(jié)點(diǎn)，這些結(jié)點(diǎn)的出現(xiàn)增加了前處理的復(fù)雜性，特別是對于大型復(fù)雜的空間問題.本節(jié)提出采用AST樣條基函數(shù)來構(gòu)造數(shù)值流形方法中的權(quán)函數(shù)，有效地避免了該問題的出現(xiàn).

圖7 (a)適合分析的T網(wǎng)格(b)擴(kuò)展T網(wǎng)格Fig.7(a)An analysis-suitable T-mesh(b)the extended T-mesh

在數(shù)值流形方法中，采用如圖8所示的T網(wǎng)格作為數(shù)學(xué)覆蓋系統(tǒng)，圖中黑色圓點(diǎn)為T樣條基函數(shù)的錨點(diǎn)，在此表示一個(gè)數(shù)學(xué)覆蓋的位置.當(dāng)p=1，2和3時(shí)，T網(wǎng)格中不存在豎直的T結(jié)點(diǎn)擴(kuò)展與水平的T結(jié)點(diǎn)擴(kuò)展相交或接觸的情況，所以此T網(wǎng)格產(chǎn)生的T樣條為AST樣條.將AST樣條的基函數(shù)作為數(shù)值流形方法中的權(quán)函數(shù)，數(shù)學(xué)覆蓋的構(gòu)造形式及范圍可按3.2節(jié)中T樣條基函數(shù)的相應(yīng)規(guī)則進(jìn)行確定.以右側(cè)裂紋尖端所在單元為例，當(dāng)T樣條基函數(shù)的階次為1時(shí)，如圖9(a)所示緊鄰的4個(gè)矩形網(wǎng)格形成一個(gè)數(shù)學(xué)覆蓋，流形單元為相應(yīng)的4個(gè)數(shù)學(xué)覆蓋形成的物理覆蓋的重疊區(qū)域.當(dāng)T樣條基函數(shù)的階次為2時(shí)，如圖9(b)所示緊鄰的9個(gè)矩形網(wǎng)格形成一個(gè)數(shù)學(xué)覆蓋，流形單元為相應(yīng)的9個(gè)數(shù)學(xué)覆蓋形成的物理覆蓋的重疊區(qū)域.類似的，當(dāng)T樣條基函數(shù)的階次為3時(shí)，如圖9(c)所示緊鄰的16個(gè)矩形網(wǎng)格形成一個(gè)數(shù)學(xué)覆蓋，流形單元為相應(yīng)的16個(gè)數(shù)學(xué)覆蓋形成的物理覆蓋的重疊區(qū)域.

p=1時(shí)，{x1,x2,x3}和{y1,y2,y3}分別為x,y方向的單調(diào)不減的實(shí)數(shù)序列，對應(yīng)的權(quán)函數(shù)可通過式(1),式(2)和式(9)推導(dǎo)，權(quán)函數(shù)可以表示為

圖8 采用AST網(wǎng)格的數(shù)值流形方法Fig.8 Numerical manifold method with AST mesh

圖9 不同階次的數(shù)學(xué)覆蓋Fig.9 Mathematical cover with dif f erent order

基于AST樣條建立的數(shù)值流形方法，權(quán)函數(shù)具有AST樣條基函數(shù)的重要性質(zhì)，滿足數(shù)值流形方法中對權(quán)函數(shù)的要求.此外，權(quán)函數(shù)具有Cp-1連續(xù)性，權(quán)值在3種不同階次數(shù)學(xué)覆蓋區(qū)域中的分布見圖10～圖12.需要指出的是，當(dāng)T網(wǎng)格中不存在T結(jié)點(diǎn)且階次為1時(shí)，AST基函數(shù)等價(jià)于一階Lagrange插值函數(shù).

圖10 p=1時(shí)權(quán)值在一個(gè)數(shù)學(xué)覆蓋中的分布Fig.10 Distribution of weights in a mathematical cover whenp=1

圖11 p=2時(shí)權(quán)值在一個(gè)數(shù)學(xué)覆蓋中的分布Fig.11 Distribution of weights in a mathematical cover whenp=2

圖12 p=3時(shí)權(quán)值在一個(gè)數(shù)學(xué)覆蓋中的分布Fig.12 Distribution of weights in a mathematical cover whenp=3

4 數(shù)值算例

4.1 Timoshenko懸臂梁

懸臂梁模型如圖13所示，左側(cè)為固定端，右側(cè)施加剪切分布力P=1.模型參數(shù)：D=1，L=5；平面應(yīng)變問題，楊氏模量E=1000，泊松比v=0.該問題的解析解為

圖13 Timoshenko懸臂梁Fig.13 Timoshenko cantilever beam

采用如圖14所示的T網(wǎng)格作為數(shù)值流形方法中的數(shù)學(xué)覆蓋系統(tǒng)，由于網(wǎng)格中不存在T結(jié)點(diǎn)，故T網(wǎng)格是適合分析的.分別采用不同階次(p=1,2,3)的AST樣條基函數(shù)構(gòu)造數(shù)值流形方法中數(shù)學(xué)覆蓋的權(quán)函數(shù).

圖14 數(shù)學(xué)覆蓋網(wǎng)格圖Fig.14 Mathematical covers

圖15和圖16分別給出了懸臂梁底邊(y=0)位移及應(yīng)力隨位置的變化圖.從圖中可看出：當(dāng)采用一次T樣條基函數(shù)時(shí)，數(shù)值流形方法(NMM based on linear AST-spline，NMM_T1)的計(jì)算結(jié)果與解析解存在一定的計(jì)算誤差；采用高階次的 T樣條基函數(shù)后，相應(yīng)的高階數(shù)值流形方法(NMM based on quadratic AST-spline，NMM T2；NMM based on cubic AST-spline，NMM_T3)與解析解吻合得很好.

圖17為不同密度T網(wǎng)格下，懸臂梁應(yīng)變能的計(jì)算誤差.應(yīng)變能的計(jì)算誤差定義為

其中U為應(yīng)變能的理論解，Uh為數(shù)值流形方法的求解結(jié)果.從圖中可以隨著T網(wǎng)格的加密，基于3種不同階次AST基函數(shù)的數(shù)值流形方法均是收斂的，但收斂速度隨著AST基函數(shù)的階次提高不斷增大.

圖15 y=0時(shí)y方向位移變化圖Fig.15ydisplacement versus position wheny=0

圖16 y=0時(shí)σx變化圖Fig.16 σxversus position wheny=0

圖17 懸臂梁應(yīng)變能誤差Fig.17 Error in strain energy of the cantilever beam

4.2 含中心邊裂紋的有限板受單向拉伸

如圖18所示有限矩形板，中部有一條長度為a的邊裂紋，矩形板單向受力σ=4，長2h=3，寬b=1.此問題的理論解可以表示為

式中

圖18 含單邊裂紋的有限板Fig.18 Finite plate with crack

圖19 數(shù)學(xué)覆蓋網(wǎng)格加密Fig.19 Mathematical mesh refinemen

取不同裂紋長度a，利用相互作用積分法[31]計(jì)算不同裂紋長度下裂紋強(qiáng)度因子 (stress intensity factors，SIFs).采用圖19所示的兩種數(shù)學(xué)加密網(wǎng)格,(a)為針對裂紋區(qū)域的數(shù)學(xué)網(wǎng)格全局加密，(b)為針對裂紋區(qū)域的局部加密.裂紋強(qiáng)度因子的計(jì)算結(jié)果見表1～表3.圖20和圖21分別為兩種不同加密策略下裂紋強(qiáng)度因子的計(jì)算誤差.對比可以得出：局部加密策略要比全局加密策略更有優(yōu)勢，即在計(jì)算精度基本一致的情況下采用局部加密數(shù)學(xué)網(wǎng)格的自由度要小于全局加密數(shù)學(xué)網(wǎng)格的自由度；在相同數(shù)學(xué)覆蓋網(wǎng)格下，隨著AST樣條基函數(shù)階次的提高，裂紋強(qiáng)度因子的計(jì)算誤差在逐步減小.

圖20 全局?jǐn)?shù)學(xué)網(wǎng)格加密時(shí)裂紋強(qiáng)度因子的計(jì)算誤差Fig.20 Error of SIFs using global refinemen mesh

圖21 局部數(shù)學(xué)網(wǎng)格加密時(shí)裂紋強(qiáng)度因子的計(jì)算誤差Fig.21 Error of SIFs using local refinemen mesh

表1 采用一次AST樣條數(shù)值流形方法計(jì)算結(jié)果Table 1 SIFs results for NMMT1

表2 采用二次AST樣條數(shù)值流形方法計(jì)算結(jié)果Table 2 SIFs results for NMMT2

表3 采用三次AST樣條數(shù)值流形方法計(jì)算結(jié)果Table 3 SIFs results for NMMT3

5 結(jié)論

本文將AST樣條引入到數(shù)值流形方法中，通過提高AST樣條基函數(shù)的階次及改變相應(yīng)的數(shù)學(xué)覆蓋的形式，建立了高階數(shù)值流形方法的分析格式.AST樣條基函數(shù)具有線性無關(guān)，單位分解，局部加密等許多重要性質(zhì)，使其非常適合用于工程設(shè)計(jì)及分析.文中將AST基函數(shù)取代傳統(tǒng)的Lagrange插值函數(shù)來建立數(shù)值流形方法中的權(quán)函數(shù)，在建立高階數(shù)值流形方法的分析格式的同時(shí)，實(shí)現(xiàn)了數(shù)學(xué)網(wǎng)格的局部加密.算例結(jié)果表明：隨著AST樣條基函數(shù)階次的提高，數(shù)值流形方法的計(jì)算結(jié)果有了明顯的改善；基于AST樣條基函數(shù)的數(shù)值流形方法在保持計(jì)算精度的前提下降低了自由度的數(shù)量.

文中僅列舉了3種不同階次的AST樣條，我們?nèi)钥赏ㄟ^提高AST樣條基函數(shù)階次及改變相應(yīng)的數(shù)學(xué)覆蓋形式建立更高階次的數(shù)值流形方法.在引入AST樣條思想之后，數(shù)值流形方法中數(shù)學(xué)覆蓋的局部加密變得更加容易，更自然.

1 Shi GH.Manifold method//The 1st Int.Forum on DDA Simulation of Discontinuous Media.Bekerley,California,USA.1996：52-204

2 石根華,裴覺民.數(shù)值流形方法與非連續(xù)變形分析.北京：清華大學(xué)出版社,1997(Shi Genhua,Pei Juemin.Numerical ManifoldMethod and Discontinue Deformation Analysis.Beijing：Tsinghua University Press,1997(in Chinese))

3 Tal Y,Hatzor YH,Feng XT.An improved numerical manifold method for simulation of sequential excavation in fractured rocks.International Journal of Rock Mechanics and Mining Sciences, 2014,65：116-128

4 焦健,喬春生,徐干成.開挖模擬在數(shù)值流形方法中的實(shí)現(xiàn).巖土力學(xué),2010,31(9)：2951-2957(Jiao Jian,Qiao Chunsheng,Xu Gancheng.Simulation of excavation in numerical manifold method.Rock and Soil Mechanics,2010,31(9)：2951-2957(in Chinese))

5 朱愛軍,鄧安富,曾祥勇.數(shù)值流形方法對巖土工程開挖卸荷問題的模擬.巖土力學(xué),2006,27(2)：179-183(Zhu Aijun,Deng Anfu, Zeng Xiangyong.Numerical manifold method for simulation of excavation unloading in geotechnical engineering.Rock and Soil Mechanics,2006,27(2)：179-183(in Chinese))

6 姜清輝，王書法.錨固巖體的三維數(shù)值流形方法.巖石力學(xué)與工程學(xué)報(bào)，2006，25(3)：528-532(Jiang Qinghui,Wang Shufa.Three dimensional numerical manifold method simulation of anchor bolt supported rock mass.Chinese Journal of Rock Mechanics and Engineering，2006,25(3)：528-532(in Chinese))

7 Wei W,Jiang QH,Peng J.New rock bolt model and numerical implementation in numerical manifold method.International Journal of Geomechanics2016,E4016004

8 王水林,葛修潤.流形方法在模擬裂紋擴(kuò)展中的應(yīng)用.巖石力學(xué)與工程學(xué)報(bào),1997,16(5)：405-410(Wang Shuilin,Ge Xiuyun.Application of manifold method in simulating crack propagation.Chinese Journal of Rock Mechanics and Engineering,1997,16(5)：405-410 (in Chinese))

9 Chiou YJ,Lee YM,Tesay RJ.Mix mode fracture propagation by manifold method.International Journal of Fracture,2002,114：327-347

10 Wu ZJ,Wong LNY.Friction crack initiation and propagation analysis using the numerical manifold method.Computers&Geotechnics,2012,39：38-53

11 Yang SK,Ma GW,Ren XH,et al.Cover refinemen of numerical manifold for crack propagation simulation.Engineering Analysis with Boundary Elements,2014,43：37-49

12 Yang YT,Tang XH,Zheng H.Three-dimensional fracture propagation with numerical manifold method.Engineering Analysis with Boundary Elements,2016,72：65-77

13 Jiang QH，Deng SS，Zhou CB.Modeling unconfine seepage fl w using three-dimensionalnumericalmanifold method.Journalof Hydrodynamics，2010，22(4)：554-561

14 姜清輝,鄧書申,周創(chuàng)兵.三維高階數(shù)值流形方法研究.巖土力學(xué),2006,27(9)：1471-1474(Jiang Qinghui,Deng Shushen,Zhou Chuangbing.Study of three-dimensional high-order numerical manifold method.Rock and Soil Mechanics,2006,27(9)：1471-1474(in Chinese))

15 蘇海東,崔建華,謝小玲.高階數(shù)值流形方法的初應(yīng)力公式.計(jì)算力學(xué)學(xué)報(bào),2010,27(2)：270-274(Su Haidong,Cui Jianhua,Xie Xiaoling.Initial stress equation for high-order manifold method.Chinese Journal of Computational Mechanics,2010,27(2)：270-274(in Chinese))

16 蘇海東.固定網(wǎng)格的數(shù)值流形方法研究.力學(xué)學(xué)報(bào),2011,43(1)：169-178(Su Haidong.Study on numerical manifold method with fi ed meshs.Chinese Journal of Theoretical and Applied Mechanics,2011,43(1)：169-178(in Chinese))

17 Wang Y,Hu MS,Zhou QL,et al.A nes second-order numerical manifold method model with an efficient scheme for analyzing free surface fl w with inner drains.Applied Mathematical Modelling, 2016,40：1427-1445

18 Zheng H,Xu DD.New strategies for some issues of numerical manifold method in simulation of crack propagation.International Journal for Numerical Methods in Engineering,2014,97：986-1010

19 郭朝旭,鄭宏.高階數(shù)值流形方法中線性相關(guān)問題的研究.工程力學(xué),2012,29(12)：228-232(Guo Chaoxu,Zheng Hong.Study of linear dependence problem in high-order numerical manifold method.Engineering Mechanics,2012,29(12)：228-232(in Chinese))

20 蔡永昌,劉高揚(yáng).基于獨(dú)立覆蓋的高階數(shù)值流形方法.同濟(jì)大學(xué)學(xué)報(bào)(自然科學(xué)版),2015,43(12)：1794-1799(Cai Yongchang, Liu Gaoyang.High-order manifold method with independent covers.Journal of Tongji University:Natural Science,2015,43(12)：1794-1799(in Chinese))

21 Fan H,Zheng H,He SM,et al.A novel numerical manifold with derivative degrees of freedom and without linear dependence.Engineering Analysis with Boundary Elements,2016,64：19-37

22 Hughes TJR,Cottrell JA,Bazilevs Y.et al.Isogeometric analysis：CAD,finit elements,NURBS,exact geometry,and mesh refine ment.Computer Methods in Applied Mechanics and Engineering, 2005,194：4135-4195

23 Piegl L,Tiller W.The NURBS Book.Springer-Verlag Berlin and Heidelberg GmbH,Berlin,1997：292-295

24 Sederberg TW,Zheng J,Bakenov A,et al.T-splines and T-nurccs.ACM Transactions on Graphic,2003,22：477-484

25 Bazilevs Y,Calo VM,Cottrell JA,et al.Isogeometric analysis using T-splines.Computer Methods in Applied Mechanics and Engineering,2010,199：229-263

26 Zhang YL,Liu DX,Tan F.Numerical manifold method based on isogeometric analysis.Science China Technological Sciences,2015, 58：1520-1532

27 Li X,Zheng J,Sederberg TW,et al.On linear independence of T-splines blending functions.Computer Aided Geometric Design, 2012,29：63-76

28 Ma GW,An XM,Zhang HH,et al.Modelling complex crack problems using the numerical manifold method.International Journal of Fracture,2009,156：21-35

29 Cox MG.The numerical evaluation of B-Splines.IMA Journal of Applied Mathematics,1972,10：134-149

30 Carl DB.On calculating with B-Splines.Journal of Approximation Theory,1972,6：50-62

31 Sukumar N,Huang ZY,Prevost JH,et al.Partition of unity enrichment for bimaterial interface cracks.International Journal for Numerical Methods in Engineering,2004,59：1075-1102

HIGHER-ORDER NUMERICAL MANIFOLD METHOD BASED ON ANALYSIS-SUITABLE T-SPLINE1)

Liu Dengxue?,?,2)Zhang Youliang?Liu Gaomin?,?

?(State Key Laboratory of Geomechanics and Geotechnical Engineering，Institute of Rock and Soil Mechanics，Chinese Academy of Sciences，Wuhan430071，China)

?(University of Chinese Academy of Sciences，Beijing100049，China)

Numerical manifold method(NMM)is a very fl xible numerical method which contains and combines finit element method(FEM)and discontinuous deformation analysis(DDA).High-order numerical manifold method can be constructed by increasing the order of the weight function.This method often needs to configur the appropriate edge nodes along the element boundary,the emergence of these nodes increase the complexity of pre-processing,especially for large and complex spatial problems.On the other hand,the level of approximation of NMM can be improved by splitting the elements into smaller ones(known as h-refinement)With regard to the h-refinement a cover refinemen strategy is necessary to overcome the singularity of the stress when simulating crack propagation in NMM.One traditional solution is to refin the entire mesh which can lead to a significan decrease in the computational efficiency.In this paper analysis-suitable T-spline is introduced into NMM and regular rectangular meshes are used as the mathematical cover system. Specificall,analysis-suitable T-spline is linearly independent,forms a partition of unity,and can be locally refine which make it meet the demands of both design and analysis.The basis function of analysis-suitable T-spline is adopted as the weight function in NMM to construct high-order NMM and make the local refinemen for feasible adaptive procedure. Two numerical examples are given to demonstrate the accuracy and efficiency of the proposed method and the results show that the higher order AS T-spline based NMM shows higher accuracies when solving both continuous and discontinuous problems.Furthermore,the local mesh refinemen using analysis-suitable T-spline reduces the number of degrees of freedom while maintaining calculation accuracy at the same time.

high-order numerical manifold method，linear independence，analysis-suitable T-spline，local refinemen

O302

A doi：10.6052/0459-1879-16-217

2016-08-01收稿，2016-11-17錄用,2016-11-23網(wǎng)絡(luò)版發(fā)表.

1)國家重點(diǎn)基礎(chǔ)研究發(fā)展計(jì)劃(973計(jì)劃)(2014CB047100)和國家自然科學(xué)基金(11272330)資助項(xiàng)目.

2)劉登學(xué)，在讀博士，主要研究方向：計(jì)算巖石力學(xué).E-mail：liudengxue123@sina.cn

劉登學(xué),張友良,劉高敏.基于適合分析T樣條的高階數(shù)值流形方法.力學(xué)學(xué)報(bào),2017,49(1)：212-222

Liu Dengxue,Zhang Youliang,Liu Gaomin.High-order numerical manifold method based on analysis-suitable T-spline.Chinese Journal of Theoretical and Applied Mechanics,2017,49(1)：212-222