NEM与IEM耦合方法在地下工程中的应用研究＊

李晓龙 王复明 徐 平 李晓楠

（郑州大学交通运输工程系1） 郑州 450002） （中原工学院计算机学院2） 郑州 450007）

自然单元法（natural element method，NEM）是一种新型的偏微分方程数值解法，由J.Braun和 M.Sambridge于1995年提出［1］.由于该方法兼具无网格法和有限元两者的优点，非常适合于求解涉及到大变形、裂隙或节理扩展、多种支护材料并且需要模拟分步施工过程的岩土及地下工程问题，因而受到了岩土工程学者的广泛关注，并展开了大量的相关研究：晏汀将NEM用于处理弹性力学问题［2］；蔡永昌研究了NEM用于地下工程计算的网格自动生成技术，并将其补充到“曙光”软件中，应用于隧道、地铁、基坑等岩土工程的线弹性分析［3－4］；朱合华实现了二维 NEM 基于Von－Mises屈服准则的弹塑性分析［5］.

尽管自然单元法在岩土工程中的应用研究已取得了不少成果，然而在处理隧洞或基坑等岩土工程无限域或半无限域问题时，仍存在与有限元方法类似的问题，即需要“人为”截取一定区域并设置相应边界条件，这种处理方式在理论上必然引起误差.无限单元法（infinite element method，IEM）是20世纪70年代为解决无限域问题而发展起来的一种数值方法，最初主要用于处理水表面波的传播问题［6］，后来被成功应用于岩土工程中［7］.为弥补NEM在处理岩土及地下工程问题时存在的不足，引入IEM模拟无穷远处边界条件，与NEM相结合形成耦合分析方法（the coupling analysis of NEM and IEM，NEM－IEM），利用深埋圆形隧洞算例对算法的正确性和对计算精度的改善作用进行检验，并探讨了计算范围选取对纯自然单元法和耦合方法分析结果的影响.

1 自然单元法的原理

自然单元法的主要实现步骤为［8］：对分析区域作结点离散；搜索积分点的自然邻接点，利用自然邻接点计算积分点的插值形函数；对计算区域积分得到结构总体刚度阵；形成等效结点荷载列阵；求解线性方程组.

1.1 自然邻接点的寻找

在搜索积分点的自然邻接点时，主要涉及到Delaunay三角化和Voronoi结构两个概念.

Delaunay三角化，即利用离散结点生成Delaunay三角形，它具有两个重要性质：（1）最大最小角性质，即在给定结点所有可能生成的三角形中，Delaunay三角形具有最大的最小内角；（2）空圆准则.通过Delaunay三角形3个顶点的外接圆里不包含其他顶点.

Voronoi结构，其定义为到点p的距离小于到其他任何结点xi的距离的集合，如图1所示，数学表达式为

式中：d（x，p）为结点x到p的距离.

将结构作结点离散后，对于分析区域上每一个按一定规则分布的计算点p，用Delaunay三角化的空圆准则求点p周围的自然邻接点：如果Delaunay三角形△l的外接圆圆心cl到p的距离d小于△l的外接圆半径R，则△l为点p的自然邻接三角形，其结点为点p的自然邻接点.

1.2 Laplace插值

Laplace插值，又称 Non－Sibson插值［9］，是一种新的自然邻点插值方式，与较早的Sibson插值相比，Laplace插值的优越性主要体现在：（1）直接利用结点的1阶Voronoi单胞的边长和点到Voronoi边的距离计算插值基函数，构造方式简单，避免了Sibson插值所采用的Waston面积叠加算法需反复进行迭代的麻烦，计算量大大降低；（2）Sibson插值在凸区域的边界是线性精确的，但对于凹区域的边界，插值并不精确；而Laplace插值对于非凸区域边界也是线性精确的，因此可以准确地施加本质边界条件.

图1 点p的自然邻接点及其一阶Voronoi结构

在图1所示的在二维空间中，点1，2，3，4为搜索得到的积分点p的所有自然邻接点.设三角形△p21，△p23，△p34，△p41的外接圆圆心分别为c21p，c23p，c34p，c41p，连接各圆心即构成点p的1阶Voronoi结构.分别定义s1，…，s4为点p的Voronoi结构的各边边长，h1，…，h4为各自然邻接点到点p的距离.则对点p的任一自然邻点i，其Laplace形函数定义为

