使用openfoam但不采用贴体网格计算圆柱绕流的方法，请各位同学老师指点一下，谢谢

OItoCFD

@李东岳李老师是的我尝试了体积力插到面上代码是这样的 f是算好的。主要是三句：

surfaceScalarField phiIB(fvc::flux(rAtU*f));

fvScalarMatrix pEqn
    (
        fvm::laplacian(rAtU(), p) == fvc::div(phiHbyA)  + fvc::div(phiIB)
    );

U = HbyA - rAtU*fvc::grad(p) + rAtU*f;

即插值得到phiIB，加到pEqn，然后下面是U = HbyA - rAtU*fvc::grad(p) + rAtU*f，这样做无法收敛，时间步会无限小下去。不知道错在哪了，这是详细的Peqn.C：

volScalarField rAU(1.0/UEqn.A());
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
surfaceScalarField phiHbyA
(
    "phiHbyA",
    fvc::flux(HbyA)
  + MRF.zeroFilter(fvc::interpolate(rAU)*fvc::ddtCorr(U, phi, Uf))
);

MRF.makeRelative(phiHbyA);

if (p.needReference())
{
    fvc::makeRelative(phiHbyA, U);
    adjustPhi(phiHbyA, U, p);
    fvc::makeAbsolute(phiHbyA, U);
}

tmp<volScalarField> rAtU(rAU);

if (pimple.consistent())
{
    rAtU = 1.0/max(1.0/rAU - UEqn.H1(), 0.1/rAU);
    phiHbyA +=
        fvc::interpolate(rAtU() - rAU)*fvc::snGrad(p)*mesh.magSf();
    HbyA -= (rAU - rAtU())*fvc::grad(p);
}
surfaceScalarField phiIB(fvc::flux(rAtU*f));



if (pimple.nCorrPiso() <= 1)
{
    tUEqn.clear();
}

// Update the pressure BCs to ensure flux consistency
constrainPressure(p, U, phiHbyA, rAtU(), MRF);

// Non-orthogonal pressure corrector loop
while (pimple.correctNonOrthogonal())
{
    Info << "!herehrehrehrehrheh" << nl;
    fvScalarMatrix pEqn
    (
        fvm::laplacian(rAtU(), p) == fvc::div(phiHbyA)  + fvc::div(phiIB)
    );

    pEqn.setReference(pRefCell, pRefValue);

    pEqn.solve();

    if (pimple.finalNonOrthogonalIter())
    {
        phi = phiHbyA - pEqn.flux();
    }
}

#include "continuityErrs.H"

// Explicitly relax pressure for momentum corrector
p.relax();

U = HbyA - rAtU*fvc::grad(p) + rAtU*f;
U.correctBoundaryConditions();
fvOptions.correct(U);

李东岳

@oitocfd

这样做无法收敛，时间步会无限小下去。

那就是压力速度耦合这面没处理好，我有的时候植入一个新算法也会这样。你这一段也应该考虑IB力：

if (pimple.finalNonOrthogonalIter())
    {
        phi = phiHbyA + fvc::div(phiIB) - pEqn.flux();
    }

OItoCFD

@李东岳好的感谢李老师我去尝试一下看看

OItoCFD

@李东岳
李老师按您说的尝试了下还是时间步无限减小

surfaceScalarField phiIB(fvc::flux(rAtU*f));
...

fvm::laplacian(rAtU(), p) == fvc::div(phiHbyA)  + fvc::div(phiIB)
...

phi = phiHbyA - pEqn.flux() + phiIB;
...

U = HbyA - rAtU*fvc::grad(p) +  rAtU*f;

然后只能算步长0.02，算2步到0.04秒就不收敛了，结果是错误的。全域一开始1m/s,入口1m/s,然后0.04s时速度和压强云图是这个样子，速度在圆柱一圈变得很大，到了20多米每秒，压强是一个很奇怪的完全一半一半对称的形状：

然后我也试着直接把f加到HbyA里。一样是不收敛的

volVectorField HbyA(constrainHbyA(rAU*UEqn.H() + rAU*f, U, p));

奇怪了，感觉加个f而已应该很简单才对呀。

李东岳

@oitocfd 你这个力是耦合在一起的，不是解耦的，所以肯定比加一个重力之类的麻烦，感觉还是方程推导/植入上不consistent，需要更深层次debug

OItoCFD

@李东岳好的李老师谢谢指点我再认真去思考下方法

李东岳

@oitocfd 你的第一种方法也是锯齿状的么

OItoCFD

@李东岳不是第一种方法用了个quadratic function插值还是比较光滑的

OItoCFD

@李东岳就是圆柱是光滑的IB点连接的速度是quadratic function从流体网格插值到光滑的一圈圆的散点上

Vortex

我也是用的direct-forcing immersed boundary method，方法跟你说的类似，不过我做的是大气边界层山脉的模拟。浸入边界交界面内外侧的应力项要处理好，我感觉有可能是你没处理好这一步。

OItoCFD

@vortex 可否请教下您的f是以哪种方式加的吗？耦合在peqn里面吗？

Vortex

@oitocfd 我没用OpenFOAM，用的求解器是lesgo，开源程序自带植入了浸没边界。它的实现我记得是没有耦合的，只在更新最后速度才加上IBM力，让固体内速度为0，详细的可参考这篇文章Modeling turbulent flow over fractal trees with renormalized numerical simulation。IBM的实现有很多变种的，有的确实是有耦合的。

OItoCFD

@vortex 是的我课题组的代码也是您说的这种是两部映射法里面中间步用f更新了速度再求压力泊松方程方程再更新最终速度然而改到openfoam发现出了点小问题

CFD中文网

使用openfoam但不采用贴体网格计算圆柱绕流的方法，请各位同学老师指点一下，谢谢