基于平衡顆粒法生成地質(zhì)體復(fù)雜節(jié)理面網(wǎng)格的方法研究

1、基于平衡顆粒法生成地質(zhì)體復(fù)雜節(jié)理面網(wǎng)格的方法研究 2012 2012 9 16 -18 顆粒材料計(jì)算力學(xué)會(huì)議 年 月 日 日National Conference on Computational Mechanics of Granular Materials CMGM-2012 湖南 ? 張家界 論文編號(hào):CMGM2012-D06ISBN :978-7-5611-7240-7 基于平衡顆粒法生成地質(zhì)體復(fù)雜節(jié)理面 *網(wǎng)格的方法研究 1? 2 2段文杰 李世海 馮春1 , 300401 河北工業(yè)大學(xué) 天津2 中國(guó)科學(xué)院力學(xué)研究所, 北京 100190 摘要: 宏觀地質(zhì)體中錯(cuò)綜復(fù)雜的節(jié)理拓?fù)潢P(guān)系給

2、節(jié)理面網(wǎng)格的劃分帶來了困難,進(jìn)而制約了對(duì)含復(fù)雜節(jié)理地質(zhì)體的力學(xué)效應(yīng)研究。文章提出了一種利用顆粒生成節(jié)理面網(wǎng)格節(jié)點(diǎn)的方法? 平衡顆粒法。該方法首先將復(fù)雜節(jié)理系統(tǒng)按照節(jié)理交線進(jìn)行分割,而后通過節(jié)理面內(nèi)顆粒的隨機(jī)排布產(chǎn)生密實(shí)、微嵌入的顆粒系統(tǒng),進(jìn)而采用全量法計(jì)算顆粒間的接觸力,同時(shí)引入指向節(jié)理面中心的牽引勢(shì)函數(shù),在粘性阻尼的控制下,通過顯式迭代,逐漸實(shí)現(xiàn)接觸力產(chǎn)生的斥力與勢(shì)函數(shù)產(chǎn)生的引力之間的平衡。算例表明,基于平衡顆粒法產(chǎn)生的節(jié)理網(wǎng)格節(jié)點(diǎn)排布均勻,在節(jié)理交線處未出現(xiàn)奇異節(jié)點(diǎn),表明了該方法的正確性。 關(guān)鍵詞:節(jié)理, 網(wǎng)格劃分, 顆粒1 前言節(jié)理裂隙是地質(zhì)體的特有屬性, 其發(fā)育狀況、 空間展布特征、

3、節(jié)理間的拓?fù)潢P(guān)系、 節(jié)理面1, 2內(nèi)的充填膠結(jié)情況,直接影響著地質(zhì)體的宏觀滲透特性及力學(xué)性能 。數(shù)值分析是研究節(jié)理裂隙系統(tǒng)復(fù) 雜物理力學(xué) 特性的有效 手段之一, 而節(jié)理單元 的網(wǎng)格剖分 是進(jìn)行數(shù)值 分析的第一步。 有限元網(wǎng) 格 生成技術(shù) 可 以按所生 成 的單元類 型 將它們劃 分 為生成結(jié) 構(gòu) 化網(wǎng)格的 方 法 和 生3成非結(jié)構(gòu)化網(wǎng)格的方法 。 生成結(jié)構(gòu)化網(wǎng)格的方法在邊界上網(wǎng)格質(zhì)量較低, 不適用于節(jié)理面網(wǎng)格的生成。 生成非結(jié)構(gòu) 化網(wǎng)格的方 法主要有布 點(diǎn)及三角化 方法,拓?fù)?分解法,幾 何分解法,基于柵格法 幾種。目前 已經(jīng)有利用 野外露頭上 采集到的裂 隙的觀測(cè)數(shù) 據(jù),通過計(jì) 算機(jī)處理昀4

