[dsemikin]

Hello.


I am trying to solve simple sample external aerodynamics problem. The body is a block with one face slightly inclined and all edges and corners rounded. Sizes: 1 x 0.5 x 0.2 m


Air velocity: 80 m/s.


From this post: How to resolve boundary layer in OF ?


I concluded, that to resolve boundary layer I just need to create more dense layers, which I tried to do.



When I use mesh with first mesh layer height ~ Y+ = 30 (~ 3.0E-4 m) or more, then solution runs smoothly.


But when I try to use mesh with 30 layers, where the first one has height of ~ Y+ = 1 (~ 1.0e-5 m), then the same simulation always diverges.


If I look at the picture of the latest iterations before divergence, then problematic faces are located in quite random (from my point of view) place often on the big flat surface. Surface mesh size at this location may be ~0.01m, so the aspect ratio of the layers is not good, but I was sure, that it should not be an issue, because there should be no velocity component normal to the wall...


The stacktrace in the log contains the following lines, so I assume, that it is where the floating point exception happens:

[13] #4  Foam::GaussSeidelSmoother::smooth(Foam::Field<double>&, Foam::Field<double> const&, unsigned char, int) const at ??:?
[6] #5  Foam::smoothSolver::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const at ??:?

Below are also the log of checkMesh and latest iterations before divergence.



In checkMesh log I don't see anything suspicious.



And in simpleFoam log one can see, that abruptly the step continuity error explodes and then floating point exception happens.


I would highly appreciate, if someone could explain, what could be the reason for divergence. In my understanding, usually it is poor quality of mesh, but in this case I don't see, how it is the case, and what property of mesh is "bad"...


Of course, also is of interest, how to overcome this. I tried to decrease relaxation factors and tolerances, but it did not help. Probably some other linear solver could giver more robust convergence? Or some other numerical settings?



Thank you in advance.


checkMesh log:

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  8
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
Build  : 8-9b73cf21a682
Exec   : checkMesh
Date   : Oct 06 2020
Time   : 11:29:53
Host   : "openfoam-vm"
PID    : 4516
I/O    : uncollated
Case   : /data/super-ka-pjt/super-ka/2020-10-05_block-01_06_resolve_boundary2
nProcs : 1
sigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).
fileModificationChecking : Monitoring run-time modified files using timeStampMaster (fileModificationSkew 10)
allowSystemOperations : Allowing user-supplied system call operations

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

Create polyMesh for time = 0

Time = 0

Mesh stats
    points:           3798462
    faces:            22123546
    internal faces:   21889090
    cells:            9373751
    faces per cell:   4.695306713
    boundary patches: 7
    point zones:      0
    face zones:       0
    cell zones:       1

Overall number of cells of each type:
    hexahedra:     0
    prisms:        6517632
    wedges:        0
    pyramids:      0
    tet wedges:    0
    tetrahedra:    2856119
    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                  
    air_y_mid           29554    23007    ok (non-closed singly connected)  
    air_side            269      160      ok (non-closed singly connected)  
    air_bottom          380      223      ok (non-closed singly connected)  
    air_top             330      194      ok (non-closed singly connected)  
    air_inlet           177      107      ok (non-closed singly connected)  
    air_outlet          70       48       ok (non-closed singly connected)  
    car_surf            203676   102087   ok (non-closed singly connected)  

Checking geometry...
    Overall domain bounding box (-2.964644661 -4.587996649e-14 -3) (10 3.25 3.2)
    Mesh has 3 geometric (non-empty/wedge) directions (1 1 1)
    Mesh has 3 solution (non-empty) directions (1 1 1)
    Boundary openness (-2.18097724e-17 7.340355336e-16 -3.994958304e-18) OK.
    Max cell openness = 1.517392784e-14 OK.
    Max aspect ratio = 585.8712006 OK.
    Minimum face area = 6.316036633e-09. Maximum face area = 0.6207478254.  Face area magnitudes OK.
    Min volume = 1.95161539e-12. Max volume = 0.1435894578.  Total volume = 261.1931449.  Cell volumes OK.
    Mesh non-orthogonality Max: 55.79104786 average: 9.668060231
    Non-orthogonality check OK.
    Face pyramids OK.
    Max skewness = 0.5989422464 OK.
    Coupled point location match (average 0) OK.

Mesh OK.

End

Latest iterations of simpleFoam:

...


Time = 53

smoothSolver:  Solving for Ux, Initial residual = 0.4824992416, Final residual = 0.0427002703, No Iterations 1
smoothSolver:  Solving for Uy, Initial residual = 0.5263585749, Final residual = 0.01438578753, No Iterations 2
smoothSolver:  Solving for Uz, Initial residual = 0.4369752154, Final residual = 0.0368871982, No Iterations 2
GAMG:  Solving for p, Initial residual = 0.06624748546, Final residual = 0.003156955786, No Iterations 5
time step continuity errors : sum local = 1.56182585, global = -0.00110296805, cumulative = -0.1549218169
smoothSolver:  Solving for omega, Initial residual = 0.08453582241, Final residual = 0.006250842107, No Iterations 1
smoothSolver:  Solving for k, Initial residual = 0.2660798304, Final residual = 0.01876795111, No Iterations 1
ExecutionTime = 245.59 s  ClockTime = 247 s

