COMSOL建模指南

翻译自COMSOL的用户手册。 对使用其他仿真软件也很有指导价值。

 

用户可以调整求解设置、使用对称或其他模型属性来获得模拟解，如果不行的话，可以中断求解来恢复部分解。

提纲：

1.       使用对称性

2.       有效的内存管理

3.       选择单元类型

4.       分析模型收敛性和准确性

5.       求解非线性方程时实现收敛

6.       避免过快的瞬变

7.       物理相关检查和指南

8.       非物理值

 

正文：

1.       使用对称性

建模时使用对称性可以减少一半或以上的计算量，对于求解大的问题这是一种有效的手段。它适用于包含对称几何结构或模型假定的情况下。对称有轴对称和对称和反对称平面或线：

a)         轴对称：通常是圆柱体或相似的3D几何。如果几何是轴对称的，那么它只能在轴向和径向变化，而在周向保持不变。

b)         对称和反对称平面或线：在2D和3D模型中比较常见。对称：分割线或平面两边是相同的。对于标量场，穿过对称线的法向流量为零。在结构力学中，对称条件是不同的。反对称：分割线或平面两边是反向平衡的。对于标量场，穿过对称平面或线的独立变量为零。结构力学使用另外的反对称。很多模型界面都支持对称边界条件。

为了使用对称平面或线，所有的几何、材料性能、边界条件都必须是对称的，任何载荷或源项也必须是对称或反对称的。然后就可以建立一个对称模型，它可以是整个几何体的1/2，1/4，1/8，然后施加合适的对称边界条件。

 

2.       有效的内存管理

特别在3D建模时，巨大的内存消耗必须采取额外的预防措施。首先，检查是否选择了线性迭代求解器。通常，用户不需要考虑选择什么求解器，因为模型界面会自动挑选合适的求解器。但是在一些情况下，用户还是需要对求解设置做适当的调整。

估算模型消耗的内存

当内存分配不足时COMSOL会提示”out-of-memory”。虽然内存大小和数组所占用的相同就足够了，但是由于内存碎片的存在，需要更大的内存。

那么求解一个特定的模型需要多大的内存，则有下面几个因素决定：

a)         节点数目

b)         独立和不独立变量数

c)         单元阶数

d)         矩阵的稀疏模式。稀疏模式取决于几何和网格的形状。例如，一个扩展的椭球比球的矩阵更稀疏。

MUMPS和PARDISO out-of-core求解器可以在内存不够时，借用磁盘存储来求解大的模型。用户可以通过观察右下角的内存使用情况来了解当前物理内存和虚拟内存使用情况。

建立一个高效内存使用的几何

首先在处理大的模型时，应该尽量减小求解几何。通常使用对称会使模型减少到原来的1/2，1/4甚至是1/8。虽然内存消耗并不是线性的，而是多项式(Cnk,k>1)，那么这意味着当模型减小一半时，内存消耗会减少一半以上。而其他的途径是：

 

a)         避免使用那些不必要的小几何体，用贝塞尔曲线代替多边形轮廓；

b)         尽量使用线性单元；

c)         保证划分的网格质量高。高质量的网格对线性迭代求解器尤为重要，这样会提高迭代的收敛性；

d)         避免几何存在尖的、狭窄的转角。网格单元会在这些地方变得很薄，使得邻近的网格质量下降。尖角同时会造成非物理现象和巨大的应力集中。

 

3.       选择单元类型

大部分模型界面默认选择一阶或二阶Lagrange单元或形函数。二阶及以上的单元类型会在网格单元的中点和内部节点上增加额外的自由度。自由度的增加会提高解的准确性，但是也会由于减小了离散系统的稀疏程度而加大了内存消耗。对于很多的应用，如结构力学中的应力分析，需要使用二阶单元来提高准确度，这是由于一阶单元的空间差分为常数。

COMSOL推荐用户使用默认的单元类型，而对于某些计算可能需要减低单元阶数，但是用户必须注意那些重要的量被求解了。

 

4.       分析模型收敛性和准确度

有限元法可以精确的获得如应力集中模型解的局部变化。在有些情况下，用户可以与参考书、测量结果或其他的原始数据做比较。而COMSOL模型库中也包含了一些计算好的结果或理论解。通过查看这些基准模型，来检查结果的准确性。

如果一个模型不能通过其他方法来验证正确性，那么只能通过收敛性测试来检查网格密度是否合适。用户通过加密网格重新计算来检查解是否可以稳定的收敛。如果加密网格后解的结果发生改变，那么解依赖于网格，该模型需要更密的网格。用户可以使用自适应网格，它基于误差准则来加密那些误差较大的区域网格。

为了收敛性，必须避免几何中的奇异点。

 

5.       在求解非线性方程时实现收敛

非线性问题通常很难被求解。在很多情况下，它的解不唯一。COMSOL使用Newton迭代法来求解非线性PDE方程组。这种方法对初始值敏感。如果初始条件偏离期望值，那么可能就无法收敛，即便该方程组使用其他初值时可能很容易收敛。

用户可以通过下面的方法来提高得到难解的非线性方程组的相关解：

a)         提供最可能的初始值；

b)         顺序求解单一方程，当获取到好的初始值后再进行全耦合计算；

c)         保证边界条件与初始值一致，相邻的边界有一个兼容条件从而不产生奇异值；

d)         加密那些梯度大的区域网格；

e)         对于对流问题，引入人工扩散来提高数值计算性能。大多数流动和化学传质建模界面将引用人工扩散作为默认设置；

f)          当一个解的成分趋于零时应采用缩放因子。在这些情况下，自动缩放可能并不能工作。

