Hi Foamers,
I am trying to calculate transient flow field in a single bifurcating tube(using OF-2.1.1 and swak4Foam) with inlet boundary condition varying with space and time:
I will start illustrating my case by giving 'checkMesh' report:
Create time
Create polyMesh for time = 0
Time = 0
Mesh stats
points: 140311
faces: 1127971
internal faces: 1083988
cells: 521531
boundary patches: 4
point zones: 0
face zones: 1
cell zones: 1
Overall number of cells of each type:
hexahedra: 0
prisms: 125835
wedges: 0
pyramids: 0
tet wedges: 0
tetrahedra: 395696
polyhedra: 0
Checking topology...
Boundary definition OK.
Cell to face addressing OK.
Point usage OK.
Upper triangular ordering OK.
Face vertices OK.
Number of regions: 1 (OK).
Checking patch topology for multiply connected surfaces ...
Patch Faces Points Surface topology
INLET 910 612 ok (non-closed singly connected)
OUTLET1 562 390 ok (non-closed singly connected)
OUTLET2 566 394 ok (non-closed singly connected)
WALL 41945 21065 ok (non-closed singly connected)
Checking geometry...
Overall domain bounding box (-0.108847 -0.161672 -0.0190498) (0.108851 0.127 0.0190498)
Mesh (non-empty, non-wedge) directions (1 1 1)
Mesh (non-empty) directions (1 1 1)
Boundary openness (2.61084e-17 -1.48172e-16 3.15641e-17) OK.
Max cell openness = 2.91142e-16 OK.
Max aspect ratio = 7.00156 OK.
Minumum face area = 1.85341e-07. Maximum face area = 3.64501e-06. Face area magnitudes OK.
Min volume = 5.00262e-11. Max volume = 2.11586e-09. Total volume = 0.000347759. Cell volumes OK.
Mesh non-orthogonality Max: 63.6578 average: 14.781
Non-orthogonality check OK.
Face pyramids OK.
Max skewness = 2.94579 OK.
Coupled point location match (average 0) OK.
Mesh OK.
End
Boundary Conditions:
0/U
Code:
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0 0 0);
boundaryField
{
INLET
{
type groovyBC;
variables "vavg=(-0.000163-0.00336*cos(time()*0.3254)-0.3416*sin(time()*0.3254));r2=(pow(pos().x,2)+pow(pos().z,2));R2=sum(area())/pi;para=-((R2-r2)/R2)*normal();";
valueExpression "2*vavg*para";
value uniform (0 0 0);
}
OUTLET1
{
type zeroGradient;
}
OUTLET2
{
type zeroGradient;
}
WALL
{
type fixedValue;
value uniform (0 0 0);
}
}
0/p
Code:
dimensions [0 2 -2 0 0 0 0];
internalField uniform 0;
boundaryField
{
INLET
{
type zeroGradient;
}
OUTLET1
{
type fixedValue;
value uniform 0;
}
OUTLET2
{
type fixedValue;
value uniform 0;
}
WALL
{
type zeroGradient;
}
}
controlDict
Code:
libs ( "libOpenFOAM.so"
"libsimpleSwakFunctionObjects.so"
"libswakFunctionObjects.so"
"libgroovyBC.so" ) ;
application pimpleFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 20;
deltaT 0.02;
writeControl adjustableRunTime;
writeInterval 0.4;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression uncompressed;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
adjustTimeStep no;
maxCo 100;
functions {
initSwak {
// needed to allow DataEntry to work
type initSwakFunctionObject;
region region0;
}
};
functions
{
volFlowInlet
{
type swakExpression;
valueType patch;
patchName INLET;
verbose true;
expression "phi";
accumulations ( sum );
}
volFlowOutlet1
{
type swakExpression;
valueType patch;
patchName OUTLET1;
verbose true;
expression "phi";
accumulations ( sum );
}
volFlowOutlet2
{
type swakExpression;
valueType patch;
patchName OUTLET2;
verbose true;
expression "phi";
accumulations ( sum );
}
}
fvSchemes
Code:
ddtSchemes
{
default CrankNicholson 1;
}
gradSchemes
{
default Gauss linear;
grad(p) Gauss linear;
grad(p_rgh) Gauss linear;
grad(U) Gauss linear;
}
divSchemes
{
default none;
div(phi,U) Gauss linearUpwind grad(U);
div(phi,k) Gauss limitedLinear 1;
div(phi,epsilon) Gauss limitedLinear 1;
div(phi,R) Gauss limitedLinear 1;
div(R) Gauss linear;
div(phi,nuTilda) Gauss limitedLinear 1;
div((nuEff*dev(T(grad(U))))) Gauss linear;
}
laplacianSchemes
{
default none;
laplacian(nuEff,U) Gauss linear corrected;
laplacian((1|A(U)),p) Gauss linear corrected;
laplacian((1|A(U)),p_rgh) Gauss linear corrected;
laplacian(rAU,p_rgh) Gauss linear corrected;
laplacian(DkEff,k) Gauss linear corrected;
laplacian(DepsilonEff,epsilon) Gauss linear corrected;
laplacian(DREff,R) Gauss linear corrected;
laplacian(DnuTildaEff,nuTilda) Gauss linear corrected;
}
interpolationSchemes
{
default linear;
interpolate(U) linear;
}
snGradSchemes
{
default corrected;
}
fluxRequired
{
default no;
p ;
p_rgh ;
}
fvSolution
Code:
solvers
{
p
{
solver GAMG;
tolerance 1e-06;
relTol 0.01;
smoother GaussSeidel;
cacheAgglomeration true;
nCellsInCoarsestLevel 10;
agglomerator faceAreaPair;
mergeLevels 1;
}
pFinal
{
solver GAMG;
tolerance 1e-06;
relTol 0.01;
smoother GaussSeidel;
cacheAgglomeration true;
nCellsInCoarsestLevel 10;
agglomerator faceAreaPair;
mergeLevels 1;
}
"(U|kl|kt|omega)"
{
solver PBiCG;
preconditioner DILU;
tolerance 1e-06;
relTol 0;
}
"(U|kl|kt|omega)Final"
{
$U;
tolerance 1e-06;
relTol 0;
}
}
PIMPLE
{
nOuterCorrectors 50;
nCorrectors 1;
nNonOrthogonalCorrectors 2;
pRefCell 0;
pRefValue 0;
residualControl
{
"(p|U)"
{
tolerance 1e-4;
relTol 0;
}
}
}
relaxationFactors
{
fields
{
}
equations
{
"U.*" 0.5;
"kl.*" 1;
"kt.*" 1;
"omega.*" 1;
}
}
The simulation blows blows at time=0.9s or little higher. The pulse lasts 20s and the reynolds number peaks(Re=2200) at 5s and 15s. If I choose adjustable time step then the it blows by showing time step in range of 10^-13 or so. If I choose a fixed time step of 0.01 or 0.02 or 0.005 it blows showing high courant number or unusually high flow rates. I have employed meshes of various kinds.....even pure tetrahedral mesh doesn't seem to work.
Here is the link to download the case directory along with the mesh:
https://www.dropbox.com/sh/u7sqlfvqtc6uynq/1MMKNryk7Z
Following are the error message snippets:
Code:
Expression volFlowInlet : sum=0.000155198
Expression volFlowOutlet1 : sum=-1.01936e+12
Expression volFlowOutlet2 : sum=1.01938e+12
Courant Number mean: 4.53531e+15 max: 9.33679e+17
Time = 1.23
Also the time step continuity errors goes exponentially high:
Code:
PIMPLE: iteration 50
DILUPBiCG: Solving for Ux, Initial residual = 0.0124649, Final residual = 3.77817e-07, No Iterations 5
DILUPBiCG: Solving for Uy, Initial residual = 0.00973289, Final residual = 3.76167e-07, No Iterations 5
DILUPBiCG: Solving for Uz, Initial residual = 0.0198487, Final residual = 4.315e-08, No Iterations 6
GAMG: Solving for p, Initial residual = 0.0537055, Final residual = 0.00027507, No Iterations 3
GAMG: Solving for p, Initial residual = 0.00744695, Final residual = 6.84879e-05, No Iterations 3
GAMG: Solving for p, Initial residual = 0.00143221, Final residual = 8.22524e-06, No Iterations 3
time step continuity errors : sum local = 3.47333e+09, global = -7.64618e+08, cumulative = -2.18106e+10
PIMPLE: not converged within 50 iterations
ExecutionTime = 55306.7 s ClockTime = 88055 s
Velocity fields have exploded and there is a flow reversal too. I tried avoiding flow reversal using outletinlet boudary condition since my outlets should behave as inlets in the beginning:
Code:
Code:
type outletInlet;
outletValue uniform 0;
value uniform 0;
but it didn't help either.
Please suggest me how to make it work .