4、終形成三維裂隙網(wǎng)絡(luò)的人工幾何模型的方法 , 也已經(jīng)有多種在節(jié)理網(wǎng)絡(luò)圖上生成有限元網(wǎng)格5, 6 7, 8的方法 以 及生成節(jié)理巖體網(wǎng)格的方法 。其 中文獻(xiàn)5 使 用的是 AFMAdvancing Front Method 法 ,首先將具有節(jié)理網(wǎng)絡(luò)的計(jì)算域轉(zhuǎn)化成單連通域, 然后在單 連通域內(nèi)按 AFM 法生成三角形單 元。這種方 法在節(jié)理面 附近的網(wǎng)格 質(zhì)量很差, 而且目前還 不能用于三 維問題。文獻(xiàn)7 采用的 是流形 切割 的方法 ,將 三維有 限元 網(wǎng)格覆 蓋到 物理模 型上 ,然后 對(duì)覆 蓋網(wǎng)格 進(jìn) 行剖分,形成昀后的網(wǎng)格。這種方法可以用于三維問題,但是切割出的網(wǎng)格質(zhì)量無法控制。 9Bagi

5、 提出了 用大小不同的顆粒填滿一個(gè)空間的方法 ,Jianfei Liu 等人提 出了用球形填充10Sphere packing 來構(gòu)造有 限元網(wǎng)格的方法 。 受這兩種方法啟發(fā), 本文提出了一種利用顆粒生*資助項(xiàng)目: 國(guó)家 自然科學(xué)基金11172063, 國(guó)家重點(diǎn)基礎(chǔ)研究發(fā)展計(jì)劃2010CB731502 E-mail: duanwenjiepkugmail - 340 -段文杰等: 基于平衡顆粒法生成地質(zhì)體復(fù)雜節(jié)理面網(wǎng)格的方法研究 成節(jié)理面網(wǎng) 格節(jié)點(diǎn)的方 法?平衡 顆粒法。該 方法首先將 復(fù)雜節(jié)理系 統(tǒng)按照節(jié)理 交線進(jìn)行分割,而后通 過節(jié)理面內(nèi) 顆粒的隨機(jī) 排布產(chǎn)生密 實(shí)、微嵌入 的顆粒系統(tǒng) ,

6、進(jìn)而采用 全量法計(jì)算顆粒間的接 觸力,同時(shí) 引入指向節(jié) 理面中心的 牽引勢(shì)函數(shù) ,在粘性阻 尼的控制下 ,通過顯式迭代,逐漸實(shí)現(xiàn)接觸力產(chǎn)生的斥力與勢(shì)函數(shù)產(chǎn)生的引力之間的平衡。 2 平衡顆粒法基本原理 節(jié)理面之間的空間關(guān)系可能很復(fù)雜, 如果針對(duì)每一種情況如何劃分網(wǎng)格進(jìn)行規(guī)定工作量難以估計(jì)。用 顆粒組成的 面代替幾何 的節(jié)理面, 然后讓顆粒 自行尋找能 量低的狀態(tài) 就成為一種可行的思路。這就需要把節(jié)理面盡可能填滿顆粒。 2.1 節(jié)理面分割 如果直接對(duì)每個(gè)節(jié)理面進(jìn)行顆粒填充, 那么在兩個(gè)相交節(jié)理面的交線處, 顆粒將出現(xiàn)重合或嵌入情況 ,從而影響 節(jié)理網(wǎng)格的 劃分質(zhì)量。 因此,在節(jié) 理面填充顆 粒之前

7、,需 要對(duì)相交節(jié)理面進(jìn)行分割 也即把相交的節(jié)理面劃分為不相交的節(jié)理面 。1圖 三組節(jié)理相交Fig.1 Three groups of intersected joint faces如圖 1 所示, 節(jié)理面的形狀很簡(jiǎn)單,而且節(jié)理面的法向分別沿著三個(gè)坐標(biāo)軸,并沒有在一定范圍內(nèi)隨機(jī),即便如此,如果不劃分,在兩個(gè)面交叉的部分就會(huì)產(chǎn)生重疊的顆粒。 首先把節(jié)理面簡(jiǎn)化為空間平面上的多邊形。 在求兩個(gè)節(jié)理面是否相交時(shí), 首先求兩個(gè)空間平面的交線 ,然后再檢 驗(yàn)這一相交 直線是否和 兩個(gè)多邊形 都相交,如 果都相交, 則說明兩個(gè)節(jié)理面相交 ,同時(shí)可以 求得交點(diǎn)。 昀后利用兩 個(gè)交點(diǎn)和第 一個(gè)節(jié)理面 的頂點(diǎn)信息

