[6863523]

Dear,
I have been facing the problem below for two geometries.

/*---------------------------------------------------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  4.1                                   |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
Build  : 4.1
Exec   : myIcoFoam -parallel
Date   : Mar 09 2017
Time   : 14:06:12
Host   : "avakas-frontend2"
PID    : 1938
Case   : /home/ywang005/run/GWDf1/test/Re15000
nProcs : 12
Slaves : 
11
(
"avakas-frontend2.1939"
"avakas-frontend2.1940"
"avakas-frontend2.1941"
"avakas-frontend2.1942"
"avakas-frontend2.1943"
"avakas-frontend2.1944"
"avakas-frontend2.1945"
"avakas-frontend2.1946"
"avakas-frontend2.1947"
"avakas-frontend2.1948"
"avakas-frontend2.1949"
)

Pstream initialized with:
    floatTransfer      : 0
    nProcsSimpleSum    : 0
    commsType          : nonBlocking
    polling iterations : 0
sigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).
fileModificationChecking : Monitoring run-time modified files using timeStampMaster
allowSystemOperations : Allowing user-supplied system call operations

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


Create mesh for time = 0


PISO: Operating solver in PISO mode

Reading transportProperties

Reading field p

Reading field U

Reading/calculating face flux field phi


Starting time loop

volRegion volumeAverage:all():
    total cells  = 8396
    total volume = 20249.7



Time = 0.0125

Courant Number mean: 0 max: 0
DILUPBiCG:  Solving for Ux, Initial residual = 1, Final residual = 3.14067e-11, No Iterations 18
DILUPBiCG:  Solving for Uy, Initial residual = 0, Final residual = 0, No Iterations 0
DICPCG:  Solving for p, Initial residual = 1, Final residual = 9.96628e-08, No Iterations 246
time step continuity errors : sum local = 2.62573e-13, global = 2.24611e-24, cumulative = 2.24611e-24
DICPCG:  Solving for p, Initial residual = 0.00826593, Final residual = 9.75861e-08, No Iterations 216
time step continuity errors : sum local = 5.1487e-13, global = 4.59247e-24, cumulative = 6.83858e-24
ExecutionTime = 0.26 s  ClockTime = 0 s

Time = 0.025

Courant Number mean: 5.34236e-06 max: 0.00488047
DILUPBiCG:  Solving for Ux, Initial residual = 0.754291, Final residual = 8.55034e-11, No Iterations 16
DILUPBiCG:  Solving for Uy, Initial residual = 0.330802, Final residual = 6.98711e-11, No Iterations 16
DICPCG:  Solving for p, Initial residual = 0.420872, Final residual = 9.52255e-08, No Iterations 265
time step continuity errors : sum local = 5.26749e-13, global = -1.48591e-23, cumulative = -8.02052e-24
DICPCG:  Solving for p, Initial residual = 0.0611023, Final residual = 9.76153e-08, No Iterations 261
time step continuity errors : sum local = 7.56513e-13, global = -8.53319e-24, cumulative = -1.65537e-23
ExecutionTime = 0.29 s  ClockTime = 0 s

Time = 0.0375

Courant Number mean: 1.02253e-05 max: 0.0113322
DILUPBiCG:  Solving for Ux, Initial residual = 0.528265, Final residual = 5.87478e-11, No Iterations 17
DILUPBiCG:  Solving for Uy, Initial residual = 0.266124, Final residual = 9.26735e-11, No Iterations 17
DICPCG:  Solving for p, Initial residual = 0.538462, Final residual = 8.96817e-08, No Iterations 287
time step continuity errors : sum local = 1.26773e-12, global = 6.13145e-22, cumulative = 5.96591e-22
DICPCG:  Solving for p, Initial residual = 0.579184, Final residual = 9.6864e-08, No Iterations 280
time step continuity errors : sum local = 3.73869e-12, global = -6.34438e-22, cumulative = -3.78469e-23
ExecutionTime = 0.33 s  ClockTime = 0 s

Time = 0.05

Courant Number mean: 4.4832e-05 max: 0.370209
DILUPBiCG:  Solving for Ux, Initial residual = 0.390239, Final residual = 6.12355e-11, No Iterations 18
DILUPBiCG:  Solving for Uy, Initial residual = 0.312606, Final residual = 1.82604e-11, No Iterations 18
DICPCG:  Solving for p, Initial residual = 0.916305, Final residual = 9.79778e-08, No Iterations 292
time step continuity errors : sum local = 3.3445e-11, global = -1.1319e-21, cumulative = -1.16975e-21
DICPCG:  Solving for p, Initial residual = 0.7854, Final residual = 9.34491e-08, No Iterations 282
time step continuity errors : sum local = 1.78264e-10, global = -2.58022e-20, cumulative = -2.6972e-20
ExecutionTime = 0.36 s  ClockTime = 0 s

Time = 0.0625

Courant Number mean: 0.00188543 max: 15.577
DILUPBiCG:  Solving for Ux, Initial residual = 0.344186, Final residual = 2.53233e-11, No Iterations 19
DILUPBiCG:  Solving for Uy, Initial residual = 0.325652, Final residual = 1.0124e-11, No Iterations 19
DICPCG:  Solving for p, Initial residual = 0.887083, Final residual = 9.42448e-08, No Iterations 295
time step continuity errors : sum local = 1.40168e-09, global = -4.72842e-18, cumulative = -4.75539e-18
DICPCG:  Solving for p, Initial residual = 0.968391, Final residual = 9.63933e-08, No Iterations 290
time step continuity errors : sum local = 1.7746e-08, global = 1.16238e-17, cumulative = 6.86843e-18
ExecutionTime = 0.39 s  ClockTime = 0 s

Time = 0.075

Courant Number mean: 0.535059 max: 1724.54
DILUPBiCG:  Solving for Ux, Initial residual = 0.642258, Final residual = 0.0046226, No Iterations 1001
DILUPBiCG:  Solving for Uy, Initial residual = 0.966405, Final residual = 0.00529744, No Iterations 1001
DICPCG:  Solving for p, Initial residual = 0.892676, Final residual = 0.036825, No Iterations 1001
time step continuity errors : sum local = 0.0900662, global = -9.51614e-17, cumulative = -8.82929e-17
DICPCG:  Solving for p, Initial residual = 0.979527, Final residual = 0.0190689, No Iterations 1001
time step continuity errors : sum local = 2.48803, global = 7.10808e-14, cumulative = 7.09925e-14
ExecutionTime = 0.66 s  ClockTime = 0 s

Time = 0.0875

Courant Number mean: 858.698 max: 2.36043e+07
DILUPBiCG:  Solving for Ux, Initial residual = 0.999976, Final residual = 1.00727, No Iterations 1001
DILUPBiCG:  Solving for Uy, Initial residual = 0.999974, Final residual = 3.69304, No Iterations 1001
DICPCG:  Solving for p, Initial residual = 0.999999, Final residual = 0.119223, No Iterations 1001
time step continuity errors : sum local = 2.24302e+07, global = -9.7959e-10, cumulative = -9.79519e-10
DICPCG:  Solving for p, Initial residual = 0.981905, Final residual = 0.335337, No Iterations 1001
time step continuity errors : sum local = 3.96954e+09, global = 1.92904e-07, cumulative = 1.91924e-07
ExecutionTime = 0.91 s  ClockTime = 1 s

Time = 0.1

Courant Number mean: 8.4947e+09 max: 1.06993e+14
DILUPBiCG:  Solving for Ux, Initial residual = 0.997632, Final residual = 51.6564, No Iterations 1001
DILUPBiCG:  Solving for Uy, Initial residual = 0.997067, Final residual = 1.91243, No Iterations 1001
DICPCG:  Solving for p, Initial residual = 1, Final residual = 12.3467, No Iterations 1001
time step continuity errors : sum local = 3.52888e+15, global = -0.0227559, cumulative = -0.0227557
DICPCG:  Solving for p, Initial residual = 0.977625, Final residual = 6.17759, No Iterations 1001
time step continuity errors : sum local = 1.0004e+18, global = 41.4258, cumulative = 41.403
ExecutionTime = 1.17 s  ClockTime = 1 s

Time = 0.1125

Courant Number mean: 2.11207e+18 max: 3.46701e+22
DILUPBiCG:  Solving for Ux, Initial residual = 0.997125, Final residual = 0.282022, No Iterations 1001
DILUPBiCG:  Solving for Uy, Initial residual = 0.995716, Final residual = 0.775586, No Iterations 1001
DICPCG:  Solving for p, Initial residual = 1, Final residual = 76.605, No Iterations 1001
time step continuity errors : sum local = 1.66192e+23, global = 9.79902e+06, cumulative = 9.79906e+06
DICPCG:  Solving for p, Initial residual = 0.980146, Final residual = 1.73321, No Iterations 1001
time step continuity errors : sum local = 1.30072e+25, global = -4.72389e+08, cumulative = -4.6259e+08
ExecutionTime = 1.44 s  ClockTime = 1 s

Time = 0.125

Courant Number mean: 2.54375e+25 max: 5.71293e+29
DILUPBiCG:  Solving for Ux, Initial residual = 0.999787, Final residual = 2.98944, No Iterations 1001
DILUPBiCG:  Solving for Uy, Initial residual = 0.995576, Final residual = 1.0834, No Iterations 1001
DICPCG:  Solving for p, Initial residual = 1, Final residual = 24.7115, No Iterations 1001
time step continuity errors : sum local = 1.00218e+30, global = 8.54027e+12, cumulative = 8.53981e+12
DICPCG:  Solving for p, Initial residual = 0.96636, Final residual = 4.26514, No Iterations 1001
time step continuity errors : sum local = 1.13572e+32, global = -8.38086e+15, cumulative = -8.37232e+15
ExecutionTime = 1.69 s  ClockTime = 2 s

Time = 0.1375

Courant Number mean: 2.13843e+32 max: 7.878e+36
DILUPBiCG:  Solving for Ux, Initial residual = 0.996796, Final residual = 1.95122, No Iterations 1001
DILUPBiCG:  Solving for Uy, Initial residual = 0.997762, Final residual = 8.63034, No Iterations 1001
DICPCG:  Solving for p, Initial residual = 1, Final residual = 5.35897, No Iterations 1001
time step continuity errors : sum local = 2.20198e+37, global = -2.42385e+21, cumulative = -2.42386e+21
DICPCG:  Solving for p, Initial residual = 0.990478, Final residual = 0.285598, No Iterations 1001
time step continuity errors : sum local = 5.67344e+38, global = 3.51519e+22, cumulative = 3.2728e+22
ExecutionTime = 1.95 s  ClockTime = 2 s

Time = 0.15

Courant Number mean: 1.17474e+39 max: 2.13026e+42
DILUPBiCG:  Solving for Ux, Initial residual = 0.999335, Final residual = 3.57986, No Iterations 1001
DILUPBiCG:  Solving for Uy, Initial residual = 0.996917, Final residual = 1.66825, No Iterations 1001
DICPCG:  Solving for p, Initial residual = 1, Final residual = 1.30652, No Iterations 1001
time step continuity errors : sum local = 2.24258e+43, global = 1.45234e+27, cumulative = 1.45237e+27
DICPCG:  Solving for p, Initial residual = 0.957264, Final residual = 1.25507, No Iterations 1001
time step continuity errors : sum local = 9.32738e+44, global = -5.20839e+27, cumulative = -3.75602e+27
ExecutionTime = 2.22 s  ClockTime = 2 s

Time = 0.1625

Courant Number mean: 2.50055e+45 max: 2.32987e+49
DILUPBiCG:  Solving for Ux, Initial residual = 0.999395, Final residual = 2.36754, No Iterations 1001
DILUPBiCG:  Solving for Uy, Initial residual = 0.999991, Final residual = 8.579, No Iterations 1001
DICPCG:  Solving for p, Initial residual = 1, Final residual = 14.3288, No Iterations 1001
time step continuity errors : sum local = 3.43893e+50, global = -3.69694e+34, cumulative = -3.69694e+34
DICPCG:  Solving for p, Initial residual = 0.986568, Final residual = 3.39767, No Iterations 1001
time step continuity errors : sum local = 7.05453e+52, global = -5.91511e+35, cumulative = -6.28481e+35
ExecutionTime = 2.48 s  ClockTime = 2 s

Time = 0.175

Courant Number mean: 1.26851e+53 max: 1.75204e+56
DILUPBiCG:  Solving for Ux, Initial residual = 0.998998, Final residual = 4.58537, No Iterations 1001
DILUPBiCG:  Solving for Uy, Initial residual = 0.999707, Final residual = 0.631888, No Iterations 1001
DICPCG:  Solving for p, Initial residual = 1, Final residual = 4.38271, No Iterations 1001
time step continuity errors : sum local = 1.06825e+57, global = -1.93055e+41, cumulative = -1.93056e+41
DICPCG:  Solving for p, Initial residual = 0.982177, Final residual = 14.2403, No Iterations 1001
time step continuity errors : sum local = 7.99535e+59, global = 1.62391e+43, cumulative = 1.60461e+43
ExecutionTime = 2.75 s  ClockTime = 3 s

Time = 0.1875

Courant Number mean: 1.63464e+60 max: 3.13789e+63
DILUPBiCG:  Solving for Ux, Initial residual = 0.998533, Final residual = 9.02987, No Iterations 1001
DILUPBiCG:  Solving for Uy, Initial residual = 0.997819, Final residual = 0.571402, No Iterations 1001
DICPCG:  Solving for p, Initial residual = 1, Final residual = 4.80908, No Iterations 1001
time step continuity errors : sum local = 9.14904e+64, global = -3.10405e+48, cumulative = -3.10404e+48
DICPCG:  Solving for p, Initial residual = 0.985965, Final residual = 0.660977, No Iterations 1001
time step continuity errors : sum local = 4.48705e+66, global = -2.57311e+50, cumulative = -2.60415e+50
ExecutionTime = 3.01 s  ClockTime = 3 s

Time = 0.2

Courant Number mean: 8.25292e+66 max: 1.12749e+71
DILUPBiCG:  Solving for Ux, Initial residual = 0.999802, Final residual = 0.740613, No Iterations 1001
DILUPBiCG:  Solving for Uy, Initial residual = 0.997251, Final residual = 0.921345, No Iterations 1001
DICPCG:  Solving for p, Initial residual = 1, Final residual = 0.720763, No Iterations 1001
time step continuity errors : sum local = 5.70309e+70, global = -2.47988e+53, cumulative = -2.48248e+53
DICPCG:  Solving for p, Initial residual = 0.961155, Final residual = 1.68471, No Iterations 1001
time step continuity errors : sum local = 6.78909e+72, global = -1.77559e+56, cumulative = -1.77807e+56
ExecutionTime = 3.27 s  ClockTime = 3 s

Time = 0.2125

Courant Number mean: 1.06959e+73 max: 2.0314e+76
DILUPBiCG:  Solving for Ux, Initial residual = 0.999582, Final residual = 1.20361, No Iterations 1001
DILUPBiCG:  Solving for Uy, Initial residual = 0.999582, Final residual = 1.32655, No Iterations 1001
DICPCG:  Solving for p, Initial residual = 1, Final residual = 1.88597, No Iterations 1001
time step continuity errors : sum local = 6.02945e+76, global = -1.78773e+60, cumulative = -1.78791e+60
DICPCG:  Solving for p, Initial residual = 0.993942, Final residual = 0.367036, No Iterations 1001
time step continuity errors : sum local = 4.62511e+78, global = 1.37298e+62, cumulative = 1.3551e+62
ExecutionTime = 3.53 s  ClockTime = 3 s

Time = 0.225

Courant Number mean: 1.55053e+79 max: 4.87871e+82
DILUPBiCG:  Solving for Ux, Initial residual = 0.999908, Final residual = 0.330781, No Iterations 1001
DILUPBiCG:  Solving for Uy, Initial residual = 0.999967, Final residual = 3.45014, No Iterations 1001
DICPCG:  Solving for p, Initial residual = 1, Final residual = 27.8568, No Iterations 1001
time step continuity errors : sum local = 5.8116e+84, global = -6.98023e+68, cumulative = -6.98023e+68
DICPCG:  Solving for p, Initial residual = 0.977871, Final residual = 10.1664, No Iterations 1001
time step continuity errors : sum local = 2.19292e+87, global = 7.59449e+70, cumulative = 7.52469e+70
ExecutionTime = 3.79 s  ClockTime = 4 s

Time = 0.2375

Courant Number mean: 3.84897e+87 max: 7.63325e+90
DILUPBiCG:  Solving for Ux, Initial residual = 0.998776, Final residual = 1.47147, No Iterations 1001
DILUPBiCG:  Solving for Uy, Initial residual = 0.999918, Final residual = 6.49108, No Iterations 1001
DICPCG:  Solving for p, Initial residual = 1, Final residual = 8.05985, No Iterations 1001
time step continuity errors : sum local = 2.34193e+92, global = 9.36871e+75, cumulative = 9.36878e+75
DICPCG:  Solving for p, Initial residual = 0.976013, Final residual = 1.47685, No Iterations 1001
time step continuity errors : sum local = 1.45774e+94, global = 6.55809e+77, cumulative = 6.65178e+77
ExecutionTime = 4.05 s  ClockTime = 4 s

Time = 0.25

Courant Number mean: 3.07094e+94 max: 1.31767e+98

I have successfully run simulations of one geometry, but when I make test on another geometry, the error shows.
The checkMesh is OK. The createPatch to create the cyclic boundaries is OK.
The timestep is 0.0125, even if I changed it to 1e-8, there is still such problems. And there is nothing new if I changed the Re to a very small value.
I suspect there is something wrong with the geometry. But the error shows up again for the third geometry. An example geometry is attached.
I wonder if you guys could give some instructions on this problem.

Thank you in advance.
Best Regards,
Bill

[piu58]

Large Courant numbers mean the time step ist too large. I see your time step is 10 ms. This is not especially small for complex geometries. Please keep in mind: Complex geometry = many nodes = small element sizes = small time for traversing an element = high Courant number.

The first steps have a small Courant number because the flow is zero initially. At least I suppose so.

If you can afford the simulation time use a much smaller time step. I guess you need far below 1 ms, rather 100 µs. If that ist not possible you have to reduce the complexity of the geometry. That may lead to coarse solutions for some regions of your model. In that case I recommend to use the coarse solution as a starting point for isolating a Region Of Interest and simulating this region isolated. This makes much more manual effort, of course.

[rdbisme]

Moreover if you use the adjustRunTime option in the controlDict together with the maxCo option you can let the solver choose the right timestep each time to stay under the MaxCo value of the Courant Number.

[6863523]

Dear,
Thank you for your replies.

Dr. Pilz,
I have tried to use timestep=1e-8s, the problem also appears. In my case, this is rather a small value and quite impossible for my simulation, since I need to run for 3e4s to ensure that I can observe the flow unsteadiness.

I have read through other threads
IcoFoam unstability, courant number gets large!
High Courant Number @ icoFoam
How to choose the right courant number?

It seems that this way does not work.

Dear Ruben,
According to what I have read, it seems that the adjustRunTime and MaxCo options are not designed for icoFoam.

I think I need to check the geometry over and over again, since one geometry is executable without this problem.

Any possible reason (except the geometry) for this problem is still welcomed.
Thank you.
Best,
Bill

[piu58]

If I understand the geometry right you try to simulate some kind of draining through a porous material like sand or rubble. If I am right I recommend:

- Use a very small clipping first
- Look at your boundary conditions. What is the b.c. for the particles?

General the first thoughts should be sacrificed to "What I want to find out" and next "How can this be simulated in an easy way". More complicated and complex conditions may added later. It would be helpful to know what you try to find out and which purpose your simulation serves.

~

> controlDict together with the maxCo option
This doesn't work with icoFoam, at least with my version of OF. pimpleFoam is able to accept these entries.

[6863523]

Dear Dr. Pilz,
Thank you for you detailed reply.
1) Reducing the time-step does not solve the problem, only delay the time for the occurance of the error.
2) You are right, I am modeling the flow in the porous media. I started from a simple case and then complex geometry. Interesting thing is that, one of the complex geometries works well, giving me the good results, while for other complex geometries, the courant number will increase.

The checkMesh function in OpenFOAM gives me the following information.

Time = 0

Mesh stats
    points:           1681124
    internal points:  0
    faces:            5732540
    internal faces:   2419290
    cells:            1630366
    faces per cell:   5
    boundary patches: 6
    point zones:      0
    face zones:       0
    cell zones:       0

Overall number of cells of each type:
    hexahedra:     0
    prisms:        1630366
    wedges:        0
    pyramids:      0
    tet wedges:    0
    tetrahedra:    0
    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
    LEFT                892      1800     ok (non-closed singly connected)
    RIGHT               892      1800     ok (non-closed singly connected)
    BOTTOM              841      1706     ok (non-closed singly connected)
    TOP                 841      1706     ok (non-closed singly connected)
    FRONTANDBACK        3260732  1681124  ok (non-closed singly connected)
    SOLID               49052    98184    ok (non-closed singly connected)

Checking geometry...
    Overall domain bounding box (0 0 0) (30 30 0.01)
    Mesh has 3 geometric (non-empty/wedge) directions (1 1 1)
    Mesh has 3 solution (non-empty) directions (1 1 1)
    Boundary openness (2.67981e-19 -1.23503e-19 3.0877e-14) OK.
    Max cell openness = 3.03166e-16 OK.
    Max aspect ratio = 7.5 OK.
    Minimum face area = 5.96e-05. Maximum face area = 0.00127061.  Face area magnitudes OK.
    Min volume = 5.96e-07. Max volume = 1.27061e-05.  Total volume = 7.53512.  Cell volumes OK.
    Mesh non-orthogonality Max: 71.1018 average: 4.27172
   *Number of severely non-orthogonal (> 70 degrees) faces: 9.
    Non-orthogonality check OK.
  <<Writing 9 non-orthogonal faces to set nonOrthoFaces
    Face pyramids OK.
    Max skewness = 0.624914 OK.
    Coupled point location match (average 0) OK.

Mesh OK.

End

Nothing special compared with the case that works well.

I wonder if there is any other possible reasons for this error, since all the settings in OpenFOAM are totally the same and there seems to be nothing wrong with the mesh.

I have looked some related threads, but no final reasons for this error.
Thank you.
Best,
Bill

[piu58]

It is hard to give an advice at this level of information.
I see two possible reasons:
1) There is something wrong with the mesh.
2) The problem gets instabil with larger complexity

You may set the end time to a point where the Courant number does not exceed 1 and look at the solution. If you have probems with the mesh you should find ares where the flow and pressure changes fast. Increasing Co number means that there is a very large flow somewhere.

If you find any mistake this way and believe in your mesh, the you could play around with underrelaxation.

Another way is slightly increasing the complexity and looking, when the problem starts to occur. In general, complex simulations are often more problematic because of numeric problems solving system of equations. If possible you may switch to an approximation technique. But I have experineces with that in OF.

But I would look at the mesh first, which is often the problematic part in such instabilities.

[laurentD]

*Number of severely non-orthogonal (> 70 degrees) faces: 9.     
Non-orthogonality check OK.   
<<Writing 9 non-orthogonal faces to set nonOrthoFaces

I have faced a similar problem. The 'checkMesh' utility checks the mesh (what a surprise !) without any flow conditions.
But the flow can be very complex in local areas. For example, in my case, there was a recirculation of the flow, with huge velocities, in the area where the severe non-orthogonal elements are detected.
Maybe you should take a look on the area concerned by this warning about the non-orthogonality, i guess that you will find a particular area with recirculations or things like that.

[piu58]

Areas with high flow need to meshed finer. A good measure of the real mesh quality is the relation between the max and mean Courant number, which should be low. You may get an idea of the flow values even from a coarse and imperfect mesh. Use that to establish a better mesh.

It is a fine idea only to use hex elements. Simulation is more stable in this case. Automated meshers don't reach the quality of a good manual mesh.

[6863523]

Dear Laurant,
Thank you for the reply.
Yes, you are right. This is most probably the reason for this problem.
I will try to figure the "severely non-orthogonal" problems.
Thank you.
Cheers,
Bill

[Mithun011]

Originally Posted by 6863523

Dear,
Thank you for your replies.

Dr. Pilz,
I have tried to use timestep=1e-8s, the problem also appears. In my case, this is rather a small value and quite impossible for my simulation, since I need to run for 3e4s to ensure that I can observe the flow unsteadiness.

I have read through other threads
IcoFoam unstability, courant number gets large!
High Courant Number @ icoFoam
How to choose the right courant number?

It seems that this way does not work.

Dear Ruben,
According to what I have read, it seems that the adjustRunTime and MaxCo options are not designed for icoFoam.

I think I need to check the geometry over and over again, since one geometry is executable without this problem.

Any possible reason (except the geometry) for this problem is still welcomed.
Thank you.
Best,
Bill

Hello Bill,

I was also facing the same problem with icoFoam, where the Courant Number increases to a very large number. But finally I resolved. It can be solved by the following:
1. Use pimpleFoam with laminar module.
2. Modify the mesh.
3. If you still want to use icoFoam and do not want to do any change in mesh, then manually calculate (roughly), the maximum velocity in steady state condition. Using the minimum mesh size that you had used, and limiting the Courant number to 1, determine the minimum time step required. With a bit of factor of safety (decrease the value a bit), put the time step in controlDict.
Hope this will work.

Thanks,
Jagannath

[sktsang]

I also have similar problems and doing searches in this forum. I currently using Abaqus CFD and runs satisfactorily. I wants to try out OpenFOAM becaujse of its capability in multiphase problems.


I tried running CFD for three intake manifold models for car engines, two out of three models have Courant numbers large right from the start of the simulations even though the checkMesh reports ok, and both simulations blew up very soon. When the simulations could still running I could see pimpleFoam do adjust the time step to even ~1e-12 s, as I turned on maxCo option (which set to 1, 0.45 or even 0.1).


What is interesting (or frustrating ... ) is that the same meshes ran perfectly in Abaqus CFD. The 3d meshes consists of tetrahedral and penta elements. All have penta elements modelling the boundary layers.

The geometry of the intake manifolds are always complex, making the mesh entirely out of brick (hexa) elements are just not possible, and because of relatively tight project schedule we cannot afford very long run time and very big problem size.


The models are intake manifolds for car engines, geometries were input, clean-up and meshed in ANSA (BetaCAD system). Snappyhexmesh was not used, as doing geometry clean-up is much easier in ANSA, and my computer runs on Window 7 (My OpenFOAM v5 is actually on Bluecore CFD).


I read the posts in this thread, if as piu58 said only hex elements are fine for OpenFOAM then I have a feeling that the problems I dealt with might be difficult for OpenFOAM ...

[rdbisme]

Originally Posted by sktsang

I also have similar problems and doing searches in this forum. I currently using Abaqus CFD and runs satisfactorily. I wants to try out OpenFOAM becaujse of its capability in multiphase problems.


I tried running CFD for three intake manifold models for car engines, two out of three models have Courant numbers large right from the start of the simulations even though the checkMesh reports ok, and both simulations blew up very soon. When the simulations could still running I could see pimpleFoam do adjust the time step to even ~1e-12 s, as I turned on maxCo option (which set to 1, 0.45 or even 0.1).


What is interesting (or frustrating ... ) is that the same meshes ran perfectly in Abaqus CFD. The 3d meshes consists of tetrahedral and penta elements. All have penta elements modelling the boundary layers.

The geometry of the intake manifolds are always complex, making the mesh entirely out of brick (hexa) elements are just not possible, and because of relatively tight project schedule we cannot afford very long run time and very big problem size.


The models are intake manifolds for car engines, geometries were input, clean-up and meshed in ANSA (BetaCAD system). Snappyhexmesh was not used, as doing geometry clean-up is much easier in ANSA, and my computer runs on Window 7 (My OpenFOAM v5 is actually on Bluecore CFD).


I read the posts in this thread, if as piu58 said only hex elements are fine for OpenFOAM then I have a feeling that the problems I dealt with might be difficult for OpenFOAM ...

Actually OpenFOAM works also with tetrahydra, just generally less smoothly.

By the way generally the problem you're facing is associated with wrong BCs. Are you sure that you translated the boundary conditions correctly from Abaqus to OF?

[sktsang]

Hi tidusuper91,

Just back, all three cases there is velocity at inlet and pressure at outlet, these are same in both Abaqus CFD and OpenFOAM runs.

Because of IP protection issue I cannot share the files, I summarise the BCs below, one from successful model, the other two from failed models:

This is from a successful model, running with Spalart Allmaras turb model:
Inlet:
p - zeroGradient
U - surfaceNormalFV, value -10 (m/s)
nut - calculated, value uniform 0
nuTilda - fixedValue, 0.01

Outlet (one of the engine cylinder port):
p - fixedValue, 0
U - zeroGradient
nut - calculated, value uniform 0
nuTilda - zeroGradient

Wall and all the other closed engine cylinder ports:
p - zeroGradient
U - noSlip
nut - nutkWallFunction, Cmu 0.09, kappa 0.41, epsilon 9.8, value uniform 0
nuTilda - zeroGradient

This is from a failed model, also running with Spalart Allmaras turb model:
Inlet:
p - zeroGradient
U - surfaceNormalFV, value -29.16 (m/s)
nut - calculated, value uniform 0
nuTilda - fixedValue, 0.01

Outlet (one of the engine cylinder port):
p - fixedValue, 0
U - zeroGradient
nut - calculated, value uniform 0
nuTilda - zeroGradient

Wall and all the other closed engine cylinder ports:
p - zeroGradient
U - noSlip
nut - nutkWallFunction, Cmu 0.09, kappa 0.41, epsilon 9.8, value uniform 0
nuTilda - zeroGradient

This is from another failed model, running with k-epsilon turb model:
Inlet:
p - zeroGradient
U - surfaceNormalFV, value -17.67(m/s)
k - fixed value, 1.
epsilon - fixedValue, 1.

Outlet (one of the engine cylinder port):
p - fixedValue, 0
U - zeroGradient
k - zeroGradient
epsilon - zeroGradient

Wall and all the other closed engine cylinder ports:
p - zeroGradient
U - noSlip
k - kqRWallFunction, value uniform 1.
epsilon - epsilonWallFunction, value uniform 1.

[Snigdha]

Hi,
I am also working with flow in the porous media (in between the gap of two spheres). I am using interFoam solver. After running for quite some time, the deltaT goes to very low value(1e-108 or something). I am using constantAlphaContactAngle condition on the sphere walls.

Interestingly, when I use zeroGradient boundary condition for alpha, the program runs well. Any idea what may be going wrong?

[rdbisme]

Originally Posted by Snigdha

Hi,
I am also working with flow in the porous media (in between the gap of two spheres). I am using interFoam solver. After running for quite some time, the deltaT goes to very low value(1e-108 or something). I am using constantAlphaContactAngle condition on the sphere walls.

Interestingly, when I use zeroGradient boundary condition for alpha, the program runs well. Any idea what may be going wrong?

Well, that hints that probably you're setting wrong BCs. Could it be?

[Snigdha]

Actually, I want to see the effect of contact angle on the liquid flow process. Hence, i want to put the solid liquid contact angle conditions.

Moreover, I made my mesh from sappy Hex mesh.

The following are my BCS for alpha, U and p_rgh respectively -





/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 7
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
object alpha.water;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions [0 0 0 0 0 0 0];

internalField uniform 0;

boundaryField
{
atmosphere
{
type zeroGradient;

}

two_sphere_bottom_sphere
{
type zeroGradient;



}

two_sphere_top_sphere
{
type zeroGradient;

}



}


/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 7
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volVectorField;
location "0";
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

internalField uniform (0 0 0);

boundaryField
{
atmosphere
{
type pressureInletOutletVelocity;
value uniform (0 0 0);
}


two_sphere_bottom_sphere
{
type noSlip;
}
two_sphere_top_sphere
{
type noSlip;
}


}
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 7
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
object p_rgh;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

internalField uniform 0;

boundaryField
{
atmosphere
{
type fixedValue;
value uniform 0;


}


two_sphere_bottom_sphere
{
type fixedFluxPressure;

}

two_sphere_top_sphere
{
type fixedFluxPressure;

}


}


Can you tell me what may be going wrong? if it is the mesh, then what shall improve?

[rdbisme]

Originally Posted by Snigdha

Actually, I want to see the effect of contact angle on the liquid flow process. Hence, i want to put the solid liquid contact angle conditions.

Moreover, I made my mesh from sappy Hex mesh.

The following are my BCS for alpha, U and p_rgh respectively -





/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 7
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
object alpha.water;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions [0 0 0 0 0 0 0];

internalField uniform 0;

boundaryField
{
atmosphere
{
type zeroGradient;

}

two_sphere_bottom_sphere
{
type zeroGradient;



}

two_sphere_top_sphere
{
type zeroGradient;

}



}


/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 7
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volVectorField;
location "0";
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

internalField uniform (0 0 0);

boundaryField
{
atmosphere
{
type pressureInletOutletVelocity;
value uniform (0 0 0);
}


two_sphere_bottom_sphere
{
type noSlip;
}
two_sphere_top_sphere
{
type noSlip;
}


}
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 7
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
object p_rgh;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

internalField uniform 0;

boundaryField
{
atmosphere
{
type fixedValue;
value uniform 0;


}


two_sphere_bottom_sphere
{
type fixedFluxPressure;

}

two_sphere_top_sphere
{
type fixedFluxPressure;

}


}


Can you tell me what may be going wrong? if it is the mesh, then what shall improve?

Just put your code within

CODE

BBcode for better readability.

What's your mesh? Showing up the mesh could actually help to understand better if you put the right BCs. Reading the name of the patches is not really enough

[Snigdha]

I have attached the wireframe of the mesh structure. The two half spheres are top sphere and bottom sphere respectively, and the other portions are atmosphere.

[rdbisme]

Originally Posted by Snigdha

I have attached the wireframe of the mesh structure. The two half spheres are top sphere and bottom sphere respectively, and the other portions are atmosphere.

As a first guess, I'll try to split the atmosphere patch in "inlet" and "outlet", and put an inlet condition for velocity and zeroGradient for pressure on "inlet", and zeroGradient for velocity and fixed pressure for "outlet".