[robbirobocop]

Dear Foamers,

in the last weeks I tried to successfully use bubbleFoam to simulation multiphase flow inside a steamdrum. After finally getting the simulation more or less stable - the timeStep and maxCo are still very low - yet another problem occurs.

I am loosing water. There is one inlet where a steam water mixture is put in. After some time that is 0.05 s in the attached picture it seems like only steam was put in and not a mixture...

And I have no idea what BC or scheme or whatever arouses this problem. It seems that steam only pushes away the water level that was initialised with setFields. No separation takes place.

The used fvSchemes and fvSolution entries are as follows:

fvSchemes

ddtSchemes
{
    default         Euler;
}
gradSchemes
{
    default         Gauss linear;
}
divSchemes
{
    default         none;
    div(phia,Ua)    Gauss upwind;
    div(phib,Ub)    Gauss upwind;
    div(phib,k)     Gauss upwind;
    div(phib,epsilon) Gauss upwind;
    div(phi,alpha)  Gauss upwind;
    div((-nuEffa*T(grad(Ua)))) Gauss linear;
    div((-nuEffb*T(grad(Ub)))) Gauss linear;
}
laplacianSchemes
{
    default         Gauss linear corrected;
}
interpolationSchemes
{
    default         linear;
}
snGradSchemes
{
    default         corrected;
}
fluxRequired
{
    default         no;
    p               ;
}

fvSolution

solvers
{
    p
    {
        solver          GAMG;
        tolerance       1e-04;
        relTol          0.1;
        smoother        DIC;
        nPreSweeps      0;
        nPostSweeps     2;
        nFinestSweeps   2;
        cacheAgglomeration true;
        nCellsInCoarsestLevel 10;
        agglomerator    faceAreaPair;
        mergeLevels     1;
    }
    pFinal
    {
        $p;
        tolerance       1e-04;
        relTol          0;
    }
       "(k|epsilon)"
    {
        solver          PBiCG;
        preconditioner  DILU;
        tolerance       1e-4;
        relTol          0.0;
    }
    "(k|epsilon)Final"
    {

        solver          PBiCG;
        preconditioner  DILU;
        tolerance       1e-4;
        relTol          0.0;
    }
    alpha
    {
        solver          PBiCG;
        preconditioner  DILU;
        tolerance       1e-5;
        relTol          0.1;
    }
    alphaFinal
    {
        solver          PBiCG;
        preconditioner  DILU;
        tolerance       1e-5;
        relTol          0;
    }
}
    nCorrectors     1; //1
    nNonOrthogonalCorrectors 1;
    nAlphaCorr      1;
    correctAlpha    yes;

The BC for the INLET is alpha of 0.6, Ua with ~6 m/s, Ub with ~4 m/s and p is of course zeroGradient. K and epsilon are calculated with the typical formulaes and set to a fixedValue. A normal water level is set with setFields. Thus, half of the geometry is filled with water.

The problem occurs 2D as well as 3D, laminar and turbulent...
The picture that is attached, illustrates the "loosing of the fluid"... (right side / patch)

Any idea or help would be very appreciated.

Thanks.

[alberto]

It would be useful to see a small case that reproduces the problem. Also, what version of OpenFOAM are you using?

Your linear solver tolerances are way too high however. A tolerance of 1.0e-8 over p, 10^-10 over alpha, and 10^-5 for the other variables is what I usually adopt. Also, if your grid is not skewed, set the nonOrthogonalCorrectors to 0, and perform 2 alpha correctors, and 2 nCorrectors.

I would say that using upwind should not be necessary, but this depend on the mesh quality :-)

[robbirobocop]

Well, if you would send me your e-mail address via this forum I can send you the case. Other users are welcome, too.

Since the grid is generated with ICEM (the snappy grid is not too good) the case directory is about 20 MB.

I am using OpenFoam 2.0.1.

I used the higher tolerances in order to decrease the computational time. The aim of this 2D case was actually to reproduce a first "quick and dirty" solution to see what problems might occur that could have major influence on the 3D case. With lower tolerances the stability is even worse.

And I found out that the picture I attached is the initialisation with setFields. So basically I am not loosing water. It is just that nothing really happens, although I have 4/6 m/s of water/steam that goes in. So after a time of 0.05 seconds (which already took a lot of time on an i7 - 8GB RAM ... system) multiplied by a speed of 4 or 6 m/s would tend to result in the water moving around 0.2 m and steam around 0.3 m. But as can be seen in the picture neither the first one nor the second one happens.

I then realized that the characteristic length is a bit smaller. Thus, epsilon is a lot higher. Now, with the new initialisation the simulation crashes even faster...

[alberto]

Originally Posted by robbirobocop

Well, if you would send me your e-mail address via this forum I can send you the case. Other users are welcome, too.

Since the grid is generated with ICEM (the snappy grid is not too good) the case directory is about 20 MB.

