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

[OItoCFD]

各位老师和同学好。之前一直使用课题组代码，最近在转向学习OpenFOAM，在做一个小东西来练习，想不采用贴体网格，而是浸入边界法或者cut cell的方法来计算圆柱绕流。

1. 方法1：浸入边界法。拟采用论文1：Ji C, Munjiza A, Williams J J R. A novel iterative direct-forcing immersed boundary method and its finite volume applications[J]. Journal of Computational Physics, 2012, 231(4): 1797-1821.

2. 方法2：cut cell方法。拟采用论文2：Lin P. A fixed-grid model for simulation of a moving body in free surface flows[J]. Computers & fluids, 2007, 36(3): 549-561.

首先采用方法1：论文1中的方法是一个direct-forcing immersed boundary method，固体的作用通过体积力施加给NS方程，固体内部相当于存在虚拟流体，虚拟流体和固体边界外流体刚好在固体边界处满足无滑移边界条件。具体的实现在论文1中是基于两步映射法，在第一步求得中间速度后将流体速度插值到固体点边界，固体速度与流体插值过来的速度应相等满足无滑移，则可计算出速度差得额外体积力f，并将该体积力通过分布函数施加给邻近的影响的流体网格，更新中间速度，然后采用中间速度求压力泊松方程得到该步压强再更新最终速度。(注：论文中增加了一个迭代循环，但是循环1次也即等同于传统直接应力的浸入边界法)


那么在OpenFOAM中，采用的piso或pimple算法。我应该在每一步将额外体积力f施加在什么部分呢？
首先我试了采用上一步计算得到的流体速度来插值到固体边界先求出f，然后施加在Ueqn里

//do something to calculate f using U of last step, 
//then :

tmp<fvVectorMatrix> tUEqn
(
    fvm::ddt(U) + fvm::div(phi, U)
  + MRF.DDt(U)
  + turbulence->divDevSigma(U)
 ==
    fvOptions(U) + f
);

用这种方法计算无法收敛。
然后我在UEqn之后，采用该中间速度算出f，施加在momentumPredictor的部分：

tmp<fvVectorMatrix> tUEqn
(
    fvm::ddt(U) + fvm::div(phi, U)
  + MRF.DDt(U)
  + turbulence->divDevSigma(U)
 ==
    fvOptions(U)
);
fvVectorMatrix& UEqn = tUEqn.ref();

UEqn.relax();

fvOptions.constrain(UEqn);

//do something to calculate f using the intermediate U , 
//then :

if (pimple.momentumPredictor())
{
    solve(UEqn == -fvc::grad(p) + f);

    fvOptions.correct(U);
}

这种方法计算是可以收敛的，而且可以正常算出涡，我计算的是一个Re=100的圆柱绕流，但是这个涡是不对的，和正确位置比滞后了很多，而且拖曳力系数是正确拖曳力系数的好几倍。并且这种施加方法，结果会和piso或pimple算法的correct步数相关了(nCorrectors)。correct步数次数在openfoam中默认2次，可以增加次数，当correct次数改变时候，我计算得到的拖曳力系数成倍增长，涡的形态也发生变化。
这是Re=100圆柱绕流正确的形态，为我使用课题组代码计算得到的：

这是当前我采用Openfoam利用我上述施加f在momentumPredictor()部分的方法算出的结果：

涡滞后了很多，而且拖曳力大了几倍甚至几十倍，随着correct次数的增加而增加。

在谷歌一些相关资料后，很多方法是将f施加在Peqn里，看起来类似interFoam里重力施加的方式，初步尝试总是不收敛，暂时我应该遗漏了某些细节。

除此之外，也尝试了方法2的cut cell方法：在论文2采用一个系数表征固体，流体和流固交界cell。固体部分系数为0，流体系数为1，相交部分为0-1，这类似VOF的方法。论文2核心其实是在当固体运动时采用一个相对速度来修正连续性方程保证质量守恒。我假设最最简单的情况，忽略部分切割(系数在0-1之间)的网格，且固体不动，就暂时不需要论文2对连续性方程的修正部分了。同样计算Re=100的圆柱绕流。按照论文2的方法，完全被固体占据的网格（系数=0）不参与计算要跳过（如果不跳过，会出现除以0的情况）。为了实现固体完全占据的部分不参与计算且速度设为0（有些等同于于扣掉网格了），我采用了这个帖子中@马乔 的方法：https://www.cfd-china.com/topic/3138/openfoam有方法能够使一部分网格不参与计算吗/28

\*---------------------------------------------------------------------------*/

