EFG1无网格法：原理、实现与工程实践

简介：无网格法是一种无需预先生成网格的数值计算方法，在处理自由边界、变形固体和非线性问题时比传统网格类方法更灵活。资源围绕 EFG1（无网格伽辽金法）的实现展开，包含完整的 MATLAB 函数脚本与配套数据，适合正在学习无网格法、需要参考代码来理解节点分布、插值构造与方程离散流程的学生或研究人员。压缩包共 4 个文件，其中 2 个 .m 脚本负责核心计算与节点定义，1 个 .asv 为 MATLAB 自动备份文件，另 1 个 .mat 用于存储边界或计算结果，整体仅 1KB，轻量但结构清晰。已有 266 人浏览学习。通过对照代码和数据，可以快速掌握 EFG 法的程序骨架，包括插值函数选择、系数矩阵组装的编码思路，并对后续扩展至 SPH、DEM 等其他无网格方法提供基础。

1. EFG1 无网格法：从 FEM 的网格枷锁里挣脱的第一代方案

EFG1（Element Free Galerkin）无网格法是计算力学领域里少有的、既能写论文又能直接落工程代码的数值方法。它不需要像 FEM 那样预先划分单元网格，而是用一组散乱分布的节点配合移动最小二乘（MLS）近似，直接构建形函数并完成 Galerkin 离散。这意味着你处理裂纹扩展、大变形、材料冲切这类拓扑不断变化的问题时，不用反复 remesh，整个迭代过程省掉的不只是建模工时，更是网格畸变带来的收敛性灾难。EFG1 不是最年轻的无网格格式，但它是最早把“无网格”这个词做成可计算方案的框架之一。本文面向已经在用 FEM、想引入无网格方法解决问题的工程师和研究者，从头梳理 EFG1 的数学骨架、代码实现、参数陷阱和实战验证路径。

2. EFG1 无网格法的理论基础：MLS 近似与 Galerkin 离散

2.1 为什么 EFG1 绕开网格还能构造形函数

EFG1 核心思路是：放弃单元，用移动最小二乘（MLS）在局部处处拟合节点场值。MLS 的关键动作是做“加权最小二乘拟合”，每算一个点的形函数，就以该点为中心划出一个影响半径范围内的邻近节点，对这些节点上的已知值做局部近似，这种局部拟合优度用带权重的 L2 范数度量，权重函数随距离衰减，导数连续性和紧支性由权重函数控制。形函数不是显式给出，而是通过求解一个小型线性方程组来数值确定。

MLS 形函数的数学表达为：

取基函数向量 ( p(x) = [1, x, y, x^2, xy, y^2] ) （二维二次基底），在每个评估点 ( x ) 处求系数向量 ( a(x) ) 使得

[ J(x) = \sum_{I} w(x-x_I) \left( p^T(x_I) a(x) - u_I \right)^2 ]

最小化。解出 ( a(x) ) 后，形函数 ( \Phi_I(x) = w(x-x_I) \cdot p^T(x) A^{-1}(x) p(x_I) )。其中矩阵 ( A(x) = \sum_I w(x-x_I) p(x_I) p^T(x_I) ) 是力矩矩阵，维度是基函数个数的平方，二维二次基下是 6x6。这个方程每个评估点都需要求解一次，计算量比 FEM 的单元形函数大一到两个数量级，这也是无网格方法最直接的代价。

2.1.1 权重函数的选择与影响

权重函数是整个 EFG1 近似质量的中枢。常用的是三次样条权重和四次样条权重。三次样条写作：

w(r) = 2/3 - 4r^2 + 4r^3 (0 ≤ r ≤ 1/2) 4/3 - 4r + 4r^2 - 4/3r^3 (1/2 ≤ r ≤ 1) 0 (r > 1)

其中 ( r = |x - x_I| / s )，s 为该影响半径。三次样条的优点是导数在 r=0 和 r=1 处连续，但二阶导在 r=1/2 处有跳变。四次样条则额外保证了二阶导连续，收敛性更好，代价是更多的浮点计算。我一般优先用四次样条，因为 EFG1 在计算刚度矩阵时要求形函数导数，导数高阶连续性直接影响应力场的平滑度。

