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

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来表示导热系数,并在固体传热接口中获得如下图所示的设置。下图同时还绘制了一组特定温度和通量边界条件下的解。

借助 COMSOL Multiphysics 中的固体传热接口来解决轴对称传热问题。
使用固体传热接口解决轴对称热传导问题。您可以从官网的“案例下载”中找到此标准模型。

现在,让我们使用一般形式偏微分方程接口来求解相同有限元网格中的相同问题。该接口是 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。稳态分析并不涉及质量系数和阻尼系数。

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

截图展示了不考虑曲线坐标系的影响时,一般形式偏微分方程接口的使用。
利用一般形式偏微分方程接口求解轴对称热传导问题。这里不考虑曲线坐标系的影响。

造成这种差异的原因可简要解释为:固体传热接口将\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。方程右边额外添加了一个源项,如下方屏幕截图所示。现在,我们得到的解便与通过固体传热接口获得的解相匹配了。

屏幕截图展示了正确的一般形式偏微分方程接口。
使用一般形式偏微分方程接口求解轴对称热传导问题。这里考虑了曲线坐标系的影响。

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

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

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

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

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

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

截图展示了无量纲化的使用。
通过无量纲化,您可以使用一个物理场接口求解具有类似数学结构,但尺寸和类型不同的问题。

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

借助 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

博客分类


评论 (24)

正在加载...
飞龙 岳
飞龙 岳
2019-03-25

对phi的弱形式二阶偏导如何实现?usys2.phisys2.phi软件显示有误

Haili Wang
Haili Wang
2019-04-15 COMSOL 员工

岳飞龙 ,您好!

感谢您的评论。

COMSOL的偏导可以通过d()算子来书写,例如d(u,phi)表示的就是u在phi方向上的一阶偏导,那么二阶偏导的表达式可以写成d(d(u,phi),phi) 。

谢谢!

嘉林 张
嘉林 张
2019-12-04

您好,请问在你验证的这个模型中,这个材料的导热系数kappa是个常数吗,还是会随温度升高而升高。如果是个(各项同性)变量又该怎么输入呢?

屹磊 金
屹磊 金
2019-12-05 COMSOL 员工

张嘉林,您好!

感谢您的评论。

这个验证模型旨在讲述如何使用自定义方程来模拟二维轴对称模型,因此为了突出目的,导热系数kappa设置为常数,当然在COMSOL中导热系数是可以设置为变量的,直接在材料属性窗口输入表达式即可,举例来说,如果导热系数为温度的函数kappa=f(T),则可以直接在导热系数的输入栏中输入f(T)。

嘉林 张
嘉林 张
2019-12-05

金屹磊,您好!
十分感谢您的回复!
这个模型的关于导热系数为变量的这个输入问题我已经理解,同样的,如果我把kappa换为另一个变量h,它不是一个材料参数,那么我应该怎么输入进去呢,因为我的理解是前面的哈密顿算子的运算法则,哈密顿算子会对h求梯度,这样我就不知道如何输入了?还是直接输入h*ur和h*uz,由或者是(hr)*(ur)和(hz)*(uz)。但补偿项应该还是h*(ur)/r吧?(PS:在软件中哈密顿算子的运算法则和数学的运算法则应该是一样的吧?)

屹磊 金
屹磊 金
2019-12-23 COMSOL 员工

张嘉林,您好!
没太理解你的问题,如果只是不希望对变量h求梯度,可以对您的方程进行整理,使h包含在f中。

嘉林 张
嘉林 张
2019-12-27

十分感谢您的回复,谢谢!

凯明 王
凯明 王
2019-12-20

您好,对于合金相场凝固模拟,∂C/∂t=▽·D[▽C+(1-k)C/(1-φ+kφ)▽φ],C为溶质场,φ为相场,k为常数,D为关于相场变量φ的关系式,此方程如何分解,如何输入comsol,

Shengyu Jia
Shengyu Jia
2020-02-04

感谢您的讲解!有两个问题想请教一下:
1)您提到的“将物理场接口用作偏微分方程模板”,具体是什么意思呢?怎样操作可以使得我仍使用这些接口中的微分算子的张量意义呢?
2)如何使用PDE写旋度?
谢谢!

翀 施
翀 施
2020-06-23 COMSOL 员工

“将物理场接口用作偏微分方程模板”的意思,就是采用已有的“物理场接口”而不是用自定义方程的形式来建模。旋度等算子算符,都需要将其展开,写成偏导的形式,然后输入COMSOL中。

文举 彭
文举 彭
2020-06-21

您好,老师!
请问,在柱坐标系系的设置栏中,原点、纵轴、轴方向φ=0,这三项代表啥意思呀?各项的设置均需要填写x、y、z值的数据,该怎么理解呢?

