考虑如下所示的关于 $\psi$ 的输运方程:

\[ \begin{equation} \frac{\partial \psi}{\partial t}+\nabla \cdot (\mathbf u \psi)-\nabla \cdot \gamma \nabla \psi=S_{\psi} \label{eq:trans} \end{equation} \]

其中,$S_\psi$ 为源项。

源项线性化

实际情况下 $S_\psi$ 是比较复杂的,可能包含非线性项。为了求解稳定,通常对其进行线性化处理

\[ \begin{equation} S_\psi = S_u + S_p \psi \end{equation} \]

其中,$S_u$ 和 $S_p$ 可与 $\psi$ 有关,也可与 $\psi$ 无关。

将输运方程离散成 $[\mathbf A] [\mathbf x] = [\mathbf b]$ 后, $S_u$ 将进入 $[\mathbf b]$,而 $S_p$ 将变成 $-S_p$ 进入 $[\mathbf A]$ 的对角元素中。若 $S_p$ 为正,则将削弱 $[\mathbf A]$ 的对角占优,很可能造成线性方程组求解发散。因此在对源项进行线性化时必须保证 $S_p$ 为负或零。

源项线性化的选择并不是唯一的,不同的选择有不同的精度和稳定性。举个例子说明,考虑如下形式的源项

\[ \begin{equation} S_\psi = 10 - 2 \psi^3 \end{equation} \]

