Navigation

    CFD中文网

    CFD中文网

    • Login
    • Search
    • 最新

    Coupled level set-VOF方法

    OpenFOAM
    2
    3
    364
    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.
    • Z
      Zhujh last edited by 李东岳

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

      ae8a8d4d-90d6-4702-a0b3-2fded1619f11-image.png

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

      72aa53d8-7890-48af-9176-ddbec5952d3a-image.png
      3e7dcaaf-8c2b-4df9-8d63-f18acf57a2c7-image.png
      4617006d-99d1-4499-9603-3b426f84ef57-image.png
      但是当我换到溃坝的算例的时候,出现了奇怪的问题,后处理结果来看,水体好像变得不连续了,最后会扩散到整个计算域。看上去像是没有了重力,但是算例里面我是加了的,想请教一下各位大佬,看看问题出在哪儿。
      1635111b-c932-47b4-a78b-16eb6c7cecd5-image.png 6fd36484-cc7d-4e10-b0ac-e8f6052c7e45-image.png
      以下是在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;
      }
      李东岳 1 Reply Last reply Reply Quote
      • 李东岳
        李东岳 管理员 @Zhujh last edited by

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

        CFD课程 改成线上了 http://dyfluid.com/class.html
        CFD高性能服务器 http://dyfluid.com/servers.html

        Z 1 Reply Last reply Reply Quote
        • Z
          Zhujh @李东岳 last edited by

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

          1 Reply Last reply Reply Quote
          • First post
            Last post

          CFD中文网 | 东岳流体 | 京ICP备15017992号-2