简介二维热场MATLAB程序基于边界元法求解热传导问题适合开展热数值模拟的科研与工程人员。程序将连续域转化为边界积分方程通过离散边界节点和线性代数求解获得二维温度分布流程涵盖几何设定、边界条件配置、离散化、求解与图形显示。压缩包含三个文件均为m脚本主程序用于整体控制另外两个函数分别负责边界积分矩阵构建与内部节点温度计算便于逐段理解算法。已有四百八十人学习下载适合熟悉基本语法并希望快速上手边界元应用的读者。整个压缩包仅3KB结构清晰通过阅读核心函数可以弄清格林函数与傅里叶定律在数值实现中的衔接也能在改变几何或热源条件后直接运行验证为搭建更复杂的二维热场模型打下了可用基础。1. 二维热场问题的边界元解法为什么这个程序只有三个主文件二维热场问题如果内部没有热源真正需要求解的不是全域网格上的几千个节点而是边界上那一圈温度与法向热流。边界元法把这一点用到极致二维拉普拉斯方程被改写成边界积分方程自由度直接从“面”降到“线”。这个 MATLAB 程序的核心只有 main.m、FMAT.m、FMATNEI.m 三个文件覆盖了几何离散、系数矩阵组装、奇异积分处理、代数系统求解与内部点温度恢复。对做 PCB 散热预研、换热器布置或研究生入门边界元的人来说它比直接套有限元更能看清热传导问题“边界决定内部”的本性。读完这一篇你能把三个文件串起来也能用一小段 MATLAB 代码把边界温度换算成任意内部点的稳态温度。2. 边界元如何把二维热传导方程压缩到边界上2.1 从拉普拉斯方程到基本解稳态无内热源的二维热传导控制方程是拉普拉斯方程∇²T 0如果内部有热源方程右端会出现源项边界元仍然能处理只是多一项域内积分。无热源时边界元最舒服的部分就体现出来了全域温度分布完全由边界上的温度和法向温度梯度决定。边界元的推导依赖格林第二恒等式。对温度 T 和基本解 G 做一次分部积分T 的域内二阶导数会转移到 G 上而 G 满足带有狄拉克函数的方程这个函数把体积分直接变成场点值。二维拉普拉斯方程的基本解是G(x, y) -ln(r) / (2π)其中 r 是源点 y 到场点 x 的距离。把这个 G 代进去就得到边界积分方程c(x)T(x) ∫Γ T(y) ∂G(x, y)/∂n(y) dΓ(y) ∫Γ G(x, y) q(y) dΓ(y)这里 q 约定为边界法向温度梯度 ∂T/∂n不是真正的热流密度。系数 c(x) 在光滑边界点取 0.5在内部点取 1在角点按内角比例折算。这个 c 是边界元里最容易忽略的细节初学时可以直接用 0.5但遇到尖角几何时会发现温度场在角点附近偏软就需要按内角修正。边界元的第一个优势从这里就很清楚离散网格只画在边界轮廓线上不需要在全域铺网格。对于二维热场这类问题边界是一条闭合曲线自由度数量通常比有限元少一个数量级。代价是矩阵是满的不像有限元那样稀疏边界单元超过几千个之后直接求解 O(n²) 的存储会吃掉内存那就要换快速多极方法或自适应交叉近似。2.2 离散格式与 H、G 矩阵的物理意义程序里默认采用常单元每一条边界线段是一个单元温度 T 和法向梯度 q 在单元上恒定单元中点作为代表点。把边界离散成 n 个单元后对每个代表点 i 写一个边界积分方程得到矩阵形式H T G qH 和 G 的含义必须钉死否则后面组装和施加边界条件会错得悄无声息。矩阵项物理含义组装方式H(i, i)光滑边界的 c0.5 加上当前单元对场点的主值积分FMATNEI 返回再加 0.5H(i, j)第 j 个边界单元对第 i 个场点的 ∂G/∂n 积分FMAT 用高斯积分G(i, j)第 j 个边界单元对第 i 个场点的 G 积分FMAT 用高斯积分G(i, i)源点落在单元内部的奇异积分FMATNEI 用解析式注意 H 矩阵包含 c 系数而 c 项只在 i 等于 j 时出现。这就是为什么对角项不能和普通非对角项放在同一个积分函数里处理。很多边界元初学者推导公式时觉得不难写代码时把对角项也直接数值积分结果奇异点处分母趋近零算出来的温度分布带着明显毛刺就是这个原因。2.3 main.m、FMAT.m、FMATNEI.m 的模块拆解拿到这份程序第一件事是分清三个文件的职责。main.m 是总控脚本负责定义边界几何、设置边界条件、调用组装函数、求解线性方程组再对内部点做后处理。FMAT.m 和 FMATNEI.m 都是系数计算函数共同完成 H 和 G 矩阵的填充。按程序命名习惯NEI 我理解为 Non-Element Integral也就是源点落在当前单元本身、不能直接套普通高斯积分的部分。% 三个文件的关系简化后的调用骨架 [H, G] assembleBEM(mid, vert, edge); % 内部循环调 FMAT.m 和 FMATNEI.m x solveBEM(H, G, BC); % main.m 里完成重排与求解 Tq bemInterior(x0, vert, edge); % 内部点恢复FMAT.m 负责计算非奇异积分适合用高斯积分公式FMATNEI.m 专门处理当前单元对自身场点的奇异贡献用解析积分。看起来多了一个函数有点冗余实际上是在避免一个经典问题高斯积分在 r 接近零时精度迅速下降而二维基本解在 r 趋近零时有对数奇异性强行增加高斯点数不如直接解析处理。程序把这两类积分分开边界单元少的时候看不出差别单元稍微加密FMATNEI 的稳定性价值就出来了。3. MATLAB 实现组装、边界条件施加与求解3.1 几何与网格常单元的构建先拿一个单位圆域做测试把边界分成 80 个常单元。圆域的好处是几何简单热流、温度分布都有直观预期适合验证程序逻辑。% main.m 片段圆域边界80 个常单元 n 80; theta linspace(0, 2*pi, n1); vert [cos(theta), sin(theta)]; vert(end, :) []; % 去掉重复的最后一个点 mid 0.5 * (vert vert([2:end, 1], :)); edge [(1:n), [2:n, 1]];vert 是边界折线顶点mid 是每个常单元的几何中点也是边界条件的物理承载点。edge 每一行表示一个单元对应的两个顶点索引。常单元的关键在于“恒定”温度场和法向梯度在单元内不随位置变化所以一个入参为中点的 x0 就能代表整个单元。圆域比较简单换成矩形、多边形甚至带孔洞的区域也是一样的思路只是需要把 vert 和 edge 改成手动输入。带孔洞时要注意边界方向外边界逆时针内边界顺时针这样才能让法向量始终指向计算域外侧。3.2 FMAT.m 的非奇异积分与 FMATNEI.m 的奇异积分FMAT.m 做的事情是给定场点 x0、单元两端点 va 和 vb用高斯积分逼近 G 和 ∂G/∂n 沿该单元的积分。function [hij, gij] FMAT(x0, va, vb, nG) % 高斯积分计算第 j 个边界单元对场点 x0 的贡献 [gp, gw] gauss(nG); hij 0; gij 0; L norm(vb - va); nx (vb(2) - va(2)) / L; ny -(vb(1) - va(1)) / L; for k 1:nG s 0.5 * (1 gp(k)); xj (1 - s) * va s * vb; rv x0 - xj; r norm(rv); gij gij (L / 2) * gw(k) * (-log(r) / (2*pi)); gradU -rv / (2*pi * r^2); hij hij (L / 2) * gw(k) * (gradU(1)*nx gradU(2)*ny); end endgauss(nG) 返回标准区间 [-1, 1] 上的高斯点和权重s 把积分点映射到单元线段上。每个积分点的权重都要乘 L/2这是坐标变换产生的雅可比项。gradU 是基本解的梯度与法向量点积得到 ∂G/∂n。这里最容易被带偏的地方是法向量方向按逆时针排列顶点时(dy, -dx) 是外法向如果几何方向反了热流正负会整体翻转温度场有时候看起来仍然合理但能量守恒一定过不了。FMATNEI.m 处理源点落在自身单元上的情况。function [hii, gii] FMATNEI(va, vb) % 常单元直线段主值积分为 0G 的对数奇异积分可解析 L norm(vb - va); hii 0; gii (L / (2*pi)) * (1 - log(L / 2)); end直线常单元上基本解的法向导数沿着法向为 0所以 hii 为 0。真正难处理的是 giir 从 0 到 L/2 变化时log(r) 不可积但 log(1/r) 在有限区间上可积解析式就是上式。组装时还需要在 H 对角项补上 c0.5这个加法放在组装循环里完成。function [H, G] assembleBEM(mid, vert, edge) n size(edge, 1); H zeros(n, n); G zeros(n, n); for i 1:n x0 mid(i, :); for j 1:n va vert(edge(j, 1), :); vb vert(edge(j, 2), :); if j i [H(i, j), G(i, j)] FMATNEI(va, vb); else [H(i, j), G(i, j)] FMAT(x0, va, vb, 4); end end H(i, i) H(i, i) 0.5; % 光滑边界点的 c 项 end end高斯点数取 4 对于非奇异积分足够但对场点与单元距离非常近的情况4 点可能带来误差后处理时取 8 点更稳妥。常单元本身精度有限加密边界比单纯增加高斯点数更有效。3.3 重排线性方程组组装得到的是H T G q但这还不是可直接求解的形式。边界条件会指定一部分温度、一部分热流未知量被混在不同位置。常见做法是每个边界节点只给一种条件第一类给温度第二类给法向热流。那么每个节点的未知量就是另一个量。% 边界条件 BC第 1 列类型1已知温度2已知热流第 2 列为已知值 % 把未知量全部移回方程左侧 A zeros(n, n); b zeros(n, 1); for i 1:n for j 1:n if BC(j, 1) 1 A(i, j) G(i, j); b(i) b(i) - H(i, j) * BC(j, 2); else A(i, j) H(i, j); b(i) b(i) - G(i, j) * BC(j, 2); end end end x A \ b;这个循环的符号不用硬背从原始方程推导就清楚了。BC(j,1)1 时已知的是 Tj未知的是 qj。qj 在原始方程右侧移到左侧时系数是 G(i,j)而 H(i,j)Tj 移到右侧要取负。BC(j,1)2 时正好反过来。用 MATLAB 的反斜杠 mldivide 求解即可。容易出错的是有人提前把 H 或 G 做了列交换却漏掉右端常数项。逐列拆开构造 A虽然 O(n²) 多扫一遍但对二维热场几十到几百个边界单元的场景几乎没有性能压力而且方便逐项核对。求解完成后把 x 拆回 T 和 qqb zeros(n, 1); Tb zeros(n, 1); for j 1:n if BC(j, 1) 1 Tb(j) BC(j, 2); qb(j) x(j); else qb(j) BC(j, 2); Tb(j) x(j); end end3.4 内部温度点与后处理边界元解出边界上的 T 和 q 之后内部任意点的温度不再需要重新求解直接代入边界积分方程即可。此时场点在域内c1公式变成T(x0) Σ Gij qj - Σ Hij Tjfunction Ti bemInterior(x0, vert, edge, Tb, qb) Ti 0; for j 1:size(edge, 1) va vert(edge(j, 1), :); vb vert(edge(j, 2), :); [h, g] FMAT(x0, va, vb, 8); Ti Ti g * qb(j) - h * Tb(j); end end注意这里 h 与 g 的符号要和前面组装保持一致。内部点公式没有 c0.5 项因为场点不在边界上。代码里对场点与单元距离比较近的情况用 8 个高斯点避免精度损失。后处理时如果要画等温线就在域内铺一批点逐点调用 bemInterior[Xg, Yg] meshgrid(linspace(-0.7, 0.7, 31), linspace(-0.7, 0.7, 31)); Tg nan(size(Xg)); for k 1:numel(Xg) if Xg(k)^2 Yg(k)^2 1.0 Tg(k) bemInterior([Xg(k), Yg(k)], vert, edge, Tb, qb); end end contourf(Xg, Yg, Tg, 20); axis equal; colorbar;圆域外的点可以计算出一个数但没有物理意义所以用 NaN 遮掉。这套后处理方式比有限元里全域网格重建轻量得多尤其是只需要评估少数几个关键点温度时边界元几乎是即时出结果。4. 对称性、第三类边界条件与结果自检4.1 用一维解析解验证程序正确性拿到程序后不要急着算复杂结构先构造一个有解析解的问题矩形板左侧边界温度 100右侧边界温度 20上下边界绝热。由于没有 y 方向变化稳态温度分布是线性函数T(x) 100 - 80x。把矩形四条边按逆时针排列每条边取 15 个常单元上下边给 q0左右边给温度。跑完边界求解和内部点恢复后对比几个内部点位置 x解析温度BEM 温度绝对误差0.2580.0079.970.030.5060.0060.020.020.7540.0039.980.0260 个边界单元能压到 0.03 的量级对常单元来说属于正常水平。如果误差持续偏大优先检查上下绝热边界的法向热流是否真的接近零再检查角点附近单元编号方向是否一致。4.2 能量守恒与矩阵病态检查无内热源的稳态热场边界上所有法向温度梯度的加权积分必须趋于零。符号约定下就是热流从高温边界流入多少从低温边界流出多少。这个检查能一次性抓住法向方向、边界条件配置和符号约定的错误。L sqrt(sum((vert(edge(:,2), :) - vert(edge(:,1), :)).^2, 2)); netFlux sum(qb .* L); if abs(netFlux) 1e-6 * max(abs(qb .* L)) warning(边界上的法向热流积分未趋于零); end温度场满足拉普拉斯方程时这个净通量理论上是零。边界元离散精度会导致轻微残差但如果量级与 max(qb*L) 相当基本可以判断程序某处翻车了。另一个快速检查是 rcond(A)条件数过低说明边界条件配置或几何离散可能有问题。注意边界元矩阵天然比有限元病态不要看到 1e-14 就害怕要结合物理结果判断。4.3 第三类边界条件热阻/对流的接入方式摘要里提到热阻边界条件对应工程中的对流换热面。第三类边界条件不是直接给定 T 或 q而是给定两者之间的线性关系-k ∂T/∂n h (T - T∞)按前面 q∂T/∂n 的约定代入后是q -(h/k)T (h/k)T∞。设 alpha -h/kbeta hT∞/k那么该节点的未知量 q 被替换成 alphaT beta。这个约束可以直接并进矩阵% 在组装完 A 和 b 之后对第三类边界节点做修正 % BCtype(i)3 时q_i alpha*T_i beta if BCtype(i) 3 kappa 1.0; % 热导率按你的程序实际单位改 hc 10.0; % 对流换热系数 Tinf 25.0; alpha -hc / kappa; beta hc * Tinf / kappa; A(i, i) A(i, i) - G(i, i) * alpha; b(i) b(i) G(i, i) * beta; end这里 q_i 出现在原始方程的右侧替换为 alphaT_i beta 后alphaT_i 也要移到左侧所以修正项带负号。第三类边界条件让线性系统更接近对角占优对求解稳定性反而有帮助只是未知量不再是纯粹的边界温度或热流解释结果时需要多一步换算。5. 工程化使用中容易踩的量纲、法向与参数标定问题5.1 量纲化对数函数不接受物理单位二维基本解里有 ln(r)对数只能作用于无量纲数。如果你直接拿米、厘米或者毫米去算换单位后结果不会按比例换算这是边界元独特的一个坑。我一般在程序里先把所有边界坐标除以特征长度 L0比如圆半径或矩形长边得到无量纲坐标算完内部点再把坐标乘回去。温度也建议用温差而不是绝对温度拉普拉斯方程只响应温差用绝对温度不会改变系数但发生第三类边界条件时T∞ 的基准必须与边界温度统一。5.2 法向方向决定结果正负外法向完全由边界顶点排列方向决定。按逆时针排列时单元法向量是 (dy/L, -dx/L)按顺时针排列则全局反号。一个非常实用的检查是任取一个单元中点沿法向向外移动一个小距离看它是否落在域外再沿法向向内移动看是否进入域内。对凹域、带孔洞的几何这个检查尤其必要因为局部法向可能不服从整体直觉。5.3 收敛性检查与条件数边界元加密网格时H 和 G 的非对角元素会逐渐趋近连续分布矩阵可能越来越病态。不要把 rcond 当成唯一标准最可靠的做法是单元数从 40 到 80 到 160 各跑一遍观察某个内部点温度是否收敛。常单元的理论收敛阶不高看到三位有效数字稳定就可以认为网格够用了再加密主要是验证不是验证精度。5.4 用优化工具箱做边界参数标定工程里经常出现这种情况外边界一部分热流未知但手上有内部几个测点温度这样的问题是典型参数反演。此时 BEM 的优势在于对每个候选边界条件内部点温度只需要重新解一个 n 阶线性系统计算量很小非常适合配合 MATLAB 优化工具箱做标定。% 反演未知边界热流 qHidden resid (qHidden) norm(Tmeas - runBEM(geometry, BC, qHidden)); qFit lsqnonlin(resid, q0, lb, ub, opts);runBEM 就是把边界元求解、内部点恢复封装成一个函数。用 lsqnonlin 而不是 fminsearch 的理由是残差向量可以同时包含多个内部测点能利用边界元对内部温度的高精度优势而且求解速度足够支撑迭代多次。参数标定前先把 5.1 的量纲和 5.2 的法向检查跑过一遍否则拟合出的热流很可能是一个符号相反的数值。本文还有配套的精品资源点击获取