[mvoss]

hi,
by backflow i mean a flapping of the flow somewhere, not just @ the inlet/outlet.
if the flow is partial unsteady you won´t achieve a steady solution.
regarding the residual check i must admit that i never did such a thing within OF... only knowing that from bigger commercial codes... but mainly it is storing the variable of the residual (for e.g. k, eps, and p,u,v,w) in a volume field (with NO_READ, AUTO_WRITE) for each iteration or at least for every 10th. or smth. similar. You would end up with a bunch of data when storing it every iteration. My idea would be :starting from an already found solution and calculate it for 2 periods of the oscillating residuals (e.g. 40 iterations).


neewbie

[vkrastev]

Mmmmh, still sounds a bit strange to me (for this case, and by my knowledge, it would be quite common to reach a steady-state solution using a standard RANS approach, whether if it is physically correct or not)...Maybe a more "immediate" check should be to run a transient case and see what happens, am I right? However, I'll try to investigate this matter further.

Apart from this, anyone that could tell me something about the setting of the case? Maybe Alberto?

Thank you all

V.

[eugene]

Not sure it will help, but putting zeroGradient boundaries on U on the side boundaries is a mistake. You will probably get inflow which will cause you all kinds of problems.

I suggest you set U to type slip on the top and side boundaries and set all other variables to the same as walls. Alternatively, you can just set the top and side boundaries to symmetry and get more-or-less the same effect.

[vkrastev]

Originally Posted by eugene

Not sure it will help, but putting zeroGradient boundaries on U on the side boundaries is a mistake. You will probably get inflow which will cause you all kinds of problems.

I suggest you set U to type slip on the top and side boundaries and set all other variables to the same as walls. Alternatively, you can just set the top and side boundaries to symmetry and get more-or-less the same effect.

Thanks for the hint, now I'm doing some runs with the slip condition for U at the top and side boundaries: I'll let you know if there will be some improvement in convergence patterns.

V.

[vkrastev]

Originally Posted by eugene

Not sure it will help, but putting zeroGradient boundaries on U on the side boundaries is a mistake. You will probably get inflow which will cause you all kinds of problems.

I suggest you set U to type slip on the top and side boundaries and set all other variables to the same as walls. Alternatively, you can just set the top and side boundaries to symmetry and get more-or-less the same effect.

Ok, after several runs, in which I have tried to set only the divScheme for U to te linearUpwindV option (the divSchemes for k and epsilon are left to be upwind), that is what I have obtained:

-Limiting both the general gradSchemes and the gradient option for the linearUpwindV scheme, a "kind of convergence" seems to be achieved.

-However, this convergence is reached only with the limited (in this case cellMDLimited) leastSquares choice (both for general gradSchemes and for the linearUpwindV <gradScheme>): using the limited Gauss linear option, the residuals keep to oscillate till the end of time...

-In addition, the "somewhat converging" runs gave a quite unphysical result for the nose-drag contribution of the body: in particular, after a few hundred of timesteps, the drag coefficient became negative (and I'm quite sure this is not a bug of the forceCoeffs function, because I've tested it before and it works properly)

For more clarity, I'll post down below the settings, residuals patterns and Cdrag patterns for the two "best converging" runs:

Run14 settings

fvSchemes:

ddtSchemes
{
default steadyState;
}

gradSchemes
{
default cellMDLimited leastSquares 1;
// grad(p) leastSquares;
// grad(U) leastSquares;
}

divSchemes
{
default none;
div(phi,U) Gauss linearUpwindV cellMDLimited leastSquares 1;
div(phi,k) Gauss upwind;
div(phi,epsilon) Gauss upwind;
div(R) Gauss linear;
div((nuEff*dev(grad(U).T()))) Gauss linear;
}

laplacianSchemes
{
default Gauss linear limited 0.5;
// laplacian(nuEff,U) Gauss linear corrected;
// laplacian((1|A(U)),p) Gauss linear corrected;
// laplacian(DkEff,k) Gauss linear corrected;
// laplacian(DepsilonEff,epsilon) Gauss linear corrected;
}

interpolationSchemes
{
default linear;
interpolate(U) linear;
}

snGradSchemes
{
default limited 0.5;
}

fluxRequired
{
default no;
p ;
}

// ************************************************** *********************** //

fvSolution:

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

