[sivakumar]

Hi There,
I am looking for your suggestion and help to improve the accuracy of my simulation results.

I am simulating axial flow fan using MRFSimpleFoam (OF-2.1.1),
I have obtained some results with 22% error. Now I want to improve the results, please give me your suggestions and advise.

Case setup:

Mesh :5.25Million Cells,(hex+tet)
volumeFlowRate :16m3/s
turbulence model :k-Epsilon
y+ :50 -70

CheckMesh output:

Create time

Create polyMesh for time = 0

Time = 0

Mesh stats
points: 4740987
faces: 15102544
internal faces: 14812112
cells: 5235464
boundary patches: 16
point zones: 0
face zones: 1
cell zones: 2

Overall number of cells of each type:
hexahedra: 4480000
prisms: 0
wedges: 0
pyramids: 12800
tet wedges: 0
tetrahedra: 742664
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
outlet 6400 6561 ok (non-closed singly connected)
inlet 6400 6561 ok (non-closed singly connected)
top2 44000 44631 ok (non-closed singly connected)
top1 11276 5850 ok (non-closed singly connected)
top0 12000 12231 ok (non-closed singly connected)
center1 44000 44631 ok (non-closed singly connected)
fan 30520 15608 ok (non-closed singly connected)
center0 12000 12231 ok (non-closed singly connected)
ILR0 12000 12231 ok (non-closed singly connected)
ILR1 12000 12231 ok (non-closed singly connected)
OLR0 44000 44631 ok (non-closed singly connected)
OLR1 44000 44631 ok (non-closed singly connected)
CLR0 3308 1766 ok (non-closed singly connected)
CLR1 3308 1766 ok (non-closed singly connected)
FCLR0 2610 1411 ok (non-closed singly connected)
FCLR1 2610 1411 ok (non-closed singly connected)

Checking geometry...
Overall domain bounding box (-0.0816836 0.156811 -1.5) (0.59502 0.771 6)
Mesh (non-empty, non-wedge) directions (1 1 1)
Mesh (non-empty) directions (1 1 1)
Boundary openness (-5.07799e-16 7.24744e-15 -8.38482e-17) OK.
Max cell openness = 3.31697e-16 OK.
Max aspect ratio = 12.0865 OK.
Minumum face area = 1.42058e-06. Maximum face area = 0.000219702. Face area magnitudes OK.
Min volume = 7.33795e-10. Max volume = 1.27874e-06. Total volume = 1.4718. Cell volumes OK.
Mesh non-orthogonality Max: 65.6246 average: 9.54892
Non-orthogonality check OK.
Face pyramids OK.
Max skewness = 0.920449 OK.
Coupled point location match (average 2.51368e-09) OK.

Mesh OK.

End

fvSchemes:

ddtSchemes
{
default steadyState;
}

gradSchemes
{
default Gauss linear;
grad(p) leastSquares;
grad(omega) Gauss linear;
grad(epsilon) Gauss linear;
grad(U) Gauss linear;
grad(k) Gauss linear;
}

divSchemes
{
default none;
div(phi,U) Gauss upwind;
div(phi,k) Gauss upwind;
div(phi,omega) Gauss upwind;
div(phi,epsilon) Gauss upwind ;
div((nuEff*dev(T(grad(U))))) Gauss linear;
}

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

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

snGradSchemes
{
default limited 0.5;
}

fluxRequired
{
default no;
p ;
}


Due to the non-orthogonality issue I am not able to use full corrected schemes. So I am using some limitations for that.
However the above schemes are under predicts the torque and pressure raise.
obtained torque -67.2 N-m,
expected is 72 N-m,

obtained pressure raise -177Pa
expected -227Pa

When I use full corrected scheme, I am facing convergence problem. Epsilon always bounding and fluctuating for fully corrected schemes.
full corrected schemes over predicts the pressure and torque.

I have attached the convergence plot for both cases, please have a look.

fvSolution

solvers
{
p
{
solver GAMG;
tolerance 1e-08;
relTol 0.01;
smoother GaussSeidel;
nPreSweeps 0;
nPostSweeps 2;
cacheAgglomeration true;
nCellsInCoarsestLevel 20;
agglomerator faceAreaPair;
mergeLevels 1;
}

U
{
solver smoothSolver;
smoother GaussSeidel;
nSweeps 2;
tolerance 1e-07;
relTol 0.01;
minIter 1;
maxIter 1000;
}

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

omega
{
solver smoothSolver;
smoother GaussSeidel;
tolerance 1e-7;
relTol 0.01;
nSweeps 1;
minIter 1;
maxIter 1000;
}

epsilon
{
solver smoothSolver;
smoother GaussSeidel;
nSweeps 2;
tolerance 1e-07;
relTol 0.005;
minIter 1;
maxIter 1000;
}
}

SIMPLE
{
momentumPredictor yes;
nNonOrthogonalCorrectors 1;
pRefCell 0;
pRefValue 0;
convergence 1e-08;

residualControl
{
p 1e-6;
U 1e-6;
nuTilda 1e-6;
}
}

potentialFlow
{
nNonOrthogonalCorrectors 10;
}

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


I have attached the 0 folder for your check, please have a look and give me some suggestions to improve my results.

I tried with different mesh sizes but there no much difference.

Thanks,
Sivakumar

[immortality]

why don't you use higher order schemes in return of Upwind like limitedLinear or linearUpwind,... to improve accuracy?does it go unstable?

[sivakumar]

Hi,
Thanks for your reply, I didn’t try the higher order schemes, for stability purpose I am using upwind.

If its possible can you please tell me which higher order scheme will be nice?

