[mwaqas]

Hello Everyone

I am using OF 2.1.1 to simulate fluidized bed (twoPhaseEulerPoam). My geometry is 3D cylinder. My simulation works fine upto 0.6 sec and courant number remains well below 1 e.g. Courant Number mean: 0.0015972 max: 0.102484 Max Ur Courant Number = 0.0864022

But after 0.6 sec, it stats diverging and suddenly increases up to 1000. see the log file. What could be the possible problem. Thank you


[QUOTE][/Courant Number mean: 0.00200133 max: 0.0987868
Max Ur Courant Number = 0.0853136
deltaT = 1.21324e-05
Time = 0.636997

DILUPBiCG: Solving for alpha, Initial residual = 6.99627e-05, Final residual = 4.09665e-12, No Iterations 2
Dispersed phase volume fraction = 0.243807 Min(alpha) = -2.1287e-05 Max(alpha) = 0.655422
DILUPBiCG: Solving for alpha, Initial residual = 1.04443e-07, Final residual = 1.24033e-11, No Iterations 1
Dispersed phase volume fraction = 0.243807 Min(alpha) = -4.54623e-07 Max(alpha) = 0.655422
kinTheory: max(Theta) = 1000
kinTheory: min(nua) = 9.16968e-11, max(nua) = 0.04
kinTheory: min(pa) = -1.92433e-08, max(pa) = 271622
GAMG: Solving for p, Initial residual = 0.0387725, Final residual = 0.00108617, No Iterations 1
time step continuity errors : sum local = 1.63039e-07, global = -3.89795e-08, cumulative = -1.47341e-06
GAMG: Solving for p, Initial residual = 0.00473594, Final residual = 6.99331e-09, No Iterations 15
time step continuity errors : sum local = 1.05671e-12, global = -3.54838e-13, cumulative = -1.47341e-06
DILUPBiCG: Solving for epsilon, Initial residual = 0.126719, Final residual = 3.52361e-08, No Iterations 1
DILUPBiCG: Solving for k, Initial residual = 0.0635536, Final residual = 1.17844e-08, No Iterations 3
ExecutionTime = 1456.39 s ClockTime = 1460 s

Courant Number mean: 0.00202643 max: 0.099791
Max Ur Courant Number = 0.0877428
deltaT = 1.21558e-05
Time = 0.637009

DILUPBiCG: Solving for alpha, Initial residual = 7.00857e-05, Final residual = 3.30245e-12, No Iterations 2
Dispersed phase volume fraction = 0.243807 Min(alpha) = -1.94546e-05 Max(alpha) = 0.655419
DILUPBiCG: Solving for alpha, Initial residual = 9.28193e-08, Final residual = 1.58776e-11, No Iterations 1
Dispersed phase volume fraction = 0.243807 Min(alpha) = -7.98219e-07 Max(alpha) = 0.655419
kinTheory: max(Theta) = 1000
kinTheory: min(nua) = 9.18938e-11, max(nua) = 0.04
kinTheory: min(pa) = -3.40219e-08, max(pa) = 262061
GAMG: Solving for p, Initial residual = 0.0629067, Final residual = 0.00327559, No Iterations 1
time step continuity errors : sum local = 5.09797e-07, global = 3.84477e-07, cumulative = -1.08893e-06
GAMG: Solving for p, Initial residual = 0.522737, Final residual = 6.48233e-09, No Iterations 22
time step continuity errors : sum local = 6.07518e-12, global = -2.00598e-12, cumulative = -1.08894e-06
DILUPBiCG: Solving for epsilon, Initial residual = 0.774542, Final residual = 2.43972e-13, No Iterations 3
DILUPBiCG: Solving for k, Initial residual = 0.999999, Final residual = 1.83465e-08, No Iterations 11
ExecutionTime = 1458.38 s ClockTime = 1462 s

Courant Number mean: 0.00243195 max: 2.85304
Max Ur Courant Number = 2.90461
deltaT = 4.18497e-07
--> FOAM Warning :
From function Time:perator++()
in file db/Time/Time.C at line 1010
Increased the timePrecision from 6 to 7 to distinguish between timeNames at time 0.637009
Time = 0.6370091

DILUPBiCG: Solving for alpha, Initial residual = 2.85712e-06, Final residual = 6.15549e-12, No Iterations 2
Dispersed phase volume fraction = 0.243807 Min(alpha) = -7.47147e-07 Max(alpha) = 0.655419
DILUPBiCG: Solving for alpha, Initial residual = 4.58626e-09, Final residual = 1.00269e-11, No Iterations 1
Dispersed phase volume fraction = 0.243807 Min(alpha) = -7.47147e-07 Max(alpha) = 0.655419
kinTheory: max(Theta) = 1000
kinTheory: min(nua) = 9.25553e-11, max(nua) = 0.04
kinTheory: min(pa) = -3.18112e-08, max(pa) = 255560
GAMG: Solving for p, Initial residual = 0.949397, Final residual = 0.0728741, No Iterations 1
time step continuity errors : sum local = 1.14736e-06, global = 8.30162e-07, cumulative = -2.58774e-07
GAMG: Solving for p, Initial residual = 0.00102808, Final residual = 7.87398e-09, No Iterations 16
time step continuity errors : sum local = 8.82252e-10, global = 2.98332e-10, cumulative = -2.58476e-07
DILUPBiCG: Solving for epsilon, Initial residual = 0.233328, Final residual = 1.09165e-13, No Iterations 1
DILUPBiCG: Solving for k, Initial residual = 0.964855, Final residual = 7.86549e-13, No Iterations 1
ExecutionTime = 1459.32 s ClockTime = 1463 s

Courant Number mean: 0.000215521 max: 0.289072
Max Ur Courant Number = 26.6629
deltaT = 1.56958e-09
--> FOAM Warning :
From function Time:perator++()
in file db/Time/Time.C at line 1010
Increased the timePrecision from 7 to 8 to distinguish between timeNames at time 0.637009
Time = 0.63700909

DILUPBiCG: Solving for alpha, Initial residual = 2.49435e-07, Final residual = 2.05948e-11, No Iterations 1
Dispersed phase volume fraction = 0.243807 Min(alpha) = -7.47016e-07 Max(alpha) = 0.655419
DILUPBiCG: Solving for alpha, Initial residual = 5.46278e-10, Final residual = 1.11052e-13, No Iterations 1
Dispersed phase volume fraction = 0.243807 Min(alpha) = -7.47016e-07 Max(alpha) = 0.655419
kinTheory: max(Theta) = 1000
kinTheory: min(nua) = 9.24831e-11, max(nua) = 0.04
kinTheory: min(pa) = -3.18415e-08, max(pa) = 1.18596e+06
GAMG: Solving for p, Initial residual = 0.999857, Final residual = 0.0729551, No Iterations 1
time step continuity errors : sum local = 0.0181995, global = -0.00461067, cumulative = -0.00461093
GAMG: Solving for p, Initial residual = 0.0615916, Final residual = 8.77243e-09, No Iterations 24
time step continuity errors : sum local = 9.23584e-09, global = -3.17977e-09, cumulative = -0.00461094
DILUPBiCG: Solving for epsilon, Initial residual = 0.87258, Final residual = 6.82741e-06, No Iterations 10
DILUPBiCG: Solving for k, Initial residual = 0.917014, Final residual = 41215.3, No Iterations 1001
ExecutionTime = 1465.41 s ClockTime = 1469 s

Courant Number mean: 0.572073 max: 260.204
Max Ur Courant Number = 13409.9
deltaT = 1.17047e-14
--> FOAM Warning :
From function Time:perator++()
in file db/Time/Time.C at line 1010
Increased the timePrecision from 8 to 9 to distinguish between timeNames at time 0.637009
Time = 0.637009086

DILUPBiCG: Solving for alpha, Initial residual = 1.13753e-07, Final residual = 3.97845e-12, No Iterations 1
Dispersed phase volume fraction = 0.243807 Min(alpha) = -7.29461e-07 Max(alpha) = 0.655419
DILUPBiCG: Solving for alpha, Initial residual = 3.31384e-10, Final residual = 4.6151e-14, No Iterations 1
Dispersed phase volume fraction = 0.243807 Min(alpha) = -7.2946e-07 Max(alpha) = 0.655419
kinTheory: max(Theta) = 1000
kinTheory: min(nua) = 9.1779e-11, max(nua) = 0.04
kinTheory: min(pa) = -1.82365, max(pa) = 1.89181e+08
GAMG: Solving for p, Initial residual = 0.961636, Final residual = 0.0531763, No Iterations 1
time step continuity errors : sum local = 1.45859e-07, global = -8.9757e-09, cumulative = -0.00461094
GAMG: Solving for p, Initial residual = 0.144315, Final residual = 1.49241e-07, No Iterations 1000
time step continuity errors : sum local = 2.90486e-06, global = 1.28668e-07, cumulative = -0.00461082]

