Turbulent wall treatment in multiphase flow regimes
 

Hello friends,

for quite some time now I am simulating several cases with multiphase flow regimes (3 phases) with openFOAM. I have found that the kEpsilon turbulence model is best for my setup and is easy to validate with a easy water-oil experiment.

We know that the dimensionless wall distance y+ is depending on the velocity, the actual wall distance and the kinematic viscosity of the fluid. Some one can easily calculate this value when it is a single phase system. But what value is the right one, if there is not one but three different values for kinematic viscosities? Problematic was also that I had a block mesh, which meant that I had large blocks on a large diameter shrinking to a very small block on a small diameter. I couldn't control my y+ very well due to the mesh layout and geometry.

I have run many test cases with my multiphase system on different meshes, finer and more coarse at the walls. I tested every different numerical scheme there is: limited, corrected, forward and backward ... I tested every turbulence model: kEpsilon, kOmegaSST ... Every time I received the same error: very high bounding values for the turbulent parameters. Mostly it was epsilon or omega. The result was always a crash of the solver.

I fixed this problem "simply" by not using any wall function for k, epsilon, nut! This is the only way I can run my cases without a crash. My question is now if there is a good value to choose for the turbulent parameters directly on the wall. Right now I am using the "fixedValue" option. Maybe there is a better alternative, I don't know. I still receive a bounding warnings, however this time, the values stay in a normal range.

P.S. If you start your case without a "nut" file and you want to use a "fixedValue" for the turbulent parameters, then openFOAM automatically creates an "nut" with a wall function. OpenFOAM then overwrites the k and epsilon or omega files and uses wall functions for these too. Therefore create a "nut" file with also a "fixedValue" at the wall in order to keep the original files from being overwritten.

Regards, Sebastian

Well, changing the BC of k and epsilon from wall functions to fixedValues didn't fix my problem. The crash happens at different physical times, which is weird because this doesn't look like a setup problem. If the BCs or the solver parameters were wrong, the crash would happen always at exactly the same time step. I am using OF 2.3 with MPI version 1.4.3. The error messages change also. I am thinking it maybe a problem with MPI, rather than OF since I could not reproduce the error in serial mode.


Regards, Sebastian


Attachment 45112

Attachment 45113

So, mpirun wasn't the problem either. Neither was the OF version or any of my BCs or ICs. I've build a smaller case where I had explicit non-orthogonal and skewed cells at the walls. CheckMesh had no problems with the mesh though. The maximum non-orthogonality was 45 and the skewness 1.3. This is more than ok.

This case didn't even run under normal circumstances. The problem was always high bounding epsilon and k values which either caused the iterations to exceed 1000 or caused a FPE which resulted in a solver crash.

I run my cases always with non-orthogonal correction because I always skewed cells in my meshes. I couldn't even believe to run this case turbulent on this mesh with multiphaseInterFoam anymore. Since I tried almost everything without success. However I found the problem in this case. What happened was that the solver tried to correct the epsilon and k values since the correction was higher than 1. The correction is however not limited when using the corrected schemes in laplacianSchemes and snGradSchemes. Using:

helped by limiting the maximum correction. The 0.333 means that the maximum correction has to be < 0.5 x the orthogonal part of the value. In the contrary the corrected schemes use explicit non-orthogonal correction.

Hope this helps someone :) Regards,

Sebastian

Attachment 45158