Menu
COMSOL 博客

利用最大值原理节省计算时间和资源

Author Image

作者Temesgen Kindo

2017年 5月 9日

当你的仿真计算量很大时，你会购买内存更大的计算机吗？如果求解耗时较长时，你会通宵运行它们吗？很多时候，你没有其他选择。但如果你有合适的工具，就可以通过利用数学原理找到更好的方法。今天，我们将展示如何在 COMSOL Multiphysics®软件中使用最大值原理节省计算资源和时间。

最大值原理简介

我们经常想知道某些量的最大值或最小值。当自由度数很大时，寻找极值的计算成本可能很高。如果使用极值作为输入，例如在边界条件或材料属性中，计算成本会更高。控制器设计就是这样一个例子。

假设我们事先掌握可能的关键位置的信息，那么，我们可以将搜索重点集中在这些位置。在瞬态问题中，如果可以仅存储可疑区域的数据，就可以节省内存。

挑战在于，对于许多问题，我们没有这样的信息。但是对于一些问题，我们可以证明极值出现在边界上。在数学中，这种性质被称为最大值原理。

COMSOL Multiphysics 提供了仅在边界上查找极值的方法，并允许我们将解存储在有限的一组域或边界上。

今天，我们将向您展示如何利用最大值原理高效求解和后处理。首先，我们将讨论最大值原理存在的条件，并简要叙述其数学证明。接下来，我们将演示 COMSOL Multiphysics 提供的使用这些原理的方法，包括边界上的最大值和最小值运算符、在表面上定义的物理场接口，以及边界元方法等。最后，我们将重点讨论最大值原理不成立的常见情况。

在开始之前，请注意，虽然最大值原理是有用的，但我们只应该在前提有效的情况下应用它们。否则，我们将会被路灯效应引入歧途。

路灯效应的示意图。
当你找不到某些东西时，就会发生路灯效果，因为你只在容易查找的地方搜索。它来自一个著名的笑话，讲一个人在路灯下寻找他的钱包，即使他知道他把它留在了附近的公园里。

推导最大值原理

让我们从一个简单的一维问题开始。在一段区间范围内，考虑二阶常微分方程

\frac{d^2u}{dx^2}=0

该方程可以代表几个一维稳态边界值问题，例如在没有平流或源/汇的情况下的热传导或化学输运。不论物理性质如何，其一般解都是一条直线

u(x) = ax+b,

其中，常数A和B可以根据边界条件确定。

由于解是一条直线，它的最大值和最小值将位于域（区间）的左端或右端。因此，如果我们对极值感兴趣，只需要检查左右两端的值。如果两个边界条件都是狄利克雷类型，指定u，我们甚至不需要求解方程来找到它的最大值和最小值。狄利克雷边界条件中较小的值为最小值，较大的值是最大值。这显然为我们节省了一些时间。

具有均匀、各向同性性质的对流扩散问题

在上面的分析中，我们利用对解的一般形式的知识得出有关极值位置的结论。对于更复杂的源项或高于一维的空间维度，这并不是一种实用的方法。与其推导出解的一般形式，不如使用初等微积分中的局部极大值或极小值的性质对以下线性偏微分方程做出类似的结论。

-\kappa\Delta u + \mathbf{v}\cdot \nabla u = f, \quad f<0, \textrm{ 在 } \Omega \textrm{中}

其中，\Delta是拉普拉斯算子\Delta u = \frac{\partial^2u}{\partial x^2}+\frac{\partial^2u}{\partial y^2}+\frac{\partial^2u}{\partial z^2}；\nabla是梯度算子\nabla u = (\frac{\partial u}{\partial x},\frac{\partial u}{\partial y},\frac{\partial u}{\partial z})；k是正标量（扩散系数或传导系数）；v是对流速度;f是源项。请注意，源项在域中的所有位置都是负数。

解在域的内部点能否达到最大值？从多变量函数的最优条件中，我们知道如果一个函数是平滑的，那么梯度在临界点应该为零。此外，在局部最大值处，所有二阶偏导数\frac{\partial^2u}{\partial x^2}, \frac{\partial^2u}{\partial y^2}, \frac{\partial^2u}{\partial z^2}都应为负数或零。这两个条件导致偏微分方程的左侧为正或零。

在内部局部最大值处，

-\kappa\Delta u + \mathbf{v}\cdot \nabla u \ge 0

这与源项相矛盾。在同一等式的右侧，

f<0

这意味着我们的解的最大值只能在边界上找到。我们刚才证明的被称为强最大值原理。

