Rewriting twoPhaseEulerFoam in conservative form
 

Hello,

I started to work on an improved version of a two-fluid solver with the goal of dealing with gas-particle flows with the traditional kinetic theory closures for the particle phase and, eventually, in future, more advanced techniques.

At first I thought to work on a two-phase code, because there still are a few questions to be addressed before adding complexity to the problem. So the steps I plan to follow are:

1. Write a two-fluid solver which does not rely on the phase-intensive formulation of the phase momentum equations. The phase intensive form is of great help in stabilizing the solution when the phase volume fraction tends to zero, but it introduces non-conservative terms, which might lead to problems in the solution under certain conditions (See for example Parks et al., Nuclear Engineering and Design, Vol 239, Issue 11, November 2009, Pages 2365-2371)
2. Implement a pressure-based compressible solver for the particle phase, to deal with the strongly non-linear behaviour of the particle pressure. After studying the existing algorithms for gas-particle flows, I think a pressure based solver might be the easiest approach to implement in OpenFOAM(r), compared to volume fraction based approaches. Moreover, looking at the functional form of the granular pressure in the dense limit, where the pressure diverges (the flow becomes incompressible), it becomes clear that solving for the granular pressure instead than for the volume fraction directly and then obtain the pressure from the equation of state is less troublesome.

Now, in doing this, there are some problems that need to be solved:

• Almost all the algorithms in the literature of multiphase flows are based on a staggered grid arrangement, while OpenFOAM uses colocated grids. This problem has been addressed already in twoPhaseEulerFoam, with Weller's pseudo-staggered approach, which relies on an appropriate use of the flux reconstruction to include certain contributions (gravity, drag, ...) in the momentum equation.
• The phase momentum equation is singular when the phase volume fraction becomes zero. This problem does not exist in the phase-intensive formulation, but it appears again when we consider the conservative form. A possible workaround is to define regions in the computational domain where the volume fraction alpha is smaller than a quantity alphaMin, and proceed to the solution only where alpha > alphaMin. To do this in OpenFOAM, one way is to define flags, and set the fvMaxtrix up so that it returns U = 0 where alpha < alphaMin. This seems to work correctly almost everywhere, but not at the interface between two phases. In a few cells accross the interface, (See MFIX numerics guide, www.mfix.org), the volume fraction of the dispersed phase is small, and the central coefficient becomes zero (or nearly zero) but the right-hand side of the momentum equation is not zero due to the presence of the fluid-phase pressure gradient and of the gravity term. This leads to a velocity peak at the interface that disturbs the solution significantly (it does not make the code crash in my implementation, and it does not cause oscillations to propagate in the solution, but the velocity peak is clearly wrong and negatively affects the volume fraction field for example). The MFIX numerical guide suggests to fix the problem solving for an approximate momentum balance in the interface cells, which are identified by looking for places where the central coefficient is very small, but I still need to figure out how to implement it "the OpenFOAM way" :). If you have suggestions to resolve this problem in some way, help is welcome. To give an idea of the problem, I attach a picture of the velocity profile along the axis of a channel with particles at the top that are falling under gravity (beginning of the run)

http://dl.dropbox.com/u/659842/VelocityPeak.png

Currently I have a rough implementation with a fictitious equation of state (I'm assuming the granular temperature to be constant), without the fix that causes the velocity peak at the interface.

P.S. Of course the code will be released, as soon as it works :D

Track this work
 

Hi Alberto,

This work sounds very attracting. How is it going on?

I was flooded by the code. twoPhaseEulerFoam is not easy to be adapted to new turbulence lib structure, especially for a roockie as me.

I look forward to your advance and wish all the best.

Thanks.

Roro

It is interesting!

I did some work to create a RANS multiphase turbulence library in the style of the single phase turbulence libraries. They generalised reasonably well, at least it makes writing new models easier. It would be interesting to compare how we've both done it.

Hi Laurence and Roro :-)

@Laurence: how is it going with your work? :)

This work is going on, a bit slowly due to other obligations, but it's not stuck. I can give a short status update, even if there still are some open questions. Please, keep in mind I'm interested in a code to compute gas-particle flows, so some of the developments cannot be directly generalized to a gas-liquid system.