从而点p的位移函数可表达为

式中：ui（i＝1，…，n）为点p周围自然邻结点i的结点位移；φi（x）为对应结点的插值形函数.

由式（2）知，形函数φi（x）不但满足单位分解条件，而且与有限元形函数一样满足

形函数φi（x）的导数φi，j（x）为

显然，对于NEM，只要分析区域的离散点给定，结点的三角化和自然邻接点的搜索可利用Delaunay准则由计算机自动完成，避免了FEM在网格划分和重构等方面的困难；同时，由于其插值形函数在边界上是线性精确的且满足Kronecker条件，因此可以像FEM那样施加准确的边界条件及模拟材料的不连续问题.可以说，NEM将无网格法［10］和FEM的优点巧妙地结合在一起，在数值计算中具有独特的优势.

1.3 平衡控制方程及其数值积分

利用自然邻接点采用Laplace插值得到近似位移表达式后，可按照有限元法同样的步骤，利用最小势能原理或Galerkin过程得到离散求解的平衡方程.对于线弹性平面问题对总势能取驻值得到系统的整体离散方程为

式中：D为平面问题的弹性矩阵；t为面力；fb为体力；Bi为积分点上的应变矩阵.

对式（6）的积分通常采用Delaunay三角形内的高斯积分，也可用与无网格伽辽金法（EFG）一样的规则矩形背景积分网格，或者用点积分［11］.

2 无限单元法的原理

无限单元法（infinite element merhod，IEM）的概念由Bettess和Zienkiewicz于1977年首次提出，其目的是为了弥补有限元在处理无界域问题上存在的不足.根据是否有局部坐标与整体坐标的映射关系，可分为映射无限元和非映射无限元.目前在岩土工程中，映射无限元的应用最普遍，其原理是利用映射函数将无限区域从几何上映射为一有限区域，然后对该区域作常规有限元类似的分析.

为了便于实现和自然单元的耦合，本文采用二维四节点映射无限元，其形式如图2所示.

图2 二维四节点无限元

坐标变换式为

式中：Ni为坐标变换的插值函数，其形式如下.

当η≤0时，有

当η＞0时，有

式中：ξi，ηi为节点i的局部坐标值.显然式（10）、（11）是连续的，当η→1时，对应的边界趋于无穷远.位移表达式为

式中：Mi为位移插值函数，为满足无限远处位移为零的边界条件，可选如下形式

式中：取η≤0时Ni的表达式；f（ri／r）称为位移衰减函数；r为衰减半径，指积分点到衰减中心的距离；ri为节点i的衰减半径，即节点i到衰减中心的距离；在衰减半径趋于无穷大时使位移衰减函数f（ri／r）趋于零，可令

可以看出，四结点无限元由于在坐标变换和位移模式中分别采用了特殊构造且形式不同的插值函数，实现了从局部坐标系中有限域到整体坐标系下无限域的映射以及无限远处位移为0的边界条件.

3 NEM与IEM的耦合

由于自然单元法的插值函数满足Kronecker条件并在边界上满足线性插值，可以很容易与线性无限元实现“无缝”耦合.设求解域Ω＝ΩNEM＋ΩIEM，则近似位移表达式为

式中，ΩNEM为自然元区域；ΩIEM为无限元区域.

可见由于两者的插值函数在交界处满足线性插值，确保了位移的连续性，因此在其交界面上无需特殊处理，不必像常规无网格法那样需要设置一个中间过渡区域或修改插值函数，实现方式十分简单.

4 算 例

某深埋地下洞室，如图3所示，开挖断面为圆形，半径为R＝3m，洞室周围的静水压力为p0＝5MPa，假设围岩为连续、均质及各向同性体，弹性模量E＝2GPa，泊松比μ＝0.25.

图3 圆形洞室及边界条件示意图

由于结构对称，取1／4区域作分析，沿环向离散为等间距的9个结点，径向结点间距按比例逐渐增大；为比较计算范围对计算精度的影响，按半径为8.9，10.8，13.2，16，19.5，23.8，29m 这7种情况分别用NEM和耦合方法作线弹性分析，并与理论值作比较.