8、就可以構(gòu)成兩個(gè)新的節(jié)理面,這兩個(gè)新節(jié)理面互不相交也和原來的第二個(gè)節(jié)理面不相交。 重復(fù)這一檢測(cè)過程, 直到所有初始的節(jié)理面和新產(chǎn)生的節(jié)理面都被檢測(cè)了為止。 這樣分割之后就沒有相交的節(jié)理面了。 2.2 初始顆粒生成 在用顆粒離散元進(jìn)行迭代計(jì)算之前, 首先要考慮的就是在節(jié)理面內(nèi)產(chǎn)生易于計(jì)算收斂的顆粒系統(tǒng)。這 是由于當(dāng)顆 粒數(shù)量過多 或者初始能 量很高時(shí), 迭代時(shí)間會(huì) 非常長(zhǎng),在 計(jì)算復(fù)雜問題時(shí)可能失去實(shí)用價(jià)值。 2.2.1 顆粒半 徑選取 基本原則是讓一個(gè)節(jié)理面內(nèi)填充的顆粒數(shù)量在 100-1000 左右 。 首先定義節(jié)理面的大小 D ,對(duì)于圓形節(jié)理面來,D 為圓的直徑, 對(duì)于多邊形節(jié)理面來說,D 為

9、昀長(zhǎng)對(duì)角線的長(zhǎng)度。 然后 定- 341 -顆粒材料計(jì)算力學(xué)研究進(jìn)展 義 D 為所 有節(jié)理面 D 中的昀大值。那么顆粒的半徑 r 滿足 grain1rD 1 grain 60實(shí)際算例中,r 往往取 和公式計(jì)算值相近的一個(gè)值, 比如 r 0.1273 時(shí)就取 r 0.1 。grain grain grain所有的顆粒都是同樣的半徑。 2.2.2 是顆粒 數(shù)量選取。 具體到每一個(gè)節(jié)理面應(yīng)該填充的顆粒數(shù)量, 需要首先計(jì)算節(jié)理面的面積 S,周長(zhǎng) l 和邊total界上的顆粒數(shù)量 n 。 edgeltotaln ? 2 edge2rgrain考慮到邊界上顆粒的面積只有一半是在節(jié)理面內(nèi),每個(gè)節(jié)理面填充的顆粒

10、數(shù)量 n 滿足: S 1nn 3 edge2r 2grain式中的中括號(hào)表示取大于其中值的昀小整數(shù)。 2.2.3 隨機(jī)方 法 對(duì)于節(jié)理面邊界,n 個(gè) 顆粒均勻的分布在邊界上, 對(duì)于節(jié)理面內(nèi)部, 隨機(jī)顆粒的中心坐edge標(biāo), 隨機(jī)數(shù)的平均值為節(jié)理面的中心, 隨機(jī)范圍為 D , 相當(dāng)于隨 機(jī)顆粒在一個(gè)立方體內(nèi)均勻分布。為了保證隨機(jī)產(chǎn)生的顆粒都在節(jié)理面內(nèi),產(chǎn)生的中心坐標(biāo)要投影到節(jié)理面所在的平面上,然后驗(yàn)證中心坐標(biāo)是否在節(jié)理面內(nèi),不滿足條件的顆粒將重新隨機(jī)產(chǎn)生。 同時(shí), 還需要檢驗(yàn)新產(chǎn)生的顆粒與節(jié)理面內(nèi)現(xiàn)存顆粒的距離, 這是由于如果隨機(jī)產(chǎn)生的顆粒距離很近 ,整個(gè)顆粒 系統(tǒng)開始運(yùn) 動(dòng)時(shí)就會(huì)有 很大的初始

11、 能量,這不 利于計(jì)算收 斂,也不利于節(jié)理面保持應(yīng)有的形狀。 定義一個(gè)隨機(jī)產(chǎn)生的顆粒之間的昀短距離 r , 當(dāng)新產(chǎn) 生的顆粒與min節(jié)理面內(nèi)現(xiàn)有的顆粒距離小于 r 時(shí),重新隨機(jī)產(chǎn)生這一顆粒。在實(shí)際程序中,r 是與 rmin min grain成正比的, 比例系數(shù)增 大時(shí),產(chǎn)生 顆粒的時(shí)間 增長(zhǎng)直至永 遠(yuǎn)不能完成 ,所以可以 根據(jù)經(jīng)驗(yàn)取值,本文中的程序均取 rr ?1.3 4 min grain2.3 顆粒運(yùn)動(dòng)計(jì)算 2.3.1計(jì)算流 程在迭代過程中,每個(gè)迭代步需要計(jì)算的內(nèi)容為: 在計(jì)算顆粒之間的接觸力時(shí),需要用鄰近空間的方法更新接觸列表。 2.3.2 接觸力 的計(jì)算 在計(jì)算過 程 中,接觸 的

