无网格伽辽金法求解固结方程的数值误差分析

1、第 23卷 第 18期岩石力学与工程学报 23(18:311731212004年 9月 Chinese Journal of Rock Mechanics and Engineering Sept . , 20042003年 4月 30日收到初稿, 2003年 5月 30日收到修改稿。 * 中国博士后基金 (2003033168和吉林大学青年教师基金资助项目。无网格伽辽金法求解固结方程的数值误差分析*张延军 1,2 王恩志 1 王思敬 1(1清华大学水利水电工程系 北京 100084 (2吉林大学建设工程学院 长春 130026摘要 作为一种新的计算方法, 无网格伽辽金法 (EFGM有着自己的

2、构成特点, 在求解固结方程时也会产生数值误 差。 讨论了 EFGM 中一些降低数值震荡误差的影响因素和解决方法, 首次给出了固结 EFGM 离散方程的误差分析 不等式,针对该公式提出了在 EFGM 计算中降低误差的参数最优选取方法。最后,通过二维条形地基理论模型的 数值试验,验证了积分 Cell 精细程度对 EFGM 求解固结方程的初始孔压精度和稳定性的影响。 关键词 土力学,无网格伽辽金法,固结方程,数值模拟,误差分析分类号 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 固结方程和 Biot 固结方程 1是土力学的标志方程,各种数值模拟技术如有限差分法、有限元法、边界元法等在不同时期都有应用,并取得 了大量的成果。在数值模拟中,最常用的方法是在 时间域上采用差分格式,在空间上采用网格单元 3118 岩石力学与工程学报 2004年离散。这种模拟思想在过去 20多年取得了很大成 功 2。在研究和应用中发现,由于该算法的特点, 时间步长的降低虽然会提高固结运

9、算精度,但时间 步长有一个下限,当低于该下限时,会出现数值振 荡现象。 文 3在数值求解 Biot 方程时首先探讨了该 问题,文 4系统研究并推导了时间步长的下限公 式, 文 5也探讨了该问题, 文 6在用 DEM-FEM 耦 合分析多孔介质的变形时, 也提到了固结精度问题。 常用的数值解法是有限元 (FEM和边界元法 等,它们是解初值或边值问题的有力工具,都将整 个研究域离散成多个网格单元,设在小的域中简化 未知量,一般假设单元间是连续的。为了提高近似 函数的精度, FEM 多采用高次插值或细分网格,但 增加了计算成本。在固结变形分析领域,对于大变 形、裂隙分析和局部应变、多相介质耦合以及求

10、解 域高梯度等需重新划分网格的问题,其前、后处理 工作量大,并造成误差,尽管已有了网格自动生成 器,但还是计算昂贵,精度降低。为了避免上述问题,人们发展了无网格技术。 从文 7提出无网格 Galerkin 方法 (EFGM开始,几 年间该法迅速在计算力学的各个领域全面推广,解 决了弹性介质中裂隙的扩展 8和 Timoshenko 梁的剪 切闭锁 9等传统 FEM 无法求解的问题。文 10用该 法分析了三维弹性和弹塑性问题,文 11, 12在 20012002年用点积分法等分析了固结大变形问 题。 从 2001年起国内学者在文 1316中对 EFGM 的计算参数选择及其在固结上的应用作了深入探

11、讨。无网格伽辽金法在固结分析时,也会遇到时 间步长有一个下限的问题。由于 EFGM 的局部化 特点,它的时间步长下限公式比传统有限元公式更 为复杂。为了发展 EFGM 在固结上的应用,必须探 讨它在与时间差分联合使用时产生误差的原因和消 除办法。 本文首次给出了固结 EFGM 离散方程的误 差分析不等式和影响因数,并通过理论模型进行了 验证分析。2 EFGM在 Terzaghi 固结方程上的 应用2.1Terzaghi 固结方程的 EFGM 离散形式对于一维饱和均质弹性多孔介质的 Terzaghi 固 结方程,文 5已经推导出了在时间间隔 t 的一维 有限元控制方程, 本文按相同的原则推出 E

12、FGM 的 公式为22ttpzpC v=(1式中:v C 为固结系数, 与两相的渗透性和压缩性有 关; p 为计算孔隙水压力; 为外力载荷。用标准 法处理时间间隔 t , 式 (1在 Galerkin 的加权残差处理后得+h hvzzqzptCzqp22dd(dddd +=hhh zqzqpq 222d(d(2式中:p , q H 1, H 0; H 1, H 0分别为一维和 零维 Sobelev 空间; 为 Lagrange 乘子; q为试探 函数; p 代表时间间隔 t 中的初值, EFGM 给出 近似域; 2h 为研究域的长度, EFGM 空间离散化后, 用 n +1等距节点, 并采用线

13、性基 1, z 定义形函数。 每个节点的影响域至少 2倍于两节点距离 (h , 本 文取 hd =5. 3m。设 N i 为 EFGM 形函数,在 Terzaghi 问题的标 准边界条件即单面排水、 双面排水情况下, Lagrange 乘子 自然满足,可以不考虑。因而,基本边界条 件能精确强加为=niniiniiipphNppN2(0(3式 (1可写成 EFGM 刚度矩阵形式:=+=+=ddd(''n1J11zNaNNHNNMapMptHCMjjjiijjiijjnjnjjijjijvij(4为了求解积分,可选择一种积分网格 (cell结 构 , 在网格上采用高斯积分。 求解式

14、 (4可得到任一 时刻各点的孔压值。根据上述原理,可以在微机上 进行分析。2.2数值近似和稳定性研究当用 FEM 计算时,精度要求每一个节点上的 p 不能超过 , i 代表了第 i 节点的孔压和与之第 23卷 第 18期 张延军等 . 无网格伽辽金法求解固结方程的数值误差分析 3119·相对的外载的相差系数。因此 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 (的转置 ,得出标量方程。 由于左边矩阵是正定的,式 (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 的形函数及其导数的特点, 上式 中的质量矩阵有以下性质:01, M >02, M >>0, k M (9而渗透矩阵 H 则正负不定,无法得出准确的规律。 由以上推导可知,决定 t 稳定的因素不同于 FEM 的 1/6网格,是一个多因素的组合,主要由边界点 的质量矩阵和渗流矩阵决定。而边界点的 M 和 H 又由 EFGM 的影响域、权函数、基向量、节点的积 分构成。 FEM 中的 t 的稳定计算公式为t V2(61C h (10在 EFGM 中,

17、由于形函数的计算中涉及到 A -1,目 前尚无法得到显函数表达式,只有当等距离 2, 3节点分布时, 节点影响域才等同 FEM 的要求。 但若 大于 3节点,则无法得出解析表达式。一般说来, 影响域中的节点愈少,则计算精度愈低。可见, 由于 EFGM 不满足边界的 Kronecker 性质,即Ij j I x N ( (11t 的稳定因素也比 FEM 复杂, t 的取值可能 大于 1/6×V C h / (2。 EFGM 的固结方程离散形式 的精度和稳定性的主要影响因素归纳起来有:(1t 大小; (2 权函数类型; (3 影响域大小; (4 基 向量次数和节点积分结构。3 瞬时加载固

18、结的 EFGM 数值误差分析在 EFGM 中常采用负指数型函数或样条函数作 为权函数。考虑到计算花费,笔者通过研究认为, 负指数型权函数和二次基向量在固结问题上精度较 高 (详见文 12和另文 , 因此, 在不考虑权函数和基 向量影响的前提下,讨论 t 大小、影响域大小、节 点积分精度对数值精度和稳定性的影响,可省时又 经济。3.1 时间因子的稳定性分析在 Terzaghi 固结方程 EFGM 数值模拟的误差和 稳定性分析中,理论上除了渗透边界点外的所有节 点的初始孔压应都等于瞬间外荷载。类似于 FEM 稳定公式 (10, t 不能过小,否则引起矩阵病态; 但过大将引起波动,文 5给出式 (1

19、0作参考。笔者也设计了一个 EFGM 影响因数 2maxV /d t C , 但 EFGM 求解的离散形式比 FEM 复杂, 涉及的参数较多。 本 文中 max d 代表节点计算域的最大值,在研究时间因 子时, 将其固定为 =max d 3.5h , h 为节点之间距离的 平均值,而差分系数 =1。图 1展示了不同时间因 子的计算误差百分比。在 FEM 中,文 5利用时间 因子 6(2V /h t C 作为收敛的准则;但在 EFGM 中,需要以 10(2maxV /d t C 作为准则, t 过大、过小 (图 1(a, c 都会引起 EFGM 计算数值较大振荡现象。 3.2 影响域的误差精度分

20、析上述的时间因子准则是在影响域大小 d 为常数 得出的,当影响域改变时,计算结果也有变化。图 2给出了计算条件为 =t 0.05 d,计算域划分为规则 等间距 21个节点,影响域参数 max d 分别为 1.1, 2.1, 4.1 h的计算结果。图 2中 p /0p 表示孔压比, 0p 为 初始孔压, H 代表计算模型厚度,本算例中取 =H2 m。 从图中可看出在靠近边界点处, 影响域愈小计 算误差也愈小,随着远离边界,影响域的影响逐步 变小乃至消失。这种现象产生的原因是当影响域小 时, 其计算接近 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 节点的积分结构的影响在计算精度和稳定性的影响因素中,影响最大 的是节点的积分结构的精细度,这从公式 (4可清晰 看出。如果 EFGM 中的 cell 积分结构足够精细,计 算误差可明显降低。下面的一系列数值试验可明显 证实这一观点。表

22、1给出了 3种 cell 积分方案, 3种情况的节 点误差百分比由图 3给出。最大误差出现在图 3(a的积分结构。随着结构的细化,误差值有了很大的 降低,可见积分结构的选择对计算精度和稳定性有 巨大的影响。图 3的计算时间间隔为 0.01 d。表 1 积分方案表 Table 1 Integral scheme方案节点数目 cell 积分结构(a 21 20(b 21 40 (c 21 100图 3 不同方案 cell 积分结构的孔压计算误差曲线 Fig.3 Calculation error of pore-pressure in differentintegration cells4 二维问

23、题的应用研究基于以上理论,对一个标准的条形基础的二维固结变形问题进行模拟验证。有关二维固结的EFGM 控制方程的推导和边界条件请参见文 14。 工程概况为均质土层,厚度取 12 m,土的天然容重=17 kN/m3,弹性模量 =E 4 MPa,泊松比 =0.3, 受宽度 6 m的条形基础的均布荷载 =q 20 kPa。 渗透系数 =k 2.5×10-6cm/s。 由于对称, 计算仅取一半进行。地基计算宽度取基础一半宽 (3 m的 5倍, 即 15 m。地基左右两边在水平方向固定,底面在水 平和竖直两个方向均固定;基础左右两边在水平方 向固定。地基在左边和底面以及基础部分不排水, 其他边

24、界可以自由排水。为了进一步说明 EFGM 的计算误差, 在确定时 间步长的基础上,改变 EFGM 的 cell 积分数量, 分为 60个 cell 结构和 30个 cell 结构这两种形式。 计算出的基础中心和基础边缘的 0.01 d的孔压随深 度变化曲线见图 4。这里需要说明的是计算中为了 发挥 EFGM 计算优势, 这两种计算情况都排列了随 机孔压节点 77个,改变的仅仅是积分网格的数量。从图 4中可清晰看出,随着积分 cell 的增加, 孔压值波动减小,在靠近边界处的孔压梯度大的地 方突变值明显变平滑。 这证明了本文节 3.3的结论。 图中基础边缘的孔压波动大于基础中心的原因是:在计算中

25、设计了基础底面不透水,而其他地面是自由透水边界,基础边缘是边界条件突变处,所以在 靠近地表处的 -2 m左右产生较大的误差。由此可见, 在 EFGM 求解固结方程的数值分析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 积分结构 (b 40个 cell 积分结构(c 100个 cell 积分结构第 23卷 第 18期 张延军等 . 无网格伽辽金法求解固结方程的数值误差分析 3121·中, cell 积分结构的精细程度对计算的误差和稳定 性有巨大影响。在有条件的情况下,可通过增加积 分结构数量,达到提高初始孔压计算精度的效果。 图 4 条形基础不同积分 cell 结构的孔压计算曲线 Fig.4 Calculation curves of pore-pressure with differentintegration

27、cells in stripe foundation5 结 论本 文 介 绍 了 EFGM 的 发 展 概 况 , 推 导 了Terzaghi 固结方程 EFGM 的离散形式,并对 EFGM 数值计算固结中初始孔压的误差和稳定性进行了研 究,得出如下结论:(1 由于 EFGM 的理论特点,其数值近似和稳 定性较传统 FEM 复杂, 本文首次给出了固结 EFGM 离散方程的误差分析不等式。(2 建立了类似于 FEM 的 EFGM 时间稳定因子 的数学公式,并给出了全新的公式指标和系数。(3 探讨了 EFGM 中影响域大小对 EFGM 计算 精度的影响规律。(4 一、二维数值试验表明,在 EFGM

28、 中积分 cell 结构的精细程度对 EFGM 计算初始孔压有巨大 影响。(5 EFGM 作为新的数值方法,其独特的计算 优势已经在计算力学的各个领域得到证实,深入研 究其计算精度和稳定性,将为该方法的推广和应用 提供理论基础。参 考 文 献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