[andrewburns]

I'm currently trying to do a simple steady state, laminar analysis of an imported 3D shape for a high Reynolds number.

The mesh is a simple tetrahedral mesh made with gmsh to the tune of around 1 million cells, the checkmesh results are as follows:

Checking geometry...
Domain bounding box: (-44.104 -16.605 -15.21) (184.676 16.605 18)
Boundary openness (5.39114e-17 -5.56627e-18 -1.00736e-17) OK.
Max cell openness = 1.75755e-16 OK.
Max aspect ratio = 590.305 OK.
Minumum face area = 4.02e-08. Maximum face area = 10.9795. Face area magnitudes OK.
Min volume = 8.9847e-10. Max volume = 12.1966. Total volume = 252267. Cell volumes OK.
Mesh non-orthogonality Max: 89.6816 average: 21.2003
*Number of severely non-orthogonal faces: 149.
Non-orthogonality check OK.
<<Writing 149 non-orthogonal faces to set nonOrthoFaces
Face pyramids OK.
Max skewness = 3.76367 OK.
*Edges too small, min/max edge length = 1e-05 5.86748, number too small: 1
<<Writing 2 points on short edges to set shortEdges
All angles in faces OK.
All face flatness OK.

Mesh OK.

End

------------------------------------------

Granted the mesh isn't exactly great but this is more due to the complexity of the object I'm trying to generate the flow around.

My scheme and solution setups are as follows:

-------------------------------------------

ddtSchemes
{
default steadyState;
}

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

divSchemes
{
default none;
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 Gauss linear corrected;
laplacian(nuEff,U) Gauss linear corrected;
laplacian((1|A(U)),p) Gauss linear limited 0.5;
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;
}


solvers
{
p GAMG
{
tolerance 1e-07;
relTol 0;
smoother GaussSeidel;
nCellsInCoarsestLevel 200;
mergeLevels 1;
agglomerator faceAreaPair;
cacheAgglomeration off;
nPreSweeps 1;
nPostSweeps 2;
nFinestSweeps 2;
scaleCorrection true;
directSolveCoarsest false;
};
U PBiCG
{
tolerance 1e-07;
relTol 0.1;
preconditioner DILU;
};
k PBiCG
{
tolerance 1e-07;
relTol 0.1;
preconditioner DILU;
};
epsilon PBiCG
{
tolerance 1e-07;
relTol 0.1;
preconditioner DILU;
};
R PBiCG
{
tolerance 1e-07;
relTol 0.1;
preconditioner DILU;
};
nuTilda PBiCG
{
tolerance 1e-07;
relTol 0.1;
preconditioner DILU;
};
}

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

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

-------------------------------------------

As things are currently, the solution diverges right from the start, after 2 iterations my lift and drag forces are in the order of 1e8 or so and my local continuity errors increase until eventually the system crashes, here's the readout from the first few iterations of the solution.

-------------------------------------------

Time = 1

DILUPBiCG: Solving for Ux, Initial residual = 1, Final residual = 0.0470618, No Iterations 1
DILUPBiCG: Solving for Uy, Initial residual = 1, Final residual = 0.0263094, No Iterations 1
DILUPBiCG: Solving for Uz, Initial residual = 1, Final residual = 0.0251062, No Iterations 1
GAMG: Solving for p, Initial residual = 1, Final residual = 8.6099e-08, No Iterations 18
GAMG: Solving for p, Initial residual = 0.0975858, Final residual = 8.43973e-08, No Iterations 16
GAMG: Solving for p, Initial residual = 0.0248457, Final residual = 9.02797e-08, No Iterations 13
GAMG: Solving for p, Initial residual = 0.00976953, Final residual = 5.16382e-08, No Iterations 13
GAMG: Solving for p, Initial residual = 0.00536754, Final residual = 6.71377e-08, No Iterations 12
GAMG: Solving for p, Initial residual = 0.00330867, Final residual = 7.83699e-08, No Iterations 11
time step continuity errors : sum local = 0.013064, global = 3.73088e-05, cumulative = 3.73088e-05

Total pressure Force = (-0.534847 -2.01937 7.96596)
Total viscous Force = (-0.00165139 -0.00937387 0.0347963)
Total turbulent Force = (0 0 0)
Total Force = (-0.536499 -2.02875 8.00076)
ExecutionTime = 116.99 s ClockTime = 117 s

Time = 2

DILUPBiCG: Solving for Ux, Initial residual = 0.33876, Final residual = 0.0240347, No Iterations 1
DILUPBiCG: Solving for Uy, Initial residual = 0.419772, Final residual = 0.0253754, No Iterations 1
DILUPBiCG: Solving for Uz, Initial residual = 0.41666, Final residual = 0.0234304, No Iterations 1
GAMG: Solving for p, Initial residual = 1, Final residual = 9.54611e-08, No Iterations 25
GAMG: Solving for p, Initial residual = 0.0179724, Final residual = 6.66377e-08, No Iterations 14
GAMG: Solving for p, Initial residual = 0.00555565, Final residual = 6.90765e-08, No Iterations 14
GAMG: Solving for p, Initial residual = 0.00274587, Final residual = 6.89308e-08, No Iterations 12
GAMG: Solving for p, Initial residual = 0.00169177, Final residual = 7.81096e-08, No Iterations 12
GAMG: Solving for p, Initial residual = 0.00125716, Final residual = 7.1114e-08, No Iterations 12
time step continuity errors : sum local = 0.102613, global = -7.33632e-05, cumulative = -3.60544e-05

