李曉龍 王復(fù)明 徐 平 李曉楠

（鄭州大學(xué)交通運(yùn)輸工程系1） 鄭州 450002） （中原工學(xué)院計(jì)算機(jī)學(xué)院2） 鄭州 450007）

自然單元法（natural element method，NEM）是一種新型的偏微分方程數(shù)值解法，由J.Braun和 M.Sambridge于1995年提出［1］.由于該方法兼具無(wú)網(wǎng)格法和有限元兩者的優(yōu)點(diǎn)，非常適合于求解涉及到大變形、裂隙或節(jié)理擴(kuò)展、多種支護(hù)材料并且需要模擬分步施工過(guò)程的巖土及地下工程問(wèn)題，因而受到了巖土工程學(xué)者的廣泛關(guān)注，并展開(kāi)了大量的相關(guān)研究：晏汀將NEM用于處理彈性力學(xué)問(wèn)題［2］；蔡永昌研究了NEM用于地下工程計(jì)算的網(wǎng)格自動(dòng)生成技術(shù)，并將其補(bǔ)充到“曙光”軟件中，應(yīng)用于隧道、地鐵、基坑等巖土工程的線彈性分析［3－4］；朱合華實(shí)現(xiàn)了二維 NEM 基于Von－Mises屈服準(zhǔn)則的彈塑性分析［5］.

盡管自然單元法在巖土工程中的應(yīng)用研究已取得了不少成果，然而在處理隧洞或基坑等巖土工程無(wú)限域或半無(wú)限域問(wèn)題時(shí)，仍存在與有限元方法類似的問(wèn)題，即需要“人為”截取一定區(qū)域并設(shè)置相應(yīng)邊界條件，這種處理方式在理論上必然引起誤差.無(wú)限單元法（infinite element method，IEM）是20世紀(jì)70年代為解決無(wú)限域問(wèn)題而發(fā)展起來(lái)的一種數(shù)值方法，最初主要用于處理水表面波的傳播問(wèn)題［6］，后來(lái)被成功應(yīng)用于巖土工程中［7］.為彌補(bǔ)NEM在處理巖土及地下工程問(wèn)題時(shí)存在的不足，引入IEM模擬無(wú)窮遠(yuǎn)處邊界條件，與NEM相結(jié)合形成耦合分析方法（the coupling analysis of NEM and IEM，NEM－IEM），利用深埋圓形隧洞算例對(duì)算法的正確性和對(duì)計(jì)算精度的改善作用進(jìn)行檢驗(yàn)，并探討了計(jì)算范圍選取對(duì)純自然單元法和耦合方法分析結(jié)果的影響.

1 自然單元法的原理

自然單元法的主要實(shí)現(xiàn)步驟為［8］：對(duì)分析區(qū)域作結(jié)點(diǎn)離散；搜索積分點(diǎn)的自然鄰接點(diǎn)，利用自然鄰接點(diǎn)計(jì)算積分點(diǎn)的插值形函數(shù)；對(duì)計(jì)算區(qū)域積分得到結(jié)構(gòu)總體剛度陣；形成等效結(jié)點(diǎn)荷載列陣；求解線性方程組.

1.1 自然鄰接點(diǎn)的尋找

在搜索積分點(diǎn)的自然鄰接點(diǎn)時(shí)，主要涉及到Delaunay三角化和Voronoi結(jié)構(gòu)兩個(gè)概念.

Delaunay三角化，即利用離散結(jié)點(diǎn)生成Delaunay三角形，它具有兩個(gè)重要性質(zhì)：（1）最大最小角性質(zhì)，即在給定結(jié)點(diǎn)所有可能生成的三角形中，Delaunay三角形具有最大的最小內(nèi)角；（2）空?qǐng)A準(zhǔn)則.通過(guò)Delaunay三角形3個(gè)頂點(diǎn)的外接圓里不包含其他頂點(diǎn).

Voronoi結(jié)構(gòu)，其定義為到點(diǎn)p的距離小于到其他任何結(jié)點(diǎn)xi的距離的集合，如圖1所示，數(shù)學(xué)表達(dá)式為

式中：d（x，p）為結(jié)點(diǎn)x到p的距離.

將結(jié)構(gòu)作結(jié)點(diǎn)離散后，對(duì)于分析區(qū)域上每一個(gè)按一定規(guī)則分布的計(jì)算點(diǎn)p，用Delaunay三角化的空?qǐng)A準(zhǔn)則求點(diǎn)p周?chē)淖匀秽徑狱c(diǎn)：如果Delaunay三角形△l的外接圓圓心cl到p的距離d小于△l的外接圓半徑R，則△l為點(diǎn)p的自然鄰接三角形，其結(jié)點(diǎn)為點(diǎn)p的自然鄰接點(diǎn).