上式中包含非线性项 $-2\psi^3$,需要对其进行线性化处理。这里有好几种方法可以选择:

  • 方法一:$S_u=10 - 2( \psi^o )^3$ , $S_p = 0$
  • 方法二:$S_u=10$, $S_p = -2( \psi^o )^2$
  • 方法三(Picard's method):对 $S_\psi(\psi)$ 在 $\psi^o$ 处进行泰勒展开,只保留一阶项

\[ \begin{equation} \begin{aligned} S_\psi &\approx S_\psi(\psi^o) + \frac{\partial S_\psi}{\partial \psi}\Bigr|_{\psi=\psi^o}(\psi-\psi^o) \\ &= 10-2(\psi^o)^3 + [-6(\psi^o)^2](\psi-\psi^o) \\ &= 10+4(\psi^o)^3 - 6(\psi^o)^2\psi \end{aligned} \end{equation} \]

即 $S_u = 10+4( \psi^o )^3$, $S_p = -6( \psi^o )^2$。

下图给出了以上三种不同源项线性化方法的示意图,可以看出,第一种方法精度最低,第二种其次,第三种方法的精度最高。

不同源项线性化方法示例

不同源项线性化方法示例

源项线性化的用途

指定网格单元的值

对某一单元 $P$,若需要指定该处的值为 $\psi_{P, desired}$,则可对该单元施加如下线性化后的源项

\[ \begin{equation} S_P = A_{large}\psi_{P, desired} - A_{large} \psi_P \end{equation} \]

其中 $A_{large}$ 为一个很大的数,如 $10^{10}$。

对于单元 $P$,将除源项外的所有其他项离散后的方程表示为

\[ \begin{equation} A_P \psi_P + \sum_{N}^{nb}A_N\psi_N = b_P \end{equation} \]

加上源项后,上式变为

\[ \begin{equation} (A_P+A_{large}) \psi_P + \sum_{N}^{nb}A_N\psi_N = b_P + A_{large}\psi_{P, desired} \end{equation} \]

若满足 $A_{large} \gg A_P$ 、 $A_{large} \gg A_N$ 且 $A_{large}\psi_{P, desired} \gg b_P$ ,则求解得到的 $\psi_P \approx \psi_{P, desired}$ 。

在需要修改某个网格单元的值时而不影响其他单元时(如重叠网格的插值单元赋值),这种方法格外有用。

保证物理量为非负

当源项为负数的时候,求解方程得到的值也可能为负数,而这对一些不可能为负数的物理量(例如温度 $T$、湍动能 $k$、特定耗散率 $\omega$)来说是没有意义的。

为了避免这种情况发生,可以将源项做如下处理:

\[ \begin{equation} S_\psi = S_{const} \approx S_{const} \frac{\psi}{\psi^o} \end{equation} \]

其中,$S_{const}$ 为负的源项。经过以上处理后,$-S_{const}/\psi^o$ 将作为正数计入到系数矩阵的对角元素中,既增强了对角占优,又可保证 $\psi$ 为非负。

OpenFOAM 的实现

方程 $\eqref{eq:trans}$ 用代码可以表示为:

1
2
3
4
5
6
7
8
fvScalarMatrix psiEqn
(
    fvm::ddt(psi)
  + fvm::div(phi, psi)
  - fvm::laplacian(gamma, psi)
 ==
    S_psi   
);

这里,S_psifvScalarMatrix 对象,一般通过 fvm::Su/fvm::Sp/fvm::SuSp、fvOptions 等得到12

== 运算优先级低于 +-,因此 == 最后计算。上面这段代码执行时,先计算 == 左右的表达式,得到两个 fvScalarMatrix 对象,再将这两个对象进行 == 运算。== 操作符已经被重载,实际为 - 操作符:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
// src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
template<class Type>
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator==
(
    const fvMatrix<Type>& A,
    const fvMatrix<Type>& B
)
{
    checkMethod(A, B, "==");
    return (A - B);
}

最后做的是类似 fvmA - fvmB 的操作,源项放在 fvmB 中。

OpenFOAM 对源项的实现比较灵活,有多种实现方式,这里介绍比较常见的 SuSp。

Su 和 Sp

对于 fvmB 而言,$S_u$ 将从 $[\mathbf b]$ 中被减去,相应代码如下:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
// src/finiteVolume/finiteVolume/fvm/fvmSup.C
template<class Type>
Foam::tmp<Foam::fvMatrix<Type>>
Foam::fvm::Su
(
    const DimensionedField<Type, volMesh>& su,
    const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    const fvMesh& mesh = vf.mesh();

    tmp<fvMatrix<Type>> tfvm
    (
        new fvMatrix<Type>
        (
            vf,
            dimVol*su.dimensions()
        )
    );
    fvMatrix<Type>& fvm = tfvm.ref();

    fvm.source() -= mesh.V()*su.field();

    return tfvm;
}

同样,对于 fvmB,$S_p$ 将加到 $[\mathbf A]$ 的对角上,相应代码如下:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
// src/finiteVolume/finiteVolume/fvm/fvmSup.C
template<class Type>
Foam::tmp<Foam::fvMatrix<Type>>
Foam::fvm::Sp
(
    const volScalarField::Internal& sp,
    const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    const fvMesh& mesh = vf.mesh();

    tmp<fvMatrix<Type>> tfvm
    (
        new fvMatrix<Type>
        (
            vf,
            dimVol*sp.dimensions()*vf.dimensions()
        )
    );
    fvMatrix<Type>& fvm = tfvm.ref();

    fvm.diag() += mesh.V()*sp.field();

    return tfvm;
}

源项自动处理

对于复杂的源项,在不同位置的正负可能不同。为了使整个线性方程组的求解更稳定,需要对正负源项分开处理。具体做法如下:

  • 若 $S_\psi<0$,则对其采用隐式离散,即 $S_u=0$,$S_p=S_\psi/\psi^o$,将其贡献计入 $[\mathbf A]$ 对角。

  • 若 $S_\psi>0$,则对其采用显式离散,即 $S_u=S_\psi$,$S_p=0$,将其贡献计入 $[\mathbf b]$ 中。

以上过程通过 fvm::SuSp 这个函数实现:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
// src/finiteVolume/finiteVolume/fvm/fvmSup.C
template<class Type>
Foam::tmp<Foam::fvMatrix<Type>>
Foam::fvm::SuSp
(
    const volScalarField::Internal& susp,
    const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    const fvMesh& mesh = vf.mesh();

    tmp<fvMatrix<Type>> tfvm
    (
        new fvMatrix<Type>
        (
            vf,
            dimVol*susp.dimensions()*vf.dimensions()
        )
    );
    fvMatrix<Type>& fvm = tfvm.ref();

    fvm.diag() += mesh.V()*max(susp.field(), scalar(0));

    fvm.source() -= mesh.V()*min(susp.field(), scalar(0))
        *vf.primitiveField();

    return tfvm;
}

这个函数可以将自动处理正负源项,对 fvmB 而言:

  • susp 大于 0 ,则执行 fvm.diag() += mesh.V()*susp.field(),这对应 $S_\psi < 0$。
  • susp 小于 0 ,则执行 fvm.source() -= mesh.V()*susp.field()*vf.primitiveField(),这对应 $S_\psi > 0$。

注意事项

fvm::Sp 在使用时,必须保证得到的对角元素为负,这样才能增强对角占优。若传递的sp为正,则需要在前面加上负号,形式大致如下:

1
fvmA == - fvm::Sp(sp, psi)

同理,fvm::SuSp 写在 == 右边时,按照约定,必须在前面加负号,形式大致如下:

1
fvmA == - fvm::SuSp(susp, psi)

应用示例

以 Spalart-Allmaras 这个一方程湍流模型为例。原始文献3 中需要求解的输运方程形式如下:

\[ \begin{equation} \frac{\partial (\rho \tilde{\nu})}{\partial t} + \underbrace{\nabla \cdot (\rho {\tilde{\nu} \mathbf u})}_{\mathrm{advection}} = \underbrace{\left[ \nabla \cdot \left(\rho D_{\tilde{\nu}} \nabla \tilde{\nu} \right) + \frac{1}{\sigma} C_{b2} \rho (\nabla \tilde{\nu})^2 \right]}_{\mathrm{diffusion}} + \underbrace{C_{b1} (1 - f_{t2}) \rho \tilde{S} \tilde{\nu}}_{\mathrm{procduction}} - \underbrace{\left(C_{w1} f_w - \frac{C_{b1}}{\kappa^2} f_{t2}\right) \rho \left( \frac{\tilde{\nu}}{d} \right)^2}_{\mathrm{destruction}} \end{equation} \]

OpenFOAM 的代码没有实现 $f_{t2}$ 项,这在 SpalartAllaras.H 文件中已经说明:

The model is implemented without the trip-term and hence the ft2 term is not needed.

于是公式简化成

\[ \begin{equation} \frac{\partial (\rho \tilde{\nu})}{\partial t} + \underbrace{\nabla \cdot (\rho {\tilde{\nu} \mathbf u})}_{\mathrm{advection}} = \underbrace{\left[ \nabla \cdot \left(\rho D_{\tilde{\nu}} \nabla \tilde{\nu} \right) + \frac{1}{\sigma} C_{b2} \rho (\nabla \tilde{\nu})^2 \right]}_{\mathrm{diffusion}} + \underbrace{C_{b1} \rho \tilde{S} \tilde{\nu}}_{\mathrm{procduction}} - \underbrace{\left(C_{w1} f_w \right) \rho \left( \frac{\tilde{\nu}}{d} \right)^2}_{\mathrm{destruction}} \end{equation} \]

代码实现如下:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
// src/TurbulenceModels/turbulenceModels/RAS/SpalartAllmaras/SpalartAllmaras.C
    tmp<fvScalarMatrix> nuTildaEqn
    (
        fvm::ddt(alpha, rho, nuTilda_)
      + fvm::div(alphaRhoPhi, nuTilda_)
      - fvm::laplacian(alpha*rho*DnuTildaEff(), nuTilda_)
      - Cb2_/sigmaNut_*alpha*rho*magSqr(fvc::grad(nuTilda_))
     ==
        Cb1_*alpha*rho*Stilda*nuTilda_
      - fvm::Sp(Cw1_*alpha*rho*fw(Stilda)*nuTilda_/sqr(y_), nuTilda_)
      + fvOptions(alpha, rho, nuTilda_)
    );

我们重点关注 production 和 destruction 项。其中,production 项大于0,用显式处理:

1
        Cb1_*alpha*rho*Stilda*nuTilda_

这里也可以使用 fvm::Su

1
        fvm::Su(Cb1_*alpha*rho*Stilda*nuTilda_, nuTilda_)

destruction 项小于0,用隐式处理:

1
      - fvm::Sp(Cw1_*alpha*rho*fw(Stilda)*nuTilda_/sqr(y_), nuTilda_)

参考资料