稳定方法 stabilization

原文参考https://www.comsol.com/blogs/understanding-stabilization-methods/，部分语句为了通顺有改动。

当输运过程由对流而不是扩散主导时，大多数的数值模拟方法包括有限元，有限体积，有限差分，就需要稳定方法。对于有限元方法来说，稳定方法即添加一个小的人工扩散，来增强健壮性并提高计算速度。本文旨在帮助理解稳定方法对数值模型影响，同时给出一种无需稳定方法的替代的数值方法。

什么时候需要稳定方法？

COMSOL采用有限元方法进行数值求解。有限元方法直接适用于固体力学和扩散控制的输运过程问题的求解，而对于对流主导的问题，往往会出现计算不稳定，比如数值解的震荡。不过，COMSOL会对这一过程自动的引入稳定方法。

对流扩散方程的常见形式：

beta表示速度；c表示扩散系数；h表示网格尺寸。

Peclet（Pe）数，衡量对流和扩散的大小关系。当Pe远超过 1 时，需要引入稳定方法。

例子

官网没给算例，可能是比较简单，这里为了校验文章的说法，我自己算了一个，COMSOL5.2：链接: http://pan.baidu.com/s/1qYsyJt6 密码: iiby

为了更好的理解，这里举一个关于时间相关质量输运的例子：3m 长的圆管，柱塞流 1m/s，流体为水，并且有一种化学物质溶解在水中，其初始浓度满足 1(x<0);0(x>0)，如图。红色虚线为 1s 时刻的解析解，下面尝试给出数值解。

作为 1D 问题，使用Transport of Dilute Species模块进行求解。设置，速度 1m/s；扩散系数 1e-9 m^2/s；进口浓度 1 mol/m^3。网格尺寸取为0.05，则 Pe=25e6，远大于1。如果不采用稳定方法，将得到下面的解

从图中可以清楚的看到，震荡不仅出现在浓度梯度大的地方，而是遍布整个计算域。出现解的震荡需要满足Pe>>1，以及下面三种情况中的任意一种：
1. 空间分布的初值，非网格求解结果；
2. Dirichlet边界条件，包含较大的梯度，形成边界层；
3. 初始扩散项很小，接近非定常源项或者非定常Dirichlet边界条件。
_PS. 不太理解这一段。

当然可以通过缩小网格使得Pe<1，这样网格会非常多，超过 1 亿。

COMSOL中的稳定方法

所有的输运模块，比如heat transfer，fluid flow，以及species transport都会自动
使用稳定方法。如下图操作可以使得COMSOL的稳定方法stabilization可视化。


接下来将展示所有的稳定方法，并给出相应的结果。总的来说，所有的稳定方法都是通过添加额外的扩散项，来减小Pe数。

Inconsistent stabilization: isotropic diffusion

最简单的方法就是定义一个人工扩散系数 c_art，如下：

并且添加它到物理扩散系数 c 上，总的扩散系数变成 c+c_art，对应的Pe数变成

为了确保Pe不会大于1，调整参数（tuning parameter）最大为0.5。其越小，计算就越稳定，所以COMSOL建议取0.25。因为这种方法在所有方向上添加了扩散，即各项同性，记为isotropic diffusion。它也是不一致的(inconsistent)，因为其添加的扩散量跟数值解与解析解的接近程度无关。

需要注意的是，这种人工扩散量的大小依赖于网格尺寸。一方面，多的网格意味着少的人工扩散量，但是会增加计算量。另一方面，少的网格需要更多的扩散量，会影响解的精度。

图中可以看出，Inconsistent相比consistent引入了更多的扩散，会有更大的误差，所以COMSOL中的默认选项为consistent，接下来介绍这种方法。

Consistent stabilization: streamline and crosswind diffusion

相比Inconsistent，consistent方法在数值解和解析解接近的地方添加更少的扩散，所以具有更小的误差。即是说，网格足够密的区域，无需引入人工扩散。

streamline方法即是，仅在流线方向上添加人工扩散，也叫上风格式（upwinding？）。从下图可以看出（蓝线），采用这种方法时，平滑区域震荡被消除，而大梯度区域震荡仍存在，为了解决这一问题，引入crosswind方法。

crosswind不仅添加扩散在cross方向，并且可以捕获非连续性区域，移除数值波动。

因此通常采用consistent方法。

除了不同稳定方法，网格精度和误差精度也会对结果造成影响。

下图给出不同网格精度下的解。

下图给出时间相关求解器下不同误差精度的解。

技巧tips

稀释和浓缩粒子的输运

Transport of Diluted Species模块提供了一个额外的选项，关于残差（residual）的计算。如图，有两个选择，approximation和full。通常full没有必要，而approximation忽略了扩散张量的导数。
另外，crosswind有两个选项Do Carmo and Galeao，Codina。都适用于各向异性网格，后者比前者添加扩散量更小，结果就可能更震荡。

流体

流体模块的求解通常不需要对稳定方法进行修改。但是，当使用默认的独立于Pe数的P1+P1元素时，NS方程计算不稳定，需要添加稳定方法。Pm+Pn，速度被求解使用m阶计算元素，压力被求解使用n阶元素。
默认的稳定方法会被移除，当速度阶数高于压力阶数（P2+P1，P3+P2）或Pe数很小时。

自定义方程的求解

如果要求解自定义的微分方程，需要做一些文献调研，了解模型需要采用的稳定方法。首次的计算可以尝试添加一个人工扩散项通过自定义的方法，具体可以参看官方案例（http://cn.comsol.com/model/flow-of-oldroyd-b-viscoelastic-fluid-4383）。

总结

没有必要修改稳定方法的默认设置！

不同于那些基于网格的数值方法（有限元，有限体积，有限差分），Lagrangian particle-based 方法对于较高Pe数问题的模拟比较有效，无需人工扩散项。粒子追踪模块就是采用这种方法（http://cn.comsol.com/particle-tracing-module）。对于基于网格方法和基于粒子方法都有相应的计算案例可以参考（http://cn.comsol.com/model/laminar-static-mixer-245；http://cn.comsol.com/model/particle-trajectories-in-a-laminar-static-mixer-10644），不过更多内容要参考COMSOL Multiphysics Reference Manual。