In particular, I dropped the pressure equation based on total continuity, and I implemented two pressure equations, one for the fluid and the other for the particle phase.

Proper interpolation practises removed the peak at the interface (it's really all already in OpenFOAM, in some of the solvers)

The algorithm to enforce the particle packing limit is what gave me more troubles, and there still are some open question on this point to improve the stability of the solver, given the level of stiffness introduced by the equation of state.

Best,

Very slowly...

I've got a multi-fluid solver alongside DQMoM solving a coupled PBE implemented now (multiple dispersed phase momentum equations so that each node is transported with a different velocity) , and it seems to work ok. There are very slight differences between the DQMoM and the QMoM solution, but I suppose that is to be expected, and they should vanish with sufficient refinement.

I'm still having problems when solving the counter-current contactor (due to the phase inversion). I can get a solution which looks reasonable, but I have a lot of validation to do and it's not very stable. There was a big problem with some of the breakup models (at least, with the way I'd implemented them!), so I hope to have that cleared up soon.

The difference between QMOM and DQMOM is, I would say, OK. Btw, an improved flavour of multi-variate QMOM will come out soon :D

About twoPhaseEulerFoam stability, how does it crash? Alpha going above 1 or drag not converging?

Best,
Alberto

It's when the coupling between the CFD and the PBE is activated. The CFD obviously does not like having too big/small diameters being returned to the drag force. I can get around it by imposing cut offs but it's rather arbitrary.

Probably a more implicit treatment of the drag would help. The semi-implicit treatment is known to work in some situation, but in general the partial elimination algorithm (PEA) or a coupled solution of the momentum equations is required for stability purposes.

The PEA is relatively easy to put into OF, the coupled solution requires much more work, since no coupled solver is available in OpenFOAM at the moment. MFIX uses PEA for n-phases, and it is very robust though.

Best,

Thanks for the tip, I'll look into it.

It looks like Fluent uses a coupled solver, I should have noticed this before. :o

Yes, FLUENT uses a partially coupled solution approach for its n-fluid segregated solver. The find the estimate of the velocities in coupled manner, and then solve for the pressure equation, correcting the velocities.

It is surely an effective solution, but probably also one of the reasons why FLUENT multi-fluid simulations are not exactly fast, at least compared to what MFIX can do.

It is surely possible to use the partial elimination with two-phases, and you can derive an incomplete version of the PEA for n phases (See MFIX - www.mfix.org documentation and Mathiesen's work (Prof. Hjertager group). I did something similar in MFIX to couple QMOM for kinetic equations to the fluid solver, and there I've 8 velocities (quadrature abscissas) per each particle phase :D).
Drop me an e-mail if you want to talk about it. I have some notes and references on the topic.

Best,
Alberto

I've just written the PEA for a two-fluid case, it seems to give much better results, I'll post a derivation and comparison up soon after I've checked it over.

Yes, if you have small bubbles, with strong coupling, the difference is very clear :cool:

Alberto & Laurence

I am very interested to know about the progress in this area. As you migth know, my group has had success in using the multiphase PEA algorithmn in both bubble column and fluidised bed applications also using the DQMOM. We have our own fortran code but would be interested in using OF for thes problems.

regards
Bjorn

Dear Bjorn,

you have been one of my teachers at a summer school in Udine in 2006, when I was doing my PhD in Torino. I really appreciated the insights in the multiphase numerics you gave. They were of great help for my work. :)

My focus is on gas-particle flows, and what I did up to now is to modify twoPhaseEulerFoam so that:

• The PIMPLE algorithm is used, so that convergence of the equations at each time step is ensured
• The conservative form of the momentum equation is used. This seems to improve the quality of the solution in some quite extreme case (complete phase separation in gas-particle flows, with packed particle phase). However there still is one open question on this point, and it is about the management of the limit of the phase fraction going to zero. Right now I use a simple cutoff value, not solving for U below that value, however there might be better approaches (Hints welcome :))
• The momentum predictor is solved
• PEA can be used for the momentum coupling through the drag (a bit messy with the flux reconstruction, since the PEA terms permeates all the fluxes, and, as a consequence, all the main equations).
• The pseudo-staggered approach used in twoPhaseEulerFoam is slightly modified to remove oscillations at phase interfaces, which are more troublesome with the conservative form of the momentum equation (probably the cut-off approach I use makes things worse).