solvers
{
// p
// {
// solver PCG;
// preconditioner DIC;
// tolerance 1e-11;
// relTol 0.001;
// }

p
{
solver GAMG;
tolerance 1e-12;
relTol 0.0001;
smoother GaussSeidel;
cacheAgglomeration true;
nCellsInCoarsestLevel 1000;
agglomerator faceAreaPair;
mergeLevels 1;
}

U
{
solver PBiCG;
preconditioner DILU;
tolerance 1e-10;
relTol 0.1;
}

k
{
solver PBiCG;
preconditioner DILU;
tolerance 1e-10;
relTol 0.1;
}

epsilon
{
solver PBiCG;
preconditioner DILU;
tolerance 1e-10;
relTol 0.1;
}

}

SIMPLE
{
nNonOrthogonalCorrectors 3;
}

relaxationFactors
{
p 0.3;
U 0.7;
k 0.7;
epsilon 0.7;
}


// ************************************************** *********************** //

0/U:

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

internalField uniform (40 0 0);

boundaryField
{
inlet
{
type fixedValue;
value uniform (40 0 0);
}

outlet
{
type zeroGradient;
}

body
{
type fixedValue;
value uniform (0 0 0);
}

nose
{
type fixedValue;
value uniform (0 0 0);
}

slant
{
type fixedValue;
value uniform (0 0 0);
}

back
{
type fixedValue;
value uniform (0 0 0);
}

floor
{
type fixedValue;
value uniform (0 0 0);
}

side
{
type slip;
}

top
{
type slip;
}

symmetry
{
type symmetryPlane;
}
}

// ************************************************** *********************** //

0/p:

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

internalField uniform 0;

boundaryField
{
inlet
{
type zeroGradient;
}

outlet
{
type fixedValue;
value uniform 0;
}

body
{
type zeroGradient;
}

nose
{
type zeroGradient;
}

slant
{
type zeroGradient;
}

back
{
type zeroGradient;
}

floor
{
type zeroGradient;
}

symmetry
{
type symmetryPlane;
}

side
{
type zeroGradient;
}

top
{
type zeroGradient;
}
}

// ************************************************** *********************** //

0/k:

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

internalField uniform 0.0096;

boundaryField
{
inlet
{
type fixedValue;
value uniform 0.0096;
}
outlet
{
type zeroGradient;
}
body
{
type kqRWallFunction;
value uniform 0.0096;
}

nose
{
type kqRWallFunction;
value uniform 0.0096;
}

slant
{
type kqRWallFunction;
value uniform 0.0096;
}

back
{
type kqRWallFunction;
value uniform 0.0096;
}

floor
{
type kqRWallFunction;
value uniform 0.0096;
}

side
{
type zeroGradient;
}

top
{
type zeroGradient;
}

symmetry
{
type symmetryPlane;
}
}

// ************************************************** *********************** //

0/epsilon:

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

internalField uniform 0.155;

boundaryField
{
inlet
{
type fixedValue;
value uniform 0.00155;
}
outlet
{
type zeroGradient;
}
body
{
type epsilonWallFunction;
value uniform 0.00155;
}
nose
{
type epsilonWallFunction;
value uniform 0.00155;
}
slant
{
type epsilonWallFunction;
value uniform 0.00155;
}
back
{
type epsilonWallFunction;
value uniform 0.00155;
}
floor
{
type epsilonWallFunction;
value uniform 0.00155;
}

side
{
type zeroGradient;
}

top
{
type zeroGradient;
}

symmetry
{
type symmetryPlane;
}
}

// ************************************************** *********************** //

0/nut:

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

internalField uniform 0;

boundaryField
{
inlet
{
type calculated;
value uniform 0;
}
outlet
{
type calculated;
value uniform 0;
}
body
{
type nutWallFunction;
value uniform 0;
}
nose
{
type nutWallFunction;
value uniform 0;
}
slant
{
type nutWallFunction;
value uniform 0;
}
back
{
type nutWallFunction;
value uniform 0;
}
floor
{
type nutWallFunction;
value uniform 0;
}

side
{
type zeroGradient;
value uniform 0;
}

top
{
type zeroGradient;
value uniform 0;
}

symmetry
{
type symmetryPlane;
}
}

// ************************************************** *********************** //

checkMesh report:

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create polyMesh for time = 0

Time = 0

Mesh stats
points: 350035
faces: 3138087
internal faces: 3025804
cells: 1488208
boundary patches: 10
point zones: 0
face zones: 0
cell zones: 0

