Navigation

    CFD中文网

    CFD中文网

    • Login
    • Search
    • 最新

    of中的多相流求解器

    OpenFOAM
    10
    33
    21547
    Loading More Posts
    • Oldest to Newest
    • Newest to Oldest
    • Most Votes
    Reply
    • Reply as topic
    Log in to reply
    This topic has been deleted. Only users with topic management privileges can see it.
    • M
      mohui last edited by

      在of中我有个疑问,关于多相流求解器。目前来看,两相流求解器,对于求解alpha方程都是第二相用1-alpha1得到,这样子虽然保证体积守恒怎么保证质量守恒,特别是有相变的求解器(interphasechangeFoam)。另外of中多相流求解器中的alpha方程是套了循环,也就是说每一相都求解了,那么最后是怎么保证每一个cell中所有的相的体积分数总和为1?

      1 Reply Last reply Reply Quote
      • C
        CFD中文网 last edited by

        $\alpha_1=1-\alpha_2$是用来保证alpha有界的,因为:
        \begin{equation}
        0<\alpha<1
        \end{equation}
        守恒体现在连续性方程中:
        \begin{equation}
        \frac{\partial \rho\alpha}{\partial t}+\nabla\cdot(\rho\mathbf{U}\alpha)=0
        \end{equation}

        of中多相流求解器中的alpha方程是套了循环

        具体指的是什么?MULES么?

        CFD中国标准用户测试帐号
        目前由徐笑笑登录

        M 1 Reply Last reply Reply Quote
        • M
          mohui @CFD中文网 last edited by

          @cfd-china OF中的多相流求解器是求解n相,所以需要套循环,alpha求解是用mules求解,然后我仔细看了多相流求解器,发现除了mules求解每一相之外,还有mules limitsum修正,这个不知道我理解对不对。理由是,我做了个测试,在compressibleinterFoam求解器中,对应的第二相也用mules求解,但是无法保证总的体积分数为1,算到一定时间算崩了。所以我推断多相流求解中的mules limitsum是保证总的体积分数守恒,关于这个理解不知道对不对,求解答。谢谢

          1 Reply Last reply Reply Quote
          • 队长别开枪
            队长别开枪 教授 last edited by

            VOF方程求解分为几何重构和纯代数方法两类,OpenFOAM中的MULES求解方法属于纯代数方法,具有效率高、易于实现、适用于任意多面体网格等优点,类似的还有CICSAM方法,代数类方法的求解重点在于计算单元表面体积流率和克服数值弥散。据我所知,在interFoam和interPhaseChangeFoam中都是只求解主相alpha1而且都只支持两相流动。

            M 1 Reply Last reply Reply Quote
            • M
              mohui @队长别开枪 last edited by

              @队长别开枪 现在我就想确定,alpha2通过对应方程推导,用MULES求解不可行的吗?

              队长别开枪 1 Reply Last reply Reply Quote
              • 队长别开枪
                队长别开枪 教授 @mohui last edited by

                @mohui 这个有点复杂,等我今天开完组会我们一起讨论一下。

                M 1 Reply Last reply Reply Quote
                • M
                  mohui @队长别开枪 last edited by

                  @队长别开枪 恩,我做了测试,如果第二相采用mules求解会导致alphasum无界。

                  李东岳 队长别开枪 2 Replies Last reply Reply Quote
                  • 李东岳
                    李东岳 管理员 @mohui last edited by

                    @mohui
                    如果alpha1和alpha2采用完全相同的数值处理,并且边界条件给定的互相对应,求解alpha1和alpha2的方程和求解alpha1方程并alpha2=1-alpha1应该是相同的。

                    可以贴出来你做的代码改动,你用什么算例测试的?

                    CFD高性能服务器 http://dyfluid.com/servers.html

                    M 1 Reply Last reply Reply Quote
                    • M
                      mohui @李东岳 last edited by 李东岳

                      @李东岳 ,首先我是看了您写的关于compressibleInterFoam求解器解析,然后推导完之后发现了个小小的错误。dgdt里面是没有alpha的。
                      然后这个是我推导的可压缩两相方程:
                      0_1493962493553_2017-05-05_13-33-57.bmp ,
                      根据这个方程结合源码,这是修改compressibleInterFoam中的alpha方程的代码:

                      for (int gCorr=0; gCorr<nAlphaCorr; gCorr++)
                          {
                              volScalarField::Internal Sp
                              (
                                  IOobject
                                  (
                                      "Sp",
                                      runTime.timeName(),
                                      mesh
                                  ),
                                  mesh,
                                  dimensionedScalar("Sp", dgdt.dimensions(), 0.0)
                              );
                      
                              volScalarField::Internal Su
                              (
                                  IOobject
                                  (
                                      "Su",
                                      runTime.timeName(),
                                      mesh
                                  ),
                                 
                                  divU*min(alpha1, scalar(1))
                              );
                      
                              forAll(dgdt, celli)
                              {
                                  if (dgdt[celli] > 0.0 && alpha1[celli] > 0.0)
                                  {
                                      Sp[celli] -= dgdt[celli]*alpha1[celli];
                                      Su[celli] += dgdt[celli]*alpha1[celli];
                                  }
                                  else if (dgdt[celli] < 0.0 && alpha1[celli] < 1.0)
                                  {
                                      Sp[celli] += dgdt[celli]*(max(1.0 - alpha1[celli],scalar(0)));
                                  }
                              }
                      
                      
                              surfaceScalarField alphaPhi1
                              (
                                  fvc::flux
                                  (
                                      phi,
                                      alpha1,
                                      alphaScheme
                                  )
                                + fvc::flux
                                  (
                                      -fvc::flux(-phir, alpha2, alpharScheme),
                                      alpha1,
                                      alpharScheme
                                  )
                              );
                      
                              MULES::explicitSolve
                              (
                                  geometricOneField(),
                                  alpha1,
                                  phi,
                                  alphaPhi1,
                                  Sp,
                                  Su,
                                  1,
                                  0
                              );
                      
                              //slove for alpha2
                              volScalarField::Internal Sp1
                              (
                                  IOobject
                                  (
                                      "Sp1",
                                      runTime.timeName(),
                                      mesh
                                  ),
                                  mesh,
                                  dimensionedScalar("Sp1", dgdt.dimensions(), 0.0)
                              );
                      
                              volScalarField::Internal Su1
                              (
                                  IOobject
                                  (
                                      "Su1",
                                      runTime.timeName(),
                                      mesh
                                  ),
                                  
                                  divU*min(alpha2, scalar(1))
                              );
                             
                              forAll(dgdt, celli)
                              {
                                  if (dgdt[celli] > 0.0 && alpha1[celli] > 0.0)
                                  {
                                      Sp1[celli] -= dgdt[celli]*alpha1[celli];
                                      
                                  }
                                  else if (dgdt[celli] < 0.0 && alpha1[celli] < 1.0)
                                  {
                                      Sp1[celli] += dgdt[celli]*alpha2[celli];
                                      Su1[celli] -= dgdt[celli]*(max(1.0 - alpha1[celli],scalar(0)));
                                  }
                              }
                               surfaceScalarField alphaPhi2
                              (
                                  fvc::flux
                                  (
                                      phi,
                                      alpha2,
                                      alphaScheme
                                  )
                                + fvc::flux
                                  (
                                      -fvc::flux(phir, alpha1, alpharScheme),
                                      alpha2,
                                      alpharScheme
                                  )
                              );
                               MULES::explicitSolve
                              (
                                  geometricOneField(),
                                  alpha2,
                                  phi,
                                  alphaPhi2,
                                  Sp1,
                                  Su1,
                                  1,
                                  0
                              );
                      
                      1 Reply Last reply Reply Quote
                      • 李东岳
                        李东岳 管理员 last edited by

                        可以把我那个错误得地方指出来么我好改一下。

                        我看了下你的方程和贴出来的代码没看出问题。能否把alpha1方程去掉,单独求解alpha2? 另外你的错误结果是什么?alpha1+alpha2>1么?还是结果严重不符合物理。

                        CFD高性能服务器 http://dyfluid.com/servers.html

                        M 1 Reply Last reply Reply Quote
                        • M
                          mohui @李东岳 last edited by

                          @李东岳 0_1493990651369_2017-05-05_21-16-12.bmp ,我说的错误是这个,里面是没有alpha1,alpha2的。对,算的结果是alpha1+alpha2最大值为2.。0_1493994510812_2017-05-05_22-06-09_副本.jpg 左边是原程序,右边是我改的程序。两者还是有点区别的。

                          李东岳 1 Reply Last reply Reply Quote
                          • 队长别开枪
                            队长别开枪 教授 @mohui last edited by

                            @mohui Hello, 不好意思现在才回复你。OpenFOAM的MULES求解器我不大熟悉,我自己本身是做多面体网格上的PLIC算法的,最后评估算法优劣的参数中就有一个mass/volume conservation error,这个值越小越好,如果算法能够保证这个参数恒等于零,那么分别求解alpha1和alpha2没有什么问题,但现实中算法误差和数值误差的存在,这个值不可能恒等于零,所以一般次相alpha2直接用1-alpha1计算,如果这个时候分别计算alpha1和alpha2,我觉得会叠加误差,以至于出现你遇到的alpha1+alpha2>1的情况,更新密度粘度场后误差再传进动量方程和压力泊松方程,最后造成求解发散。这些是我的看法。

                            M 3 Replies Last reply Reply Quote
                            • M
                              mohui @队长别开枪 last edited by

                              @队长别开枪 恩,确实是这样子的。

                              1 Reply Last reply Reply Quote
                              • M
                                mohui @队长别开枪 last edited by 李东岳

                                @队长别开枪 这个是compressiblemultiphaseFoam的关于求解alpha的源码,从源码看,除了用explicitsolve,在这前还使用了limit以及limitsum 限制,我在想是不是这样子才会使得求解每一相同时还能保证alpha之和为1.

                                void Foam::multiphaseMixtureThermo::solveAlphas
                                (
                                    const scalar cAlpha
                                )
                                {
                                    static label nSolves=-1;
                                    nSolves++;
                                
                                    word alphaScheme("div(phi,alpha)");
                                    word alpharScheme("div(phirb,alpha)");
                                
                                    surfaceScalarField phic(mag(phi_/mesh_.magSf()));
                                    phic = min(cAlpha*phic, max(phic));
                                
                                    PtrList<surfaceScalarField> phiAlphaCorrs(phases_.size());
                                    int phasei = 0;
                                
                                    forAllIter(PtrDictionary<phaseModel>, phases_, phase)
                                    {
                                        phaseModel& alpha = phase();
                                
                                        phiAlphaCorrs.set
                                        (
                                            phasei,
                                            new surfaceScalarField
                                            (
                                                fvc::flux
                                                (
                                                    phi_,
                                                    alpha,
                                                    alphaScheme
                                                )
                                            )
                                        );
                                
                                        surfaceScalarField& phiAlphaCorr = phiAlphaCorrs[phasei];
                                
                                        forAllIter(PtrDictionary<phaseModel>, phases_, phase2)
                                        {
                                            phaseModel& alpha2 = phase2();
                                
                                            if (&alpha2 == &alpha) continue;
                                
                                            surfaceScalarField phir(phic*nHatf(alpha, alpha2));
                                
                                            phiAlphaCorr += fvc::flux
                                            (
                                                -fvc::flux(-phir, alpha2, alpharScheme),
                                                alpha,
                                                alpharScheme
                                            );
                                        }
                                
                                        MULES::limit
                                        (
                                            1.0/mesh_.time().deltaT().value(),
                                            geometricOneField(),
                                            alpha,
                                            phi_,
                                            phiAlphaCorr,
                                            zeroField(),
                                            zeroField(),
                                            1,
                                            0,
                                            3,
                                            true
                                        );
                                
                                        phasei++;
                                    }
                                
                                    MULES::limitSum(phiAlphaCorrs);
                                
                                    rhoPhi_ = dimensionedScalar("0", dimensionSet(1, 0, -1, 0, 0), 0);
                                
                                    volScalarField sumAlpha
                                    (
                                        IOobject
                                        (
                                            "sumAlpha",
                                            mesh_.time().timeName(),
                                            mesh_
                                        ),
                                        mesh_,
                                        dimensionedScalar("sumAlpha", dimless, 0)
                                    );
                                
                                
                                    volScalarField divU(fvc::div(fvc::absolute(phi_, U_)));
                                
                                
                                    phasei = 0;
                                
                                    forAllIter(PtrDictionary<phaseModel>, phases_, phase)
                                    {
                                        phaseModel& alpha = phase();
                                
                                        surfaceScalarField& phiAlpha = phiAlphaCorrs[phasei];
                                        phiAlpha += upwind<scalar>(mesh_, phi_).flux(alpha);
                                
                                        volScalarField::DimensionedInternalField Sp
                                        (
                                            IOobject
                                            (
                                                "Sp",
                                                mesh_.time().timeName(),
                                                mesh_
                                            ),
                                            mesh_,
                                            dimensionedScalar("Sp", alpha.dgdt().dimensions(), 0.0)
                                        );
                                
                                        volScalarField::DimensionedInternalField Su
                                        (
                                            IOobject
                                            (
                                                "Su",
                                                mesh_.time().timeName(),
                                                mesh_
                                            ),
                                            // Divergence term is handled explicitly to be
                                            // consistent with the explicit transport solution
                                            divU*min(alpha, scalar(1))
                                        );
                                
                                        {
                                            const scalarField& dgdt = alpha.dgdt();
                                
                                            forAll(dgdt, celli)
                                            {
                                                if (dgdt[celli] < 0.0 && alpha[celli] > 0.0)
                                                {
                                                    Sp[celli] += dgdt[celli]*alpha[celli];
                                                    Su[celli] -= dgdt[celli]*alpha[celli];
                                                }
                                                else if (dgdt[celli] > 0.0 && alpha[celli] < 1.0)
                                                {
                                                    Sp[celli] -= dgdt[celli]*(1.0 - alpha[celli]);
                                                }
                                            }
                                        }
                                
                                        forAllConstIter(PtrDictionary<phaseModel>, phases_, phase2)
                                        {
                                            const phaseModel& alpha2 = phase2();
                                
                                            if (&alpha2 == &alpha) continue;
                                
                                            const scalarField& dgdt2 = alpha2.dgdt();
                                
                                            forAll(dgdt2, celli)
                                            {
                                                if (dgdt2[celli] > 0.0 && alpha2[celli] < 1.0)
                                                {
                                                    Sp[celli] -= dgdt2[celli]*(1.0 - alpha2[celli]);
                                                    Su[celli] += dgdt2[celli]*alpha[celli];
                                                }
                                                else if (dgdt2[celli] < 0.0 && alpha2[celli] > 0.0)
                                                {
                                                    Sp[celli] += dgdt2[celli]*alpha2[celli];
                                                }
                                            }
                                        }
                                
                                        MULES::explicitSolve
                                        (
                                            geometricOneField(),
                                            alpha,
                                            phiAlpha,
                                            Sp,
                                            Su
                                        );
                                1 Reply Last reply Reply Quote
                                • M
                                  mohui last edited by

                                  后来我又做了测试,把limit以及limitsum注释掉了,重新编译运算,发现alpha的总和仍然是守恒的,按道理我注释掉了这两个,说白了也是只有explicitsolve求解,运算的时候结果很好的保证了有界性,还请各位大神能否解释这是为什么。

                                  1 Reply Last reply Reply Quote
                                  • M
                                    mohui last edited by

                                    不好意思,昨天编译的仓促,实际上并没有真正的编译我所改掉的程序,所以运行的还是原来的程序。直到刚才我才正确重新编译,运算发现马上浮点溢出。说明limit和limtsum确实有作用。

                                    1 Reply Last reply Reply Quote
                                    • M
                                      mohui @队长别开枪 last edited by 李东岳

                                      @队长别开枪 ,最近做了一些测试,发现多相流求解器为什么求解每一相的同时保证alpha的总和有界。这里最关键的地方代码是limitsum通量限制器。这个是limitsum的源码,

                                      void Foam::MULES::limitSum(UPtrList<scalarField>& phiPsiCorrs)
                                      {
                                          forAll(phiPsiCorrs[0], facei)
                                          {
                                              scalar sumPos = 0;
                                              scalar sumNeg = 0;
                                      
                                              for (int phasei=0; phasei<phiPsiCorrs.size(); phasei++)
                                              {
                                                  if (phiPsiCorrs[phasei][facei] > 0)
                                                  {
                                                      sumPos += phiPsiCorrs[phasei][facei];
                                                  }
                                                  else
                                                  {
                                                      sumNeg += phiPsiCorrs[phasei][facei];
                                                  }
                                              }
                                      
                                              scalar sum = sumPos + sumNeg;
                                      
                                              if (sum > 0 && sumPos > VSMALL)
                                              {
                                                  scalar lambda = -sumNeg/sumPos;
                                      
                                                  for (int phasei=0; phasei<phiPsiCorrs.size(); phasei++)
                                                  {
                                                      if (phiPsiCorrs[phasei][facei] > 0)
                                                      {
                                                          phiPsiCorrs[phasei][facei] *= lambda;
                                                      }
                                                  }
                                              }
                                              else if (sum < 0 && sumNeg < -VSMALL)
                                              {
                                                  scalar lambda = -sumPos/sumNeg;
                                      
                                                  for (int phasei=0; phasei<phiPsiCorrs.size(); phasei++)
                                                  {
                                                      if (phiPsiCorrs[phasei][facei] < 0)
                                                      {
                                                          phiPsiCorrs[phasei][facei] *= lambda;
                                                      }
                                                  }
                                              }
                                          }
                                      }
                                      

                                      看完代码后,可以比较清晰的看出,它其实做的工作量就是统计每一个面上的每一相的alpha的通量,然后进行了修正,得到了修正系数再乘上了alpha的通量,输出。

                                      surfaceScalarField& phiAlpha = phiAlphaCorrs[phasei];
                                              phiAlpha += upwind<scalar>(mesh_, phi_).flux(alpha)
                                      

                                      这两行代码是指,最后这个修正的通量加上原来的通量,作为alpha的通量计算。

                                      队长别开枪 李东岳 D 3 Replies Last reply Reply Quote
                                      • 队长别开枪
                                        队长别开枪 教授 @mohui last edited by

                                        @mohui OpenFOAM的另一个非官方版本(foam-extend)里alpha方程的源代码比官方的看起来简洁一点,比如求解器interFoam,extend版本(4.0)中的alphaEqn.H为

                                        {
                                            word alphaScheme("div(phi,alpha)");
                                            word alpharScheme("div(phirb,alpha)");
                                        
                                            surfaceScalarField phic = mag(phi/mesh.magSf());
                                            phic = min(interface.cAlpha()*phic, max(phic));
                                            surfaceScalarField phir = phic*interface.nHatf();
                                        
                                            for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
                                            {
                                                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;
                                            }
                                        
                                            Info<< "Liquid phase volume fraction = "
                                                << alpha1.weightedAverage(mesh.V()).value()
                                                << "  Min(alpha1) = " << min(alpha1).value()
                                                << "  Max(alpha1) = " << max(alpha1).value()
                                                << endl;
                                        }
                                        

                                        官方版本(4.x)中的为

                                        {
                                            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)
                                                {
                                                    FatalErrorInFunction
                                                        << "Sub-cycling is not supported "
                                                           "with the CrankNicolson ddt scheme"
                                                        << exit(FatalError);
                                                }
                                        
                                                ocCoeff =
                                                    refCast<const fv::CrankNicolsonDdtScheme<scalar>>(ddtAlpha())
                                                   .ocCoeff();
                                            }
                                            else
                                            {
                                                FatalErrorInFunction
                                                    << "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));
                                            }
                                        
                                            surfaceScalarField::Boundary& phicBf =
                                                phic.boundaryFieldRef();
                                        
                                            // 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;
                                                }
                                            }
                                        
                                            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();
                                        
                                                if (alphaApplyPrevCorr && talphaPhiCorr0.valid())
                                                {
                                                    Info<< "Applying the previous iteration compression flux" << endl;
                                                    MULES::correct(alpha1, alphaPhi, talphaPhiCorr0.ref(), 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.ref(), 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;
                                        }
                                        

                                        具体区别我不大清楚,但看起来还是有区别的。

                                        M 1 Reply Last reply Reply Quote
                                        • 李东岳
                                          李东岳 管理员 @mohui last edited by

                                          @mohui
                                          有点迷惑。
                                          如果你要验证求解一个相方程和两个相方程的区别,你就不需要研究MULES::limiter(),
                                          如果你要验证MULES的限制器,你可以直接采用求解单alpha方程来仔细研究通量限制器。
                                          哪个是研究的重点?:confused:

                                          CFD高性能服务器 http://dyfluid.com/servers.html

                                          1 Reply Last reply Reply Quote
                                          • 李东岳
                                            李东岳 管理员 @mohui last edited by

                                            @mohui 在 of中的多相流求解器 中说:

                                            我说的错误是这个,里面是没有alpha1,alpha2的

                                            是说 compreisblInterFoam解析 中下面这个方程有误?

                                            \begin{equation}
                                            \mathrm{dgdt}= \left( \frac{\alpha_2}{\rho_2}\frac{\mathrm{D}\rho_2}{\mathrm{D}t} - \frac{\alpha_1}{\rho_1}\frac{\mathrm{D}\rho_1}{\mathrm{D}t} \right)
                                            \end{equation}

                                            CFD高性能服务器 http://dyfluid.com/servers.html

                                            M 1 Reply Last reply Reply Quote
                                            • M
                                              mohui @李东岳 last edited by

                                              @李东岳对,这个我推导里面是没有alpha的,源程序里也是没有的。

                                              李东岳 1 Reply Last reply Reply Quote
                                              • M
                                                mohui @队长别开枪 last edited by

                                                @队长别开枪 这两个区别在于后者采用了CrankNicolson格式。现在我最新给的代码是解释多相流求解器如何保证alpha之和有界的,你可以看下

                                                1 Reply Last reply Reply Quote
                                                • 李东岳
                                                  李东岳 管理员 @mohui last edited by

                                                  @mohui 非常感谢!原文已更新:cheeky:

                                                  有关MULES的验证(暂时说不上研究,因为基金会已经植入了),我这里分为俩部分验证,第一部分验证已经完毕并且投稿,第二部分验证还在进行中,不过进度很慢毕竟还有别的课题在进行,我只能建议你参考Flux Correct Transport部分数值算法内容,MULES的思想来源于此,只不过基金会将此方法用于多相流

                                                  有关求解一个方程或者两个方程的验证,应该不是很难,并且结果应该一致。在你之前发的相图上,我没有看出太大差别,是否有曲线图。另外,是否对比过1D Riemann问题的解析解?

                                                  CFD高性能服务器 http://dyfluid.com/servers.html

                                                  M M yhdthu 3 Replies Last reply Reply Quote
                                                  • M
                                                    mohui @李东岳 last edited by

                                                    @李东岳 东岳老师,你可能还对我提的问题存在误区,我并不是对mules做验证,仅是分析OF中多相流求解器和两相流求解器对与两相流求解的区别。以我目前的测试检验,认为多相流求解器每一相求解并能够同时保证总的alpha之和等于的1的关键在于limitsum这个限制器,关于这个限制器的代码我也上传了,通过代码可以分析的出来,就是在mules求解之前对alphaphi进行了修正,以保证流入流出单元的alphaphi总和为0。其次,就如我前面所言,如果单纯的mules求解每一相是无法保证alpha1+alpha2=1的,虽然我给的图上没有很大区别,但是我输出alpha1+alpha2最大值为2,这是不对的。以上是我这几天的发现。

                                                    1 Reply Last reply Reply Quote
                                                    • M
                                                      maoyanjun_dut @李东岳 last edited by

                                                      @李东岳

                                                      @李东岳 在 of中的多相流求解器 中说:

                                                      @mohui 非常感谢!原文已更新:cheeky:

                                                      有关MULES的验证(暂时说不上研究,因为基金会已经植入了),我这里分为俩部分验证,第一部分验证已经完毕并且投稿,第二部分验证还在进行中,不过进度很慢毕竟还有别的课题在进行,我只能建议你参考Flux Correct Transport部分数值算法内容,MULES的思想来源于此,只不过基金会将此方法用于多相流

                                                      有关求解一个方程或者两个方程的验证,应该不是很难,并且结果应该一致。在你之前发的相图上,我没有看出太大差别,是否有曲线图。另外,是否对比过1D Riemann问题的解析解?

                                                      岳哥,公式20更新了,公式19却没更新。建议更新一下,否则看着怪怪的。

                                                      Blog:http://maoyanjun.top
                                                      厚德如海,纳川有藏。上善若水,利物不张。大哉斯言,勤勉勿忘。

                                                      赵 1 Reply Last reply Reply Quote
                                                      • 赵
                                                        赵一铭 @maoyanjun_dut last edited by

                                                        @maoyanjun_dut

                                                        确实!好眼力。
                                                        @李东岳

                                                        1 Reply Last reply Reply Quote
                                                        • yhdthu
                                                          yhdthu 讲师 @李东岳 last edited by

                                                          @李东岳 前辈,compressibleInterFoam解析的(19)式忘记改了,括号内分子没有alpha:laughing:

                                                          长风破浪会有时,直挂云帆济沧海

                                                          1 Reply Last reply Reply Quote
                                                          • 李东岳
                                                            李东岳 管理员 last edited by

                                                            谢谢,刚改了,interPhaseChangeFoam一直没时间写,估计得10月份开会回来之后能写...
                                                            :big_mouth:

                                                            CFD高性能服务器 http://dyfluid.com/servers.html

                                                            yhdthu 1 Reply Last reply Reply Quote
                                                            • yhdthu
                                                              yhdthu 讲师 @李东岳 last edited by

                                                              @李东岳 好的👌,到时候跟前辈讨论学习一下~

                                                              长风破浪会有时,直挂云帆济沧海

                                                              1 Reply Last reply Reply Quote
                                                              • J
                                                                Jacobian last edited by

                                                                @mohui 没太看懂limitSum为何能保证体积分数之和为1,可否讲解一下?

                                                                1 Reply Last reply Reply Quote
                                                                • D
                                                                  dzw05 教授 @mohui last edited by

                                                                  @mohui 其实也没太懂为什么limitsum可以保证alpha总和有界。limitsum限制的是所谓的修正通量,即“高阶通量”减去“迎风通量”的部分,至于为什么这部分通量总和等于零就可以保证alpha总和有界,这个还是没看懂。

                                                                  自主匠心,普惠仿真。

                                                                  1 Reply Last reply Reply Quote
                                                                  • 小
                                                                    小龙 last edited by

                                                                    @李东岳 @mohui @赵一铭 我最近再做空化模拟,OF中只有多相流VOF模型,VOF做空化不太适合,想在OF里编译一下mixture模型,前辈们觉得实现的可行性高吗。实现的难度大吗?我初步的想法是基于interPhaseChangeFoam来修改。前辈们有什么建议?

                                                                    1 Reply Last reply Reply Quote
                                                                    • 李东岳
                                                                      李东岳 管理员 last edited by 李东岳

                                                                      如果单纯的mules求解每一相是无法保证alpha1+alpha2=1的

                                                                      现在我就想确定,alpha2通过对应方程推导,用MULES求解不可行的吗?

                                                                      MULES的出发点是保证变量的有界。所以如果你用MULES求解alpha1,那么理论上alpha1是有界的。alpha2的求解没有调用MULES,他是通过1-alpha1算出来的,如果alpha1有界,alpha2也有界。

                                                                      如果你用算法同时求解alpha1和alpha2,如何处理耦合?这俩个变量是耦合在一起的,就像速度和压力。你不能单独的去分离求解,如果分离求解就就需要迭代。迭代就导致计算速度变慢。因此现存大厂据我所知都是只求解alpha1,然后alpha2=1-alpha1。

                                                                      我在去你年底验证了MULES,理论上是可以保证有界,但是真实计算的时候,还是会越界。尤其是物理模型比较复杂的时候。个人觉得这方面内容搞出来,绝对是个好文。

                                                                      CFD高性能服务器 http://dyfluid.com/servers.html

                                                                      1 Reply Last reply Reply Quote
                                                                      • First post
                                                                        Last post

                                                                      CFD中文网 | 东岳流体 | 京ICP备15017992号-2
                                                                      论坛登录问题反馈可联系 li.dy@dyfluid.com