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年的暑假整的，但是当时没整出来，今年才发了小论文，下午注册的账号，我可太激动了！马上就发了一贴。
补充一下就是气泡那个模拟，两个求解器是用的同一套网格。

[Shihang Chen]

@Zhujh 您好，也采用了与您基本相同的CLSVOF方法进行计算，然而对于毛细张力主导的问题，这个方法表现出了更强的虚假流动（寄生流动）的问题，如无重力状态下的水滴。请问您遇到过类似的问题么？同时我也参考fluent里面的两种抑制虚假流动方法（密度和H函数），但是收效甚微，想问下您有什么建议吗

[EricLiu]

@Zhujh @Shihang-Chen 两位大佬好，请问你们用到的CLSVOF 方法有相关的代码资源吗，我最近也打算结合CLS和VOF方法，但是没有找到合适的切入点，openfoam这块属于刚开始学

[Shihang Chen]

@EricLiu 不是大佬。这个讨论顶部文章附带一些代码，或者直接github搜索clsvof也能有很多收获。

[EricLiu]

@Shihang-Chen 十分感谢，我去查一下

[capillaryFix]

@Shihang-Chen 关于spurious currents/velocities，我们此前在OpenFOAM里植入了一个非常简单但是可行的方法，代码量与植入难度不大，具体请参考：https://www.sciencedirect.com/science/article/pii/S1877750323002557

[Shihang Chen]

@capillaryFix 太好了，感谢感谢，我会好好看看。我目前一直使用的消除寄生流的方法是基于sharpen surface model来的，这个效果也是很好。

[李东岳]

@capillaryFix 这个流弊，老铁