简介无网格法是一种无需预先生成网格的数值计算方法在处理自由边界、变形固体和非线性问题时比传统网格类方法更灵活。资源围绕 EFG1无网格伽辽金法的实现展开包含完整的 MATLAB 函数脚本与配套数据适合正在学习无网格法、需要参考代码来理解节点分布、插值构造与方程离散流程的学生或研究人员。压缩包共 4 个文件其中 2 个 .m 脚本负责核心计算与节点定义1 个 .asv 为 MATLAB 自动备份文件另 1 个 .mat 用于存储边界或计算结果整体仅 1KB轻量但结构清晰。已有 266 人浏览学习。通过对照代码和数据可以快速掌握 EFG 法的程序骨架包括插值函数选择、系数矩阵组装的编码思路并对后续扩展至 SPH、DEM 等其他无网格方法提供基础。1. EFG1 无网格法从 FEM 的网格枷锁里挣脱的第一代方案EFG1Element 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 为该影响半径。三次样条的优点是导数在 r0 和 r1 处连续但二阶导在 r1/2 处有跳变。四次样条则额外保证了二阶导连续收敛性更好代价是更多的浮点计算。我一般优先用四次样条因为 EFG1 在计算刚度矩阵时要求形函数导数导数高阶连续性直接影响应力场的平滑度。权重函数影响半径 s 的取值原则二维问题中每个节点的影响半径通常取该节点到最近邻节点距离的 2.5 到 4.0 倍。太小则覆盖的邻居节点不够、力矩矩阵 A 奇异太大则拟合过于光滑局部特征被磨平。2.2 Galerkin 离散与刚度矩阵装配流程EFG1 的离散过程和 FEM 同构把试探函数和检验函数都投影到 MLS 形函数张成的空间里。以二维弹性静力学为例位移场 ( u(x) ) 近似为[ u(x) \sum_{I1}^{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 形函数不满足插值性边界节点的近似值是该节点邻域内的拟合值并非精确节点值。直接施加本质边界条件会产生系统性误差误差大小和权重函数的边界截断有关。常见处理方案方法实现复杂度精度稳定性拉格朗日乘子法高高需小心求解器罚函数法低中罚系数敏感修正变分原理中高稳定Nitsche 法中高稳定工程中时间受限时用罚函数法起步罚系数取材料弹性模量的 (10^4\sim10^6) 倍。研究精度优先时上拉格朗日乘子或者 Nitsche 法。罚函数法的致命弱点是病态条件数罚系数过大刚度矩阵条件数从 (10^8) 量级飙升到 (10^{14}) 以上线性求解器直接失准。正确做法是用双精度求解并用迭代法时做残差监控。3. 用 Python 实现一个最小可运行的 EFG1 悬臂梁求解器3.1 节点离散与背景网格生成求解对象选经典悬臂梁长 L2.0 m高 H1.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], dtypefloat) 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[i1] yb, yt y_vals[j], y_vals[j1] # 每个背景单元细分 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, basisquadratic): # 稀疏局部列表不需要全量搜索 diff nodes - pt dist np.linalg.norm(diff, axis1) 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_radius0.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*ia2, 2*ib:2*ib2] Ke[a*2:a*22, b*2:b*22] # 施加左端本质边界条件 by 罚函数法 penalty E * 1e6 for nid, node_pos in enumerate(nodes): if node_pos[0] 1e-10: for dof in range(2): K[2*niddof, 2*niddof] penalty force[2*niddof] 0.0 u np.linalg.solve(K, force)本质边界条件用罚函数法罚系数取弹性模量的 (10^6) 倍组装进对角线。代码段中的support_radius0.35是近似值实际应按节点间距的倍率动态计算设节点间距为 h则半径取 2.8h~3.2h均匀网格下 h0.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 中单纯增加高斯点不会无限制提升精度因为每个高斯点本身也是独立评估点。我用标准悬臂梁做收敛性测试时的经验数据背景单元内高斯点数位移相对误差计算时间比2x26.2%1.0x4x40.8%3.4x6x60.3%7.1x均匀背景下 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, axis1) 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 树把这部分降下来整个求解时间能直接缩短一个数量级。本文还有配套的精品资源点击获取 EFG1无网格法：原理、实现与工程实践
简介无网格法是一种无需预先生成网格的数值计算方法在处理自由边界、变形固体和非线性问题时比传统网格类方法更灵活。资源围绕 EFG1无网格伽辽金法的实现展开包含完整的 MATLAB 函数脚本与配套数据适合正在学习无网格法、需要参考代码来理解节点分布、插值构造与方程离散流程的学生或研究人员。压缩包共 4 个文件其中 2 个 .m 脚本负责核心计算与节点定义1 个 .asv 为 MATLAB 自动备份文件另 1 个 .mat 用于存储边界或计算结果整体仅 1KB轻量但结构清晰。已有 266 人浏览学习。通过对照代码和数据可以快速掌握 EFG 法的程序骨架包括插值函数选择、系数矩阵组装的编码思路并对后续扩展至 SPH、DEM 等其他无网格方法提供基础。1. EFG1 无网格法从 FEM 的网格枷锁里挣脱的第一代方案EFG1Element 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 为该影响半径。三次样条的优点是导数在 r0 和 r1 处连续但二阶导在 r1/2 处有跳变。四次样条则额外保证了二阶导连续收敛性更好代价是更多的浮点计算。我一般优先用四次样条因为 EFG1 在计算刚度矩阵时要求形函数导数导数高阶连续性直接影响应力场的平滑度。权重函数影响半径 s 的取值原则二维问题中每个节点的影响半径通常取该节点到最近邻节点距离的 2.5 到 4.0 倍。太小则覆盖的邻居节点不够、力矩矩阵 A 奇异太大则拟合过于光滑局部特征被磨平。2.2 Galerkin 离散与刚度矩阵装配流程EFG1 的离散过程和 FEM 同构把试探函数和检验函数都投影到 MLS 形函数张成的空间里。以二维弹性静力学为例位移场 ( u(x) ) 近似为[ u(x) \sum_{I1}^{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 形函数不满足插值性边界节点的近似值是该节点邻域内的拟合值并非精确节点值。直接施加本质边界条件会产生系统性误差误差大小和权重函数的边界截断有关。常见处理方案方法实现复杂度精度稳定性拉格朗日乘子法高高需小心求解器罚函数法低中罚系数敏感修正变分原理中高稳定Nitsche 法中高稳定工程中时间受限时用罚函数法起步罚系数取材料弹性模量的 (10^4\sim10^6) 倍。研究精度优先时上拉格朗日乘子或者 Nitsche 法。罚函数法的致命弱点是病态条件数罚系数过大刚度矩阵条件数从 (10^8) 量级飙升到 (10^{14}) 以上线性求解器直接失准。正确做法是用双精度求解并用迭代法时做残差监控。3. 用 Python 实现一个最小可运行的 EFG1 悬臂梁求解器3.1 节点离散与背景网格生成求解对象选经典悬臂梁长 L2.0 m高 H1.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], dtypefloat) 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[i1] yb, yt y_vals[j], y_vals[j1] # 每个背景单元细分 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, basisquadratic): # 稀疏局部列表不需要全量搜索 diff nodes - pt dist np.linalg.norm(diff, axis1) 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_radius0.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*ia2, 2*ib:2*ib2] Ke[a*2:a*22, b*2:b*22] # 施加左端本质边界条件 by 罚函数法 penalty E * 1e6 for nid, node_pos in enumerate(nodes): if node_pos[0] 1e-10: for dof in range(2): K[2*niddof, 2*niddof] penalty force[2*niddof] 0.0 u np.linalg.solve(K, force)本质边界条件用罚函数法罚系数取弹性模量的 (10^6) 倍组装进对角线。代码段中的support_radius0.35是近似值实际应按节点间距的倍率动态计算设节点间距为 h则半径取 2.8h~3.2h均匀网格下 h0.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 中单纯增加高斯点不会无限制提升精度因为每个高斯点本身也是独立评估点。我用标准悬臂梁做收敛性测试时的经验数据背景单元内高斯点数位移相对误差计算时间比2x26.2%1.0x4x40.8%3.4x6x60.3%7.1x均匀背景下 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, axis1) 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 树把这部分降下来整个求解时间能直接缩短一个数量级。本文还有配套的精品资源点击获取