[hani]

I have OpenFOAM converging a cyclic water turbine runner computation. I am however not pleased with the convergence rate, and I think that my setup of the case could be improved (and I could of course use your ideas)

My case is shown in the following figure. It shows the cyclic 1/5 of the runner which is cyclic about the y-axis. The geometry is colored gray. The cyclic boundaries cut through the center of the blades, and they are colored by the turbulent kinetic energy.

The flow comes in at the top, where I set steady highly swirling profiles for U,k,epsilon and dp/dn=0. At the walls I use wall functions, and I use the standard k-epsilon model. At the outlet I have tried different alternatives, with the following common settings: p=0, dU/dn=0, and the following alternatives: dk/dn=0, depsilon/dn=0 or inletOutlet on k and epsilon, with an inlet value similar to the average at the true inlet.

The computation converges nicely in the upper half of the computational domain, but below the runner (in the diffusing part) the convergence is extremely slow, particularly close to the rotational axis (in the region colored red below the hub). In part of this region there is a recirculation. There are several possible resons for this behaviour:
1) the flow in this region is not steady in reality (however, there does not seem to be any oscillations in the solution but the flow seems to be extremely slowly converging towards a solution),
2) the outlet boundary condition is not good enough (this is likely, I don't like the concept of constant pressure for outlets with swirling flow),
3) there is something wrong with my cyclics (in another thread I have described my problems using this mesh in OpenFOAM, and I had to change the too strict geometrical tolerance to make it work. I have however checked that the cyclic coupling looks alright),
4) there is a problem with the prism cells at the axis (they are generated due to the collapse of one of the faces in the hex-cells, the rest of the grid cells are hexas),
5) etc. etc. (please give me your ideas)

The residuals at the end of the computation are shown in the following figure:


The flow is computed with a modified simpleFoam, where the Coriolis and centrifugal forces are added as:
dimensionedVector Omega
(
"Omega",
dimensionSet(0, 0, -1, 0, 0),
vector(0, 2*M_PI/60*(-595), 0)
);
volVectorField Fcent = (Omega ^ (Omega ^ mesh.C()));
--------------------------------------
tmp<fvvectormatrix> UEqn
(
fvm::div(phi, U)
+ turbulence->divR(U)
+ (2*Omega ^ U)
+ Fcent
);

I hope that someone out there has any idea on what might cause this slow convergence below the hub.

Below you can see my fvSchemes, fvSolution

Håkan.

fvSchemes:
ddtSchemes
{
default steadyState;
}

gradSchemes
{
default Gauss linear;
grad(p) Gauss linear;
grad(U) Gauss linear;
}

divSchemes
{
default none;
div(phi,U) Gauss GammaV 1;
div(phi,k) Gauss Gamma 1;
div(phi,epsilon) Gauss Gamma 1;
// div(phi,U) Gauss upwind;
// div(phi,k) Gauss upwind;
// div(phi,epsilon) Gauss upwind;
div(phi,R) Gauss upwind;
div(R) Gauss linear;
div(phi,nuTilda) Gauss upwind;
div((nuEff*dev(grad(U).T()))) Gauss linear;
}

laplacianSchemes
{
default none;
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;
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;
}
----------------------------------------------
fvSolution:
solvers
{
// p ICCG 1e-06 0.01;
p AMG 1e-07 0.01 100;
U BICCG 1e-06 0.1;
k BICCG 1e-06 0.1;
epsilon BICCG 1e-06 0.1;
R BICCG 1e-06 0.1;
nuTilda BICCG 1e-06 0.1;
}

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

relaxationFactors
{
p 0.2;
U 0.7;
k 0.5;
epsilon 0.5;
R 0.5;
nuTilda 0.5;
}

[eugene]

Unless you have high non-orthogonality, your settings look fine. Might want to play with your convection schemes.

At a guess I would say your problem is due the erronous constant outlet pressure. At least that is the first thing I would check.

