在 COMSOL Multiphysics® 软件的大多数功能中，都包含 弱约束 选项。这篇博客，我们将深入探讨什么是弱约束、为什么要使用弱约束，以及使用弱约束时需要特别考虑的事项。
章节内容
如何在 COMSOL Multiphysics® 中实施弱约束？
几乎在所有情况下，COMSOL Multiphysics® 中默认的约束类型都是逐点约束。逐点约束直接应用于自由度 (DOF)，通常是网格中的一组节点。
在模型开发器的大多数标准约束功能的设置中，都有一个 约束设置 栏。在这个栏，您可以在两个不同的约束之间选择。一些功能的约束列表中可能包含第三个选项，即 Nitsche 约束。
约束设置 栏示例。
不过，默认情况下该栏并不显示。弱约束属于 “高级”功能。您可以在 显示更多选项 对话框中选择 高级物理场选项 来启用该栏。
拉格朗日乘子
由于弱约束是基于 拉格朗日乘子，我们有必要先从该主题的一些基本概念讲起。
拉格朗日乘子的概念是由数学家约瑟夫•路易•拉格朗日（Joseph-Louis Lagrange）在变分法研究中提出。在我们的博客诞辰快乐，约瑟夫•路易•拉格朗日中，您可以了解更多关于他的生平和工作。
一般的带约束最小化问题可表述如下：在一组约束条件 g_i(\mathbf x) = 0 的约束下，求函数 f(\mathbf x) \;\text{,} \; \mathbf x \in \Omega 的最小值。接下来的问题是，如何以一般的方式满足这些约束条件。对于简单约束（例如线性约束），可以显式地反解约束表达式并将结果代入 f(\mathbf x)，从而减少未知数的数量。但这只是例外情况，而非普遍规律。
现在可以通过求解以下函数的最小值来解决问题：
式中，\lambda_i 是一组新的未知变量，即拉格朗日乘子。由于 g_i(\mathbf x) = 0 为真解，因此很明显，\mathcal{L} 的最小值等于受约束 f 的最小值。但为什么会这样呢？
要找到函数的极值，需要计算所有变量的偏导数，并将它们设为零。在这个示例中，
和
因此，通过对方程组求关于拉格朗日乘子的偏导就可以恢复约束方程，也就是说，约束条件仍然有效。但是，第一组方程更为精妙。要理解对于特定值，为什么它们能确保获得正确的最小值，需要更详细的论证。值得注意的是，新方程现在还涉及约束函数的梯度。
如需完整且详细的解释和举例说明，请参阅：拉格朗日乘子。
有限元约束处理
考虑一个简单的线性静态有限元问题。最终离散矩阵形式可以写为
式中，\mathbf{K} 是刚度矩阵，\mathbf {u} 是自由度向量，\mathbf {f} 是节点载荷向量。为了能够求解此方程组（使刚度矩阵非奇异），必须已知至少一定数量的自由度值。例如，对于传热问题，至少需要一个温度，而对于三维固体力学问题，至少需要知道 6 个位移自由度。最常见的方法是直接输入该自由度的数值。此外还有其他选项，如热量传递分析中的对流边界条件或结构力学中的弹簧条件。
实际应用中往往会施加远超稳定性所需的约束，例如在整个边界上应用约束。只要存在一些受约束的自由度，就可以将方程组按如下形式进行划分：
式中，下标 a 表示激活的自由度，下标 c 表示受约束的自由度。矢量 \mathbf {r} 包含受约束节点上未知的反作用力。这一术语来源于固体力学反作用力的概念。例如，在热量传递中，反作用力就是热通量。大多数情况下，约束节点上不施加任何载荷\mathbf f_{\mathrm{c}} ，因为它们对解没有影响。
\mathbf u_{\mathrm{c}} 的值是已知的，例如 \mathbf u_{\mathrm{c}}= \mathbf u_0 。在求解有限元方程时，最基本的处理方法是只求解激活的自由度。约束值的影响可以移到方程右侧，得到
求解简化后的方程组后，就可以使用以下方程在结果估计步骤计算反作用力
实际计算时，只需组装 \mathbf k_{\mathrm{aa}} 矩阵，其余矩阵乘法运算可以在单元层面更高效地完成。
这种表述反作用力的方法假定问题是线性的，因此矩阵是常数。一种更通用的表述方式是，反作用力是残差计算的结果。对于激活的自由度，收敛解的残差很小。在受约束节点，残差较大，可以解释为反作用力。
在 COMSOL Multiphysics® 中，基本上 后验 是可以做到的使用 reacf() 算子计算与节点相关的反作用力。请注意，得到的值是集中的节点值而非连续分布场。根据单元中使用的形函数，这种节点值的分布可能不是很直观。不过，其求和结果在数值舍入误差范围内是精确的。
这种方法的缺点是，如果需要在计算中使用反作用力，它们就无法使用。
注：上述描述基于一般的有限元公式，是对 COMSOL Multiphysics® 内部实际处理方式的极大简化。实际上约束处理要复杂得多。例如，可能涉及非线性约束，或是连接不同物理场变量的跨场约束。
引入拉格朗日乘子法
现在，继续讨论离散化形式。为了不那么抽象，我们举一个固体力学的例子。对于线性问题，可以证明势能为
将 W 与位移最小化，就能得到之前使用过的普通方程组。现在，将所有约束条件与各自对应的拉格朗日乘子相乘后加入势能表达式：
为清晰起见，自由度向量分为无约束 DOF ( \mathbf u_{\mathrm{a}}) 和有约束 DOF ( \mathbf u_{\mathrm{c}}) 。关键区别在于：新公式中将受约束自由度同样视为未知量。
为简单起见，假设受约束的自由度无外载作用。对各组变量求导后得到如下方程组：
利用第二行方程可以很容易地验证，拉格朗日乘子与上文介绍的节点反作用向量 \mathbf {r} 完全相同。最后一行方程简单地说明了\mathbf u_{\mathrm{c}} = \mathbf u_0 。
如果原始刚度矩阵是对称的（通常是这种情况），那么这个新方程组也是对称的。但此表述存在以下局限：
- 与只求解激活的自由度不同，这种计算方法需要求解激活的自由度、受约束的自由度和拉格朗日乘子。
- 矩阵对角线存在零元素，部分线性方程组解法无法处理
- 出于数值计算考虑，必须对方程组进行适当的比例缩放。原刚度矩阵元素的量级可能与单位值相差甚远。
尽管如此，该方法仍具显著优势：
- 反作用力值本身构成问题表述的一部分，不能仅视作结果量
- 使用拉格朗日乘子时，强非线性约束条件的收敛效果会更好。
- 在逐点约束中无法使用时间导数。如果需要，使用拉格朗日乘子法是唯一的选择。
弱约束
弱约束条件同样基于拉格朗日乘子概念。不过，在离散化之前的数学描述中已经考虑了约束条件。
有限元公式的基础是使用基本方程的弱形式。在固体力学中，这也被称为虚功原理。计算公式如下
式中， \sigma是应力张量，\mathbf {t} 是应变张量， 是牵引向量， \mathbf {u} 是位移场。\delta 符号表示变分算子。在 COMSOL Multiphysics® 中，它由 test() 和 var() 算子表示。符号约定与软件保持一致。
为简单起见，假设载荷（牵引力）只作用于边界。在更一般的情况下，也可能存在体积、边和点载荷。
为了使公式更加完整，需要在边界的某些部分设置位移的指定值，
在用选定的形函数近似位移时，该数学表达式将转化为有限元方程的离散化形式。
基于拉格朗日乘子的思想，现在可以通过添加一个额外项，将约束条件纳入弱表达式中：
式中, \lambda(\mathbf x) 是一个定义在约束边界上的拉格朗日乘子。假设 \mathbf u_0 只是一个指定值（与解无关），那么最后一项可以扩展为三项，即
当使用有限元法将弱形式表达式转换为离散方程时，与位移一样，\lambda(\mathbf x) 使用形函数在单元上近似。原则上，两个场的形函数可以相互独立地选择，但这样会失去刚度矩阵的对称性（甚至会使矩阵变得奇异）。因此，通常使用与位移相同的形函数。组装后的方程组如下所示
矩阵 \mathbf L 和 \mathbf L^T 源自上述弱表达式中的第三项和第四项，只要位移和拉格朗日乘子使用相同的形函数，它们就具有对称性。
补充说明：矩阵 \mathbf L 实际上与同一边界上单位面积质量分布的质量矩阵贡献相同。在 固体力学 接口中，可以使用 附加质量 节点给出这种质量分布。
为便于比较，我们重复上述对离散化系统应用拉格朗日乘子得到的方程组：
可以看出，结构类似，但单位矩阵被 \mathbf L 所取代，右侧现在包含 \tilde {\mathbf u}_0 ，它是 \mathbf u_0 的加权形式。
那么，这种修改会产生什么影响？指定自由度的节点值将不再与指定值完全一致。另一方面，节点之间的值在平均意义上将更接近给定函数 \mathbf u_{0}(\mathbf x)。
为了说明这一点，我们来看一个简单的传热示例，比较用标准逐点约束法和弱约束法得到的结果。
示例一：传热
使用二阶拉格朗日单元在单位正方形上以 2×2 网格求解一个二维传热问题。在最右边的边界上，温度设定为 sin(2*Y)。在相反的边界上，温度设为零。在其余边界上，没有任何边界条件，对应于无通量。温度选择指定正弦函数分布是为了避免二次形函数对其进行精确描述。
如果我们沿着右边界绘制温度曲线，两种方法之间的差异几乎重合，如下图所示。
在下图中，我们可以看到温度与其设定值之间的比较，这是一个更有趣的曲线图。
从上图可以看出，在节点点上，逐点约束与指定函数完全一致。而使用弱约束条件时，情况并非如此。
还存在一个或其他更好的解吗？第一个尝试是计算与给定函数相比的平均误差。结果如下表所示。
| 约束 | T-sin(2*Y) 的平均值 |
|---|---|
| 逐点约束 | 2.5*10-4 |
| 弱约束 | 3.6*10-7 |
事实证明，使用弱约束条件时，平均误差要小三个数量级。通过添加拉格朗日乘子场，指定值将在平均意义上尽可能地一致。为了证明这不是巧合，我们进行了一次参数扫描，将边的网格密度从 1 个单元改为 20 个单元。
需要强调的是，这些误差量级本身已经非常微小。对于任何合理细化的网格而言，上述对比主要具有理论研究价值。
在接下来的案例中，约束类型的选择将产生直观可见的影响。
示例二：杆模型
此例中，我们使用的是一根横截面为正方形的杆，它受到单轴拉伸。这是一个有两个域的装配体。一个域使用六面体单元网格划分，另一个域使用四面体单元网格划分。由于这是一个装配体，因此使用 连续性 条件将两个域连接起来。其结果是，网格在共同边界处 不兼容。
默认情况下，连续性是通过逐点约束来实现的，因此 目标 侧的每个节点都必须具有与 源 侧相应位置完全相同的位移。根据选择哪一侧作为目标，将得到不同的结果。
由于两个域交界处的形函数不一致，应力场中存在明显的局部扰动。随着与连接处距离的增加，扰动会迅速消失。COMSOL博客：“圣-维南原理的应用和解释 ”中详细讨论了这一现象。
如果我们现在改用弱约束，结果会好得多。在下图中，应力采用了另一种标度。误差大约小了一个数量级。
结论表明，弱约束法能有效缓解非匹配网格连接时的应力扰动问题。
您可能会问，既然在传热示例中的影响如此之小，为什么在这个示例中两个公式之间的差异如此之大？答案是，指定温度是平滑的，可以用形函数进行合理近似。而在网格不匹配的情况下，每个单元面上的位移场都由各自的形函数表示，并不具有连续导数。现在，弱约束的平均效应变得更加重要。
Nitsche 约束
COMSOL Multiphysics® 中的一些约束功能允许使用第三种类型的约束实施方法：Nitsche 方法。本文不展开理论细节，它也是一种弱约束类型，但不依赖于拉格朗日乘子的使用。它没有增加额外的自由度。以下是将该方法应用于同一个杆示例时的结果。
Nitsche方法有多种变体可供选择。这里使用的是默认的对称方法。可以看出，源点和终点的选择不再重要，误差甚至比使用弱约束时更小。Nitsche 约束不仅涉及连接边界上的节点，还涉及与边界上一个面的单元相连的所有节点，因此提供了更多的灵活性。
Nitsche 方法的缺点是，在其默认（也是最稳定的）的实现中，会产生一个非对称刚度矩阵，这可能会大大延长求解时间。如果问题中还有其他效应导致非对称贡献，这一缺陷则可忽略，因为对称刚度矩阵的优势已不复存在。
术语说明
在 COMSOL Multiphysics® 中，当选择使用弱约束时，该功能表示”使用拉格朗日乘子进行约束”。在少数情况下，例如在约束单点或常微分方程 (ODE) 自由度时，拉格朗日乘子只是一个数值而非场量。在这种情况下，情况与上述离散情况相同：没有新的近似值，唯一的影响是可以直接获取约束力。(不过，你可能并不认为这是真正的弱约束）。
在某些情况下，特别是在结构力学接口中，调用拉格朗日乘子的公式时并没有明确提及弱约束。例如，刚性连接件 功能有一个名为 计算反作用力 的选项。
你可能会问为什么需要一个选项来获取反作用力。原因是对于该功能和其他类似功能，除了使用拉格朗日乘子外，没有其他方法可以计算反作用力。另一方面，总是使用这种计算方法可能会导致一些不明显的问题，这将在后文中讨论。因此，需要用户进行交互控制。
拉格朗日乘子的解释
上文提到，对于离散情况，拉格朗日乘子可直接解释为节点反作用力。这一特性更为普遍，并且不局限于有限元方法或固体力学。拉格朗日乘子代表执行约束所需的某种作用。通常，拉格朗日 \mathcal {L } 代表一种能量。在这种情况下，拉格朗日乘子将与约束量在能量上共轭。
拉格朗日乘子的实际单位还取决于受约束物体的尺寸。对于固体力学，采用弱约束边界的拉格朗日乘子单位为 N/m2 = Pa。拉格朗日乘子场可以解释为为了保持指定位移而需要施加到边界上的牵引力场。但是，不应期望该域能准确表示约束处的局部应力场。不过，它在综合意义上是非常准确的。
为了说明反作用力的牵引力场，我们对下图中的短悬臂梁进行了研究。在固定端施加了弱约束。
接下来研究受约束端部的剪应力。根据分析解法，剪应力在横截面上呈抛物线分布。泊松比选择为零，以尽量减少约束效应。在给定的坐标系中，外加载荷会产生剪应力 σxy 和弯曲正应力 σxx。由于固定表面的法线位于 x 方向，牵引力分量 ty = σxy。因此，它应该用 y 方向上约束条件的拉格朗日乘子来表示。下图对这些结果进行了比较。
在所使用的网格分辨率下，应力和拉格朗日乘子均未精确反映真实解。存在一个显著差异：积分后，Y 方向总反作用力的相对误差在使用应力计算时为 3%，而使用拉格朗日乘子时为 2·10-12。
如果自行编写弱约束，拉格朗日乘数的物理意义可能因约束公式的不同而有强有弱。例如，假设需约束某点，使其变形后位于以原始位置为圆心、半径为 R 的圆周上。
以下展示了两种输入此类约束的不同方式：
在 弱约束 节点的设置中输入同一约束的不同方法。
这两种表达式都能以相似的迭代次数得到相同的解。然而，在第一种情况下，拉格朗日乘子的值很难解释。从量纲角度来看，它的单位是 N/m（因为它乘以单位为 m2 的约束条件）。在第二种公式中，它实际上是对径向位移的约束。计算出的拉格朗日乘子将是作用在支撑结构上的力。反作用力的方向并不明显，但会沿径向作用。因此，在物理上，该结构由一个无摩擦的圆形边界支撑。
不过，这个问题有两种可能的正确解—— 该点可能附着于圆周上的两个位置，本质上代表圆内侧或外侧的接触，最终解取决于初始条件。
使用弱约束时应注意的事项
约束冲突
在使用弱约束时最常见的问题是，若将其与标准逐点约束应用于同一自由度，二者无法共存。这种冲突可能导致”奇异矩阵””发现 NaN 或 Inf”等错误，甚至是完全错误的解。这种情况通常并非用户有意为之，但也可能在不经意间发生。
例如，如果在一条边界上使用弱约束，而在相邻边界上使用逐点约束。那么在公共边上就会出现冲突。通常，最简单的解决方法是将所有相连的约束条件转换为相同的表述（可能是弱约束条件，因为首先在某处选择它是有原因的）。有些约束的设置中还有 排除的边 和 排除的点 的部分。通过将冲突边添加到此类选择中，也可以解决问题。
在某些情况下，发生冲突的可能性并不明显。可能有一些约束条件是你没有想到的，因为它们或多或少是自动添加的。连续性 条件就是这样一种情况。
更微妙的是壳 接口中的情况，在 壳 接口中，每个节点都对旋转自由度（围绕壳表面法线的旋转）具有隐式约束。因此，如果在 壳 接口的任何对象上添加涉及旋转的弱约束（例如 固定约束 ），就会产生冲突。所有旋转约束（包括隐式约束）都必须改为弱约束。
拉格朗日乘子的单位
在大多数情况下，COMSOL Multiphysics® 中的拉格朗日乘子不带单位。这使得它们与几乎所有其他量不同。单位是隐含的，只能通过反作用力场的物理意义推断。如果需要在结果评估中频繁使用拉格朗日乘子，建议创建包含单位的中间变量。
添加一个带有单位的变量来表示拉格朗日乘子场。
变量缩放
有时，加入拉格朗日乘子会使非线性问题的迭代次数增加，甚至无法收敛。出现这种情况的最可能原因是容差处理不当。解决这一问题并确保容差得到正确处理的最佳方法是对拉格朗日乘子的自由度进行手动缩放。为此，请转到求解器序列中的 因变量 节点，并添加适当的缩放比例。缩放应提供拉格朗日乘子的数量级，即其代表的反作用力通量。
迭代求解器
如果使用迭代求解器求解线性方程组，拉格朗日乘子可能会引发一些问题。刚度矩阵的对角线上有零点，因此它不再是正定矩阵，大多数标准有限元公式都是如此。
因此，某些线性求解器和前置条件器不能用于求解弱约束问题，即共轭梯度迭代求解器和 SOR 类前置条件器和平滑器。您可以尝试另一种迭代求解器，并使用以拉格朗日乘数作为 Vanka 变量的 Vanka 算法，或者使用不完全 LU 因子化算法作为预处理器。
对于多物理场问题，一种可行方法是采用分离式求解策略；若模型规模允许，还可对涉及拉格朗日乘数的场使用直接求解器。
拉格朗日乘子自由度的奇异性
将标准约束条件转换为弱约束条件时，在数值上会表现良好。但是，如果编写自己的非线性约束条件，或使用约束变量的某些表达式，那么刚度矩阵的拉格朗日乘子部分可能会出现奇异性。
在上文使用的非线性约束示例中，当节点被约束移动到其原始位置周围一圈时，其原始配置实际上会失效。
这里的约束表达式实际上是一个以节点为中心的抛物面，因此当 u = v = 0 时，其梯度为零。
作为替代方案，您可以在拉格朗日乘子的自由度中添加一个人为的”刚度”。在上面的例子中，可以将表达式扩展为：u^2+v^2-R^2+1e-6*lm*(u^2+v^2 < 0.01*R^2)
拉格朗日乘子变量的名称为 lm。由于采用布尔表达式，额外的贡献只对小位移有效，因此不会影响真正的解。其作用是在刚度矩阵的对角线上放置一个小数值，以避免初始奇异性。
我们在下一篇博客的"弱约束条件"部分也使用了同样的方法：如何在仿真中设置边界条件。在该部分中，出现奇异性的原因是边界中只有随时间变化的部分具有指定值。在边界的其余部分，拉格朗日乘子的自由度仍然存在，但没有定义它们的方程。
扩展资源
COMSOL Multiphysics® 中的弱约束是强大的模拟技术，但要成功运用它们，需从数值角度理解其底层机制。
如果您想了解更多有关弱约束的信息，也可以查看以下博客:
评论 (0)