OK, I will send you a message. Hopefully your mesh does not have millions of cells (I'd prefer a 2D case to debug!)

I used the higher tolerances in order to decrease the computational time. The aim of this 2D case was actually to reproduce a first "quick and dirty" solution to see what problems might occur that could have major influence on the 3D case. With lower tolerances the stability is even worse.

This just indicates that the problem is not set up correctly. Lower tolerances on the linear solvers should lead to better convergence.

However, you said you are using meshes generated in ICEM. I recently had a similar problem with a 2D mesh in ICEM, but I did not investigate it too much, and I used another mesh/mesher. To try to understand if it's a mesh conversion problem, is it possible to use another mesh?

Best,

[robbirobocop]

Dear alberto,

by using another mesh that is a snappy grid in this case, the solver crashes earlier. So I rather used the finer grid (ICEM) "without" any errors. The mesh is very fine (160k cells) in order to lock out an error due to the mesh.

I also setup a small 3D case that only has 100k cells instead of 400k (the "normal" 3D case) to test around... Actually, something happens. But the inlet patches are not meshed appropriate since some cells look "bad". So I might test around with snappy some more time...

But I think something is wrong in the setup. Usually, I get problems with "k" rather than anything else. Hopefully, you will find what's wrong... Because compared to interFoam the computational time is like x1000 right now...I tried to act in accordance to the bubbleFoam description inside the OpenFOAM wiki. But it was only partly helpful.

I sent you the case along with a few notes a few minutes ago.
Thank you very much for your help.

[alberto]

Originally Posted by robbirobocop

Dear alberto,

by using another mesh that is a snappy grid in this case, the solver crashes earlier. So I rather used the finer grid (ICEM) "without" any errors. The mesh is very fine (160k cells) in order to lock out an error due to the mesh.

I also setup a small 3D case that only has 100k cells instead of 400k (the "normal" 3D case) to test around... Actually, something happens. But the inlet patches are not meshed appropriate since some cells look "bad". So I might test around with snappy some more time...

But I think something is wrong in the setup. Usually, I get problems with "k" rather than anything else. Hopefully, you will find what's wrong... Because compared to interFoam the computational time is like x1000 right now...I tried to act in accordance to the bubbleFoam description inside the OpenFOAM wiki. But it was only partly helpful.

I sent you the case along with a few notes a few minutes ago.
Thank you very much for your help.

The converted mesh from ICEM has boundary conditions set incorrectly. If you take a look at constant/boundary, all the BC's are walls or empty. Inlets and outlets should be of type "patch". You can change this manually in the boundary file, or setting the BC's in ICEM/FLUENT.

Two considerations: you are dealing with a system involving phase separation. The implemented k-epsilon model is not well defined in such a situation. I have a version of twoPhaseEulerFoam with a turbulence model based on the mixture assumption which I can share if you want. To check if the rest of the setup is OK, switch off the turbulence model.

Second consideration: If you are interested in resolving interfaces, you might want to consider multiphaseEulerFoam, which is a multi-fluid solver, with the capability of reconstructing interfaces too. The algorithm it implements seems to be more robust. The turbulence model is LES-based however (Smagorinsky + vanDriest damping).

Best,

[robbirobocop]

Originally Posted by alberto

The converted mesh from ICEM has boundary conditions set incorrectly. If you take a look at constant/boundary, all the BC's are walls or empty. Inlets and outlets should be of type "patch". You can change this manually in the boundary file, or setting the BC's in ICEM/FLUENT.

I did not know that this might be a problem. I always used "walls" for inlets and outlets as well. For example interFoam never had a problem with the boundary file setup like this...

Originally Posted by alberto

Two considerations: you are dealing with a system involving phase separation. The implemented k-epsilon model is not well defined in such a situation. I have a version of twoPhaseEulerFoam with a turbulence model based on the mixture assumption which I can share if you want. To check if the rest of the setup is OK, switch off the turbulence model.

It would be nice if you could share your turbulence model with a "little" explanation on how it has to be implemented and how to calculate the values (if there is a difference to the usual bubbleFoam turbulence parameters)...

Originally Posted by alberto

Second consideration: If you are interested in resolving interfaces, you might want to consider multiphaseEulerFoam, which is a multi-fluid solver, with the capability of reconstructing interfaces too. The algorithm it implements seems to be more robust. The turbulence model is LES-based however (Smagorinsky + vanDriest damping).

Best,

I will keep that mind for future consideration.
Thanks a lot for the information.

[robbirobocop]

Well, alberto.

Even when I change the boundaries from "wall" to "patch" nothing really happens at the in- or outlets. Could you check my BC's once again? Right now, something happens inside the drum. But the in- and out flows are not considered...

[robbirobocop]

Dear alberto,

finally the case works at least laminar...
I therefore set up a small 3D case that has a smaller mesh then the 2D mesh I sent you via e-mail... My problem with the turbulence is the turbulent kinetic energy k... Therefore, I can only use time steps around 1e-07 / 1e-06... Thus, even the simulation of the small 3D test case would take until 2013... That leads me to two questions...

What turbulence model - the one you'd like to share - have you implemented into bubbleFoam?
How should I set up k and epsilon to have a stable simulation? The OpenFOAM wiki tells me to set up k and epsilon for the continous phase... But with these values and initialisation I wouldn't get far to a result...

Thanks for your help.

Best,
Rob