#include "fvCFD.H"
#include "pisoControl.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])
{
    //Info<< "\nStarting time loop\n" << endl;
    #include "setRootCaseLists.H"
    #include "createTime.H"
    #include "createMesh.H"

    pisoControl piso(mesh);

    #include "createFields.H"
    #include "initContinuityErrs.H"

    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

    Info<< "\nStarting time loop\n" << endl;
    
    DynamicList<label> lowerFaceCells(mesh.C().size());
    DynamicList<label>  upperFaceCells(mesh.C().size());
    //mesh LDU address
    const labelUList& lowAddr = mesh.lduAddr().lowerAddr();
    const labelUList& upAddr = mesh.lduAddr().upperAddr();
    //label refCell = 4;
    //labelList cellBlock (20);
    //DynamicList<label> cellBlock;
    if(mesh.cellZones().size())
    {
	const labelList& cellBlock = mesh.cellZones()[0];
	Info << "cell block size: " << cellBlock.size() << endl;
	//Info << cellBlock << endl;

	//low efficient
	forAll(lowAddr, faceI)
	{
	    for(auto cellI : cellBlock)
	    {
		if(lowAddr[faceI] == cellI)
		{
		    lowerFaceCells.append(faceI);
		}

		if(upAddr[faceI] == cellI)
		{
		    upperFaceCells.append(faceI);
		}
	    }
	}
    }

    lowerFaceCells.shrink();
    upperFaceCells.shrink();

    Info << "low Cells size: " << lowerFaceCells.size() << ", low Cells size: " << upperFaceCells.size() << endl;

    while (runTime.loop())
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;

        #include "CourantNo.H"

        // Momentum predictor

        fvVectorMatrix UEqn
        (
            fvm::ddt(U)
          + fvm::div(phi, U)
          - fvm::laplacian(nu, U)
        );
        
        /*label refCell = 4;
        labelList refCells (1,refCell);
        vectorField refValues (1, vector(0,0,0));
        UEqn.setValues(refCells, refValues);*/
        
        
//modify matrix
	if(mesh.cellZones().size())
	{
	    const cellZone& cellBlock = mesh.cellZones()[0]; //should select cellZone accor name

	    scalarField& lower = UEqn.lower();
	    scalarField& upper = UEqn.upper();
	    scalarField& diag = UEqn.diag();


	    for(auto cellI : cellBlock)
	    {
		diag[cellI] = 1.;
		UEqn.source()[cellI] = vector::zero; //vector(0.,0.,0.); you can set your deserved value here, or read from dict, ref. fvoptions
	    }

	    for(auto faceI : lowerFaceCells)
	    {
	 	lower[faceI] = 0.;	
	    }

	    for(auto faceI : upperFaceCells)
	    {
	 	upper[faceI] = 0.;	
	    }
	}
	else
	{
	    Info << "There is no cell selected" << endl;
	}
        

        if (piso.momentumPredictor())
        {
            solve(UEqn == -fvc::grad(p));
        }

        // --- PISO loop
        while (piso.correct())
        {
            volScalarField rAU(1.0/UEqn.A());
            volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
            surfaceScalarField phiHbyA
            (
                "phiHbyA",
                fvc::flux(HbyA)
              + fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)
            );

            adjustPhi(phiHbyA, U, p);

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

            // Non-orthogonal pressure corrector loop
            while (piso.correctNonOrthogonal())
            {
                // Pressure corrector

                fvScalarMatrix pEqn
                (
                    fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
                );

                pEqn.setReference(pRefCell, pRefValue);

		//modify p matrix
		if(mesh.cellZones().size())
		{
		    const cellZone& cellBlock = mesh.cellZones()[0];
		    //p matrix symmetry
		    scalarField& offDiag = pEqn.hasUpper() ? pEqn.upper() : pEqn.lower();
		    scalarField& diag = pEqn.diag();

		    for(auto cellI : cellBlock)
		    {
			diag[cellI] = 1.;
			pEqn.source()[cellI] = 0.; 
		    }

		    //for(auto faceI : upperFaceCells) //should find faces again?
		    //{
		    //    offDiag[faceI] = 0.;	
		    //}
		    const labelUList& offAddr = pEqn.hasUpper() ? pEqn.lduAddr().upperAddr() : pEqn.lduAddr().lowerAddr();

		    forAll(offAddr, faceI)
		    {
			for(auto cellI : cellBlock)
			{
			    if(offAddr[faceI] == cellI)
			    {
				offDiag[faceI] = 0.;	
			    }
			}
		    }
		}
                pEqn.solve();

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

            #include "continuityErrs.H"

            U = HbyA - rAU*fvc::grad(p);
            U.correctBoundaryConditions();
        }

        runTime.write();
        
        //const vectorField& UCells = U.internalField();
        
        
        
        //dimensionedVector v1 ("v1", dimVelocity, vector (1.0,0.0,0.0));
        
        //U.internalField()[4] = v1;
        
        /*Info << U << nl << endl;
        Info << mesh.C() << nl << endl;*/

        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
            << nl << endl;
    }

    Info<< "End\n" << endl;

    return 0;
}

确实可以计算，但是会一开始出现一点点涡，同样也非常滞后看起来不正常。然后随着计算时间的增加，涡完全消失不产生绕流。不知是否在这种情况还需要对跳过cell处的边界做一些特殊处理？下面的图是刚计算了几秒的时候，看起来是正常的，算一段时间后产生涡，同样滞后。再算一段时间，涡消失完全没涡了。

小弟接触openfoam时间不长，希望各位老师和同学能帮忙给点建议或者方向。解决这个采用IBM或cut cell来计算圆柱绕流的问题。不采用贴体网格是因为我后续要做的东西更适合和IBM或cut cell相结合。谢谢大家！

[李东岳]

@oitocfd 你算的圆柱不是个圆，是一个带锯齿的东西？

[OItoCFD]

@李东岳 李老师 方法2我为了先试简单的 就没把固体部分占据流体的地方考虑 不考虑部分切割 就等于锯齿状的 因为想先看效果 就算是锯齿状的 但是整体状态上应该大致有对的样子才对 大致对了再把部分切割地方的系数乘到方程里 现在等于算系数全为0的部分即完全被固体占据的地方和系数为1也就是只有水未被固体占据的部分

[李东岳]

@oitocfd

在谷歌一些相关资料后，很多方法是将f施加在Peqn里，看起来类似interFoam里重力施加的方式，初步尝试总是不收敛，暂时我应该遗漏了某些细节。

我应该在每一步将额外体积力f施加在什么部分呢？

一般是简单的把体积力插值到面上，然后加到phiHbyA里面就好。你是怎么加的

[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发现出了点小问题