权重函数影响半径 s 的取值原则：二维问题中，每个节点的影响半径通常取该节点到最近邻节点距离的 2.5 到 4.0 倍。太小则覆盖的邻居节点不够、力矩矩阵 A 奇异，太大则拟合过于光滑，局部特征被磨平。

2.2 Galerkin 离散与刚度矩阵装配流程

EFG1 的离散过程和 FEM 同构：把试探函数和检验函数都投影到 MLS 形函数张成的空间里。以二维弹性静力学为例，位移场 ( u(x) ) 近似为：

[ u(x) = \sum_{I=1}^{N} \Phi_I(x) u_I ]

把该近似代入虚位移原理或者势能泛函极小化条件，得到的离散线性系统为：

[ K_{IJ} = \int_{\Omega} B_I^T D B_J , d\Omega ]

其中 ( B_I ) 是应变矩阵，由 ( \Phi_I ) 的导数组装而成，D 是材料本构矩阵。这不是单元刚度矩阵的叠加，因为每个点上的积分域是重叠的影响域。计算中，求解域内用高斯积分背景网格积分，边界上的牵引力条件通过对称罚函数法或拉格朗日乘子法施加。

2.2.1 背景积分网格的必要性

EFG1 只是“不划分单元”，但仍需要一个辅助性的背景积分网格来数值积分刚度矩阵。常见做法有三种：规则矩形网格法、四叉树自适应网格法和 Voronoi 图法。规则矩形网格最容易实现但精度受网格对齐影响；四叉树方法能自适应加密高梯度区域，适合断裂力学中的局部高应力问题；Voronoi 图法最精确但需要额外的几何计算。对大多数梁板问题，规则矩形背景网格配上自适应细分就足够了。背景网格密度：每个节点在各方向至少被 3~4 个高斯积分点覆盖，否则会出现空间振荡。二维情况下每个积分单元用 4x4 高斯点，这是个稳妥的起点。

2.3 EFG1 与 FEM 的边界条件处理差异

FEM 中本质边界条件（固定位移、强制位移）只需把节点自由度直接约束，因为形函数满足 Kronecker delta 性质，节点值就是精确位移。而 EFG1 的 MLS 形函数不满足插值性，边界节点的近似值是该节点邻域内的拟合值，并非精确节点值。直接施加本质边界条件会产生系统性误差，误差大小和权重函数的边界截断有关。常见处理方案：

工程中时间受限时用罚函数法起步，罚系数取材料弹性模量的 (10^4\sim10^6) 倍。研究精度优先时上拉格朗日乘子或者 Nitsche 法。罚函数法的致命弱点是病态条件数：罚系数过大，刚度矩阵条件数从 (10^8) 量级飙升到 (10^{14}) 以上，线性求解器直接失准。正确做法是用双精度求解并用迭代法时做残差监控。

3. 用 Python 实现一个最小可运行的 EFG1 悬臂梁求解器

3.1 节点离散与背景网格生成

求解对象选经典悬臂梁，长 L=2.0 m，高 H=1.0 m，右端受向下集中力 P，左端固定。先把计算域铺上均匀散乱节点，再叠加背景网格。

import numpy as np # 节点密度控制参数 nx = 21 # x 方向节点数 ny = 11 # y 方向节点数 L, H = 2.0, 1.0 # 产生均匀节点 x_vals = np.linspace(0, L, nx) y_vals = np.linspace(0, H, ny) nodes = np.array([(x, y) for y in y_vals for x in x_vals], dtype=float) nnode = len(nodes) # 背景积分网格：每个 1x1 单元内再细分 4 个子网格 bg_cells = [] for i in range(nx - 1): for j in range(ny - 1): xl, xr = x_vals[i], x_vals[i+1] yb, yt = y_vals[j], y_vals[j+1] # 每个背景单元细分 2x2 以提升积分精度 for ci in range(2): for cj in range(2): x1 = xl + (xr - xl) * ci / 2 x2 = xl + (xr - xl) * (ci + 1) / 2 y1 = yb + (yt - yb) * cj / 2 y2 = yb + (yt - yb) * (cj + 1) / 2 bg_cells.append((x1, y1, x2, y2))

