twoPhaseEulerFoam convergence problem
 

Hello to everyone,
I would gladly accept any help or hint with my diverging solution using twoPhaseEulerFoam.

I am trying to simulate the separated inflow of air on top of water.
Internal and boundary are modeled using codedValue:
alpha =0; if y <= h/2 - deltaH/2,
alpha =0; if y >= h/2 + deltaH/2,
alpha =0; if h/2 - deltaH/2 <=y <= h/2 + deltaH/2,
Where h is height of domain, deltaH is the height of smooth step function.
k = 1 / ( f( h/2.0 + deltaH/2.0 ) - f( h/2.0 - deltaH/2.0 ) );
b = 1 - k * f( h/2.0 + deltaH/2.0 );
f( y ) = exp( y ) / ( exp( y )+1 );

Initial fields, boundary conditions, and alpha results for time 0.0052 sec are attached.
Another question: why we cannot set alpha = 0 in BC?
Thank You,

Best regards,
Tonnykz


Initial fields are:
alpha:
dimensions [0 0 0 0 0 0 0];
internalField #codeStream
{
codeInclude
#{
#include "fvCFD.H"
#};

codeOptions
#{
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude
#};

//libs needed to visualize BC in paraview
codeLibs
#{
-lmeshTools \
-lfiniteVolume
#};

code
#{
const IOdictionary& d = static_cast<const IOdictionary&>(dict);
const fvMesh& mesh = refCast<const fvMesh>(d.db());
scalarField alpha(mesh.nCells(), 0.);

const scalar k = 2000.01; //4000.03;
const scalar b = -1002.01; //-2004.51;
const scalar h = 0.01;
const scalar deltaH = 0.002; //0.001;

forAll( alpha, i )
{
const scalar x = mesh.C()[i][0];
const scalar y = mesh.C()[i][1];

if ( y <= h/2.0 - deltaH/2.0 )
{
alpha[i] = 0;
}
else if ( y >= h/2.0 + deltaH/2.0 )
{
alpha[i] = 1;
}
else
{
alpha[i] = k * ( exp( y ) / ( exp( y ) + 1 ) ) + b;
}
};

alpha.writeEntry("", os);
#};
};
boundaryField
{
inlet
{
type codedFixedValue;
value uniform 0;
name inletProfile1;

codeInclude #{
#include "fvCFD.H"
#include <cmath>
#include <iostream>
#};

code #{
const fvPatch& boundaryPatch = patch();
const vectorField& Cf = boundaryPatch.Cf();
scalarField& field = *this;

field = patchInternalField(); //take value from initialization

const scalar k = 2000.01; //4000.03;
const scalar b = -1002.01; //-2004.51;
const scalar h = 0.01;
const scalar deltaH = 0.002; //0.001;

forAll(Cf, faceI)
{
if ( Cf[faceI].y() <= h/2.0 - deltaH/2.0 )
{

field[faceI] = 0.01;
}
else if (
Cf[faceI].y() >= h/2.0 + deltaH/2.0
)
{
field[faceI] = 0.99;
}
else
{
field[faceI] = k * ( exp( Cf[faceI].y() )/( exp( Cf[faceI].y() ) + 1 ) ) + b;
}
}
#};

codeOptions #{
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude
#};
}

outlet
{
type zeroGradient;
}

walls
{
type zeroGradient;
}

frontAndBackPlanes
{
type empty;
}
}

p:

dimensions [1 -1 -2 0 0 0 0];
internalField uniform 1e5;
boundaryField
{
inlet
{
type calculated;
value $internalField;
}

outlet
{
type calculated;
value $internalField;
}

walls
{
type calculated;
value $internalField;
}

frontAndBackPlanes
{
type empty;
}
}


p_rgh:
dimensions [1 -1 -2 0 0 0 0];
internalField uniform 1e5;
boundaryField
{
inlet
{
type fixedFluxPressure;
value $internalField;
}

outlet
{
type prghPressure;
p $internalField;
value $internalField;
}

walls
{
type fixedFluxPressure;
value $internalField;
}

frontAndBackPlanes
{
type empty;
}
}


