[Chrisstiapis]

Dear al,
My name is Christos and I'm currently trying to simulate a distribution system for my PhD.
The distributions system has an inlet and 4 outlets. And I'm using simpleFoam.
The inlet is at the bottom and the outlets on the right.


I meshed the Geometry with the use of Salome and I got this mesh (on a branch) which is not bad on my opinion


The problem is that when I;m using simpleFoam with kEpsilon the time-continuity explodes. When I'm using realizablekE timestep error increases but but slowly. When I'm using k-Omega again the timestep explodes.
I'm attaching the 0 folder and the system folder so if anyone wants may have a look.
Mesh orthogonality 50 and y+37.

The velocity boundary conditions
boundaryField
{
inlet_extruded
{
type fixedValue;
value uniform (0 0 0.44211);
}
outlet_extruded
{
type inletOutlet;
inletValue uniform (0 0 0);
value uniform (0 0 0);
}
walls_extruded
{
type fixedValue;
value uniform (0 0 0);
}
frontAndBack
{
type empty;
}
} //

And the pressure
boundaryField
{
inlet_extruded
{ type zeroGradient;
}
outlet_extruded
{
type fixedValue;
value uniform 44500;
}
walls_extruded
{
type zeroGradient;
}
frontAndBack
{
type empty;
}
} //

The fvSolution
solvers
{
p
{
solver GAMG;
tolerance 1e-08;
relTol 0.00005;
smoother GaussSeidel;
cacheAgglomeration true;
nCellsInCoarsestLevel 20;
agglomerator faceAreaPair;
mergeLevels 1;
}
U
{
solver smoothSolver;
smoother GaussSeidel;
nSweeps 2;
tolerance 1e-07;
relTol 0.00005;
}

k
{
solver smoothSolver;
smoother GaussSeidel;
nSweeps 2;
tolerance 1e-07;
relTol 0.01;
}

epsilon
{
solver smoothSolver;
smoother GaussSeidel;
nSweeps 2;
tolerance 1e-07;
relTol 0.1;
}

}

SIMPLE
{
nNonOrthogonalCorrectors 1;
pRefCell 0;
pRefValue 0;
}

relaxationFactors
{
fields
{
p 0.3;
}
equations
{
U 0.5;
k 0.5;
epsilon 0.5;
omega 0.5;
}
}

I really don't uderstand what is happening ????????


I get

GAMG: Solving for p, Initial residual = 0.276687, Final residual = 8.57431e-06, No Iterations 14
smoothSolver: Solving for Uy, Initial residual = 0.00637261, Final residual = 5.2594e-08, No Iterations 8
smoothSolver: Solving for Uz, Initial residual = 0.0119267, Final residual = 6.9414e-08, No Iterations 8
GAMG: Solving for p, Initial residual = 0.0189897, Final residual = 4.52346e-07, No Iterations 15
time step continuity errors : sum local = 0.255507, global = 0.00466411, cumulative = 6.20251
smoothSolver: Solving for epsilon, Initial residual = 0.0277328, Final residual = 0.00100213, No Iterations 2
bounding epsilon, min: -1.56859e+12 max: 3.84711e+13 average: 3.64094e+10
GAMG: Solving for p, Initial residual = 0.0189897, Final residual = 4.52346e-07, No Iterations 15
time step continuity errors : sum local = 0.255507, global = 0.00466411, cumulative = 6.20251
smoothSolver: Solving for k, Initial residual = 3.51264e-12, Final residual = 3.51264e-12, No Iterations 0
ExecutionTime = 2561.74 s ClockTime = 3104 s

[alexeym]

Hi,

and what are ICs and BCs for k, epsilon, and omega? Can you also post fvSchemes?

Also, looking at your geometry I can't grasp why do you need triangular mesh.

[Chrisstiapis]

Alexey thanks for replying,

