简介：围绕COMSOL求解一维固体热传导方程，这份docx文档面向已有一定PDE与传热基础、希望用数值软件复现温度场演化过程的高校学生和工程技术人员。内容以“固体传热”物理场和瞬态研究为主线，涉及几何区间设置、导热系数、密度、恒压热容等用户定义参数、左右边界温度约束、网格与时间步range(0,0.01,1)配置，以及计算后导出data.txt并用Matlab作后处理的完整流程，有助于理解热传导方程ρcp∂T/∂t=k∂²T/∂x²的离散求解与结果可视化。压缩包共1个docx文件，约45KB，篇幅紧凑，适合作为COMSOL入门练习和课堂演示的参考材料。已有3685人学习，说明其在传热仿真初学者中有一定参考价值。读者可据此梳理PDE建模步骤、边界条件设置与数据导出后处理思路，减少自行摸索成本。
1. 一根 14 米长的杆：这个一维热传导仿真要解决什么
用 COMSOL 解热传导方程，卡住人的往往不是方程本身，而是几何长度为什么填 14。一个 0 到 1 秒、导热系数、密度、恒压热容全取 1 的一维固体传热瞬态模型，把计算域拉到 14 米，看着像手滑多敲了个 4，其实是刻意的：热扩散长度 √(αt) 在 t=1 时只有 1 米，域长取 14 米意味着右边界的影响在 1 秒内还没传回来，模型等价于半无限大介质，数值解才有资格和 erfc 解析解逐点对比。这套流程——建物理场、设材料参数、划网格、跑瞬态研究、导出 data.txt、丢进 MATLAB 后处理——属于典型的验证型仿真：目的不是算出某个工程温度，而是证明你能用 COMSOL 的传热接口或 PDE 接口把一个答案已知的问题算准，并且把误差说清楚。适合正在做传热数值模拟课设、需要交一份从建模到误差分析完整报告的人。
2. 物理场与参数落地：把热传导方程落到 COMSOL 的三个输入框
2.1 方程里每个符号对应哪个设置
先把要解的方程摆出来，一维、无内热源、无流动的固体传热：
ρ·cp·∂T/∂t = k·∂²T/∂x²
其中 ρ 是密度、cp 是恒压热容、k 是导热系数、T 是温度、t 是时间、x 是空间坐标。把两边同除以 ρcp，得到 ∂T/∂t = α·∂²T/∂x²，α = k/(ρcp) 就是热扩散系数。这个 α 决定了整个算例的时间尺度：温度扰动向前推进的特征距离是 √(αt)，t = 1 秒时大约是 1 米，这就是后面域长取值、网格密度、输出时间步全部要对齐的那个数。
COMSOL 的固体传热接口内部写的是更一般的形式 ρcp(∂T/∂t + u·∇T) + ∇·q = Q，q = −k∇T。一维、u = 0、Q = 0 之后自动退化成上面那个式子，所以材料节点里那三个输入框——导热系数、密度、恒压热容——正好一一对应方程里的 k、ρ、cp，不多不少。这也是选固体传热接口而不是从零搭 PDE 的核心理由：物理量和边界条件的名字是现成的，报告里写起来也和教材对得上。
把三个参数全部设成用户定义并且都等于 1，本质上是做一次无量纲化。令 x* = x/L、t* = αt/L²，方程变成 ∂T/∂t* = ∂²T/∂x*²，一个参数都不剩。取 α = 1 等于把长度单位和时间单位绑死在一起，这样算出来的曲线可以直接和标准解对比，不用在单位换算上反复折腾。答辩时被问“参数取 1 有什么物理意义”，标准答法是：这是一次无量纲验证，考察的是离散方案和边界处理，不是某种具体材料的物性；真要换成铝或钢，只需要把三个参数替换成实际值，流程一个字都不用改。
2.2 固体传热接口还是系数型 PDE
摘要里提到用 PDE 模块模拟热传导，这条路确实通，但两条路的操作量和出错概率差别不小。下面的对照表是我在两种建法之间来回切过之后整理的，课程设计报告里可以直接拿来说明选型理由：
| 对比项 | 固体传热接口 | 系数型 PDE |
|---|---|---|
| 方程形式 | ρcp∂T/∂t + ∇·q = 0，q = −k∇T | ea∂²u/∂t² + da∂u/∂t + ∇·(−c∇u − αu + γ) + β·∇u + au = f |
| 待求变量 | 温度 T，单位 K | 通用变量 u，无量纲 |
| 参数入口 | 导热系数、密度、恒压热容三个具名框 | c、da、a、f、α、β、γ 一堆系数框 |
| 左端 1000 K | 右键“温度”边界，填 1000 | Dirichlet 边界，填 1000 |
| 初始值 | 具名“初始值”节点，单位 K | 具名“初始值”节点，填数值 |
| 后处理 | 温度、热流、梯度都有内置变量 | 只有 u 和它的导数 |
系数型 PDE 的映射关系是 c = k = 1、da = ρcp = 1、a = 0、f = 0，其余系数全零，这样两个接口解的是同一个方程。选 PDE 的好处是能把弱形式、系数含义讲得更透，适合做方法学部分的展开；选固体传热的优势是边界条件、单位、后处理变量都是物理量语义，出图时横轴直接是米、纵轴直接是 K，省掉一层换算。我的习惯是主模型用固体传热跑结果，另建一个 PDE 版本做交叉验证，两者温度场最大偏差应该在 1e-6 K 量级，如果对不上，多半是某个系数忘了清零。
2.3 初始值与温度约束：阶跃边界的正确写法
初始值这一步最容易翻车。固体传热接口新建之后，初始值节点的默认温度是 293.15 K，不是 300 K。直接跑下去，全场从 293.15 K 起步，看上去只差 6.85 K，但在误差评估里这 6.85 K 会被算成模型误差，后面无论怎么加密网格都消不掉。所以第一步就是把初始值手动改成 300 K，改完回头看一眼节点上显示的数字，确认不是默认值。
边界条件是两个温度约束。右键物理场添加“温度”节点，选择左边界，填 1000 K；再加一个，选右边界，填 300 K。注意这是硬约束，左端在 t = 0 那一刻从 300 K 直接跳到 1000 K，数学上是个阶跃，解析解在 t → 0⁺ 时刻在 x = 0 处有奇性，温度梯度趋于无穷。这一点直接决定了后面两件事：第一，t = 0 那一列输出数据不能拿去和解析解比；第二，网格必须在左端加密，否则阶跃附近的数值振荡会非常难看。
关于“针对边界 1 固体 1”和“针对边界 2 固体 1”这两组重复的设置，如果几何里存在两个固体域节点（比如为了分区控制网格把 0–14 拆成两段），材料参数和初始值要在每个固体节点下各设一遍，只改第一处、漏掉第二处，跑出来的结果会是“半截模型”：左半段参数正确，右半段还在用默认的空气物性，温度曲线在中段出现一个不合理的拐点。检查方法是点开每个固体节点，确认三个参数后面都显示“用户定义”并且值为 1。
3. 几何、网格与时间步：三个互相牵制的离散参数
3.1 域长 14 的定量依据与截断误差估算
有限域和半无限域的差别，可以精确写成一项误差。初始 300 K、左端阶跃到 1000 K 的半无限大介质解是：
T(x,t) = 300 + 700·erfc( x / (2√(αt)) )
其中 erfc 是余误差函数。如果计算域只取到 L，右边界又被强制固定在 300 K，那么真实解在 x = L 处本来会升高的那部分温度被强行压回去，误差量级就是 700·erfc(L / (2√(αt)))。把 α = 1、t = 1 代进去：L = 1 时 erfc(0.5) ≈ 0.4795，右端温度本该升到 600 K 以上，误差接近 50%，这个模型根本不能用解析解对比；L = 7 时 erfc(3.5) ≈ 7.4e-7；L = 14 时 erfc(7) 已经在 1e-23 量级，比双精度浮点数的分辨率还小。
反推更直观：想让截断误差压到 1e-6 以下，需要 L / (2√(αt_max)) ≥ erfcinv(1e-6) ≈ 3.46，也就是 L ≥ 6.9 米。取 14 米是留了一倍余量，顺便让整条曲线在 1 秒内完全看不到右边界的存在。这个推导值得原样写进报告的“几何与计算域”一节，它把一个看起来随意的数字变成了有依据的设计选择。实际工程里做长杆淬火、半无限地基温度响应这类题目，也是同样的思路：域长取到特征扩散长度的 7 倍以上，右边界类型（固定温度还是绝热）就不再影响结果。
3.2 网格：把单元堆在 x = 0 附近
一维几何用“间隔”建，左端点 0、右端点 14，单位保持默认的米。网格这一步，均匀划分 140 个单元看上去 h = 0.1 挺合理，但 t = 0.01 时热扩散长度只有 0.1 米，整个温度变化都挤在 x ∈ [0, 0.3] 这一小段里，h = 0.1 意味着这一段只有两三个单元，梯度根本描述不出来。典型表现是左端出现超过 1000 K 的超调，或者温度曲线在头部鼓出一个包。
我的做法是改用分布型网格，在“分布”节点里设置单元数和单元比，让左端单元尺寸明显小于右端。粗略的经验是首单元尺寸取 √(α·t_min) 的十分之一，t_min 是最早需要看的输出时刻：t = 0.01 时 √(0.01) = 0.1，首单元尺寸 0.01 米比较稳，单元总数 200 左右。想更省事就直接给单元数 200 加一个单元比 10，然后靠网格无关性检查兜底。
下面这组量级是我在一维线性单元上反复调参时常见的结果，可以当作判断网格够不够的参照：
| 方案 | 单元数 | 首单元尺寸 | t=0.01 最大超调 | t=0.5 相对误差 |
|---|---|---|---|---|
| A | 50（均匀） | 0.28 | 约 +40 K | 3e-2 |
| B | 200（分布） | 0.01 | 约 +6 K | 2e-3 |
| C | 800（分布） | 0.002 | 约 +0.5 K | 3e-4 |
超调这一列特别值得盯：物理上温度不可能超过 1000 K，只要出现 T > 1000.5 K 的节点，就说明网格还不到位。顺带说一句，单元阶次从线性换到二次，超调会明显下降，代价是同样的单元数下自由度翻倍，一维问题这点开销完全可以接受。
3.3 时间步：range(0,0.01,1) 是输出时间，不是求解步长
瞬态研究里把时间设成range(0,0.01,1)，含义是 0 到 1 秒、间隔 0.01 秒，一共 101 个输出时刻。很多人到这里会以为求解器就按 0.01 秒一步往前走，于是在“时间步进”里改成用户控制、填 0.01，结果第一个步长就跨过了阶跃后的剧烈变化段，曲线头部出现明显的数值振荡，甚至温度跌到 300 K 以下。
正确的理解是：这 101 个点是“结果保存时刻”，求解器内部默认用自由步进，BDF 格式会自己把步长压到很小，尤其在 t 刚过 0 的那几个毫秒。所以常规配置就是输出时间写 range 函数、求解器保持“自由”、勾选在输出时间保存结果，中间的计算步长交给自适应。如果因为报告需要必须固定步长，那么步长要小到能跟上温度前沿的推进速度，起步阶段 1e-4 秒量级比较安全，之后可以再用一个分段的时间列表放宽。
还有一条从显式差分课上传下来的经验在这里不适用：α·Δt/h² ≤ 0.5 那个稳定性条件针对的是显式格式，COMSOL 默认的隐式 BDF 是无条件稳定的，限制步长的是精度而不是稳定性。答辩时如果有人拿这个条件来问网格和时间步的关系，能把这个区别讲清楚，基本就说明你真跑过而不是抄的。
4. 从 data.txt 到误差曲线：MATLAB 后处理的完整链路
4.1 导出设置：列的含义必须先定死
计算跑完之后走“结果 → 导出 → 数据”，数据集选一维绘图组或者直接用网格，表达式填 T，时间步选全部，格式选表格。导出之前一定要在界面上确认两件事：一是导出的数据集是不是整个域（选成“边界”就只剩两个端点，选成“点”就只有几个采样点）；二是时间步列表是不是完整的 101 个，漏选的话后面 MATLAB 里的列数就对不上。
data.txt 的常见结构是第一列 x 坐标、后面每一列对应一个输出时刻的温度，但也可能因为格式选项不同变成“每行一组 (x, t, T)”的长表。两种都能处理，关键是先打开文件看一眼前几行，别想当然。文件开头通常有几行以 % 开头的注释，记录导出时间和表达式，读的时候要跳过。
另一个容易忽略的点：x 列不一定单调递增。多域建模时导出顺序按单元排列，可能出现分段拼接，排序没做的话画出来的曲线会来回折返，看着像数值发散，其实是数据顺序问题。
4.2 读文件与解析解对比
读取和自检放在同一段里，列数对不上就直接报错停下，比画出一张错图再去倒查省时间：
% 读取 COMSOL 导出的一维瞬态温度数据 % 约定：第 1 列为 x 坐标（米），第 2 列起依次为 t = 0, 0.01, ..., 1.0 的温度 raw = readmatrix('data.txt', 'CommentStyle', '%'); x = raw(:,1); % 节点坐标 T = raw(:,2:end); % 每列一个时刻的温度场 t = 0:0.01:1; % 与导出列一一对应 % 排序：防止多域导出时 x 非单调导致曲线折返 [x, idx] = sort(x); T = T(idx, :); % 列数自检：时间点数必须等于数据列数 assert(size(T,2) == numel(t), ... '时间步数与数据列数不一致，检查导出时的时间步选择'); fprintf('读入 %d 个节点，%d 个时刻\n', numel(x), numel(t));readmatrix的CommentStyle参数负责跳过 % 开头的注释行；sort那两行是给多域模型兜底的；assert看着多余，但它拦住的是后面所有分析的地基错误——列错位一列，误差曲线会整体平移，很容易被误判成求解器精度问题。
接下来算解析解和误差。归一化用温差 700 K，不要用 300 K 或者绝对温度，否则相对误差会被大数掩盖：
alpha = 1; % 热扩散系数 k/(rho*cp)，三个参数均取 1 T0 = 300; % 初始温度 Tw = 1000; % 左端温度约束 dT = Tw - T0; % 归一化基准 % 取 t = 0.1 / 0.5 / 1.0 三个时刻（跳过 t=0，该时刻解析解有奇性） idx = round([0.1 0.5 1.0] / 0.01) + 1; err = zeros(1, numel(idx)); for k = 1:numel(idx) tk = t(idx(k)); T_exact = T0 + dT * erfc(x / (2*sqrt(alpha*tk))); % 半无限大介质阶跃解 err(k) = max(abs(T(:,idx(k)) - T_exact)) / dT; % 最大偏差 / 温差 end disp(table(t(idx)', err', 'VariableNames', {'t','relMaxErr'}));erfc是 MATLAB 内置的余误差函数，和前面推导用的是同一个；round(.../0.01)+1把物理时刻换算成列号；误差用最大绝对偏差除以 700，比 L2 范数更直观，报告里也好写“全场最大偏差不超过温差的千分之一”。
还有一点必须提前想明白：热流 q = −k∂T/∂x 在 x = 0 处随时间按 1/√t 发散，t 越小越离谱。如果用热流做验证指标，只能取 t ≥ 0.05 之后的时刻，早期数据纯属自找麻烦。
4.3 出图：温度剖面、时空云图和误差趋势
报告里的图不用多，三张能讲清问题就够。第一张是多时刻温度剖面：横轴 x、纵轴 T，把 t = 0.05、0.1、0.3、0.5、1.0 五条曲线画在一起，实线数值解、圆圈解析解，两条线重合程度一眼可见。第二张是时空云图，用imagesc(x, t, T')或surf，横轴位置、纵轴时间、颜色温度，能直观看到温度前沿向右推进的速度，也能看出右边界有没有被“污染”。第三张是误差随时间的曲线，横轴 t、纵轴最大偏差，正常情况下这条曲线在 t 很小时偏高（阶跃奇性导致），之后迅速回落并稳定在网格精度对应的水平上。
画图时把axis的范围固定住，别用自动缩放，否则不同网格方案生成的图比例不一致，放在一起没法比。颜色条统一用同一个范围 300–1000 K，云图之间的差异才看得出来。
5. 避坑与排查：这套流程里最容易翻车的五个地方
5.1 几何单位被改成厘米，整条曲线都对不上
现象：温度在 0.02 秒左右就“撞”到右边界，剖面图从斜坡变成拱形，和 erfc 曲线完全不像。
原因：几何节点里单位下拉框被改成了 cm，右端点填 14 实际只有 0.14 米，比 t = 1 时的扩散长度还短，半无限假设直接失效。
解决：几何、材料参数、输出时间三处的单位必须在建模开始前统一成 SI。改完之后重新检查一遍几何节点上显示的“14 m”，确认不是“14 cm”。这个坑我在帮人看模型时遇到过不止一次，而且因为数值不报错，非常隐蔽。
5.2 初始值没改，全场从 293.15 K 起步
现象：数值解和解析解整体差一个接近 7 K 的常数偏移，误差曲线水平一条线，加密网格、减小时间步都没用。
原因：固体传热接口的初始值默认是 293.15 K，建模时看到“初始值 1”这个节点名就跳过去了，没点开改数值。
解决：点开初始值节点，把温度改成 300 K，改完看节点摘要里的数字。凡是出现“误差不随离散加密下降”的情况，第一反应就该是初始条件或边界条件设错了，而不是网格不够。
5.3 把输出时间当成求解步长，阶跃处出现振荡
现象：t 很小时的几条曲线在左端附近上下抖动，甚至出现低于 300 K 或高于 1000 K 的点。
原因：在时间步进里选了用户控制并填 0.01，求解器第一步就跨过了阶跃后的剧烈变化段；或者网格在左端太粗，陡梯度下 Galerkin 格式产生超调。
解决：先恢复自由步进，让 BDF 自适应压步长；再把左端网格加密到首单元 0.01 米以下。两步做完还有超调，就把单元阶次从线性换成二次。
5.4 导出数据的列错位，图看着像发散
现象：画出来的温度剖面在两段之间跳变，或者时间轴对不上。
原因：导出时数据集选成了边界，只剩两个端点；或者 x 列因为多域拼接不单调，直接plot(x, T)会来回折返。
解决：导出前确认数据集覆盖整个域、时间步全选；读入后先sort再画图；用列数自检断言把错误拦在入口。
5.5 只跑一套网格就下结论
现象：报告里写了“数值解与解析解吻合良好”，但拿不出网格无关性的证据，答辩时被追问一句就答不上来。
原因：一维问题跑得快，容易只算一遍就收工，缺少收敛性验证这一环。
解决：至少跑三档网格（比如 50、200、800 单元）和两档时间步，把最大偏差列成表。判据是网格加密一倍后误差变化小于 5%，时间步减半后结果变化小于 2%，两条都满足才叫“网格无关”。
6. 验收与答辩：用三条曲线和一张表证明解是对的
验收这件事，别等到报告写完才想起来。我的做法是模型一跑通就先搭一张验收表，每跑一版填一行，最后直接把表贴进报告，比事后补写“结果分析”省力得多：
| 检查项 | 合格判据 | 在哪里看 |
|---|---|---|
| 边界一致性 | 左端恒为 1000 K，右端恒为 300 K，全程无漂移 | 一维绘图组的两个端点值 |
| 无超调 | 全场 300 K ≤ T ≤ 1000 K（容差 0.5 K） | 表面图的最大值最小值 |
| 与解析解偏差 | t ≥ 0.05 时最大偏差 / 700 < 1e-3 | MATLAB 误差表 |
| 网格无关性 | 单元加密一倍，最大偏差变化 < 5% | 两档网格结果对比 |
| 时间步无关性 | 输出步长减半，同一时刻温度变化 < 2% | 两档时间步结果对比 |
| 能量合理性 | 热流在 x = 0 处随时间单调下降 | 热流导出曲线 |
这张表里最容易被跳过、又最容易被问的是后两项。答辩时常见的三个问题是“为什么域长取 14”“参数取 1 有什么意义”“网格够不够”，对应的答法其实前面都推导过：第一问用 erfc(7) 的量级回答，第二问用无量纲化回答，第三问直接把三档网格的误差表甩出来——有数就有底气，没数就只能说“看起来挺光滑的”，那句话在答辩场上没有任何分量。
如果时间还宽裕，加一个交叉验证会更稳：用系数型 PDE 接口按 c = 1、da = 1、a = 0、f = 0 再建一遍，两个模型的温度场最大偏差应该在 1e-6 K 量级。这一手能让“我确实理解方程和接口的对应关系”这句话变成可展示的证据，而不是一句自述。
从那以后我每次做这类验证型仿真，都强制自己先写验收表再点计算，跑完第一版就把三档网格排上，绝不允许“只跑一遍就写结论”。这套习惯救过我一次——某次报告都快定稿了，网格无关性一测才发现 200 单元根本没收敛，改到 800 才压住超调，要是等到答辩现场被问出来，场面会很难看。这份从建模、参数设置、网格时间步到 MATLAB 后处理和验收表的完整流程，我整理成了一套可直接复现的文件，需要照着做课设或毕设的话，拿去当模板改参数就行。希望帮到你。
本文还有配套的精品资源，点击获取