T.air:
dimensions [0 0 0 1 0 0 0];
internalField uniform 300;
boundaryField
{
inlet
{
type fixedValue;
value uniform 300;
}

outlet
{
type inletOutlet;
phi phi.air;
inletValue uniform 300;
value uniform 300;
}

walls
{
type zeroGradient;
}

frontAndBackPlanes
{
type empty;
}
}


T.water:
dimensions [0 0 0 1 0 0 0];
internalField uniform 300;
boundaryField
{
inlet
{
type zeroGradient;
}

outlet
{
type inletOutlet;
phi phi.water;
inletValue uniform 300;
value uniform 300;
}

walls
{
type zeroGradient;
}

frontAndBackPlanes
{
type empty;
}
}


U.air:
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0.05 0 0);
boundaryField
{
inlet
{
type interstitialInletVelocity;
inletVelocity uniform (0.05 0 0);
alpha alpha.air;
value $internalField;
}

outlet
{
type pressureInletOutletVelocity;
phi phi.air;
value $internalField;
}

walls
{
type noSlip;
}

frontAndBackPlanes
{
type empty;
}
}


U.water:
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0.05 0 0);
boundaryField
{
inlet
{
type interstitialInletVelocity;
inletVelocity uniform (0.05 0 0);
alpha alpha.water;
value $internalField;
/* type fixedValue;
value uniform (0 0 0);*/
}

outlet
{
type pressureInletOutletVelocity;
phi phi.water;
value $internalField; /*fixedValue;
value uniform (0 0 0);*/
}

walls
{
type noSlip;/*fixedValue;
value uniform (0 0 0);*/
}

frontAndBackPlanes
{
type empty;
}
}


fvSchemes:


ddtSchemes
{
default Euler;
}

gradSchemes
{
default Gauss linear;
}

divSchemes
{
default none;

"div\(phi,alpha.*\)" Gauss vanLeer;
"div\(phir,alpha.*\)" Gauss vanLeer;

"div\(alphaRhoPhi.*,U.*\)" Gauss limitedLinearV 1;
"div\(phi.*,U.*\)" Gauss limitedLinearV 1;

"div\(alphaRhoPhi.*,(h|e).*\)" Gauss limitedLinear 1;
"div\(alphaRhoPhi.*,K.*\)" Gauss limitedLinear 1;
"div\(alphaPhi.*,p\)" Gauss limitedLinear 1;

"div\(\(\(\(alpha.*\*thermo:rho.*\)\*nuEff.*\)\*de v2\(T\(grad\(U.*\)\)\)\)\)" Gauss linear;
}

laplacianSchemes
{
default Gauss linear uncorrected;
// bounded Gauss linear uncorrected;
}

interpolationSchemes
{
default linear;
}

snGradSchemes
{
default uncorrected;
// bounded uncorrected;
}

fvSolution:

solvers
{
"alpha.*"
{
nAlphaCorr 1;
nAlphaSubCycles 2;
/*
smoothLimiter 0.1;

implicitPhasePressure yes;
solver smoothSolver;
smoother symGaussSeidel;
tolerance 1e-9;
relTol 0;
minIter 1;*/
}

p_rgh
{
solver GAMG;
smoother DIC;
tolerance 1e-8;
relTol 0;
}

p_rghFinal
{
$p_rgh;
relTol 0;
}

"U.*"
{
solver smoothSolver;
smoother symGaussSeidel;
tolerance 1e-5;
relTol 0;
minIter 1;
}

"e.*" //"(h|e).*"
{
solver smoothSolver;
smoother symGaussSeidel;
tolerance 1e-8; //1e-6;
relTol 0;
minIter 1;
}
/*
"Theta.*"
{
solver smoothSolver;
smoother symGaussSeidel;
tolerance 1e-6;
relTol 0;
minIter 1;
}

"(k|epsilon).*"
{
solver smoothSolver;
smoother symGaussSeidel;
tolerance 1e-5;
relTol 0;
minIter 1;
}*/
}