I think some relatively minor cleanup is still needed in the code for what concerns the kinetic theory/frictional models.

About QMOM, a colleague of mine in Dr Fox group, who will graduate next week has been working on a new flavour of QMOM for multi-variate distributions (Conditional QMOM - CQMOM) applied to nano-precipitation. She coupled it with a single-phase solver in OpenFOAM. At the moment, the implementation is not done "the OpenFOAM way", but the code is functional.

Some extra updates. QMOM has also been extended to deal with kinetic equation (Boltzmann-like, to be clear), with the goal of directly approximating the kinetic equation to describe poly-disperse flows and overcome the limitations of multi-fluid models. In this method QMOM-like techniques are used to solve for the dispersed phase velocity, without making restrictive assumptions to derive closures for the phase viscosity as in KTGF. This development is not done in OpenFOAM yet, since the main initial goal was to develop the approach without too many complications. However we have a functional code based on MFIX, and the procedure is co-located by nature, so porting it should not be too painful (we are thinking about it, since now the base multiphase code seems to be quite robust).

Best,

Dear Bjorn,

I see you are at Stavanger. I lived there and in Bergen until I was 7. I can't remember much apart from the snow though. :)

Yes, I've read the work of Stefano Bove in particular and enjoyed reading Neils Deen's work.

I've been working on multiphase DQMoM for gas-liquid/liquid-liquid applications, where each quadrature node has its own velocity and hence a separate momentum equation is solved for. I also want to increase the number of internal coordinates in the population balance to simulate mass transfer in a real application.

I've extended bubbleFoam to solve for N dispersed phases and coded multiphase PEA (hope it's correct - it does give equivalent results to the semiImplicit method with better residuals when the drag force is large). I also have some multiphase turbulence models and some surface forces implemented in a runtime selectable framework.

I have problems with solution stability when solving a momentum equation for each quadrature node. The residuals are very oscillatory and some cases diverge completely when coalescence and breakup in the population balances are on. I am not sure how to proceed from here: I may look at using a coupled approach. Any suggestions would be appreciated.

Alberto: Regarding the cut off for solving the momentum equations. We had convergence problems when imposing arbitrary cut-offs for reactive flows using DQMoM, and found it better to filter/smooth the source terms towards zero. I'm not sure if it would be relevant (some limited discretisation schemes probably do a similar thing?) but the paper is here.

Dear Laurence,

We also had some problems in solving the PBM using DQMOM in Fluent. You migth look at the PhD thesis of Jesper Madsen where he explains what he did:
http://hugin.aue.auc.dk/publ.php?id=2007

Regards
Bjorn

Hello,

So you found an operator splitting was necessary?

I used that for the reacting flow 'version' of DQMoM, where the chemical time scales are much smaller than that of the flow, we needed to have an operator splitting.

Interesting, this should be quick to try. Did you use RADAU5 or CVODE for the time integration?

Regards,

We do something similar to what Bjorn told you. Split and use an ODE (if necessary, stiff) solver for the source term part.

About my problem with the cuttoff, I do not have instabilities related to it, but a certain sensitivity to the cut-off value. In other words the "small" value used as cutoff value cannot be too small (it is typically between 10^-4 and 10^-6). In some cases 10^-6, which is what I would like to use, is too small and gives some trouble at interfaces (wrinkles), while 10^-4 works fine, but it might be too high in some case, depending on the application.

Thank you for your reference. I'll take a look and let you know :)

Hello there.

I've been working with multiphase simulations for a while now.

One of the new goals of the taskforce I work with is the evaluate OpenFOAM as a possible substitute ofr other simualtions packages.

Would one og guys mind telling me how does the code used in the tutorial example "DamBreak4phase" ensures that the sum of all components fractions will be equal to 1???

After reading your discussion I thought one of you could help me.

Thanks a lot

Hi, the limitation is included directly in the equation for alpha, following the same approach used in interFoam, where the convective flux is re-written in terms of the mixture and relative velocity.

You can find the implementation of this into multiphaseInterFoam in the solveAlphas method in multiphaseMixture.C.

Best,