[hani]

Thanks for your response, Eugene.

Hex-meshes in these kinds of geometries tend to have both non-orthogonal and skew cells, which is the case also in the present mesh. Are you proposing that I should increase nNonOrthogonalCorrectors? I haven't yet compared the reults quantitatively with experiments and previous computations, but the results look fine in the upper part of the domain, which also has non-orthogonal and skew cells. The convergence problem is only below the hub.

I have tried both upwind and Gamma convection schemes with the same problem, but I actually don't think that the convergence problem is due to the convection schemes. What I do know is that upwind is likely to give bad results for this kind of flow.

I agree that there is probably a problem with the outlet boundary condition. What boundary conditions would you recommend at the outlet for this case with highly swirling flow? In previous research with another code I have used Neumann for all variables at the outlet.

I saw the newly discussed comparison between Fluent and OpenFOAM for the cyclinder in cross-flow and that there was also a problem with the outlet boundary condition for that case in OpenFOAM. It is really important to have good outlet boundary conditions unless you can afford to extend your computational domain, so if anyone has some good suggestion I would be very happy.

Håkan.

[eugene]

If you run checkmesh and your non-orthogonality is above say 70 degrees, you will want to put a limiter on the corrector of the laplacian schemes like:
laplacian(1|A(U),p) Gauss linear limited 0.3;

But as you say, this is unlikely to be your problem. As mentioned in the other thread, you can simulate an avergae pressure outlet boundary by updating the outlet each iteration in the top level code provided the outlet is on a single processor.

Basically something like:
//before solver starts
scalar pMeanOut = 0;

//check if this processor has an outlet boundary
label outletPatchID = mesh.boundary().patch().findPatchID("yourOutletBou ndaryNameHere")
if(outletPatchID != -1)
{
pMeanOut = mean(p.boundaryfield()[outletPatchID]);
}

//inside the time loop
if(outletPatchID != -1)
{
scalarField newP = p.boundaryField()[outletPatchID];

const labelList& outletPatchCells
= mesh.boundary()[outletPatchID]
.faceCells();

forAll(newP, faceI)
{
newP[faceI] = p[outletPatchCells[faceI]];
}
newP /= pMeanOut/mean(newP);
p.boundaryField()[outletPatchID]
== newP;
}

This hasnt been tested, so take it only as a suggestion. And before you try this just check if zero gradient on outlet pressure doesnt produce a stable solution.

[hani]

I will try the limiter on the corrector on the laplacian scheme in the future. I have many cells with non-orthogonality above 70 degrees. I think that this will have impact on the quality on the results, but to start with the issue is to get rid of the convergence problem, which I don't think will be dependent on this limiter.

About setting a mean pressure outlet boundary condition as you suggest, what is the difference numerically between that and a pure Neumann bc? As I understand it the only difference is that the mean pressure at the outlet is set to the specified reference pressure instead of specifying the reference pressure in a point when using pure Neumann. This would only result in a difference in the level of the pressure of the final result. I think that the way the reference pressure is set is mostly important when doing unsteady simulations, where it should be set in a way that it does not affect the overall pressure level, i.e. don't set the reference pressure in a point where the pressure is likely to oscillate. If you do your whole pressure field will start oscillating and you will probably have slow computations due to the imposed global fluctuating pressure level. In my case the computation is steady, so it should not matter how I specify the reference pressure (as long as I don't say that the pressure should be constant at the outlet). My conclusion is that a Neumann outlet boundary condition for the pressure (together with a reference pressure point) should be sufficient, and I'm testing that right now. The reason why I didn't test this already is that I had problems with this for another case. Maybe it is too unstable to use from scratch? Now I use it when restarting the solution I show above, and it is still converging after some 700 iterations.

I really appreciate your help Eugene!

Håkan.

[sshyu]

Håkan,

Normally I will look at the residual distribution and how it changes during the iterations. I guess that will give you more feeling about your problem.