PIMPLE
{
nOuterCorrectors 3;
nCorrectors 1;
nNonOrthogonalCorrectors 0;
}

relaxationFactors
{
equations
{
".*" 1;
}
}

controlDict:


application reactingTwoPhaseEulerFoam;

startFrom startTime;

startTime 0;

stopAt endTime;

endTime 1e-2;

deltaT 1e-4;

writeControl runTime; //adjustableRunTime;

writeInterval 1e-4;

purgeWrite 0;

writeFormat ascii;

writePrecision 6;

writeCompression off;

timeFormat general;

timePrecision 6;

runTimeModifiable on;

adjustTimeStep no; //true;

maxCo 0.5;

maxDeltaT 1;//1e-05;

allowSystemOperations 1;

Hello Tonnykz,



I use multiphaseEulerFoam and I figured out that my simulations only run when I use PISO (fvSolutions - no nOuterCorrectors). I allow a Co-Number of 0.99 (ControlDict).

I am sorry that I do not know why you get an error regarding alpha but maybe the hint with Piso helps you.



Kind regards,
Stefanie

Hello Stefanie,
Thank You for Your suggestion.
Actually, I have made progress with this issue. There were couple of things to fix:
1. For 0/U.* Instead of putting interstitialInletVelocity one should use fixedValue with specified velocity value. Because the former one should be used for phase that is dispersed in continuous phase, e.g. bubbles in column.
2. Max Courant number is 0.3. Although, I still vaguely understand the role of this number in implicit solver.
3. AdjustableTimeStep yes.
4. I have disabled virtual mass force in my phase properties as well as aspect ratio. Because, I was testing this BC and solver for enclosed cavity case with all walls, whereas lower half of it is filled with water and top is air. The test case diverged. And after taking a look into governing equations from http://powerlab.fsb.hr/ped/kturbo/Op...chePhD2002.pdf, one can see that even if we put velocity field equal to zero and expect current state to be steady we see movement. The reason is virtual mass force in formulation. After disabling it physically correct solution was obtained. Picture enclosed.
Sorry for poor formulation and lack of formulas.

I hope that will help someone with similar problems.
Thank You,
Tonnykz

Hello Tonnykz,



thank you for sharing your findings and especially your case. I recognized how you define p and p-rgh different from me. This helped me to fix my BCs. My simulation runs since yesterday evening without obvious errors now and it seems like I have normal air pressure and a freestream outlet with my fixedVelocity/Groovy inlets.

( https://www.cfd-online.com/Forums/op...eulerfoam.html)

Thank you for your help!


Kind regards, Stefanie

Hello Stefanie,
I am happy to help.
Regarding my current simulation:
1. Velocity of air in water is much higher than imposed (by the power of 10!). Same for water in air. What seems not physical to me since fractions of one phase should have the same speed as of continuous phase. Probably that was the reason of imposing interstitial inlet velocity. I will read a thesis of H.Rusche regarding this solver. ( http://powerlab.fsb.hr/ped/kturbo/Op...chePhD2002.pdf ).
2. Although I impose uniform temperature everywhere, its field have weird distribution. Upper wall seems to be cooling phases. Although zero gradient condition was imposed.
Will work on it.
Best regards,
Tonnykz

Update on progress
 

Hello everyone,
Purpose is to provide update on current simulation, maybe someone going the same path will find it useful.
1. High velocity in dispersed phase is fictional since the phase concentration is small (~e-10).
2. I have moved from reactingTwoPhaseEulerFoam to interFoam. Since I could not obtain convergence in my simulations with sharp interface. Simulation with same BC and initial conditions works fine with interFoam.
3. I have added energy equation to interFoam and plan to add species transport as well as chemical reactions in it.
Best regards,
Tonnykz