Overall number of cells of each type:
hexahedra: 0
prisms: 211059
wedges: 0
pyramids: 0
tet wedges: 0
tetrahedra: 1277149
polyhedra: 0

Checking topology...
Boundary definition 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
symmetry 22584 12567 ok (non-closed singly connected)
top 1814 994 ok (non-closed singly connected)
side 2602 1393 ok (non-closed singly connected)
outlet 372 213 ok (non-closed singly connected)
inlet 384 219 ok (non-closed singly connected)
floor 14174 7364 ok (non-closed singly connected)
nose 8710 4464 ok (non-closed singly connected)
body 54795 27714 ok (non-closed singly connected)
back 2882 1513 ok (non-closed singly connected)
slant 3966 2067 ok (non-closed singly connected)

Checking geometry...
Overall domain bounding box (-1.26 0 0.1945) (6.36 1.46 1.1945)
Mesh (non-empty, non-wedge) directions (1 1 1)
Mesh (non-empty) directions (1 1 1)
Boundary openness (-1.30596e-19 3.09526e-19 7.90637e-19) OK.
Max cell openness = 2.7631e-16 OK.
Max aspect ratio = 5.38151 OK.
Minumum face area = 1.25535e-06. Maximum face area = 0.0132756. Face area magnitudes OK.
Min volume = 6.89797e-10. Max volume = 0.000546044. Total volume = 11.0702. Cell volumes OK.
Mesh non-orthogonality Max: 54.7518 average: 14.569
Non-orthogonality check OK.
Face pyramids OK.
Max skewness = 2.0193 OK.

Mesh OK.

Additional notes:

OpenFOAM version: 1.6
turbulence model: realizableKE


Run15 settings

The same as Run14, but:

-with URF's for k and epsilon set to 0.6 instead of 0.7
-with relTol parameters for U, k and epsilon lowered to 10e-03 instead of 10e-01 (and with the smooth solver option instead of thr DILU-PBiCG solver one)

Regards

V.

PS-The red-dot horizontal pattern corresponds to the experimental nose-drag measured by Ahmed et al. (1984)

[eugene]

This
default cellMDLimited leastSquares 1;

is too restrictive in my opinion. Especially since it is applied to the pressure gradient, which affects continuity. I have also heard that the leastSquares implementation in standard OPENFOAM has some problems.

In any case, it is very odd that your result is almost exactly the experimental value times -1. I would check the forces FO and the flow in the vicinity of the object a few more times.

[alberto]

Originally Posted by eugene

This
default cellMDLimited leastSquares 1;

is too restrictive in my opinion. Especially since it is applied to the pressure gradient, which affects continuity. I have also heard that the leastSquares implementation in standard OPENFOAM has some problems.

What problems for the chronicles? :-)
It would be useful to have this reported as a bug report, if there is some way of showing it.

P.S. With restrictive you mean the limiter should not be applied to grad(p)?

[eugene]

Hi Alberto,

Check the dev version for the modified least-squares implementation. I have not studied the issue in detail. For now it is just an unsubstantiated comment, as I can't remember the specifics.

Correct, the limiter should not be applied to grad(p) if you can help it. (Sometimes it is better to get a less accurate result than no result at all.)

[alberto]

I have to look into that more carefully, however the least-squares implementations for scalars identical (a minor difference in a constructor), but not for vectors, where the weighting factors are different (OpenCFD release uses a surface-based weighting).

In this post Henry discussed the implementation:

http://www.cfd-online.com/Forums/ope...stsquares.html

The two implementations of extendedLeastSquares seem consistent in both 1.7.x and 1.5-dev (no surface-based weights. There are a few other differences, unrelated).

Best,

[eugene]

Hmm, I never noticed it. Thanks.

[vkrastev]

Originally Posted by eugene

This
default cellMDLimited leastSquares 1;

is too restrictive in my opinion. Especially since it is applied to the pressure gradient, which affects continuity.

I've removed the limiter on grad(p), but as you can see from the results posted below now the solution doesn't even seem to converge (the Run18 settings are the same as Run15 but with only leastSquares scheme on grad(p))

Originally Posted by eugene

I have also heard that the leastSquares implementation in standard OPENFOAM has some problems.

Should I use extendedLeastSquares or Gauss linear instead? And should I do it also for the gradScheme option in the linearUpwindV div scheme?

Originally Posted by eugene

In any case, it is very odd that your result is almost exactly the experimental value times -1. I would check the forces FO and the flow in the vicinity of the object a few more times.

