用于上限有限元的非結(jié)構(gòu)化網(wǎng)格重劃加密方法研究

楊峰，陽軍生

(中南大學(xué) 土木工程學(xué)院，湖南 長沙，410075)

目前，與有限元技術(shù)相結(jié)合的極限分析上限法即上限有限元法用于巖土工程穩(wěn)定性分析已成為研究熱點(diǎn)[1-5]。該法假定巖土體為理想彈塑性材料，利用上限原理將巖土穩(wěn)定性問題轉(zhuǎn)化為線性或非線性規(guī)劃模型進(jìn)行求解[6-8]。上限有限元法具有數(shù)值求解穩(wěn)定、快速、無需假定破壞模式、適應(yīng)復(fù)雜模型等諸多優(yōu)勢。然而，即便對(duì)于二維問題，目前該法還存在如下困難：1) 上限有限元形成的線性或非線性規(guī)劃模型含有大量的決策變量和約束條件；2) 常采用三節(jié)點(diǎn)或六節(jié)點(diǎn)三角形單元，其他類型單元應(yīng)用難度大；3) 常需在相鄰單元間設(shè)速度間斷線，網(wǎng)格數(shù)據(jù)格式與共節(jié)點(diǎn)的有限元不同；4) 巖土材料失穩(wěn)破壞往往形成局部化的剪切帶和塑性區(qū)，獲得較優(yōu)上限解和破壞形態(tài)常需劃分密集的網(wǎng)格；5) 目前尚無成熟的上限有限元軟件可用。對(duì)于剪切帶和塑性區(qū)的精確模擬問題，李大鐘等[9]通過衡量屈服準(zhǔn)則殘余和等效應(yīng)變，對(duì)網(wǎng)絡(luò)進(jìn)行剖分加密以實(shí)現(xiàn)極限分析有限元網(wǎng)絡(luò)自適應(yīng)。為改進(jìn)上限有限元研究，本文作者提出1 種二維上限有限元非結(jié)構(gòu)化網(wǎng)格數(shù)據(jù)轉(zhuǎn)化方法并編制數(shù)據(jù)處理和計(jì)算程序，將MESH 2D[10]網(wǎng)格劃分工具用于上限有限元分析；提出基于耗散能密度評(píng)判參數(shù)的網(wǎng)格重劃加密技術(shù)，通過3 次網(wǎng)格重劃加密和上限有限元計(jì)算，達(dá)到不顯著增加計(jì)算規(guī)模并獲得較優(yōu)上限解的目的；實(shí)現(xiàn)用于上限有限元的三節(jié)點(diǎn)三角形單元+速度間斷線和六節(jié)點(diǎn)三角形單元無間斷線這2 種網(wǎng)格模式。最后，通過對(duì)比已有文獻(xiàn)中關(guān)于邊坡穩(wěn)定性算例，驗(yàn)證上述方法的有效性。

1 上限有限元與六節(jié)點(diǎn)三角形單元

三節(jié)點(diǎn)三角形單元+速度間斷線模式用于上限有限元的相關(guān)理論見文獻(xiàn)[6，8]。對(duì)于六節(jié)點(diǎn)三角形單元用于上限有限元，Yu 等[11]以及Makrodimopoulos等[12-13]均進(jìn)行了研究，前者研究了內(nèi)摩擦角φ 為0°且模型為線性規(guī)劃的情況，后者采用二階錐規(guī)劃模型進(jìn)行研究。本文旨在建立基于線性規(guī)劃模型和六節(jié)點(diǎn)三角形單元無間斷線模式的上限有限元并用于c - φ型土體，其基本理論[14]如下。

1.1 六節(jié)點(diǎn)三角形單元

單元內(nèi)部任一點(diǎn)的速度可表示為

式中：ui和vi分別為節(jié)點(diǎn)x 和y 向速度分量；Ni為單元形函數(shù)，若用面積坐標(biāo)L1，L2和L3表示，則

圖1 六節(jié)點(diǎn)三角形單元Fig.1 Six nodal triangular element

1.2 單元內(nèi)部約束條件

其中：

