CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Running, Solving & CFD

kEpsilon diverge

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree6Likes
  • 1 Post By piu58
  • 1 Post By piu58
  • 1 Post By Ahyar
  • 3 Post By HPE

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   March 24, 2022, 11:46
Post kEpsilon diverge
  #1
Member
 
saidc
Join Date: Feb 2020
Location: Türkiye
Posts: 60
Rep Power: 6
saidc. is on a distinguished road
Hi,

I'm trying to simulate natural convection in packed-bed geometry. As seen in the geometry figure1, it consists of 300 spheres inside the cylinder. The bottom of the cylinder, its side-face and the spheres are defined as walls. Only the top of the cylinder is defined as the pressure outlet (no inlet). I use buoyantBoussinesqPimpleFoam as the solver and kEpsilon as the turbulence model. This case worked on laminar and SpalartAllmaras compositions and it converged, but the kEpsilon model diverges after a few iterations.

Things I've tried:
  • Reduced the time-step to 1e-5.
  • Reduced the relaxation factors.
  • Reduced turbulence schemes to first order (QUICK --> upwind).
The geometry has approximately 25 million cells and I don't want to change that number unless I have to. So refining the mesh is the last thing I would try.

Mesh:
Code:
 Mesh stats
     faces:     79630130

      cells:       24519096

 Checking geometry...
    Max aspect ratio = 7.23740009 OK.
    Minimum face area = 5.64402614e-11. Maximum face area = 3.4520711e-07.  Face area magnitudes OK.
    Min volume = 2.52506288e-15. Max volume = 1.66835424e-10.  Total volume = 0.000298417221.  Cell volumes OK.
    Mesh non-orthogonality Max: 55.0000107 average: 14.5384019
    Non-orthogonality check OK.
    Face pyramids OK.
 ***Max skewness = 4.7217643, 18 highly skew faces detected which may impair the quality of the results
  <<Writing 18 skew faces to set skewFaces
     Coupled point location match (average 0) OK.

Failed 1 mesh checks.

End
I guess the boundary conditions are not physical. I set the outlet as zeroGradient for all variable except pressure. The other 3 walls (side-face, bottom, spheres):
  • k, epsilon and nut: wallfunction
  • U: noSlip
  • p_rgh: fixedFluxPressure
  • T: spheres fixedValue, walls adiabatic
  • alphat: alphatJayatillekeWallFunction
fvSchemes and fvSolutions dictionaries attached. I can share more information if needed. I'm open to any suggestions.

Kind regards,
Said.
Attached Images
File Type: png figure1.png (12.9 KB, 12 views)
Attached Files
File Type: txt fvSchemes.txt (1.3 KB, 2 views)
File Type: txt fvSolution.txt (1.9 KB, 2 views)
saidc. is offline   Reply With Quote

Old   March 24, 2022, 13:34
Default
  #2
Senior Member
 
piu58's Avatar
 
Uwe Pilz
Join Date: Feb 2017
Location: Leipzig, Germany
Posts: 744
Rep Power: 15
piu58 is on a distinguished road
There may be problems with the mesh. A cylinder packed with spheres may arise a lopt of problems. I think of: Are all parts of the fluid connected to each other? No isolated regions?
What does checkMesh say?

Next: I don't recommend your setting with an output but no input. That may arise problems. What happens if (mathematical) the volume shrinks? I guess your spheres are hot and this is not physical possible. But when the simulation starst, some slightly strange thinks may happen until it is stabilized?

I recommend inletoutlet boundary condition. This prevents an unstable situation.
saidc. likes this.
__________________
Uwe Pilz
--
Die der Hauptbewegung überlagerte Schwankungsbewegung ist in ihren Einzelheiten so hoffnungslos kompliziert, daß ihre theoretische Berechnung aussichtslos erscheint. (Hermann Schlichting, 1950)
piu58 is offline   Reply With Quote

Old   March 24, 2022, 14:25
Default
  #3
Senior Member
 
piu58's Avatar
 
Uwe Pilz
Join Date: Feb 2017
Location: Leipzig, Germany
Posts: 744
Rep Power: 15
piu58 is on a distinguished road
Just another thing: Your setup is slightly strange. Yo u don't get a stable flow with this setup. You get a warming of the fluid with a slight stream upwards. It finishes when the all temperature and the fluid temperature is equal.

This is not a technical application. Form practical point of view, an inlet at the lower side and an output at the upper side gives some kind of reactor or similar.
saidc. likes this.
__________________
Uwe Pilz
--
Die der Hauptbewegung überlagerte Schwankungsbewegung ist in ihren Einzelheiten so hoffnungslos kompliziert, daß ihre theoretische Berechnung aussichtslos erscheint. (Hermann Schlichting, 1950)
piu58 is offline   Reply With Quote

