基于方程的轴对称组件建模指南

Temesgen Kindo 2016年 10月 5日

柱坐标系对于高效求解和后处理旋转对称问题而言很有用。COMSOL Multiphysics® 软件为轴对称物理场接口中的柱坐标系提供了内置支持。当您使用数学接口对定制的偏微分方程(partial differential equation,简称 PDE)进行定义时,请务必仔细辨明它们的意义。偏微分方程接口默认采用笛卡尔坐标系来描述偏微分方程,因此您需要手动将笛卡尔坐标系变换为柱坐标系。在下文中,我们将介绍当使用自定义的偏微分方程时如何实现这种坐标变换。

一个轴对称的热传导问题

今天的案例涉及到求解圆柱体上的轴对称传热问题。该模型取自 NAFEMS 标准模型集并在 COMSOL Multiphysics 中求解,同时还使用了固体传热 接口。我们将采集模型的尺寸和材料属性,并借助一般形式偏微分方程 接口重现结果。实质上,我们正在通过比较手动偏微分方程接口与固体传热 接口(该接口会自动处理轴对称性),从而验证手动偏微分方程接口的方法是否可以解决上述问题。

刚性固体上的稳态温度分布 T 的方程为

\rho C_p\mathbf{u}\cdot \nabla T + \nabla \cdot (-\kappa \nabla T)=Q,

如果平移速度 u 始终为零,则热容 Cp 不会产生影响,导热系数 κ 便是我们唯一需要指定的材料属性。我们使用参数 kappa 来表示导热系数,并在固体传热 接口中获得如下图所示的设置。下图同时还绘制了一组特定温度和通量边界条件下的解。

Heat Transfer in Solids interface 基于方程的轴对称组件建模指南
使用 固体传热接口解决轴对称热传导问题。您可以从官网的“案例下载”中找到此标准模型。

现在,让我们使用一般形式偏微分方程 接口来求解相同有限元网格中的相同问题。该接口是 COMSOL Multiphysics 中的数学接口之一,可用于求解定制方程。当您要求解的方程并非建立在任何一个物理场接口中时,您便可以使用数学接口求解代数方程、常微分方程及偏微分方程。新添加的方程与线性或非线性无关,既可以自行求解,也可以通过与其他预置的物理场接口耦合进行求解。因内置接口本身支持传热方程,所以可便捷地将使用传热方程得到的结果与使用我们自己的偏微分方程得到的结果进行比较,这会是一个非常良好的出发点。

一般形式偏微分方程 接口的模板为

e_a\frac{\partial^2 u}{\partial t^2} + d_a\frac{\partial u}{\partial t} + \nabla \cdot \Gamma = f,

为了定义我们的数学模型,我们必须提供质量系数 ea,阻尼系数 da,通量 Γ 以及源 f

如果我们将该模板与固体传热 接口中的方程进行比较,其中因变量 u 代表温度,则我们应该将通量定义为 \Gamma = -\kappa \nabla u。为了实现这一点,我们将偏微分方程接口中通量的 r 分量和 z 分量分别指定为 -kappa*ur-kappa*uz。稳态分析并不涉及质量系数和阻尼系数。

如下图所示,使用上述设置获得的解不同于使用固体传热 接口获得的解。原因何在?

General Form PDE interface 基于方程的轴对称组件建模指南
利用 一般形式偏微分方程接口求解轴对称热传导问题。这里不考虑曲线坐标系的影响。

造成这种差异的原因可简要解释为:固体传热 接口将 \nabla \cdot (-\kappa \nabla T) 理解为热通量的散度,然而偏微分方程接口将 \nabla \cdot \Gamma 理解为

\frac{\partial \Gamma_r}{\partial r} + \frac{\partial \Gamma_z}{\partial z}.

因此,该表达式并不是柱坐标系中 Γ 的散度的表达式。关于如何修正这一问题,我们将在本文的后半部分进行演示,接下来我们再更为详细地解释一下为何会造成这种差异。

守恒定律

在 COMSOL Multiphysics 中,各类物理场接口求解的方程均是物理定律的数学抽象化表达。这些定律通常表示为守恒定律或原理解释,用于描述发生在某个域或跨域边界的活动所引起的某种量的变化。虽然守恒定律适用于所有材料,但是不同材料对域力和边界通量的响应程度也有所不同。材料响应由所谓的本构方程或状态方程进行指定。我们必须确保本构方程不违反热力学第二定律。守恒定律与有效的本构方程配合使用,可提供获取适定性数学模型所需的充足信息。

