简介无网格法是一种无需预先生成网格的数值计算方法尤其适合自由边界、大变形及高度非线性问题的求解。压缩包聚焦无网格法的基础实现适合正在学习数值方法、需要相关示例代码的本科生或工程师参考使用。压缩包共4个文件包括2个MATLAB脚本、1个MATLAB矩阵数据文件以及1个自动保存文件整体体积仅1KB。脚本主要对应节点分布处理与插值权重构造数据文件可能存储边界条件或计算中间结果共同构成一个极简的无网格法求解流程示例。目前已有266人学习下载。资源体量虽小却覆盖了从节点选取、插值函数构造到离散控制方程求解的完整代码线索便于读者对照移动最小二乘法、径向基函数等经典插值方法理解实现细节也可以作为快速验证无网格思路的轻量参考甚至为后续学习光滑粒子流体动力学等扩展方法打下基础。1. 无网格法不是黑匣子EFG1压缩包到底装了什么无网格法这三个字听起来像是离工程很远但一套几百行的 MATLAB 代码就能把核心流程跑通。EFG1 就是这种适合练手的资源解压后四个文件其中两个是 .m 脚本、一个是权函数、一个是 .mat 数据文件刚好覆盖节点分布、形函数构造、离散求解和结果保存这条完整链路。它解决的问题是有限元遇到网格畸变、裂纹扩展、大变形时卡壳的那一类场景。适合正在读研做数值模拟、或者从有限元转无网格法但不想先啃一堆公式的同学。把这份资源完整跑一遍你对 EFG 的认知会比看十篇综述更牢。2. 无网格法与 EFG 方法为什么不用网格反而更麻烦很多人第一次听说无网格法时以为它是把有限元的网格删掉就完事。真上手会发现节点怎么撒、支持域怎么定、积分怎么做比画网格更微妙。EFG1 这套资源选的是 EFGElement-Free Galerkin无单元 Galerkin方法它在无网格法家族里属于“精度高但需要背景积分”的那一派。这一章先把它的理论位置和核心步骤讲清楚后面几章再落到文件和参数上。2.1 从有限元到无网格问题出在网格身上用过有限元的人都知道模型里最脆弱的部分往往是网格本身。二维问题还好三维问题一旦遇到大变形单元扭曲到一定程度雅可比行列式变号求解直接发散。更常见的是裂纹扩展每往前一步裂纹尖端的网格都要重新划分一次前一步的解还要映射到新网格上一不留神误差就把结果吃掉。做动态裂纹模拟的人经常把一半时间花在网格重划分和场量映射上而不是物理问题本身。无网格法换了个思路只要在计算域里撒一批节点每个节点自带坐标和自由度然后通过局部拟合构造形函数不需要单元连通关系。节点可以随变形自由移动也可以重新分布裂纹尖端加节点、删除节点都很方便。EFG1 文件名里的 del.m 和 del.asv我倾向于理解成与“删除/更新节点”或“分布”有关的前缀这也正好对得上无网格法里动态调整节点分布的需求。但“无需网格”并不意味着“无需积分”。EFG 这类基于弱形式的方法最后还是要把偏微分方程在计算域上积分于是引入背景积分单元或背景网格。这个背景网格只承担积分功能不承担插值功能节点之间依然不需要连通。换句话说EFG 是“无网格插值有网格积分”这也是它区别于 SPH 的一个关键点。2.2 EFG 的核心移动最小二乘与背景积分EFG 的形函数由移动最小二乘Moving Least Squares, MLS近似构造。给一个评估点 x在它的影响域内取若干节点 x_i用一组基函数做加权最小二乘拟合让近似函数在所有采样点上误差最小。MLS 的形函数写成φ_i(x) p^T(x) A^{-1}(x) B_i(x)其中 A(x) Σ w_i(x) p(x_i) p^T(x_i)B_i(x) w_i(x) p(x_i)p 是多项式基向量w_i 是权函数。这里最关键的就是权函数它决定了每个节点对评估点的影响强度和作用范围。EFG1 里的 cubwgt.m按照命名习惯大概率是三次样条权函数cubic spline weight“cub”指三次。这个函数写出来常见是这样function w cubwgt(r, dm) % 三次样条权函数EFG 中最常用的紧支撑权函数 % r - 评估点到节点的距离 % dm - 影响半径 if r 0 w 1; elseif r 0.5 * dm w 2/3 - 4*(r/dm)^2 4*(r/dm)^3; elseif r dm w 4/3 - 4*(r/dm) 4*(r/dm)^2 - 4/3*(r/dm)^3; else w 0; end end这段代码把权函数分成三段距离为 0 时权重取 1保证节点自身信息有最大权重中间段用三次多项式过渡超过影响半径后权重直接归零实现紧支撑。参数 dm 本质上是每个节点的“势力范围”选小了会导致局部矩阵奇异选大了会让形函数过于平滑反而抹掉局部梯度。后面第四章会专门讲怎么调。有了形函数EFG 的求解流程和有限元很像先建立弱形式把控制方程乘上试探函数再积分然后代入 MLS 形函数得到关于节点位移的线性方程组 K u F。区别在于有限元的 K 是单元刚度矩阵装配出来的EFG 的 K 是每个背景积分点对全局矩阵直接累加矩阵通常是满阵没有有限元那种窄带稀疏结构。2.3 从弱形式到代数方程组EFG1 里 K u F 怎么组装以二维线弹性问题为例平衡方程和边界条件都喜欢省略但离散过程是绕不开的。对每个背景积分点 x_q先计算该点影响域内所有节点的形函数向量 φ 和导数矩阵 dφ然后构造应变-位移矩阵 B最后累加K B^T D B * det(J) * w_q这里的 D 是材料本构矩阵w_q 是高斯积分权重det(J) 是背景积分单元的雅可比行列式。EFG1 里的 del.m 如果按照常见教学代码结构来写主体就是对这个循环的封装。我一般拿到这类代码会先看它的积分循环长什么样而不是急着运行% del.m 的典型循环逻辑示意实际以源码为准 for e 1:size(elems,1) [gp, wq] gauss_rule(elems(e,:), nodes); % 背景积分单元上的高斯点 for q 1:length(wq) [phi, dphi] mls_shape(gp(q,:), nodes, dm); % 调 cubwgt 和基函数 B build_B(dphi(1), dphi(2)); K K B * D * B * wq(q) * detJ(e); end end u Kff \ ff;这段循环说明了两个容易被新手忽略的点。第一背景积分单元和节点是两套独立的结构积分点周围如果没有足够节点B 矩阵算出来是零K 就会奇异。第二K 是全局稠密矩阵规模不大时直接用反斜杠求解不需要像有限元那样做带宽优化。EFG1 里如果自由度在几百到几千这个范围高斯消元法就是最稳的选择。3. 拆解 EFG1 压缩包四个文件在解决什么下载资源后别急着双击运行。先解压到独立文件夹用记事本或 MATLAB 编辑器打开每个文件看一眼。多花五分钟能省掉后面两小时的报错排查。这一章把四个文件的角色和运行路径讲明白。3.1 文件清单del.asv、cubwgt.m、del.m、bbb.mat 各管哪一段文件名类型在我看来的作用运行/使用优先级del.mMATLAB 主脚本节点分布、离散、组装、求解的主流程直接运行cubwgt.mMATLAB 函数三次样条权函数供形函数计算调用被 del.m 调用del.asvMATLAB 自动保存文件del.m 的自动备份可能比 del.m 旧不要作为入口bbb.matMAT 数据文件节点坐标、边界条件或背景积分信息由 del.m 读取这里最坑的就是 del.asv。MATLAB 编辑器默认会为正在编辑的 .m 文件生成一个 .asv 自动存档文件文件名相同、后缀不同。原作者打包时如果没清理这个文件就会混进来。它的内容往往是不完整的中间版本直接当脚本运行大概率报错或者结果不对。正确的做法是忽略它或者用movefile del.asv del_old.m把它备份成另一个名字再对比两个版本之间有没有有效内容。bbb.mat 是另一个容易迷糊的点。很多人一上来就load(bbb.mat)然后发现工作区没有预期变量。正确的检查方式是在 MATLAB 命令窗口执行whos -file bbb.mat先看里面存了哪些变量名、维度多少。常见情况下它可能存着 nodes、elems、dbc 之类的变量但作者也可能用 bbb 这个结构体把全部信息包起来调用时要写bbb.nodes而不是nodes。3.2 用 MATLAB 跑通 EFG1从解压到看到结果的完整步骤我建议在 MATLAB 里用绝对路径进入目录避免当前目录被切换导致找不到函数。下面几步是比较保险的起手式cd(D:\work\EFG1); % 换成你的解压路径别带中文和空格 ls % 列出文件确认四个文件都在 whos -file bbb.mat % 查看数据文件里的变量先摸底 run(del.m); % 执行主脚本执行完whos -file bbb.mat后你会看到变量名和尺寸。比如如果看到nodes是 100×2 的矩阵说明计算域里有 100 个节点每个节点存 x、y 坐标如果看到elems是 80×4 的矩阵说明背景积分单元用的是四边形。这些信息有助于你判断 del.m 内部是否对数据结构有硬编码假设。如果run(del.m)报错先看错误信息里提到哪个文件名九成是路径或变量名问题。跑通之后预期输出会因为具体算例不同而不同。常见的是在命令窗口打印节点数、方程自由度数、最大位移有的版本会直接画出一张位移云图。如果运行后工作区里出现了 u 或 U 变量那说明求解完成。没有输出图表也没关系用whos看一下工作区变量只要 K 和 u 存在后面怎么可视化都可以自己补。3.3 数据流bbb.mat 怎么和 del.m 配合很多教学代码喜欢把预处理结果直接存成 .mat 文件省得每次运行都重新生成节点。EFG1 里的 bbb.mat 大概率就是这个角色。del.m 在运行开始时用load(bbb.mat)把节点和边界信息读进来之后每个背景积分点上的 MLS 近似都会用到这些节点坐标。如果你要改成自己的算例有两种做法。第一种是直接改 bbb.mat在 MATLAB 里重新定义 nodes、边界条件然后用save(bbb.mat,nodes,elems,dbc)覆盖原文件再运行 del.m。第二种是改 del.m把它开头的 load 部分注释掉改成用meshgrid或linspace直接在脚本里生成节点。我更喜欢第二种因为每次运行前都清楚当前节点是什么样不会被旧的 .mat 文件误导。4. 关键参数对照权函数、影响半径与求解器选择无网格法跑出结果不难跑出“正确”的结果难。十个无网格代码里有八个能算出曲线但曲线对不对取决于三个参数权函数形式、影响半径、求解方式。EFG1 里的 cubwgt.m 已经定死了权函数剩下两件事需要你根据实际问题手动调。4.1 权函数cubwgt.m的“立方”到底指什么cubwgt.m 名字里的 cub我判断是 cubic 的意思也就是三次样条权函数。这类权函数的特点是二阶连续MLS 形函数的光滑性也能跟上去。在 MATLAB 里调用它时注意输入参数的单位要保持一致。如果你把节点坐标换算成毫米影响半径 dm 也要用毫米否则权函数会在整个域里都失效。一个容易犯的错误是把权函数当成形函数本身。cubwgt 只是 MLS 近似里的权重最终形函数是基函数、A 矩阵和 B 矩阵的组合。你在 del.m 里会看到类似w cubwgt(r, dm)的调用随后紧接着A A w * p * p.这后面的累加才是 ML S 的核心。改 cubwgt.m 里的三次多项式系数会改变形函数的光滑性但不会直接改变物理方程所以调试时先别动它。4.2 影响半径与节点密度精度和稳定性怎么平衡影响半径 dm 与节点平均间距 h 的比值我习惯记成 rho dm / h。rho 太小每个评估点周围节点不够A 矩阵接近奇异位移场会出现不自然的振荡rho 太大权函数覆盖太多节点MLS 拟合过于平均局部大梯度被抹平。对线弹性问题rho 取 1.52.5 是比较稳的区间。rho dm / h常见现象建议 1.2矩阵奇异、空域存在脚本报错增大 dm 或加密节点1.52.5收敛稳定精度良好保持 3.5精度下降计算量明显增加减小 dm或增大节点间距让节点间距均匀很重要。有些算例里为了模拟局部突变会在局部加聚节点这时候 dm 不能全模型用同一个值。我一般会在 del.m 里对每个节点单独计算它周围另一个节点的平均距离再乘一个系数得到该节点的影响半径。代码可以写成% 按节点间距动态设置影响半径 for i 1:size(nodes,1) dist sqrt(sum((nodes - nodes(i,:)).^2, 2)); dist(dist 0) inf; % 去掉自身 h_i min(dist); % 最近邻距离 dm(i) 1.8 * h_i; % 常见系数取 1.8 end这段代码的逻辑是先找每个节点的最近邻距离再乘以 1.8 作为影响半径。动态影响半径的好处是能适应非均匀节点分布坏处是计算量比固定值大一些。对 EFG1 几百个节点的规模完全不用在意这点开销。4.3 线性方程组求解为什么高斯消元在小规模 EFG 里够用EFG 的刚度矩阵是满阵因为每个背景积分点上的形函数可能覆盖很多节点两个相距甚远的自由度之间也有非零耦合。有限元里的稀疏带状优势在 EFG 里不存在。所以求解时MATLAB 的A\b是最省事的方案它内部会基于列主元消去法或自适应的 LDL 分解对几千阶的满阵也能秒级求解。如果问题规模超过几万自由再想用A\b就有点吃力了。这时候我一般不会在 EFG1 这种教学资源里做大规模扩展而是转去考虑预处理共轭梯度法或者做区域分解。但工程手算和课程作业阶段直接法完全够。你在 del.m 里看到u Kff \ ff;那行别手痒把它改成inv(Kff) * ff反斜杠更快也更稳inv不仅慢还容易让人忽略奇异矩阵警告。5. 避坑手册五个让 EFG1 翻车的常见问题无网格法的调错和有限元不太一样有限元报错通常能定位到某个单元无网格法报错往往是全局性的。这里写五个我在拿这类资源练手时真遇到过的问题每一条都是现象、原因、解决三步走。5.1 现象矩阵奇异或“Matrix is close to singular”警告运行 del.m 到求解步骤时MATLAB 弹出警告说矩阵接近奇异或者干脆算出几百亿的位移。原因通常有两个影响半径 dm 太小导致某些背景积分点周围没有足够节点或者节点坐标里有重复点让矩阵出现线性相关。解决方法是先检查节点生成部分用unique去掉重复坐标再把 dm 增大到最近邻距离的 1.5 倍以上。我一般会临时在代码里输出每个积分点影响域内的节点数量如果发现小于 3 个立刻就能定位问题区域。5.2 现象边界处位移严重偏离理论解如果你已经能跑出曲线但发现边界处数值解跟解析解差很多先别怀疑位移边界条件是写错。EFG 的 MLS 形函数不满足 Kronecker delta 性质也就是说边界节点上的近似位移并不等于节点自由度值。直接像有限元那样把边界自由度强行改成给定值边界附近必然出误差。解决方法是改用罚函数法或者拉格朗日乘子法施加本质边界条件。EFG1 如果没处理这一步你需要自己在组装 K 之后处理边界节点。罚函数法里取一个很大的惩罚系数 alpha通常比材料模量大 6 个量级边界条件就能以近似的形式进入方程。5.3 现象cubwgt.m 报错说维度不匹配错误信息可能类似“Matrix dimensions must agree”。原因是调用 cubwgt 时传入的 r 和 dm 维度不一致或者 r 是一个向量而代码里写了r 0这种标量判断。解决方法是先把 r 统一成列向量再把 dm 展开成和 r 相同长度的向量。我习惯在 cubwgt 的第一行加一句r r(:);强制转成列向量避免外界调用时传出行向量。这个改动不影响原逻辑但能让脚本更抗造。5.4 现象把 del.asv 当成主脚本运行结果就是个错误练习MATLAB 编辑器在保存 del.m 时自动生成了 del.asv时间一长你会忘了哪个是最新版。用type del.asv看一眼经常发现它只有前半段代码连求解和画图都没写。解决方法是养成看修改时间的习惯在文件浏览器里把 del.m 和 del.asv 的“修改日期”列出来用最新的 .m。如果你确实需要在 del.asv 里找回一段误删的代码把它复制出来贴在 del.m 里不要直接运行 .asv 文件。5.5 现象结果震荡或发散但矩阵没有奇异如果你把影响半径调得很大比如 rho 超过 3.5矩阵不奇异但结果开始出现锯齿状振荡。原因在于权函数覆盖范围太大MLS 拟合把局部的曲率变化过度平滑同时背景积分误差也被放大。解决方法是减小 dm并且检查背景积分单元的高斯点数。EFG 对积分点的要求比有限元严格很多实现里每个背景单元默认只有一个积分点这对二次变化的问题是不够的。常见做法是用 2×2 高斯积分即四个积分点再乱就换 3×3。6. 从 EFG1 延伸用一维杆问题校准你的无网格代码跑通 EFG1 只是第一步。真正要把无网格法用到自己的问题上你需要一个能快速判断“代码没写错”的方法。我最常用的是拿一维杆问题做收敛性验证因为一维问题有精确解析解节点少计算快任何形函数或积分上的问题都能暴露出来。6.1 一维杆问题的解析解与 EFG 收敛性代码取一根长度 L 1 的杆左端固定右端施加均布拉力 F 1材料弹性模量 E 1解析解为 u(x) F x / (E)。这一步可以直接用 EFG1 里的思路只保留 x 方向把二维形函数退化成一维。下面这段校核代码可以放到 del.m 之后运行L 1; E 1; F 1; exact (x) F * x / E; for n [5 9 17 33] x linspace(0, L, n); u solve_efg_1d(x, E, F); % 自行实现的一维 EFG 求解函数 err max(abs(u - exact(x))); fprintf(n%3d, Loo_err%1.2e\n, n, err); end这段代码的用途是观察误差随节点数的收敛情况。如果从 n5 到 n33误差在持续下降说明你的无网格实现是自洽的如果误差不降或突然变大问题多半在影响半径或边界条件处理上。一维问题里影响半径取节点间距的 1.5 倍左右线性基函数就可以了。6.2 从 bbb.mat 里挖出你需要的结果并保存很多算例跑完后结果只存在工作区一关 MATLAB 就没了。我会建议在 del.m 结尾追加一段保存代码把节点、位移和网格统一存成一个结构化数据文件方便后续画图或者写论文时重新调用% 保存结果 result.nodes nodes; result.u u; result.parameters struct(dm, dm, nGP, nGP); save(result.mat, -struct, result); load(result.mat); scatter(result.nodes(:,1), result.nodes(:,2), 20, result.u(:,1), filled); colorbar; title(Displacement u_x);这段代码的关键是save(result.mat, -struct, result)它把 result 结构体里的每个字段变成 .mat 文件中的独立变量下次读取时可以直接load出 nodes 和 u。而后面的 scatter 则画出一张节点位移云图。只要位移数据不是零矩阵云图颜色就应该是连续过渡的如果出现明显的斑马纹说明某些节点的位移异常回头找影响半径的问题。从那以后我每次拿到一份无网格法源码都先不做复杂算例而是跑一遍一维杆或者悬臂梁的解析解收敛性验证。这个习惯帮我避开过不少看着正常、实际精度被边界条件吃掉的翻车案例。记住无网格法的参数不是玄学但需要你亲自把每个量测一遍。希望帮到你。本文还有配套的精品资源点击获取