Code of the Week #2: ddt(U)

本周，我们讨论一下这一个语句：

fvm::ddt(U);

这一项是动量方程 $$ \frac{\partial U}{\partial t} + \nabla\cdot(UU) = -\nabla p + \nabla\cdot\tau $$ 的左边第一项。

求解这一方程时，需要对 $p$ 和 $U$ 进行解耦。但是，本周，我们不讨论解耦的问题，而讨论一下时间项是如何离散的。为方便说明，设 $$ \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)] \Delta V/\Delta T = S\Delta V $$ 则在新的时间步 $$ \varphi(t+\Delta t)\frac{\Delta V}{\Delta t} = S/\Delta t + \varphi(t)\frac{\Delta V}{\Delta t} $$

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

template<class Type>                                        // 类模板
tmp<fvMatrix<Type> >                                        // 返回类型: fvMatrix 矩阵
EulerDdtScheme<Type>::fvmDdt                      // 类名::函数名，定义
(
    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();    // 把旧值作为源项，vf*V/deltaT
    }

    return tfvm;
}

更详细的讨论及版本请至：Code of the Week #2: ddt(U)