关于interFoam 中的alphaEqn里面的两个系数cAlpha和icAlpha



  • 版本3.0.x中的interFoam里面的alphaEqn.H line46-54

    // Standard face-flux compression coefficient
    surfaceScalarField phic(mixture.cAlpha()*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));
    }
    

    查看了一下源文件, 这两个系数都是用户可以指定的,小弟不太明白,这两个系数的作用的什么?有知道的朋友告诉一下,谢谢!



  • 麻烦把所有alphaeqn.H里面的内容发一下:confused:



  • @cfd-china alphaeqn.H里面的内容如下

    {
        word alphaScheme("div(phi,alpha)");
        word alpharScheme("div(phirb,alpha)");
    
        tmp<fv::ddtScheme<scalar> > ddtAlpha
        (
            fv::ddtScheme<scalar>::New
            (
                mesh,
                mesh.ddtScheme("ddt(alpha)")
            )
        );
    
        // Set the off-centering coefficient according to ddt scheme
        scalar ocCoeff = 0;
        if
        (
            isType<fv::EulerDdtScheme<scalar> >(ddtAlpha())
         || isType<fv::localEulerDdtScheme<scalar> >(ddtAlpha())
        )
        {
            ocCoeff = 0;
        }
        else if (isType<fv::CrankNicolsonDdtScheme<scalar> >(ddtAlpha()))
        {
            if (nAlphaSubCycles > 1)
            {
                FatalErrorIn(args.executable())
                    << "Sub-cycling is not supported "
                       "with the CrankNicolson ddt scheme"
                    << exit(FatalError);
            }
    
            ocCoeff =
                refCast<fv::CrankNicolsonDdtScheme<scalar> >(ddtAlpha()).ocCoeff();
        }
        else
        {
            FatalErrorIn(args.executable())
                << "Only Euler and CrankNicolson ddt schemes are supported"
                << exit(FatalError);
        }
    
        scalar cnCoeff = 1.0/(1.0 + ocCoeff);
    
        // Standard face-flux compression coefficient
        surfaceScalarField phic(mixture.cAlpha()*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 = phic.boundaryField()[patchi];
    
            if (!phicp.coupled())
            {
                phicp == 0;
            }
        }
    
        tmp<surfaceScalarField> phiCN(phi);
    
        // Calculate the Crank-Nicolson off-centred volumetric flux
        if (ocCoeff > 0)
        {
            phiCN = cnCoeff*phi + (1.0 - cnCoeff)*phi.oldTime();
        }
    
        if (MULESCorr)
        {
            fvScalarMatrix alpha1Eqn
            (
                (
                    LTS
                  ? fv::localEulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
                  : fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
                )
              + fv::gaussConvectionScheme<scalar>
                (
                    mesh,
                    phiCN,
                    upwind<scalar>(mesh, phiCN)
                ).fvmDiv(phiCN, alpha1)
            );
    
            alpha1Eqn.solve();
    
            Info<< "Phase-1 volume fraction = "
                << alpha1.weightedAverage(mesh.Vsc()).value()
                << "  Min(" << alpha1.name() << ") = " << min(alpha1).value()
                << "  Max(" << alpha1.name() << ") = " << max(alpha1).value()
                << endl;
    
            tmp<surfaceScalarField> talphaPhiUD(alpha1Eqn.flux());
            alphaPhi = talphaPhiUD();//based on updated alpha from alpha1Eqn.solve() to correct alphaPhi
    
            if (alphaApplyPrevCorr && talphaPhiCorr0.valid())
            {
                Info<< "Applying the previous iteration compression flux" << endl;
                MULES::correct(alpha1, alphaPhi, talphaPhiCorr0(), 1, 0);
    
                alphaPhi += talphaPhiCorr0();
            }
    
            // Cache the upwind-flux
            talphaPhiCorr0 = talphaPhiUD;
    
            alpha2 = 1.0 - alpha1;
    
            mixture.correct();
        }
    
    
        for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
        {
            surfaceScalarField phir(phic*mixture.nHatf());
    
            tmp<surfaceScalarField> talphaPhiUn
            (
                fvc::flux
                (
                    phi,
                    alpha1,
                    alphaScheme
                )
              + fvc::flux
                (
                   -fvc::flux(-phir, alpha2, alpharScheme),
                    alpha1,
                    alpharScheme
                )
            );
    
            // Calculate the Crank-Nicolson off-centred alpha flux
            if (ocCoeff > 0)
            {
                talphaPhiUn =
                    cnCoeff*talphaPhiUn + (1.0 - cnCoeff)*alphaPhi.oldTime();
            }
    
            if (MULESCorr)
            {
                tmp<surfaceScalarField> talphaPhiCorr(talphaPhiUn() - alphaPhi);
                volScalarField alpha10("alpha10", alpha1);
    
                MULES::correct(alpha1, talphaPhiUn(), talphaPhiCorr(), 1, 0);
    
                // Under-relax the correction for all but the 1st corrector
                if (aCorr == 0)
                {
                    alphaPhi += talphaPhiCorr();
                }
                else
                {
                    alpha1 = 0.5*alpha1 + 0.5*alpha10;
                    alphaPhi += 0.5*talphaPhiCorr();
                }
            }
            else
            {
                alphaPhi = talphaPhiUn;
    
                MULES::explicitSolve(alpha1, phiCN, alphaPhi, 1, 0);
            }
    
            alpha2 = 1.0 - alpha1;
    
            mixture.correct();
        }
    
        if (alphaApplyPrevCorr && MULESCorr)
        {
            talphaPhiCorr0 = alphaPhi - talphaPhiCorr0;
        }
    
        if
        (
            word(mesh.ddtScheme("ddt(rho,U)"))
         == fv::EulerDdtScheme<vector>::typeName
        )
        {
            rhoPhi = alphaPhi*(rho1 - rho2) + phiCN*rho2;
        }
        else
        {
            if (ocCoeff > 0)
            {
                // Calculate the end-of-time-step alpha flux
                alphaPhi = (alphaPhi - (1.0 - cnCoeff)*alphaPhi.oldTime())/cnCoeff;
            }
    
            // Calculate the end-of-time-step mass flux
            rhoPhi = alphaPhi*(rho1 - rho2) + phi*rho2;
        }
    
        Info<< "Phase-1 volume fraction = "
            << alpha1.weightedAverage(mesh.Vsc()).value()
            << "  Min(" << alpha1.name() << ") = " << min(alpha1).value()
            << "  Max(" << alpha1.name() << ") = " << max(alpha1).value()
            << endl;
    }
    


  • 在论文DYNAMIC MESHING AROUND FLUID-FLUID INTERFACES WITH APPLICATIONS TO DROPLET TRACKING IN CONTRACTION GEOMETRIES第58页最上部有写cAlpha的作用。

    The cAlpha keyword specified in the file <case>/system/fvSolution is
    a factor that controls the compression of the interface, where 0 corresponds to no
    compression. After solving the volume fluid equation, the density and viscosity are
    modified using the new values of α.



  • @supersoldier thanks



  • @supersoldier 非常有用的回复,直接帮我解决了很多疑问。谢谢!!



  • @supersoldier 你好,可以再详细解释一下不?


登录后回复
 

与 CFD中文网 的连接断开,我们正在尝试重连,请耐心等待