式(4)為節(jié)點(diǎn)i 上的塑性流動(dòng)約束條件，單元頂點(diǎn)對(duì)應(yīng)的節(jié)點(diǎn)(i=1，2，3)均需施加上述約束。

1.3 速度邊界條件

其中：

1.4 自重功率

在如圖1 所示坐標(biāo)系中，每個(gè)單元自重功率Pt為

其中：

γ為土體容重；A 為三角形單元面積。

1.5 單元內(nèi)部耗散功率

每個(gè)單元內(nèi)部耗散功率Pc為

其中：

c 為土體黏聚力；φ為內(nèi)摩擦角。

1.6 上限有限元的線性規(guī)劃模型

式(8)中各量均為全局變量，與式(4)～(7)中的局部變量對(duì)應(yīng)，如X1為所有節(jié)點(diǎn)速度向量，X2為所有單元塑性乘子向量。

2 非結(jié)構(gòu)化網(wǎng)格重劃加密方法

2.1 六節(jié)點(diǎn)三角形單元無間斷線模式

本文非結(jié)構(gòu)化網(wǎng)格重劃加密方法系指由前次計(jì)算結(jié)果重新劃分網(wǎng)格且變形較大的區(qū)域依網(wǎng)格尺寸權(quán)重加密，其中網(wǎng)格尺寸權(quán)重h 與單元邊長直接關(guān)聯(lián)。以下以六節(jié)點(diǎn)三角形單元為例，將耗散能密度(這里與塑性乘子一致)作為網(wǎng)格尺寸權(quán)重的評(píng)判參數(shù)進(jìn)行網(wǎng)格重劃與計(jì)算。

1) 初次計(jì)算。初次計(jì)算可獲得模型破壞區(qū)域的大致信息。以均勻網(wǎng)格劃分模型，網(wǎng)格尺寸權(quán)重h 為定值。然后轉(zhuǎn)換網(wǎng)格數(shù)據(jù)，再進(jìn)行上限有限元計(jì)算。三角形單元數(shù)目合理選取，若過大則將使計(jì)算變緩甚至不可行，若過小則計(jì)算結(jié)果精度降低。

為便于設(shè)置評(píng)判閥值Δ，將每個(gè)單元的λ～歸一化，于是，評(píng)判參數(shù)η為

網(wǎng)格重劃時(shí)，對(duì)η ≥ Δ的區(qū)域，劃分較密的網(wǎng)格，即權(quán)重h 取較小值；而對(duì)η＜Δ的區(qū)域，劃分稀疏的網(wǎng)格，即h 取較大值。之后利用該網(wǎng)格進(jìn)行上限有限元二次計(jì)算。

3) 最終計(jì)算。最終計(jì)算前，網(wǎng)格重劃所需尺寸權(quán)重h 依據(jù)二次計(jì)算所得塑性乘子值，權(quán)重計(jì)算點(diǎn)仍為三角形單元形心。評(píng)判參數(shù)η的計(jì)算與前面的相同。所不同的是：此次網(wǎng)格重劃時(shí)，由0≤η≤1 取值范圍設(shè)多級(jí)閥值Δ1，Δ2，…，Δn，不同的閥值區(qū)間遞減設(shè)置權(quán)重(h1，h2，…，hn)，以此劃分網(wǎng)格并進(jìn)行上限有限元最終計(jì)算。

2.2 三節(jié)點(diǎn)三角形單元+間斷線模式

對(duì)于三節(jié)點(diǎn)三角形單元+間斷線模式，其網(wǎng)格重劃加密方法和步驟與六節(jié)點(diǎn)三角形單元的一致，不同的是評(píng)判參數(shù)η需考慮速度間斷線耗散能。

仍選取耗散能密度作為網(wǎng)格權(quán)重h 的評(píng)判參數(shù)，計(jì)算點(diǎn)為三角形單元形心，將單元相鄰間斷線的耗散能的一半疊加到該三角形單元耗散能上。

每個(gè)三角形單元耗散能密度 ρe為[8]

該單元相鄰間斷線耗散能密度為[8]

設(shè)ρd1，ρd2和ρd3分別為該單元三條邊所在間斷線耗散能密度，將其一半疊加到該單元可得