Time = 54

smoothSolver:  Solving for Ux, Initial residual = 0.3337720154, Final residual = 0.01002013169, No Iterations 5
smoothSolver:  Solving for Uy, Initial residual = 0.3340450882, Final residual = 0.01110384285, No Iterations 4
smoothSolver:  Solving for Uz, Initial residual = 0.3731827088, Final residual = 0.009254611244, No Iterations 6
GAMG:  Solving for p, Initial residual = 0.8799070613, Final residual = 0.03313636654, No Iterations 2
time step continuity errors : sum local = 893.6400296, global = 0.03382824822, cumulative = -0.1210935687
smoothSolver:  Solving for omega, Initial residual = 0.318287741, Final residual = 0.002645842889, No Iterations 1
smoothSolver:  Solving for k, Initial residual = 0.1526440356, Final residual = 0.00615373916, No Iterations 1
ExecutionTime = 250.03 s  ClockTime = 252 s

Time = 55

smoothSolver:  Solving for Ux, Initial residual = 0.7455203745, Final residual = 0.009375546185, No Iterations 2
smoothSolver:  Solving for Uy, Initial residual = 0.7133473046, Final residual = 0.005211766964, No Iterations 2
smoothSolver:  Solving for Uz, Initial residual = 0.5544917958, Final residual = 0.0123354264, No Iterations 2
GAMG:  Solving for p, Initial residual = 0.3172217772, Final residual = 0.0135137624, No Iterations 4
time step continuity errors : sum local = 3838.404266, global = 3.380422244, cumulative = 3.259328675
smoothSolver:  Solving for omega, Initial residual = 0.1870470368, Final residual = 0.01061101629, No Iterations 2
smoothSolver:  Solving for k, Initial residual = 0.1152355589, Final residual = 0.008423114327, No Iterations 2
ExecutionTime = 254.61 s  ClockTime = 256 s

Time = 56

smoothSolver:  Solving for Ux, Initial residual = 0.445060718, Final residual = 0.02737065041, No Iterations 16
smoothSolver:  Solving for Uy, Initial residual = 0.5242167792, Final residual = 0.02489138409, No Iterations 16
smoothSolver:  Solving for Uz, Initial residual = 0.4304852584, Final residual = 0.02229787888, No Iterations 16
GAMG:  Solving for p, Initial residual = 0.9922608279, Final residual = 0.03727616882, No Iterations 33
time step continuity errors : sum local = 4.971918037e+10, global = -76473931.53, cumulative = -76473928.27
smoothSolver:  Solving for omega, Initial residual = 0.2511483143, Final residual = 0.007224566153, No Iterations 2
smoothSolver:  Solving for k, Initial residual = 0.08664252862, Final residual = 0.006499034169, No Iterations 3
ExecutionTime = 266.15 s  ClockTime = 268 s

regIOobject::readIfModified() :
    Re-reading object controlDict from file "/data/super-ka-pjt/super-ka/2020-10-05_block-01_06_resolve_boundary2/system/controlDict"
Time = 57

smoothSolver:  Solving for Ux, Initial residual = 0.3679538748, Final residual = 0.03369544687, No Iterations 25
smoothSolver:  Solving for Uy, Initial residual = 0.5404464887, Final residual = 0.04559830362, No Iterations 22
smoothSolver:  Solving for Uz, Initial residual = 0.3205770382, Final residual = 0.02557449322, No Iterations 29
GAMG:  Solving for p, Initial residual = 0.9999999995, Final residual = 1.102454466, No Iterations 1000
time step continuity errors : sum local = 3.032927413e+23, global = 1.375129909e+21, cumulative = 1.375129909e+21
smoothSolver:  Solving for omega, Initial residual = 0.3376810441, Final residual = 0.01272912993, No Iterations 2
smoothSolver:  Solving for k, Initial residual = 0.04744437071, Final residual = 0.003486061951, No Iterations 2
ExecutionTime = 433.7 s  ClockTime = 436 s



... (stack trace)

[dsemikin]

Minor update:


Now I also tried to map solution from coarser mesh (Y+ = 30) to the finer mesh. And then run the simulation with the fine mesh. But the result diverged on a first iteration...


Still no luck .

[Tobermory]

I haven't come across this problem before, but in response to your question: "which aspect of the mesh is bad?", one answer could be that your aspect ratio is 585, which is pretty high (I know that checkMesh says this is "okay", but it is still pretty high). Alternatives could be that there is a large jump in mesh size between your prism layers and the background core mesh (I can't tell this from just the checkMesh output .... maybe you could upload a picture?).

For fixes, try looking at the gradient and diffusion schemes that you have used in fvSchemes (perhaps upload your fvSchemes file?). For hints on how to make this more stable, take a look at some of the online content by e.g. Guerrero
http://www.wolfdynamics.com/training...ession2020.pdf

or there's a youtube video of a lecture by Prof Jasak where he gives some advice on solver settings:
https://www.youtube.com/watch?v=a4B_oXR5Kzs

Good luck, and let us know how you solved it in the end!

[dsemikin]

Hello Tobermory,


Thank you for the links. I will check it out.


Meanwhile, here is my fvSChemes and screenshots of the body and mesh attached.


with kind regards
Dmitrii.