关于icoFoam、pisoFoam等求解过程的连续性误差的疑惑
-
各位老师好,我关于piso算法进行速度修正有一点疑惑,想请教一下。
- U = HbyA - rAU*fvc::grad(p);修正后的速度场准确满足连续性条件吗?
- 在速度压力耦合计算结束后,如何构建满足连续性条件的phi?
谢谢!
以下是具体问题:
这是icoFoam中最后一部分构建压力泊松方程并且迭代求解的代码:while (piso.correctNonOrthogonal()) { // Pressure corrector fvScalarMatrix pEqn ( fvm::laplacian(rAU, p) == fvc::div(phiHbyA) ); pEqn.setReference(pRefCell, pRefValue); pEqn.solve(mesh.solver(p.select(piso.finalInnerIter()))); if (piso.finalNonOrthogonalIter()) { phi = phiHbyA - pEqn.flux(); } } #include "continuityErrs.H" U = HbyA - rAU*fvc::grad(p); U.correctBoundaryConditions();
下面是 #include "continuityErrs.H"
{ volScalarField contErr(fvc::div(phi)); scalar sumLocalContErr = runTime.deltaTValue()* mag(contErr)().weightedAverage(mesh.V()).value(); scalar globalContErr = runTime.deltaTValue()* contErr.weightedAverage(mesh.V()).value(); cumulativeContErr += globalContErr; Info<< "time step continuity errors : sum local = " << sumLocalContErr << ", global = " << globalContErr << ", cumulative = " << cumulativeContErr << endl; }
使用cavity进行测试发现此时的
time step continuity errors : sum local = 9.66354e-09,在1e-8这个量级,满足连续性条件。理论上最后根据(U = HbyA - rAUfvc::grad(p);)得到的流场应该满足连续性条件才对,那为什么是在U修正之前进行连续性的计算呢?如果在(U = HbyA - rAUfvc::grad(p);)之后再次构建phi,进行连续性判断会发现局部连续性误差显著增大。
... U = HbyA - rAU*fvc::grad(p); surfaceScalarField phi2 = fvc::flux(U); volScalarField contErr(fvc::div(phi2)); scalar sumLocalContErr = runTime.deltaTValue()* mag(contErr)().weightedAverage(mesh.V()).value(); scalar globalContErr = runTime.deltaTValue()* contErr.weightedAverage(mesh.V()).value(); cumulativeContErr += globalContErr; Info<< "time step continuity errors : sum local = " << sumLocalContErr << ", global = " << globalContErr << ", cumulative = " << cumulativeContErr << endl;
使用cavity进行测试发现此时的
time step continuity errors : sum local = 00222664,在2e-3左右,误差较大。
因此,是否是这样的phi的构造方式有问题呢?如何得到速度修正之后,最终满足连续性条件的phi呢 。 -
可以参考这个帖子 pisoFoam计算的U,求div(U)结果为何不是严格等于0
-
@李东岳
@coolhhh
感谢两位老师,我提出这个问题的主要原因是在现有的CFD-DEM流固耦合计算中,一种基于FD/IBM的细网格方法是需要对速度场U进行改变,并针对连续性方程进行修正的。而这个连续性修正的过程就需要再次显式构建连续性条件:div(U)=0。当然这是理论上的条件。
$\nabla \cdot (\nabla \phi) = \nabla \cdot U$
而下面代码中的连续性修正因子的phiIB理论上只在颗粒存在的条件下才有意义,同样才有值。但是如果div(U)在颗粒不存在得情况就不为0了,那么这样的修正就是错误的。虽然当前这个方面的主要算法就是下面展示的这样。目前看来在PISO外部,U和phi是分离的,但是CFD-DEM的修正是针对U的,这样操作应该是有问题的,误差大小有待考察。而我原本的目标是希望能够在PISO外,得到满足连续性的U或者phi,不过目前看来是很难实现的。
上述过程具体代码如下:
void Foam::cfdemCloud::calcVelocityCorrection ( volScalarField& p, volVectorField& U, volScalarField& phiIB, volScalarField& voidfraction//颗粒体积分数(颗粒内部为1 ) { void Foam::cfdemCloudIB::calcVelocityCorrection ( volScalarField& p, volVectorField& U, volScalarField& phiIB, volScalarField& voidfraction ) { setParticleVelocity(U);//改变颗粒所在区域流体的速度 // 修改速度会使其不满足连续性,因此使用一个phiIB进行修正 fvScalarMatrix phiIBEqn ( fvm::laplacian(phiIB) == fvc::div(U) + fvc::ddt(voidfraction) ); if(phiIB.needReference()) { phiIBEqn.setReference(pRefCell_, pRefValue_); } phiIBEqn.solve(); U=U-fvc::grad(phiIB); U.correctBoundaryConditions(); // correct the pressure as well p=p+phiIB/U.mesh().time().deltaT(); p.correctBoundaryConditions(); }