Regards

[sharonyue]

Hello dude,

From my experience, this solver is mesh sensitive, So I suspect this might be the problem. could u pls tell us which kind of mesh do u use? tet or hex?

Best,

[mwaqas]

Hello sharonyue

Thank you for your reply. I am using hexa mesh generated by using SHM (snappyHexMesh). Here is the my checkMesh result.


Mesh stats
points: 42241
faces: 122452
internal faces: 116332
cells: 40144
boundary patches: 4
point zones: 0
face zones: 0
cell zones: 0

Overall number of cells of each type:
hexahedra: 38064
prisms: 2080
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
allBoundary 0 0 ok (empty)
inlet_inlet 772 797 ok (non-closed singly connected)
outlet_outlet 772 797 ok (non-closed singly connected)
wall_wall 4576 4664 ok (non-closed singly connected)

Checking geometry...
Overall domain bounding box (-0.11499987 -0.1148842 0) (0.11499986 0.1148842 0.38)
Mesh (non-empty, non-wedge) directions (1 1 1)
Mesh (non-empty) directions (1 1 1)
Boundary openness (3.4623065e-17 7.0086455e-16 -1.2106427e-16) OK.
Max cell openness = 3.0660809e-16 OK.
Max aspect ratio = 3.499081 OK.
Minumum face area = 1.5619455e-05. Maximum face area = 8.1890858e-05. Face area magnitudes OK.
Min volume = 9.1210546e-08. Max volume = 5.4318036e-07. Total volume = 0.015750168. Cell volumes OK.
Mesh non-orthogonality Max: 24.633834 average: 3.0930508
Non-orthogonality check OK.
Face pyramids OK.
Max skewness = 0.7049643 OK.
Coupled point location match (average 0) OK.

