kEpsilon diverge
 

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

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.

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.

Dear Uwe,


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

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

First you should repair your mesh:
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.

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.

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.

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.