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

2017年 5月 9日

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

最大值原理简介

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

想象一下,事先掌握可能的关键位置的信息。现在,我们可以将寻找重点放在这些位置。在瞬态问题中,如果我们可以仅存储可疑区域的数据,就可以节省内存。

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

COMSOL Multiphysics 提供了仅在边界上查找极值的工具。该软件还允许我们将解存储在有限的一组域或边界上。

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

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

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

推导最大值原理

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

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

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

u(x) = ax+b,

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

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

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

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

-\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σ 是热导率和电导率,TV 是温度和电势。

源项 \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 的一些标量函数 Pu)来满足泊松方程。这样的标量函数被称为 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} 是应变张量,而应变张量又通过应变-位移关系

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

与位移 ui 有关

如果材料属性 λμ 是恒定的,并且体力 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,i}+2\mu\epsilon_{ij,j}=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 博客