如何设置特殊边界条件

假设你正在模拟这样一种情况：载荷横跨不同网格单元和边界移动。如果你希望只对一部分几何边界或只在特定条件下施加边界条件，该怎么处理？在这篇博客中，我们将讨论如何利用 COMSOL Multiphysics 灵活处理这种特殊情况。

边界条件的分类

对偏微分方程进行数学处理时，会碰到三种边界条件：狄利克雷、诺伊曼和罗宾。狄利克雷条件用于指定要求解的变量；诺伊曼条件用于指定通量，也就是因变量的梯度；罗宾条件结合了前两类边界条件，用于指定变量及其梯度之间的关系。

下表列举了这三种边界条件在一些不同物理场中的应用示例，以及对应的物理解释。

在使用有限元方法时，这三类边界条件会对所求解问题的矩阵结构产生不同的影响。

诺伊曼条件

诺伊曼条件是“载荷”，位于方程组右侧。在 COMSOL Multiphysics 的方程视图中，这类边界条件显示为弱贡献。诺伊曼条件纯粹是方程组右侧附加的贡献，因此包含以下变量的任何函数：时间、坐标或参数值。

我们来看一个传热问题示例。在这个示例中，半径为R的圆形热源沿x轴方向以速度v移动。其强度呈抛物线分布，峰值为q_0。该载荷的数学描述为

很明显，移动载荷不可能有域边界，甚至不可能存在一个始终适合载荷分布的网格。

我们可以在该表达式中直接输入载荷分布。因为有两处会用到径向坐标变量r，所以将其定义成变量是一个不错的想法。该移动热源的完整输入如下图所示。

描述移动热源的参数。

描述从移动热源的当前中心触发的局部径向坐标的变量。

输入热通量。

下方动画显示的是采用上述数据进行瞬态仿真后的结果。假设该问题关于yz平面对称，那么载荷实际上被施加到一个半圆形移动热源上。

热源沿块移动时的温度分布动画。

狄利克雷条件

给定狄利克雷条件意味着指定了因变量，所以无须求解因变量。我们可以从问题中删除这类自由度方程。因此狄利克雷条件会改变刚度矩阵的结构。在 COMSOL Multiphysics 的方程视图中，这类条件将显示为约束。

假定要将移动点的温度指定为 450 K，这或许有点刻意，但这能表现出诺伊曼条件和狄利克雷条件之间的一个重要区别。假如要添加一个温度节点并输入一个类似的表达式（if(r < R,450[K],0)），这意味着没有被热源点覆盖的那部分边界的温度将被设定为绝对零度。

不过，我们的目标是在热源点之外停用狄利克雷条件。为此可以使用一个小窍门：如果将if(r < R,450[K],ht.Tvar)作为指定值输入，就能实现目标（如下方动画所示）。

含条件约束的狄利克雷条件设置。

指定温度的热源点沿块移动时的温度分布动画。

为了理解其运行原理，我们启用了方程视图，来查看狄利克雷条件的实现（本示例中为指定温度）：

启用方程视图。

方程视图中的温度节点 。

约束用公式表示为ht.T0-ht.Tvar，代表ht.T0-ht.Tvar = 0。第一项是作为输入的指定温度。第二项是转换成变量的温度自由度。这个约束使温度等于给定值，除非给定值恰好是ht.Tvar。如果发生这种情况，符号代数在分配过程中会将表达式简化为ht.Tvar-ht.Tvar，并进一步计算为零。当约束表达式为0时，表示没有约束。

弱约束

在 COMSOL Multiphysics 中，狄利克雷条件实际上有两种可能的实现方式：默认情况下为上文提到的逐点约束，还可以使用弱约束。在弱约束中，我们要添加方程，而不是移除方程。热通量需要被添加为额外的自由度（拉格朗日乘子），以强制指定温度值。

只要对上面的温度表达式稍作更改，就可以使用本质上相同的小窍门使弱约束含条件限制。要使用弱约束，必须先启用高级物理场接口选项。

启用高级物理场接口选项。

在物理场接口内对某个节点选中弱约束后，就会为拉格朗日乘子添加额外的自由度。这个示例中，该自由度名为T_lm。

如果应用了上文中的温度表达式，那么在停用狄利克雷条件的边界部分，额外的自由度不会得到任何刚度矩阵的贡献，因此刚度矩阵变为奇异矩阵。为了避免这种情形，我们将if(r < R,450[K],ht.Tvar)更改为if(r < R,450[K],ht.Tvar-T_lm*1e-2)。不同的模型和物理场，用于T_lm的乘子也不同，因此可能需要进行一些调整。

至于其中的工作原理，我们会像教科书中一样“作为练习留给读者”。

使用弱约束时，含条件限制的狄利克雷条件设置。

罗宾条件

罗宾条件通常对刚度矩阵和方程右侧都会造成影响。虽然刚度矩阵的结构不会受到影响，但会在现有位置上的值添加值。在方程视图中，罗宾条件同样显示为弱贡献。将这类条件转换为关于时间、空间和其他变量的函数与使用诺伊曼条件时的做法一致。