背景网格的关键变量是单元坐标四元组(x1, y1, x2, y2)。每个背景单元内部积分点用于数值积分，单元尺寸不能超过节点间距的 1.5 倍，否则高斯点周围没有足够的节点权重覆盖。

3.2 MLS 形函数计算与导数输出

这是 EFG1 的核心函数，所有性能瓶颈都在这里。下面代码实现了 2D 一次基和二次基可切换的 MLS 求解。

def mls_shape_and_derivatives(pt, nodes, support_radius, basis='quadratic'): # 稀疏局部列表，不需要全量搜索 diff = nodes - pt dist = np.linalg.norm(diff, axis=1) mask = dist <= support_radius idx = np.where(mask)[0] d = dist[idx] d[d < 1e-15] = 1e-15 w = spline_weight(d, support_radius) xi, yi = nodes[idx, 0], nodes[idx, 1] if basis == 'linear': p = np.array([np.ones_like(xi), xi, yi]).T # 3xM dp = np.array([[0]*len(idx), np.ones(len(idx)), np.zeros(len(idx)), np.zeros(len(idx)), np.zeros(len(idx)), np.ones(len(idx))]).T elif basis == 'quadratic': p = np.array([np.ones_like(xi), xi, yi, xi*xi, xi*yi, yi*yi]).T # 6xM dp = np.array([np.zeros(len(idx)), np.ones(len(idx)), np.zeros(len(idx)), xi, yi, np.zeros(len(idx)), np.zeros(len(idx)), np.ones(len(idx)), np.zeros(len(idx)), np.zeros(len(idx)), xi, yi]).T.reshape(len(idx), 2, 6) A = (p.T * w) @ p # 力矩矩阵 # 求逆再加物理扰动避免奇异 A_inv = np.linalg.inv(A + 1e-12 * np.eye(A.shape[0])) # D = A^{-1} * p^T * W D = A_inv @ (p.T * w) phi = D.T @ p[0] # 形函数值（实际上 phi 就是 D 的第一列以外部分） # 更标准的是 phi_I = sum_j p_j(x) * D[j, I] print(phi.shape) return phi, D, idx

形函数计算结果实际使用时的组织方式和 FEM 不同。FEM 里形函数矩阵是稀疏的且和单元关联，EFG1 中每个评估点产生的是一组稠密的局部向量phi和权重矩阵D。刚度矩阵组装时要把局部贡献散射到全局自由度。注意A_inv中用1e-12扰动项是必要的数值稳定措施，但扰动太大（超过 1e-6）会破坏精度。

3.2.1 样条权重函数的实现细节

def spline_weight(r, s): # r: 已归一化前的距离数组, s: 影响半径 xi = r / s w = np.zeros_like(xi) mask1 = xi <= 0.5 mask2 = (xi > 0.5) & (xi <= 1.0) w[mask1] = 2.0/3.0 - 4.0*xi[mask1]**2 + 4.0*xi[mask1]**3 w[mask2] = 4.0/3.0 - 4.0*xi[mask2] + 4.0*xi[mask2]**2 - (4.0/3.0)*xi[mask2]**3 return w

边界附近需要特殊处理：当评估点靠近计算域边界时，影响域被截断，力矩矩阵可能接近奇异。这是 EFG1 最经典的稳定性陷阱。解决方案是遇到奇异时自动扩充半径到原来 1.2 倍并重算。

3.3 刚度矩阵组装与线性求解

