
假设你在模拟这样的情况:载荷移动时横跨不同的网格单元和边界。这时你只希望对一部分几何边界施加边界条件或只在特定条件下施加。在本篇博客文章中,我们将讨论如何利用 COMSOL Multiphysics 灵活处理这类特殊情况。
边界条件的分类
对偏微分方程进行数学处理时,会遇到狄氏、纽曼和洛平这三类边界条件。狄氏条件用于指定要求解的变量;纽曼条件用于指定通量,也就是因变量的梯度;洛平条件结合了前两类边界条件,用于指定变量及其梯度之间的关系。
下表列举了这三类边界条件用于不同物理场的几个示例以及对应的物理解释。
物理场 | 狄氏条件 | 纽曼条件 | 洛平条件 |
---|---|---|---|
固体力学 | 位移 | 牵引(应力) | 弹簧 |
传热 | 温度 | 热通量 | 对流 |
压力声学 | 声压 | 法向加速度 | 阻抗 |
电流 | 固定电势 | 固定电流 | 阻抗 |
在使用有限元方法时,这三类边界条件会对求解问题的矩阵结构产生不同的影响。
纽曼条件
纽曼条件是“载荷”,出现在方程组右侧。在 COMSOL Multiphysics 的方程视图 中,这类边界条件显示为弱贡献。纽曼条件纯粹是方程组右侧附加的贡献,因此可以包含以下变量的任何函数:时间、坐标或参数值。
我们来看一个传热问题,其中半径为 R 的圆形热源沿 x 轴方向以速度 v 移动。其强度分布呈抛物线形,峰值为 q_0。此载荷的数学描述可以是
很明显,移动载荷不可能有域边界,甚至不可能存在一个始终适合载荷分布的网格。
我们可以在该表达式中直接输入载荷分布本身。因为有两处会用到径向坐标变量 r ,所以将其定义成变量是一个好方法。移动热源的完整输入如下图所示。
描述移动热源的参数。
描述移动热源的局部径向坐标相对于当前中心的变量。
输入热通量。
下方动画显示的是采用上述数据进行瞬态仿真后的结果。假定该问题关于 yz 平面对称,则载荷实际施加到了半圆形的移动热源上。
热源沿条块移动时温度分布的动画。
狄氏条件
当给定狄氏条件时,因变量就指定了,所以无须对其求解。我们可以从问题中删除这一类自由度方程。因此狄氏条件会改变刚度矩阵的结构。在 COMSOL Multiphysics 的方程视图 中,这类条件显示为约束。
假定要将移动点的温度指定为刚好 450 K,这或许有点刻意,但是能表现出纽曼条件和狄氏条件之间的一个重要区别。假如要添加一个温度 节点并输入类似表达式( if(r < R,450[K],0)
),这意味着将热点不会覆盖的那部分边界的温度设定为绝对零度。
不过,我们的目的是在热点之外停用狄氏条件。为此可以使用一个小窍门:输入 if(r < R,450[K],ht.Tvar)
作为指定值,就能获得所期望的停用(如下方动画所示)。
含条件限制的狄氏条件设置。
温度指定的热点沿条块移动时的温度分布动画。
为了理解其中的工作原理,我们启用方程视图,查看狄氏条件(本示例中为指定温度)的实现:
启用 方程视图。
方程视图中的 温度节点 。
约束公式为 ht.T0-ht.Tvar
,暗指 ht.T0-ht.Tvar = 0
。第一项是作为输入的指定温度。第二项是转换成变量的温度自由度,这一约束使温度等于给定值,不过给定值恰好是字符串 ht.Tvar
的情况除外。如果发生这种情况,符号代数在汇编期间会将表达式简化为 ht.Tvar-ht.Tvar
,并进一步计算为零。约束表达式为 0
,即没有约束。
弱约束
在 COMSOL Multiphysics 中,狄氏条件实际上有两种可能的实现方式:默认情况下为上文所提到的逐点约束,还可以使用弱约束。在弱约束中,我们要添加方程,而不是移除方程。热通量添加为额外的自由度(拉格朗日乘子),以强制指定温度值。
只要对上面的温度表达式稍作更改,就可以使用本质上相同的窍门使弱约束含条件限制。要使用弱约束,须先启用高级物理场接口选项。
启用 高级物理场接口选项。
在物理场接口内对某个节点选中弱约束后,拉格朗日乘子上就会添加额外的自由度。这个示例中,该自由度名为 T_lm
。
如果应用了上文中的温度表达式,那么在停用狄氏条件的边界部分,此额外自由度不会得到任何刚度矩阵的贡献,刚度矩阵由此变为奇异矩阵。为了避免这种情形,我们将 if(r < R,450[K],ht.Tvar)
更改为 if(r < R,450[K],ht.Tvar-T_lm*1e-2)
。不同的模型和物理场, T_lm
所用的乘子也不同,因此可能需要进行一些调整。
至于其中的工作原理,就像教科书所讲的那样,“留给读者作为练习”。
使用弱约束时含条件限制的狄氏条件设置。
洛平条件
洛平条件通常都会影响刚度矩阵和方程右侧。虽然刚度矩阵的结构不会受到影响,但现有位置上会添加值。在方程视图 中,洛平条件同样显示为弱贡献。将这类条件转换为关于时间、空间和其他变量的函数,这与使用纽曼条件时的做法一致。
不过有趣的是,选择合适的值确实可以转换洛平条件,使之近似为狄氏条件或纽曼条件。如果仿真期间你希望在这两类边界条件之间切换,那么这一点十分重要。
要创建狄氏条件,需要对“刚度”指派一个高值,例如弹簧常数或传热系数。在数学术语中,这实质上是狄氏条件的罚 实现。刚度越高,自由度的指定值就越精确。但这里需要注意:刚度过高会影响刚度矩阵的数值稳定性。而在传热问题中,要选择“高”的传热系数 \alpha,起始值可以是
其中 \lambda 表示导热系数,h 表示特征单元尺寸。
用适当的材料属性(也就是固体力学中的杨氏模量)替换 \lambda,即可以在其他物理场实现相同的计算。将因子设为 1000 只是一个建议,可以替换成 104 或 105。
如果要使用对流模拟上一个示例中移动温度为 450 K 的热点,则可以采用下图中的设置。单元尺寸的内置变量 h
就应用到了表达式中。
利用对流条件指定温度。
当年刚刚开始使用有限元分析时,我有时无法在结构力学的有限元程序中指定非零位移。这一限制是由于编程越来越复杂造成的。遇到这种情况时,最佳方法是添加预变形的刚性弹簧来使用罚函数法。但是刚度不能过大;因为那时还在采用用单精度运算!
下面,我们转而讨论近似为纽曼条件。我们希望热通量与表面温度无关。在传热这一示例中,洛平条件声明向内热通量 q 为
其中 \alpha 表示传热系数,T 表示边界温度,T_{\textrm{ext}} 表示外部温度。
所以如果 T_{\textrm{ext}} 比表面预期温度高很多,那么 q \approx \alpha T_{\textrm{ext}}。则采用的策略可以是选择一个任意大的 T_{\textrm{ext}},然后计算一个合适的传热系数,如下图高亮处所示。
利用对流条件指定热通量。
在实际情况中,设计人员利用这一原理将以下固定力引入实际机械部件中:又长又软的预应力弹簧。如果弹簧的预变形远大于与弹簧连接的部件的位移,则弹簧的力几乎保持不变。
解决离散误差
当边界条件受制于布尔表达式(类似 if(r < R,...
)时,施加了此边界条件的区域边界很可能不会跟随网格单元的边,这将引入离散误差。
使用纽曼条件或洛平条件时,我们要对每个有限元执行数值积分,并计算该单元中大量离散高斯点的函数值。如果网格单元的尺寸大于载荷的几何尺寸,那么施加有载荷的确切高斯点数会显著影响总载荷。由此,施加有载荷的小块区域上应该总是存在几个单元。
载荷位置稍微移动一点就可能改变所含积分点的数量。(在实际情况中,所改变的积分点数会更多。)
因此,至少从离散点的意义上讲,要对网格节点施加狄氏条件。下图显示,在模拟指定温度为 450 K 的圆形点的移动时,某一时间点上的温度分布和热通量。在热点的前端,一个温度为 260 K 的暗影清晰可见。由于模拟的初始温度和环境温度均为 293 K,因此这一暗影不应出现。出现这一数值差异,是基于这样一个事实:并非每个单元上的所有节点都设定了狄氏条件。在狄氏条件的不连续性处,就会出现奇异点。上一篇博客文章曾探讨过这一话题。细化网格可以减少这一类问题。
下图中的绿色箭头表示指定温度后产生热通量的节点。模型中的网格密度太小,使半圆的近似变得相当困难。
指定温度的半圆周围的温度分布和热通量。
解对边界条件的依赖性
要将解包含入边界条件有多种方法。这样做往往会引入非线性,COMSOL Multiphysics 可以自动检测到这样的非线性。
我们以一个梁为例,梁的稍下方有一个支撑,其作用是在梁发生一定挠曲后阻止梁的进一步移动。在梁 接口的指定位移/旋转 节点中,设置一个含条件限制的狄氏条件可以实现这一模拟。
具有挠曲、控制支撑和分布载荷的梁。
指定梁应在挠曲为 2 厘米时停止移动的设置。
这一分析表明希望挠曲为 2 厘米时梁停止移动。载荷较小时,挠曲形状呈对称结构,而载荷较大时,梁上受额外支撑的点会停止移动。载荷达到最大时,梁上的曲率甚至会由正值变为负值。变形图可以显示这一变化,不过弯矩图中更加明显。
梁在支撑点处移动 2 厘米后停止的梁位移。
载荷不同时梁上的弯矩。
此处重点讨论的方法还相当不成熟,而且迭代解的收敛特性可能不够好。更可靠的实现方式是在支撑点处使用高度非线性弹簧,由此反作用力就是关于位移的连续可微函数。这种做法实际上类似于在固体力学 接口中实现罚接触。
在 COMSOL Multiphysics 中使用边界条件的结论和思考
对于非标准边界条件的指定,COMSOL Multiphysics 提供了强大的机制。今天,我们探讨了几个使用此类条件的案例。
如果你对分析含有移动载荷的模型感兴趣,请看一下“App 库”中的移动载荷教程模型。
如果你在建模过程中对如何指定非标准边界条件有其他疑问,请立即联系我们。
评论 (22)
恭正 陈
2018-11-13comsol怎么模拟一个正方形的不稳定性运动?
恭正 陈
2018-11-13请教下COMSOL如何模拟正方形的不稳定运动?
Hehua Xu
2018-11-15在求解地球上的区域热结构的问题中,通常是已知表层的温度,同时通过测试也可以知道表明的热流(热通量),而某个深度的温度和热流都未知。在选定求解区域后,最好的办法是在上边界同时加载温度和热流边界,但软件中无法实施,不知专家有什么建议。谢谢。
Yong Jiang
2019-03-14不同的物理场,如何选择合适的节点?
Haili Wang
2019-04-15Hehua ,您好!
感谢您的评论。
对于您的问题,建议您联系系我们的技术支持团队:
在线支持中心:cn.comsol.com/support
Email: support@comsol.com
谢谢!
Haili Wang
2019-04-15Yong Jiang ,您好!
感谢您的评论。
您的问题有点过于笼统,建议您将问题细化之后,联系COMSOL的技术支持团队:
在线支持中心:cn.comsol.com/support
Email: support@comsol.com
谢谢!
义奎 温
2019-06-06关于动载荷的施加和变化的电场是否计算没有太大的差别呢
Jiawei Kao
2020-01-29您好!comsol在固体力学中,建立了较为复杂的几何模型后,在选择其中一些边界施加相应条件时,发现这些边界的相应编号不是按着建模顺序的前后进行赋予的,请问这个有没有相关的详细说明?
Jing Yuan
2020-07-12是否有人可以解释一下,这里面洛平边界条件是否对应的是 k*dT/dx=T0-T 这种边界条件?我的问题是,如何在Comsol中设置形如 k.dT/dx=T0-T这样的边界条件呢?我尝试用convective heat flux, q0=h(Text-T), 可是计算结果q0和conductive heat flux (k*dT/dx) 并不相等,这是为什么呢?谢谢
翀 施
2020-07-13 COMSOL 员工洛平(罗宾)边界条件在COMSOL的传热模块中对应的条件应该是-k*dT/dx=h*(T0-T) ,需要添加一个负号,另外这个表达式是相对与1D模型而言,如果是二维甚至是三维还存在T关于y和z的梯度项。COMSOL中的convective heat flux, q0=h(Text-T),就是一个现成的洛平(罗宾)边界条件。
hang zhang
2020-09-21请问我的激光光源是条状的不是圆形的,我该怎么定义条状的激光光源
Ying (Grace) Xu
2020-09-24感谢您的评论。
定义条状的激光光源需要有描述条状光源强度分布的数学表达式,然后替换博客中圆形强度分布的表达式。
如还有其他问题联系COMSOL的技术支持团队:
在线支持中心:cn.comsol.com/support
Email: support@comsol.com
谢谢!
shelly 赖
2021-05-10请问用comsol如何得到方形板四周受均布压力时的位移云图。总是会提示相对误差(7.4)大于相对容差。comsol中有像ansys中的弱弹簧平衡微小力的设置吗?
Yang Zhang
2021-05-11 COMSOL 员工您这是由于模型欠约束导致的报错,您可以在边或域上添加“弹簧基础(spring foundation)”实现其平衡。
shelly 赖
2021-05-13你好,请问如果在四周受压的方形板中间挖一个洞,并给定洞的边界上的点的位移,计算结果显示方形性板四个角的位移极大且超出合理范围。如果只指定洞的边界的位移,计算结果合理。若是加上四周的面力,四个角的位移极大。请问这个要怎么处理?
Lei Cao
2021-05-19 COMSOL 员工shelly 赖, 您好!
感谢您的评论。
仅凭描述无法准确判断您的问题,目前猜测可能是施加的力的量级存在问题,也就是设置的力过大。或是对于大变形问题未考虑几何非线性,计算都按照线弹性材料模型进行考虑。
建议您将问题和相关模型发送至COMSOL的技术支持团队:
在线支持中心:cn.comsol.com/support
Email: support@comsol.com
谢谢!
柳春 杨
2021-05-17我想请教一下边界条件怎么设置不会反弹超声波
Lei Cao
2021-05-19 COMSOL 员工柳春 杨, 您好!
感谢您的评论。
您指的应该是声波的反射。在材料声阻抗变化的边界上都会产生一定的反射情况。针对不同的建模情况有不同的设置方式,可尝试辐射边界、阻抗边界、完美匹配层、低反射边界(固体瞬态)等方式。
如果有进一步问题,建议您联系COMSOL的技术支持团队:
在线支持中心:cn.comsol.com/support
Email: support@comsol.com
谢谢!
Yuxi Zhang
2021-09-24您好,请教下:无限大的空间,建模时只截取有限空间(比如一立方体)进行仿真,那如何在立方体的边界设置连续边界条件呢?及立方体内的物理量=立方体外的物理量。 谢谢
Liwen Yang
2021-09-27 COMSOL 员工您好,截取无限域的有限空间进行有限元仿真,对于截断区域有三种处理方法:无限元、完美匹配层、吸收层。对于不同类型的问题需要选择不同的处理方式,详细内容可以参考:http://cn.comsol.com/support/knowledgebase/1272 。
xiaoling chen
2023-03-20高压可压缩流体的入口条件如何同时设置为压力和速度边界条件?
越 赵
2023-03-22 COMSOL 员工您可以使用高马赫数流接口中的入口边界条件来设置,压力体现在入口静压上,速度可以转换为马赫数进行边界设置。