Total pressure Force = (-5.78598e+06 -8.6209e+06 2.9444e+07)
Total viscous Force = (-23.9807 -11.6406 33.9165)
Total turbulent Force = (0 0 0)
Total Force = (-5.78601e+06 -8.62091e+06 2.94441e+07)
ExecutionTime = 233.26 s ClockTime = 234 s

Time = 3

DILUPBiCG: Solving for Ux, Initial residual = 0.278798, Final residual = 0.00263389, No Iterations 2
DILUPBiCG: Solving for Uy, Initial residual = 0.274962, Final residual = 0.0194582, No Iterations 1
DILUPBiCG: Solving for Uz, Initial residual = 0.281969, Final residual = 0.0142537, No Iterations 1
GAMG: Solving for p, Initial residual = 0.149281, Final residual = 8.71556e-08, No Iterations 23
GAMG: Solving for p, Initial residual = 0.00644011, Final residual = 8.54485e-08, No Iterations 11
GAMG: Solving for p, Initial residual = 0.00191251, Final residual = 7.56476e-08, No Iterations 12
GAMG: Solving for p, Initial residual = 0.000863805, Final residual = 6.18423e-08, No Iterations 11
GAMG: Solving for p, Initial residual = 0.00051502, Final residual = 7.79148e-08, No Iterations 11
GAMG: Solving for p, Initial residual = 0.000364454, Final residual = 7.7788e-08, No Iterations 11
time step continuity errors : sum local = 0.245931, global = 6.00056e-06, cumulative = -3.00538e-05

Total pressure Force = (-3.12704e+07 1.20367e+08 4.65635e+08)
Total viscous Force = (-28.0773 -10.4569 33.8212)
Total turbulent Force = (0 0 0)
Total Force = (-3.12704e+07 1.20367e+08 4.65635e+08)
ExecutionTime = 340.78 s ClockTime = 341 s

Time = 4

DILUPBiCG: Solving for Ux, Initial residual = 0.256524, Final residual = 0.00287552, No Iterations 2
DILUPBiCG: Solving for Uy, Initial residual = 0.212373, Final residual = 0.0138832, No Iterations 1
DILUPBiCG: Solving for Uz, Initial residual = 0.209775, Final residual = 0.0157668, No Iterations 1
GAMG: Solving for p, Initial residual = 0.0282568, Final residual = 9.9138e-08, No Iterations 33
GAMG: Solving for p, Initial residual = 0.00410752, Final residual = 8.15701e-08, No Iterations 17
GAMG: Solving for p, Initial residual = 0.00116876, Final residual = 7.89225e-08, No Iterations 13
GAMG: Solving for p, Initial residual = 0.000591451, Final residual = 8.9654e-08, No Iterations 9
GAMG: Solving for p, Initial residual = 0.000380289, Final residual = 7.98343e-08, No Iterations 9
GAMG: Solving for p, Initial residual = 0.000275326, Final residual = 9.12885e-08, No Iterations 8
time step continuity errors : sum local = 0.286709, global = -0.000247643, cumulative = -0.000277697

Total pressure Force = (-1.14291e+08 4.25839e+07 7.23494e+08)
Total viscous Force = (-31.085 -10.3115 52.0479)
Total turbulent Force = (0 0 0)
Total Force = (-1.14291e+08 4.25839e+07 7.23494e+08)
ExecutionTime = 459.04 s ClockTime = 460 s


-------------------------------------------

It's my understanding that OpenFoam is more sensitive to bad meshes than commercial solvers, however this is probably the best mesh I'm going to be able to generate with gmsh, so is there any way I can increase the stability of my solution?

[andrewburns]

So nobody has any ideas?

[msha]

FOAM is senssitive to mesh as much as others or even less.
instead you must focus on cfd method to stabilize solution by means of relaxation or so.

[andrewburns]

That's all well and good but I've run down to a pressure relaxation factor of 0.01 and 10 orthagonal correctors and still my solution explodes to infinity in a few iterations. Clearly something is wrong because I can get other meshes to work just fine, the only difference is this time gmsh can't really cope with the geometry and so I have a fairly poor mesh.

[hjasak]

Sorry for being stupid, but where is your turbulence model being updated? There's no solver calls for k-epsilon, RSM or similar...

If it's not that, you probably messed up the boundary conditions. Your mesh is criminally bad, but the solver should still run.

Which turbulence model are you using? Did you try running potentialFoam and looking at the velocity field? This will give you an idea if the mesh is usable and give you a quick check on the boundary conditions.

Hrv

[andrewburns]

The turbulence is set to laminar and switched off in simplefoam. I ran a potential foam solution and it too could not reach a stable solution however I don't think my boundary conditions are incorrect, I made sure to open the mesh in paraview and check that all the patches are what I want them to be.

[hjasak]

If you cannot get a stable solution with potentialFoam, your mesh is simply not good enough. What does the velocity feel look like? Is it horrid + where?

[andrewburns]

Yeah I don't think gmsh is going to let me mesh this thing. If you catch a look at the velocity field before it blows up it's horribly splotchy, after the solution blows up it doesn't make any sense anyway.

Thanks for the help anyway, I guess I'll have to look into other ways of meshing.

[mighelone]

Why are you using laminar flow with high Re number flow? Probably your flow is not laminar, but turbulent, try to use a turbulence model in order to dump the instabilities of the flow.