COMSOL Multiphysics® 中的固体瞬态加热建模介绍

COMSOL Multiphysics^®软件经常被用来模拟固体的瞬态加热。瞬态加热模型很容易建立和求解，但它们在求解时也不是没有困难。例如，对瞬态加热结果的插值甚至会使高级 COMSOL^®用户感到困惑。在这篇博客中，我们将探讨一个简单的瞬态加热问题的模型，并利用它来深入了解这些细微差别。

一个简单的瞬态加热问题

图1显示了本文所讨论主题的建模场景。在这个场景中，将一个空间上均匀分布的热载荷施加在一个具有均匀初始温度的圆柱体材料顶面的圆形区域内。最开始载荷很高，但在一段时间后会逐渐下降。除了施加热载荷外，还添加了一个边界条件来模拟整个顶面的热辐射，它使零件重新冷却。假设材料属性(热导率、密度和比热)和表面辐射率在预期温度范围内保持不变，并且假设没有其他作用的物理场。我们的建模目标是用它来计算圆柱体材料内随时间变化的温度分布。

在 COMSOL 案例库中的硅晶片激光加热教程模型中，有一个类似的建模场景，但请记住，本文讨论的内容适用于任何涉及瞬态加热的情况。

图1.顶面有一个热源的圆柱体材料几何模型。

尽管我们很想通过绘制图1中所示的精确几何结构开始建立模型，但我们可以从一个更简单的模型开始。在图1中，可以看到几何体和载荷是围绕中心线轴向对称的，所以我们可以合理地推断，解也将是轴向对称的。因此，我们可以将模型简化为二维轴对称建模平面。(点击此处，了解如何使用对称性来减小模型尺寸。)

在中间的圆形区域内，热通量是均匀的。最简单的建模方法是通过在二维域的边界上引入一个点来修改几何形状。这个点将边界划分为受热和未受热的部分。在几何形状上增加这个点，可以确保所产生的网格与热通量的变化完全一致。考虑到这些，我们可以创建一个等效于三维模型的二维轴对称计算模型(图2)。

图2.相当于三维模型的二维轴对称模型。显示的是默认网格。

此外，我们还考虑了施加的热通量大小的瞬时变化的情况；在 t=0.25s 时，它的值变得较低。载荷的这种阶梯式变化应该通过使用事件接口来解决，如 COMSOL 知识库中关于求解包含时变载荷阶跃变化的模型一文所述。简单来说，事件接口会准确地告诉求解器载荷的变化什么时候发生，求解器将相应地调整时间步长。我们可能也想知道求解器采取的时间步长，这可以通过修改求解器的设置，按求解器的步长输出结果，然后就可以绘制零件顶部中心点的温度，如图3所示。

图3.某一点的温度随时间变化的曲线图，各点显示了求解器在载荷突然变化的附近采取的步长较短。

接下来，我们用不同的求解器相对容差值重新运行该模型，并在图中进行比较(图4)。这类图表明，像预期的那样，随着公差变小，解迅速向同一个值收敛。

图4. 用不同的相对容差求解出的随时间变化的某一点的温度图。

另一个可以计算的量是进入该域的总能量。我们可以对通过边界的总热通量的表达式ht.nteflux进行积分，使用timeint()算子对时间进行积分，得到总能量。积分的结果在下面的表格中列出，用于增加时间步长的相对容差。(提示：你可以在 COMSOL知识库中了解更多关于计算空间和时间积分的信息，在这篇关于如何计算质量守恒和能量平衡的博客中了解更多关于计算能量平衡的信息)。

求解器相对容差	通入域的通量的时间积分(J)
1e-2	32.495
1e-3	32.469
1e-4	32.463

从数据中我们可以观察到，进入系统的总能量实际上几乎与时间步长容差无关。乍一看，这似乎是对我们模型的一个奇妙的验证。然而，需要指出的是，我们在这里观察到的是有限元法(FEM)的基本数学特性。简单说，就是总能量总是会很好地平衡。这并不意味着模型中没有错误，错误只是出现在不同的地方……接下来，我们就去寻找错误。

错误：很容易产生，但很难定义

我们应该在这里暂停一下，来非常谨慎地处理上文中提到的一个词，即错误这个提示，它在建模和仿真的世界中经常被使用，但没有固定出现的场合。在本节的后面部分，我们将对各种建模案例中可能出现的不同错误进行一些详细描述。(如果你想直接跳到与模型中的错误有关的部分，请点击这里)。

输入错误

输入错误，顾名思义，是指模型输入中的错误，如材料属性输入不正确或几何形状绘制错误。最有危害的一个输入错误就是遗漏错误，例如忘记添加一个边界条件。输入错误与输入中的不确定性是不同的，例如，当不知道确切的材料属性时，就会出现不确定性。前者输入错误只能通过仔细检查来解决，而后者输入中的不确定性可以通过 COMSOL 软件的不确定性量化模块来解决。对于我们的例子，我们确定没有输入错误或不确定因素。