1.2 Laplace插值

Laplace插值，又稱 Non－Sibson插值［9］，是一種新的自然鄰點(diǎn)插值方式，與較早的Sibson插值相比，Laplace插值的優(yōu)越性主要體現(xiàn)在：（1）直接利用結(jié)點(diǎn)的1階Voronoi單胞的邊長(zhǎng)和點(diǎn)到Voronoi邊的距離計(jì)算插值基函數(shù)，構(gòu)造方式簡(jiǎn)單，避免了Sibson插值所采用的Waston面積疊加算法需反復(fù)進(jìn)行迭代的麻煩，計(jì)算量大大降低；（2）Sibson插值在凸區(qū)域的邊界是線性精確的，但對(duì)于凹區(qū)域的邊界，插值并不精確；而Laplace插值對(duì)于非凸區(qū)域邊界也是線性精確的，因此可以準(zhǔn)確地施加本質(zhì)邊界條件.

圖1 點(diǎn)p的自然鄰接點(diǎn)及其一階Voronoi結(jié)構(gòu)

在圖1所示的在二維空間中，點(diǎn)1，2，3，4為搜索得到的積分點(diǎn)p的所有自然鄰接點(diǎn).設(shè)三角形△p21，△p23，△p34，△p41的外接圓圓心分別為c21p，c23p，c34p，c41p，連接各圓心即構(gòu)成點(diǎn)p的1階Voronoi結(jié)構(gòu).分別定義s1，…，s4為點(diǎn)p的Voronoi結(jié)構(gòu)的各邊邊長(zhǎng)，h1，…，h4為各自然鄰接點(diǎn)到點(diǎn)p的距離.則對(duì)點(diǎn)p的任一自然鄰點(diǎn)i，其Laplace形函數(shù)定義為

從而點(diǎn)p的位移函數(shù)可表達(dá)為

式中：ui（i＝1，…，n）為點(diǎn)p周?chē)匀秽徑Y(jié)點(diǎn)i的結(jié)點(diǎn)位移；φi（x）為對(duì)應(yīng)結(jié)點(diǎn)的插值形函數(shù).

由式（2）知，形函數(shù)φi（x）不但滿足單位分解條件，而且與有限元形函數(shù)一樣滿足

形函數(shù)φi（x）的導(dǎo)數(shù)φi，j（x）為

顯然，對(duì)于NEM，只要分析區(qū)域的離散點(diǎn)給定，結(jié)點(diǎn)的三角化和自然鄰接點(diǎn)的搜索可利用Delaunay準(zhǔn)則由計(jì)算機(jī)自動(dòng)完成，避免了FEM在網(wǎng)格劃分和重構(gòu)等方面的困難；同時(shí)，由于其插值形函數(shù)在邊界上是線性精確的且滿足Kronecker條件，因此可以像FEM那樣施加準(zhǔn)確的邊界條件及模擬材料的不連續(xù)問(wèn)題.可以說(shuō)，NEM將無(wú)網(wǎng)格法［10］和FEM的優(yōu)點(diǎn)巧妙地結(jié)合在一起，在數(shù)值計(jì)算中具有獨(dú)特的優(yōu)勢(shì).

1.3 平衡控制方程及其數(shù)值積分

利用自然鄰接點(diǎn)采用Laplace插值得到近似位移表達(dá)式后，可按照有限元法同樣的步驟，利用最小勢(shì)能原理或Galerkin過(guò)程得到離散求解的平衡方程.對(duì)于線彈性平面問(wèn)題對(duì)總勢(shì)能取駐值得到系統(tǒng)的整體離散方程為

式中：D為平面問(wèn)題的彈性矩陣；t為面力；fb為體力；Bi為積分點(diǎn)上的應(yīng)變矩陣.

對(duì)式（6）的積分通常采用Delaunay三角形內(nèi)的高斯積分，也可用與無(wú)網(wǎng)格伽遼金法（EFG）一樣的規(guī)則矩形背景積分網(wǎng)格，或者用點(diǎn)積分［11］.

2 無(wú)限單元法的原理

無(wú)限單元法（infinite element merhod，IEM）的概念由Bettess和Zienkiewicz于1977年首次提出，其目的是為了彌補(bǔ)有限元在處理無(wú)界域問(wèn)題上存在的不足.根據(jù)是否有局部坐標(biāo)與整體坐標(biāo)的映射關(guān)系，可分為映射無(wú)限元和非映射無(wú)限元.目前在巖土工程中，映射無(wú)限元的應(yīng)用最普遍，其原理是利用映射函數(shù)將無(wú)限區(qū)域從幾何上映射為一有限區(qū)域，然后對(duì)該區(qū)域作常規(guī)有限元類似的分析.

