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=∑∣U⃗f⋅S⃗f∣(1)sumPhi=\sum|\vec{U}_f\cdot\vec{ S}_f| \tag{1} sumPhi=∑∣Uf⋅Sf∣(1)
CoNum=12max(∣U⃗f⋅S⃗f∣VΔT)(2)CoNum=\frac 12 max(\frac{|\vec{U}_f\cdot\vec{ S}_f| }{V}\Delta T )\tag 2CoNum=21max(V∣Uf⋅Sf∣ΔT)(2)
meanCoNum=12sum(∣U⃗f⋅S⃗f∣VΔT)(3)meanCoNum=\frac 12 sum(\frac{|\vec{U}_f\cdot\vec{ S}_f| }{V}\Delta T )\tag 3meanCoNum=21sum(V∣Uf⋅Sf∣ΔT)(3)
注：在 OpenFOAM 中，若为六面体网格，则其对应的面通量值ϕ\phiϕ加和加了两次，因此公式 (2) 和(3) 需要除 2。

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

上述代码用数学公式可表述为：
sumPhi=(ϕ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)∑∣Uf⋅Sf∣,0.01<ϕ1<0.99(5)
则两相库朗数计算结果为：
sumPhi≈ϕ1phi2∑∣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∑∣Uf⋅Sf∣,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=∑∣ρfU⃗f⋅S⃗f∣2∗Comax⋅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⋅ρ∑∣ρfUf⋅Sf∣(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∗Coαmax⋅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∑∣Uf⋅Sf∣,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-局部时间步，自动调节步长技术相关推荐

非稳态计算时间步和最大迭代数的设定（分享）
fluent中max iteration pertime step.time step size.number of time steps的设定来源于:https://zhuanlan.zhihu. ...
降低百倍时间步，精度媲美传统神经网络：上交等机构提出ANN-SNN转换框架
©作者 | 刘方鑫,赵文博,蒋力来源 | 机器之心脉冲神经网络(Spiking Neural Network, SNN)被誉为第三代的神经网络,以其丰富的时空领域的神经动力学特性.多样的编码机制. ...
我业余时间如何学习多门技术以及开发业余项目的一些心得
我的学习方法在开始介绍我开发业余项目的经验前,先和大家分享一下我的学习方法吧,不过我认为学习这件事情因人而异,每个人都有适合自己的学习方式,所以这是作为一个参考,希望大家也都能先找到适合自己的学习方 ...
[MATLAB粒子模拟笔记]初始化半个时间步的位置
function particle = position(particle,prm) % Update the position in one stepslx = prm.slx p = partic ...
NLP：自然语言处理技术领域的代表性算法概述(技术迭代路线图/发展时间路线)、四大技术范式变迁概述(统计时代→大模型时代)、四个时代的技术方法论探究(少数公司可承担的训练成本原因)之详细攻略
NLP:自然语言处理技术领域的代表性算法概述(技术迭代路线图/发展时间路线).四大技术范式变迁概述(统计时代→大模型时代).四个时代的技术方法论探究(少数公司可承担的训练成本原因)之详细攻略目录一 ...
源码学习之LAMMPS的一个时间步是如何工作的
从本推文开始将开始学习LAMMPS的源码.这部分我也不太了解,所以一边学一边分享.学习源码的最好方法当然后阅读手册的PROGRAMMER GUIDE.为了督促自己学习,我首先将PROGRAMMER G ...
matlab找到非定常涡流的每个时间步的涡的涡核位置和这个涡环量(以及重叠网格扑翼流场的涡动力学参数求解的解决方案)
如果提前已经算好fluent每个时间步的cas和dat文件,那么可以用jou文件导出每个节点的u,v,dudy,dvdx,overset-cell-type.fluent中只能直接导出涡量幅值,想要计 ...
【Fluent】导出瞬态计算过程每一秒或每一个时间步的各个坐标/节点的物理量-温度场-压力场,TUI命令/file/export和/file/transient-export
一.功能需求如果你进行的是稳态计算,你需要将物理场中的每一个节点上的物理量数据(例如温度.压力)导出成类似txt或Excel表格的文件. 文件里的内容形式是:每一行中有节点ID.节点的XYZ坐标.物 ...
怀旧服服务器维护后断牙刷新吗,魔兽世界怀旧服：猎人抓断牙方法介绍，时间和位面技术缺一不可...
原标题:魔兽世界怀旧服:猎人抓断牙方法介绍,时间和位面技术缺一不可如果问魔兽世界怀旧服猎人最想获得的宠物是什么,相信很多猎人玩家都会脱口而出"断牙".断牙作为猫科宠物,虽然不是唯 ...

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

LTS-局部时间步，自动调节步长技术相关推荐

最新文章

热门文章