NEM與IEM耦合方法在地下工程中的應用研究＊

李曉龍 王復明 徐 平 李曉楠

（鄭州大學交通運輸工程系1） 鄭州 450002） （中原工學院計算機學院2） 鄭州 450007）

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

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

1 自然單元法的原理

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

1.1 自然鄰接點的尋找

在搜索積分點的自然鄰接點時，主要涉及到Delaunay三角化和Voronoi結構兩個概念.

Delaunay三角化，即利用離散結點生成Delaunay三角形，它具有兩個重要性質：（1）最大最小角性質，即在給定結點所有可能生成的三角形中，Delaunay三角形具有最大的最小內角；（2）空圓準則.通過Delaunay三角形3個頂點的外接圓里不包含其他頂點.

Voronoi結構，其定義為到點p的距離小于到其他任何結點xi的距離的集合，如圖1所示，數(shù)學表達式為

式中：d（x，p）為結點x到p的距離.

將結構作結點離散后，對于分析區(qū)域上每一個按一定規(guī)則分布的計算點p，用Delaunay三角化的空圓準則求點p周圍的自然鄰接點：如果Delaunay三角形△l的外接圓圓心cl到p的距離d小于△l的外接圓半徑R，則△l為點p的自然鄰接三角形，其結點為點p的自然鄰接點.

1.2 Laplace插值

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

圖1 點p的自然鄰接點及其一階Voronoi結構

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

從而點p的位移函數(shù)可表達為

式中：ui（i＝1，…，n）為點p周圍自然鄰結點i的結點位移；φi（x）為對應結點的插值形函數(shù).

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

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

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

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

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

式中：D為平面問題的彈性矩陣；t為面力；fb為體力；Bi為積分點上的應變矩陣.

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

2 無限單元法的原理

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

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

圖2 二維四節(jié)點無限元

坐標變換式為

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

當η≤0時，有

當η＞0時，有

式中：ξi，ηi為節(jié)點i的局部坐標值.顯然式（10）、（11）是連續(xù)的，當η→1時，對應的邊界趨于無窮遠.位移表達式為

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

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

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

3 NEM與IEM的耦合

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

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

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

4 算 例

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

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

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

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

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

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

圖6 徑向應力沿半徑的分布（NEM）

圖7 徑向應力沿半徑的分布（NEM－IEM）

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

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

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

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

表1 徑向位移值的比較

5 結 論

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

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

［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］晏 汀，沈成武.自然單元法在彈性力學中的應用［J］.武漢理工大學學報：交通科學與工程版，2006，30（6）：1052－1054.

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

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

［5］朱合華，楊寶紅，蔡永昌，等.無網(wǎng)格自然單元法在彈塑性分析中的應用［J］.巖石力學，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］周世良，胡 曉，王 江.無限元在巖土工程數(shù)值分析中的應用［J］.重慶交通學院學報，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ǎng)格方法［J］.力學進展，2000，30（3）：55－56.

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