LTS-USB迷|专注于互联网分享

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;

s u m P h i = ∑ ∣ U ⃗ f ⋅ S ⃗ f ∣ (1) sumPhi=\sum|\vec{U}_f\cdot\vec{ S}_f| \tag{1} sumPhi=∑∣U f⋅S f∣(1)
C o N u m = 1 2 m a x ( ∣ U ⃗ f ⋅ S ⃗ f ∣ V Δ T ) (2) CoNum=\frac 12 max(\frac{|\vec{U}_f\cdot\vec{ S}_f| }{V}\Delta T )\tag 2 CoNum=21max(V∣U f⋅S f∣ΔT)(2)
m e a n C o N u m = 1 2 s u m ( ∣ U ⃗ f ⋅ S ⃗ f ∣ V Δ T ) (3) meanCoNum=\frac 12 sum(\frac{|\vec{U}_f\cdot\vec{ S}_f| }{V}\Delta T )\tag 3 meanCoNum=21sum(V∣U f⋅S f∣ΔT)(3)
注：在 OpenFOAM 中，若为六面体网格，则其对应的面通量值 ϕ \phi ϕ加和加了两次，因此公式 (2) 和(3) 需要除 2。

库朗数一般定义如下：
C o = ∣ U ⃗ f ⋅ Δ T ∣ Δ x (4) Co=\frac {|\vec U_f \cdot \Delta T| } {\Delta x}\tag 4 Co=Δx∣U f⋅ΔT∣(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;

上述代码用数学公式可表述为：
s u m P h i = ( ϕ 1 − 0.01 ) ( 0.99 − ϕ 1 ) ∑ ∣ U ⃗ f ⋅ S ⃗ f ∣ , 0.01 < ϕ 1 < 0.99 (5) sumPhi=(\phi_1-0.01)(0.99-\phi_1)\sum|\vec{U}_f\cdot\vec{ S}_f|,0.01<\phi_1<0.99 \tag5 sumPhi=(ϕ1−0.01)(0.99−ϕ1)∑∣U f⋅S f∣,0.01<ϕ1<0.99(5)
则两相库朗数计算结果为：
s u m P h i ≈ ϕ 1 p h i 2 ∑ ∣ U ⃗ f ⋅ S ⃗ f ∣ , 0.01 < ϕ 1 < 0.99 (6) sumPhi\approx \phi_1phi_2\sum|\vec{U}_f\cdot\vec{ S}_f|,0.01<\phi_1<0.99 \tag6 sumPhi≈ϕ1phi2∑∣U f⋅S f∣,0.01<ϕ1<0.99(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()));

若第二项的值相对较大，则其数学公式表述如下：
Δ T = ∑ ∣ ρ f U ⃗ f ⋅ S ⃗ f ∣ 2 ∗ C o m a x ⋅ V ⋅ ρ (7) \Delta T =\frac {\sum|\rho_f \vec{U}_f \cdot \vec{S}_f|} {2*Co_{max} \cdot V \cdot \rho} \tag 7 ΔT=2∗Comax⋅V⋅ρ∑∣ρfU f⋅S f∣(7)

 if (maxAlphaCo < maxCo){// Further limit the reciprocal time-step// in the vicinity of the interfacevolScalarField alpha1Bar(fvc::average(alpha1));rDeltaT.ref() = max(rDeltaT(),pos(alpha1Bar() - alphaSpreadMin)*pos(alphaSpreadMax - alpha1Bar())*fvc::surfaceSum(mag(phi))()()/((2*maxAlphaCo)*mesh.V()));}

这里，我们有
Δ T = ( ϕ 1 − 0.01 ) ( 0.99 − ϕ 1 ) ∑ ∣ U ⃗ f ⋅ S ⃗ f ∣ 2 ∗ C o α m a x ⋅ V , 0.01 < ϕ 1 < 0.99 (8) \Delta T =(\phi_1-0.01)(0.99-\phi_1) \frac {\sum|\vec{U}_f\cdot\vec{ S}_f|}{2*Co_{\alpha _{max}}\cdot V},0.01<\phi_1<0.99 \tag8 ΔT=(ϕ1−0.01)(0.99−ϕ1)2∗Coαmax⋅V∑∣U f⋅S f∣,0.01<ϕ1<0.99(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主要解决稳态问题

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

LTS