[GerhardHolzinger]

I have a bubble column which simulates just fine with twoPhaseEulerFoam-2.3.0. However, twoPhaseEulerFoam-2.3.x breaks the case.

I attached some images to illustrate my problem.

With OF-2.3.0 all turns out just fine. The screenshots of this case are at t=3s.

With OF-2.3.x, however, the simulation crashes at some point after t=0.2s.

I attached also the case for you to play around.


So, to repeat my problem or question:

How can I get my case running in twoPhaseEulerFoam-2.3.x with the fully conservative formulation (i.e. fvm::ddt(alpha1, rho1, U1) instead of fvm::ddt(alpha1, U1))?

[cutter]

Hi,

I'm currently investigating twoPhaseEulerFoam too. I have no idea what the reason for the crash is, but I can confirm that it's also crashing after 0.1301399s with my version of OF 2.3.x.

Just an idea: You could use 'git bisect' to find the patch that is causing the problem. This will take some time (depending on the number of git commits and your compile times), but it might help you finding the problem or maybe filing a bug report.

BTW: nice talk at the conference!

Cutter


http://git-scm.com/book/en/v2/Git-To...gging-with-Git

[alberto]

The new formulation couples the system of equations in a much tighter way. You must ensure that the whole system converges well.
The modification was introduced to ensure the conservation of the disperse phase, and to guarantee energy conservation when there are large gradients of the dispersed phase fraction.

In short, use multiple (20 - 40) outer correctors, and a sufficiently strict Courant condition.

[GerhardHolzinger]

Thanks for your reply and your note at my bug report (http://www.openfoam.org/mantisbt/view.php?id=1441#c3281).

Sorry, for being a whining baby, but is this the only way to get the cases to run?

Making the time step small by prescribing a small Co number and performing tens of outer iteration causes the case to take ages.

I know generalisation is a good thing, but I don't want my isothermal, incompressible cases to take x times as long just because I use a newer version of OpenFOAM.

Although I like the features introduced in OF-2.3.x such as the fvOptions mechanism in twoPhaseEulerFoam, I highly dislike the pain it causes me to get any case running.

[GerhardHolzinger]

I still don't manage to get the case running.

I use maxCo = 0.2 and the following residualControls.

To get to 0.1s simulated time, I had to wait for 2h and the result is still the same unphysical mess.

Does anybody get this bubble column to work?

PIMPLE
{
    nOuterCorrectors 40;
    nCorrectors      3;
    nNonOrthogonalCorrectors 0;
    
    residualControl
    {
        "(U|k|epsilon)"
        {
            relTol          0.0001;
            tolerance       0.0001;
        }
    }

    turbOnFinalIterOnly off;

}

[alberto]

I'm giving it a try now :-)

[alberto]

One question: is the inlet made by a central square and a small square in a corner?

That is what topoSet + createPatch seems to generate here.

[alberto]

Ok. Not sure how that happened. But I answered my question... just a square in the center.

[alberto]

Originally Posted by GerhardHolzinger

Although I like the features introduced in OF-2.3.x such as the fvOptions mechanism in twoPhaseEulerFoam, I highly dislike the pain it causes me to get any case running.

2.3.x is a development version. I would stick to 2.3.0 for production, if you don't need the features in 2.3.x.

Also, your case has some inconsistency which are independent from the problem you are encountering: you fill the column with water, water cannot go out (slip condition), but you add mass to the system. Also, Uwater at the inlet should be zero, being alpha.water zero.

[alberto]

I investigated a bit more, and run your simulation with Schiller-Naumann drag, obtaining reasonable results. Could you check that? Maybe it is a problem in the implementation of the drag law.

[GerhardHolzinger]

Thank you for having a look on my problem. I apologize for posting a messy case, that's clearly my homework to do.

This bubble column is based on the thesis and other publications of N. Deen, e.g. [1,2]


The problem with the water not being able to leave the domain is a tricky one. How could I model the upper boundary without including the free surface into my simulation domain?
Including the free surface introduces additional problems, e.g. phase-inversion.


I cleaned up and attached the case.

[1] N. G. Deen, T. Solberg, and B. H. Hjertager. Large eddy simulation of the gas-liquid flow in
a square cross-sectioned bubble column. Chemical Engineering Science, 56:6341–6349, 2001.

[2] Niels G. Deen. An experimental and computational study of fluid dynamics in gas-liquid
chemical reactors. PhD thesis, Aalborg University Esbjerg, 2001.

[alberto]

Thanks for the updated case. I will look into it.

Originally Posted by GerhardHolzinger

The problem with the water not being able to leave the domain is a tricky one. How could I model the upper boundary without including the free surface into my simulation domain?
Including the free surface introduces additional problems, e.g. phase-inversion.

The problem of phase inversion is well addressed in twoPhaseEulerFoam, so I would not worry about it. For stability, it is probably a good idea to actually include the free surface (it won't be resolved in detail, so it won't add computational cost) because it allows recirculation of water at the boundary to be avoided. Additionally, it lets you specify an easier boundary condition (constant pressure).

In any case, I will come back with additional comments once I look at the updated case :-)