將模型耗散能密度ρ歸一化，網(wǎng)格權(quán)重h 的評(píng)判參數(shù)η為

2.3 MESH 2D 網(wǎng)格劃分與數(shù)據(jù)轉(zhuǎn)換

MESH 2D 是Engwirda[10]編制的一個(gè)二維非結(jié)構(gòu)網(wǎng)格劃分工具，可嵌入MATLAB。網(wǎng)格劃分所得數(shù)據(jù)為節(jié)點(diǎn)編號(hào)與坐標(biāo)、三角形單元編號(hào)及內(nèi)部節(jié)點(diǎn)編號(hào)。對(duì)于六節(jié)點(diǎn)三角形單元無間斷線模式，上限有限元所需的網(wǎng)格數(shù)據(jù)轉(zhuǎn)化思路如下：

1) 原型每個(gè)節(jié)點(diǎn)含x 和y 向速度分量。

2) 獲取單元的邊編號(hào)，每條邊中點(diǎn)附加1 個(gè)節(jié)點(diǎn)，該節(jié)點(diǎn)含x 和y 向速度分量。

3) 塑性乘子歸屬于每個(gè)單元的3 個(gè)頂點(diǎn)對(duì)應(yīng)的節(jié)點(diǎn)，分屬不同單元的同位置頂點(diǎn)節(jié)點(diǎn)處塑性乘子不共有。

4) 單元節(jié)點(diǎn)及附加量順序按圖1 所示排序。

5) 依模型邊界標(biāo)識(shí)搜尋并記錄邊界節(jié)點(diǎn)、單元及局部坐標(biāo)法向等信息，以便施加邊界條件。

對(duì)于三節(jié)點(diǎn)三角形單元+間斷線模式，上限有限元所需的網(wǎng)格數(shù)據(jù)轉(zhuǎn)化思路為：

1) 速度分量歸屬于單元節(jié)點(diǎn)，相鄰單元同位置處的節(jié)點(diǎn)不共有，塑性乘子歸屬于單元。

2) 所有單元的邊依邊界標(biāo)識(shí)分為速度間斷線和邊界邊2 種。速度間斷線附加輔助變量，邊界邊記錄所屬單元、所含節(jié)點(diǎn)及局部坐標(biāo)法向等信息。

3) 搜索、判斷并建立速度間斷線兩側(cè)單元編號(hào)、間斷線節(jié)點(diǎn)在兩側(cè)單元節(jié)點(diǎn)的位置等信息，以便施加速度間斷線塑性流動(dòng)約束。

2.4 上限有限元計(jì)算及網(wǎng)格重劃加密流程

為便于編程，將上限有限元線性規(guī)劃模型決策變量之一(節(jié)點(diǎn)速度分量)轉(zhuǎn)化為非負(fù)量。上限有限元主程序采用MATLAB 編制，網(wǎng)格劃分利用MESH 2D[10]，線性規(guī)劃模型采用SeDuMi[15]求解。

當(dāng)采用網(wǎng)格重劃加密和六節(jié)點(diǎn)三角形單元無間斷線模式時(shí)，上限有限元的計(jì)算流程如圖2 所示。

圖2 上限有限元計(jì)算流程圖Fig.2 Flow chart for finite element upper bound solution

3 算例驗(yàn)證與討論

選取文獻(xiàn)[12-13]中的邊坡穩(wěn)定性算例進(jìn)行上限有限元計(jì)算驗(yàn)證。該算例求解不同內(nèi)摩擦角φ 對(duì)應(yīng)的邊坡穩(wěn)定性系數(shù)Ns。模型的初次網(wǎng)格劃分、邊界條件、幾何尺寸如圖3 所示。邊坡坡度為70°，左右和下邊界均約束2 個(gè)方向的速度，即u=0 m/s，v=0 m/s。邊坡穩(wěn)定性系數(shù)Ns(無量綱)定義為

圖3 邊坡穩(wěn)定性計(jì)算模型網(wǎng)格劃分及邊界條件Fig.3 Mesh and boundary condition of calculating model with slope stability

