無(wú)網(wǎng)格伽遼金法求解固結(jié)方程的數(shù)值誤差分析

1、第 23卷 第 18期巖石力學(xué)與工程學(xué)報(bào) 23(18:311731212004年 9月 Chinese Journal of Rock Mechanics and Engineering Sept . , 20042003年 4月 30日收到初稿, 2003年 5月 30日收到修改稿。 * 中國(guó)博士后基金 (2003033168和吉林大學(xué)青年教師基金資助項(xiàng)目。無(wú)網(wǎng)格伽遼金法求解固結(jié)方程的數(shù)值誤差分析*張延軍 1,2 王恩志 1 王思敬 1(1清華大學(xué)水利水電工程系 北京 100084 (2吉林大學(xué)建設(shè)工程學(xué)院 長(zhǎng)春 130026摘要 作為一種新的計(jì)算方法, 無(wú)網(wǎng)格伽遼金法 (EFGM有著自己的

2、構(gòu)成特點(diǎn), 在求解固結(jié)方程時(shí)也會(huì)產(chǎn)生數(shù)值誤 差。 討論了 EFGM 中一些降低數(shù)值震蕩誤差的影響因素和解決方法, 首次給出了固結(jié) EFGM 離散方程的誤差分析 不等式,針對(duì)該公式提出了在 EFGM 計(jì)算中降低誤差的參數(shù)最優(yōu)選取方法。最后,通過(guò)二維條形地基理論模型的 數(shù)值試驗(yàn),驗(yàn)證了積分 Cell 精細(xì)程度對(duì) EFGM 求解固結(jié)方程的初始孔壓精度和穩(wěn)定性的影響。 關(guān)鍵詞 土力學(xué),無(wú)網(wǎng)格伽遼金法,固結(jié)方程,數(shù)值模擬,誤差分析分類號(hào) TU 43 文獻(xiàn)標(biāo)識(shí)碼 A 文章編號(hào) 1000-6915(200418-3117-05NUMERICAL ERROR ANALYSIS FOR CONSOLIDATIO

3、NEQUATION BY ELEMENT-FREE GALERKIN METHODZhang Yanjun1,2, Wang Enzhi1, Wang Sijing1(1Department of Hydraulic and Hydropower Engineering, Tsinghua University, Beijing 100084 China (2College of Environment and Construction Engineering, Jilin University, Changchun 130026 .ChinaAbstract The new numerica

4、l method based on element-free Galerkin method (EFGM and finite element method (FEM is a promising method to solve the consolidation problem by using element background mesh and shape function from the moving least square approximation. And it may also produce some numerical errors in solving the co

5、nsolidation equation. In this paper, the FEM-EFGM coupling method is developed and numerically implemented. An error inequality formula for EFGM is presented for consolidation problem based on the Terzaghis theory. Then, some influence factors are discussed to reduce the oscillatory errors in EFGM c

6、alculation , such as the time factor, the influence of domain and the integral cell structure. In the end, through numerical experiments of the model of two-dimensional stripe foundation, the effect of the integral refined degree of cell structure is validated on both accuracy and stability of the i

7、nitial pore-pressure in solving the consolidation equation with EFGM. The work of the paper will help to enhance the possibility of the application of EFGM to geotechnical engineering and also provide a new numerical analysis tool for solving the solid-fluid coupling problem.Key words soil mechanics

8、, element-free Galerkin method, consolidation equation, numerical simulation, error analysis1 前 言Terzaghi 固結(jié)方程和 Biot 固結(jié)方程 1是土力學(xué)的標(biāo)志方程,各種數(shù)值模擬技術(shù)如有限差分法、有限元法、邊界元法等在不同時(shí)期都有應(yīng)用,并取得 了大量的成果。在數(shù)值模擬中,最常用的方法是在 時(shí)間域上采用差分格式,在空間上采用網(wǎng)格單元 3118 巖石力學(xué)與工程學(xué)報(bào) 2004年離散。這種模擬思想在過(guò)去 20多年取得了很大成 功 2。在研究和應(yīng)用中發(fā)現(xiàn),由于該算法的特點(diǎn), 時(shí)間步長(zhǎng)的降低雖然會(huì)提高固結(jié)運(yùn)

9、算精度,但時(shí)間 步長(zhǎng)有一個(gè)下限,當(dāng)?shù)陀谠撓孪迺r(shí),會(huì)出現(xiàn)數(shù)值振 蕩現(xiàn)象。 文 3在數(shù)值求解 Biot 方程時(shí)首先探討了該 問(wèn)題,文 4系統(tǒng)研究并推導(dǎo)了時(shí)間步長(zhǎng)的下限公 式, 文 5也探討了該問(wèn)題, 文 6在用 DEM-FEM 耦 合分析多孔介質(zhì)的變形時(shí), 也提到了固結(jié)精度問(wèn)題。 常用的數(shù)值解法是有限元 (FEM和邊界元法 等,它們是解初值或邊值問(wèn)題的有力工具,都將整 個(gè)研究域離散成多個(gè)網(wǎng)格單元,設(shè)在小的域中簡(jiǎn)化 未知量,一般假設(shè)單元間是連續(xù)的。為了提高近似 函數(shù)的精度, FEM 多采用高次插值或細(xì)分網(wǎng)格,但 增加了計(jì)算成本。在固結(jié)變形分析領(lǐng)域,對(duì)于大變 形、裂隙分析和局部應(yīng)變、多相介質(zhì)耦合以及求

10、解 域高梯度等需重新劃分網(wǎng)格的問(wèn)題,其前、后處理 工作量大,并造成誤差,盡管已有了網(wǎng)格自動(dòng)生成 器,但還是計(jì)算昂貴,精度降低。為了避免上述問(wèn)題,人們發(fā)展了無(wú)網(wǎng)格技術(shù)。 從文 7提出無(wú)網(wǎng)格 Galerkin 方法 (EFGM開(kāi)始,幾 年間該法迅速在計(jì)算力學(xué)的各個(gè)領(lǐng)域全面推廣,解 決了彈性介質(zhì)中裂隙的擴(kuò)展 8和 Timoshenko 梁的剪 切閉鎖 9等傳統(tǒng) FEM 無(wú)法求解的問(wèn)題。文 10用該 法分析了三維彈性和彈塑性問(wèn)題,文 11, 12在 20012002年用點(diǎn)積分法等分析了固結(jié)大變形問(wèn) 題。 從 2001年起國(guó)內(nèi)學(xué)者在文 1316中對(duì) EFGM 的計(jì)算參數(shù)選擇及其在固結(jié)上的應(yīng)用作了深入探

11、討。無(wú)網(wǎng)格伽遼金法在固結(jié)分析時(shí),也會(huì)遇到時(shí) 間步長(zhǎng)有一個(gè)下限的問(wèn)題。由于 EFGM 的局部化 特點(diǎn),它的時(shí)間步長(zhǎng)下限公式比傳統(tǒng)有限元公式更 為復(fù)雜。為了發(fā)展 EFGM 在固結(jié)上的應(yīng)用,必須探 討它在與時(shí)間差分聯(lián)合使用時(shí)產(chǎn)生誤差的原因和消 除辦法。 本文首次給出了固結(jié) EFGM 離散方程的誤 差分析不等式和影響因數(shù),并通過(guò)理論模型進(jìn)行了 驗(yàn)證分析。2 EFGM在 Terzaghi 固結(jié)方程上的 應(yīng)用2.1Terzaghi 固結(jié)方程的 EFGM 離散形式對(duì)于一維飽和均質(zhì)彈性多孔介質(zhì)的 Terzaghi 固 結(jié)方程,文 5已經(jīng)推導(dǎo)出了在時(shí)間間隔 t 的一維 有限元控制方程, 本文按相同的原則推出 E

12、FGM 的 公式為22ttpzpC v=(1式中:v C 為固結(jié)系數(shù), 與兩相的滲透性和壓縮性有 關(guān); p 為計(jì)算孔隙水壓力; 為外力載荷。用標(biāo)準(zhǔn) 法處理時(shí)間間隔 t , 式 (1在 Galerkin 的加權(quán)殘差處理后得+h hvzzqzptCzqp22dd(dddd +=hhh zqzqpq 222d(d(2式中:p , q H 1, H 0; H 1, H 0分別為一維和 零維 Sobelev 空間; 為 Lagrange 乘子; q為試探 函數(shù); p 代表時(shí)間間隔 t 中的初值, EFGM 給出 近似域; 2h 為研究域的長(zhǎng)度, EFGM 空間離散化后, 用 n +1等距節(jié)點(diǎn), 并采用線

13、性基 1, z 定義形函數(shù)。 每個(gè)節(jié)點(diǎn)的影響域至少 2倍于兩節(jié)點(diǎn)距離 (h , 本 文取 hd =5. 3m。設(shè) N i 為 EFGM 形函數(shù),在 Terzaghi 問(wèn)題的標(biāo) 準(zhǔn)邊界條件即單面排水、 雙面排水情況下, Lagrange 乘子 自然滿足,可以不考慮。因而,基本邊界條 件能精確強(qiáng)加為=niniiniiipphNppN2(0(3式 (1可寫(xiě)成 EFGM 剛度矩陣形式:=+=+=ddd(''n1J11zNaNNHNNMapMptHCMjjjiijjiijjnjnjjijjijvij(4為了求解積分,可選擇一種積分網(wǎng)格 (cell結(jié) 構(gòu) , 在網(wǎng)格上采用高斯積分。 求解式

14、 (4可得到任一 時(shí)刻各點(diǎn)的孔壓值。根據(jù)上述原理,可以在微機(jī)上 進(jìn)行分析。2.2數(shù)值近似和穩(wěn)定性研究當(dāng)用 FEM 計(jì)算時(shí),精度要求每一個(gè)節(jié)點(diǎn)上的 p 不能超過(guò) , i 代表了第 i 節(jié)點(diǎn)的孔壓和與之第 23卷 第 18期 張延軍等 . 無(wú)網(wǎng)格伽遼金法求解固結(jié)方程的數(shù)值誤差分析 3119·相對(duì)的外載的相差系數(shù)。因此 p 可用下式表達(dá):1( i i p =; i 0 (=i 1, 2, 1n (5 由式 (4和 (5可得+11V 1111V 111V 11V 111V 1111V 11 n n n n n n n n k n k k k n n tH C M tH C MtH C M t

15、H C M tH C M tH C M , , , , , , , , , , , L M M L M M L +=n n n n k k n k tH C M tH C M tH C M , , , , , 1V 10V 001V 0111 MM M M (6 將式 (6兩邊都乘以 T (的轉(zhuǎn)置 ,得出標(biāo)量方程。 由于左邊矩陣是正定的,式 (6右邊也總是正定的。 另外,由于對(duì)稱,有 11=n , 因此,可得到如下EFGM 誤差分析不等式:+ ( (02V 02201V 011, , , , tH C M tH C M (0V 0, , k k k tH C M +L 0 (7t 0V 02V

16、 201V 10022011 , , , , , , k k k k H C H C H C M M M +L (8考慮到 EFGM 的形函數(shù)及其導(dǎo)數(shù)的特點(diǎn), 上式 中的質(zhì)量矩陣有以下性質(zhì):01, M >02, M >>0, k M (9而滲透矩陣 H 則正負(fù)不定,無(wú)法得出準(zhǔn)確的規(guī)律。 由以上推導(dǎo)可知,決定 t 穩(wěn)定的因素不同于 FEM 的 1/6網(wǎng)格,是一個(gè)多因素的組合,主要由邊界點(diǎn) 的質(zhì)量矩陣和滲流矩陣決定。而邊界點(diǎn)的 M 和 H 又由 EFGM 的影響域、權(quán)函數(shù)、基向量、節(jié)點(diǎn)的積 分構(gòu)成。 FEM 中的 t 的穩(wěn)定計(jì)算公式為t V2(61C h (10在 EFGM 中,

17、由于形函數(shù)的計(jì)算中涉及到 A -1,目 前尚無(wú)法得到顯函數(shù)表達(dá)式,只有當(dāng)?shù)染嚯x 2, 3節(jié)點(diǎn)分布時(shí), 節(jié)點(diǎn)影響域才等同 FEM 的要求。 但若 大于 3節(jié)點(diǎn),則無(wú)法得出解析表達(dá)式。一般說(shuō)來(lái), 影響域中的節(jié)點(diǎn)愈少,則計(jì)算精度愈低?？梢?jiàn), 由于 EFGM 不滿足邊界的 Kronecker 性質(zhì),即Ij j I x N ( (11t 的穩(wěn)定因素也比 FEM 復(fù)雜, t 的取值可能 大于 1/6×V C h / (2。 EFGM 的固結(jié)方程離散形式 的精度和穩(wěn)定性的主要影響因素歸納起來(lái)有:(1t 大小; (2 權(quán)函數(shù)類型; (3 影響域大小; (4 基 向量次數(shù)和節(jié)點(diǎn)積分結(jié)構(gòu)。3 瞬時(shí)加載固

18、結(jié)的 EFGM 數(shù)值誤差分析在 EFGM 中常采用負(fù)指數(shù)型函數(shù)或樣條函數(shù)作 為權(quán)函數(shù)?？紤]到計(jì)算花費(fèi),筆者通過(guò)研究認(rèn)為, 負(fù)指數(shù)型權(quán)函數(shù)和二次基向量在固結(jié)問(wèn)題上精度較 高 (詳見(jiàn)文 12和另文 , 因此, 在不考慮權(quán)函數(shù)和基 向量影響的前提下,討論 t 大小、影響域大小、節(jié) 點(diǎn)積分精度對(duì)數(shù)值精度和穩(wěn)定性的影響,可省時(shí)又 經(jīng)濟(jì)。3.1 時(shí)間因子的穩(wěn)定性分析在 Terzaghi 固結(jié)方程 EFGM 數(shù)值模擬的誤差和 穩(wěn)定性分析中,理論上除了滲透邊界點(diǎn)外的所有節(jié) 點(diǎn)的初始孔壓應(yīng)都等于瞬間外荷載。類似于 FEM 穩(wěn)定公式 (10, t 不能過(guò)小,否則引起矩陣病態(tài); 但過(guò)大將引起波動(dòng),文 5給出式 (1

19、0作參考。筆者也設(shè)計(jì)了一個(gè) EFGM 影響因數(shù) 2maxV /d t C , 但 EFGM 求解的離散形式比 FEM 復(fù)雜, 涉及的參數(shù)較多。 本 文中 max d 代表節(jié)點(diǎn)計(jì)算域的最大值,在研究時(shí)間因 子時(shí), 將其固定為 =max d 3.5h , h 為節(jié)點(diǎn)之間距離的 平均值,而差分系數(shù) =1。圖 1展示了不同時(shí)間因 子的計(jì)算誤差百分比。在 FEM 中,文 5利用時(shí)間 因子 6(2V /h t C 作為收斂的準(zhǔn)則;但在 EFGM 中,需要以 10(2maxV /d t C 作為準(zhǔn)則, t 過(guò)大、過(guò)小 (圖 1(a, c 都會(huì)引起 EFGM 計(jì)算數(shù)值較大振蕩現(xiàn)象。 3.2 影響域的誤差精度分

20、析上述的時(shí)間因子準(zhǔn)則是在影響域大小 d 為常數(shù) 得出的,當(dāng)影響域改變時(shí),計(jì)算結(jié)果也有變化。圖 2給出了計(jì)算條件為 =t 0.05 d,計(jì)算域劃分為規(guī)則 等間距 21個(gè)節(jié)點(diǎn),影響域參數(shù) max d 分別為 1.1, 2.1, 4.1 h的計(jì)算結(jié)果。圖 2中 p /0p 表示孔壓比, 0p 為 初始孔壓, H 代表計(jì)算模型厚度,本算例中取 =H2 m。 從圖中可看出在靠近邊界點(diǎn)處, 影響域愈小計(jì) 算誤差也愈小,隨著遠(yuǎn)離邊界,影響域的影響逐步 變小乃至消失。這種現(xiàn)象產(chǎn)生的原因是當(dāng)影響域小 時(shí), 其計(jì)算接近 FEM , 而 FEM 滿足解答 Kronecker式 (11,所以邊界處誤差變小。 3120

21、 巖石力學(xué)與工程學(xué)報(bào) 2004年圖 1 不同時(shí)間因子的 EFGM 計(jì)算誤差圖 Fig.1 The calculation errors of EFGM with different timefactor圖 2 不同影響域的孔壓消散曲線Fig.2 Pore-pressure isochrones with different influencedomain3.3 節(jié)點(diǎn)的積分結(jié)構(gòu)的影響在計(jì)算精度和穩(wěn)定性的影響因素中,影響最大 的是節(jié)點(diǎn)的積分結(jié)構(gòu)的精細(xì)度,這從公式 (4可清晰 看出。如果 EFGM 中的 cell 積分結(jié)構(gòu)足夠精細(xì),計(jì) 算誤差可明顯降低。下面的一系列數(shù)值試驗(yàn)可明顯 證實(shí)這一觀點(diǎn)。表

22、1給出了 3種 cell 積分方案, 3種情況的節(jié) 點(diǎn)誤差百分比由圖 3給出。最大誤差出現(xiàn)在圖 3(a的積分結(jié)構(gòu)。隨著結(jié)構(gòu)的細(xì)化,誤差值有了很大的 降低,可見(jiàn)積分結(jié)構(gòu)的選擇對(duì)計(jì)算精度和穩(wěn)定性有 巨大的影響。圖 3的計(jì)算時(shí)間間隔為 0.01 d。表 1 積分方案表 Table 1 Integral scheme方案節(jié)點(diǎn)數(shù)目 cell 積分結(jié)構(gòu)(a 21 20(b 21 40 (c 21 100圖 3 不同方案 cell 積分結(jié)構(gòu)的孔壓計(jì)算誤差曲線 Fig.3 Calculation error of pore-pressure in differentintegration cells4 二維問(wèn)

23、題的應(yīng)用研究基于以上理論,對(duì)一個(gè)標(biāo)準(zhǔn)的條形基礎(chǔ)的二維固結(jié)變形問(wèn)題進(jìn)行模擬驗(yàn)證。有關(guān)二維固結(jié)的EFGM 控制方程的推導(dǎo)和邊界條件請(qǐng)參見(jiàn)文 14。 工程概況為均質(zhì)土層,厚度取 12 m,土的天然容重=17 kN/m3,彈性模量 =E 4 MPa,泊松比 =0.3, 受寬度 6 m的條形基礎(chǔ)的均布荷載 =q 20 kPa。 滲透系數(shù) =k 2.5×10-6cm/s。 由于對(duì)稱, 計(jì)算僅取一半進(jìn)行。地基計(jì)算寬度取基礎(chǔ)一半寬 (3 m的 5倍, 即 15 m。地基左右兩邊在水平方向固定,底面在水 平和豎直兩個(gè)方向均固定;基礎(chǔ)左右兩邊在水平方 向固定。地基在左邊和底面以及基礎(chǔ)部分不排水, 其他邊

24、界可以自由排水。為了進(jìn)一步說(shuō)明 EFGM 的計(jì)算誤差, 在確定時(shí) 間步長(zhǎng)的基礎(chǔ)上,改變 EFGM 的 cell 積分?jǐn)?shù)量, 分為 60個(gè) cell 結(jié)構(gòu)和 30個(gè) cell 結(jié)構(gòu)這兩種形式。 計(jì)算出的基礎(chǔ)中心和基礎(chǔ)邊緣的 0.01 d的孔壓隨深 度變化曲線見(jiàn)圖 4。這里需要說(shuō)明的是計(jì)算中為了 發(fā)揮 EFGM 計(jì)算優(yōu)勢(shì), 這兩種計(jì)算情況都排列了隨 機(jī)孔壓節(jié)點(diǎn) 77個(gè),改變的僅僅是積分網(wǎng)格的數(shù)量。從圖 4中可清晰看出,隨著積分 cell 的增加, 孔壓值波動(dòng)減小,在靠近邊界處的孔壓梯度大的地 方突變值明顯變平滑。 這證明了本文節(jié) 3.3的結(jié)論。 圖中基礎(chǔ)邊緣的孔壓波動(dòng)大于基礎(chǔ)中心的原因是:在計(jì)算中

25、設(shè)計(jì)了基礎(chǔ)底面不透水,而其他地面是自由透水邊界,基礎(chǔ)邊緣是邊界條件突變處,所以在 靠近地表處的 -2 m左右產(chǎn)生較大的誤差。由此可見(jiàn), 在 EFGM 求解固結(jié)方程的數(shù)值分析0.0 -0.4 +0.4 / % 0.0 -0.4 +0.4 0.0-0.4+0.4 / % / %(a 10C V t /d max= 0.26 2(b 10C V t /d max = 1 2 (c 10C V t /d max = 2.620 0.2 0.41.01.23/41/21/4p /p 0z /Hd max = 1.1 hd max = 2.1 h d max = 4.1 h0.0-0.4+0.4 / %0

26、.0 -0.4+0.4 0.0-0.4+0.4 / % / %(a 20個(gè) cell 積分結(jié)構(gòu) (b 40個(gè) cell 積分結(jié)構(gòu)(c 100個(gè) cell 積分結(jié)構(gòu)第 23卷 第 18期 張延軍等 . 無(wú)網(wǎng)格伽遼金法求解固結(jié)方程的數(shù)值誤差分析 3121·中, cell 積分結(jié)構(gòu)的精細(xì)程度對(duì)計(jì)算的誤差和穩(wěn)定 性有巨大影響。在有條件的情況下,可通過(guò)增加積 分結(jié)構(gòu)數(shù)量,達(dá)到提高初始孔壓計(jì)算精度的效果。 圖 4 條形基礎(chǔ)不同積分 cell 結(jié)構(gòu)的孔壓計(jì)算曲線 Fig.4 Calculation curves of pore-pressure with differentintegration

27、cells in stripe foundation5 結(jié) 論本 文 介 紹 了 EFGM 的 發(fā) 展 概 況 , 推 導(dǎo) 了Terzaghi 固結(jié)方程 EFGM 的離散形式,并對(duì) EFGM 數(shù)值計(jì)算固結(jié)中初始孔壓的誤差和穩(wěn)定性進(jìn)行了研 究,得出如下結(jié)論:(1 由于 EFGM 的理論特點(diǎn),其數(shù)值近似和穩(wěn) 定性較傳統(tǒng) FEM 復(fù)雜, 本文首次給出了固結(jié) EFGM 離散方程的誤差分析不等式。(2 建立了類似于 FEM 的 EFGM 時(shí)間穩(wěn)定因子 的數(shù)學(xué)公式,并給出了全新的公式指標(biāo)和系數(shù)。(3 探討了 EFGM 中影響域大小對(duì) EFGM 計(jì)算 精度的影響規(guī)律。(4 一、二維數(shù)值試驗(yàn)表明,在 EFGM

28、 中積分 cell 結(jié)構(gòu)的精細(xì)程度對(duì) EFGM 計(jì)算初始孔壓有巨大 影響。(5 EFGM 作為新的數(shù)值方法,其獨(dú)特的計(jì)算 優(yōu)勢(shì)已經(jīng)在計(jì)算力學(xué)的各個(gè)領(lǐng)域得到證實(shí),深入研 究其計(jì)算精度和穩(wěn)定性,將為該方法的推廣和應(yīng)用 提供理論基礎(chǔ)。參 考 文 獻(xiàn)1Boit M A. General theory of three-dimensional consolidationJ. J. App. Physics, 1941, 12(2:155164 2Desai C S, Christan J T. Numerical Methods in Geotechnical Engineering M.New Yor

29、k:McGraw-Hill , 1977 3Booker J R, Small J C.An investigation of the stability of numerical solution of Biots equations of consolidationJ. Int. J. Solid Struct., 1975, 11(6:9079174 Vermeer P , Verruijt A. An accuracy condition for consolidation byfinite elementJ. Int. J. Num. Analy. Mech. in Geomech.

30、, 1981, 5(1:114 5朱百里,沈珠江 . 計(jì)算土力學(xué) M. 上海:上?？萍汲霭嫔? 1990, 1031056 Hormoz Modaressi , Philippe Aubert. A diffuse element-finite elementtechnique for transient coupled analysisJ. Int. J. for Num. Methods in Engng. , 1996, 39(5:3 8093 8387 Ted Belytschko , Lu Y Y, Gu L. Element-free Galerkin methodsJ. Int.J. for Num. Methods in Engng., 1994, 37(3:229256.8 Yu Xu , Sunil Saigal. An element free Galerkin formulation for stablecrack growth in an elastic solidJ. Computer M