[TommyCP]

Dear fellow FOAMers, I’ve run into a problem with an interFoam simulation that I would like to request your help with.

I’m trying to simulate the filling of a tank with a high-viscosity liquid. The vessel has both the inlet and outlet placed on top. At the start, the vessel is filled half with air, half with liquid. The vessel is then filled with liquid until full.

The expectation of the simulation is that the vessel would be filled up almost completely with some air pockets remaining. However, it currently will not fill for more than 80-85 %. Once it is filled up this far, the liquid starts getting ‘sucked’ into the outlet in an unphysical way. Oddilly enough, it does not come into contact with the ceiling of the vessel nor with the walls of the outlet.

I’ve tried many things to create a more realistic behavior, all to no avail. Things I’ve tried:

Geometry-related:
- tilting the whole mesh at a very slight angle
- increasing the height of the outlet (see attached images)

Mesh-related:
- adding more/less layers on the walls
- increasing/decreasing mesh fineness
- using a more fine mesh only around the inlet and outlet
- note: checkMesh never reports any errors

transportProperties-related:
- varying sigma from 0.07 down to 0

BC-related:
- trying various boundary conditions for inlet, outlet, and walls (see below for current BCs), based on the capillaryRise and damBreak examples as well as various hints on these forums. To name a few of the important ones:
o for alpha: * inlet: inletOutlet; * walls: zeroGradient or constrantAlphaContactAngle (with various values for theta0); * outlet: inletOutlet or zeroGradient
o for p_rgh: * inlet: fixedFluxPressure; * walls: fixedFluxPressure or fixedValue; * outlet: fixedValue or totalPressure
o for U: * inlet: flowRateInletVelocity; * walls: fixedValue; * outlet: pressureInletOutletVelocity or zeroGradient.

Solver-related:
- trying interpolationSchemes ‘linear’ and ‘midPoint’
- trying divSchemes ‘Gauss limitedLinearV 1’ and ‘Gauss upwind’ for div(rhoPhi,U)
- trying divSchemes ‘Gauss linear’ and ‘Gauss interfaceCompression’ for div(phirb,alpha)
- playing with nAlphacorr, nAlphaSubCycles, cAlpha
- enabling/disabling the momentumPredictor
- reducing maxCo from 0.95 to 0.2

Currently, I’m at a loss so I would like to ask if anyone has an idea how to get more realistic results. Below are some of the settings I last used. I also attached a few screenshots of the simulation. Thank you for any help you can provide!


alpha

internalField   uniform 0; 
   
  boundaryField
  {
      inlet
      {
          type            inletOutlet;
          inletValue      uniform 1;
          value           uniform 1;
      }
      main
      {
          type            constantAlphaContactAngle;
          theta0          45;
          limit           gradient;
          value           nonuniform 0;
      }
      outlet
      {
          type            inletOutlet;
          inletValue      uniform 0;
          value           uniform 0;
      }
  }

p_rgh

internalField   uniform 0;
   
  boundaryField
  {
      inlet
      {
          type            fixedFluxPressure;
          value           uniform 0;
      }
      main
      {
          type            fixedFluxPressure;
          value           uniform 0;
      }
      outlet
      {
          type            totalPressure;
              p0          uniform 0;
              U           U;
              phi         phi;
              rho         rho;
              psi         none;
              gamma       1;
              value       uniform 0;
      }
  }

U

internalField   uniform (0 0 0);
   
  boundaryField
  {
      inlet
      {
          type            flowRateInletVelocity;
          massFlowRate    constant 0.1;
          value           uniform (0 0 -0.1);
      }
      main
      {
          type            fixedValue;
          value           uniform (0 0 0);
      }
      outlet
      {
          type            zeroGradient;
      }
  }

fvSolution