举例来说,刚性物体中的热传导遵循热能守恒定律。热能变化率等于域中热源供热速率加上跨边界的热通量速率。根据经验观察我们可以看出,固体中的热通量与温度梯度成正比,并且从高温区域向低温区域传播。这就是本构方程。由此,我们可以使用多元微积分,将各向同性材料的传热方程写为

\rho C_p\frac{\partial T}{\partial t}+\rho C_p\mathbf{u}\cdot \nabla T + \nabla \cdot (-\kappa \nabla T)=Q,

其中,T 表示温度,是我们须主要追踪的变量;材料参数 pCpK 分别用于描述密度、热容及导热系数;u 为平移速度;Q 为单位体积热源。边界条件也可用于描述边界处发生的情况。

坐标系

在推导出与上述方程相似的数学模型后,下一步就是求解主要因变量和其他您希望获取的物理量,例如热传导过程中温度随时间和空间的分布。我们选择坐标系来对空间中的点进行辨别。合适的坐标系有助于后续分析,而一个不合适的坐标系会给我们的工作增加不必要的复杂性。无论选择什么坐标系,都需要保证方程的物理意义始终不变。

以包含温度梯度 \nabla T 的热传导方程为例,它在笛卡尔坐标系中的表达式为

\nabla T = \frac{\partial T}{\partial x}\mathbf{e}_x + \frac{\partial T}{\partial y}\mathbf{e}_y+\frac{\partial T}{\partial z}\mathbf{e}_z,

而在柱坐标系中的表达式为

\nabla T = \frac{\partial T}{\partial r}\mathbf{e}_r + \frac{1}{r}\frac{\partial T}{\partial \phi}\mathbf{e}_{\phi}+\frac{\partial T}{\partial z}\mathbf{e}_z.

在第一种情况下,温度标量关于自变量的偏导数直接表述了温度梯度的各个分量,而在像圆柱坐标系这样的曲线坐标系中则不是这样的。根据上文提到的温度梯度的物理意义,该矢量指向温度增幅最大的方向,且大小等于增长率,这是应该保持不变的。为了确保这种恒定性,我们必须使用协变导数来取代常规偏导数。在笛卡尔坐标系中,协变导数和偏导数一致;但在曲线坐标系中,它们并不一致。

另一个微分算子为散度算子 \nabla \cdot。在传热方程中,该算子作用于通量 -\kappa \nabla T。使用笛卡尔坐标系时,为了计算通量的散度,只有通量的分量 -\kappa \frac{\partial T}{\partial x}-\kappa \frac{\partial T}{\partial y}-\kappa \frac{\partial T}{\partial z} 需要微分。基矢 exeyez 在从一点到另一点的微分过程中保持不变。在柱坐标系及其他曲线坐标系中,基矢会发生变化。因此,求解散度还需要求解基矢的导数。而协变导数考虑到了这一点,在使用曲线坐标系时,将偏导数简单相加会丢失散度的物理意义,这反映出了矢量从特定点发散量的多少。

\nabla \cdot (\kappa \nabla T) = \frac{\partial }{\partial x}(\kappa \frac{\partial T}{\partial x}) + \frac{\partial }{\partial y}(\kappa \frac{\partial T}{\partial y}) + \frac{\partial }{\partial z}(\kappa \frac{\partial T}{\partial z}) = \frac{1}{r}\frac{\partial }{\partial r}(r \kappa \frac{\partial T}{\partial r}) + \frac{1}{r}\frac{\partial }{\partial \phi}( \frac{\kappa}{r}\frac{\partial T}{\partial \phi})+\frac{\partial }{\partial z}(\kappa \frac{\partial T}{\partial z}).

同样地,点积 \mathbf{u} \cdot \nabla T 应该始终与两个矢量的大小和二者夹角的余弦的乘积保持一致。在笛卡尔坐标系中,这仅仅是两个矢量的对应分量的乘积之和。而在曲线坐标系中,坐标系的度规张量也混于其中。

由于曲线坐标系会增加数学计算的复杂程度,所以您可能会十分好奇,为何我们宁可使用其他坐标系,也不想使用笛卡尔坐标系。不过,正如我们在这里强调的,曲线坐标系对于某些应用来说非常有用。

应在什么情况下使用曲线坐标系