本例邊坡坡高H=25 m，黏結(jié)力c=1 kPa，內(nèi)摩擦角φ=35°。以Ns為目標(biāo)函數(shù)，令單位容重(γ=1)時(shí)自重功率為1，利用式(6)施加自重約束條件。摩爾-庫侖屈服準(zhǔn)則線性化參數(shù)p 取32。

采用網(wǎng)格重劃加密的上限有限元法，按六節(jié)點(diǎn)三角形單元無間斷線和三節(jié)點(diǎn)三角形單元+間斷線2 種模式所得邊坡穩(wěn)定性系數(shù)Ns分別如表1 和表2 所示。表1 和表2 中也列出了文獻(xiàn)[12]和[16]中的計(jì)算結(jié)果，并以文獻(xiàn)[16]中的計(jì)算結(jié)果作為評(píng)判相對(duì)誤差的參考值。

從表1 可看出：采用六節(jié)點(diǎn)三角形單元進(jìn)行上限有限元計(jì)算，初次和二次計(jì)算結(jié)果Ns的相對(duì)誤差分別為14.90%和5.50%，Ns值均與文獻(xiàn)[12]中的相應(yīng)值14.97 和14.45 接近，不過文獻(xiàn)[12]中的第2 種網(wǎng)格單元數(shù)約為本文二次計(jì)算網(wǎng)格單元數(shù)的2.5 倍。此外，本文最終計(jì)算的單元數(shù)為4 006，Ns相對(duì)誤差為1.86%，而文獻(xiàn)[12]中采用的網(wǎng)格單元數(shù)為28 864，Ns相對(duì)誤差為2.38%。本文結(jié)果相對(duì)誤差小且模型網(wǎng)格單元數(shù)少，這是由于采用了網(wǎng)格重劃加密方法。

表1 邊坡穩(wěn)定性系數(shù)Ns 上限有限元計(jì)算結(jié)果(六節(jié)點(diǎn)三角形單元無間斷線)Table 1 Results of slope stability factor Ns by finite element upper bound solution with six nodal triangular elements

表2 邊坡穩(wěn)定性系數(shù)Ns 上限有限元計(jì)算結(jié)果(三節(jié)點(diǎn)三角形單元+間斷線)Table 2 Results of slope stability factor Ns by finite element upper bound solution with three nodal triangular elements and velocity discontinuities

從表2 可以看出：采用三節(jié)點(diǎn)三角形單元+間斷線進(jìn)行上限有限元計(jì)算，二次和最終計(jì)算所得Ns的相對(duì)誤差分別為11.75%和8.74%，均小于文獻(xiàn)[12]中的13.13%和11.54%，而本文采用的網(wǎng)格單元數(shù)少，特別是最終計(jì)算時(shí)采用的網(wǎng)格單元數(shù)更少。這說明基于三節(jié)點(diǎn)三角形單元+間斷線模式的網(wǎng)格重劃加密方法對(duì)降低計(jì)算規(guī)模也具有優(yōu)勢。

對(duì)比表1 和表2 最終計(jì)算的網(wǎng)格參數(shù)可知：表1中單元、節(jié)點(diǎn)和約束條件數(shù)目均比表2 中的少，僅變量數(shù)目比表2 中的大，這說明二者計(jì)算規(guī)模差異不大。然而，所得Ns相對(duì)誤差分別為1.86%和8.74%，說明采用網(wǎng)格重劃加密方法時(shí)，六節(jié)點(diǎn)三角形單元無間斷線模式比三節(jié)點(diǎn)三角形單元+間斷線模式精度更高。

此外，本文認(rèn)為對(duì)該邊坡算例，文獻(xiàn)[16]采用基于假定剛性滑體破壞模式所得上限解更優(yōu)，這是由于其模型和破壞模式均簡單明確，且形成明顯剪切帶。對(duì)于難以預(yù)先找到合理破壞模式的復(fù)雜模型，假定剛體滑塊破壞模式的上限法精度優(yōu)勢將喪失。

