简介一份基于边界元方法BEM的二维热场MATLAB计算程序面向涉及热传导数值模拟的工程师、科研人员及相关专业学生。程序将二维热传导问题转化为边界积分方程通过格林函数离散化边界配合线性代数方程求解与温度场可视化适用于复杂几何形状或非均匀边界条件的热分析场景。资源压缩包共3个文件全部为m脚本主程序、边界元核心函数与后处理模块相互分离整体仅3KB结构简洁便于阅读和二次开发适合用作边界元初学者理解算法流程的入门样例。已有480人学习下载。通过这份代码使用者可快速掌握边界元法求解热场的关键步骤并针对实际工况调整边界条件或几何参数节省自主编程与调试时间结合描述中的方向也可进一步扩展为工程热管理设计与理论验证的实用参考。1. 二维热场为什么要选边界元而不是有限元做热分析的第一反应通常是有限元。但二维稳态热场有个特殊性质控制方程是拉普拉斯方程基本解是解析可写的于是整个问题可以等价地改写成边界上的积分方程。真实未知量从区域内的温度场退化成边界上的温度和热流离散维度少一维。网格只画边界一圈算完之后再用积分公式反推任意内部点温度这就是边界元法BEM。边界元最大的优势在无限域或半无限域问题里体现得最明显。散热器外周的热量扩散到无穷远处有限元必须在远处截断并处理截断边界反射边界元的基本解本身就满足无穷远条件费半天劲画的大区域网格可以直接扔掉。对二维热场问题边界元通常能把有限元需要几千上万个单元的模型压缩到几百个边界单元代价是系数矩阵从稀疏变成稠密。适合用它的人是做热设计验证、埋地管道、电缆载流量、电子器件散热模型这类以“看边界温度/热流分布、算内部几个关键点温度”为主的工程问题。2. 二维热场边界元的数学底子基本解、边界积分方程与符号习惯2.1 稳态热传导方程与二维基本解不含内热源的二维稳态温度场满足拉普拉斯方程k(∂²T/∂x² ∂²T/∂y²) 0均匀介质里 k 可以约掉剩下的问题是纯几何的。边界元把偏微分方程转成积分方程的钥匙是基本解Green 函数二维拉普拉斯算子的基本解写作φ -1/(2π) · ln r其中 r 是源点与场点之间的距离坐标形式是 r √((x-x₀)² (y-y₀)²)。注意这里 ln r 前的负号是人为约定的很多教材写成 1/(2π)·ln(1/r)两者等价只是后面的边界积分每一项都会跟着变符号。这是边界元程序最容易翻车的地方之一建议整个代码统一采用 φ -ln r / (2π)后文所有矩阵组装都以这个符号约定为准。这个基本解本身对应的是“二维空间内单位点热源产生的温度场”它是圆的。边界元正是靠它把任意形状边界上的温度和热流联系起来。2.2 从格林第二恒等式到边界积分方程把拉普拉斯方程和基本解做加权余量经过格林第二恒等式变换可以得到如下边界积分方程c(ξ)·T(ξ) ∫Γ T(y)·∂φ/∂n(y) dΓ ∫Γ φ(y)·∂T/∂n(y) dΓ其中 ξ 是场点y 是边界上的积分点n 是边界外法线。系数 c(ξ) 由场点位置决定场点在计算域内部时 c1在光滑边界上时 c1/2在域外时 c0。这个式子把二维区域的偏微分方程严格转化成了沿边界曲线 Γ 的一维积分这就是“降维”的数学来源。等式左边包含边界温度 T 和基本解法向导数 ∂φ/∂n 的乘积右边包含边界法向热流 ∂T/∂n 和基本解的乘积。实际应用时边界上每个点要么给定温度第一类边界条件、要么给定法向热流第二类边界条件不可能同时给两者。边界积分方程妙就妙在通过联立所有边界点可以把未知的那一侧反解出来。2.3 边界热流的方向约定工程上热流密度 q 通常定义为 q -k·∂T/∂n即沿外法线方向为正表示流出计算域。边界元公式里直接用的是 ∂T/∂n它和 q 相差一个负号。下面的矩阵系统统一用 ∂T/∂n 作为求解量后处理如果需要物理热流记得翻转符号。提示在验证程序时先检查符号而不是先怀疑网格。一个最简单的检查方法是均匀温度场 Tconst此时 ∂T/∂n0边界积分方程应退化为 c·T ∫∂φ/∂n·T dΓ 0可以据此逐项核对 H 矩阵。3. 二维热场边界元离散化常数单元、配点与矩阵组装3.1 为什么常数单元是入门首选边界元里最常用的单元是常数单元、线性单元和二次单元。常数单元把每条直线段内的温度 T 和法向导数 ∂T/∂n 都视为常数取单元中点作为配点collocation point。此时场点落在边界上任意单元中点处时c(ξ)1/2 恒定不需要处理角和边的不连续问题代码写起来最干净。线性单元在相邻单元交点上共享温度自由度精度更高但在角点处法向导数本身可能不连续需要额外处理。对二维热场做初步仿真、验证方法、快速估计热流分布常数单元是最好的选择。它的收敛阶在内部点上接近二阶边界热流是一阶对工程估算足够。对于要求高精度的场合通常做法是先跑通常数单元再根据同样的配点框架换成线性单元。3.2 边界离散与法线朝向规则以圆环域为例计算域是 R1 r R2 的环形区域。离散分两步先在圆周上按角度等分生成单元再取每个单元的中点为配点。边界元有一条铁律计算域必须在边界走向的左侧法线指向计算域外部。圆环有两条边界。外圆边界 R2 按逆时针方向划分法线指向圆外内圆边界 R1 按顺时针方向划分法线指向圆心即指向孔洞内部这是计算域外部。下表总结了两种单元的法线与流向关系边界位置几何走向外法线方向计算域侧内圆 R1顺时针指向孔心法线反向外圆 R2逆时针指向远离圆心法线同向这个表中只要有一个法线方向反了最后的接出来的温度场会在某个边界附近出现明显的漏热或吸热假象收敛性也完全破坏。3.3 矩阵 H 与 G 的数值组装把边界离散成 N 个常数单元后对每个配点 i 写积分方程得到 N×N 线性系统H·T G·q其中 q 表示 ∂T/∂n 的边界节点值。矩阵元素定义为H(i,j)在单元 j 上对 ∂φ/∂n 做积分当 i≠j 时用高斯积分ij 时取 H(i,i)0.5G(i,j)在单元 j 上对 φ 做积分当 i≠j 时用高斯积分ij 时使用奇异积分解析值对于长度 L 的常数单元自作用奇异积分的解析结果是G(i,i) L/(2π) · (1 - ln(L/2))这里的 L 是第 i 个单元的长度。这个解析公式在边界元里几乎必写它避免了在高斯积分时被零距离除掉。如果不处理自作用项H 矩阵对角线缺 0.5、G 矩阵对角线缺奇异项线性系统解出来全是错的而且错误不会随网格加密消失。高斯积分用于 i≠j 的远场项时2 点高斯对这个二维热场问题已经足够。原因是常数单元上的积分核在非奇异情况下很光滑2 点高斯可以精确积分三次以下多项式继续加密到 4 点或者 8 点收敛速度不会有实际改善只会增加组装时间。3.4 边界条件重整把未知量摆到左边原始系统 H·T G·q 中每个节点上 T 和 q 只有一个已知、一个未知。设第 j 个节点上如果给的是 T第一类边界则未知量是 q_j如果给的是 q第二类边界则未知量是 T_j。于是构造一个统一的代数系统 A·x b当 j 节点的 T 未知时A(i,j) 对应的列取 H(i,j)x 分量是 T_j当 j 节点的 q 未知时A(i,j) 对应的列取 -G(i,j)x 分量是 q_j已知项全部移动到右侧 b这类重组在 MATLAB 里用逻辑索引即可实现不需要物理重排矩阵行。组装完成后求解 x A\b再把 x 里解出的值放回 T 和 q 向量。4. matlab实现二维热场边界元环域温度场的最小可运行程序4.1 生成圆环边界单元内圆半径 R1 和外圆半径 R2 之间夹着的区域就是计算域。下面这个函数生成所有边界单元的配点坐标、外法线方向、单元长度和内外圈标记function [xm, ym, nx, ny, len, isInner] ringMesh(R1, R2, N1, N2) th_i linspace(0, 2*pi, N11); th_i th_i(1:end-1) pi/N1; % 内圈单元中点角度 th_o linspace(0, 2*pi, N21); th_o th_o(1:end-1) pi/N2; % 外圈单元中点角度 xm [R1*cos(th_i), R2*cos(th_o)]; ym [R1*sin(th_i), R2*sin(th_o)]; nx [-cos(th_i), cos(th_o)]; % 内圈法线指向孔心 ny [-sin(th_i), sin(th_o)]; len [R1*2*pi/N1*ones(1,N1), R2*2*pi/N2*ones(1,N2)]; isInner [true(1,N1), false(1,N2)]; end这里每个单元用一个点中点代表,所以没有显式存储端点坐标。为单元数 N1, N2 分别控制内外圈离散密度如果板与板间距较大而内外圈半径差较大时N1和N2可以不同这个参数在收敛性分析中会反复修改。内圈法线取 -cos、-sin 就是指向孔心外圈取 cos、sin 指离圆心整个边界的外法线系统是闭合的。4.2 组装 H 矩阵与 G 矩阵function [H, G] assembleBEM(xm, ym, len, nx, ny) N length(xm); H zeros(N, N); G zeros(N, N); gp [-0.577350269189626, 0.577350269189626]; gw [1, 1]; for i 1:N xi xm(i); yi ym(i); for j 1:N if i j H(i,j) 0.5; G(i,j) len(i)/(2*pi) * (1 - log(len(i)/2)); else for g 1:2 t 0.5 * len(j) * gp(g); xj xm(j) - t*nx(j); % 从配点沿法线反向走到单元端点附近 yj ym(j) - t*ny(j); rx xj - xi; ry yj - yi; r2 rx^2 ry^2; r sqrt(r2); gradrn (rx*nx(j) ry*ny(j)) / r2; % d(ln r)/dn phi -log(r) / (2*pi); dphidn -gradrn / (2*pi); H(i,j) H(i,j) 0.5*len(j)*gw(g)*dphidn; G(i,j) G(i,j) 0.5*len(j)*gw(g)*phi; end end end end end这段代码有三个关键点。第一自作用项 H(i,i)0.5 来自常数单元配点在光滑边界的几何关系不需要积分G(i,i) 用解析式直接写。第二远场积分时被积点坐标用xm(j) - t*nx(j)表示意思是沿法线反向偏移半个单元长度再乘以高斯积分点的比例系数。这块的几何含义是单元中点向两侧延伸到单元端点正好覆盖整个单元长度。第三二维基本解的法向导数展开为-gradrn/(2*pi)其中 gradrn 是 ln r 沿法线的方向导数符号来自 φ-ln r/(2π) 这一约定和前面章节保持一致。4.3 组装代数方程并求解边界未知量R1 0.5; R2 1.0; N1 48; N2 64; [xm, ym, nx, ny, len, isInner] ringMesh(R1, R2, N1, N2); N length(xm); [H, G] assembleBEM(xm, ym, len, nx, ny); Tbc zeros(N, 1); Tbc(~isInner) 20; % 外圈温度 20 度 Tbc(isInner) 100; % 内圈温度 100 度 % 圆环内外都给定的是温度未知量全部是 q q G \ (H * Tbc);这一步直接用左除G \ (H * Tbc)求解因为两条边界都是第一类边界条件系数矩阵就是 G。如果有第二类边界条件混在里面需要按照 3.4 节的重组方式把未知量排列到左边。求解完成后 q 里存的是每个边界单元的 ∂T/∂n 值想换算成热流要乘 -k。4.4 计算内部任意点温度并画出二维热场nxq 120; nyq 120; xx linspace(-R2, R2, nxq); yy linspace(-R2, R2, nyq); [X, Y] meshgrid(xx, yy); inside (X.^2 Y.^2 R1^2 eps) (X.^2 Y.^2 R2^2 - eps); Tp NaN(size(X)); for ii 1:nxq for jj 1:nyq if inside(ii,jj) xp X(ii,jj); yp Y(ii,jj); sumT 0; for k 1:N dx xm(k) - xp; dy ym(k) - yp; r2 dx^2 dy^2; r sqrt(r2); gradrn (dx*nx(k) dy*ny(k)) / r2; phi -log(r) / (2*pi); dphidn -gradrn / (2*pi); sumT sumT G_int(k)*q(k) - H_int(k)*Tbc(k); end Tp(ii,jj) sumT; end end end等温线图画法很不讲究contourf(X, Y, Tp, 20); colorbar;就能直接出图。内部点温度积分所用的公式和组装 H、G 时完全相同只是配点在域内时 c1不再有 0.5 的自作用项。常见做法是另写一个反演函数先对每个内部点循环全部边界单元累加G_int·q - H_int·T。这段代码慢但胜在直观。4.5 完整主脚本与参数快速对照上面四个段落合起来就是完整的二维热场边界元程序。把边界单元数、半径和温度值改成实际问题参数即可直接运行。各参数的作用与调试关注点如下表参数作用调试关注点N1 / N2内外圈边界离散密度加密时看 q 的收敛趋势R1 / R2计算域几何范围检查是否满足 R1 R2Tbc已知边界温度确认 isInner 方向和温度一一对应nxq / nyq内部绘图网格密度只影响后处理不影响边界解5. 二维热场边界元的验证套路与最常见错误二维热场边界元程序写完后最有效的验证不是直接拿复杂工程模型去算而是用解析解做收敛性对照。圆环问题恰好有精确解。径向稳态温度分布满足T(r) T1 (T2 - T1) · ln(r/R1) / ln(R2/R1)取 T1100、T220、R10.5、R21.0对任意边界单元数 N内部点数值解应该与上式几乎重合。实际验证时这样测固定一个内部点 r0.75在 N16、32、64、128 四组网格下分别求该点温度与解析解做差。常数单元的局部误差应该随边界单元长度 h 线性到二次之间衰减画 log-log 图时斜率保持在 1 到 2 之间就说明矩阵组装无误。最常见的错误有一个固定套路法线方向。形容一下症状如果内圈法线方向写反计算出的 q 符号全部翻转温度场会出现内圈附近温度极端、等温线不对称的假象。可以加一个热流守恒检查稳态无热源问题中边界净热流应为零net_flux sum(q .* len); fprintf(边界净热流 %e\n, net_flux);圆环内外圈如果都是等温边界理论上 net_flux 为零实际计算时由于离散误差通常是一个 1e-12 量级的小数。如果这个值是 1e-1 量级基本可以断定法线或符号约定出了问题。还有一个容易被忽略的细节G 矩阵自作用项的解析式依赖于单元长度。内圈外圈单元长度不同每个单元的自作用项都不同不能用一个全局值代替。用常数单元时H 矩阵对角线恒为 0.5 与单元无关但一旦改用线性单元或二次单元对角线处的基本解法向导数积分不再是 0.5需要重新推导角点系数。这个差别是很多人在把常数单元程序扩展成高阶单元时卡住的地方。进阶应用中如果边界条件里有纯热流边界最后得到的 q 在边界上可能震荡这是配点法的固有现象尤其是角点附近。对一个正方形计算域加上两个相邻边绝热时角点处场解不唯一常数单元会在角点产生小的伪热流。缓解办法是把角点处的绝热边界拆成两个不同单元并忽略角点配点或者改用线性单元并在角点做双节点处理。先用圆环验证过符号和组装逻辑再往带角点的几何扩展排查范围就小得多。本文还有配套的精品资源点击获取 二维热场边界元法MATLAB实现：从基本解到环域温度场
简介一份基于边界元方法BEM的二维热场MATLAB计算程序面向涉及热传导数值模拟的工程师、科研人员及相关专业学生。程序将二维热传导问题转化为边界积分方程通过格林函数离散化边界配合线性代数方程求解与温度场可视化适用于复杂几何形状或非均匀边界条件的热分析场景。资源压缩包共3个文件全部为m脚本主程序、边界元核心函数与后处理模块相互分离整体仅3KB结构简洁便于阅读和二次开发适合用作边界元初学者理解算法流程的入门样例。已有480人学习下载。通过这份代码使用者可快速掌握边界元法求解热场的关键步骤并针对实际工况调整边界条件或几何参数节省自主编程与调试时间结合描述中的方向也可进一步扩展为工程热管理设计与理论验证的实用参考。1. 二维热场为什么要选边界元而不是有限元做热分析的第一反应通常是有限元。但二维稳态热场有个特殊性质控制方程是拉普拉斯方程基本解是解析可写的于是整个问题可以等价地改写成边界上的积分方程。真实未知量从区域内的温度场退化成边界上的温度和热流离散维度少一维。网格只画边界一圈算完之后再用积分公式反推任意内部点温度这就是边界元法BEM。边界元最大的优势在无限域或半无限域问题里体现得最明显。散热器外周的热量扩散到无穷远处有限元必须在远处截断并处理截断边界反射边界元的基本解本身就满足无穷远条件费半天劲画的大区域网格可以直接扔掉。对二维热场问题边界元通常能把有限元需要几千上万个单元的模型压缩到几百个边界单元代价是系数矩阵从稀疏变成稠密。适合用它的人是做热设计验证、埋地管道、电缆载流量、电子器件散热模型这类以“看边界温度/热流分布、算内部几个关键点温度”为主的工程问题。2. 二维热场边界元的数学底子基本解、边界积分方程与符号习惯2.1 稳态热传导方程与二维基本解不含内热源的二维稳态温度场满足拉普拉斯方程k(∂²T/∂x² ∂²T/∂y²) 0均匀介质里 k 可以约掉剩下的问题是纯几何的。边界元把偏微分方程转成积分方程的钥匙是基本解Green 函数二维拉普拉斯算子的基本解写作φ -1/(2π) · ln r其中 r 是源点与场点之间的距离坐标形式是 r √((x-x₀)² (y-y₀)²)。注意这里 ln r 前的负号是人为约定的很多教材写成 1/(2π)·ln(1/r)两者等价只是后面的边界积分每一项都会跟着变符号。这是边界元程序最容易翻车的地方之一建议整个代码统一采用 φ -ln r / (2π)后文所有矩阵组装都以这个符号约定为准。这个基本解本身对应的是“二维空间内单位点热源产生的温度场”它是圆的。边界元正是靠它把任意形状边界上的温度和热流联系起来。2.2 从格林第二恒等式到边界积分方程把拉普拉斯方程和基本解做加权余量经过格林第二恒等式变换可以得到如下边界积分方程c(ξ)·T(ξ) ∫Γ T(y)·∂φ/∂n(y) dΓ ∫Γ φ(y)·∂T/∂n(y) dΓ其中 ξ 是场点y 是边界上的积分点n 是边界外法线。系数 c(ξ) 由场点位置决定场点在计算域内部时 c1在光滑边界上时 c1/2在域外时 c0。这个式子把二维区域的偏微分方程严格转化成了沿边界曲线 Γ 的一维积分这就是“降维”的数学来源。等式左边包含边界温度 T 和基本解法向导数 ∂φ/∂n 的乘积右边包含边界法向热流 ∂T/∂n 和基本解的乘积。实际应用时边界上每个点要么给定温度第一类边界条件、要么给定法向热流第二类边界条件不可能同时给两者。边界积分方程妙就妙在通过联立所有边界点可以把未知的那一侧反解出来。2.3 边界热流的方向约定工程上热流密度 q 通常定义为 q -k·∂T/∂n即沿外法线方向为正表示流出计算域。边界元公式里直接用的是 ∂T/∂n它和 q 相差一个负号。下面的矩阵系统统一用 ∂T/∂n 作为求解量后处理如果需要物理热流记得翻转符号。提示在验证程序时先检查符号而不是先怀疑网格。一个最简单的检查方法是均匀温度场 Tconst此时 ∂T/∂n0边界积分方程应退化为 c·T ∫∂φ/∂n·T dΓ 0可以据此逐项核对 H 矩阵。3. 二维热场边界元离散化常数单元、配点与矩阵组装3.1 为什么常数单元是入门首选边界元里最常用的单元是常数单元、线性单元和二次单元。常数单元把每条直线段内的温度 T 和法向导数 ∂T/∂n 都视为常数取单元中点作为配点collocation point。此时场点落在边界上任意单元中点处时c(ξ)1/2 恒定不需要处理角和边的不连续问题代码写起来最干净。线性单元在相邻单元交点上共享温度自由度精度更高但在角点处法向导数本身可能不连续需要额外处理。对二维热场做初步仿真、验证方法、快速估计热流分布常数单元是最好的选择。它的收敛阶在内部点上接近二阶边界热流是一阶对工程估算足够。对于要求高精度的场合通常做法是先跑通常数单元再根据同样的配点框架换成线性单元。3.2 边界离散与法线朝向规则以圆环域为例计算域是 R1 r R2 的环形区域。离散分两步先在圆周上按角度等分生成单元再取每个单元的中点为配点。边界元有一条铁律计算域必须在边界走向的左侧法线指向计算域外部。圆环有两条边界。外圆边界 R2 按逆时针方向划分法线指向圆外内圆边界 R1 按顺时针方向划分法线指向圆心即指向孔洞内部这是计算域外部。下表总结了两种单元的法线与流向关系边界位置几何走向外法线方向计算域侧内圆 R1顺时针指向孔心法线反向外圆 R2逆时针指向远离圆心法线同向这个表中只要有一个法线方向反了最后的接出来的温度场会在某个边界附近出现明显的漏热或吸热假象收敛性也完全破坏。3.3 矩阵 H 与 G 的数值组装把边界离散成 N 个常数单元后对每个配点 i 写积分方程得到 N×N 线性系统H·T G·q其中 q 表示 ∂T/∂n 的边界节点值。矩阵元素定义为H(i,j)在单元 j 上对 ∂φ/∂n 做积分当 i≠j 时用高斯积分ij 时取 H(i,i)0.5G(i,j)在单元 j 上对 φ 做积分当 i≠j 时用高斯积分ij 时使用奇异积分解析值对于长度 L 的常数单元自作用奇异积分的解析结果是G(i,i) L/(2π) · (1 - ln(L/2))这里的 L 是第 i 个单元的长度。这个解析公式在边界元里几乎必写它避免了在高斯积分时被零距离除掉。如果不处理自作用项H 矩阵对角线缺 0.5、G 矩阵对角线缺奇异项线性系统解出来全是错的而且错误不会随网格加密消失。高斯积分用于 i≠j 的远场项时2 点高斯对这个二维热场问题已经足够。原因是常数单元上的积分核在非奇异情况下很光滑2 点高斯可以精确积分三次以下多项式继续加密到 4 点或者 8 点收敛速度不会有实际改善只会增加组装时间。3.4 边界条件重整把未知量摆到左边原始系统 H·T G·q 中每个节点上 T 和 q 只有一个已知、一个未知。设第 j 个节点上如果给的是 T第一类边界则未知量是 q_j如果给的是 q第二类边界则未知量是 T_j。于是构造一个统一的代数系统 A·x b当 j 节点的 T 未知时A(i,j) 对应的列取 H(i,j)x 分量是 T_j当 j 节点的 q 未知时A(i,j) 对应的列取 -G(i,j)x 分量是 q_j已知项全部移动到右侧 b这类重组在 MATLAB 里用逻辑索引即可实现不需要物理重排矩阵行。组装完成后求解 x A\b再把 x 里解出的值放回 T 和 q 向量。4. matlab实现二维热场边界元环域温度场的最小可运行程序4.1 生成圆环边界单元内圆半径 R1 和外圆半径 R2 之间夹着的区域就是计算域。下面这个函数生成所有边界单元的配点坐标、外法线方向、单元长度和内外圈标记function [xm, ym, nx, ny, len, isInner] ringMesh(R1, R2, N1, N2) th_i linspace(0, 2*pi, N11); th_i th_i(1:end-1) pi/N1; % 内圈单元中点角度 th_o linspace(0, 2*pi, N21); th_o th_o(1:end-1) pi/N2; % 外圈单元中点角度 xm [R1*cos(th_i), R2*cos(th_o)]; ym [R1*sin(th_i), R2*sin(th_o)]; nx [-cos(th_i), cos(th_o)]; % 内圈法线指向孔心 ny [-sin(th_i), sin(th_o)]; len [R1*2*pi/N1*ones(1,N1), R2*2*pi/N2*ones(1,N2)]; isInner [true(1,N1), false(1,N2)]; end这里每个单元用一个点中点代表,所以没有显式存储端点坐标。为单元数 N1, N2 分别控制内外圈离散密度如果板与板间距较大而内外圈半径差较大时N1和N2可以不同这个参数在收敛性分析中会反复修改。内圈法线取 -cos、-sin 就是指向孔心外圈取 cos、sin 指离圆心整个边界的外法线系统是闭合的。4.2 组装 H 矩阵与 G 矩阵function [H, G] assembleBEM(xm, ym, len, nx, ny) N length(xm); H zeros(N, N); G zeros(N, N); gp [-0.577350269189626, 0.577350269189626]; gw [1, 1]; for i 1:N xi xm(i); yi ym(i); for j 1:N if i j H(i,j) 0.5; G(i,j) len(i)/(2*pi) * (1 - log(len(i)/2)); else for g 1:2 t 0.5 * len(j) * gp(g); xj xm(j) - t*nx(j); % 从配点沿法线反向走到单元端点附近 yj ym(j) - t*ny(j); rx xj - xi; ry yj - yi; r2 rx^2 ry^2; r sqrt(r2); gradrn (rx*nx(j) ry*ny(j)) / r2; % d(ln r)/dn phi -log(r) / (2*pi); dphidn -gradrn / (2*pi); H(i,j) H(i,j) 0.5*len(j)*gw(g)*dphidn; G(i,j) G(i,j) 0.5*len(j)*gw(g)*phi; end end end end end这段代码有三个关键点。第一自作用项 H(i,i)0.5 来自常数单元配点在光滑边界的几何关系不需要积分G(i,i) 用解析式直接写。第二远场积分时被积点坐标用xm(j) - t*nx(j)表示意思是沿法线反向偏移半个单元长度再乘以高斯积分点的比例系数。这块的几何含义是单元中点向两侧延伸到单元端点正好覆盖整个单元长度。第三二维基本解的法向导数展开为-gradrn/(2*pi)其中 gradrn 是 ln r 沿法线的方向导数符号来自 φ-ln r/(2π) 这一约定和前面章节保持一致。4.3 组装代数方程并求解边界未知量R1 0.5; R2 1.0; N1 48; N2 64; [xm, ym, nx, ny, len, isInner] ringMesh(R1, R2, N1, N2); N length(xm); [H, G] assembleBEM(xm, ym, len, nx, ny); Tbc zeros(N, 1); Tbc(~isInner) 20; % 外圈温度 20 度 Tbc(isInner) 100; % 内圈温度 100 度 % 圆环内外都给定的是温度未知量全部是 q q G \ (H * Tbc);这一步直接用左除G \ (H * Tbc)求解因为两条边界都是第一类边界条件系数矩阵就是 G。如果有第二类边界条件混在里面需要按照 3.4 节的重组方式把未知量排列到左边。求解完成后 q 里存的是每个边界单元的 ∂T/∂n 值想换算成热流要乘 -k。4.4 计算内部任意点温度并画出二维热场nxq 120; nyq 120; xx linspace(-R2, R2, nxq); yy linspace(-R2, R2, nyq); [X, Y] meshgrid(xx, yy); inside (X.^2 Y.^2 R1^2 eps) (X.^2 Y.^2 R2^2 - eps); Tp NaN(size(X)); for ii 1:nxq for jj 1:nyq if inside(ii,jj) xp X(ii,jj); yp Y(ii,jj); sumT 0; for k 1:N dx xm(k) - xp; dy ym(k) - yp; r2 dx^2 dy^2; r sqrt(r2); gradrn (dx*nx(k) dy*ny(k)) / r2; phi -log(r) / (2*pi); dphidn -gradrn / (2*pi); sumT sumT G_int(k)*q(k) - H_int(k)*Tbc(k); end Tp(ii,jj) sumT; end end end等温线图画法很不讲究contourf(X, Y, Tp, 20); colorbar;就能直接出图。内部点温度积分所用的公式和组装 H、G 时完全相同只是配点在域内时 c1不再有 0.5 的自作用项。常见做法是另写一个反演函数先对每个内部点循环全部边界单元累加G_int·q - H_int·T。这段代码慢但胜在直观。4.5 完整主脚本与参数快速对照上面四个段落合起来就是完整的二维热场边界元程序。把边界单元数、半径和温度值改成实际问题参数即可直接运行。各参数的作用与调试关注点如下表参数作用调试关注点N1 / N2内外圈边界离散密度加密时看 q 的收敛趋势R1 / R2计算域几何范围检查是否满足 R1 R2Tbc已知边界温度确认 isInner 方向和温度一一对应nxq / nyq内部绘图网格密度只影响后处理不影响边界解5. 二维热场边界元的验证套路与最常见错误二维热场边界元程序写完后最有效的验证不是直接拿复杂工程模型去算而是用解析解做收敛性对照。圆环问题恰好有精确解。径向稳态温度分布满足T(r) T1 (T2 - T1) · ln(r/R1) / ln(R2/R1)取 T1100、T220、R10.5、R21.0对任意边界单元数 N内部点数值解应该与上式几乎重合。实际验证时这样测固定一个内部点 r0.75在 N16、32、64、128 四组网格下分别求该点温度与解析解做差。常数单元的局部误差应该随边界单元长度 h 线性到二次之间衰减画 log-log 图时斜率保持在 1 到 2 之间就说明矩阵组装无误。最常见的错误有一个固定套路法线方向。形容一下症状如果内圈法线方向写反计算出的 q 符号全部翻转温度场会出现内圈附近温度极端、等温线不对称的假象。可以加一个热流守恒检查稳态无热源问题中边界净热流应为零net_flux sum(q .* len); fprintf(边界净热流 %e\n, net_flux);圆环内外圈如果都是等温边界理论上 net_flux 为零实际计算时由于离散误差通常是一个 1e-12 量级的小数。如果这个值是 1e-1 量级基本可以断定法线或符号约定出了问题。还有一个容易被忽略的细节G 矩阵自作用项的解析式依赖于单元长度。内圈外圈单元长度不同每个单元的自作用项都不同不能用一个全局值代替。用常数单元时H 矩阵对角线恒为 0.5 与单元无关但一旦改用线性单元或二次单元对角线处的基本解法向导数积分不再是 0.5需要重新推导角点系数。这个差别是很多人在把常数单元程序扩展成高阶单元时卡住的地方。进阶应用中如果边界条件里有纯热流边界最后得到的 q 在边界上可能震荡这是配点法的固有现象尤其是角点附近。对一个正方形计算域加上两个相邻边绝热时角点处场解不唯一常数单元会在角点产生小的伪热流。缓解办法是把角点处的绝热边界拆成两个不同单元并忽略角点配点或者改用线性单元并在角点做双节点处理。先用圆环验证过符号和组装逻辑再往带角点的几何扩展排查范围就小得多。本文还有配套的精品资源点击获取