Mesh OK.

Regards

[sharonyue]

Hello,

From your log, looks like this is caused of your k and epsilon fields. have u ever tried with running it in laminar? Or maybe u can modify fvscheme of k and epsilon. Also u can check your latest result in paraview to see if epsilon is very big. and see if this thing happens in some bad cells.

Best,
Dongyue

[mwaqas]

Hello Dongyue

I also tried with laminar flow, the result is almost same. With laminar flow simulation didn't crash but delta t (time step) has decreased up to 10*e-11 after 0.4 sec of simulation. Initially it was taken 0.00001 and courant number was well below 0.1.
Secondly I also tried with 2D mesh (simple rectangular geometry), all cells are of same size and no bad cell is present. Again my simulation crashed after 0.4 sec

PS: for low velocity my simulation result is fine and it didn't crash but at higher velocity it crashes (even at high velocity courant number is smaller than 0.1 till 0.4 sec)

my fvScheme is as follows. can you give some suggestions please about it

ddtSchemes
{
default Euler;
}

gradSchemes
{
default Gauss linear;
}

divSchemes
{
default none;
div(phia,Ua) Gauss limitedLinearV 1;
div(phib,Ub) Gauss limitedLinearV 1;
div(phib,k) Gauss limitedLinear 1;
div(phib,epsilon) Gauss limitedLinear 1;
div(phi,alpha) Gauss limitedLinear01 1;
div(phir,alpha) Gauss limitedLinear01 1;
div(phi,Theta) Gauss limitedLinear 1;
div(Rca) Gauss linear;
div(Rcb) Gauss linear;
}

laplacianSchemes
{
default none;
laplacian(nuEffa,Ua) Gauss linear corrected;
laplacian(nuEffb,Ub) Gauss linear corrected;
laplacian((rho*(1|A(U))),p) Gauss linear corrected;
laplacian(alphaPpMag,alpha) Gauss linear corrected;
laplacian(Galphaf,alpha) Gauss linear corrected;
laplacian(DkEff,k) Gauss linear corrected;
laplacian(DepsilonEff,epsilon) Gauss linear corrected;
laplacian(kappa,Theta) Gauss linear corrected;
}

interpolationSchemes
{
default linear;
}

snGradSchemes
{
default corrected;
}

fluxRequired
{
default no;
p ;
}

Regards

[BlnPhoenix]

Hi, did u manage to stabalize your simulation?

I'm also experiencing crashes with twoPhaseEulerFoam and i am not sure what the cause could be.
I'm running in laminar also. Do you have some experience with this solver u can share?

Thanks!

[sharonyue]

Hi guys,

I just found that in his Checkmesh log, there are prisms. I dont think twoPhaseEulerFoam can handle prisms. From my experience it can only paly with hex.

Best,

[BlnPhoenix]

Hi sharonyue,

this is a very interesting comment! I have indeed various types of cells in my geometry. I'm importing the mesh form ccm+ and actually i am confused i have all these types of cells, i thougth i would only have polyhedras. So i make from your post i should switch to 100% hexahedras. I will try it!

Thanks




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

Create polyMesh for time = 0

Time = 0

Mesh stats
points: 2062262
faces: 2403400
internal faces: 2350866
cells: 336222
faces per cell: 14.1403
boundary patches: 4
point zones: 0
face zones: 0
cell zones: 0

Overall number of cells of each type:
hexahedra: 317
prisms: 531
wedges: 3
pyramids: 1
tet wedges: 3
tetrahedra: 0
polyhedra: 335367
Breakdown of polyhedra by number of faces:
faces number of cells
5 4
6 22
7 513
8 2349
9 7081
10 16375
11 22209
12 29420
13 44478
14 57475
15 57773
16 45517
17 28370
18 14311
19 6169
20 2235
21 753
22 239
23 60
24 10
25 4