12、發(fā)現(xiàn),顆 粒 的接觸力 的 計(jì)算,阻 尼 的引入和 顯 示迭代都 和 傳統(tǒng)離散元沒有區(qū)別。 顆粒的平動(dòng)方程為 - 342 -段文杰等: 基于平衡顆粒法生成地質(zhì)體復(fù)雜節(jié)理面網(wǎng)格的方法研究ki?dVimF+F 5ic,ijdijdtj1式中,m 與 V 分別為顆 粒 i 的質(zhì)量和 速度。t 為時(shí) 間,F 與 F 分別為顆 粒 i 與 j 的接 觸力 與i i c,ij d,ij粘性接觸阻尼力,k 為所 有與顆粒 i 接 觸的顆?？倲?shù)。 由于顆粒只是用于生成網(wǎng)格結(jié)點(diǎn), 所以計(jì)算中不計(jì)重力,不計(jì)算轉(zhuǎn)動(dòng),于是也就不需要計(jì)算切向的摩擦力和力矩。 法向接觸力的計(jì)算采用 Cundall 模型 ,有:FF-k

13、?n 6 ci, j cnij n n?Fc-Vnn 7dn,ij n ijk n c 式中, 為 法向彈簧剛度, 為單位 法向向量。 為阻尼系數(shù) 。n n R RR-R8 ni j j ii j V R為顆粒 與 接觸時(shí)的侵入深度, 為兩個(gè)處于接觸顆粒接觸點(diǎn)的相對(duì)速度, 為顆粒ij ini 的坐標(biāo)。 采用的接觸發(fā)現(xiàn)算法為 Static Cell 法 , 也就是靜態(tài)盒子法。 將顆粒體系所占據(jù)的空間劃分mmm?成規(guī)則的網(wǎng)格, 將三維空間劃分為 個(gè)立方體單元, 粒子就當(dāng)前的位置被分配到某個(gè)單元,只有在同一個(gè)單元或直接相鄰單元內(nèi)的粒子間才可能發(fā)生相互作用。 2.3.3 控制 方程和邊界條件為了保證顆

14、粒在運(yùn)動(dòng)之后依然能代表節(jié)理面,需要引入控制方程和邊界條件。 1 控制方程 節(jié)理面中心對(duì)屬于節(jié)理面的粒子施加引力,方向指向節(jié)理面中心,大小為 F 。 g?Fk? R R 9 gi gravity pc其中 R 為顆 粒中心坐標(biāo),k 為引力系數(shù)。這樣當(dāng)顆粒接近中心時(shí),引力就會(huì)減小。 pc gravity為了確定一個(gè)合適的 k , 保證顆粒不會(huì)因?yàn)橐^大而擠在一起, 產(chǎn)生很大的變形, 需要gravity用下面的公式確定 k 。 gravity Fk D 10 edge gravity kDgravity ?/ Fk11 edge nknkD?gravity 0.01 12 22rkrgrain n

15、 grain2kr rn grain graink 0.01 0.02k13 gravity nDD - 343 -顆粒材料計(jì)算力學(xué)研究進(jìn)展 其中 F 是 節(jié)理面邊緣處顆粒所受的引力,? 是引力引起的顆粒形變, 是相對(duì)形變,edge2rgrain用節(jié)理面邊緣處顆粒相對(duì)形變?yōu)?1% 為條件可以解出 k 。 gravity2 邊界條件 對(duì)于顆粒 的 運(yùn)動(dòng)有兩 個(gè) 限制條件 , 一個(gè)是運(yùn) 動(dòng) 范圍限制 , 一個(gè)是運(yùn) 動(dòng) 一致性限 制 。運(yùn)動(dòng)范圍限制是指顆粒在運(yùn)動(dòng)時(shí), 在初始節(jié)理面法向上, 偏離的距離不能超過 2r , 在初始節(jié)理grain10r面沿面方向上,偏離的距離不能超過 ,其表達(dá)式為 grai