K = np.zeros((nnode*2, nnode*2)) # 2D 问题每个节点 2 个自由度 force = np.zeros(nnode*2) E, nu = 1e5, 0.3 Dmat = E / (1 - nu**2) * np.array([[1, nu, 0], [nu, 1, 0], [0, 0, (1-nu)/2]]) for (x1, y1, x2, y2) in bg_cells: # 2x2 高斯积分 gauss_pts, gauss_wts = gauss2d(x1, y1, x2, y2) for gp, gw in zip(gauss_pts, gauss_wts): pt = np.array(gp) phi, D, idx = mls_shape_and_derivatives(pt, nodes, support_radius=0.35) # 组装局部应变矩阵 B (每个点扫描全套 D) # B = [dphi/dx 0; 0 dphi/dy; dphi/dy dphi/dx] # 注意此处用克尔系数 D 矩阵中的导数结构 B = build_local_B(phi_derivs, nnode10, idx) Ke = B.T @ Dmat @ B * gw * (x2-x1) * (y2-y1) for a, ia in enumerate(idx): for b, ib in enumerate(idx): K[2*ia:2*ia+2, 2*ib:2*ib+2] += Ke[a*2:a*2+2, b*2:b*2+2] # 施加左端本质边界条件 by 罚函数法 penalty = E * 1e6 for nid, node_pos in enumerate(nodes): if node_pos[0] < 1e-10: for dof in range(2): K[2*nid+dof, 2*nid+dof] += penalty force[2*nid+dof] = 0.0 u = np.linalg.solve(K, force)

本质边界条件用罚函数法，罚系数取弹性模量的 (10^6) 倍组装进对角线。代码段中的support_radius=0.35是近似值，实际应按节点间距的倍率动态计算：设节点间距为 h，则半径取 2.8h~3.2h，均匀网格下 h=0.1，则半径约为 0.3。固定半径的问题在于非均匀网格下部分节点会孤立。

4. EFG1 无网格法的门槛参数与三个常见失败模式

4.1 影响半径的定量调整策略

影响半径 s 是 EFG1 里最重要的超参数，比 FEM 里网格密度还敏感。s 过小，力矩矩阵接近奇异，位移场出现棋盘式振荡；s 过大，形函数趋于恒定，应力集中被完全抹掉。推荐计算公式：

[ s_I = d_{I,\min} \times \alpha ]

其中 ( d_{I,\min} ) 是节点 I 到最近邻节点的距离，( \alpha ) 取值范围 2.5~4.0。二维均匀网格时 ( \alpha ) 建议 3.0，非均匀网格下要让每个节点影响域内至少覆盖 8 到 12 个邻居节点，同时避免影响域跨过几何边界。

怎么验证当前半径是否合适？看位移解的高频振荡幅度。具体操作：对比线性基和二次基计算的同一问题结果，如果两者相对误差超过 5%，优先怀疑半径偏小；如果解的梯度场出现规则性波纹，则半径偏大需要缩小。

4.2 背景积分点数与精度之间的权衡

背景网格高斯积分的点数决定了求解精度的上限。EFG1 中单纯增加高斯点不会无限制提升精度，因为每个高斯点本身也是独立评估点。我用标准悬臂梁做收敛性测试时的经验数据：

均匀背景下 4x4 高斯点性价比最优。继续加密的时候误差下降变缓，但计算时间线性增长。值得注意，背景网格也需要跟随节点分布变化：如果节点数量多但背景网格粗，误差主导项是积分误差；反过来背景网格过细则只会增加求解耗时。

4.3 非均匀节点分布下的矩阵病态问题的处理与排错

工程里非均匀节点不可避免：裂纹尖端要加密，远离应力集中处可以放稀。非均匀分布会导致两个问题：一是影响半径不一致导致刚度矩阵对角占比浮动大；二是力矩矩阵在稀疏区域的二阶矩偏差大。

# 排错第一步：计算每个节点的邻居数量，输出低于阈值的区域 python check_neighbors.py --input nodes.vtk --radius-ratio 3.0 --min-neighbors 8

没写这个脚本前，先在实现里加一段诊断代码，输出二维直方图显示每个节点影响域覆盖的节点数。如果发现局部不足，把 alpha 临时从 3.0 调到 4.0 再看覆盖数是否达到门槛。这比反复试算位移结果快得多。