不过有趣的是，选择合适的值确实可以转换罗宾条件，使之近似为狄利克雷条件或诺伊曼条件。如果模拟过程中你希望在这两类边界条件之间切换，那么这一点十分重要。

要创建狄利克雷条件，需要对“刚度”指派一个高值，例如弹簧常数或传热系数。在数学术语中，这实质上是狄利克雷条件的罚实现。刚度越高，自由度的指定值就越精确。但这里需要注意：刚度过高会影响刚度矩阵的数值稳定性。而在传热问题中，要选择“高”的传热系数\alpha，起始值可以是

其中，\lambda表示导热系数，h表示特征单元尺寸。

用适当的材料属性（也就是固体力学中的杨氏模量）替换\lambda，就可以在其他物理场实现相同的计算。将因子设为 1000 只是一个建议，您也可以替换成 10⁴或 10⁵。

如果要使用上一个对流模拟示例中的移动温度为 450 K 的热点，则可以采用下图中的设置。表达式中使用了单元大小的内置变量h。

利用对流条件指定温度。

我刚开始使用有限元分析时，有时无法在结构力学的有限元程序中指定非零位移。这一限制是由于编程越来越复杂造成的。遇到这种情况时，最佳方法是添加预变形的刚性弹簧来使用罚函数法。但是刚度不能过大，因为那时还在采用单精度运算！

下面，我们来讨论近似为诺伊曼条件。我们希望热通量与表面温度无关。在这个传热示例中，罗宾条件表明向内热通量q为

其中，\alpha表示传热系数，T表示边界温度，T_{\textrm{ext}}表示外部温度。

所以，如果T_{\textrm{ext}}比表面预期温度高很多，那么q \approx \alpha T_{\textrm{ext}}。这时就是选择一个任意大的T_{\textrm{ext}}，然后计算一个合适的传热系数，如下图高亮处所示。

利用对流条件指定热通量。

在实际情况中，设计人员利用这个原理在一个实际机械组件引入了一个固定力：一根预应力的又长又软的弹簧。如果弹簧的预变形远大于与弹簧连接部件的位移，则弹簧的力几乎保持不变。

解决离散误差

当边界条件受制于布尔表达式（类似if(r < R,...）时，施加了此边界条件的区域边界很可能不会跟随网格单元的边，这将引入离散误差。

使用诺伊曼条件或罗宾条件时，要对每个有限元进行数值积分，并计算该单元中大量离散高斯点的函数值。如果网格单元的大小大于载荷的几何大小，那么施加有载荷的确切高斯点数会显著影响总载荷。由此，施加有载荷的小块区域上应该总是存在几个单元。

载荷位置稍微移动一点就可能改变所含积分点的数量。（在实际情况中，改变的积分点数会更多。）

因此，至少从离散点的意义上讲，要对网格节点施加狄利克雷条件。下图显示，在模拟指定温度为 450 K 的圆形点的移动时，某一时间点上的温度分布和热通量。在热点的前端，一个温度为 260 K 的暗影清晰可见。由于模拟的初始温度和环境温度均为 293 K，因此这一暗影不应出现。出现这一数值差异，是基于这样一个事实：并非每个单元上的所有节点都设定了狄利克雷条件。在狄利克雷条件的不连续性处，就会出现奇异点。这一话题我们在上一篇博客中探讨过。细化网格可以减少这类问题发生。

下图中的绿色箭头表示指定温度后产生热通量的节点。模型中的网格密度太小，使得近似半圆变得相当困难。

指定温度的半圆周围的温度分布和热通量。

边界条件中解的依赖性

要将解作为边界条件的输入有多种方法。这样做往往会引入非线性，但COMSOL Multiphysics 可以自动检测到非线性。

我们以一个梁为例。梁的下方有一个支撑，其作用是在梁发生一定挠曲后阻止梁的进一步移动。在梁接口的指定位移/旋转节点中，设置一个含条件限制的狄利克雷条件可以实现模拟。

具有挠曲、控制支撑和分布载荷的梁。

指定梁应在挠曲为 2 厘米时停止移动的设置。

分析表明预期的梁行为发生了。当载荷较小时，挠曲形状呈对称结构，而载荷较大时，梁上受额外支撑的点会停止移动。载荷达到最大时，梁上的曲率甚至会由正值变为负值。变形图可以显示这一变化，不过弯矩图中更加明显。

在支撑点处移动 2 厘米后停止的梁位移。

载荷不同时梁上的弯矩。

文中重点讨论的方法其实还相当不成熟，而且迭代解的收敛特性可能还不够好。更可靠的实现方式是在支撑点处使用高度非线性弹簧，由此反作用力就是关于位移的连续可微函数。这种做法实际上相当于在固体力学接口中实现罚接触。

在 COMSOL Multiphysics 中使用边界条件的结论和思考

COMSOL Multiphysics 允许您访问各种强大的原理，用于指定非标准边界条件。今天，我们只讨论了几个使用这类条件的案例。

如果你对分析含有移动载荷的模型感兴趣，请查看 COMSOL App 库中的移动载荷教程模型。

如果你在模拟过程中对于如何指定非标准边界条件有其他问题，请立即联系我们。