以下對(duì)初次、二次和最終三階段的網(wǎng)格塑性變形及最終所得耗散能密度分布進(jìn)行討論。圖4 所示為采用六節(jié)點(diǎn)三角形單元無間斷模式所得網(wǎng)格塑性變形圖(僅顯示了模型局部范圍)。初次計(jì)算采用均勻網(wǎng)格(見圖4(a))。二次計(jì)算對(duì)耗散能較大的區(qū)域網(wǎng)格加密，加密區(qū)域網(wǎng)格相對(duì)均勻，其網(wǎng)格變形在坡腳處較劇烈(見圖4(b))。最終計(jì)算采用的網(wǎng)格在坡腳附近相對(duì)密集，向上至坡頂形成1 條密集的帶狀區(qū)域，稍遠(yuǎn)處網(wǎng)格稀疏(見圖4(c))。圖5 所示為最終計(jì)算所得耗散能密度分布圖。耗散能密度以塑性乘子替換并進(jìn)行歸一化，使之成為無量綱參數(shù)，在坡腳處其值最大，向上呈圓滑、狹窄的帶狀直至地表。耗散能密度所反映的帶狀區(qū)域即為邊坡滑動(dòng)面。

圖4 不同計(jì)算階段網(wǎng)格塑性變形圖(六節(jié)點(diǎn)三角形單元)Fig.4 Plastic deformation of mesh for various calculating stages with six nodal triangular elements

圖5 最終計(jì)算所得耗散能密度分布圖(六節(jié)點(diǎn)三角形單元)Fig.5 Distribution of dissipated energy density for last calculating stage with six nodal triangular elements

圖6 不同計(jì)算階段網(wǎng)格塑性變形圖(三節(jié)點(diǎn)三角形單元+間斷線)Fig.6 Plastic deformation of mesh for various calculating stage with three nodal triangular elements and velocity discontinuities

圖6 所示為采用三節(jié)點(diǎn)三角形單元+間斷線模式所得網(wǎng)格塑性變形圖。其初次、二次和最終網(wǎng)格塑性變形圖與圖4 所示的基本類似。不同之處在于由于存在間斷線，坡腳處單元之間錯(cuò)動(dòng)明顯，顯示了間斷線模擬劇烈變形區(qū)的優(yōu)勢。然而，三節(jié)點(diǎn)單元模擬劇烈變形區(qū)效果差，最終計(jì)算前網(wǎng)格重劃加密形成的密集區(qū)較寬，計(jì)算結(jié)果精度也較低。圖7 所示為采用該模式最終計(jì)算所得耗散能密度分布圖，耗散能密度同樣進(jìn)行了無量綱處理，該圖顯示的耗散能密度分布形狀與圖5 所示的類似，不過其表示的剪切帶更寬泛。

圖7 最終計(jì)算所得耗散能密度分布圖(三節(jié)點(diǎn)三角形單元+間斷線)Fig.7 Distribution of dissipated energy density for last calculating stage with three nodal triangular elements and velocity discontinuities

4 結(jié)論

1) 通過網(wǎng)格單元數(shù)據(jù)轉(zhuǎn)換，可將MESH 2D 作為上限有限元的網(wǎng)格非結(jié)構(gòu)化網(wǎng)格劃分工具，從而減小復(fù)雜模型建模的工作量。

2) 基于耗散能密度為評(píng)判參數(shù)的網(wǎng)格重劃加密方法，能依據(jù)前次結(jié)果信息重劃網(wǎng)格對(duì)應(yīng)變率較大的區(qū)域網(wǎng)格加密，達(dá)到降低計(jì)算規(guī)模且不影響計(jì)算精度的目的。

3) 配合網(wǎng)格重劃加密方法，采用六節(jié)點(diǎn)三角形單元無間斷線模式比三節(jié)點(diǎn)三角形單元+間斷線模式更能有效提高計(jì)算精度并降低計(jì)算規(guī)模。

[1] Augarde C E, Lyamin A V, Sloan S W.Stability of an undrained plane strain heading revisited[J]. Computers and Geotechnics,2003,30：419-430.

[2] 王均星, 王漢輝, 吳雅峰. 土坡穩(wěn)定的有限元塑性極限分析上限法研究[J]. 巖石力學(xué)與工程學(xué)報(bào), 2004, 23(11)：1867-1873.WANG Junxing, WANG Hanhui, WU Yafeng. Stability analysis of soil slope by finite element method with plastic limit upper bound[J]. Chinese Journal of Rock Mechanics and Engineering,2004,23(11)：1867-1873.