4.4 无网格法与有限元耦合时的界面网格过渡做法

在实际项目里很少全程用 EFG1，更常见的做法是用 FEM 处理规则区域，用 EFG1 处理裂纹、大变形或材料破坏区域，两者之间需要界面耦合。常见方案有三种。

桥接节点法：界面处保留 FEM 节点，EFG1 形函数在界面附近强制还原为 FEM 形函数。实现简单但精度在耦合界面明显下降，因为 MLS 在边界处不再保持一致性。

界面拉格朗日乘子：在界面上额外引入 Lagrange 乘子场，分别约束 FEM 位移和 EFG1 位移在各个高斯点到同一值。精度高，但系统矩阵出现鞍点结构，需要专用求解器，复杂度提升明显。

杂交函数的精确覆盖法：用类似 PU（单位分解）的方式在界面区域把 FEM 和 MLS 形函数做加权组合。这是稳健性和实现简便性的平衡点，工程中我用得最多。

耦合问题的关键检查项是界面处的应力连续性，不连续性超过 2%~3% 就要检查界面节点是否与两端网格都保持了合适的影响域重叠。

5. 应力计算中的 EFG1 后处理陷阱与自适应加密应用

EFG1 后处理和 FEM 有本质差异：FEM 中节点应力是单元应力的平均，而 EFG1 中的位移解本身就用 MLS 近似，直接对该近似求导得到的是光滑连续的应变场。看似省事了，但也带来特效问题：如果直接对 MLS 形函数求空间导数来计算应力，那么在支持域边缘，由于权重函数截断，应力会出现微小波动。实际做法是对应力做第二重 MLS 平滑，即把参考点的应力值投射回节点上做加权平均，再在节点间插值。这个平滑本质上是一个低通滤波，对于应力集中区域要格外小心。

# MLS 应力平滑示例 node_stress = np.zeros(nnode) # 在所有节点上评估应力场，然后做移动最小二乘重投影 for i, nd in enumerate(nodes): diffs = nodes - nd dist2 = np.sum(diffs**2, axis=1) weight = np.exp(-dist2 / (support**2)) numerator = np.sum(weight * raw_stress) denom = np.sum(weight) node_stress[i] = numerator / denom

这个简单的高斯核平滑在均匀网格上效果不错，但核宽度要根据局部节点密度自适应。核宽过大，裂纹尖端的应力峰值被削弱；核宽过小，平滑作用不明显。做法是用每个节点到第 k 个最近邻的距离作为核的局部尺度。

工程里 EFG1 的主要价值在断裂力学：裂纹扩展不需要 remesh，因为节点一直就在那里，裂纹面两侧的连续性靠修改权重函数的可见性消掉。具体做法是在计算权重时做可见性检查，若评估点和节点连线穿过裂纹面，则该节点权重清零。这是 EFG1 相对 FEM 最自然的优势。

自适应加密上，常用指示因子是应变能密度的后验误差：

[ \eta_I = \sqrt{\sum_{J \in N_I} (| \varepsilon_J^{raw} - \varepsilon_J^{smooth} | \cdot \Omega_J)} ]

逐节点算出该值后，按从大到小排序，将不均匀度最大的 10%~20% 节点附近新增节点，重新离散后再求解。由于无网格法的形函数完全由节点位置决定，加节点后不需要任何连通性更新，代码层面只是往节点列表里 append 几个点，整个流程比 FEM 重划网格简单得多。

最后提一条实战经验：EFG1 的效率和精度很大程度上取决于局部支撑域搜索的实现。做三维瞬态问题时，k-d 树的构建开销能占到整个求解时间 30% 以上，可以把影响域搜索和积分循环融合。如果在弹性静力学里你的 EFG1 代码比同等自由度 FEM 慢超过 50 倍，先检查是不是在每个高斯点全量扫描了所有节点——用空间哈希或 k-d 树把这部分降下来，整个求解时间能直接缩短一个数量级。

本文还有配套的精品资源，点击获取