不同情况下NEM和耦合方法所得结点径向位移沿半径分布曲线分别如图4和5所示.从中可以看出，自然单元法得到的结点径向位移计算值随计算范围的增大逐渐趋近于理论解，而当计算范围减小时，两者的偏差迅速变大；与自然单元法不同的是，耦合方法的分析结果虽然也受到计算区域选取的影响，但影响程度很小，随着计算范围的变化，径向位移只发生轻微波动，始终与理论解非常接近.

图4 径向位移沿半径的分布（NEM）

图5 径向位移沿半径的分布（NEM－IEM）

图6 径向应力沿半径的分布（NEM）

图7 径向应力沿半径的分布（NEM－IEM）

图8 环向应力沿半径的分布（NEM）

图9 环向应力沿半径的分布（NEM－IEM）

NEM和耦合方法得到的径向和环向应力沿半径分布曲线分别如图6～9所示，可以看出，两种方法得到的径向及环向应力计算值与位移解具有相似特点：即耦合方法的分析结果受计算范围影响小，求解精度很高，计算曲线与理论曲线几乎重合；而纯自然单元法的精度较低，当计算区域减小时，误差增加很快.

为了更具体地比较两种方法的精度差别，选取3种典型情况下耦合与非耦合方法得到的结点径向位移列于表1中，可以发现，耦合方法的计算精度远高于非耦合方法，当计算范围均取29m时，前者得到的隧洞内表面径向位移误差仅为0.34%，而非耦合方法误差则高达3.3%，两者计算精度几乎相差10倍；值得注意的是，对于耦合方法，当计算范围仅取8.9m时，其精度仍然远高于非耦合方法29m时的计算结果，而此时前者的结点总数仅为54，后者却达到108，那么据此结果，在满足工程要求的情况下，采用耦合方法显然可以大大缩小结构刚度矩阵的规模，从而减小计算工作量，这无论对于岩土工程正向分析还是反演都具有重要意义.

表1 径向位移值的比较

5 结 论

1）在处理无限域或半无限域问题时，采用自然元与无限元耦合分析方法，能减少纯自然单元法因“人为”截取分析区域及规定边界条件而引起的误差，显著提高计算结果的精度，从而进一步提高了自然单元法处理岩土工程问题的能力.

2）分析范围的大小对耦合方法所得结果影响不显著，在满足工程要求的前提下，采用耦合方法只需选取较小的计算区域，就能得到令人满意的结果，因而可以大大降低结构刚度矩阵的规模，减小计算工作量，这对于复杂岩土工程尤其是其反演分析具有非常重要的意义.

［1］Braun J，Sambridge M.A numerical method for solving partial differential equations on highly irregular evolving grids［J］.Nature，1995，376：655－660.

［2］晏 汀，沈成武.自然单元法在弹性力学中的应用［J］.武汉理工大学学报：交通科学与工程版，2006，30（6）：1052－1054.

［3］蔡永昌，朱合华.岩土工程数值计算中的无网格法及其全自动布点技术［J］.岩土力学，2003，24（1）：21－24.

［4］蔡永昌，朱合华，夏才初.无网格自然邻接点法及其在岩土工程数值模拟中的应用［J］.岩石力学与工程学报，2005，24（11）：1888－1924.

［5］朱合华，杨宝红，蔡永昌，等.无网格自然单元法在弹塑性分析中的应用［J］.岩石力学，2004，25（4）：671－674.

［6］Bettess P，Zienkiewicz O C.Diffraction and refraction of surface waves using finite and infinite elements［J］.International Journal for Numerical Methods in Engineering，1977，11：1271－1290.

［7］周世良，胡 晓，王 江.无限元在岩土工程数值分析中的应用［J］.重庆交通学院学报，2004，23（增刊）：61－64.

［8］Sukumar N，Moran B，Belytschko T.The nature element method in solid mechamics［J］.Int J Numer Meth Engng，1998，43：839－887.

［9］Sukumar N，Moran B，Semenov A Y，et al.Natural neighbor Galerkin method［J］.Int.J.Numer.Meth.Engng.，2001，50：1－27.

［10］宋康祖，陆明万，张 雄.固体力学中的无网格方法［J］.力学进展，2000，30（3）：55－56.

［11］卢 波.自然单元法的发展及其应用［D］.武汉：中国科学院武汉岩土力学所，2005.