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

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

目前，結(jié)合有限技術(shù)的限制分析上限法，即限制分析法用于巖石工程的穩(wěn)定性分析，已成為研究的重點(diǎn)。1上限偏差和六個節(jié)點(diǎn)的三角形單元三節(jié)點(diǎn)三角形單元+速度間斷線模式用于上限有限元的相關(guān)理論見文獻(xiàn)[6,8]。對于六節(jié)點(diǎn)三角形單元用于上限有限元,Yu等1.1元每邊節(jié)點(diǎn)節(jié)點(diǎn)六節(jié)點(diǎn)三角形單元如圖1所示,六節(jié)點(diǎn)三角形單元每邊中點(diǎn)增設(shè)1個節(jié)點(diǎn)。每個節(jié)點(diǎn)含2個速度分量,單元頂點(diǎn)對應(yīng)的節(jié)點(diǎn)1,2和3還分別包含p個塑性乘子單元內(nèi)部任一點(diǎn)的速度可表示為式中:u1.2塑性流動約束條件對摩爾-庫侖屈服函數(shù)線性化時,單元內(nèi)部應(yīng)服從相關(guān)聯(lián)流動法則,于是,節(jié)點(diǎn)i(i=1,2,3)處塑性應(yīng)變率其中:式(4)為節(jié)點(diǎn)i上的塑性流動約束條件,單元頂點(diǎn)對應(yīng)的節(jié)點(diǎn)(i=1,2,3)均需施加上述約束。1.3速度邊界條件設(shè)節(jié)點(diǎn)i位于與x軸夾角為θ的邊界上,其切向和法向速度分別為u和v,則其中:1.4重量損失在如圖1所示坐標(biāo)系中,每個單元自重功率P其中:γ為土體容重;A為三角形單元面積。1.5單元內(nèi)部的耗散能量每個單元內(nèi)部耗散功率P其中:c為土體黏聚力;uf066為內(nèi)摩擦角。1.6生成目標(biāo)函數(shù)根據(jù)上限定理,由耗散能、外力功率和虛功率平衡方程即可建立目標(biāo)函數(shù)。其決策變量為節(jié)點(diǎn)速度分量u式(8)中各量均為全局變量,與式(4)~(7)中的局部變量對應(yīng),如X2非結(jié)構(gòu)網(wǎng)格重復(fù)計加密方法2.1網(wǎng)格重劃與計算本文非結(jié)構(gòu)化網(wǎng)格重劃加密方法系指由前次計算結(jié)果重新劃分網(wǎng)格且變形較大的區(qū)域依網(wǎng)格尺寸權(quán)重加密,其中網(wǎng)格尺寸權(quán)重h與單元邊長直接關(guān)聯(lián)。以下以六節(jié)點(diǎn)三角形單元為例,將耗散能密度(這里與塑性乘子一致)作為網(wǎng)格尺寸權(quán)重的評判參數(shù)進(jìn)行網(wǎng)格重劃與計算。1)初次計算。初次計算可獲得模型破壞區(qū)域的大致信息。以均勻網(wǎng)格劃分模型,網(wǎng)格尺寸權(quán)重h為定值。然后轉(zhuǎn)換網(wǎng)格數(shù)據(jù),再進(jìn)行上限有限元計算。三角形單元數(shù)目合理選取,若過大則將使計算變緩甚至不可行,若過小則計算結(jié)果精度降低。2)二次計算。二次計算前的網(wǎng)格重劃需設(shè)置網(wǎng)格權(quán)重h,以特征點(diǎn)塑性乘子作為權(quán)重h評判參數(shù)。特征點(diǎn)為三角形單元形心,該點(diǎn)塑性乘子uf06c(4)為為便于設(shè)置評判閥值(35),將每個單元的uf06c(4)歸一化,于是,評判參數(shù)η為式中:網(wǎng)格重劃時,對3)最終計算。最終計算前,網(wǎng)格重劃所需尺寸權(quán)重h依據(jù)二次計算所得塑性乘子值,權(quán)重計算點(diǎn)仍為三角形單元形心。評判參數(shù)η的計算與前面的相同。所不同的是:此次網(wǎng)格重劃時,由0≤η≤1取值范圍設(shè)多級閥值Δ2.2速度間斷線耗散能密度h對于三節(jié)點(diǎn)三角形單元+間斷線模式,其網(wǎng)格重劃加密方法和步驟與六節(jié)點(diǎn)三角形單元的一致,不同的是評判參數(shù)η需考慮速度間斷線耗散能。仍選取耗散能密度作為網(wǎng)格權(quán)重h的評判參數(shù),計算點(diǎn)為三角形單元形心,將單元相鄰間斷線的耗散能的一半疊加到該三角形單元耗散能上。每個三角形單元耗散能密度uf072該單元相鄰間斷線耗散能密度為式中:u設(shè)uf072將模型耗散能密度ρ歸一化,網(wǎng)格權(quán)重h的評判參數(shù)η為2.3節(jié)點(diǎn)三角形單元+間斷線模式MESH2D是Engwirda1)原型每個節(jié)點(diǎn)含x和y向速度分量。2)獲取單元的邊編號,每條邊中點(diǎn)附加1個節(jié)點(diǎn),該節(jié)點(diǎn)含x和y向速度分量。3)塑性乘子歸屬于每個單元的3個頂點(diǎn)對應(yīng)的節(jié)點(diǎn),分屬不同單元的同位置頂點(diǎn)節(jié)點(diǎn)處塑性乘子不共有。4)單元節(jié)點(diǎn)及附加量順序按圖1所示排序。5)依模型邊界標(biāo)識搜尋并記錄邊界節(jié)點(diǎn)、單元及局部坐標(biāo)法向等信息,以便施加邊界條件。對于三節(jié)點(diǎn)三角形單元+間斷線模式,上限有限元所需的網(wǎng)格數(shù)據(jù)轉(zhuǎn)化思路為:1)速度分量歸屬于單元節(jié)點(diǎn),相鄰單元同位置處的節(jié)點(diǎn)不共有,塑性乘子歸屬于單元。2)所有單元的邊依邊界標(biāo)識分為速度間斷線和邊界邊2種。速度間斷線附加輔助變量,邊界邊記錄所屬單元、所含節(jié)點(diǎn)及局部坐標(biāo)法向等信息。3)搜索、判斷并建立速度間斷線兩側(cè)單元編號、間斷線節(jié)點(diǎn)在兩側(cè)單元節(jié)點(diǎn)的位置等信息,以便施加速度間斷線塑性流動約束。2.4速度分量非負(fù)量轉(zhuǎn)化為非負(fù)量b轉(zhuǎn)化為便于編程,將上限有限元線性規(guī)劃模型決策變量之一(節(jié)點(diǎn)速度分量)轉(zhuǎn)化為非負(fù)量。上限有限元主程序采用MATLAB編制,網(wǎng)格劃分利用MESH2D當(dāng)采用網(wǎng)格重劃加密和六節(jié)點(diǎn)三角形單元無間斷線模式時,上限有限元的計算流程如圖2所示。3算例2:某選取文獻(xiàn)[12-13]中的邊坡穩(wěn)定性算例進(jìn)行上限有限元計算驗(yàn)證。該算例求解不同內(nèi)摩擦角uf066對應(yīng)的邊坡穩(wěn)定性系數(shù)N本例邊坡坡高H=25m,黏結(jié)力c=1kPa,內(nèi)摩擦角采用網(wǎng)格重劃加密的上限有限元法,按六節(jié)點(diǎn)三角形單元無間斷線和三節(jié)點(diǎn)三角形單元+間斷線2種模式所得邊坡穩(wěn)定性系數(shù)N從表1可看出:采用六節(jié)點(diǎn)三角形單元進(jìn)行上限有限元計算,初次和二次計算結(jié)果N從表2可以看出:采用三節(jié)點(diǎn)三角形單元+間斷線進(jìn)行上限有限元計算,二次和最終計算所得N對比表1和表2最終計算的網(wǎng)格參數(shù)可知:表1中單元、節(jié)點(diǎn)和約束條件數(shù)目均比表2中的少,僅變量數(shù)目比表2中的大,這說明二者計算規(guī)模差異不大。然而,所得N此外,本文認(rèn)為對該邊坡算例,文獻(xiàn)[16]采用基于假定剛性滑體破壞模式所得上限解更優(yōu),這是由于其模型和破壞模式均簡單明確,且形成明顯剪切帶。對于難以預(yù)先找到合理破壞模式的復(fù)雜模型,假定剛體滑塊破壞模式的上限法精度優(yōu)勢將喪失。以下對初次、二次和最終三階段的網(wǎng)格塑性變形及最終所得耗散能密度分布進(jìn)行討論。圖4所示為采用六節(jié)點(diǎn)三角形單元無間斷模式所得網(wǎng)格塑性變形圖(僅顯示了模型局部范圍)。初次計算采用均勻網(wǎng)格(見圖4(a))。二次計算對耗散能較大的區(qū)域網(wǎng)格加密,加密區(qū)域網(wǎng)格相對均勻,其網(wǎng)格變形在坡腳處較劇烈(見圖4(b))。最終計算采用的網(wǎng)格在坡腳附近相對密集,向上至坡頂形成1條密集的帶狀區(qū)域,稍遠(yuǎn)處網(wǎng)格稀疏(見圖4(c))。圖5所示為最終計算所得耗散能密度分布圖。耗散能密度以塑性乘子替換并進(jìn)行歸一化,使之成為無量綱參數(shù),在坡腳處其值最大,向上呈圓滑、狹窄的帶狀直至地表。耗散能密度所反映的帶狀區(qū)域即為邊坡滑動面。圖6所示為采用三節(jié)點(diǎn)三角形單元+間斷線模式所得網(wǎng)格塑性變形圖。其初次、二次和最終網(wǎng)格塑性變形圖與圖4所示的基本類似。不同之處在于由于存在間斷線,坡腳處單元之間錯動明顯,顯示了間斷線模擬劇烈變形區(qū)的優(yōu)勢。然而,三節(jié)點(diǎn)單元模擬劇烈變形區(qū)效果差,最終計算前網(wǎng)格重劃加密形成的密集區(qū)較寬,計算結(jié)果精度也較低。圖7所示為采用該模式最終計算所得耗散能密度分布圖,耗散能密度同樣進(jìn)行了無量綱處理,該圖顯示的耗散能密度分布形狀與圖5所示的類似,不過其表示的剪切帶更寬泛。4網(wǎng)格重劃加密的數(shù)學(xué)模型1)通過網(wǎng)格單元數(shù)據(jù)轉(zhuǎn)換,可將MESH2D作為上限有限元的網(wǎng)格非結(jié)構(gòu)化網(wǎng)格劃分工具,