还有弱最大值原理，当允许源项为零时，该原理成立；即，f\le0是条件。在这种情况下，最大值要么出现在边界处，要么解始终是恒定的。无论哪种方式，如果要寻找最大值，只需寻找边界就足够了。我们在这里跳过弱最大值原理的证明。当f\geq0，我们可以推导出最小值的类似原理。

最后，如果没有源或汇，我们可以结合弱最大值和最小值原理，得出解的最大和最小值都应该在边界上的结论。习惯上，我们使用最大值原理一词来指代最大值原理和最小值原理，具体指代哪个应结合上下文理解。

请注意，我们需要解具有平滑度才能得出上述结论。

各向异性，非均质相对流扩散问题

到目前为止，我们的结论可用于稳态对流扩散型物理问题，如传热、电流或化学传输，前提是电导率（扩散率）是均匀的和各向同性的，并且对流速度v也是均匀的。现在让我们取消这些限制。

通常，稳态对流扩散问题的形式为

\nabla \cdot (-\mathcal{D} \nabla u + \mathbf{v}u) = f,

其中，D是电导率（扩散）张量，可以是各向异性和异构性的；v是潜在的不均匀对流速度;f是热（化学/电流）源。

建立平滑解的弱最大值原理的要求如下：

  1. 扩散张量D是正定的
  2. 对流速度v是无发散的；即\nabla \cdot \mathbf{v} = 0
  3. 源项为非正数；即f(\mathbf{x}) \leq 0, \quad \forall \mathbf{x} \in \Omega

从物理上讲，这涵盖了不可压缩对流速度的传输问题，例如传热、反应流和电荷传输问题。无论如何，如果f \geq 0，可以得出最小值应该在边界上的结论。

考虑电加热器、保险丝和其他导体中的焦耳加热问题(详见这个视频)，其中电气部件中的电阻损耗会导致加热。固体传热方程由下式给出

\nabla \cdot (-\mathcal{\kappa} \nabla T ) = \sigma|\nabla V|^2,

其中，k和σ是热导率和电导率，T和V是温度和电势。

源项\sigma|\nabla V|^2始终为非负数。根据我们到目前为止的讨论，可以得出结论，最低温度必须在模型的边界上。请参阅下图，其中标注了体积和表面最小值的位置。这两个位置重叠。

在 COMSOL Multiphysics 中模拟焦耳热。
标注了最小整体温度的COMSOL模型。
标注了最小值边界温度的 COMSOL 模型。

焦耳热问题中的温度分布(左)、最低总温度的位置(中心)和最低边界温度的位置(右)。

瞬态问题

瞬态问题大致可分为抛物线问题和双曲线问题。第一组包含扩散型现象，例如传热，化学反应和地下水流。第二组包含波型问题，如声学、电磁波和结构动力学。

对于抛物线问题，我们可以证明解的总体最大值或最小值是在初始时间或边界处遇到的。与稳态问题一样，源项的符号在这里起着决定性的作用。双曲线问题通常不遵循最大值/最小值原理。

考虑以下具有适当的初始和边界条件的抛物线方程

\frac{\partial u}{\partial t}+\nabla \cdot (-\mathcal{D} \nabla u +\mathbf{v}u) = f

如果扩散张量D是正定的，那么对流速度v是无发散的，源项f \leq 0。可以证明解在所有时间和点的整体最大值（而非任意给定时刻的最大值）要么出现在边界处，要么在初始时刻出现。

为了简化这里的讨论，让我们考虑上述抛物线方程的一维版本的强最大值原理。

\frac{\partial u}{\partial t}+\frac{\partial}{\partial x}(-D \frac{\partial u}{\partial x} +vu) = f, \quad x \in [0,L], \quad t\in[0,T], u(x,0) = u_0(x)

用导数的乘积法则扩展此等式，得到

\frac{\partial u}{\partial t}-\frac{\partial D}{\partial x}\frac{\partial u}{\partial x} -D \frac{\partial^2 u}{\partial x^2} +v\frac{\partial u}{\partial x}+u\frac{\partial v}{\partial x} = f, \quad x \in [0,L], \quad t\in[0,T], u(x,0) = u_0(x)

为了证明强极大值原理，我们将假设u的最大值位于时空域的内部点[0,L] \times [0,T]，并说明这一假设如何导致矛盾。

在内部最大值处，平滑解将具有\frac{\partial u}{\partial x}=\frac{\partial u}{\partial t}=0。此外，对于一维问题，不可压缩的平流速度将处处都有\frac{\partial v}{\partial x}=0。因此，在内部最大值下，上述等式简化为