几何体的离散化错误

当通过有限元网格离散几何体时，特别是在对非平面边界进行网格划分时，会产生一个几何体的离散错误。这些错误随着网格细化程度的增加而减少，并且可以在不实际求解有限元模型的情况下进行计算。本文示例中的二维轴对称建模域没有弯曲的边界，不必担心这种类型的错误。

解的离散化错误

解离散错误是由于有限元基函数不能完全代表真实的解场及其在此域内的导数。它从根本上存在于有限元方法中。这种与几何离散误差有内在联系的误差总是存在的，对于任何良好的有限元问题来说，它总是随着网格的细化而减少。

时间步长误差

了解时域模型中的误差传播是相当复杂的。这篇博客，我们只要说在任何一个时间步长中引入的或已经存在的任何误差都会向前传播就足够了，但对于文中讨论的扩散类问题，它们会逐渐衰减。这种类型的误差总是存在的，而且这些误差的大小是由瞬态求解器容差和网格控制的。

插值错误

还有一种类型的错误是比较定性的，那就是插值错误。这些错误发生在对结果的意义和产生方式没有准确理解的情况下。其中最著名的是尖角处的奇异性，这种情况经常出现在结构力学以及电磁场建模中。当存在输入错误时，插值错误尤其经常出现。因此，如果你对你的结果有任何不确定的地方，一定要回去仔细检查(甚至三番五次检查！)你模型的所有输入。

上面列举的错误清单并不完整。例如，我们还可以谈一谈由于线性系统求解器的有限精度算术、非线性系统求解器和数值积分误差而产生的数值误差。然而，这些以及其他类型的误差，基本上规模都小得多的。

有了上述的这组定义，现在准备回到我们的模型了。

追踪空间和时间中的错误

到目前为止，我们已经观察了模型中某一点的解，并观察到随着我们完善瞬态求解器的相对容差，解似乎收敛得很好，所以我们应该已经理解了收紧瞬态求解器的相对容差将减少时间步长误差的想法。现在，我们来看看空间温度分布。我们将从沿中心线的温度开始，对于最宽松的容差 1e-2，看初始时间的解以及求解器采取的第一个时间步长，如下图所示。

图5.初始时间和第一个时间步长的中心线温度图。

从初始值图中，我们可以看到，沿中心线的温度与规定的初始温度不一致——有些地方甚至低于初始值。这是由于 COMSOL Multiphysics 使用了所谓的一致初始化，即调整初始时间的解场，使其与初始时间的边界条件和初始值一致。一致初始化包括采取一个额外的非常小的人工时间步长，我们可以认为是在零时刻发生的。一致性初始化可以在求解器的初始时间步长设置中关闭，也可以在显式事件和隐式事件功能中关闭，但是这样做的时候应该谨慎。在常见的的多物理场模型中，尤其是涉及到流体流动的模型，默认情况下要让它处于可能更稳健的启用状态，所以我们在这里将讨论这种情况。

在这种情况下，考虑一致初始化的方式是，调整温度场使之与施加的载荷和边界条件相一致。由于施加的热载荷最初是不为零的，温度场的梯度与热通量成正比，最初也必须不为零。我们还需要考虑，这个场是用有限元基函数离散的。沿着中心线，这些基函数是多项式，但多项式不可能完全匹配真实的解；因此，在一致初始化步骤之后，我们最终得到的是一个会略微超过或低于预期结果的解。从第一个时间步长的解中，我们还可以看到，零件远端的温度已经在上升，这是意料之外的。虽然这些与预期的变化幅度非常小，但我们还是希望将它们降到最低。

在通过修改模型来减少这些解的离散化误差之前，让我们应用一点物理上的直觉来解决这个问题。在模拟时间开始时，沿中心线的温度分布将与涉及模拟通过一维板块的传热建模方案相当类似。对于这种类型的建模方案，已经存在一个分析解，这在许多关于传热分析的书籍中经常会谈到(实际上，这个例子被用作我桌上的一本教科书Fundamentals of Heat and Mass Transfer的封面插图)。

为了简洁起见，我们将跳过分析解法直接陈述结果。当对表面零件施加热量时，表面的温度将开始上升，最终内部区域也会变暖。请注意，离边界较远的点需要更多的时间来加热。板块内的温度是不会均匀变化的。在内部更远的点，与靠近表面的点相比，需要更长的时间才能使温度开始变化。值得注意的是，由于传热方程的扩散性质，空间温度变化将随着时间的推移而趋于平稳。有了这种认识，让我们再回到我们的模型，看看如何改进它。