Checking topology...
Boundary definition OK.
Cell to face addressing OK.
Point usage OK.
<<Found 93 neighbouring cells with multiple inbetween faces.
Upper triangular ordering OK.
<<Writing 186 unordered faces to set upperTriangularFace
Face vertices OK.
Number of regions: 1 (OK).

Checking patch topology for multiply connected surfaces...
Patch Faces Points Surface topology
Body2.Wall 43291 84241 ok (non-closed singly connected)
Body2.Outlet 1705 3416 ok (non-closed singly connected)
Body2.Edges 2645 5290 ok (non-closed singly connected)
Body2.Inlet 4893 9786 ok (non-closed singly connected)

Checking geometry...
Overall domain bounding box (0 0 -1.04157e-015) (1.34 0.31 0.5)
Mesh (non-empty, non-wedge) directions (1 1 1)
Mesh (non-empty) directions (1 1 1)
Boundary openness (-1.6881e-016 3.15279e-017 3.68234e-016) OK.
Max cell openness = 2.53791e-016 OK.
Max aspect ratio = 8.16626 OK.
Minimum face area = 4.03175e-010. Maximum face area = 0.000356587. Face area magnitudes OK.
Min volume = 3.27702e-012. Max volume = 1.14262e-005. Total volume = 0.159208. Cell volumes OK.
Mesh non-orthogonality Max: 73.5329 average: 14.3827
*Number of severely non-orthogonal (> 70 degrees) faces: 4.
Non-orthogonality check OK.
<<Writing 4 non-orthogonal faces to set nonOrthoFaces
Face pyramids OK.
***Max skewness = 5.75398, 1 highly skew faces detected which may impair the quality of the results
<<Writing 1 skew faces to set skewFaces
Coupled point location match (average 0) OK.

Failed 1 mesh checks.

End

[BlnPhoenix]

Hi sharonyue,

its me again. I'm having a hard time getting rid of the prism cells in my geometry. I thought if i use snappyhex mesh it would only create hexadra cells but i also get these prism cells u said don't go well with twophaseeulerfoam. Since u worked with this solver can u give me a hint how i can avoid these bad cell types? is there any way around the solver crashes?

thanks a lot for ur suggestions!

[sharonyue]

Originally Posted by BlnPhoenix

Hi sharonyue,

its me again. I'm having a hard time getting rid of the prism cells in my geometry. I thought if i use snappyhex mesh it would only create hexadra cells but i also get these prism cells u said don't go well with twophaseeulerfoam. Since u worked with this solver can u give me a hint how i can avoid these bad cell types? is there any way around the solver crashes?

thanks a lot for ur suggestions!

Hi BlnPhoenix

Yeah, I use this solver to simulate stirred tanks. From my experience and the bug reports, twoPhaseEulerFoam and compressibleTwoPhaseEulerFoam need a perfect mesh, especially in turbulent situation.

Here are a bug reports:

http://www.openfoam.org/mantisbt/view.php?id=1379

Heney address this problem in OpenFOAM-dev. I did not try it, but I tried it with OpenFOAM 22x and 23x. Did not work.

Another example for me is this:



For this kind of irregular geometry, at first I tried with snappyHexMesh as of its simplicity. But it will diverge when running in turbulent simulation. It will tell you the epsilon (or k) field is quite large, then you will see the deltaT will be below of 1e-7 like that. So for me a absolutely perfect mesh is the only way.

Even for the hex mesh above, I need to pay lots of time to make it better near the stick region.

So I think at first you can try OpenFOAM-dev to try the new algorithm.

In the mean time, you can write the bad fields in some time step then check it in paraview. I guess the velocity is quite large. Check if this quite large velocity locates in the bad mesh cells.

Why dont u upload your geometry? Lets see if its possible for a hex mesh.

[masoudsh]

Hi Dongyue

I see the errrors you me mentions:

"It will tell you the epsilon (or k) field is quite large, then you will see the deltaT will be below of 1e-7 like that. So for me a absolutely perfect mesh is the only way."

Do you solve it?

Best Regards
Masoud

[hooman.4028]

Hello BlnPhoenix and Dongyue,

Hope all is well with you. I am trying to use twoPhaseEulerFoam and reactingTwoPhaseEulerFoam. I tried reactingTwoPhaseEulerFoam because I read that it is more stable that the twoPhaseEulerFoam. But unfortunately both blow up. My temperatures go to huge number even though I do not solve for h/e equations (i simply put minIter and maxIter to zero). I used SHM to mesh my geometry and I do have prims mesh elements. I also consider a turbulent flow (k-E). Can you please let me know how you solved your problem?

I am using OF 4.1!

Thank you,