[3] 殷建華, 陳健, 李焯芬. 巖土邊坡穩(wěn)定性的剛體有限元上限分析法[J]. 巖石力學(xué)與工程學(xué)報(bào),2004,23(6)：898-905.YIN Jianhua, CHEN Jian, LI Zhuofen. Upper bound limit analysis of stability of rock and soil slopes using rigid finite elements[J]. Chinese Journal of Rock Mechanics and Engineering,2004,23(6)：898-905.

[4] 楊峰, 陽軍生, 張學(xué)民, 等. 黏土不排水條件下淺埋隧道穩(wěn)定性上限有限元分析[J]. 巖石力學(xué)與工程學(xué)報(bào), 2010, 29(S2)：3952-3959.YANG Feng, YANG Junsheng, ZHANG Xuemin, et al. Finite element analysis of upper bound solution of shallow-buried tunnel stability in undrained clay[J]. Chinese Journal of Rock Mechanics and Engineering,2010,29(S2)：3952-3959.

[5] YANG Feng,YANG Junsheng. Stability of shallow tunnel using rigid blocks and finite element upper bound solutions[J].International Journal of Geomechanics,2010,10(6)：242-247.

[6] Sloan S W.Upper bound limit analysis using finite elements and linear programming[J]. International Journal for Nummerical and Analytical Methods in Geomechanics, 1989, 13(3)：263-282.

[7] Sloan S W, Kleeman P W. Upper bound limit analysis using discontinuous velocity fields[J]. Computer Methods in Applied Mechanics and Engineering,1995,127(1)：293-314.

[8] 楊峰, 陽軍生, 張學(xué)民. 基于線性規(guī)劃模型的極限分析上限有限元的實(shí)現(xiàn)[J]. 巖土力學(xué),2011,32(3)：914-921.YANGFeng, YANGJunsheng,ZHANGXuemin.Implementation of finite element upper bound solution of limit analysis based on linear programming model[J]. Rock and Soil Mechanics,2011,32(3)：914-921.

[9] 李大鐘, 鄭榕明, 王金安, 等. 自適應(yīng)有限元極限分析及巖土工程中的應(yīng)用[J]. 巖土工程學(xué)報(bào),2013,35(5)：922-929.LI Dazhong, ZHENG Rongming, WANG Jinan, et al.Application of finite-element-based limit analysis with mesh adaptation in geotechnical engineering[J]. Chinese Journal of Geotechnical Engineering,2013,25(5)：922-929.

[10] Engwirda D. MESHZD-automatic mesh generation[EB/OL].http：//www. mathworks. com/ matlabcentral/ fileexchange/25555-mesh2d-automatic-mesh-generation,2009-10-12.

[11] Yu H S,Sloan S W,Kleeman P W.A quadratic element for upper bound limit analysis[J]. Engineering Computations, 1994, 11(3)：195-212.

[12] Makrodimopoulos A, Martin C M. Upper bound limit analysis using simplex strain elements and second-order cone programming[J]. International Journal for Numerical and Analytical Methods in Geomechanics,2007,31(6)：835-865.

[13] Makrodimopoulos A, Martin C M. Upper bound limit analysis using discontinuousquadraticdisplacementfields[J].Communications in Numerical Methods in Engineering, 2008,24(11)：911-927.

[14] 楊峰, 陽軍生, 李昌友, 等. 基于六節(jié)點(diǎn)三角形單元和線性規(guī)劃模型的上限有限元研究[J]. 巖石力學(xué)與工程學(xué)報(bào), 2012,31(12)：2556-2563.YANG Feng,YANG Junsheng,LI Changyou, et al.Investigation of finite element upper bound solution based on six nodal triangular elements and linear programming modal[J]. Chinese Journal of Rock Mechanics and Engineering, 2012, 31(12)：2556-2563.

[15] Sturm J F. Using SeDuMi 1.02, a MATLAB toolbox for optimization over symmetric cones[J]. Optimization Methods and Software,1999,11(12)：625-653.

[16] Chen W F. Limit analysis and soil mechanics[M]. New York：Elsevier Scientific Publishing Company,1975：254-255.