[alberto]

I have come to the conclusion that the major culprit is the TomiyamaAnalytic drag. I ran the case with both Schiller-Naumann and Tomiyama-Kataoka-Zun-Sakaguchi drag law (I implemented it), and I don't see the very strange velocity profile that comes out with the TomiyamaAnalytic implementation.

Is the grid you are using here the same you had when you ran with 2.3.0 and obtained good results, or did you coarsen it? I don't see a lot of fluctuations in the bubble plume with the Schiller-Naumann drag. The case with Tomiyama-Kataoka-Zun-Sakaguchi is still running, so I will know soon.

[GerhardHolzinger]

Hi,

it is the same the mesh. I should definitely try other drag laws.

I will look into this issue.

Thanks for taking so many looks on my stuff.

It was definitely good to have a different pair of eyes on a problem I am stuck with.

[alberto]

Originally Posted by GerhardHolzinger

Hi,

it is the same the mesh. I should definitely try other drag laws.

I will look into this issue.

Thanks for taking so many looks on my stuff.

It was definitely good to have a different pair of eyes on a problem I am stuck with.

Any improvements?

[demichie]

Originally Posted by GerhardHolzinger

I have a bubble column which simulates just fine with twoPhaseEulerFoam-2.3.0. However, twoPhaseEulerFoam-2.3.x breaks the case.

I attached some images to illustrate my problem.

With OF-2.3.0 all turns out just fine. The screenshots of this case are at t=3s.

With OF-2.3.x, however, the simulation crashes at some point after t=0.2s.

I attached also the case for you to play around.


So, to repeat my problem or question:

How can I get my case running in twoPhaseEulerFoam-2.3.x with the fully conservative formulation (i.e. fvm::ddt(alpha1, rho1, U1) instead of fvm::ddt(alpha1, U1))?

I also have some problem with twoPhaseEulerFoam that seem related with the drag. In particular, I have noticed that in version 2.2 there was an equation for u1 where the drag term was formulated as

alfa2/rho1*k*(u2-u1).

In version 2.3.0 the equation for the momentum is written for alfa1*u1 and the drag term is:

k/rho1*(u2-u1).

In version 2.3.x the equation for the momentum is written for alfa1*rho1*u1 (conservative form) and the drag term is:
k*(u2-u1).

The terms k are different in the different versions of OpenFoam, but it seems to me that they are not consistent.
In particular, looking at 2.3.0, k is defined as

0.75*Cd*Re*mu_c/(alfa_c *d_p^2),

where the index f and p are for the continuous and dispersed phases.

In 2.3.x k is defined as

0.75 * Cd*Re*alfa_p*mu_c/( d_p^2).

So, there is a difference in the way the the volume fractions are considered in the drags, and to me the right formulation should be:

0.75 * Cd*Re*alfa_c*alfa_p*mu_c/( d_p^2).

In addition, it is not clear how the blending for he drag coefficients is implemented. In BlendedInterfacialModel.C there are two coefficients f1 and f2, functions of the volumetric fractions of the phases and the values given in the dictionary phaseproperties. From what I have understood when one phase is fully dispersed we should have f1=0 and f2=1, or vice versa.
But then the two drags are added in the following way:

x() += model1In2_->K()*(1 - f1);
x() += model2In1_->K()*f2;

So it seems that the drags evaluated for the two cases (1 dispersed in 2 and 2 dispersed in 1) are summed in any case, because when f1=0 and f2=1 =>(1-f1)=1 and f2=1. Is this correct?

Thank you for any help!
Mattia

[GerhardHolzinger]

Hi,

you can find here [1] some of my thoughts and findings on twoPhaseEulerFoam.

In OF-2.2 you could choose the drag phase to be "1", "2" or blended. In OF-2.3 the blending does this job.

You are correct with your summary. Since, twoPhaseEulerFoam-2.3 is formulated in a very general way, you have phase 1 dispersed in phase 2 and vice versa.




[1] https://github.com/OSCCAR-PFM/OSCCAR-DOC-PUBLIC