I'm using triangular mesh cause this geometry and mesh is generated by an automated routine in Salome, where i use the salome functionality to batch generate geometry and meshes.
I tried and experimented with many schemes and now I'm using that one that is considered very stable although diffusive.
I took a look in the bibliography and it states that in case of multiple outlets only a velocity boundary at the inlet needs to be specified and pressure at the outlet.



The velocity U

boundaryField
{
inlet_extruded
{
type fixedValue;
value uniform (0 0 0.44);
}
outlet_extruded
{
type zeroGradient;
}
walls_extruded
{
type fixedValue;
value uniform (0 0 0);
}
frontAndBack
{
type empty;
}

And the pressure P
dimensions [ 0 2 -2 0 0 0 0 ];

internalField uniform 44900;

boundaryField
{
inlet_extruded
{ type zeroGradient;
}
outlet_extruded
{
type fixedValue;
value uniform 44900;
}
walls_extruded
{
type zeroGradient;
}
frontAndBack
{
type empty;
}
}

Also the fVSchemes

ddtSchemes
{
default steadyState;
}

gradSchemes
{
default cellLimited leastSquares 1.0;
}

divSchemes
{
default none;
div(phi,U) bounded Gauss linearUpwind cellLimited Gauss linear 1;
div(phi,k) bounded Gauss upwind;
div(phi,epsilon) bounded Gauss upwind;
div(phi,R) bounded Gauss upwind;
div(R) Gauss linear;
div(phi,nuTilda) bounded Gauss upwind;
div((nuEff*dev(T(grad(U))))) Gauss linear;
div(phi,omega) bounded Gauss upwind;
}

laplacianSchemes
{
default Gauss linear limited 0.5;
laplacian(1,p) Gauss linear corrected;
}

interpolationSchemes
{
default linear;
}

snGradSchemes
{
default limited 0.333;
}

fluxRequired
{
default no;
p ;
}

[alexeym]

According to your log file, epsilon equation is a source of the problems:

time step continuity errors : sum local = 0.255507, global = 0.00466411, cumulative = 6.20251
smoothSolver: Solving for epsilon, Initial residual = 0.0277328, Final residual = 0.00100213, No Iterations 2
bounding epsilon, min: -1.56859e+12 max: 3.84711e+13 average: 3.64094e+10

check ICs and BCs for this variable. Well, in fact just ICs as boundary conditions seem to be reasonable.

Also if you know pressure are the outlet, why not use pressureDirectedInletOutletVelocity BC for velocity? And are you sure all your outlets will have the same pressure?

[Chrisstiapis]

Dear Alexey, this distribution system is in the bottom of a vesssel(not pressured) so the pressure on the outlet is the hydrostatic pressure and as long as the outlets share the same height i suppose that the pressure will be equal.

I decided to neglect the turbulence Effects till i get my case converged in laminar.
I did as u suggested and I'm using the pressureDirectedInletOutletVelocity condition for velocity and the the total pressure. I read somewhere in this forum that this coupling is stable.
I'm getting again a explosion in time-step continuity ?????

pressureDirectedInletOutletVelocity

smoothSolver: Solving for Uy, Initial residual = 0.243622, Final residual = 1.00739e-08, No Iterations 13
smoothSolver: Solving for Uz, Initial residual = 0.271367, Final residual = 1.00003e-08, No Iterations 13
GAMG: Solving for p, Initial residual = 0.410951, Final residual = 4.08045e-06, No Iterations 413
time step continuity errors : sum local = 2.12258, global = 0.445773, cumulative = -0.296711.
The velocity boundary now is
dimensions [ 0 1 -1 0 0 0 0 ];

internalField uniform (0 0 0);

boundaryField
{
inlet_extruded
{
type pressureDirectedInletOutletVelocity;
inletDirection uniform (0 0 1);
value uniform (0 0 0.44);

}
outlet_extruded
{
type zeroGradient;

}
walls_extruded
{
type fixedValue;
value uniform (0 0 0);
}
frontAndBack
{
type empty;
}
} //

and the pressure
dimensions [ 0 2 -2 0 0 0 0 ];

internalField uniform 44900;

boundaryField
{
inlet_extruded
{ type totalPressure;
U U;
phi phi;
rho none;
psi none;
gamma 1.4;
p0 uniform 244900;
}
outlet_extruded
{
type fixedValue;
value uniform 44900;
}
walls_extruded
{
type zeroGradient;
}
frontAndBack
{
type empty;
}
} //

?????

[alexeym]

Well,

If your case is laminar, why even bother about turbulence model? If your case is turbulent, it won't converge without turbulence model.

We can continue to play in a game called "guess real physical conditions in your system" (to suggest right boundary conditions), or you can try to describe in more details what you're trying to model. Right now your BCs don't make sense.

[Chrisstiapis]

Thank you Alexey,
My case is turbulent (RE>50000).
What i was saying earlier was to minimize the mesh sensitivity by minimizing the velocity so it will be laminar in order to test the boundary conditions without worrying about the mesh sensitivity.
Now, my case.
There is a distribution system in the bottom of a vessel. Therefore the pressure on the outlet is known from the hydrostatic pressure.
The pressure on the inlet is known, cause there is a pump. The velocity at the inlet is known because of the flowrate of the pump.
Can I ask you, what do u think is wrong on the previously described boundary conditions???
My epsilon BCS

internalField uniform 0.00000001;

boundaryField
{
inlet_extruded
{
type fixedValue;
value $internalField;
}
outlet_extruded
{
type zeroGradient;
}
walls_extruded
{
type epsilonWallFunction;
value $internalField;
}

}
frontAndBack
{
type empty;
}

And My K boundary conditions
internalField uniform 0.0000001;

boundaryField
{
inlet_extruded
{
type fixedValue;
value uniform 0.00750571558131;
}
outlet_extruded
{
type zeroGradient;
}
walls_extruded
{
type kqRWallFunction;
value $internalField;
}
frontAndBack
{
type empty;
}
}

although at some point converges due to the very loose tolerances

[RodriguezFatz]

Chris, you use the wrong syntax for the under-relaxation factors. Accordingly, OpenFOAM doesn't recognize them. You need to put

relaxationFactors
{
    fields
    {
        "(p)"           0.3;
    }
    equations
    {
        "(U)"           0.7;
        "(k)"           0.8;
        "(epsilon)"       0.8;
    }
}

into your fvSolution. Simple without under-relaxation will rarely converge, as you probably know.
If you post code, hit the "Go advanced" button below the text box and use the "#" sign, to wrap "code" tags. This is much easier to read. Also some tabs will help to understand the code.

Try this and post again, if it worked or not.

[Chrisstiapis]

Thank you guys for replying, really appreciated.
I changed it and it doesn't seem to make a difference. I think that in OpenFoam 2.3.0 the previous systax is allowed as well as the one suggested.
No I'm doing the simulation and it seems to converge very slowly,

#Time = 159

smoothSolver: Solving for Uy, Initial residual = 0.000998984, Final residual = 5.55782e-10, No Iterations 14
smoothSolver: Solving for Uz, Initial residual = 0.000555796, Final residual = 7.87531e-10, No Iterations 13
GAMG: Solving for p, Initial residual = 0.434494, Final residual = 9.3188e-09, No Iterations 1000
time step continuity errors : sum local = 9.16547e-08, global = 2.3667e-10, cumulative = 2.67703e-09
smoothSolver: Solving for epsilon, Initial residual = 0.000423026, Final residual = 4.34165e-10, No Iterations 9
smoothSolver: Solving for k, Initial residual = 8.5116e-10, Final residual = 8.5116e-10, No Iterations 0
ExecutionTime = 44.03 s ClockTime = 47 s

Time = 160

smoothSolver: Solving for Uy, Initial residual = 0.00092316, Final residual = 4.05633e-10, No Iterations 14
smoothSolver: Solving for Uz, Initial residual = 0.000517306, Final residual = 7.15796e-10, No Iterations 13
GAMG: Solving for p, Initial residual = 0.434494, Final residual = 9.3188e-09, No Iterations 1000
time step continuity errors : sum local = 9.16547e-08, global = 2.3667e-10, cumulative = 2.67703e-09
smoothSolver: Solving for epsilon, Initial residual = 0.000423026, Final residual = 4.34165e-10, No Iterations 9
GAMG: Solving for p, Initial residual = 0.434494, Final residual = 9.3188e-09, No Iterations 1000
time step continuity errors : sum local = 9.16547e-08, global = 2.3667e-10, cumulative = 2.67703e-09
smoothSolver: Solving for k, Initial residual = 8.5116e-10, Final residual = 8.5116e-10, No Iterations 0
ExecutionTime = 45.38 s ClockTime = 48 s
#

I beievee that this convergence is due to the very stict tolerances that I have imposed in fvSolution

#p
{
solver GAMG;
tolerance 6e-10;
relTol 0.00000001;
smoother GaussSeidel;
nPreSweeps 0;
nPostSweeps 2;
cacheAgglomeration on;
agglomerator faceAreaPair;
nCellsInCoarsestLevel 10;
mergeLevels 1;
}

"(U|k|epsilon|omega|R|nuTilda)"
{
solver smoothSolver;
smoother symGaussSeidel;
tolerance 1e-09;
relTol 0.0000001;
}#

I really don't know what is going wrong. I neverthought that this case will get so complicated .
The turbulence mdoel I'm using is realizableKE.

[RodriguezFatz]

"relTol" normally doesn't need to be lower than 0.1 or 0.01 in simple. The non-linear error of the outer iterations is much larger anyway. In any case, 0.00000001 is insane
I think for SIMPLE, you should try to have just a few iterations in all equations, such as maximum 5.

Also, you don't actually need the non-orthogonal correction, as you mesh isn't that bad. Hrvoje Jasak stated somewhere in this forum, that you should only use them, if your mesh orthogonal quality is larger than 70°.

In your fvSchemes:

• start with "Gauss linear" gradient scheme
• with "bounded Gauss upwind" for all convective schemes, except the "div((nuEff*dev(T(grad(U))))) Gauss linear;".
• set "default Gauss linear uncorrected;" for all laplacian schemes. This is the most stable.

If these settings diverge, your initial settings are probably bad.

[RodriguezFatz]

Originally Posted by Chrisstiapis

I changed it and it doesn't seem to make a difference. I think that in OpenFoam 2.3.0 the previous systax is allowed as well as the one suggested.

Edit: Ok, you are right, my fault.

[alexeym]

Hi,

Just as PoC I've made similar case and it converges quite quickly. My BCs (I've decided not to mess with turbulence, so it's just p and U):

inlet:
U -> fixedValue
p -> zeroGradient

outlet:
U -> zeroGradient (as I don't expect back-flow, otherwise something more fancy can be used)
p -> fixedValue (as you've said)

walls:
U -> non-slip
p -> zeroGradient

You can find a case attached to the message.

Concerning your turbulence BCs - where did you get these constants for k and epsilon? RNG? Experimental values? Estimations from turbulence intensity and hydraulic radius?

[RodriguezFatz]

I can highly recommend to use Alexey's settings. These are the most basic / common settings for "p" and "u" and thus (let's hope) the most stable ones. I use them for all my pipe simulations, also in Fluent. I would not try to use any others, just because some person in the forum told you so, except you really know what you are doing.

Regarding my upper post, I just found this: http://www.openfoam.org/mantisbt/view.php?id=1410#c3246
A statement of henry: Also when running the case I improved the schemes, in particular it is a VERY bad idea to limit all gradients: the pressure gradient should NEVER be limited.

[Chrisstiapis]

For the inlet
turbulent intensity
I=0.16Re−1/8. The turbulent length scale can be estimated as


k=32(UI)2,
ϵ=cμ3/4k3/2l−1.andl=0.07L, with L a characteristic length.
For internal flows this may take the value of the inlet duct (or pipe) width (or diameter) or the hydraulic diameter.
This is what I'm using in order to calculate the
l=0.07L

And the k file

 internalField uniform 0.075;
boundaryField
{
  inlet_extruded
  {
    type fixedValue;
    value uniform 0.075;
    intensity 0.160000;

for epsilon

internalField uniform 0.02;
boundaryField
{
  inlet_extruded
  {
            type            turbulentIntensityKineticEnergyInlet;
        intensity       0.16;       // 16% turbulent intensity
        value           $internalField;
  }

internalField uniform 1;

boundaryField
{
inlet_extruded
{
type turbulentMixingLengthDissipationRateInlet;
mixingLength 0.02; // 0.02m - half pipe diameter
value $internalField;
}

The weird thing is that it converges although the time-step continuity is high???

Time = 3072

smoothSolver:  Solving for Uy, Initial residual = 0.000147551, Final residual = 8.79075e-07, No Iterations 7
smoothSolver:  Solving for Uz, Initial residual = 0.00014121, Final residual = 7.46354e-07, No Iterations 7
GAMG:  Solving for p, Initial residual = 0.00998099, Final residual = 9.77039e-05, No Iterations 7
time step continuity errors : sum local = 0.000324409, global = 4.5041e-06, cumulative = -0.0117465
smoothSolver:  Solving for epsilon, Initial residual = 3.33489e-05, Final residual = 1.33563e-07, No Iterations 4
smoothSolver:  Solving for k, Initial residual = 8.77387e-10, Final residual = 8.77387e-10, No Iterations 0
ExecutionTime = 584.71 s  ClockTime = 586 s


SIMPLE solution converged in 3072 iterations

[Chrisstiapis]

The problem that I'm suspecting now,as alexey was saying, maybe there is an issue with the IC.
I'm a starter now in CFD so I'm a little bit lost

[alexeym]

Well, if ICs affect your steady state solution, I guess there's something wrong with the solution

If you take a look at src/finiteVolume/cfdTools/incompressible/continuityErrs.H:

{
    volScalarField contErr(fvc::div(phi));

    scalar sumLocalContErr = runTime.deltaTValue()*
        mag(contErr)().weightedAverage(mesh.V()).value();

    scalar globalContErr = runTime.deltaTValue()*
        contErr.weightedAverage(mesh.V()).value();
    cumulativeContErr += globalContErr;

    Info<< "time step continuity errors : sum local = " << sumLocalContErr
        << ", global = " << globalContErr
        << ", cumulative = " << cumulativeContErr
        << endl;
}

So time step continuity errors depend on fvc::div(phi) and also phi depends on cell->face interpolation, guess if you've got ugly mesh contErr can be rather large.

[RodriguezFatz]

1) Did you set the schemes as I suggested in post # 10?
2) You can increase the convergence criteria? Now it says it converged just because all initial residuals are below 1e-3. Set them to 1e-12 and see how low they can get.

[Chrisstiapis]

I did as u suggested, that's why took so long to answer,many iterations,
Although, the continuity has been greatly reduced the cumulative error remains high

Time = 13146
smoothSolver:  Solving for Uz, Initial residual = 3.10449e-06, Final residual = 2.44747e-09, No Iterations 8
GAMG:  Solving for p, Initial residual = 8.14567e-05, Final residual = 7.23971e-08, No Iterations 10
time step continuity errors : sum local = 1.39972e-07, global = 4.68047e-09, cumulative = -0.000204088
smoothSolver:  Solving for Uy, Initial residual = 2.03736e-06, Final residual = 1.46767e-09, No Iterations 8
smoothSolver:  Solving for Uz, Initial residual = 3.13649e-06, Final residual = 2.42425e-09, No Iterations 8
smoothSolver:  Solving for epsilon, Initial residual = 5.74286e-07, Final residual = 2.65494e-10, No Iterations 6
smoothSolver:  Solving for Uz, Initial residual = 3.09894e-06, Final residual = 2.41098e-09, No Iterations 8
smoothSolver:  Solving for Uy, Initial residual = 2.05149e-06, Final residual = 1.42138e-09, No Iterations 8
smoothSolver:  Solving for Uz, Initial residual = 3.06732e-06, Final residual = 2.39964e-09, No Iterations 8
smoothSolver:  Solving for epsilon, Initial residual = 5.66788e-07, Final residual = 2.62298e-10, No Iterations 6
smoothSolver:  Solving for Uy, Initial residual = 2.02899e-06, Final residual = 1.45363e-09, No Iterations 8
smoothSolver:  Solving for Uy, Initial residual = 2.091e-06, Final residual = 1.47467e-09, No Iterations 8
smoothSolver:  Solving for k, Initial residual = 5.72958e-07, Final residual = 4.05752e-10, No Iterations 6
smoothSolver:  Solving for Uz, Initial residual = 3.08764e-06, Final residual = 2.3985e-09, No Iterations 8

Additionally, what do u think about the boundary conditions in k and epsilon???

[RodriguezFatz]

1) Chris, to get the inlet profiles, I usually run a case with just a short pipe / channel with periodic boundary conditions with the shape of the inlet of my second case. Then, I use "mapFields" to map the outlet of my first case (U, k, omega, epsilon, nu) to the inlet of my second case. Now, I have physically meaningful inlet conditions. See this post for a sketch:
http://www.cfd-online.com/Forums/openfoam-pre-processing/141196-mapfields-question.html#post508830

2) For the other problem: I have never seen, that the cumulative error remains high, while all others decrease. What pressure and velocity inlet and outlet b.c. do you use right now?

3) If it converges, you can also start to change convective schemes to "linearUpwind". If you get numerical problems you can set a limiter, but just for the upwind gradient scheme:

gradSchemes
{
    default             Gauss linear;
    grad(linUpwind)    faceLimited Gauss linear 1.0;
}

divSchemes
{
       div((nuEff*dev(T(grad(U))))) Gauss linear;
        div(phi,U)        Gauss linearUpwind grad(linUpwind);
        div(phi,k)        Gauss linearUpwind grad(linUpwind);
        div(phi,epsilon)  Gauss linearUpwind grad(linUpwind);
}

[Chrisstiapis]

My velocity and Pressure

U

FoamFile
{
 version 2.0;
 format ascii;
 class volVectorField;
 object U;
}

dimensions [ 0 1 -1 0 0 0 0 ];

internalField uniform (0 0 0);

boundaryField
{
  inlet_extruded
  {
    type fixedValue;
    value uniform (0 0 0.44211);
  }
  outlet_extruded
  {
    type zeroGradient;
  }
  walls_extruded
  {
    type fixedValue;
    value uniform (0 0 0);
  }
  frontAndBack
  {
    type empty;
  }
}     // **************************************

and P

dimensions      [0 2 -2 0 0 0 0];

internalField   uniform 49000;

boundaryField
{
    inlet_extruded
    {
        type            zeroGradient;
    }

    outlet_extruded
    {
        type            fixedValue;
        value           uniform 49000;
    }

    walls_extruded
    {
        type            zeroGradient;
    }


    frontAndBack
    {
        type            empty;
    }
}

A problem that I was facing was that although i know the pressure on the outlet cannot use it cause whenever I was putting fixed value or total pressure the solution was becoming divergent
This procedure with the map-fields looks promising although i have never used this function before.