g)         将一个稳态非线性问题转化为瞬态问题。通过这种方法可以使问题易于收敛。当瞬态计算趋于稳态时，再去求解稳态问题；

h)         采用参数化扫描求解器，调整一个材料的属性或调整PDE系数的起始值从而减小方程组的非线性。用这种方法，可以得到一系列非线性解，使得非线性问题容易计算。

 

6.       避免过快的瞬变

在求解不一致初始条件的瞬态问题时，或着其边界条件、源项在一段时间内发生巨大的变化，那么求解系统引入了很强的瞬态信号。这时时间步进算法会选用很小的步长来求解这个变化，因此求解时间变得很长，有的时候求解会终止。稳态问题则可能会遇上网格分辨率的问题，诸如无线通量问题造成的过冲或下冲。

 

除非你想知道这些瞬变的细节，一般以可以获得一致解达到稳态的条件开始。然后再改成模型在真实条件下一个时间段内的边界值、源项或驱动流。

在多数情况下，可以通过一个平滑函数来引入源项。使得阶跃或转变不那么突然。

 

7.       物理相关检查和指南

 

流体流动和传质现象

a)         如果所有的边界条件都不含压力条件，那么需要在流体域内某个点上指定压强。如果不指定它，求解就会缺乏约束从而难以收敛；

b)         确保网格足够细，例如在狭窄的通道内至少保证有4-6个网格单元。

c)         确保边界条件和初始条件与瞬态问题相吻合。例如，在流体域内初始速度为零时，就不要使用全速度，而是采用一个平滑或斜坡函数将初始速度缓缓的加上去，这样才能使流场的初始值与边界条件相匹配，最终达到全速度。

d)         对于流动模型，需要通过Reynolds数来判断是层流还是湍流。当流动为湍流时，应采用湍流模型。

 

声学、结构、电磁波传播

对于描述波传播的模型，需要考虑它的时间和空间性。这意味着，最大的网格单元约提供每个波长下10个线性或5个两阶单元。对于瞬态问题，固定时间步长需要足够小。

 

结构力学

a)         确保模型是全约束的。约束模型以避免所有的刚性运动，对于3D实体力学模型需要6个约束，即3个平移和3个转动。否则模型就没有定义好，计算无法收敛。不可能把6个约束条件加载到一个点上，在它上面最多只能限制3个平移自由度。对于一个3D实体力学模型，可以使用“3-2-1准则”来约束3个自由度在一个点（一个固定约束）上，2个自由度在另外一个点上，1个自由度在第三个点上。可以选择三个分开的点（顶点）来操作。然后在三个方向上固定第一个点。限制第二个点在两个正交的方向上，并确保两点间的连线没有变形的限制。最后约束第三个点在由三个点组成的平面法线方向上。为了验证这种方法，模型一个温度载荷下可以膨胀或收缩，整个过程中应力没有集中。2D模型则是在一个点上约束2个平动自由度，在另外一个点上约束转动自由度。

b)         在变形量小的时候考虑使用线弹性材料模型，否则使用非线性材料模型；

c)         避免有尖的转角，否则会带来非物理解或带来无限大的应力集中。

 

8.       非物理解

哪里和为什么会出现非物理解？

在一些模型中，由于数值计算会带来一些少量的非物理解。如：

a)         质量扩散中的负浓度；

b)         瞬态热传导中温度略大于初始条件；

c)         结构力学中在没荷载的方向上出现少量反作用力；

d)         接触分析中少量负间隙；

e)         少量负的有效塑性应变；

f)          实体力学中理想塑性材料的应力大于屈服应力。

 

一些原因：

a)         数值噪声造成的。当一个独立变量的值趋于零，那么数值噪声的影响就相对比较大，从而引起了少量的非物理的负值。

b)         插值和外推造成的。例如，结果在有限元中的一个积分点（Gauss点）对于一个弹塑性材料是正确的（在一定的误差内），但是当外推到单元边界上时变成了非物理值。

c)         模型不连续造成的。例如，由于初值不连续产生负浓度。对于对流模型，沿着某个边界初值为零，例如，物理插值是一个初始尖锐，远离边界则是逐渐扩散的，然而，有限元模型默认的形函数（二阶Lagrange单元）只允许采用连续的函数。因此程序会在迭代前首先修改这些不连续的初始条件，这通常会在起始时间在解中引入一些“小坑”。如下图所示，浓度在t=0时出现负值：

 

d)         缺少足够的网格分辨率。如在高的对流区域（高的Peclet数）和大的反应项或快速反应动力学中（高的Damkohler数）会产生负浓度从而遇到收敛性问题。

e)         不合理的物理模型。例如对于质量输运，只能在高浓度时使用定常减少的反应项，当浓度接近零时，反应项继续消耗组分，从而产生负浓度。

 

避免非物理值

a)         在一些情况下，也可以对独立变量使用基线，从而在独立变量接近零时数值噪音不会影响求解。例如加任意一个基线溶度使浓度大于0。这种缩放方法对于某些情况无效，例如，一个反应项它与浓度有关。

b)         避免使用不连续值，用合适的平滑阶跃函数代替。

c)         使用对数变量可以消除网格分辨率和负值问题。这是因为线性变换的网格有时无法描述浓度变化的指数行为。另外，使用对数变量可以保证求解过程中原变量不会出现负值。

d)         在绘图时避免显示数值噪音造成的非物理值。例如使用c*(c>0)来代替c，这样在c小于零的部分都强制等于0。同时还可以通过调整绘图的范围来去掉它们。

e)         研究细化网格时非物理值是否变好或变坏。如果变好，那么就可以继续加密网格。如果变坏，那么就应检查物理模型是否正确。