如何设置特殊边界条件

Henrik Sönnerlind 2016年 6月 27日

假设你在模拟这样的情况:载荷移动时横跨不同的网格单元和边界。这时你只希望对一部分几何边界施加边界条件或只在特定条件下施加。在本篇博客文章中,我们将讨论如何利用 COMSOL Multiphysics 灵活处理这类特殊情况。

边界条件的分类

对偏微分方程进行数学处理时,会遇到狄氏纽曼洛平这三类边界条件。狄氏条件用于指定要求解的变量;纽曼条件用于指定通量,也就是因变量的梯度;洛平条件结合了前两类边界条件,用于指定变量及其梯度之间的关系。

下表列举了这三类边界条件用于不同物理场的几个示例以及对应的物理解释。

物理场 狄氏条件 纽曼条件 洛平条件
固体力学 位移 牵引(应力) 弹簧
传热 温度 热通量 对流
压力声学 声压 法向加速度 阻抗
电流 固定电势 固定电流 阻抗

在使用有限元方法时,这三类边界条件会对求解问题的矩阵结构产生不同的影响。

纽曼条件

纽曼条件是“载荷”,出现在方程组右侧。在 COMSOL Multiphysics 的方程视图 中,这类边界条件显示为弱贡献。纽曼条件纯粹是方程组右侧附加的贡献,因此可以包含以下变量的任何函数:时间、坐标或参数值。

我们来看一个传热问题,其中半径为 R 的圆形热源沿 x 轴方向以速度 v 移动。其强度分布呈抛物线形,峰值为 q_0。此载荷的数学描述可以是

q(x,y,t)=q_0\left(1-\left(\frac{r}{R}\right)^2\right), \quad r=\sqrt{(x-vt)^2+y^2}, \quad r < R

很明显,移动载荷不可能有域边界,甚至不可能存在一个始终适合载荷分布的网格。

我们可以在该表达式中直接输入载荷分布本身。因为有两处会用到径向坐标变量 r ,所以将其定义成变量是一个好方法。移动热源的完整输入如下图所示。

Heat Source Parameters 如何设置特殊边界条件
描述移动热源的参数。
Local radial coordinate variable 如何设置特殊边界条件
描述移动热源的局部径向坐标相对于当前中心的变量。
Heat flux input 如何设置特殊边界条件
输入热通量。

下方动画显示的是采用上述数据进行瞬态仿真后的结果。假定该问题关于 yz 平面对称,则载荷实际施加到了半圆形的移动热源上。

 

热源沿条块移动时温度分布的动画。

狄氏条件

当给定狄氏条件时,因变量就指定了,所以无须对其求解。我们可以从问题中删除这一类自由度方程。因此狄氏条件会改变刚度矩阵的结构。在 COMSOL Multiphysics 的方程视图 中,这类条件显示为约束。

假定要将移动点的温度指定为刚好 450 K,这或许有点刻意,但是能表现出纽曼条件和狄氏条件之间的一个重要区别。假如要添加一个温度 节点并输入类似表达式( if(r < R,450[K],0)),这意味着将热点不会覆盖的那部分边界的温度设定为绝对零度。

不过,我们的目的是在热点之外停用狄氏条件。为此可以使用一个小窍门:输入 if(r < R,450[K],ht.Tvar) 作为指定值,就能获得所期望的停用(如下方动画所示)。

Conditional Dirichlet condition settings 如何设置特殊边界条件
含条件限制的狄氏条件设置。

 

温度指定的热点沿条块移动时的温度分布动画。

为了理解其中的工作原理,我们启用方程视图,查看狄氏条件(本示例中为指定温度)的实现:

Enabling Equation View 如何设置特殊边界条件
启用 方程视图
Equation View for Temperature node 如何设置特殊边界条件
方程视图中的 温度节点 。

