企业官网建设流程全解析

简介：一份基于边界元方法（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(ξ) 由场点位置决定：场点在计算域内部时 c=1，在光滑边界上时 c=1/2，在域外时 c=0。这个式子把二维区域的偏微分方程严格转化成了沿边界曲线 Γ 的一维积分，这就是“降维”的数学来源。

等式左边包含边界温度 T 和基本解法向导数 ∂φ/∂n 的乘积，右边包含边界法向热流 ∂T/∂n 和基本解的乘积。实际应用时，边界上每个点要么给定温度（第一类边界条件）、要么给定法向热流（第二类边界条件），不可能同时给两者。边界积分方程妙就妙在通过联立所有边界点，可以把未知的那一侧反解出来。

2.3 边界热流的方向约定

工程上热流密度 q 通常定义为 q = -k·∂T/∂n，即沿外法线方向为正表示流出计算域。边界元公式里直接用的是 ∂T/∂n，它和 q 相差一个负号。下面的矩阵系统统一用 ∂T/∂n 作为求解量，后处理如果需要物理热流，记得翻转符号。

提示：在验证程序时，先检查符号而不是先怀疑网格。一个最简单的检查方法是均匀温度场 T=const，此时 ∂T/∂n=0，边界积分方程应退化为 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 按顺时针方向划分，法线指向圆心（即指向孔洞内部，这是计算域外部）。下表总结了两种单元的法线与流向关系：

这个表中只要有一个法线方向反了，最后的接出来的温度场会在某个边界附近出现明显的漏热或吸热假象，收敛性也完全破坏。

3.3 矩阵 H 与 G 的数值组装

把边界离散成 N 个常数单元后，对每个配点 i 写积分方程，得到 N×N 线性系统：

H·T = G·q

其中 q 表示 ∂T/∂n 的边界节点值。矩阵元素定义为：

• H(i,j)：在单元 j 上对 ∂φ/∂n 做积分，当 i≠j 时用高斯积分，i=j 时取 H(i,i)=0.5
• G(i,j)：在单元 j 上对 φ 做积分，当 i≠j 时用高斯积分，i=j 时使用奇异积分解析值

对于长度 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, N1+1); th_i = th_i(1:end-1) + pi/N1; % 内圈单元中点角度 th_o = linspace(0, 2*pi, N2+1); 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 时完全相同，只是配点在域内时 c=1，不再有 0.5 的自作用项。常见做法是另写一个反演函数，先对每个内部点循环全部边界单元，累加G_int·q - H_int·T。这段代码慢，但胜在直观。

4.5 完整主脚本与参数快速对照

上面四个段落合起来就是完整的二维热场边界元程序。把边界单元数、半径和温度值改成实际问题参数即可直接运行。各参数的作用与调试关注点如下表：

5. 二维热场边界元的验证套路与最常见错误

二维热场边界元程序写完后，最有效的验证不是直接拿复杂工程模型去算，而是用解析解做收敛性对照。圆环问题恰好有精确解。

径向稳态温度分布满足：

T(r) = T1 + (T2 - T1) · ln(r/R1) / ln(R2/R1)

取 T1=100、T2=20、R1=0.5、R2=1.0，对任意边界单元数 N，内部点数值解应该与上式几乎重合。实际验证时这样测：固定一个内部点 r=0.75，在 N=16、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 在边界上可能震荡，这是配点法的固有现象，尤其是角点附近。对一个正方形计算域加上两个相邻边绝热时，角点处场解不唯一，常数单元会在角点产生小的伪热流。缓解办法是把角点处的绝热边界拆成两个不同单元并忽略角点配点，或者改用线性单元并在角点做双节点处理。先用圆环验证过符号和组装逻辑，再往带角点的几何扩展，排查范围就小得多。

本文还有配套的精品资源，点击获取