翀 施
翀 施
2020-06-23 COMSOL 员工

原点在这里指的是在笛卡尔坐标系中对应的柱坐标系的原点位置,默认与笛卡尔坐标系的原点相同;纵轴即柱坐标系下的a,默认对应笛卡尔坐标系的z轴;初始方位角phi=0默认与笛卡尔坐标系的x轴一致。所以其实这三项是用于确定柱坐标系的位置,一般按照默认即可。

叶 柳
叶 柳
2021-06-03

你好!请问如果是实心的轴对称问题,r=0处,f/r无意义,散度如何计算呢?

Qihang Lin
Qihang Lin
2021-06-09 COMSOL 员工

对于一个矢量场 F 而言,散度有两种不同的定义方式。由散度的定义可知,div F 表示在某点处的单位体积内散发出来的矢量F的通量,所以 div F 描述了通量源的密度。

泽华 汪
泽华 汪
2021-07-07

你好,请问关于力学方面应力应变怎么用分量形式表示该张量。

Lei Cao
Lei Cao
2021-07-09 COMSOL 员工

泽华 汪, 您好!

感谢您的评论。
力学问题求解的是位移分量,应力应变有多种类型各有定义公式。应力应变张量本身就是矩阵形式,其中各元素即为正向和切向分量。
若使用结构力学相关接口建模,可在后处理中找到相应变量。若使用 PDE 建模需进行推导,将矩阵分解整合到方程组。

如果有进一步问题,建议您联系COMSOL的技术支持团队:
在线支持中心:www.denkrieger.com/support
Email:support@comsol.com
谢谢!

泽华 汪
泽华 汪
2021-07-11

就是pde接口的时候不太会分解,比如εij,就说是矢量,可以教一下嘛。谢谢

Lei Cao
Lei Cao
2021-07-13 COMSOL 员工

泽华 汪, 您好!

建议您先学习下 PDE 建模的基本操作方法,再结合固体力学接口中的方程以及内置变量定义。PDE 相关视频资源如下:
//www.denkrieger.com/video/equation-based-modeling-webinar-cn
//www.denkrieger.com/video-training/advanced-skills-and-tricks-3

如果有进一步问题,建议您联系COMSOL的技术支持团队:
在线支持中心:www.denkrieger.com/support
Email:support@comsol.com
谢谢!

婧 高
婧 高
2022-04-27

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

Qihang Lin
Qihang Lin
2022-04-29 COMSOL 员工

您好,文中已经写明了使用 “源” 这个功能节点添加补偿,同时在下方公式推导中我们可以发现原文将方程左边的一项移至右侧,代表将该项当作方程的源项进行处理。如果您不清楚在软件中如何进行具体设置,可以结合文章第一段提到的此博客使用的案例模型进行学习://www.denkrieger.com/model/steady-state-2d-axisymmetric-heat-transfer-with-conduction-453

远 张
远 张
2022-06-15

你好,我在上边将组件中坐标系设定为柱坐标系(r,phi,z)后,pde一般形式微分方程中还是∇.Γ=(∂/∂r,∂/∂phi,∂/∂z),怎么变成我想要的∇.Γ=(∂/∂r,1/r*∂/∂phi,∂/∂z)呢?可以解答一下吗,非常感谢

Qihang Lin
Qihang Lin
2022-06-22 COMSOL 员工

这部分需要您进行手动的坐标转换。具体您可以看高级培训第三天的PDE建模坐标系换算部分://www.denkrieger.com/video-training/advanced-skills-and-tricks-1

秋莉 牛
秋莉 牛
2023-07-13

您好,柱坐标系中第一轴是r,第二轴是phi,,第三轴是a,那么如何进行量纲化成我想要的量纲一化坐标系 ( θ,η,ζ)呢?其中θ = x /R,η = y /h,ζ = z/( L /2)。请帮我解答一下疑惑,非常谢谢您。

Anran Wei
Anran Wei
2023-07-18 COMSOL 员工

您好,若您是想把柱坐标系 (r, phi, a) 下的方程转换成无量纲坐标系 (r’, phi’, a’) 下的形式,其中 r’=r/R, phi’=phi/h, a’=a/(L/2), 则无量纲坐标系下的PDE方程中 ∂/∂r’, ∂/∂phi’, ∂/∂a’ 可表示为 (1/R)*∂/∂r, (1/h)*∂/∂phi, (2/L)*∂/∂z。在一般形式偏微分方程接口的设置中,需要在输入Γ的表达式时对应添加前面的系数。

浏览 COMSOL 博客
Baidu
map