為了便于實(shí)現(xiàn)和自然單元的耦合，本文采用二維四節(jié)點(diǎn)映射無(wú)限元，其形式如圖2所示.

圖2 二維四節(jié)點(diǎn)無(wú)限元

坐標(biāo)變換式為

式中：Ni為坐標(biāo)變換的插值函數(shù)，其形式如下.

當(dāng)η≤0時(shí)，有

當(dāng)η＞0時(shí)，有

式中：ξi，ηi為節(jié)點(diǎn)i的局部坐標(biāo)值.顯然式（10）、（11）是連續(xù)的，當(dāng)η→1時(shí)，對(duì)應(yīng)的邊界趨于無(wú)窮遠(yuǎn).位移表達(dá)式為

式中：Mi為位移插值函數(shù)，為滿足無(wú)限遠(yuǎn)處位移為零的邊界條件，可選如下形式

式中：取η≤0時(shí)Ni的表達(dá)式；f（ri／r）稱為位移衰減函數(shù)；r為衰減半徑，指積分點(diǎn)到衰減中心的距離；ri為節(jié)點(diǎn)i的衰減半徑，即節(jié)點(diǎn)i到衰減中心的距離；在衰減半徑趨于無(wú)窮大時(shí)使位移衰減函數(shù)f（ri／r）趨于零，可令

可以看出，四結(jié)點(diǎn)無(wú)限元由于在坐標(biāo)變換和位移模式中分別采用了特殊構(gòu)造且形式不同的插值函數(shù)，實(shí)現(xiàn)了從局部坐標(biāo)系中有限域到整體坐標(biāo)系下無(wú)限域的映射以及無(wú)限遠(yuǎn)處位移為0的邊界條件.

3 NEM與IEM的耦合

由于自然單元法的插值函數(shù)滿足Kronecker條件并在邊界上滿足線性插值，可以很容易與線性無(wú)限元實(shí)現(xiàn)“無(wú)縫”耦合.設(shè)求解域Ω＝ΩNEM＋ΩIEM，則近似位移表達(dá)式為

式中，ΩNEM為自然元區(qū)域；ΩIEM為無(wú)限元區(qū)域.

可見(jiàn)由于兩者的插值函數(shù)在交界處滿足線性插值，確保了位移的連續(xù)性，因此在其交界面上無(wú)需特殊處理，不必像常規(guī)無(wú)網(wǎng)格法那樣需要設(shè)置一個(gè)中間過(guò)渡區(qū)域或修改插值函數(shù)，實(shí)現(xiàn)方式十分簡(jiǎn)單.

4 算 例

某深埋地下洞室，如圖3所示，開(kāi)挖斷面為圓形，半徑為R＝3m，洞室周?chē)撵o水壓力為p0＝5MPa，假設(shè)圍巖為連續(xù)、均質(zhì)及各向同性體，彈性模量E＝2GPa，泊松比μ＝0.25.

圖3 圓形洞室及邊界條件示意圖

由于結(jié)構(gòu)對(duì)稱，取1／4區(qū)域作分析，沿環(huán)向離散為等間距的9個(gè)結(jié)點(diǎn)，徑向結(jié)點(diǎn)間距按比例逐漸增大；為比較計(jì)算范圍對(duì)計(jì)算精度的影響，按半徑為8.9，10.8，13.2，16，19.5，23.8，29m 這7種情況分別用NEM和耦合方法作線彈性分析，并與理論值作比較.

不同情況下NEM和耦合方法所得結(jié)點(diǎn)徑向位移沿半徑分布曲線分別如圖4和5所示.從中可以看出，自然單元法得到的結(jié)點(diǎn)徑向位移計(jì)算值隨計(jì)算范圍的增大逐漸趨近于理論解，而當(dāng)計(jì)算范圍減小時(shí)，兩者的偏差迅速變大；與自然單元法不同的是，耦合方法的分析結(jié)果雖然也受到計(jì)算區(qū)域選取的影響，但影響程度很小，隨著計(jì)算范圍的變化，徑向位移只發(fā)生輕微波動(dòng)，始終與理論解非常接近.

圖4 徑向位移沿半徑的分布（NEM）

圖5 徑向位移沿半徑的分布（NEM－IEM）

圖6 徑向應(yīng)力沿半徑的分布（NEM）

圖7 徑向應(yīng)力沿半徑的分布（NEM－IEM）

圖8 環(huán)向應(yīng)力沿半徑的分布（NEM）

圖9 環(huán)向應(yīng)力沿半徑的分布（NEM－IEM）

