[saidc.]

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:

 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.

[piu58]

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.

[piu58]

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.]

Dear Uwe,


checkMesh

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.

[Ahyar]

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

[piu58]

First you should repair your mesh:

***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.

[saidc.]

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.

[HPE]

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.

[saidc.]

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.