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

2、構(gòu)成特點, 在求解固結(jié)方程時也會產(chǎn)生數(shù)值誤 差。 討論了 EFGM 中一些降低數(shù)值震蕩誤差的影響因素和解決方法, 首次給出了固結(jié) EFGM 離散方程的誤差分析 不等式,針對該公式提出了在 EFGM 計算中降低誤差的參數(shù)最優(yōu)選取方法。最后,通過二維條形地基理論模型的 數(shù)值試驗,驗證了積分 Cell 精細程度對 EFGM 求解固結(jié)方程的初始孔壓精度和穩(wěn)定性的影響。 關(guān)鍵詞 土力學,無網(wǎng)格伽遼金法,固結(jié)方程,數(shù)值模擬,誤差分析分類號 TU 43 文獻標識碼 A 文章編號 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是土力學的標志方程,各種數(shù)值模擬技術(shù)如有限差分法、有限元法、邊界元法等在不同時期都有應用,并取得 了大量的成果。在數(shù)值模擬中,最常用的方法是在 時間域上采用差分格式,在空間上采用網(wǎng)格單元 3118 巖石力學與工程學報 2004年離散。這種模擬思想在過去 20多年取得了很大成 功 2。在研究和應用中發(fā)現(xiàn),由于該算法的特點, 時間步長的降低雖然會提高固結(jié)運

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

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

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

12、FGM 的 公式為22ttpzpC v=(1式中:v C 為固結(jié)系數(shù), 與兩相的滲透性和壓縮性有 關(guān); p 為計算孔隙水壓力; 為外力載荷。用標準 法處理時間間隔 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 代表時間間隔 t 中的初值, EFGM 給出 近似域; 2h 為研究域的長度, EFGM 空間離散化后, 用 n +1等距節(jié)點, 并采用線

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

14、 (4可得到任一 時刻各點的孔壓值。根據(jù)上述原理,可以在微機上 進行分析。2.2數(shù)值近似和穩(wěn)定性研究當用 FEM 計算時,精度要求每一個節(jié)點上的 p 不能超過 , i 代表了第 i 節(jié)點的孔壓和與之第 23卷 第 18期 張延軍等 . 無網(wǎng)格伽遼金法求解固結(jié)方程的數(shù)值誤差分析 3119·相對的外載的相差系數(shù)。因此 p 可用下式表達: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)置 ,得出標量方程。 由于左邊矩陣是正定的,式 (6右邊也總是正定的。 另外,由于對稱,有 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ù)及其導數(shù)的特點, 上式 中的質(zhì)量矩陣有以下性質(zhì):01, M >02, M >>0, k M (9而滲透矩陣 H 則正負不定,無法得出準確的規(guī)律。 由以上推導可知,決定 t 穩(wěn)定的因素不同于 FEM 的 1/6網(wǎng)格,是一個多因素的組合,主要由邊界點 的質(zhì)量矩陣和滲流矩陣決定。而邊界點的 M 和 H 又由 EFGM 的影響域、權(quán)函數(shù)、基向量、節(jié)點的積 分構(gòu)成。 FEM 中的 t 的穩(wěn)定計算公式為t V2(61C h (10在 EFGM 中,

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

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

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

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

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

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

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

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

25、設計了基礎底面不透水,而其他地面是自由透水邊界,基礎邊緣是邊界條件突變處,所以在 靠近地表處的 -2 m左右產(chǎ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個 cell 積分結(jié)構(gòu) (b 40個 cell 積分結(jié)構(gòu)(c 100個 cell 積分結(jié)構(gòu)第 23卷 第 18期 張延軍等 . 無網(wǎng)格伽遼金法求解固結(jié)方程的數(shù)值誤差分析 3121·中, cell 積分結(jié)構(gòu)的精細程度對計算的誤差和穩(wěn)定 性有巨大影響。在有條件的情況下,可通過增加積 分結(jié)構(gòu)數(shù)量,達到提高初始孔壓計算精度的效果。 圖 4 條形基礎不同積分 cell 結(jié)構(gòu)的孔壓計算曲線 Fig.4 Calculation curves of pore-pressure with differentintegration

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

28、 中積分 cell 結(jié)構(gòu)的精細程度對 EFGM 計算初始孔壓有巨大 影響。(5 EFGM 作為新的數(shù)值方法,其獨特的計算 優(yōu)勢已經(jīng)在計算力學的各個領(lǐng)域得到證實,深入研 究其計算精度和穩(wě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朱百里,沈珠江 . 計算土力學 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