– D \frac{\partial^2 u}{\partial x^2} = f

由于D为正，f为负，这意味着在内部最大值\frac{\partial^2 u}{\partial x^2} >0，这与微积分的通常要求相矛盾，即在内部最大值时，\frac{\partial^2 u}{\partial x^2} \leq 0。

这就使我们把时空域的边界作为最大解的可能位置。这个复合域的边界包括初始时间 t = 0、结束时间t=T以及 x=0 和x=L处的空间边界。

时空域图表。
一维瞬态偏微分方程的时空域。

我们已经将最大解的位置限定在上图所示的四个边界之内。我们最初的观点是，最大解可能出现在空间边界或者初始条件所定义的时刻。这就排除了那些在空间内部但在终端时间t=T的点；也就是上述时域中矩形的顶部。我们来说明为什么会这样。

假设最大值位于时间边界t=T的一个空间内部点\hat{x} \in (0,L)。如果我们从那个点向左或向右移动，即在空间中移动而不是在时间上移动，解不应该增加。因此\frac{\partial u}{\partial x}=0, \frac{\partial^2 u}{\partial x^2} \leq 0，然而，第一个时间偏导数不必为零，因为我们只能从终端时间向一个方向移动。这里的限制是，要想在终端时间发生最大值，时间偏导数在那里不能为负。也就是说，我们需要\frac{\partial u}{\partial t} \ge 0。这与偏微分方程一致吗？将第一个空间导数设置为零，我们得到

\frac{\partial u}{\partial t}- D \frac{\partial^2 u}{\partial x^2} = f \Rightarrow \frac{\partial u}{\partial t}= f+D \frac{\partial^2 u}{\partial x^2}

由此，我们可以看到，关于源项的起始假设和对二空间偏导数的限制导致了一个负时间导数。这个基于偏微分方程的结论与我们根据局部最大值的要求所作的陈述相矛盾。相反，在初始时间使用类似的推理表明，在初始条件下可能的内部最大值处的时间导数是可能的。

在计算工作中，这意味着如果我们的目标是寻找最大的解，应该将搜索限制在边界内，但除了初始时刻，也应该考虑内部点。但是请注意，这种说法仅适用于所有时间和积分范围内的最大值。给定时间的最大值仍然可能在内部点。

系统和高阶方程

最大值原理通常不适用于高于 2 阶的方程，例如双谐波方程或二阶系统（例如二维和三维应力分析）。寻找最大值原理的基本思想是得到解向量u的一些标量函数P（u）来满足泊松方程。这样的标量函数被称为 P 函数，以 L. E. Payne 的名字命名，他对这一思想做出了重大贡献。工程上的挑战是找到具有实际意义的 P 函数。

例如，考虑稳态应力分析。平衡方程由方程组给出

\sigma_{ij,j}+f_i = 0, i,j=1,2,3,

式中，σij和fi分别是每单位质量的应力张量和体力的分量。

在线性弹性的情况下，应力和应变通过本构关系相关

\sigma_{ij}=\lambda \epsilon_{kk}\delta_{ij}+2\mu \epsilon_{ij}

其中，\epsilon_{ij}是应变张量的分量，且与位移ui密切相关，这种关系通过应变-位移关系来定义

\qquad \epsilon_{ij}=\frac{1}{2}(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i})

如果材料属性λ和μ是恒定的，并且体力f是螺线型的（无发散），例如非相对论应用中的引力情况，我们可以证明，体积变化\epsilon_{kk}=\epsilon_{11}+\epsilon_{22}+\epsilon_{33}，也称为膨胀，同时满足最大值和最小值原理。

我们可以通过代入平衡方程中的本构关系来证明这一点，得到

\lambda \epsilon_{kk,i}+2\mu\epsilon_{ij,j}=0, \quad i,j=1,2,3

并取上述向量方程的散度得到

\lambda \epsilon_{kk,ii}+2\mu\epsilon_{ij,ji}=0, \quad i,j=1,2,3.

从微分阶的应变-位移关系和微分次序互换性可以看出\epsilon_{ij,ji}=\epsilon_{jj,ii}。这最终导致

\epsilon_{kk,ii}=\Delta \epsilon_{kk} = 0

由于膨胀满足拉普拉斯方程，因此其最大值和最小值都必须在边界处实现，除非膨胀在物体上是均匀的。对本构关系的追踪表明，在线性弹性、各向同性材料属性和无发散体积力假设下，平均应力p=\frac{1}{3}\sigma_{kk}也是如此。

