OpenFOAM中matrix relax的bug

李东岳

背景：

为了方便google搜索，一些用英文来表示。

OpenFOAM中的松弛分为matrix relax，以及field relax。在field relax中，理论上新的时间步的值与旧的时间步之间的差不会那么大，就是说$T^{t+\dt}=T^{t}+\alpha(T^{t+\dt}-T^t)$。在OpenFOAM的代码中，也是这样植入的。

现在讨论matrix relax，在matrix relax中，理论上，

给定松弛因子为1（$\alpha=1$），不会影响计算结果，不会改变矩阵结构。
~~给定任意的松弛因子，不会影响计算结果~~。这句话不对。更新：给定任意的松弛因子，同样影响计算结果。这样的话，OpenFOAM跟理论结果就一样了。只不过就是松弛的计算方法不一样

目前在OpenFOAM代码中，如果对matrix relax，存在以下问题：

给定松弛因子为1（$\alpha=1$），会改变矩阵结构，基本不影响计算结果
给定非常小的松弛因子（$\alpha=0.01$），改变矩阵结构，影响瞬态计算结果，无时间推进特征。

测试方法：

在瞬态pisoFoam求解器中，在fvSolution中对速度方程增加matrix relax，

relaxationFactors
{
    equations
    {
        U 0.01;
    }
}

查看结果。会发现所有的场都没有按照时间步进行推进。

原因：

fvMatrix.C中的relax方法并没有按照传统的松弛方法来进行。OpenFOAM中的松弛方法非常简单。考虑下面的离散矩阵，假定非对角线系数比对角线系数的绝对值大：
$$
-a_NU_N^{t+\Delta t}+\left(\frac{1}{\Delta t}+a_P\right )U_P^{t+\Delta t}=\frac{1}{\Delta t}U_P^t
$$
如果矩阵松弛因子给1。OpenFOAM中的松弛方法为：
$$
-a_NU_N^{t+\Delta t}+a_NU_P^{t+\Delta t}=\frac{1}{\Delta t}U_P^t+\left(a_N-\left(\frac{1}{\Delta t}+a_P\right ) \right )U_P^{t}
$$
理论上就是将对角线系数增加，然后源项（采用已知量）也进行相应增加来互相抵消。这种方法如果不考虑边界条件，最差最差会创立一个对角相等矩阵。一般会创立一个对角占优矩阵。下面是Henry的解释

Correct, relaxation is applied to a matrix which is first made at least diagonally equal, otherwise it would not improve diagonal dominance. Thus relaxation with a factor of 1 makes the matrix at least diagonally equal but no more.

如果矩阵松弛因子给非常小的值，比如0.01，变为：

$$
-a_NU_N^{t+\Delta t}+\frac{a_N}{0.01}U_P^{t+\Delta t}=\frac{1}{\Delta t}U_P^t+\left(\frac{a_N}{0.01}-\left(\frac{1}{\Delta t}+a_P\right ) \right )U_P^{t}
$$

这会导致方程近似为

$$
\frac{a_N}{0.01}U_P^{t+\Delta t}=\frac{a_N}{0.01}U_P^{t}
$$

这会创立一个非常强的对角相等矩阵。但下一个时间步的值，基本等于上一个时间步的值，因此不会实现真实的时间步推进。

关联问题：

pisoFoam中的UEqn.relax();也进行方程松弛？

关于亚松弛的疑问

关于压力松弛因子的疑问

李东岳

另外要注意的是。瞬态求解器没有任何tutorials对matrix做了relax，因此对计算结果不会产生影响。但是自己调试的算例一定要注意一下。这个OpenFOAM代码跟我们relax的理论有区别。

xpqiu

@李东岳在 OpenFOAM中matrix relax的bug 中说：

给定任意的松弛因子，不会影响计算结果。

这个要成立，前提条件是需要引入额外的迭代。不管是 matrix relax 还是 field relax，本质上都是把上个步骤的值跟当前步的值做一个混合之后来作为当前步的值，也就是上面提到的
$$
T^{n+1}=T^{n}+\alpha(T^{n+1}-T^n)
$$
matrix relax 虽然具体实现方式不同，但是本质不变。为避免混淆，这里使用 n+1 表示当前步，n表示上一步，n 更多的时候是内迭代步，而不是时间。
那么要想实现任意松弛因子对最终结果没有影响，唯一可能的就是迭代到一个稳定的值，也就是
$$
T^{n+1} = T^{n}
$$

如果对 PISO 这样的非迭代式算法（指的是一个时间步内动量方程只求解一次）使用 matrix relax，而且设置很小的值，那么肯定会影响计算结果的啊，比如你设置U 的松弛为0.01，那么就相当于说每次都是使用99%的上一步的值加上1%的当前步算出来的值，来作为当前步的值，这样肯定就会表现为时间上基本无推进了。

OpenFOAM 里面的 PIMPLE 算法的 matrix relax 是可以使用小于1的松弛因子的，因为 PIMPLE 算法可以理解为在时间推进步里面还有一层内迭代步，也是PIMPLE可以算比较大的 Co 数的主要原因。姑且认为从这一个时间步推进到下一个时间步都存在准确解，那么 PISO 相当于是通过一次求解动量方程，就把下一步的解算出来，PIMPLE 使用 relax，则相当于是把从上一个时间步到下一步时间步推进这个过程，使用一种类似稳态方法的方式来进行，使用小于1的松弛因子，同时增加迭代步数，这样计算会更稳定一些尤其是大时间步情况下。通过这样的迭代（多次求解动量方程和压力方程），可以逼近达到
$$
T^{n+1} = T^{n}
$$
这里的n就代表内迭代步了。从而实现虽然使用了松弛因子，但是松弛因子对结果基本无影响这个效果。

另外，回到 OpenFOAM，或者说跟 OpenFOAM 类似的使用 collocated 网格的程序，都需要使用 Rhie-Chow 插值来解决棋盘压力问题。但是原始的 Rhie-Chow 方法都做不到仿真结果不受松弛因子影响。OpenFOAM 里面，可以尝试用最简单的二维 cavity 算例，用 simpleFoam 求解，一直算下去直到 initial residual 小于 1e-15，换不同的松弛因子来算，会发现虽然都能收敛，但是最终会收敛到不同的结果，虽然可能差异不会很大。

现在也有一些对 Rhie-Chow 的修正，比如

Cubero, A., & Fueyo, N. (2007). A compact momentum interpolation procedure for unsteady flows and relaxation. Numerical Heat Transfer, Part B: Fundamentals, 52(6), 507–529. https://doi.org/10.1080/10407790701563334

和

Cubero, A., Sánchez-Insa, A., & Fueyo, N. (2014). A consistent momentum interpolation method for steady and unsteady multiphase flows. Computers and Chemical Engineering, 62, 96–107. https://doi.org/10.1016/j.compchemeng.2013.12.002

李东岳

我思考了下，

matrix relax 虽然具体实现方式不同，但是本质不变。

你说的是对的。一楼我说的应该有问题（我更正一下）。不仅仅field relax，matrix relax，也会影响计算结果。我当时把matrix relax理解成不会影响计算结果了。

matrix relax是使得计算结果更慢的趋向于真实值，但换来一个对角占优矩阵
field relax是使得计算结果更慢的趋向于真实值

两个都使得计算结果更慢的趋向于真实值。

当时看这个公式的时候：$\phi^{n+1}=\phi^n+\beta(\phi^{n+1}-\phi^n)$，以为$\beta$变小时候，$\phi^{n+1}$还是原来那个值。然而并不是好尴尬