CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Running, Solving & CFD

twoPhaseEulerFoam convergence problem

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree2Likes
  • 1 Post By Stefanie.S.W.
  • 1 Post By tonnykz

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   April 4, 2019, 15:34
Default twoPhaseEulerFoam convergence problem
  #1
New Member
 
Join Date: Oct 2016
Posts: 29
Rep Power: 6
tonnykz is on a distinguished road
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;
Attached Images
File Type: jpg alpha.jpg (29.5 KB, 24 views)
File Type: jpg mesh.jpg (101.5 KB, 19 views)
tonnykz is offline   Reply With Quote

Old   April 9, 2019, 12:09
Default
  #2
Member
 
Stefanie Wolf
Join Date: Nov 2018
Location: Aachen
Posts: 32
Rep Power: 4
Stefanie.S.W. is on a distinguished road
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
tonnykz likes this.
Stefanie.S.W. is offline   Reply With Quote

Old   April 9, 2019, 12:24
Default
  #3
New Member
 
Join Date: Oct 2016
Posts: 29
Rep Power: 6
tonnykz is on a distinguished road
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
Attached Images
File Type: jpg alpha.jpg (27.5 KB, 29 views)
Stefanie.S.W. likes this.
tonnykz is offline   Reply With Quote

Old   April 10, 2019, 07:50
Default
  #4
Member
 
Stefanie Wolf
Join Date: Nov 2018
Location: Aachen
Posts: 32
Rep Power: 4
Stefanie.S.W. is on a distinguished road
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.

( non-physical Results with multiphaseEulerFoam)

Thank you for your help!


Kind regards, Stefanie
Stefanie.S.W. is offline   Reply With Quote

Old   April 13, 2019, 10:02
Default
  #5
New Member
 
Join Date: Oct 2016
Posts: 29
Rep Power: 6
tonnykz is on a distinguished road
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
tonnykz is offline   Reply With Quote

Old   April 22, 2019, 15:43
Default Update on progress
  #6
New Member
 
Join Date: Oct 2016
Posts: 29
Rep Power: 6
tonnykz is on a distinguished road
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
tonnykz is offline   Reply With Quote

Reply

Tags
multiphaseeulerfoam, twophaseeulerfoam

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Convergence problem in Fluent for quenching process kaeran FLUENT 4 December 1, 2014 02:14
Rotate frame reference convergence problem! wjy-c CFX 2 September 26, 2014 06:03
Centrifugal pump OpenFOAM, convergence problem, ANSA model RDD OpenFOAM Running, Solving & CFD 0 July 5, 2014 09:12
Convergence Problem in Axisymmetric Periodic Flow atheresia FLUENT 3 February 10, 2014 03:00
Convergence of CFX field in FSI analysis nasdak CFX 2 June 29, 2009 01:17


All times are GMT -4. The time now is 18:33.