设想这样一个三维问题:模型的几何结构、材料属性、边界条件及热源都是轴对称,因而得到的解也是围绕此轴旋转对称的。如果我们将坐标系中的 z 轴方向用作对称轴,那么所有与 Φ 相关的偏导数都会消失,所以留给我们的只有一个 rz 平面内的二维问题。这带来的好处是,二维问题中的方程会比在笛卡尔坐标系中的方程更为简单,这是因为笛卡尔坐标系中的方程包含全部三个空间偏导数。对这类问题进行有限元建模时,轴对称公式的使用让我们可以使用二维网格,而非三维网格,这将会极大地节省内存和时间。

在材料属性(例如导热系数 κ)并非各向同性的情况下,曲线坐标系依然十分有帮助。如果在某些方向为各向异性,我们可以使坐标系优先与这些方向对齐,从而简化材料属性输入。注意,COMSOL Multiphysics 中的曲线坐标系 接口可生成非标准坐标系。您可以在COMSOL 博客中阅读更多关于如何使用曲线坐标系 接口如何在曲线坐标系中求解各向异性问题的信息。

在下文的示例中,我们将重点介绍使用标准的柱坐标系。

比较轴对称组件中的物理场接口与偏微分方程接口

在 COMSOL Multiphysics 中,有多个基于物理场的接口可用于求解由一个或多个守恒定律推导的方程。例如,在上文提到的刚性固体传热的例子中,我们遵循了热能守恒定律。在进行等温应力分析时,我们遵循的是质量守恒、线性动量守恒及角动量守恒定律。当您使用这些物理场接口中的任意一个时,软件都会通过使用对应笛卡尔坐标系或曲线坐标系的表达式来保持微分算子的物理意义。也就是说,在物理场接口中我们使用的是协变导数而非偏导数,正是这一原因才使不同坐标系下的物理场保持不变。

但是,当您使用系数型偏微分方程一般形式偏微分方程 接口时,软件会采用偏导数。

换句话说,物理场接口将 \nabla u\nabla \cdot \Gamma\nabla \times \Gamma 分别作为物理标量 u 或高阶张量 Γ 的梯度、散度及旋度。另一方面,在偏微分方程接口中,这些张量的意义不再适用,并且关于自变量的偏导数替代了这些算子。举个例子,自变量或许不能表示物理坐标。

现在我们来分析一下 \nabla \cdot \Gamma。在物理场接口中,它代表散度;但是在偏微分方程接口中,它代表三维组件中的 [\frac{\partial}{\partial x}, \frac{\partial}{\partial y}, \frac{\partial}{\partial z}]\cdot[\Gamma_x,\Gamma_y,\Gamma_z] 及二维轴对称组件中的 [\frac{\partial}{\partial r}, \frac{\partial}{\partial z}]\cdot[\Gamma_r, \Gamma_z]。第一个表达式可推出 \frac{\partial \Gamma_x}{\partial x} + \frac{\partial \Gamma_y}{\partial y} + \frac{\partial \Gamma_z}{\partial z},它与散度一致,第二个表达式可推出 \frac{\partial \Gamma_r}{\partial r} + \frac{\partial \Gamma_z}{\partial z},但这不是 Γ 的散度。

柱坐标系中散度的表达式为

\nabla \cdot \Gamma = \frac{1}{r}\frac{\partial (r \Gamma)}{\partial r} + \frac{1}{r}\frac{\partial \Gamma_{\phi}}{\partial \phi}+\frac{\partial \Gamma_z}{\partial z}.

针对轴对称问题,此总和中的第二项会消失,仅剩下

(1)

\nabla \cdot \Gamma = \frac{1}{r}\frac{\partial (r \Gamma_r)}{\partial r} +\frac{\partial \Gamma_z}{\partial z} = \frac{\partial \Gamma_r}{\partial r} + \frac{\partial \Gamma_z}{\partial z}+\frac{\Gamma_r}{r}.

偏微分方程接口中的表达式 \nabla \cdot \Gamma 并不包括最后一项,因为算子 \nabla \cdot 可被简单地表述为偏导数之和。为了得到柱坐标系中的物理量的散度,我们必须对方程 (1) 中的最后一项进行补偿。接下来,我们将借助示例来展示一个通过使用源项实现补偿的方法。

轴对称热传导问题的解决方案

现在,让我们回到最初的问题,我们要解决的是一个稳态问题

\nabla \cdot \Gamma = f,