NEM和耦合方法得到的徑向和環(huán)向應(yīng)力沿半徑分布曲線分別如圖6～9所示，可以看出，兩種方法得到的徑向及環(huán)向應(yīng)力計(jì)算值與位移解具有相似特點(diǎn)：即耦合方法的分析結(jié)果受計(jì)算范圍影響小，求解精度很高，計(jì)算曲線與理論曲線幾乎重合；而純自然單元法的精度較低，當(dāng)計(jì)算區(qū)域減小時(shí)，誤差增加很快.

為了更具體地比較兩種方法的精度差別，選取3種典型情況下耦合與非耦合方法得到的結(jié)點(diǎn)徑向位移列于表1中，可以發(fā)現(xiàn)，耦合方法的計(jì)算精度遠(yuǎn)高于非耦合方法，當(dāng)計(jì)算范圍均取29m時(shí)，前者得到的隧洞內(nèi)表面徑向位移誤差僅為0.34%，而非耦合方法誤差則高達(dá)3.3%，兩者計(jì)算精度幾乎相差10倍；值得注意的是，對(duì)于耦合方法，當(dāng)計(jì)算范圍僅取8.9m時(shí)，其精度仍然遠(yuǎn)高于非耦合方法29m時(shí)的計(jì)算結(jié)果，而此時(shí)前者的結(jié)點(diǎn)總數(shù)僅為54，后者卻達(dá)到108，那么據(jù)此結(jié)果，在滿足工程要求的情況下，采用耦合方法顯然可以大大縮小結(jié)構(gòu)剛度矩陣的規(guī)模，從而減小計(jì)算工作量，這無(wú)論對(duì)于巖土工程正向分析還是反演都具有重要意義.

表1 徑向位移值的比較

5 結(jié) 論

1）在處理無(wú)限域或半無(wú)限域問(wèn)題時(shí)，采用自然元與無(wú)限元耦合分析方法，能減少純自然單元法因“人為”截取分析區(qū)域及規(guī)定邊界條件而引起的誤差，顯著提高計(jì)算結(jié)果的精度，從而進(jìn)一步提高了自然單元法處理巖土工程問(wèn)題的能力.

2）分析范圍的大小對(duì)耦合方法所得結(jié)果影響不顯著，在滿足工程要求的前提下，采用耦合方法只需選取較小的計(jì)算區(qū)域，就能得到令人滿意的結(jié)果，因而可以大大降低結(jié)構(gòu)剛度矩陣的規(guī)模，減小計(jì)算工作量，這對(duì)于復(fù)雜巖土工程尤其是其反演分析具有非常重要的意義.

［1］Braun J，Sambridge M.A numerical method for solving partial differential equations on highly irregular evolving grids［J］.Nature，1995，376：655－660.

［2］晏 汀，沈成武.自然單元法在彈性力學(xué)中的應(yīng)用［J］.武漢理工大學(xué)學(xué)報(bào)：交通科學(xué)與工程版，2006，30（6）：1052－1054.

［3］蔡永昌，朱合華.巖土工程數(shù)值計(jì)算中的無(wú)網(wǎng)格法及其全自動(dòng)布點(diǎn)技術(shù)［J］.巖土力學(xué)，2003，24（1）：21－24.

［4］蔡永昌，朱合華，夏才初.無(wú)網(wǎng)格自然鄰接點(diǎn)法及其在巖土工程數(shù)值模擬中的應(yīng)用［J］.巖石力學(xué)與工程學(xué)報(bào)，2005，24（11）：1888－1924.

［5］朱合華，楊寶紅，蔡永昌，等.無(wú)網(wǎng)格自然單元法在彈塑性分析中的應(yīng)用［J］.巖石力學(xué)，2004，25（4）：671－674.

［6］Bettess P，Zienkiewicz O C.Diffraction and refraction of surface waves using finite and infinite elements［J］.International Journal for Numerical Methods in Engineering，1977，11：1271－1290.

［7］周世良，胡 曉，王 江.無(wú)限元在巖土工程數(shù)值分析中的應(yīng)用［J］.重慶交通學(xué)院學(xué)報(bào)，2004，23（增刊）：61－64.

［8］Sukumar N，Moran B，Belytschko T.The nature element method in solid mechamics［J］.Int J Numer Meth Engng，1998，43：839－887.

［9］Sukumar N，Moran B，Semenov A Y，et al.Natural neighbor Galerkin method［J］.Int.J.Numer.Meth.Engng.，2001，50：1－27.

［10］宋康祖，陸明萬(wàn)，張 雄.固體力學(xué)中的無(wú)網(wǎng)格方法［J］.力學(xué)進(jìn)展，2000，30（3）：55－56.

［11］盧 波.自然單元法的發(fā)展及其應(yīng)用［D］.武漢：中國(guó)科學(xué)院武漢巖土力學(xué)所，2005.