简介这套MATLAB/Octave代码包围绕各向异性网格自适应、拓扑优化与有限元方法展开面向固体力学、热传导、结构优化等方向的科研人员与工程师特别适合需要借助数值仿真探索最优材料布局的读者。代码将自适应网格技术与拓扑优化求解相结合能在二维和三维模型中高效处理最佳导热路径设计也可用于线性弹性最小柔度优化问题方法的理论背景可追溯到WCSMO12会议论文具备较强的可复现性。压缩包共5个文件主体为3个m脚本分别承担热传导与弹性优化等经典算例另附1个rst格式说明文档和1个md格式开源许可整体大小仅约173KB结构紧凑便于快速阅读和二次开发。已有323人浏览/学习对于希望掌握拓扑优化实现细节、开展算例复现或扩展算法功能的MATLAB/Octave用户来说是一套轻量且实用的工具包。1. TrulleKrul能做什么把网格自适应从“手动调”变成“自动长”做拓扑优化时网格自适应常常是最容易被当成“分支任务”的一环。TrulleKrul这套MATLAB/Octave脚本核心是把各向异性网格自适应、拓扑优化和有限元方法串成一条自动链路。它的价值不是给你一个更复杂的求解器而是让你在求解过程中自动调整单元的疏密、形状和拉伸方向让应力集中区、边界层和优化结构界面获得足够分辨率又不必让整体自由度爆炸。适合谁手里已经有一套SIMP或Level Set拓扑优化代码、但网格一直靠手动加密和重划分的开发者也包括刚接触拓扑优化、想用轻量脚本理解“误差估计网格重构场插值”这套完整机制的研究生。它不是银弹但如果你愿意花一个下午把最小流程跑通后面省下的调网格时间会相当可观。2. 各向异性网格自适应原理TrulleKrul为什么要用“长条单元”而不是单纯加密2.1 各向异性网格和普通加密的本质区别普通自适应加密遇到应力集中就把单元切成更小的同类三角形单元长宽比保持接近1。这样做有一个麻烦在弯曲边界层或薄板结构里应力往往只在法线方向变化剧烈切向方向变化平缓。用等向单元去细分切线方向本质上是在浪费自由度。各向异性网格的思路是允许单元在一个方向拉长、另一个方向压缩让单元形状追踪解的曲率。TrulleKrul走的正是这条路线所以它适合的问题是明确的细长边界、局部接触、壳结构里的弯曲层以及拓扑优化里从灰色过渡到黑白的那条窄带边界。从实现角度看各向异性网格自适应不是单纯“加密某些单元”它包含两个动作。第一是估计判断哪个位置需要更细的单元以及哪个方向需要更细第二是重构在给定“理想单元形状”约束下重新生成网格。TrulleKrul把这两件事拆成了独立模块用户只改误差指示器不影响网格生成部分。也是因为这个拆分它作为教学代码比很多商业化黑匣子更友好——你可以看到密度场或应力场如何一步一步变成网格尺寸场。需要澄清一个常见误解各向异性不等于边界层加密。边界层加密只是在一个固定层数里把单元压扁而各向异性自适应是在每个点用误差估计器自主决定“哪个方向扁、扁多少”。前者需要人工指定层数后者只需要给一个目标误差或拉伸上限网格自己会找位置。这也是TrulleKrul和普通网格细化工具最大的差别。2.2 误差估计与黎曼度量TrulleKrul的核心逻辑TrulleKrul这类脚本普遍采用基于解的Hessian矩阵来构造度量张量。设有限元解为u插值误差沿某个方向h的估计可以写成|u - u_h| ≤ c * max_s (hᵀ |H| h)这里的H是由节点值恢复出的二阶导数矩阵也就是解的曲率。如果某个方向曲率大h就必须收缩曲率平缓的方向h可以放大。把H做特征值分解得到三个正交方向和期望尺寸再由目标误差算出整体拉伸比例。TrulleKrul内部把这些信息整理成每个网格点上的度量张量M网格生成器读到的不是“这个地方要加密”而是“在这个点附近沿方向D₁边长0.02、沿方向D₂边长0.18”。这套思路的边界也很清晰如果解本身不够光滑Hessian恢复出来全是噪声度量张量就会乱。所以几乎所有实现里都有一步“平滑”或“截断”。TrulleKrul一般通过设定hmin、hmax以及最大拉伸比来控制度量张量的特征值范围。换句话说你要调的不是某个神秘的“网格密度系数”而是这三个有物理意义的数字。我习惯在调用自适应之前先把网格点上的度量特征值分布画出来确认最大拉伸比在合理区间。最常遇到的问题是最小特征值被某个噪声点拉得很小导致那一小片区域单元被极限压扁网格规模直接跳涨一个量级。遇到这种情况不要急着改网格生成器先看Hessian恢复处的光滑性或者把hmin上调两倍问题通常会缓解。2.3 一套最小自适应流程的可执行骨架不需要先读源码也能把TrulleKrul的思路跑通。我建议把流程写成下面这样的封装solve → estimate → adapt。之后无论是换求解器还是换误差指示器都只改其中一个函数。function [p, tri, info] run_trullekrul_adapt(problem, opt) % 最小自适应循环solve - estimate - adapt p problem.p0; tri problem.tri0; for it 1:opt.numit u solve_fe(p, tri, problem); % 1. 当前网格上求解 H recover_hessian(p, tri, u); % 2. 恢复节点二阶导数 metric metric_from_hessian(H, opt); % 3. 由H构造度量张量 [p, tri, info(it)] adapt_by_metric(p, tri, metric); % 4. 重构网格 end end关键点在于solve_fe可以是任意的有限元求解函数recover_hessian可以用超收敛补片恢复也可以简单地对解求两次梯度后用L2投影光滑metric_from_hessian负责把H的特征值映射到目标边长adapt_by_metric对应TrulleKrul内部的重划分调用。把这个骨架钉死后你能在未知内部实现的前提下先手动把每一步打出来看结果。参数opt里我一般只设四个字段hmin、hmax、stretch_max、numit。hmin是全局最小单元边长hmax是最大边长stretch_max是最长边与最短边的比值上限numit是自适应循环次数。第一次调试时只动这四个其他参数一律不动否则问题会被参数之间的相互作用掩盖。需要注意这个封装里的adapt_by_metric是抽象函数。TrulleKrul不同版本的函数名或有差异但骨架不变。你接手一个更新版本时唯一要做的就是去解析包里自带demo找到那行真正调用自适应模块的位置再用if/else把它替换到adapt_by_metric里。这个替换过程不超过半小时却能替后面省下大量排查时间。3. 在MATLAB/Octave里跑通最小示例一条命令跑完三分钟看结果3.1 环境准备路径、工具箱和Octave兼容性TrulleKrul是纯MATLAB/Octave脚本意味着不需要编译。第一件事是解压后把整个目录放进路径。我习惯在项目根目录放一个startup.m内容如下% startup.m addpath(fullfile(pwd, trullekrul)); addpath(fullfile(pwd, solvers)); addpath(fullfile(pwd, meshers));这里把三个子目录都加进来核心脚本、有限元求解器、网格初始化和重划分工具。不同版本子目录名可能不同但思路一致。Octave下设置路径的方式一样只是解压目录后要先确认没有需要编译的MEX文件。如果包里带了MEXOctave需要对应的编译环境否则建议找到同功能的纯M文件版本避免一开始就踩编译坑。性能方面要提前有心理预期Octave的循环解释开销比MATLAB大TrulleKrul这种反复做网格重构的脚本在Octave下容易慢五倍以上。如果你是Octave用户先跑通二维小规模算例不要一上来就三维。三维各向异性网格重构的内存占用比二维高一个量级Octave在单元数超过几万后会有明显卡顿。后面第5章会具体说怎么缓解。3.2 最小示例二维泊松问题的自适应网格不要用包内demo作为第一个跑的项目demo往往叠加了后处理和绘图逻辑混淆视线。最好自己写一个脚本二维泊松方程源项集中在中心一个很小的峰值上跑三次自适应。% demo_minimal.m % 2D Poisson problem with TrulleKrul-style adaptive loop clear; clc; % 1. build initial mesh on unit square, deliberately coarse [p, tri] init_rectangle(0, 1, 0, 1, 0.2); % 2. problem definition problem.f (x,y) 20 * exp(-((x-0.5).^2 (y-0.5).^2) / 0.002); problem.g (x,y) 0; % Dirichlet boundary % 3. adaptation parameters opt.hmin 0.005; opt.hmax 0.2; opt.stretch_max 15; opt.numit 3; % 4. run [p, tri, info] run_trullekrul_adapt(problem, opt); % 5. visualize trimesh(tri, p(:,1), p(:,2)); axis equal;这段代码最值得注意的是hmin0.005。假设有人为了“精度”把hmin调到1e-4中心峰周围单元数量就会爆炸矩阵组装和网格重构都会卡在几分钟级别。正确做法是先从粗网格开始让自适应自己决定是否细化。stretch_max15是各向异性的上限太低会退回普通自适应太高则可能出现畸形单元。二维问题从10到20开始调是稳定的。另一处容易误用是把hmax设成与hmin相等。那样所有单元的尺寸都被限制在同一个区间各向异性网格的优势完全发挥不出来。hmax应该至少比hmin大一个数量级网格才有空间在平缓区域保留大单元。3.3 三个必调参数hmin、hmax和stretch_maxTrulleKrul里有一堆参数但真正影响结果的是下面三个它们直接决定计算量和可靠性。参数表如下。参数含义参考范围调参建议hmin最小单元边长防止局部过细目标特征尺寸的1/101/20先固定不要一开始就取极端小值hmax最大单元边长控制整体规模结构宏观尺寸的1/51/10想整体加密就调小它stretch_max最大拉伸比限制各向异性程度二维820三维38出现负Jacobian就先调低还有一个隐藏参数常被忽视度量张量的截断因子你可以把它理解成自适应“敏感度”。建议不要上来就碰它。先把上面三个数调到结果合理再根据“哪些区域还在过度加密”去调截断系数。否则参数之间的相互作用会让排查变得很痛苦。如何判断调参是否合理一个简单检验是在stretch_max10和stretch_max20下各跑一遍记录最终单元数和求解误差。如果单元数减少30%但误差反而变大说明拉伸比已经越界。这时候不要继续增大要么调低stretch_max要么提高hmin让整体单元变细。我见过很多人在这个点上反复折腾最后发现是hmin设得太小导致度量张量截断完全失效。跑通第一个算例后我建议立刻做一组对照把stretch_max改成1重新跑一遍。这样你能直观看到等向网格和各项异性网格的单元数量差别。很多时候这一个对照实验比读十遍文档更能理解TrulleKrul到底在干什么。4. 把TrulleKrul嵌进SIMP拓扑优化三步封装避免“网格一改、目标函数就跳变”4.1 为什么网格自适应能帮拓扑优化收敛拓扑优化过程中密度场从均匀分布进化到0/1分布结构边界从模糊区域逐渐清晰。如果固定网格为了捕捉边界通常要把整个设计域画得很细但细网格意味着每次有限元求解成本大增。各向异性自适应能让边界附近的单元沿边界方向拉长、沿法线方向压薄边界分辨率高但自由度少。另一个层面上网格自适应还能让优化中期密度过渡区更平滑。基于当前密度场重新剖分后灰色区域的单元变细灵敏度过滤和密度投影结果更干净。但自适应不能每迭代步都做否则密度场在旧网格上积累的演化记忆被网格重构打断目标函数曲线会出现锯齿状跳变。这是TrulleKrul接入优化时最需要控制的一点。4.2 最小嵌入代码从固定网格改成自适应网格下面这段代码是把经典的SIMP拓扑优化改成TrulleKrul自适应版的骨架省略了具体板结构和载荷细节重点看网格重构部分。% simp_adapt.m % SIMP topology optimization with periodic TrulleKrul remeshing nelx 60; nely 20; volfrac 0.4; penal 3; rmin 1.5; maxiter 80; adapt_every 20; % initialize density and design variables x repmat(volfrac, nely, nelx); xPhys x; mesh build_uniform_mesh(nelx, nely); % initial uniform mesh for iter 1:maxiter % solve on current mesh [u, mesh] fe_solve(mesh, xPhys); % sensitivity from SIMP dc simp_sensitivity(mesh, u, xPhys, penal); dc sensitivity_filter(mesh, dc, rmin); % OC update xNew oc_update(x, dc, volfrac); xPhys density_projection(xNew, penal); if mod(iter, adapt_every) 0 % build metric from density gradient and stress metric metric_from_fields(mesh, xPhys, u); % remesh meshNew trullekrul_adapt(mesh, metric); % transfer density from old mesh to new mesh xPhys transfer_volume_weighted(mesh, meshNew, xPhys); % reset design variable to match new discretization x xPhys; mesh meshNew; end end最关键的是transfer_volume_weighted这一步。很多人忽略它直接把旧网格密度当成新网格密度结果自适应一次后结构轮廓全花掉。正确做法是对每个新网格单元找到它与旧网格单元的覆盖关系用面积加权平均得到密度值。线性插值为什么不行因为拓扑优化的密度场经过Heaviside投影后大量值接近0或1线性插值会把边界上的0/1抹成0.5的灰色云团。rmin过滤半径也要留意。如果过滤半径以物理长度为单位网格变化后影响不大如果以单元数为单位网格加密后同一个半径覆盖的单元数量会完全不同灵敏度分布会被悄悄改变。用TrulleKrul时建议把rmin改写成物理长度才不会被网格尺寸变化干扰。4.3 自适应节奏与密度传输每隔多少步重构一次我的习惯是每20到40步重构一次而不是每步都做。更稳的节奏是“三段式”前40步固定网格让拓扑基本成型中间每20步自适应一次让边界变清晰最后20步固定新网格只细化边界附近的单元。如果最后的细化不需要改变拓扑甚至可以把TrulleKrul的度量输入从“全场密度梯度”切换成“边界距离场”让加密只发生在结构边界。密度传输的实现可以单独封装下面是一个面积加权版本的核心逻辑function xNew transfer_volume_weighted(meshOld, meshNew, xOld) % volume-weighted density transfer from old to new mesh xNew zeros(size(meshNew.elem,1), 1); for e 1:numel(xNew) cen meshNew.elemCenter(e,:); % find covering old element by centroid location oldElems pointLocation(meshOld.tri, cen); if isnan(oldElems) % fallback: nearest old element oldElems nearestNeighbor(meshOld.dt, cen); end xNew(e) xOld(oldElems); end end这段代码在二维三角网格里用pointLocation找覆盖单元找不到时回退到最近邻。这里只写了面积加权的简化版本实际更稳健的处理是遍历新网格单元计算它被哪些旧网格单元覆盖再按交集面积权重平均。后者复杂度更高但对边界上的密度保持更友好。从工程实践看如果传输后的密度分布与传输前相比总体积偏差超过1%说明网格变化太剧烈。TrulleKrul并不是每步都要把网格彻底改头换面它也可以只做局部重划分。在优化中后期把自适应范围限制在密度介于0.1到0.9的过渡带附近能明显减少插值误差。5. TrulleKrul踩坑实录网格爆炸、负Jacobian与灵敏度跳变5.1 现象初始网格很细程序卡在第一次自适应上我第一次接手TrulleKrul时用了0.01的初始网格尺寸觉得这样“精度够”结果自适应循环跑了二十分钟还没结束。原因是初始网格过细第一次误差估计后峰值周围又加密一轮单元数量直接翻了几十倍。原因hmin和hmax是全局范围参数初始网格尺寸不在这个范围内就会被度量张量继续强制加密。解决从比hmax更粗的网格起步让自适应自己去“生长”细网格而不是一开始就把全局网格铺细。我后来把初始尺寸设为hmax的两到三倍程序通常在几次迭代内就能稳定下来。5.2 现象自适应后单元畸形求解器报负Jacobian报错信息里出现negative Jacobian时几乎可以断定是某个单元发生了翻转。各向异性强度过大时单元在拉伸方向被压得太扁两个节点交叉导致面积变成负数。原因stretch_max设得过高或度量张量没有做足够的正则化。解决先把stretch_max降到8试二维算例再逐步往上加。同时检查Hessian恢复是否引入了负特征值。如果recover_hessian没有做特征值截断负特征值会让度量张量出现“负尺寸”区域这比单纯的大拉伸比更容易让网格翻车。5.3 现象拓扑优化迭代几步后密度场变成一片灰色结构“融化”SIMP优化在固定网格上能收敛到0/1但接入TrulleKrul后每20步网格重构一次结果密度场越来越灰最后整个设计域呈半透明状。原因网格重构后使用了线性插值传输密度0/1边界被抹成中间值。解决改用上一节说的面积加权传输并且把密度限制在[0,1]。如果边界仍然模糊检查自适应频率是否过密。网格重构后需要给OC算法至少5到10步去重新锐化边界20步一次通常不会出大问题10步一次就会偏灰色。5.4 现象Octave下跑同样算例比MATLAB慢五倍以上TrulleKrul在Octave里的性能落差主要是两类开销一是循环解释慢二是稀疏矩阵求解器不如MATLAB稳定。原因Octave对大量小循环和动态增长数组不友好而TrulleKrul这类自适应代码正好充满这类操作。解决提前预分配所有数组尽量用向量化替换for循环把自适应频率降低从每10步一次改成每30步一次求解器部分使用lu显式分解并复用符号分解。实测中仅把矩阵组装的稀疏索引从动态增长改成预分配Octave性能就能提升一倍以上。5.5 现象相邻两次自适应的目标函数发生“跳水”优化曲线断裂自适应后的第一步目标函数突然比上一步低了15%然后又开始缓慢上升整个优化曲线像断了一截。原因新旧网格离散误差不同目标函数基于两套不同网格计算数值基准发生变化。解决自适应后先固定网格跑若干步再记录曲线。更严格的做法是把每次自适应后的目标函数用旧网格插值重新估算得到一个校正系数让曲线连续。工程师更常用的是“不比较网格重构前后的目标函数绝对值只比较重构后5步内的变化趋势”。6. 用两个小习惯验证TrulleKrul的结果质量直方图和灵敏度一致性测试6.1 每次都画网格质量直方图TrulleKrul跑完不要只看结果图就收工。我习惯先把所有单元的边长比算出来画一张直方图。正常情况下分布会集中在1到stretch_max之间越靠近上限说明各向异性利用越充分。如果直方图尾部出现个别长宽比超过100的单元说明度量张量里有尖峰没被截断。这时不要去找网格生成器的问题先回头看hmin和Hessian恢复是否产生了异常极值。质量直方图是一个极好的早期预警器。6.2 用灵敏度一致性测试判断该不该继续信任新网格另一个验证技巧是灵敏度一致性测试。做法是在自适应前计算一次目标函数和灵敏度然后把密度传输到新网格在新网格上重新计算目标函数和灵敏度。如果两者相对偏差小于2%说明网格重构没有剧烈改变离散误差如果偏差超过5%说明这次自适应太激进需要降低stretch_max或减少自适应覆盖范围。这个测试不需要额外写很复杂的代码只需要在自适应函数前后各调用一次求解器和灵敏度计算。长期观察下来TrulleKrul的结果可信度和这个偏差大小高度相关。偏差小的算例后处理出来的边界提取和应力校验都干净偏差大的算例往往会在后续强度校核里翻车。我现在的固定习惯是每次跑完一个TrulleKrul算例先出质量直方图再做灵敏度一致性测试两个关卡都过了才敢把优化结果交给下游做几何重建。看起来多花了十分钟实际上省掉了后面反复校验的半天时间。希望帮到你。本文还有配套的精品资源点击获取