16、n- RRn ?2r 14 origin originFace grain?- RRR -R n n 10r15 origin origin originFace originFace grainR R n 式中, 為顆?，F(xiàn)在的位置, 為顆粒的初始位置, 為初始節(jié) 理面法向。origin originFacearcsin1/15 4這一邊界條件保證了經(jīng)過迭代,節(jié)理面的法向變化不會(huì)超過 ,約等于 度,1/6D 節(jié)理面的中心位置不會(huì)變化超過 。運(yùn)動(dòng)一致 性 限制是指 屬 于一個(gè)節(jié) 理 面的顆粒 在 每一個(gè)時(shí) 間 步結(jié)束后 仍 然需要在 同 一 個(gè) 平面內(nèi),所以需要在一步迭代結(jié)束后調(diào)整顆粒的位置。R

17、RR?R ?n n 16 average average averageR n其中 和 表 示一個(gè)節(jié)理面的顆粒在運(yùn)動(dòng)后的平均位置和任選兩個(gè)顆粒和平均位置average average組成的平面 的平均法向 。經(jīng)過這一 調(diào)整之后, 一個(gè)節(jié)理面 的顆粒都被 調(diào)整到它們 的位置決定的一個(gè)平均平面上,維持著在同一平面的狀態(tài)。2.3.4 收斂判 據(jù) 由于顆粒只是用于生成網(wǎng)格節(jié)點(diǎn), 所以不需要像普通離散元一樣計(jì)算到能量和不平衡力收斂。只需要判斷顆粒之間的昀大嵌入深度 0.2 是否滿 足條件即可。在程序中,也就是把nFcn,ij 0.2作為收斂判據(jù)。 實(shí)際程序中, 采用一次運(yùn)行 10 步, 然后監(jiān)視 F 的

18、cn, ijkrn grain方法來實(shí)現(xiàn)。 3 算例 3.1 雙節(jié)理 兩個(gè)互 相垂 直而且 相交 的矩形 節(jié)理 面,第 一個(gè) 節(jié)理面 法向 為 x 方向, 四個(gè)頂 點(diǎn)坐 標(biāo) 為0,3,30,-3,30,-3,-30,3,-3 ;第二 個(gè)節(jié)理面 法向?yàn)?z 方向,四 個(gè)頂點(diǎn)坐 標(biāo)為2,0,0-2,0,0-2,4,02,4,0 AFM。兩個(gè)節(jié)理面都沒有完全被對(duì)方分割,如果按照 法生成網(wǎng)格,在交線端點(diǎn)處的網(wǎng)格質(zhì)量不能保證。2 3 3采用平衡顆粒法生成的網(wǎng)格節(jié)點(diǎn)如圖 、 所示 。如圖 所 示,通過該方法產(chǎn)生的顆粒排2布較為均勻,顆粒間的空隙不足以再容納一個(gè)顆粒。如圖 所示,三個(gè)顏色表示原來的兩個(gè)3 1

19、 2 3 節(jié)理面被分成了 個(gè)不 同的節(jié)理面, 粉色和綠色面 號(hào)面 、 號(hào)面 之間被紫色面 號(hào)面 擠出- 344 -段文杰等: 基于平衡顆粒法生成地質(zhì)體復(fù)雜節(jié)理面網(wǎng)格的方法研究 了一道縫, 這就保證了 粉色面和紫 色面相交的 位置的網(wǎng)格 不會(huì)是畸形 網(wǎng)格。每個(gè) 節(jié)理面中的顆粒數(shù)分別為藍(lán)色面 551 ,粉色面 652 ,青色面 600 。3 號(hào)面 1 號(hào)面 2 號(hào)面圖 2 兩個(gè)節(jié)理面處理之后的整體視圖 Fig.2 A holistic view of the two joint faces after the procedure 3 號(hào)面 1 號(hào)面2 號(hào)面 3圖 三個(gè)面顆粒分布情況的獨(dú)立特寫圖Fi