简单来说，为了使这个解的离散化误差最小化，我们需要在场变化剧烈的位置划分更精细的网格。根据我们的经验（或者分析解，如果我们想查的话），我们知道，在非常接近表面和边界的法线方向上，场的变化非常大，但在内部则变得更加平滑。这正是需要边界层网格划分的情况，如图6所示，它在边界的法线上创建薄的单元。

图6.通过沿着晶圆顶部的一个边界添加一个边界层来修改网格划分序列。

现在，我们可以重新运行模拟，并绘制出初始时间和下一个时间步长的解。

图7.使用边界层网格时，初始和第一个时间步长的中心线温度图。

在图7中，我们可以观察到，在初始值时，温度的下调在空间上更为局部。事实证明，使用更精细的网格也会导致随时间变化的求解器采取更小的时间步长。因此，通过这种细化的网格，我们减少了空间离散和时间步长误差。

我们还可以看一下沿建模域顶部边界的结果，代表暴露表面的温度分布。图8中显示的是使用 1e-2 的容差绘制的初始时间和第一个时间步长。在这些图中，我们可以观察到空间中相当剧烈的震荡场。这是空间离散化的一个表现。请记住，热载荷沿着径向轴经历了大小的阶跃变化，我们在这里观察到的有点类似于吉布斯现象。

图8.使用边界层网格绘制的在初始值和第一个时间步长中沿顶面的温度图。

解与之前类似，但现在我们必须在过渡的位置细化网格。对于这个问题，可以对划线点应用更加精细的大小设置，从而得到如下所示的网格。

图9.网格设置和网格的图示，在划定热载荷分布的点上应用较小的网格大小。

从图10的温度结果中，我们可以看到，现在解中的振荡已经减少了，在空间和时间上的传播也没有那么多。即使求解器的相对容差为 1e-2，求解的结果也已经有了很大的改善。

图10.在细化网格后，使用0.01的相对容差，在初始值和第一个时间步沿顶面的温度要准确得多。

我们可以使用更多的网格和求解器容差细化来继续这个练习。但通过我们目前所做的细化，已经开始看到误差迅速减少–由于瞬态传热方程的扩散性质，即使是仍然存在的误差也会在空间和时间上被平滑掉。实际上，我们应该花同样多的精力来研究模型输入中的不确定性的影响。

还有什么可以发挥作用？

在这个例子中，施加在边界上的热载荷并没有及时移动，所以划分边界的方法是合理的。如果热载荷分布要移动，那么受热面的整个网格就需要更加精细。(你可以在这篇博客中了解在 COMSOL^®中对移动载荷和约束进行建模的 3 种方法)。

在这篇博客的前面部分，我们假设材料属性随温度的变化而保持不变，并且不依赖于任何其他物理场。这是一个重大的简化，因为所有的材料属性都会随温度变化。材料甚至可以经历相变，如融化。模拟相变可以用几种不同的方法，包括表观热容法，它是使用高度非线性的比热来说明相变的潜热。我们也可以很容易地预见到这是一个多物理场问题，如涉及热固化的方程式，甚至是材料非线性电磁加热问题。在这种情况下，我们不仅需要监测温度场的收敛性，还需要监测所有其他正在求解的场变量，甚至可能是它们的时间和空间导数。这些情况可能都需要在建模域的所有地方采用非常精细的网格，所以从本例这种简单情况中得到的经验可能并不适用。然而，即使在对更复杂的模型进行网格划分和求解时，记住最简单的情况总是好的(即使它只是作为一个概念性的起点)。

此外，我们应该强调的是，这篇文章只是关于固体材料的时变传热。如果有一个移动的流体，控制方程将发生重大变化；流体流动模型的网格划分是一个单独的、相对更复杂的话题。对于波动型问题，网格的选择和求解器的设置就变得相当简单了。

结束语

在这篇博客中，我们复习了一个典型的传热建模问题。我们注意到，在空间和时间上，在载荷突然变化的情况下，解会出现某些错误。读者现在应该对这些错误类型有所了解，并知道它们是有限元方法的固有结果，就像所有的数值方法一样，只是对现实的一种近似。尽管这些误差看起来很大，但由于瞬态传热方程的扩散性，它们的大小在空间和时间上都会衰减。我们已经表明，网格细化将减少空间离散误差，同时隐含着减少时间步长误差的效果。最后，我们讨论了如何通过求解器相对容差细化来进一步减少时间步长误差。

还值得做一个更简短的总结：如果你主要对一个相对长的时间后的解感兴趣，使用相当粗的网格和默认的求解器相对容差是完全可以接受的。另一方面，如果你对短时和小范围的温度变化感兴趣，那么必须研究求解器相对容差和网格细化。

理解了这些，我们就可以避免犯解释错误。这将使我们能够自信和轻松地从简单的模型中建立更复杂的模型。

下一步

点击下面的按钮，你将进入 COMSOL 案例库，可以尝试自己模拟这篇博客中提到的示例模型。