[demichie]

Thank you for the link to your document. I just want to point that both 2.3 and 2.3.x should use the same blending where both the phases can be the continuous phase, but there are some differences in the way they consider the volume fractions in the drag term.


This is the momentum equation in 2.3.0:

U1Eqn =
(
fvm::ddt(alpha1, U1)
+ .... ==
- fvm::Sp(dragCoeff/rho1, U1)

with the dragCoeff term evaluated here:

0.75
*CdRe()
*swarmCorrection_->Cs()
*pair_.continuous().rho()
*pair_.continuous().nu()
/(
max(pair_.continuous(), residualAlpha_)
*sqr(pair_.dispersed().d())
);

and this is the formulation in 2.3.x:


U1Eqn =
(
fvm::ddt(alpha1, rho1, U1) + fvm::div(alphaRhoPhi1, U1)
.... ==
);
U1Eqn.relax();
U1Eqn += fvm::Sp(dragCoeff, U1);

with its dragCoeff:

0.75
*CdRe()
*max(pair_.dispersed(), residualAlpha_)
*swarmCorrection_->Cs()
*pair_.continuous().rho()
*pair_.continuous().nu()
/sqr(pair_.dispersed().d());

Please note that the difference in the difference in the momentum equation is that in 2.3 there is a non-conservative form, without the density in the derivative with respect to time, and for this reason dragCoeff is divided by rho1 in the drag term.
But in version 2.3.x dragCoeff is defined multiplying for the volume fraction of the dispersed phase ( *max(pair_.dispersed(), residualAlpha_) ), while in 2.3 the volume fraction is at the denominator:

/(
max(pair_.continuous(), residualAlpha_)
*sqr(pair_.dispersed().d())
);


To me these problems appeared when there was the change from the original formulation in 2.2 where the equation was written for U, and the momentum equation was divided by both the density and the volume fraction to obtain the desired formulation. But when we write the equation for the conservative form both the formulations used in 2.3 and 2.3.x are wrong, because in dragCoeff there should be both the volume fractions at the numerator. This is for example what I think is done in Fluent:

http://aerojet.engr.ucdavis.edu/flue...heory-exchange

Ciao
Mattia


P.S. Unfortunately the presence of the volume fractions at the numerator or denominator does not changes the dimensions of the drag term, so it is not possible to check the correctness of the term by looking at the dimensions only.

[demichie]

Hi all,
I have attached a .pdf describing how the drag terms K_1in2 and K_2in1 contribute to the total drag term K. In addition to there terms there is also the drag computed when neither phase is truly dispersed. Here there is the code taking into account this "segregated drag" (model_->K()):


if (model_.valid())
{
x() += model_->K()*(f1() - f2());
}

My problem is that the term f1-f2 can be negative...
I think instead for the segregated case the following formulation should be used (always positive, null when one phase is fully dispersed, and maximum when none of the phases is dispersed):


if (model_.valid())
{
Info<< "FIRST f1 min.internalField()" << endl;
x() += model_->K()*(1 - f2()) * f1();
}

Ciao
Mattia

[alberto]

Originally Posted by demichie

The terms k are different in the different versions of OpenFoam, but it seems to me that they are not consistent.
In particular, looking at 2.3.0, k is defined as

0.75*Cd*Re*mu_c/(alfa_c *d_p^2),

where the index f and p are for the continuous and dispersed phases.

In 2.3.x k is defined as

0.75 * Cd*Re*alfa_p*mu_c/( d_p^2).

So, there is a difference in the way the the volume fractions are considered in the drags, and to me the right formulation should be:

0.75 * Cd*Re*alfa_c*alfa_p*mu_c/( d_p^2).

If you look at the reliable literature on two-fluid models (Simonin, Troshko, Lahey, Enwald, ...), you will see that only the disperse phase volume fraction is included in the definition of K:

K = 3/4 * alpha_disperse * rho_continuous * C_D * Ur / dp

You find this formulation also in CFX. Some authors (Wen and Yu, for example) modify C_D and Re to include the effect of the presence of multiple particles, but K must remain unchanged, or you would have an incorrect form of some of the drag laws (Schiller and Naumann, for example, does not have a dependency on alpha_continuous anywhere).

If you carefully examine for example the implementation of the Wen and Yu drag model in OpenFOAM, you will see that the dependency on alpha_continuous is introduced correctly when needed (they multiply by alpha_continuous^-2.65). As a consequence, adding alpha_continuous in the function that calculates K would not be correct.

I hope this helps.

Best,