虽然上述结果在理论上很有吸引力，但工程师很少对寻找膨胀或平均应力的极值感兴趣。另一方面，最大 von Mises 应力和最大剪切应力是重要的标量。例如，弹性极限（屈服标准）通常根据这两个量中的任何一个来确定。不幸的是，对于这些和其他实际利益量，我们没有广泛适用的最大值原则。但针对特殊情况仍然可以推导出相应的最大原则。

一个很好的例子是对扭转的均匀梁分析。如果我们选择z轴作为梁的轴，应力分析问题可以根据在横截面 Ω 上定义的应力函数\phi(x,y)重新表述。在没有体积力的情况下，应力函数满足

\frac{\partial^2 \phi}{\partial x^2}+\frac{\partial^2 \phi}{\partial y^2} =-2, \quad (x,y) \in \Omega

应力的非零分量为

\sigma_{xz}=\mu\theta\frac{\partial \phi}{\partial y}, \quad \sigma_{yz}=-\mu\theta\frac{\partial \phi}{\partial x}

其中，θ是梁每单位长度的扭转角度。

COMSOL Multiphysics 的一张屏幕截图，显示了一个使用梁横截面接口的模型。
使用梁横截面接口计算的扭转应力分布。

我们对最大剪切应力感兴趣。请注意，由于扭转问题不存在法向应力，因此剪切应力的大小与 von Mises 应力成正比。对于剪切应力的大小τ，我们有

\tau = \sqrt{\sigma_{xz}^2+\sigma_{yz}^2} \Rightarrow \tau^2= \sigma_{xz}^2+\sigma_{yz}^2=(\mu\theta)^2|\nabla\phi|^2

我们以|\nabla\phi|^2的散度来研究最大值原理的可能性。使用指数表示法

div(|\nabla\phi|^2)=\frac{\partial^2}{\partial x_j\partial x_j}(\frac{\partial \phi}{\partial x_i}\frac{\partial \phi}{\partial x_i})=2\frac{\partial}{\partial x_j}(\frac{\partial \phi}{\partial x_i}\frac{\partial^2\phi}{\partial x_i\partial x_j})=2\frac{\partial^2\phi}{\partial x_i\partial x_j}\frac{\partial^2\phi}{\partial x_i\partial x_j}+2\frac{\partial\phi}{\partial x_i}\frac{\partial^3\phi}{\partial x_j\partial x_i\partial x_j}

如果我们交换上述等式最后一项中的微分顺序，可以得到

div(|\nabla\phi|^2)=2\frac{\partial^2\phi}{\partial x_i\partial x_j}\frac{\partial^2\phi}{\partial x_i\partial x_j}+2\frac{\partial\phi}{\partial x_i}\frac{\partial}{\partial x_i}(\frac{\partial^2\phi}{\partial x_j\partial x_j})

括号中的项是应力函数的散度。根据我们最初使用的等式，这个值是 -2。因此，上面的第二个右手项为零。第一个右手项是非负的，因为它是\phi的 Hessian（二阶偏导数矩阵）的分量的平方和。此外，它不能为零，因为 Hessian 的对角线加起来是一个非零值 -2。

因此，我们有

div(\tau^2)>0 \quad \forall (x,y) \in \Omega

这意味着剪切应力的大小在梁横截面的边界处达到最大值。许多关于线弹性的教科书都包含了针对简单横截面和边界条件的 Saint-Venant 扭转问题的解析解。这些显式解表明，最大的剪切应力出现在边界处。我们在这里所做的是针对任意横截面和边界条件得出这一结论，而没有明确解决边界值问题。请参阅下面承受扭转的工字梁的比较。

一个梁模型，标注了整体最大剪应力。
标注了边界最大剪应力的梁模型。

比较扭转作用下梁中整体最大剪应力(左)和边界最大剪应力(右)的位置。

如何利用 COMSOL Multiphysics®软件中的最大值原理

假设感兴趣的物理量存在最大值原理，或者你只想关注边界处的物理量，那么你可能就只对在问题边界上的极值感兴趣。COMSOL Multiphysics 有哪些功能可以帮助你做到这一点呢？

后处理阶段的最大值和最小值计算

如果从模型生成器或 GUI 的菜单中进入结果>派生值，就会找到在体积、面或直线上定义的最大值和最小值运算符。当你对边界上的极值感兴趣时，可以使用如下所示的表面或线算子，具体取决于问题的空间维度。