其中,\nabla \cdot 为散度算子。对于轴对称问题来说,我可以得出

\frac{\partial \Gamma_r}{\partial r} + \frac{\partial \Gamma_z}{\partial z}+\frac{\Gamma_r}{r}=f.

这相当于

\frac{\partial \Gamma_r}{\partial r} + \frac{\partial \Gamma_z}{\partial z}=-\frac{\Gamma_r}{r}+f.

方程左边在轴对称组件中的一般形式偏微分方程 接口中可被理解为 \nabla \cdot \Gamma。方程右边额外添加了一个源项,如下方屏幕截图所示。现在,我们得到的解便与通过固体传热 接口获得的解相匹配了。

Proper General Form PDE interface 基于方程的轴对称组件建模指南
使用 一般形式偏微分方程接口求解轴对称热传导问题。这里考虑了曲线坐标系的影响。

当然,源项中添加的项并非物理热源。它仅仅是求解散度所需的协变微分和利用偏微分方程接口求解的偏微分之间缺失项的补偿。一种较好的解决方法是在“偏微分方程”设置窗口中添加该项,在“模型开发器”中添加物理热源的操作方法为:右键单击偏微分方程接口并选择“源”。

在该示例中,我们只使用了散度算子。使用其他微分算子(例如旋度)求解物理场问题时,您也同样需要执行类似的补偿。

这些调整适用于嵌在欧氏空间中的曲线坐标系。如果您使用的底层空间并非欧氏空间,那么便根本不存在笛卡尔坐标系。在这类情况下,即使您使用的是非轴对称的二维或三维组件,仍需要注意这一点。COMSOL Multiphysics 的内置物理场接口可解决这类问题,这些接口包括结构力学领域的 接口和 接口、流体力学领域的薄膜流动 接口、电磁学领域的电流,壳 接口等。您可以查看我们官网的产品规格表,获取更为完整的列表。

将物理场接口用作偏微分方程模板

另一种基于方程的建模方法是将物理场接口用作处理具有相似数学结构问题的偏微分方程参考模板。举例来说,为了求解具有对流-扩散-反应特征的物理问题,我们可以使用传热 接口或稀物质传递 接口。这些接口能够保持轴对称组件中微分算子的张量意义。您只需使用在数学意义上与想添加的项相匹配的边界条件和域条件即可。

这种方法的唯一缺点就是变量的单位可能与已有的单位不匹配。在这种情况下,您可以使用无量纲化 法。在获得无量纲方程之后,您可以前往“模型开发器”的根节点并将单位系统 设置为,从而使您的 COMSOL 模型无量纲化。

nondimensionalization 基于方程的轴对称组件建模指南
通过无量纲化,您可以使用一个物理场接口求解具有类似数学结构,但尺寸和类型不同的问题。

关于基于方程的轴对称组件建模的总结性思考

借助 COMSOL Multiphysics 中的系数型偏微分方程一般形式偏微分方程 接口,您可以建立偏微分方程来求解一些软件中尚未包含的新问题。这些偏微分方程可能推导自某个物理问题,也可能不是。因此,偏微分方程接口中的微分算子被刻意设置得比较简单,并且不会自动转化为张量算子。当使用曲线坐标系求解某个物理问题(例如想要利用轴对称)时,您必须确保输入到软件模型样板中的项能准确表述您的物理问题。COMSOL Multiphysics 的基于物理场的接口便可处理这一点。但在偏微分方程接口中,软件无法为您的方程赋予任何意义,这需要您自己完成。

解决这个问题的方法之一是通过添加额外的源项来平衡偏导数和协变导数之间的差异。另一种方法是使用与方程具有相同数学结构的现有物理场接口。如果您对上述解决方案或此话题有任何疑问,请立即联系我们

如需更多关于曲线坐标系中的微分算子及表面偏微分方程的信息,请参阅有关张量运算或微分几何的书籍。以下列出的是我个人非常喜欢的书,排名不分先后:

  • Rutherford Aris, Vectors, Tensors, and the Basic Equations of Fluid Mechanics, Dover Publications, Inc., 1989
  • Pavel Grinfeld, Introduction to Tensor Analysis and the Calculus of Moving Surfaces, Springer Science+Business Media, 2013
  • Mikhail Itskov, Tensor Algebra and Tensor Analysis for Engineers: with Applications to Continuum Mechanics, Springer-Verlag, 2013

博客标签

技术资料
加载评论……