Code of the Week #2: ddt(U)

本周,我们讨论一下这一个语句:

fvm::ddt(U);

这一项是动量方程

\frac{\partial U}{\partial t} + \nabla\cdot(UU) = -\nabla p + \nabla\cdot\tau

的左边第一项。

求解这一方程时,需要对 pU 进行解耦。但是,本周,我们不讨论解耦的问题,而讨论一下时间项是如何离散的。为方便说明,设

\frac{\partial U}{\partial t} = S

右边的 S 为源项。

大家都知道,FVM 方法是建立在离散的空间网格上的,对这一项进行离散(积分)的方法就是

\iiint_V\int_t^{t+\Delta t} \frac{dU}{dt} = [\varphi(t+\Delta t) -\varphi(t)] \Delta V = S\Delta V \Delta t

变换一下,

[\varphi(t+\Delta t) -\varphi(t)] \frac{\Delta V}{\Delta t} = S\Delta V

则在新的时间步

\varphi(t+\Delta t)\frac{\Delta V}{\Delta t} = S\Delta V + \varphi(t)\frac{\Delta V}{\Delta t}

S=0 ,则有

\varphi(t+\Delta t)\frac{\Delta V}{\Delta t} = \varphi(t)\frac{\Delta V}{\Delta t}

下面,以 \varphi^{n+1} 代替 \varphi(t+\Delta t) , 以 \varphi^{n} 代替 \varphi(t)
如果以一均分网格作为示例,


\left[ \begin{array}% \frac{V}{\Delta t} & 0 & 0 & 0\\ 0 & \frac{V}{\Delta t} & 0 & 0 \\ 0 & 0 & \frac{V}{\Delta t} & 0\\ 0 & 0 & 0 & \frac{V}{\Delta t} \end{array} \right] \times \left[ \begin{array}% \varphi_1^{n+1}\\ \varphi_2^{n+1}\\ \varphi_3^{n+1}\\ \varphi_4^{n+1}\\ \end{array} \right] = \left[ \begin{array}% \varphi_1^n\frac{V}{\Delta t}\\ \varphi_2^n\frac{V}{\Delta t}\\ \varphi_3^n\frac{V}{\Delta t}\\ \varphi_4^n\frac{V}{\Delta t}\\ \end{array} \right]

获得这一矩阵后,再与动量方程的其他项离散后得到的矩阵相加,如上一期的 Code of the Week #1: div(phi,U) ,就可以形成动量方程离散后的完整矩阵方程。

以下程序段为 foam-extend-3.1 内, 源文件 src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C
从第 260 行起:

template<class Type>                          // 类模板
tmp<fvMatrix<Type> >                          // 返回类型: fvMatrix 矩阵
EulerDdtScheme<Type>::fvmDdt                  // 类名::函数名,定义
(                                             // 对 vf 进行时间项离散
    GeometricField<Type, fvPatchField, volMesh>& vf  //  vf 是参数
)
{
    tmp<fvMatrix<Type> > tfvm                        // 矩阵定义及初始化,返回的变量
    (
        new fvMatrix<Type>
        (
            vf,
            vf.dimensions()*dimVol/dimTime
        )
    );

    fvMatrix<Type>& fvm = tfvm();

    scalar rDeltaT = 1.0/mesh().time().deltaT().value(); // 1/deltaT

    fvm.diag() = rDeltaT*mesh().V();                     // 对角线 V/deltaT

    if (mesh().moving())          // 网格是否有移动
    {
        fvm.source() = rDeltaT*vf.oldTime().internalField()*mesh().V0();
    }
    else
    {
        fvm.source() = rDeltaT*vf.oldTime().internalField()*mesh().V(); // RHS,vf*V/deltaT
    }

    return tfvm;
}

在程序中,返回的矩阵里面 0 元素位置上是空的,以节约内存。也就是说,返回的值中,包括矩阵对角线上的值,和等式右边的值 (源项)。 S 在 ddt() 中没有体现,将另外说明。

6 Likes

没太看明白,矩阵方程等候右边的phi上标表示什么意思?如果是表达时间推进,这样做是不是更好点:
%E5%9B%BE%E7%89%87

2 Likes

如果使用 \varphi_i^{n+1} = S\Delta t + \varphi_i^n 这一类方式,通过前一时刻的值来直接计算后一时刻的值,就是时间项的显式离散。如果通过解矩阵方程求解的话,则是隐式离散。

1 Like

OpenFoam 提供三种方式:

  1. Euler
\frac{\partial \varphi}{\partial t} = \frac{\varphi^{n+1} - \varphi^n}{\Delta t}
  1. Crank-Nicolson

  2. Backward