Old   March 24, 2022, 14:38
Default
  #4
Member
 
saidc
Join Date: Feb 2020
Location: Türkiye
Posts: 60
Rep Power: 6
saidc. is on a distinguished road
Dear Uwe,


checkMesh
Code:
Mesh stats
    points:           32216753
    faces:            79630130
    internal faces:   72897652
    cells:            24519096
    faces per cell:   6.22077511
    boundary patches: 4
    point zones:      0
    face zones:       0
    cell zones:       2

Overall number of cells of each type:
    hexahedra:     18338817
    prisms:        893152
    wedges:        0
    pyramids:      0
    tet wedges:    7187
    tetrahedra:    0
    polyhedra:     5279940
    Breakdown of polyhedra by number of faces:
        faces   number of cells
            4   1439273
            5   959668
            6   1023468
            7   12187
            8   13110
            9   787164
           10   2128
           11   7550
           12   574076
           13   472
           14   5533
           15   413901
           16   206
           17   1587
           18   39279
           19   88
           20   108
           21   141
           22   1

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
                  wall_0   109487   112105  ok (non-closed singly connected)
                 z_min_0    15643    16084  ok (non-closed singly connected)
                 z_max_0    13300    13441  ok (non-closed singly connected)
           STL_spheres_0  6594048  8290445      ok (closed singly connected)

Checking faceZone topology for multiply connected surfaces...
    No faceZones found.

Checking basic cellZone addressing...
                CellZone        Cells       Points       VolumeBoundingBox
                 square      9379283     12413686 9.27568804e-05 (-0.016302401 -0.0163072004 0) (0.0162903377 0.0163140554 0.151)
            innerCircle     15139813     19983201 0.00020566034 (-0.030500095 -0.0305000242 0) (0.0305001003 0.030500022 0.151)

Checking geometry...
    Overall domain bounding box (-0.030500095 -0.0305000242 0) (0.0305001003 0.030500022 0.151)
    Mesh has 3 geometric (non-empty/wedge) directions (1 1 1)
    Mesh has 3 solution (non-empty) directions (1 1 1)
    Boundary openness (2.97417334e-16 -7.74866121e-17 3.17896828e-16) OK.
    Max cell openness = 4.69920191e-16 OK.
    Max aspect ratio = 7.23740009 OK.
    Minimum face area = 5.64402614e-11. Maximum face area = 3.4520711e-07.  Face area magnitudes OK.
    Min volume = 2.52506288e-15. Max volume = 1.66835424e-10.  Total volume = 0.000298417221.  Cell volumes OK.
    Mesh non-orthogonality Max: 55.0000107 average: 14.5384019
    Non-orthogonality check OK.
    Face pyramids OK.
 ***Max skewness = 4.7217643, 18 highly skew faces detected which may impair the quality of the results
  <<Writing 18 skew faces to set skewFaces
    Coupled point location match (average 0) OK.

Failed 1 mesh checks.
End
I think there are no isolated regions between spheres. For making geometry same as reality I used bridge (local) contact point approximation so there are cylinders on every contact point. Therefore meshing process is not overwhelming for snappyhexmesh and due to this the geometry must have been snapped realistic.

I'll try inletOutlet BC at outlet for all variables. However, mostly it will apply zeroGradient because Re Number is around 300, there is no big velocities and eddies (not sure) so as you said without any inlet patch it may be explode again. What if I just change bottom wall with inlet patch? But I'm not sure which boundary condition composition will be physical for my problem because there should be a wall. Is this combination poor?
  • p: zeroGradient
  • U: fixedValue (0) or outletInlet
  • k, epsilon, nut: outletInlet with 0 value (or fixedValue)
Kind regards,
Said.
saidc. is offline   Reply With Quote

Old   March 24, 2022, 21:13
Default
  #5
Member
 
Muhammad Ahyar
Join Date: Mar 2020
Posts: 30
Rep Power: 6
Ahyar is on a distinguished road
i used to face the similar problem where the kEpsilon diverge. What i did to solve that problem was changing the 3D model format from .stl to .obj (I used the snappyHexMesh tool to create the mesh), and then my simulation ran well. Until now, i don't know why that method worked
saidc. likes this.
Ahyar is offline   Reply With Quote

Old   March 25, 2022, 02:11
Default
  #6
Senior Member
 
piu58's Avatar
 