约束公式为 ht.T0-ht.Tvar,暗指 ht.T0-ht.Tvar = 0。第一项是作为输入的指定温度。第二项是转换成变量的温度自由度,这一约束使温度等于给定值,不过给定值恰好是字符串 ht.Tvar 的情况除外。如果发生这种情况,符号代数在汇编期间会将表达式简化为 ht.Tvar-ht.Tvar,并进一步计算为零。约束表达式为 0,即没有约束。

弱约束

在 COMSOL Multiphysics 中,狄氏条件实际上有两种可能的实现方式:默认情况下为上文所提到的逐点约束,还可以使用弱约束。在弱约束中,我们要添加方程,而不是移除方程。热通量添加为额外的自由度(拉格朗日乘子),以强制指定温度值。

只要对上面的温度表达式稍作更改,就可以使用本质上相同的窍门使弱约束含条件限制。要使用弱约束,须先启用高级物理场接口选项

Enabling Advanced Physics Options 如何设置特殊边界条件
启用 高级物理场接口选项

在物理场接口内对某个节点选中弱约束后,拉格朗日乘子上就会添加额外的自由度。这个示例中,该自由度名为 T_lm

如果应用了上文中的温度表达式,那么在停用狄氏条件的边界部分,此额外自由度不会得到任何刚度矩阵的贡献,刚度矩阵由此变为奇异矩阵。为了避免这种情形,我们将 if(r < R,450[K],ht.Tvar) 更改为 if(r < R,450[K],ht.Tvar-T_lm*1e-2)。不同的模型和物理场, T_lm 所用的乘子也不同,因此可能需要进行一些调整。

至于其中的工作原理,就像教科书所讲的那样,“留给读者作为练习”。

conditional Dirichlet condition when using weak constraints 如何设置特殊边界条件
使用弱约束时含条件限制的狄氏条件设置。

洛平条件

洛平条件通常都会影响刚度矩阵和方程右侧。虽然刚度矩阵的结构不会受到影响,但现有位置上会添加值。在方程视图 中,洛平条件同样显示为弱贡献。将这类条件转换为关于时间、空间和其他变量的函数,这与使用纽曼条件时的做法一致。

不过有趣的是,选择合适的值确实可以转换洛平条件,使之近似为狄氏条件或纽曼条件。如果仿真期间你希望在这两类边界条件之间切换,那么这一点十分重要。

要创建狄氏条件,需要对“刚度”指派一个高值,例如弹簧常数或传热系数。在数学术语中,这实质上是狄氏条件的实现。刚度越高,自由度的指定值就越精确。但这里需要注意:刚度过高会影响刚度矩阵的数值稳定性。而在传热问题中,要选择“高”的传热系数 \alpha,起始值可以是

\alpha=1000 \frac{\lambda}{h}

其中 \lambda 表示导热系数,h 表示特征单元尺寸。

用适当的材料属性(也就是固体力学中的杨氏模量)替换 \lambda,即可以在其他物理场实现相同的计算。将因子设为 1000 只是一个建议,可以替换成 104 或 105

如果要使用对流模拟上一个示例中移动温度为 450 K 的热点,则可以采用下图中的设置。单元尺寸的内置变量 h 就应用到了表达式中。

Convection condition for temperature 如何设置特殊边界条件
利用对流条件指定温度。

当年刚刚开始使用有限元分析时,我有时无法在结构力学的有限元程序中指定非零位移。这一限制是由于编程越来越复杂造成的。遇到这种情况时,最佳方法是添加预变形的刚性弹簧来使用罚函数法。但是刚度不能过大;因为那时还在采用用单精度运算!

下面,我们转而讨论近似为纽曼条件。我们希望热通量与表面温度无关。在传热这一示例中,洛平条件声明向内热通量 q

q=\alpha(T_{\textrm{ext}}-T)

其中 \alpha 表示传热系数,T 表示边界温度,T_{\textrm{ext}} 表示外部温度。

所以如果 T_{\textrm{ext}} 比表面预期温度高很多,那么 q \approx \alpha T_{\textrm{ext}}。则采用的策略可以是选择一个任意大的 T_{\textrm{ext}},然后计算一个合适的传热系数,如下图高亮处所示。