Thanks,
Sivakumar

[immortality]

test vanLeer for div(phi,...( all except U)) and vanLeerV for div(phi,U)
in one run and in another gamma and gammaV respectively(check for its correct capitals,i don't remember) and limitedLinear 1; and limitedLinearV 1;
if they are stable or not.

[sivakumar]

Hi,
Now I changed my fvSchemes as follows,

ddtSchemes
{
default steadyState;
}

gradSchemes
{
default leastSquares;
}

divSchemes
{
default none;
div(phi,U) Gauss vanLeerV;
div(phi,k) Gauss vanLeer;
div(phi,omega) Gauss vanLeer;
div(phi,epsilon) Gauss vanLeer;
div((nuEff*dev(T(grad(U))))) Gauss linear;
}

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

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

snGradSchemes
{
default limited 0.5;
}

fluxRequired
{
default no;
p ;
}

epsilon and k both started bounding

Time = 1

smoothSolver: Solving for Ux, Initial residual = 1, Final residual = 0.00119926, No Iterations 4
smoothSolver: Solving for Uy, Initial residual = 1, Final residual = 0.000995463, No Iterations 4
smoothSolver: Solving for Uz, Initial residual = 1, Final residual = 0.00217619, No Iterations 4
GAMG: Solving for p, Initial residual = 1, Final residual = 0.00988547, No Iterations 176
GAMG: Solving for p, Initial residual = 0.0371835, Final residual = 0.000146369, No Iterations 3
time step continuity errors : sum local = 0.164183, global = -0.000512996, cumulative = -0.000512996
smoothSolver: Solving for epsilon, Initial residual = 0.0464226, Final residual = 0.000107885, No Iterations 4
bounding epsilon, min: -22.5128 max: 335.166 average: 2.43186
smoothSolver: Solving for k, Initial residual = 1, Final residual = 0.00464138, No Iterations 4
bounding k, min: -15.3328 max: 108.947 average: 0.465543
ExecutionTime = 92.65 s ClockTime = 100 s

Time = 2

smoothSolver: Solving for Ux, Initial residual = 0.353353, Final residual = 0.00311845, No Iterations 4
smoothSolver: Solving for Uy, Initial residual = 0.365416, Final residual = 0.0034596, No Iterations 4
smoothSolver: Solving for Uz, Initial residual = 0.875299, Final residual = 0.00135475, No Iterations 6
GAMG: Solving for p, Initial residual = 0.0259538, Final residual = 0.000255184, No Iterations 38
GAMG: Solving for p, Initial residual = 0.0297287, Final residual = 0.000223586, No Iterations 4
time step continuity errors : sum local = 0.113884, global = -0.0134757, cumulative = -0.0139887
smoothSolver: Solving for epsilon, Initial residual = 0.133081, Final residual = 0.000566866, No Iterations 4
bounding epsilon, min: -1001.87 max: 130408 average: 17.0609
smoothSolver: Solving for k, Initial residual = 0.384752, Final residual = 0.00320007, No Iterations 4
bounding k, min: -66.4828 max: 700.014 average: 0.775518
ExecutionTime = 120.93 s ClockTime = 130 s

Time = 3

smoothSolver: Solving for Ux, Initial residual = 0.703536, Final residual = 0.00279422, No Iterations 4
smoothSolver: Solving for Uy, Initial residual = 0.814474, Final residual = 0.00330137, No Iterations 4
smoothSolver: Solving for Uz, Initial residual = 0.672708, Final residual = 0.00104579, No Iterations 6
GAMG: Solving for p, Initial residual = 0.0775572, Final residual = 0.000750949, No Iterations 80
GAMG: Solving for p, Initial residual = 0.0451604, Final residual = 0.000448057, No Iterations 11
time step continuity errors : sum local = 0.121181, global = 0.016846, cumulative = 0.00285731
smoothSolver: Solving for epsilon, Initial residual = 0.201837, Final residual = 0.000159865, No Iterations 6
bounding epsilon, min: -37133.4 max: 1.16809e+06 average: 77.2635
smoothSolver: Solving for k, Initial residual = 0.428879, Final residual = 0.000492869, No Iterations 6
bounding k, min: -2619.3 max: 26710.7 average: 1.564
ExecutionTime = 170.6 s ClockTime = 181 s

it seems the calculation will be stable only for first order schemes (upwind).

any Idea?

Thanks,
Sivakumar

[immortality]

what do you mean by full corrected schemes?did you use unCorrected scheme in laplacianScheme so far?
i thought when there is non-orthogonality in mesh its better to use corrected scheme to calculate the gradients normal to faces.whats the correct usage of it?

[sivakumar]

Hi,
Till now I am using limited 0.5 for LaplacianScheme, because my mesh has orthogonal angle of 65, so I need to use it, otherwise the calculation will diverge.
I don't know, is there any other way to get rid of this issue?

Just have a look of my first post, if you know any other way I am ready to use it.

Thanks,
Sivakumar

[immortality]

i don't understand why should corrected scheme lead to diverge when the mesh has non orthogonality?could you tell me?
thanks.

[sivakumar]

I have seen in the training materials, when your non-orthogonal angle is

below 50 you can use full corrected scheme,
if it has angle of 50 to 70 you need to use limited correction
70 to 80 stability possible accuracy compromised
above 80 stability very difficult to attain.

based on this information I am using limited 0.5, I tried 0.7 as well but I am not able to get the stability. I feel still there are some tricky things to get good accuracy, that is the information I am looking from you people.

Thanks,
Sivakumar

[immortality]

then when uncorrected is applicable?when mesh is orthogonal.correct?