模型开发器的截图显示了最大值计算选项。
可以在不同的几何实体级别上完成最大值和最小值计算。

最大值和最小值组件耦合算子

上述算子的使用仅限于后处理。如果你必须在物理场中使用极值，例如实现控制器，则可以在模型的定义部分中定义最大值或最小值组件耦合算子。这些运算符还允许在体积、面、曲线或点之间选择几何实体。添加算子后，在设置窗口中选择几何实体级别，如下所示。

显示COMSOL软件中组件耦合运算符的几何实体选择的拼图。
组件耦合运算符的几何实体选择。

除了由于自由度较低而节省时间和内存外，边界耦合运算符对所涉及的矩阵结构造成的干扰也比域耦合运算符小。如果你知道要求解的方程具有最大值或最小值原理，并且需要在问题表述中使用极值，那么 COMSOL Multiphysics 提供了一种经济的替代方案，即使用边界耦合算子而不是域耦合算子。

将解存储在几何体的选定零件上

在一个瞬态的问题中，你只对边界数据的后处理感兴趣，当模型较大且存在内存限制时，只存储几何体的关键部分的解可能是一种既经济又必要的方法。我们在之前的博客中讨论过，在 COMSOL Multiphysics 中，有不止一种方法可以做到这一点，。

表面疲劳

如果你正在进行疲劳分析，并且知道临界行为位于几何体的边界处，COMSOL Multiphysics 将为你提供仅限于边界的疲劳分析功能。

在COMSOL Multiphysics中突出域和边界的疲劳评估的注解屏幕截图。
域与边界的疲劳分析。

边界元法

COMSOL Multiphysics 主要基于有限元方法。对于某些问题，它支持边界元方法。在最新版本的软件(5.3 版)中，边界元方法可用于静电、腐蚀和用户定义的方程。

边界元方法解决的问题通常不含源项，因此很可能同时满足最大值和最小值原理。不过，利用最大值原理并非选择边界元方法的充分理由。但当采用边界元方法时，问题在边界上极可能存在极值。

注意事项

最大值原理在应用于其前提适用的问题时非常有用。然而，我们应该避免将其应用于这些前提条件不适用的问题。例如，源项符号的不连续性和变化就不复合这些前提。

不连续性

在我们的数学分析中，我们假设系数的平滑度。如果存在不连续性，则这些原则不适用于整个域。然而，我们可以将域细分为若干连续的区域，并在每个子域上应用最大值原理。实际上，这意味着你必须包括内部边界和外部边界。具有不同材料属性的两种材料之间的界面就是一个例子。

我们之前提到过，在 COMSOL Multiphysics 中，疲劳分析可以在域或面上进行。在许多问题中，疲劳失效始于表面。但是在接触力学中，已知疲劳失效始于次表面点。有关这方面的更多信息，包括有限元网格划分策略的含义，请参阅以下有关接触疲劳建模的博客文章。

更改源项的符号

在非零源项的问题上，最大值原理要求源项具有非负性或非正性。因此，如果源项在域内改变符号，最大值原理将不再适用。例如，在焦耳热问题中，电磁热源始终为正。另一方面，在化学反应中，相同的化学物质在域的一部分区域产生，而在其他区域消耗。除非我们确定反应在整个域内都是单一地产生或消耗某种化学物质，否则不应使用最大值原理。

离散化的影响

我们对最大值原理的推导直接使用偏微分方程。离散化过程可能会破坏这一原理的适用性。因此，当发现某个理论上应遵循最大值原理的问题在数值求解中未能体现时，我们应检查所使用的数值方法是否遵循离散最大值原理。此外，还需仔细检查是否存在数值误差。

关于使用最大值原理的结论性思考

最大值原理在偏微分方程的分析中占据重要地位，常用于证明解的存在性与唯一性、确定误差界限以及实现其他类似目标。在这篇博客中，我们展示了如何利用这些原理进行高效的计算分析。

请注意，本文第一部分中提供的数学讨论针对更广泛的受众。有关更严格的讨论，请参阅偏微分方程分析中的标准文本。以下是我个人推荐的有关这个主题的一些参考：

  • Lawrence C. Evans,Partial Differential Equations, American Mathematical Society, 2010
  • Rene Sperb,Maximum Principles and Their Applications, Academic Press, Inc., 1981

如果你对本主题或使用 COMSOL Multiphysics 有任何疑问，请联系我们。

联系我们

延伸阅读

博客分类

相关主题

评论 (0)

正在加载...
浏览 COMSOL 博客
类别
标记
Baidu
map