基于蒙特卡罗积分的无网格伽辽金法及其在渗流分析中的应用.pdf

第 29 卷 第 12 期 岩 土 工 程 学 报 Vol.29 No.12 2007 年 12 月 Chinese Journal of Geotechnical Engineering Dec., 2007 基于蒙特卡罗积分的无网格伽辽金法及其在渗流 分析中的应用 介玉新 1, 2，刘 岩1, 2 （1. 清华大学水利水电工程系，北京 100084；2. 水沙科学与水利水电工程国家重点实验室清华大学，北京 100084） 摘 要无网格法在处理有限元法难以解决的问题时具有显著的优势，且前后处理比较简单。无网格伽辽金法 （Element-free Galerkin ，EFG）是无网格法的一种，它采用滑动最小二乘法近似场函数，计算精度较高，且有 较好的稳定性，在结构和渗流分析中受到欢迎。但 EFG 方法需要背景积分网格进行高斯积分，离真正的“无网格”方 法还有一定距离。对其积分方法进行改进，采用蒙特卡罗方法（Monte Carlo）进行积分运算，由此提出了基于蒙特卡 罗积分的无网格法（MCEFG）。从而摆脱 EFG 方法对背景积分网格的依赖，使之成为真正意义上的无网格方法。基于 MCEFG 方法的渗流分析程序对有自由面的渗流问题进行了较好模拟。 该方法可以推广应用到其他需要背景积分网格的 方法中，从而使无网格法真正与网格脱离。 关键词无网格法；无单元法；蒙特卡罗方法；数值积分；渗流 中图分类号O302；O211 文献标识码A 文章编号1000–4548200712–1794–06 作者简介介玉新1970– ，男，副教授，博士，主要从事土工数值分析、基础工程、土工合成材料等方面的教学和 科研工作。E-mail jieyx。 Element-free Galerkin based on Monte Carlo integration and its application in seepage analysis JIE Yu-xin1, 2，LIU Yan1, 2 1. Department of Hydraulic Engineering, Tsinghua University, Beijing, 100084, China; 2. State Key Laboratory of Hydroscience and Engineering, Tsinghua University, Beijing 100084, China Abstract Meshless s were preponderant in dealing with problems when it was difficult to employ finite element . And they could simplify pre- and post-processings. Element-free Galerkin EFG, based on moving least squares MLS, was one of the meshless s. Owing to its accuracy and stability, EFG was widely applied in structure and seepage analysis. However, EFG was not a “truly” meshless , since mesh was still used to conduct numerical integration. A Monte Carlo-based EFG MCEFG was proposed to do the integration with Monte Carlo . The mesh or cell was not used in MCEFG . So it was truly a meshless . A program based on MCEFG was compiled and applied to seepage analysis. It was shown that the was suitable for simulating seepage with a free surface. The idea in this study could be generally extended to other s with integral mesh, and transed them to be truly meshless. Key words meshless ; element-free ; Monte Carlo ; numerical integration; seepage 0 引 言 有限元法是目前工程问题数值模拟的主要手段。 但在解决大变形、变动边界等问题时，无网格法具有 显著的优势。另外，无网格法摆脱了常规有限元方法 中将域划分为单元的限制，能够灵活地布置节点，前 后处理比较简单。随着计算机硬件和软件的发展，工 程计算的发展方向是只需给出结构尺寸、材料性质和 边界条件，就可以直接进行计算，无需用户进行网格 剖分。与有限元方法相比，无网格法在适应这方面需 要上也更有竞争力。 目前已提出了很多种无网格方法。它们之间的主 要区别在于所使用的试探函数（trial function，如滑动 ─────── 基金项目国家重点基础研究（973）资助项目2007CB714102 收稿日期2007–01–17 第 12 期 介玉新，等. 基于蒙特卡罗积分的无网格伽辽金法及其在渗流分析中的应用 1795 最小二乘、重构核函数、单位分解、径向基函数、点 插值等）和微分方程的等效形式。在积分方案上，有 的可以通过节点直接积分， 但这种积分方案不太稳定， 有的依赖背景网格进行积分。 无网格法的研究历史可追溯到 20 世纪 70 年代， 一般认为无网格法的发展开始于光滑粒子流体动力学 （smoothed particle hydrodynamics，SPH）方法[1-3]的 提出。1992 年 Nayroles 等[4]将滑动最小二乘近似 （moving least square，MLS）引入 Galerkin 法中，提 出了漫射元法（diffuse element ，DEM） 。 Belytschko 等[5]对 DEM 进行了改进，提出了无单元 Galerkin 法（element-free Galerkin ，EFG） 。随 后，许多学者开始了无网格法的研究。 EFG 方法是无网格法的一种。在国内，寇晓东对 EFG 方法进行了较为系统的研究，并将其用于追踪结 构开裂计算[6]。黄岩松将其推广用于三维结构计算[7]。 葛锦宏[8]首先将其用于渗流分析，计算有自由面的渗 流问题。周晓杰将其用于管涌问题的模拟[9]。 EFG 方法在国内有不同的翻译名称，主要有“无 单元法” 、 “无单元伽辽金法” 、 “无网格伽辽金法”等。 由于该方法仍然需要划分积分网格，所以称为“无单 元伽辽金法”可能更恰当。但很多中文文献中将其称 为 “无网格伽辽金法” ， 以及考虑到与无网格法这一大 类方法的衔接关系， 所以在此也将 “无网格伽辽金法” 作为 EFG 的中文名称。 EFG 方法主要利用背景网格进行数值积分。由于 网格和节点相互独立，积分网格的构造比有限元方法 更为方便。但它毕竟没有摆脱对网格的依赖，虽然可 以通过节点直接积分，但这种积分方案不太稳定，难 以推广应用。 为了解决 EFG 方法对积分网格的依赖， 本文采用 蒙特卡罗方法（Monte Carlo）进行积分运算。蒙特卡 罗方法即随机模拟方法（random simulation） ，亦称作 随机抽样技术（random sampling）或统计试验方法 （statistical testing） 。由于随机模拟的原型来自博采， 人们就用博采之都 Monte Carlo 作为其别称。蒙特卡 罗方法计算手段及程序结构简单，在数学上主要用于 高维多重积分等常规方法难以实现的积分计算[10-11]。 在有限元计算中由于高斯积分等的成功而使其他方法 受到忽视。本文将蒙特卡罗方法引入 EFG 的积分计 算，由此提出了基于蒙特卡罗积分的无网格法 （MCEFG） 。该方法相当于用一系列随机分布的点集 代替积分域进行积分计算。只要将求解域用两套点集 进行离散就能进行数值计算。一套为节点，用于构造 形函数；一套为随机分布的积分样本点，用于数值积 分。这种方法彻底摆脱了对背景积分网格的依赖，从 而使之成为真正意义上的无网格方法。 笔者编写了基于 MCEFG 方法的计算程序，并将 其用于渗流分析，对有自由面的渗流问题进行了数值 计算。计算表明，该方法能够满足岩土工程计算精度 的要求。 本文的研究在思路具有普遍性，可以用于其他需 要背景网格的无网格方法中。即使对于无需背景网格 的方法，如果考虑将积分点与节点分离，即节点只用 于构造插值函数，积分计算则通过另外构造的一套点 集进行，也可能会为解决计算不稳定等问题提供新的 思路。 1 EFG 方法简介 假定场函数为u X，其中 T , x yX为场点。 ,1,2,..., ii inu Xu为场函数在n个给定场点上 的已知值。滑动最小二乘的插值函数为[6,8,12] 1 n ii i n ∑Gu XX u ， 1 其中， i n X为i节点的形函数在X点的取值。 1 m iiji j nwc ∑ XXX ， 2a ,, jji ji j qq c b X XXX X X ， 2b 2 1 , n jiji i bwq ∑XXXX ， 2c 其 中 ， i w X为 权 函 数 ， j q为 正 交 基 函 数 1,2,...,jm。实践表明用正交基不但可以减少计算 工作量，而且可以改善插值计算的精度。 EFG方法与常规有限元方法的主要不同之处是 形函数的不同。在积分方案上通常都采用高斯积分的 形式。 2 渗流分析中的 EFG 方法 在渗流分析中，自由面的计算是难点。用有限元 法处理这类问题比较困难。因为自由面位置事先并不 知道，需要根据求解结果来确定。传统上一般是先假 定一个自由面，然后根据计算结果进行调整。一般需 要通过迭代来实现，有限元网格随自由面位置的变动 而调整[13]。该方法的缺点是每次迭代都要更新有限元 网格，如果假定的自由面与实际相去较远，可能会导 致单元畸形。 尤其在自由面附近单元材料类型较多时， 处理起来就更为困难。为此人们发展了固定网格法， 如Desai于1976年提出了剩余流量法[14]。 也有采用饱 和－非饱和渗流理论来处理问题的方法[15]。固定网格 法改善了计算效果，但也存在迭代计算上的困难。笔 1796 岩 土 工 程 学 报 2007 年 者也曾尝试用适体坐标变换方法求解自由面问题 [16-17]，但它是一种差分方法，对流量边界的处理比较 困难，对水位变动情况下自由面的计算精度较低，所 以也不太适合这类问题。由于无网格法能够灵活地布 置节点，特别适用于边界条件随计算结果而调整的情 况，所以在自由面的计算中是有优势的。 利用变分原理或加权余量（残值）法可以写出二 维饱和稳定渗流的有限元格式 [ ]{ } { }0KHF ， 3 其中 T [][ ] [ ][ ]d∫∫ Ω KBk BΩ ， 4a Ω Ω d y n y n k x n x n kK j i y j i xij ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∫∫ ，4b d ii Fqn ∫ 2 Γ Γ 。 5 上式中的形函数可以通过式（2）构造。 3 EFG 的积分方案 虽然与有限元方法采用相似的网格进行积分计 算，但EFG积分网格的布置比较灵活，只要能够将域 离散即可。一种方案是采用如图1所示的矩形网格进 行积分[6,8]；一种是采用常规与有限元类似的网格进 行，此时网格仅用于积分计算，与式（1）中的插值节 点则没有必然联系。 图 1 背景网格积分 Fig. 1 Background integration mesh 图2是采用有限元网格进行积分的例子[8]。由于 网格只是用来离散积分区域，所以不必象有限元那样 严格要求网格之间的邻接关系。 还有一种方法是利用节点进行积分，但其计算稳 定性较差[1]。 4 蒙特卡罗积分 蒙特卡罗方法是一种利用统计抽样理论近似求解 数学问题或物理问题的方法。它的主要理论基础是概 图 2 有限元网格积分 Fig. 2 Integration mesh following FEM 率论中的大数定律，主要手段为随机变量的抽样[18]。 它的主要优点是方法及程序简单，适应性强，收敛的 概率性和收敛速度与问题维数无关。在数学上主要用 于高维多重积分的计算。而在有限元计算中，由于高 斯积分的成功使蒙特卡罗方法被忽视。另外，早期由 于计算机运行速度限制，不能进行大量随机抽样，其 实用性受到制约。随着计算机硬件的发展，这种方法 应当重新受到重视。 基于蒙特卡罗的积分方法有频率法 （随机投点法） 和期望法 （平均值法） ， 以及一些通过降低方差提高精 度的方法。本文采用期望法来计算积分。这种方法比 较简单，易于编程。 期望法的核心是把积分看成某个随机变量的期 望。 以一维积分为例， 最自然的是看成积分域[ , ]a b上 均匀随机变量的期望。设 U[ , ]a bη（[ , ]a b上的均 匀分布） ，则 d b a If xxba Efη− ∫ 。 6 于是对于N个独立的[ , ]a b上的均匀随机数，可 以用矩估计 1 N ff Iba N ηη − L 7 作为 d b a If xx∫的估计。显然，它是无偏估计[11]。 上述蒙特卡罗求积方法是以随机变量θ的简单子 样 12 ,,..., N θ θθ的算术平均值 1 1 N i i N θθ ∑ 8 作为所求积分I的近似值。如果随机变量序列 12 ,,..., N θ θθ相互独立、同分布、期望值存在，则由大 数定理可知，当N → ∞时，θ以概率1收敛于I。同 时，根据中心极限定理，若有限标准差0σ≠，那么 当N → ∞时，随机变量//YINθσ−渐进地服 从标准正态分布0,1N，即[10] 第 12 期 介玉新，等. 基于蒙特卡罗积分的无网格伽辽金法及其在渗流分析中的应用 1797 2 2 1 ed 2π x t r P Ytx α α − −∞ ， 有 r P Ytα r t PI N ασ θ− 2 2 0 2 ed1 2π x t x α α − ≈ − ∫ ，这表示不等式 t I N ασ θ− 。 10 成立的概率近似等于1α−。 其中， 当α很小时，α称 为显著水平，1α−称为置信水平。因此，蒙特卡罗方 法的误差ε可以写为 t N ασ ε 。 11 该误差是在概率意义下的统计估值。 模拟结果误差 ε 与标准差σ成正比，与模拟次数的平方根N成反 比。在许多情况下，无法知道理论值σ，因此对于本 文采用的平均值法，在实际计算中可以用标准样本方 差来代替σ。从式（11）可以看出，只要N足够大， 能够保证误差ε满足精度要求。 为了验算上述方法，首先对三个算例进行计算 （a） 1 0 dx x ∫ ，（b） 11 00 d dxy x y ∫ ∫ ，（c） 22 1xy Ω −− ∫∫ d dx y（积分域Ω为单位圆） 。 蒙特卡罗数值积分计算结果及其与精确解的相对 误差列于表1中。计算采用的随机样本点个数为 320000。 可以通过Fortran自带的随机数生成器先在区 间[0, 1]上生成样本点，然后变换到相应域上即可。表 中计算值为10次计算的平均值。从计算结果可以看 出，蒙特卡罗数值积分的计算精度能够得到保证。 表 1 蒙特卡罗数值积分计算结果 Table 1 Integration results with Monte Carlo 蒙特卡罗积分 积分式 精确解 计算值 误差/ 1 0 dx x ∫ 2 3 0.666689980.003 11 00 d dxy x y ∫ ∫ 0.25 0.24983280.067 22 1d dxyx y Ω −− ∫∫ 2.0943951 2.09501520.030 5 基于蒙特卡罗积分的 EFG 方法 （MCEFG）及其应用 采用蒙特卡罗方法计算式（4）所示的积分，即可 得到基于蒙特卡罗数值积分的无网格Galerkin方法 （MCEFG） 。在实际应用中，如图3所示，为了避免 计算不规则积分域Ω的面积， 可直接在包含积分域Ω 的规则矩形域内布点， 假设积分点的总数为M， 落入 积 分 域 内 的N个 积 分 点 对 应 的 函 数 值 为 1122 ,,,,...,, NN f x yf xyf xy，则有 ∑∫∫ Ω ≈ N i ii yxf M S yxyxfI 1 ,dd, ， 12 其中，S为包含积分域的规则矩形域的面积。 这种处理对于不规则积分域是方便的，如果积分 域为梯形等容易计算其面积的情况，可以直接进行面 积计算。 图 3 积分域的处理 Fig. 3 Integration domain 为了验证MCEFG方法的适用性，首先对如图4 （a）所示的均质矩形土坝进行计算，渗透系数 0.1 m/d xy kk，上游水位6 m，下游水位1 m， 计算时假定的初始自由面位置如图4（a）中的直线。 采用MCEFG法对其进行渗流计算，最终得到的自由 面位置如图4（b）所示（图中的黑点是计算结束时节 点的位置） 。采用该方法与传统EFG法的结果对比如 表2。 可以看出，MCEFG法的计算结果能够满足要求。 图 4 MCEFG法计算矩形土坝 Fig. 4 Calculation of a rectangular dam by MCEFG 图2是传统EFG方法采用的积分网格[8,19]。需要 1798 岩 土 工 程 学 报 2007年 在计算前事先完成。MCEFG方法无需积分网格，只 需在计算域内通过Fortran自带的随机数生成器布置 随机分布的积分点集即可。它相当于用大量的点来离 散积分区域。 表 2 自由面节点水头 Table 2 Water head of the node on free surface 单位m x 坐标 1.0 2.0 3.0 4.0 试验值 5.63 5.10 4.38 3.25 有限元 5.65 5.22 4.44 3.53 EFG 5.65 5.14 4.45 3.34 MCEFG 5.65 5.12 4.42 3.36 进一步对如图5所示的均质土坝进行计算，土坝 坡度1∶2， 渗透系数0.1 m/d xy kk， 上游水位6 m，下游水位处于基准面位置。假定的初始自由面位 置如图5（a）所示。采用MCEFG法对其进行渗流计 算，最终得到的自由面位置如图5（b） 。计算得到的 溢出点坐标与采用传统EFG方法的计算结果对比见 表3。 图 5 MCEFG法计算均质土坝 Fig. 5 Calculation of a homogeneous dam by MCEFG 表 3 溢出点坐标的比较 Table 3 Calculated location of seepage exit 单位 m 计算方法 溢出点坐标 EFG （32.05，1.98） MCEFG （31.89，2.06） 可以看出，在工程意义上，MCEFG的计算精度 是能够满足要求的。MCEFG方法最大的优点是前处 理比较简单，只要给出结构或求解域的几何特征、材 料参数和边界条件，采用合适的自动布点技术[19]就可 以直接进行计算，无需划分网格、检查网格特性等繁 杂的前处理工作。 用蒙特卡罗方法进行积分的主要缺点是收敛速 度慢，计算精度与N/1有关，所以为提高精度需要 进行大量运算。但随着计算机运行速度的提高，大规 模运算已逐渐成为可能。此外，在数学上已有一些通 过降低方差来提高计算精度的方法，但如何用于EFG 方法有待进一步研究。相信随着研究的深入，结合无 网格法工程计算特点构造合适的改善精度的算法，并 采用相应的权函数，提高计算效率是不难实现的。 6 结 论 （1）EFG方法计算精度和稳定性较高，但它需 要背景网格进行积分计算。本文对其积分方法进行改 进，采用蒙特卡罗方法（Monte Carlo）进行积分运算， 提出了基于蒙特卡罗积分的无网格法（MCEFG） 。从 而摆脱EFG方法对背景积分网格的依赖， 使之成为真 正意义上的无网格方法。通过算例计算表明MCEFG 方法能够满足工程计算精度的需要。 （2）MCEFG方法在思路上的特点在于采用了两 套点集来离散求解域一套是插值节点，用于构造形 函数；一套是积分样本点，用于离散积分区域。两者 彼此独立。 该方法对点集的质量要求较低， 性质稳定， 而且自动布点技术的实现也比划分有限元网格容易得 多。这就使得MCEFG方法编程简单，自动化程度较 高。 （3）这种“以点代面”的思路可以推广用于其 他无网格方法。对于需要背景网格的无网格法，采用 Monte Carlo积分能够使其真正与网格脱离。 即使对于 无需背景网格的方法， 如果考虑将积分点与节点分离， 即节点只用于构造插值函数，积分计算则通过另外构 造的一套点集进行，也会为解决计算不稳定问题提供 新的可能。 （4）Monte Carlo方法的缺点是收敛速度慢，计 算量偏大。不过，随着技术进步，计算机硬件运行速 度的迅速提高能够抵消计算方法在效率方面的不足。 此时人们会更关注方法本身的简便性和稳定性。在简 便性方面，Monte Carlo方法是具有竞争力的。结合无 网格法计算特点对算法进行优化，通过降低方差等构 造精度更高的Monte Carlo积分方法，也能够改善其 计算效率。 （5）随着计算机技术以及以有限元为代表的各 种计算手段的发展，数值计算已深入人心。在一些特 定问题上，无网格法将会是有限元法的重要补充。另 一方面，计算机发展也使一些原来被人们忽视的算法 重新显现出活力，这也为科研工作提供了新的思考问 题的方式。 参考文献 [1] 张 雄, 刘 岩. 无网格法[M]. 北京 清华大学出版社, 2004. ZHANG Xiong, LIU Yan. Meshless [M]. Beijing Tsinghua University Press, 2004. in Chinese 第12期 介玉新，等. 基于蒙特卡罗积分的无网格伽辽金法及其在渗流分析中的应用 1799 [2] LIU G R, LIU M B. 光滑粒子流体动力学 一种无网格粒子 法[M]. 韩 旭, 杨 刚, 强洪夫, 译. 长沙 湖南大学出 版社, 2005. LIU G R, LIU M B. Smoothed particle hydrodynamics a meshless particle [M]. HAN Xu, YANG Gang, QIANG Hong-fu, trans. Changsha Hunan University Press, 2005. in Chinese [3] LUCY L B. A numerical approach to the testing of the fission hypothesis[J]. The Astron J 1977, 812 1013–1024. [4] NAYROLES B, TOUZOT G, VILLON P. Generalizing the finite element diffuse approximation and diffuse elements[J]. Comput Mech, 1992, 10 307–318. [5] BELYTSCHKO T, LU Y Y, GU L. Element-free Galerkin s[J]. Int J Num Meth Engng, 1994, 37 229–256. [6] 寇晓东. 无单元法追踪结构开裂及拱坝稳定研究[D]. 北京 清华大学, 1998. KOU Xiao-dong. Element-free tracing structure fracture and arch dam stability analysis[D]. Beijing Tsinghua University, 1998. in Chinese [7] 黄岩松. 三维无单元伽辽金法及其在坝工破坏分析中的应 用[D]. 北 京 清 华 大 学, 2005. HUANG Yan-song. Threedimensional element-free and its application in dam analysis[D]. Beijing Tsinghua University, 2005. in Chinese [8] 葛锦宏. 无单元法渗流计算及其在长江堤防中的应用[D]. 北京 清华大学, 2002. GE Jin-hong. Seepage analysis with element-free and its application in embankment of Yangtze River[D]. Beijing Tsinghua University, 2002. in Chinese [9] 周晓杰. 堤防的渗透变形及其发展的研究[D]. 北京 清华 大学, 2006. ZHOU Xiao-jie. Research on the generation and evolution of seepage deation in levee[D]. Beijing Tsinghua University, 2006. in Chinese [10] 徐钟济. 蒙特卡罗方法[M]. 上海 上海科学技术出版社, 1985. XU Zhong-ji. Monte Carlo [M]. Shanghai Shanghai Science Technology Press, 1985. in Chinese [11] 龚光鲁, 钱敏平. 应用随机过程教程及在算法和智能计算 中的随机模型[M]. 北京 清华大学出版社, 2004. GONG Guang-lu, QIAN Min-ping. Applied random process and random model in algorithm and intelligent computation[M]. Beijing Tsinghua University Press, 2004. in Chinese [12] LU Y Y, BELYTSCHKO T, GU L. A new implementation of the element-free Galerkin [J]. Computer s in Applied Mechanics and Engineering, 1994, 113 397–414. [13] DESAI C S. Finite element s for flow in porous media[M]. UK John Wiley and Sons, 1975. [14] DESAI C S. Finite element residual schemes for unconfined flow[J]. International Journal for Numerical s in Engineering, 1976, 10 1415–1418. [15] LAM L, FREDLUND D G. Saturated-unsaturated transient finite element seepage model for geotechnical engineering[J]. Advances in Water Resources, 1984, 7 132–136. [16] 介玉新, 揭冠周, 李广信. 用适体坐标变换方法求解渗流 [J]. 岩土工程学报, 2004, 261 52–56. JIE Yu-xin, JIE Guan-zhou, LI Guang-xin. Seepage analysis by boundary- fitted coordinate transation [J]. Chinese Journal of Geotechnical Engineering, 2004, 261 52–56. in Chinese [17] 揭冠周. 基于适体坐标变换的有限差分法在渗流分析中 的应用[D]. 北京 清华大学, 2002. JIE Guan-zhou. Application of boundary fitted coordinate transation in seepage analysis[D]. Beijing Tsinghua University, 2002. in Chinese [18] 高明坤. 实用概率统计学（概率及其应用）[M]. 北京 国 防工业出版社, 1988. GAO Ming-kun. Practical probability statistics probability and its application[M]. Beijing Defence Industry Press, 1988. in Chinese [19] 刘 岩. 基于模拟退火的约束布点技术及其在渗流分析 中的应用[D]. 北京 清华大学, 2006. LIU Yan. Point sets generation with prescribed boundary points based on simulated annealing and its application in seepage analysis[D]. Beijing Tsinghua University, 2006. in Chinese