Convection Condition for heat flux 如何设置特殊边界条件
利用对流条件指定热通量。

在实际情况中,设计人员利用这一原理将以下固定力引入实际机械部件中:又长又软的预应力弹簧。如果弹簧的预变形远大于与弹簧连接的部件的位移,则弹簧的力几乎保持不变。

解决离散误差

当边界条件受制于布尔表达式(类似 if(r < R,...)时,施加了此边界条件的区域边界很可能不会跟随网格单元的边,这将引入离散误差。

使用纽曼条件或洛平条件时,我们要对每个有限元执行数值积分,并计算该单元中大量离散高斯点的函数值。如果网格单元的尺寸大于载荷的几何尺寸,那么施加有载荷的确切高斯点数会显著影响总载荷。由此,施加有载荷的小块区域上应该总是存在几个单元。

Change in number of contributing integration points 如何设置特殊边界条件
载荷位置稍微移动一点就可能改变所含积分点的数量。(在实际情况中,所改变的积分点数会更多。)

因此,至少从离散点的意义上讲,要对网格节点施加狄氏条件。下图显示,在模拟指定温度为 450 K 的圆形点的移动时,某一时间点上的温度分布和热通量。在热点的前端,一个温度为 260 K 的暗影清晰可见。由于模拟的初始温度和环境温度均为 293 K,因此这一暗影不应出现。出现这一数值差异,是基于这样一个事实:并非每个单元上的所有节点都设定了狄氏条件。在狄氏条件的不连续性处,就会出现奇异点。上一篇博客文章曾探讨过这一话题。细化网格可以减少这一类问题。

下图中的绿色箭头表示指定温度后产生热通量的节点。模型中的网格密度太小,使半圆的近似变得相当困难。

Temperature distribution and heat flux around semicircular temperature 如何设置特殊边界条件
指定温度的半圆周围的温度分布和热通量。

解对边界条件的依赖性

要将解包含入边界条件有多种方法。这样做往往会引入非线性,COMSOL Multiphysics 可以自动检测到这样的非线性。

我们以一个梁为例,梁的稍下方有一个支撑,其作用是在梁发生一定挠曲后阻止梁的进一步移动。在接口的指定位移/旋转 节点中,设置一个含条件限制的狄氏条件可以实现这一模拟。

Beam with a deflection 如何设置特殊边界条件
具有挠曲、控制支撑和分布载荷的梁。
Prescribing location for stop of beam deflection 如何设置特殊边界条件
指定梁应在挠曲为 2 厘米时停止移动的设置。

这一分析表明希望挠曲为 2 厘米时梁停止移动。载荷较小时,挠曲形状呈对称结构,而载荷较大时,梁上受额外支撑的点会停止移动。载荷达到最大时,梁上的曲率甚至会由正值变为负值。变形图可以显示这一变化,不过弯矩图中更加明显。

Beam displacement at 2 cm 如何设置特殊边界条件
梁在支撑点处移动 2 厘米后停止的梁位移。
Bending moment of beam for various loads 如何设置特殊边界条件
载荷不同时梁上的弯矩。

此处重点讨论的方法还相当不成熟,而且迭代解的收敛特性可能不够好。更可靠的实现方式是在支撑点处使用高度非线性弹簧,由此反作用力就是关于位移的连续可微函数。这种做法实际上类似于在固体力学 接口中实现罚接触。

在 COMSOL Multiphysics 中使用边界条件的结论和思考

对于非标准边界条件的指定,COMSOL Multiphysics 提供了强大的机制。今天,我们探讨了几个使用此类条件的案例。

如果你对分析含有移动载荷的模型感兴趣,请看一下“App 库”中的移动载荷教程模型

如果你在建模过程中对如何指定非标准边界条件有其他疑问,请立即联系我们


博客分类

博客标签

技术资料
加载评论……