LTS-局部时间步，自动调节步长技术

在查看interFoam求解器时，会发现其中引入了LTS：

if (LTS)
        {
            #include "setRDeltaT.H"
        }

那么它的含义和功能是什么呢？

• LTS（locall time-step）是一种局部时间步求解器。该求解器建立于局部时间步下，它会提取pimple算法控制字典t里面的maxCo与maxAlphaCo，并以此为依据，根据本地计算条件，计算当前网格的时间步长deltaT。
• 由于每个网格有不同的 deltaT，这会使质量失衡（瞬态问题每个时间步的质量守恒），因此需要迭代求解。在最终稳态的情况下，时间项的导数为 0 ，这时输出的结果和最终瞬态达到稳定状况下的结果相同。

LTS与库朗数是息息相关的，为此，我们先分析interFoam 中的 CourantNo.H ，了解其中库朗数的定义：

scalar CoNum = 0.0;    //定义标量 CoNum 并赋初值
scalar meanCoNum = 0.0;//定义标量 meanCoNum 并赋初值

if (mesh.nInternalFaces())
{
    scalarField sumPhi
    (
        fvc::surfaceSum(mag(phi))().primitiveField()
    );

    CoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();

    meanCoNum =
        0.5*(gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaTValue();
}

Info<< "Courant Number mean: " << meanCoNum << " max: " << CoNum << endl;

$$sumPhi=\sum |{\vec{U}}_f\cdot {\vec{S}}_f|\text{\hspace{1em}(1)}$$
$$CoNum=\frac{1}{2}max(\frac{|{\vec{U}}_f\cdot {\vec{S}}_f|}{V}\Delta T)\text{\hspace{1em}(2)}$$
$$meanCoNum=\frac{1}{2}sum(\frac{|{\vec{U}}_f\cdot {\vec{S}}_f|}{V}\Delta T)\text{\hspace{1em}(3)}$$
注：在 OpenFOAM 中，若为六面体网格，则其对应的面通量值$\varphi$加和加了两次，因此公式 (2) 和(3) 需要除 2。

库朗数一般定义如下：
$$Co=\frac{|{\vec{U}}_f\cdot \Delta T|}{\Delta x}\text{\hspace{1em}(4)}$$
获得了库朗数后，进一步求取alpha 库郎数，参见 alphaCourantNo.H：

scalar maxAlphaCo
(
	//查找controlDict字典文件中的maxAlphaCo并进行读取
    readScalar(runTime.controlDict().lookup("maxAlphaCo")) 
);
//声明变量赋初值
scalar alphaCoNum = 0.0;
scalar meanAlphaCoNum = 0.0;
//求解公式
if (mesh.nInternalFaces())
{
    scalarField sumPhi
    (
        mixture.nearInterface()().primitiveField()
       *fvc::surfaceSum(mag(phi))().primitiveField()
    );

    alphaCoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();

    meanAlphaCoNum =
        0.5*(gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaTValue();
}

Info<< "Interface Courant Number mean: " << meanAlphaCoNum
    << " max: " << alphaCoNum << endl;

上述代码用数学公式可表述为：
$$sumPhi=({\varphi }_1-0.01)(0.99-{\varphi }_1)\sum |{\vec{U}}_f\cdot {\vec{S}}_f|,0.01<{\varphi }_1<0.99\text{\hspace{1em}(5)}$$
则两相库朗数计算结果为：
$$sumPhi\approx {\varphi }_{1}ph{i}_2\sum |{\vec{U}}_f\cdot {\vec{S}}_f|,0.01<{\varphi }_1<0.99\text{\hspace{1em}(6)}$$
同理可求得alphaCoNum和meanAlphaCoNum。

在了解了库朗数和两相库朗数后，我们进入setDeltaT.H文件

if (adjustTimeStep)//如果调节时间步
{
    scalar maxDeltaTFact =
        min(maxCo/(CoNum + SMALL), maxAlphaCo/(alphaCoNum + SMALL));
	//在没有达到设定的库朗数值的时候，deltaTFact每次扩大1.2倍
    scalar deltaTFact = min(min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2);
	//与controlDict字典中的 maxDeltaT比较，将DeltaT设置成较小的那个
    runTime.setDeltaT
    (
        min
        (
            deltaTFact*runTime.deltaTValue(),
            maxDeltaT
        )
    );

    Info<< "deltaT = " <<  runTime.deltaTValue() << endl;
}

下面是setRDeltaT.H,这里面声明了计算所需要的场量，以及如何计算局部时间步.由于该文件较长，这里只将其中的计算局部时间步的代码提取出来，

   rDeltaT.ref() = max
    (
        1/dimensionedScalar("maxDeltaT", dimTime, maxDeltaT),
        fvc::surfaceSum(mag(rhoPhi))()()
       /((2*maxCo)*mesh.V()*rho())
    );

若第二项的值相对较大，则其数学公式表述如下：
$$\Delta T=\frac{\sum |{\rho }_f{\vec{U}}_f\cdot {\vec{S}}_f|}{2\ast C{o}_{max}\cdot V\cdot \rho }\text{\hspace{1em}(7)}$$

 if (maxAlphaCo < maxCo)
    {
        // Further limit the reciprocal time-step
        // in the vicinity of the interface

        volScalarField alpha1Bar(fvc::average(alpha1));

        rDeltaT.ref() = max
        (
            rDeltaT(),
            pos(alpha1Bar() - alphaSpreadMin)
           *pos(alphaSpreadMax - alpha1Bar())
           *fvc::surfaceSum(mag(phi))()()
           /((2*maxAlphaCo)*mesh.V())
        );
    }

这里，我们有
$$\Delta T=({\varphi }_1-0.01)(0.99-{\varphi }_1)\frac{\sum |{\vec{U}}_f\cdot {\vec{S}}_f|}{2\ast C{o}_{{\alpha }_{max}}\cdot V},0.01<{\varphi }_1<0.99\text{\hspace{1em}(8)}$$

以下，通过fvc::smooth、fvc::spread、fvc::sweep继续光顺局部时间步。

alpha方程求解：

fvScalarMatrix alpha1Eqn
        (
            (
                LTS
              ? fv::localEulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
              : fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
            )
          + fv::gaussConvectionScheme<scalar>
            (
                mesh,
                phiCN,
                upwind<scalar>(mesh, phiCN)
            ).fvmDiv(phiCN, alpha1)
        );
alpha1Eqn.solve();

从这段代码我们可以看出，interFoam采用的是欧拉法进行时间离散，而LTS，则对时间项采用局部欧拉法进行离散。

注：

• LTS是一种可用可不用的自动调节步长方式，它是一种单纯的计算加速技术，能够加快收敛进程，但准确性还有待考究。
• LTS主要解决稳态问题

本文内容参考了东岳流体，关于求解器的相关知识可以查阅该网站。