简介对于计算力学与数值仿真中界面演化问题的研究者而言这套将Level Set方法与有限元法FEM相结合的MATLAB程序包提供了从理论到代码的完整参考适合偏微分方程基础扎实的学生、科研人员及工程师学习与复用。资源共32个文件以27个.m脚本为主体辅以PDF说明文档、README、版本控制配置文件压缩包仅约1MB紧凑而便于快速下载和本地复现。代码实现覆盖Level Set函数初始化与重新初始化、Eikonal方程数值求解、有限元刚度/质量/载荷矩阵装配、边界条件处理、界面耦合以及区域误差计算等关键环节并将水平集对流更新、保守型和标准型重新初始化等高级算法封装为独立模块。借助Level Set函数的隐式表达程序在界面移动时无需重新剖分网格可自然处理拓扑合并与分裂显著简化了传统界面追踪实现。资源还配有演示算例、基准测试和仿真动画脚本可直接运行观察界面运动与拓扑变化。目前已有156人学习过该资源对于想掌握LSM与FEM耦合实现、开展流体力学、固体力学或相变问题数值模拟的读者是一份难得的实践参考。1. 当水平集方法撞上有限元MATLAB 里的一类移动边界求解器第一次见到dekespo-Level_Set_Method_FEM_MATLAB.zip这种命名我猜不少人和我一样先想到的是研究生课程作业或者论文复现包。名字里三个词其实已经交代了技术路线用 level set 隐式描述移动界面用有限元法在背景网格上求解物理场最后用 MATLAB 把这两套系统时序耦合在一起。这类代码通常解决的是相变传热、拓扑优化、液滴运动、裂纹扩展等既需要跟踪界面又需要求解场变量的问题。适合正在写数值方法课作业、做科研复现或者想把水平集-FEM 耦合框架快速跑通再做二次开发的工程师。关键不在代码本身而在理解两个方程之间的数据传递节奏。2. 先把理论摆正水平集方程、Galerkin 离散与两套变量的耦合顺序2.1 水平集函数如何描述移动界面水平集方法的核心是用一个高维函数 φ(x,y,t) 的零等值面隐式表示界面。所谓“高维”对于二维问题就是“在平面内定义函数用 φ0 这条线代表界面”。φ 通常初始化为符号距离函数界面内为负、界面外为正函数值等于到界面的距离。这种表示最大的好处是拓扑变化比如两个液滴合并或一个界面分裂不需要刻意处理节点连接关系因为界面是被“切”出来的。界面运动满足哈密顿-雅可比型方程∂φ/∂t F|∇φ| 0其中 F 是界面沿外法向的速度。这个方程是一个双曲型守恒律空间离散需要用迎风格式而不是 FEM 那种中心差分。熟悉计算流体的人都明白直接对 |∇φ| 做中心差分会引起振荡导致界面扭曲。这一点是后面写 MATLAB 时最容易踩的坑。F 的来源正是有限元解。在纯几何问题里F 可以是常量比如圆扩大但耦合问题中F 依赖温度场、应力场或浓度场。这个依赖关系把两个物理系统绑在一起。2.2 有限元那边的控制方程以移动边界热传导为例假设我们要模拟一个固体-液体界面在温度场中移动的过程。整个计算域内热导率在界面两侧不同控制方程为稳态热传导∇·(k(x,y) ∇T) 0这里的 k 取决于坐标而坐标所在的区域由 φ 的正负决定。标准做法是用平滑 Heaviside 函数 H_ε(φ) 把两种材料参数连续地“混合”起来k(φ) k_solid (k_liquid - k_solid) * H_ε(φ)其中 H_ε 是 φ 的连续近似常用 tanh 函数ε 控制过渡带的宽度。将上式代入弱形式得到∫Ω k(φ) ∇T · ∇v dΩ 0对 v 取标准线性有限元形函数就能得到刚度矩阵 K求解 K T 0加上边界条件后。这个矩阵因 φ 的变化而随时间更新所以每推进一个水平集时间步就要重新组装一次 K。注意这里用“稳态”而非瞬态纯属省事。界面移动速度远小于热扩散速度时这种准静态假设是合理的也是这类代码最常用的简化方式。如果想做瞬态还需要质量矩阵和时间积分格式框架一样只是多一个矩阵。2.3 两个系统的数据交换谁传给谁传的是什么耦合的顺序要非常明确用当前 φ 更新材料参数 k。组装有限元矩阵求解温度场 T。由 T 计算界面法向速度 F。把 F 带入水平集方程推进 φ 一个时间步。每隔若干步重新初始化 φ让它恢复符号距离函数性质。其中第 3 步是耦合的关键也是最容易被跳过或写错的一步。F 在节点上的值一般不自直接相邻于温度场比如在 Stefan 问题中界面速度取决于界面两侧的温度梯度突跳。在 MATLAB 实现里常见做法是在界面附近提取温度梯度再投影到一周的节点上。这类代码里通常有一个get_velocity_from_temperature函数输入 T 和 φ输出 F。为了让你直观看到材料参数的耦合下面这段 MATLAB 代码演示如何从 φ 生成热导率分布function kappa phi_to_kappa(phi, k_solid, k_liquid, epsilon) % phi: N-by-1每个节点上的水平集函数值 % epsilon: 界面过渡带半宽一般取 1.5~3 倍网格尺寸 H 0.5 * (1 tanh(phi / epsilon)); kappa k_solid (k_liquid - k_solid) * H; end这里的tanh平滑函数让热导率在界面两侧连续过渡避免有限元组装时因材料参数跳变而出现局部振荡。epsilon 如果取得比网格尺寸还小过渡带内会出现一个单元的突变导致刚度矩阵条件数变差取得太大则界面被抹得模糊影响 F 的计算精度。这是第一个要调的参数。下面用一个表格总结两套系统的分工子系统方程类型离散方法输出给另一半的变量水平集双曲型哈密顿-雅可比有限差分迎风材料参数 k(φ)几何位置 φ0有限元椭圆型热传导线性三角形 Galerkin节点法向速度 F时间耦合—显式 Eulerφ 推进 dt 后的新值这个框架不是唯一选择。有人会在水平集更新时重新生成贴合界面的有限元网格有人会使用边界面上的积分来光滑 F。但最稳妥、也最适合 MATLAB 快速复现的还是同一套背景网格、同一组节点用有限差分处理水平集方程用有限元处理物理方程。这种单网格做法虽然不炫但它把耦合调试成本降到最低。3. 在 MATLAB 里搭一个能运行的最小实现从网格生成到主循环3.1 矩形域三角网格与节点编号MATLAB 里做二维有限元最常用的基础网格是矩形域上的均匀三角剖分。每个小矩形分成两个三角形节点编号按列优先排列。这种网格的好处是单元面积相同刚度矩阵元素规律性强方便验证坏处是模板死板不适用于复杂几何。不过对水平集-FEM 演示算例来说已经足够。网格生成代码可以这样写nx 40; ny 40; % 每个方向的单元数 x linspace(0, 1, nx 1); y linspace(0, 1, ny 1); [X, Y] meshgrid(x, y); nodes [X(:), Y(:)]; % 第 i 个节点的坐标 % 节点编号与 (j,i) 网格点对应j 是 y 方向索引 id reshape(1:((nx1)*(ny1)), ny1, nx1); elements zeros(2 * nx * ny, 3); k 1; for j 1:ny for i 1:nx n1 id(j,i); n2 id(j,i1); n3 id(j1,i); n4 id(j1,i1); elements(k, :) [n1 n2 n3]; elements(k1, :) [n2 n4 n3]; k k 2; end end mesh.nodes nodes; mesh.elements elements;这段代码把每个矩形切成左下和右上两个三角形。这里节点序号的确定非常关键因为后面的有限元装配要不断用elements索引节点坐标如果节点顺序乱了三角形面积会出现负值。调试时可以在打开第一个单元的三条边看看是否零面积。3.2 线性三角形单元的刚度矩阵组装组装刚度矩阵是有限元中最核心的部分。在 MATLAB 中用稀疏矩阵存储可以避免全矩阵的内存爆炸。下面这个函数针对每个单元计算局部刚度矩阵再累加到全局矩阵。function K assemble_stiffness(mesh, kappa) N size(mesh.nodes, 1); Ne size(mesh.elements, 1); K sparse(N, N); for e 1:Ne idx mesh.elements(e, :); xe mesh.nodes(idx, 1); ye mesh.nodes(idx, 2); % 三角形面积负值说明节点顺序反了 area 0.5 * (xe(1)*(ye(2)-ye(3)) ... xe(2)*(ye(3)-ye(1)) ... xe(3)*(ye(1)-ye(2))); if area 0 error(单元 %d 的面积非正请检查节点顺序, e); end % 形函数梯度矩阵每列对应一个节点的 dN/dx, dN/dy grad zeros(2, 3); grad(:,1) [ye(2)-ye(3); xe(3)-xe(2)] / (2*area); grad(:,2) [ye(3)-ye(1); xe(1)-xe(3)] / (2*area); grad(:,3) [ye(1)-ye(2); xe(2)-xe(1)] / (2*area); % 单元热导率取三个节点平均值 ke mean(kappa(idx)); Ke ke * grad * grad * area; K(idx, idx) K(idx, idx) Ke; end end这个函数把三角形单元视为常系数区域ke取节点平均。如果你想要更高精度可以用高斯积分但在水平集耦合算例中因为材料参数本身有过渡带平均值的误差是可以接受的。面积area的正负判断是一种简单自查万一你手工录入节点顺序错了程序会直接报错而不是给出一个静默错误的结果。边界条件不在这个函数里处理。常见做法是在主循环中先组装完整矩阵然后找出边界节点用“置换行法”或“罚函数法”强制温度值。这里推荐罚函数法简单且不破坏稀疏结构function [K, F] apply_dirichlet(K, F, fixed_nodes, fixed_values) penalty 1e10 * max(max(abs(K))); % 足够大的数 K(fixed_nodes, :) 0; K(:, fixed_nodes) 0; K(fixed_nodes, fixed_nodes) penalty; F(fixed_nodes) penalty * fixed_values; end可以看到强制固定温度就是把对应节点上的对角元变成一个很大的数右端项设为penalty * fixed_values。这相当于用一个刚度很大的“弹簧”把该节点拉到目标值比直接删行保留列要更容易写进循环。3.3 水平集演化中的迎风差分与时间步进水平集方程是双曲型必须用迎风格式。一个稳定的选择是 Godunov 格式只对法向速度分量做单边差分。这里为了简洁用一阶差分并限制 F 的正负符号function phi_new upwind_levelset(phi, F, dx, dy, dt) [phix, phiy] gradient(phi, dx, dy); % 中心差分仅用于检测方向 % 使用简单迎风根据 F 的符号选择一侧差分 phixl gradient(phi, dx) ; % 默认前向 phixb gradient(phi, dx); % 占位 % 实际实现应如下 phi_x_p zeros(size(phi)); phi_x_m phi_x_p; phi_x_p(:, 2:end) (phi(:, 2:end) - phi(:, 1:end-1)) / dx; phi_x_m(:, 1:end-1) (phi(:, 2:end) - phi(:, 1:end-1)) / dx; phi_x_p(:, end) phi_x_p(:, end-1); phi_x_m(:, end) phi_x_m(:, end-1); % 同理 y 方向处理省略实际按 dx,dy 正确索引 grad_upwind sqrt(max(phi_x_p, 0).^2 min(phi_x_m, 0).^2 ... max(phi_y_p, 0).^2 min(phi_y_m, 0).^2); phi_new phi - dt * F .* grad_upwind; end上面这段代码是为了展示迎风的基本思想但其中 y 方向没有写完整完整版本应该同样处理phi_y_pphi_y_m。实际工程里我一般直接用文献中成熟的汉密尔顿-雅可比求解器比如 MIT 的ToolboxLS或者自己写一个 Lax-Friedrichs 通量function phi_new hamilton_jacobi_step(phi, F, dx, dy, dt) % 一阶 Lax-Friedrichs 格式速度 F 是标量函数 grad_x (phi([2:end, end], :) - phi([1, 1:end-1], :)) / (2*dx); grad_y (phi(:, [2:end, end]) - phi(:, [1, 1:end-1])) / (2*dy); laplacian_x (phi([2:end, end], :) - 2*phi phi([1, 1:end-1], :)) / dx^2; laplacian_y (phi(:, [2:end, end]) - 2*phi phi(:, [1, 1:end-1])) / dy^2; phi_new phi - dt * F .* sqrt(grad_x.^2 grad_y.^2) ... 0.5 * dt * max(F(:)) * dx * (laplacian_x laplacian_y); endLF 格式加入人工扩散项系数0.5 * dt * max(F)能有效抑制振荡代价是界面略微模糊。这是新手最容易理解也最容易改的格式。如果你想做高精度换成 WENO5但一维轮子已经够演示耦合。这个函数里有一个关键点对 F 的取值要么在节点上给定常数要么在每一步用有限元解插值得到。如果你的 F 在界面两侧不连续需要先做平滑否则 LF 格式的人工扩散会显得不足。3.4 主循环骨架与初始参数表主循环把上述所有函数串联起来。一个典型的主循环在 MATLAB 里长这样% 初始化 phi sqrt((nodes(:,1)-0.5).^2 (nodes(:,2)-0.5).^2) - 0.2; kappa phi_to_kappa(phi, 1.0, 0.2, epsilon); T zeros(size(phi)); for n 1:nt kappa phi_to_kappa(phi, k_solid, k_liquid, epsilon); K assemble_stiffness(mesh, kappa); [K, Fv] apply_dirichlet(K, zeros(size(phi)), fixed_nodes, fixed_values); T K \ Fv; % 求解稳态温度场 F get_velocity_from_temperature(T, phi, dx, dy, dt); phi hamilton_jacobi_step(phi, F, dx, dy, dt); if mod(n, reinit_freq) 0 phi reinitialize(phi, dx, dy); end if mod(n, 10) 0 fprintf(step %d, interface area %.4f\n, n, sum(phi 0)*dx*dy); end endget_velocity_from_temperature和reinitialize这两个函数每个写起来都有一页纸。reinitialize 最简单的实现是迭代求解 φ 的稳态方程直到收敛。但为了让你先跑通可以先用 MATLAB 自带的implicit_distance插件或者干脆跳过重新初始化只在你发现界面距离函数性质丢失时再补。这里给出初始参数的典型值参数值含义nx, ny40, 40单元数调大网格更密但时间成本上升dx, dy1/nx, 1/ny网格步长dtCFL * dx / max(F)时间步长epsilon1.5 * dx水平集平滑带宽reinit_freq5每 5 步重新初始化一次k_solid1.0界内热导率k_liquid0.2界外热导率CFL 数建议小于 0.5这样才能保证界面在一个时间步内不会跨过一个单元。这一步是最容易忽略的水平集方程是显式时间推进稳定性限制比有限元椭圆问题要严格得多。4. 算例调参与常见翻车现场从界面偏移到面积漂移4.1 一个具体算例初始圆在温差场中的演化在单位正方形内初始界面是一个圆心(0.5,0.5)、半径0.2的圆。设定左边界温度 1.0右边界温度 0.0上下边界绝热。材料热导率在圆内部取 1.0外部取 0.2。这样高温从左侧向右侧扩散靠近左侧的界面受热界面随温度场逐渐向左移动还是向右具体方向取决于你定义的 F 正负。这里我们约定界面外法向指向 phi0 区域F 为正表示界面收缩即圆的面积变小。运行nt 200步后你可能在figure; contour(reshape(phi,nx1,ny1), [0 0]); axis equal;里看到一个仍然近似圆但形状慢慢变的界面。这符合预期因为温度梯度在左右两侧不同界面移动速度不均匀。这个算例的价值在于验证程序能否稳定运行。如果 200 步内界面面积变化超过 10%说明你的 dt 太大或 reinit 过于频繁。下面就要讲如何调参。4.2 三个必调参数CFL、重新初始化频率、带宽第一个是 CFL。经验公式dt 0.4 * dx / max(abs(F))只是起点。如果你的速度场里有尖峰或者界面接近边界需要降到 0.2。MATLAB 里检查稳定性的最直接方法是在每个时间步输出max(abs(F)) * dt / dx如果这个数大于 0.5立刻把 dt 折半重跑。第二个是重新初始化频率。频繁重新初始化会让界面位置发生人为偏移因为距离函数重建过程本身也有数值误差。我一般先用 10 步一次再对比 5 步一次的结果。如果两者差异很小就用低频如果差异大说明 F 计算过度依赖 φ 的梯度需要降低 reinit 周期。判断指标是界面所围面积A sum(phi 0) * dx * dy理想情况 A 随时间缓慢单调变化。第三个是 epsilon。epsilon 是 Heaviside 的光滑宽度直接影响热导率过渡带。如果 epsilon 太小比如0.5*dx那么节处材料参数会接近阶跃有限元求解温度场时容易出现局部尖峰进而导致 F 在界面附近也不光滑如果 epsilon 太大比如5*dx界面被抹了半根指头宽F 的计算就失去物理意义。我的习惯是epsilon 1.5*dx做网格收敛性研究时同时按比例缩放 epsilon。下面用表格总结一下故障现象和调参方向现象可能原因先调哪个参数调完还不解决再看界面在 20 步内严重变形dt 太大CFL 降到 0.2检查 F 是否平滑面积单调漂移但形状正常reinit 太频繁或太少将 reinit_freq 改为 5/10 对比ε 是否合适温度场等值线在界面附近锯齿状ε 太小增大 epsilon 到 2.5*dx检查网格质量整体计算发散边界条件错误检查 apply_dirichlet 中罚因子检查 K 对角线是否非负4.3 界面模糊、面积漂移、迭代发散分别怎么调界面模糊是最常见的问题直观表现是contour(phi,[0 0])画出的界面线越来越粗或者界面附近的等值线不再重合。原因通常是 LF 格式的数值扩散太大。左上角的修法是降低人工扩散系数即把 LF 项中的0.5*dt*max(F)换成0.1*dt*max(F)。如果界面出现锯齿则说明人工扩散不足需要提高一点。面积漂移一般不是时间步造成的而是水平集方程和重新初始化之间的平衡破了。你可以做一个自检把 F 设为零只演化水平集跑 100 步观察初始圆形面积是否仍然守恒。如果 F 为零面积还在变化那一定是水平集求解器自身的守恒性问题。此时检查边界处理是否正确尤其是四个角点因为水平集方程的迎风差分在边界上需要外插。迭代发散的最常见原因是温度场求解失败。请你在主循环里加一条断言assert(isfinite(norm(T, inf)), 温度场出现 NaN 或 Inf);如果触发断言先检查边界条件是否挂载正确、K\F是否有警告。罚函数法的罚因子太大可能导致条件数爆炸通常1e10已经足够继续加并不会更准确反而会丢失精度。5. 验证和加速水平集-FEM 程序的两条捷径5.1 用守恒性检查做自洽验证任何数值代码在深入研究之前都应该先跑自洽验证。对水平集-FEM 程序来说最容易做的验证是“零速度守恒”。把 F 设为常数零更新水平集方程界面应完全不动。如果动说明时间积分或边界处理有 bug。第二个验证是圆界面在均匀外法向速度下保持圆形。令 F 1初始圆半径 R00.2理论半径 R(t)R0t。程序运行 20 步后计算实际半径用phi0的节点面积等效相对误差小于 1% 即为正常。这里有现成命令R_measured sqrt(sum(phi(:) 0) * dx * dy / pi); R_expected sqrt(pi * 0.2^2 t * pi) ; % 注意公式取决于面积增长率不要只看半径速度场为零时面积守恒是最严格的测试比半径更敏感。5.2 用 parfor 和稀疏重构加快调试迭代MATLAB 程序中耗时最多的是每步重新组装刚度矩阵。如果你的网格是均匀矩形你可以提前把所有单元的几何信息面积、梯度矩阵缓存下来在循环里只更新kappa并重新计算Ke ke * grad*grad*area。这样能省掉最耗时的坐标提取和面积计算。在装配循环里如果把for e 1:Ne换成parfor注意每个e迭代必须完全独立。上面代码里的K(idx, idx) K(idx, idx) Ke是一种 accumulate 操作在 parfor 里不能直接累加。常见做法是预先分配三个向量row,col,val在每个迭代里填充本单元贡献的三行再用一次sparse(row, col, val, N, N)构造整体矩阵。这样既快又符合 parfor 要求。如果你想走得更远可以把assemble_stiffness里最耗时的循环改写成向量化矩阵运算。对于均匀网格所有单元的面积和梯度矩阵是相同的可以直接把每个单元的节点索引做成一个Ne x 3的矩阵用accumarray一次性把局部刚度矩阵的 9 个分量累加进全局索引。这种写法能把 400x400 网格的装配时间降低一个数量级但代码可读性差建议只在验证通过之后做优化。最后一招在调用K\F之前先对 K 做一次符号分解Ks decomposition(K, lower); T Ks \ F;如果边界条件固定且材料参数每步都变分解矩阵无法复用但如果 kappa 变化很小可以每 10 步才更新一次预分解。这是一个工程取舍不是通用准则。另外水平集初始化的距离函数重建在 MATLAB 中常调用dist2p函数但它在大型网格上很慢。更快的做法是调用scipy.ndimage或 C 扩展但既然已经锁定在 MATLAB 环境我建议对中等规模网格100x100 以内坚持用迭代法重新初始化并记录每次迭代的 max 差量差量稳定后提前退出能节约一半时间。本文还有配套的精品资源点击获取