Don't get misleaded about this value: the forceCoeffs function has always been working perfectly before these cases (the results for the round cylinder are in perfect agreement with the experimental data). Also, for the same case, but with Gauss upwind divSchemes, it gave absolutely reasonable results for the drag coefficients (both the total one and the partial contribution for the nose), wich are quite higher than the experimental ones but very close to other RANS calculations from the literature.

Originally Posted by alberto

What problems for the chronicles? :-)

For me, the problems ere that still I cannot find a way to get my runs working on with higher order divSchemes

Originally Posted by alberto

The two implementations of extendedLeastSquares seem consistent in both 1.7.x and 1.5-dev (no surface-based weights. There are a few other differences, unrelated).

Again, should it help if I use extendedLeastSquares (OF 1.6) instead of leastSquares?

Thank you

V.

[alberto]

Originally Posted by vkrastev

About to check where the residuals are higher in the domain, can you give me a hint on how it could be done? Thank you

The easy, but a bit messy, way is to

1. Define a volScalarField for p, called say pRes, and a volVectorField for U, uRes.
2. Set pRes = fvc::.... your pressure equation.
3. Set uRes = fvc::... your momentum equation.

With fvc::... I mean that all the terms in the momentum and pressure equation must be computed explicitly. The fields pRes and uRes will contain the residuals.

Best,

[alberto]

Originally Posted by vkrastev

I've removed the limiter on grad(p), but as you can see from the results posted below now the solution doesn't even seem to converge (the Run18 settings are the same as Run15 but with only leastSquares scheme on grad(p))

My nose smells a mesh problem somewhere

Should I use extendedLeastSquares or Gauss linear instead? And should I do it also for the gradScheme option in the linearUpwindV div scheme?

Yes! The schemes used for gradients must be consistent.

For me, the problems ere that still I cannot find a way to get my runs working on with higher order divSchemes

Again, should it help if I use extendedLeastSquares (OF 1.6) instead of leastSquares?

Only if your mesh is very, very bad.

Best,

[mvoss]

hi,
thanks for the answer.

Originally Posted by alberto

The easy, but a bit messy, way is to

1. Define a volScalarField for p, called say pRes, and a volVectorField for U, uRes.
2. Set pRes = fvc::.... your pressure equation.
3. Set uRes = fvc::... your momentum equation.

With fvc::... I mean that all the terms in the momentum and pressure equation must be computed explicitly. The fields pRes and uRes will contain the residuals.

I don´t get it.
fake-numerically spoken --> fvc:: means to calculate U from smth. known, in this case the momentum eq. with already known U´s.
So if i fvc smth. i end up in the explicitly calculated value for each node from the iterated (fvm)-solution. I understand that if i calculate smth. from a known rhs i end up with my "failure"(smth. like a residual) on the lhs.
Is this correct so fare ?
A residual is smth. like res=abs(U_new-U_old)/abs(U_old) and all i have to do is to normalize the Ures values ?

Sorry for the messy english.

neewbie

[alberto]

The uRes you compute with fvc is the "distance of UEqn from zero", meaning:

UEqn(U) = uRes

where UEqn contains all its terms.

Best,

[mvoss]

hi,

by "distance" from zero you mean lhs-rhs= distance ?

I tried the approach above and got extremely high values in my uRes field, that´s why i was thinking about normalization.

neewbie

[Eren10]

I am using k-epsilon model. I have set everything to second order scheme, unless the k and epsilon terms. When I make them also second order the simulation blows up.

Why can't k and epsilon be seconder order ? ( I am not sure if both can't be second order)

[makaveli_lcf]

They can, but since 2nd order produces numerical oscillations you can try to use underelaxation for k and epsilon to stabilize it.

[s.m]

Originally Posted by marine

OK that explain why the linear scheme doesn't work, my Peclet number is too big.

I don't obtain the same results for the viscous forces on my hull with the linearUpwind scheme or the limitedLinear scheme (15% difference) . As they say in the book (Peric) that linearUpwind is unbounded I think is more accurate than limitedLinear ( I still don't know how this one works) and is comparable to the second order upwinded scheme we can find in Starccm+ or Fluent. Do you think I'm wrong?

regards,

Marine

Hi dear Marine,

how we should calculate the "local (cell) Peclet number" ?
thank you

[akidess]

Use the "Pe" tool that comes with OF-2.2.x (earlier versions were buggy).