20、g.3 Separate view of the three joint faces 3.2 多節(jié)理 本算例模型區(qū)域范圍為 10m×10m×10m ,原點(diǎn)為0,0,0 。共有 3 組 128 個(gè) 節(jié)理面,每個(gè)節(jié)理面均呈圓 盤狀。輸入 程序的原始 數(shù)據(jù)包括: 裂隙編號(hào), 中心坐標(biāo), 半徑,傾向 和傾角等,且同一組節(jié)理面的傾向和傾角將在一定范圍內(nèi)隨機(jī)。 如圖 4 所示 , 給出了模型的三維結(jié)構(gòu)圖, 每組節(jié)理系統(tǒng)由不同的顏色表示。 由于程序通過空間多邊形來描述節(jié)理面,故本算例中采用正六邊形來表示圓盤狀節(jié)理。圖 4 未處理的節(jié)理面的三維結(jié)構(gòu)圖 Fig.4 Structure of

21、the multiple joint faces before the procedure通過節(jié)理分割操作及邊界截?cái)嗖僮骱?節(jié)理模型中將不可避免地出現(xiàn)若干極細(xì)小節(jié)理面,考慮到這些 節(jié)理面對(duì)地 質(zhì)體性質(zhì)的 影響與其他 節(jié)理面相比 可以忽略, 故在處理過 程中刪去了劃分后邊的長(zhǎng)度小于 2 倍顆粒直徑的面。 采用平衡 顆 粒法進(jìn)行 節(jié) 理網(wǎng)格節(jié) 點(diǎn) 的均勻化 布 設(shè),計(jì)算 完 畢后的節(jié) 理 網(wǎng)格節(jié)點(diǎn) 如 圖 5 所示。 如圖 5.3 所示, 在周圍多個(gè)節(jié)理面上的顆粒影響下, 節(jié)理面上的顆粒仍然保持了節(jié)理面的- 345 -顆粒材料計(jì)算力學(xué)研究進(jìn)展 形態(tài),兩個(gè)節(jié)理面之間沒有發(fā)現(xiàn)有穿透的現(xiàn)象。如圖 5.1

22、 和 5.4 所示, 能觀察到原本屬于同一組的節(jié)理面 仍然保持基 本平行,并 沒有因?yàn)轭w 粒運(yùn)動(dòng)而使 節(jié)理面的法 向發(fā)生大的 變化。如圖5.2 所示,由 于程序使用空間多邊形表述節(jié)理面,所以切割后會(huì)產(chǎn)生三角形節(jié)理面,與模型中的圓盤狀節(jié)理面差別較大。 5.1 5.2 5.3 5.4圖 5 幾個(gè)不同角度的處理后結(jié)果圖 Fig.5 Different views of the result4 結(jié)論 平衡顆粒法的核心是在節(jié)理面內(nèi)隨機(jī)生成顆粒, 進(jìn)而通過顆粒離散元的計(jì)算獲取高質(zhì)量的節(jié)理網(wǎng)格節(jié) 點(diǎn)。其中的 關(guān)鍵是在計(jì) 算中保持節(jié) 理面的位置 和法向基本 不變,也就 是牽引勢(shì)函數(shù)和位移邊界條件的引入。 雙節(jié)理和多節(jié)理的算例表明, 該方法產(chǎn)生的節(jié)理網(wǎng)格節(jié)點(diǎn)排布均勻,在節(jié)理交線處無奇異節(jié)點(diǎn),這為后續(xù)高質(zhì)量節(jié)理網(wǎng)格的劃分提供了基礎(chǔ)。 平衡顆粒法的后續(xù)工作包括: 基于平衡顆粒產(chǎn)生高質(zhì)量的節(jié)理網(wǎng)格, 計(jì)算過程中節(jié)理面中心和方向的 調(diào)整對(duì)節(jié)理 系統(tǒng)物理力 學(xué)特性的影 響,及節(jié)理 分割操作中 相交節(jié)理面 的連通性問題處理等。 參考文獻(xiàn) 1 耿克勤, 吳永 平. 拱壩 和壩 肩巖體的力學(xué) 和滲流的耦合 分析實(shí)例J. 巖石力學(xué)與工 程學(xué)報(bào), 1997, 162: 125-1312 陶連金, 蘇生瑞, 張倬元. 節(jié)理 巖體邊坡的動(dòng)力穩(wěn)定性分析J. 工程地質(zhì)學(xué)報(bào),