Coupled level set-VOF方法

[Zhujh]

各位大佬们好，我是研究波浪与结构物相互作用的，之前看到了下面这篇文章 https://www.sciencedirect.com/science/article/pii/S0301932206001832?via%3Dihub ，因为用VOF捕捉自由液面会存在钝化的问题，特别是想研究破碎波的时候，水气交界面的网格需要非常细，于是就想用一下这个方法。

根据这篇文章网上有公开的一部分代码，我就照着植入到了interFoam里面，取了个名字叫clsinterFoam，下面几张是我完成后对比两个求解器气泡上升的结果，可以看到clsinterFoam确实对交界面的捕捉更到位（左边是interFoam，右边是clsinterFoam）。

但是当我换到溃坝的算例的时候，出现了奇怪的问题，后处理结果来看，水体好像变得不连续了，最后会扩散到整个计算域。看上去像是没有了重力，但是算例里面我是加了的，想请教一下各位大佬，看看问题出在哪儿。

以下是在interFoam上改的部分代码：

	#include "mappingPsi.H"
	#include "solveLSFunction.H"
	#include "calcNewCurvature.H"

    turbulence->validate();

    if (!LTS)
    {
        #include "CourantNo.H"
        #include "setInitialDeltaT.H"
    }

    Info<< "\nStarting time loop\n" << endl;

    while (runTime.run())
    {
        #include "readDyMControls.H"

        if (LTS)
        {
            #include "setRDeltaT.H"
        }
        else
        {
            #include "CourantNo.H"
            #include "alphaCourantNo.H"
            #include "setDeltaT.H"
        }

        runTime++;

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

        // --- Pressure-velocity PIMPLE corrector loop
        while (pimple.loop())
        {
            if (pimple.firstIter() || moveMeshOuterCorrectors)
            {
                mesh.update();

                if (mesh.changing())
                {
                    // Do not apply previous time-step mesh compression flux
                    // if the mesh topology changed
                    if (mesh.topoChanging())
                    {
                        talphaPhi1Corr0.clear();
                    }

                    gh = (g & mesh.C()) - ghRef;
                    ghf = (g & mesh.Cf()) - ghRef;

                    MRF.update();

                    if (correctPhi)
                    {
                        // Calculate absolute flux
                        // from the mapped surface velocity
                        phi = mesh.Sf() & Uf();

                        #include "correctPhi.H"

                        // Make the flux relative to the mesh motion
                        fvc::makeRelative(phi, U);

                        mixture.correct();
                    }

                    if (checkMeshCourantNo)
                    {
                        #include "meshCourantNo.H"
                    }
                }
            }

            #include "alphaControls.H"
            #include "alphaEqnSubCycle.H"

			#include "mappingPsi.H"
			#include "solveLSFunction.H"
			#include "calcNewCurvature.H"
			#include "updateFlux.H"

			mixture.correct();

            if (pimple.frozenFlow())
            {
                continue;
            }

            #include "UEqn.H"

            // --- Pressure corrector loop
            while (pimple.correct())
            {
                #include "pEqn.H"
            }

            if (pimple.turbCorr())
            {
                turbulence->correct();
            }
        }

		// reInitialise the alpha equation 
		if (runTime.outputTime())
		{
			Info << "Overwriting alpha" << nl << endl;
			alpha1 = H;
			volScalarField& alpha10 = const_cast<volScalarField&>(alpha1.oldTime());
			alpha10 = H.oldTime();
			//const_cast<volScalarField&>(alpha1.storeOldTime()()) = H.oldTime();
		}

        runTime.write();

        runTime.printExecutionTime(Info);
    }

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

    return 0;

//mappingPsi.H
// mapping alpha value to psi0
   psi0 == (double(2.0)*alpha1-double(1.0))*gamma;

//solvevelSFunction.H
// solve Level-Set function as the re-initialization equation
   Info<< "solve the reinitialization equation"     
       << nl << endl;

   psi == psi0;

   for (int corr=0; corr<int(epsilon.value()/deltaTau.value()); corr++)
   {
      psi = psi + psi0/mag(psi0)*(double(1)-mag(fvc::grad(psi)*dimChange))*deltaTau;
      psi.correctBoundaryConditions();
   }


// update Dirac function
   forAll(mesh.cells(),celli)
   {
      if(mag(psi[celli]) > epsilon.value())
         delta[celli] = double(0);
      else
         delta[celli] = double(1.0)/(double(2.0)*epsilon.value())*(double(1.0)+Foam::cos(M_PI*psi[celli]/epsilon.value()));
   };

// update Heaviside function
   forAll(mesh.cells(),celli)
   {
      if(psi[celli] < -epsilon.value())
         H[celli] = double(0);
      else if(epsilon.value() < psi[celli])
         H [celli] = double(1);
      else
         H[celli] = double(1.0)/double(2.0)*(double(1.0)+psi[celli]/epsilon.value()+Foam::sin(M_PI*psi[celli]/epsilon.value())/M_PI);
   };

//calcNewCurvature.H
// calculate normal vector
volVectorField gradPsi(fvc::grad(psi));
surfaceVectorField gradPsif(fvc::interpolate(gradPsi));
surfaceVectorField nVecfv(gradPsif/(mag(gradPsif)+scalar(1.0e-6)/dimChange));
surfaceScalarField nVecf(nVecfv & mesh.Sf());

// calculate new curvature based on psi (LS function)
   C == -fvc::div(nVecf);

//updataFlux.H
{
    word alphaScheme("div(phi,alpha)");
    word alpharScheme("div(phirb,alpha)");
	// Standard face-flux compression coefficient
    surfaceScalarField phic(mixture.cAlpha()*mag(phi/mesh.magSf()));
	surfaceScalarField::Boundary& phicBf = phic.boundaryFieldRef();
    //surfaceScalarField phic(mag(phi/mesh.magSf()));

    // Add the optional isotropic compression contribution
    if (icAlpha > 0)
    {
        phic *= (1.0 - icAlpha);
        phic += (mixture.cAlpha()*icAlpha)*fvc::interpolate(mag(U));
    }

    // Do not compress interface at non-coupled boundary faces
    // (inlets, outlets etc.)
    forAll(phic.boundaryField(), patchi)
    {
        fvsPatchScalarField& phicp = phicBf[patchi];

        if (!phicp.coupled())
        {
            phicp == 0;
        }
    }

    surfaceScalarField phir(phic*nVecf);

    surfaceScalarField phiAlpha
    (
        fvc::flux
        (
            phi,
            alpha1,
            alphaScheme
        )
      + fvc::flux
        (
            -fvc::flux(-phir, scalar(1) - alpha1, alpharScheme),
            alpha1,
            alpharScheme
        )
    );

        MULES::explicitSolve(alpha1, phi, phiAlpha, 1, 0);

        rhoPhi = phiAlpha*(rho1 - rho2) + phi*rho2;
}

[李东岳]

@zhujh 看起来非常有意思。你那个算的比VOF强很多。不过我要收拾行李去杭州了。只能国庆回来能详细看看了。但愿其他大佬能帮到你。比如 @队长别开枪

[Zhujh]

@李东岳 谢谢谢谢李老师！！这个是我20年的暑假整的，但是当时没整出来，今年才发了小论文，下午注册的账号，我可太激动了！马上就发了一贴。
补充一下就是气泡那个模拟，两个求解器是用的同一套网格。