solvers
  {
      alpha.water
      {
          nAlphaCorr      1;
          nAlphaSubCycles 2;
          cAlpha          1;
      }
   
      pcorr
      {
          solver          PCG;
          preconditioner  DIC;
          tolerance       1e-10;
          relTol          0.01;
      }
   
      p_rgh
      {
          solver          GAMG;
          smoother        GaussSeidel
          tolerance       1e-07;
          relTol          0.05;
      }
   
      p_rghFinal
      {
          $p_rgh;
          tolerance       1e-08;
          relTol          0.01;
      }
   
      U
      {
          solver          smoothSolver;
          smoother        symGaussSeidel;
          tolerance       1e-06;
          relTol          0;
      }
   
      UFinal
      {
          $U;
          tolerance       1e-08;
          relTol          0;
      }
  }
   
  PIMPLE
  {
      momentumPredictor no;
      nCorrectors     3;
      nOuterCorrectors   1;
      nNonOrthogonalCorrectors 0;
  }

fvSchemes

ddtSchemes
  {
      default            Euler;
  }
   
  gradSchemes
  {
      default         Gauss linear;
  }
   
  divSchemes
  {
      div(rhoPhi,U) Gauss limitedLinearV 1;
      div(phi,alpha) Gauss vanLeer;
      div(phirb,alpha) Gauss interfaceCompression;
      div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
  }
   
  laplacianSchemes
  {
      default         Gauss linear corrected;
  }
   
  interpolationSchemes
  {
      default         midPoint;
  }
   
  snGradSchemes
  {
      default         corrected;
  }
   
  fluxRequired
  {
      default         no;
      p_rgh;
      pcorr;
      alpha;
  }

controlDict

application     interFoam;
  startFrom       startTime;
  startTime       0;
  stopAt          endTime;
  endTime         125;
  deltaT          2e-3;
  adjustTimeStep  yes;
  writeControl    adjustableRunTime;
  writeInterval   5;
  purgeWrite      0;
  writeFormat     ascii;
  writePrecision  6;
  writeCompression uncompressed;
  timeFormat      general;
  timePrecision   6;
  runTimeModifiable yes;
  adjustTimeStep  yes;
  maxCo           0.95;
  maxAlphaCo      0.95;
  maxDeltaT       1;

[TommyCP]

I've done some additional experimentation with my case. Using water (instead of the viscous liquid used earlier) did not lead to better results -- again, the liquid phase got 'sucked' into the outlet in an unphysical way. I'm currently running out of ideas on how to get a more realistic behavior.

If anyone has any ideas on this, I'd really appreciate it!

[TobM]

Try to use pressureInletOutletValocity BC on outlet (refer to point 5.2.3.2 in the user guide).
Use "normal" schemes as in damBreak or so.
Your mesh seems pretty coarse at the outlet, only 10 or so cells is not enough.

[TommyCP]

Thank you for your suggestion. I've increased mesh fineness 2x in general and 4x around the inlet and outlet. Also, I changed my schemes back to those used in the damBreak example. Unfortunately, it did not change the behavior of my system. As for pressureInletOutletVelocity, I was already using this on the outlet.

My current BCs are as follows:

alpha
- inlet: inletOutlet
- outlet: inletOutlet
- wall: constantAlphaContactAngle

p_rgh
- inlet: fixedFluxPressure
- outlet: totalPressure
- wall: fixedFluxPressure

U
- inlet: flowRateInletVelocity
- outlet: pressureInletOutletVelocity
- wall: fixedValue

If you have any further suggestions, I'd much appreciate it.

Edit: just upgraded to OpenFOAM 4.0. Will try the new noSlip BC in U and make some changes to fvSolution.

[TobM]

You haven't read the user guide .
pressureInletOutletVelocity for U is used together with totalPressure for p_rgh

[TommyCP]

Sorry, that is an error in my post. I do indeed use totalPressure for p_rgh and pressureInletOutletVelocity for U.

[TobM]

Please compare your entire settings with settings from a similar tutorial, and check if your mesh is ok --> checkMesh.

[TommyCP]

Meshes always check out OK with checkMesh. I'm currently testing my case with settings taken from the angledDuct tutorial. I will let you know how it works out.

Thank you for your suggestions thus far.

[TommyCP]

Update: settings from the angledDuct tutorial result in a more physically correct behavior. Unfortunately, the case blows up when the tank is around 85 % full. I've found no way to avoid this.

I think it's time to give up and move. I've been trying to a lot of things to make this works, but apparently there's more to this seemingly simple case than meets the eye.

Again, thank you for your help.

[mwmalkawi]

Hi Tommy

Have yo umnaged to find the soluiton fo your problem

Thakns