Uwe Pilz
Join Date: Feb 2017
Location: Leipzig, Germany
Posts: 744
Rep Power: 15
piu58 is on a distinguished road
First you should repair your mesh:
Code:
***Max skewness = 4.7217643, 18 highly skew faces detected which may impair the quality of the results
Failed 1 mesh checks.
Probably there are way simplifying your geometry at large. What is your real world application? Have you ever thought if it is sufficient simulation only one or a few spheres?

This type of simulation is thought for boundaries produce by gravity. Because you have spheres, there is no continuous boundary, there are always free flow regions between them. Looking what happens at a single sphere should give much insight.
__________________
Uwe Pilz
--
Die der Hauptbewegung überlagerte Schwankungsbewegung ist in ihren Einzelheiten so hoffnungslos kompliziert, daß ihre theoretische Berechnung aussichtslos erscheint. (Hermann Schlichting, 1950)
piu58 is offline   Reply With Quote

Old   March 25, 2022, 05:06
Default
  #7
Member
 
saidc
Join Date: Feb 2020
Location: Türkiye
Posts: 60
Rep Power: 6
saidc. is on a distinguished road
Hi Uwe,

Yes, firstly we did simulate bcc-fcc geometries with periodic bcs, but now we need to simulate all geometry (full-core).

After I get the results of inlet patch and inletOutlet/outletInlet bcs I will post the results here. I hope it is about the boundary conditions because I don't know how to improve this mesh any more. If it diverge again I'll try skewCorrected schemes as my last chance. If it diverges again, I will look for new methods to improve the Mesh.

Kind regards,
Said.

Last edited by saidc.; March 25, 2022 at 07:17.
saidc. is offline   Reply With Quote

Old   March 26, 2022, 15:48
Default
  #8
HPE
Senior Member
 
HPE's Avatar
 
Herpes Free Engineer
Join Date: Sep 2019
Location: The Home Under The Ground with the Lost Boys
Posts: 932
Rep Power: 12
HPE is on a distinguished road
I wouldn't use the standard kEpsilon for wall-bounded flow applications, since it was derived as a high-Re number model. Either consider RNGkEpsilon or realizableKE, or just the kOmegaSST which involves the kEpsilon in the freestream anyway. If you insist on it, avoid the y+ being below 30-40 (that may require remeshing.)

Also, unlike ANSYS, epsilon-based models are a bit unstable in OpenFOAM. I speculate that the software that provide stable-epsilon models use the "homogeneous" epsilon field in the background (i.e. epsilonH = epsilon - nonHomogeneousEpsilon) to compute epsilon or use implicit treatments in the epsilon equation (e.g. nextFoam). The standard OpenFOAM does not have any of the two to stabilise epsilon iterations.
piu58, hogsonik and saidc. like this.
HPE is offline   Reply With Quote

Old   March 31, 2022, 07:09
Post
  #9
Member
 
saidc
Join Date: Feb 2020
Location: Türkiye
Posts: 60
Rep Power: 6
saidc. is on a distinguished road
Hi again,

I reduced the skewness, now checkMesh does not give any fail. However, there is something strange in my results. Before I try two equation turbulence models I have to find the reason why I couldn't achive correct results with laminar case and one equation model. Lets say, Experimental and LBM results are reaches 1 m/s velocity but I only reach 0.25 m/s and there is really a big difference (Measurements are made at a height where flow can develop). Also, their heat transfer occours so fast like 3-4 seconds in real time, but on my case have to run the simulation 60 seconds for just reach 0.1 m/s because the heat transfer not occours fast and correctly. Although the walls are adiabatic, the temperature values decrease to lower than the initial value.


With all this in mind, I think I'm doing something non-physical. I'll edit here if I can find what I missed.


Edit: The problem is the residual control tolerance was too low for converge. I tried same test case with different residual tolerance (1e-2 and 1e-4) and with 1e-2 the convergence time has doubled.


Kind regards,
Said.

Last edited by saidc.; April 19, 2022 at 10:07.
saidc. is offline   Reply With Quote

Reply

Tags
boussinesq, buoyantboussinesqpimple, natural conection

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
kEpsilon vs. kOmegaSST Tobi OpenFOAM Running, Solving & CFD 4 February 14, 2020 16:04
simulation supersonic external flow arround the missile diverge!! raminostadi FLUENT 2 April 14, 2018 13:13
Make wall-function read a newly defined field inside modified kEpsilon model Radunz OpenFOAM Programming & Development 2 July 18, 2017 23:39
pressure diverge for pipeflow using rhoSimpleFoam waq OpenFOAM Running, Solving & CFD 1 June 25, 2015 09:38
kEpsilon divergence s.m OpenFOAM Running, Solving & CFD 0 May 27, 2013 10:30


All times are GMT -4. The time now is 04:55.