[Aurelien Thinat]

Hi everybody,

I need some advice in adding a source term for the dispersed phase in the bubbleFoam solver.

At first I put a volscalarfield named "source" in my domain and I transformed the alpha equation :

fvScalarMatrix alphaEqn
(
fvm::ddt(alpha)
+ fvm::div(phi,alpha,scheme)
+ fvm::div(-fvc::flux(-phir,beta,scheme),alpha,scheme)
=
source
)

But I had a loss of "beta" in my calculations :
A column of water (z<z0)with a volumic source term of air bubbles at the bottom. The top of the column is filled with air (z>z0).
After a few time, the water column collapses on itself.

I guess that's because of the definition beta = scalar(1) - alpha. If I create some alpha, I delete some beta right ?

So I would need some way to add this source term of alpha keeping the same mass of beta in the whole domain. Is that possible without rewriting a new solver ?

Thanks.

[Aurelien Thinat]

Ok, I got the answer.

I have modified the alpha equation the right way, but I didn't do the same on the pressure equation...

Thanks Cyp.

[Aurelien Thinat]

Hello everybody,

I'm still missing something.

I have modified my phase continuity equation as written above :
fvScalarMatrix alphaEqn
(
fvm::ddt(alpha)
+ fvm::div(phi,alpha,scheme)
+ fvm::div(-fvc::flux(-phir,beta,scheme),alpha,scheme)
=
source
)

So I need to modify my pressure equation right ?
To obtain the pressure equation, one should force div(phi)=0 in the case with no source term, which becomes div(phi) = source in my case, am I still right ?
Then the final pressure equation equation should be :
pEqn
(
fvm::laplacian(Dp,p) == fvc::div(phi)-source
)

EDIT :
Well in fact I get lost in Henrik Rusche Thesis at p.112 :
He wrote that the mixture continuity equation is the equation 3.39 by combining the 2 phase continuity equation with weightning coefficients wi.
In my case the 0 should become "rhoa*source/wa".

But then he decides to solve the volumetric mixture continuity equation which is div(U)=0.
Since he got this equation by combining the 2 phase continuity equations, I guess I still should write :
div(U)= div(alpha*Ua + beta*Ub)=source
And that leads to the pressure equation :
pEqn
(
fvm::laplacian(Dp,p) == fvc::div(phi)-source
)
But when I code this, I have an infinite creation of water in my domain : the water level is rising and the mean volumic ratio of the dispersed phase is decreasing...

So I'm missing something but I don't know what. If someone has any clue...Thanks.

[Aurelien Thinat]

Hi again,

I continue on my problem. After few tries, I got that I need a coefficient next to my source term in the pressure equation : I have test manually few coefficients for a given value of "source" and I found one that allow mass conservation of the dispersed/continue phase.


In his thesis, Henrik Rusche wrote :
"we chose to solve for the mixture pressure by recasting the volumetric mixture continuity equation (3.3) into one for pressure. To achieve this, it is formulated at the cell faces, such that:
div(alpha_af * phi_a + alpha_bf * phi_b) = 0 "

So, the pressure equation is written at cell's faces only. Does it mean I have to use some kind of "phi_source", a projection of the source term on each face ?


EDIT :

I'm trying to interpolate my "volScalarField source" at the cell's face. The volScalarField alpha is interpolate at the cell's face by the command :
"surfaceScalarField alphaf(fvc::interpolate(alpha));"

So I defined :
"surfaceScalarField sourcef(fvc::interpolate(source));"


But then, when I try to include sourcef in the pEqn equation, I have an issue about their type with the operator "+" :

pEqn.H:54: error: no match for 'operator+' in
'Foam::fvc::div(const Foam::GeometricField<Type, Foam::fvsPatchField, Foam::surfaceMesh>&) [with Type = double]()
+
Foam::fvc::interpolate(const Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh>&) [with Type = double]()'


If anyone has an idea...

[Aurelien Thinat]

No one has ever try to add a source term of alpha in this solver ?

A little summary (and a up of the post) :

1) fvc::div(phi) is a volScalarField,
2) fvm::laplacian(Dp,p) is a volScalarField too,
3) The initial pressure equation is

fvm::laplacian(Dp,p)=fvc::div(phi)

as the continuity equation at the face center is :

div(alpha_f*phia + beta_f*phib)=0

4) I add a source term in alphaEqn :

fvm::ddt(alpha)

+ fvm::div(phi,alpha,scheme)

+ fvm::div(-fvc::flux(-phir,beta,scheme),alpha,scheme)

==
source

5) My continuity equation at the face is so :

div(alpha_f*phia + beta_f*phib)=source_f (if anyone can confirm I am right)

I defined sourcef as fvc::interpolate(source). It is a surfaceScalarField.

6) How can I add it in the fvScalarMatrix pEqn ?

[alberto]

Hi,

consider the two continuity equations:

ddt(alpha) + div(alpha*Ua) = Sa
ddt(beta) + div(beta*Ub) = 0

Since they are volumetric already (we assumed constant density, so the equations were divided by rhoa and rhob respectively), there is no need of weighting again.

Sum the equations to obtain:

div(alpha*Ua) + div(beta*Ub) = Sa

Replace the velocities with face fluxes, and you have what you are looking for.

Best,

[Aurelien Thinat]

Thank you for answering Alberto.

So if I understand what you are saying I will have :
div(alpha_f*phia + beta_f*phib) - laplacian(Dp,p)=Sa (I take your notation)

That's the equation I have coded at first but with this one, the mean volume fraction in the whole domain is not stable. It leads to an endless creation of continuous phase.

[alberto]

Originally Posted by Aurelien Thinat

Thank you for answering Alberto.

So if I understand what you are saying I will have :
div(alpha_f*phia + beta_f*phib) - laplacian(Dp,p)=Sa (I take your notation)

That's the equation I have coded at first but with this one, the mean volume fraction in the whole domain is not stable. It leads to an endless creation of continuous phase.

That's the equation that comes out from what you said. However the formulation of the model seems a bit inconsistent. If the mass leaves phase a, where does it go? In deriving the pressure equation you are imposing that the total mass (volume, since rho = const) is conserved.

[Aurelien Thinat]

I'm trying to simulate the case of bubble plume. So I have a sharp interface between the dispersed phase a and the continuous phase b.
At t= 0 :
- In the top part of the domain I have alpha = 1,
- In the bottom part, alpha = 0.
The dispersed phase a can go through an outlet at the top of the domain. The creation should be ensured by the source term I'm trying to implement (I'd like to simulate an air bubble creation from a chemical reaction).

Normally, my interface should be moved to the top and my volume fraction a should be constant after few seconds (well, less than seconds).

[alberto]

Then I think your model is missing a term. If you transfer mass from a to b, then b has a sink of mass in its continuity equation.

Also, keep in mind that transferring mass, you transfer momentum. As a consequence you will have also a source term in the momentum equations.

[Aurelien Thinat]

I have edited the wrong post :

"I (and you) may have solved my problem. I wasn't confident in my pressure equation but it seems the problem came from the boundary condition in pressure I fixed : buoyantPressure, which I'm used for VOF models.
I have defined a zeroGradient boundary condition to the walls and the volume fraction is now constant."

Thank you again for your time.

Best regards,

Aurélien

[alberto]

Hi,

buoyantPressure is identical to zeroGradient if p is used instead of p_rgh, so it is a bit strange it was the source of your problems.

Best,
Alberto

[Aurelien Thinat]

Hi everyone,

Me again. But with another "source term" problem. I have to implement a turbulent dispersion force in Euler/Euler solvers. So let's start with bubbleFoam.

Starting with Lopez de Bertodano Model :
Force of turbulent dispersion, phase a = Ftd,a = -Ftd,b = C * rhob * k * grad(alpha)
with k the kinetic energy of the continuous phase.

Since the momemtum equations are divided by rhoa * alpha (rhob * beta for the second one), I have got :

Mtd,a = C * rhob * k * grad(alpha) / ( rhoa * alpha)

So I'm adding this term in the UaEqn and it's opposite in UbEqn (divided by rhob * beta instead). And...at the first pressure equation it blows up because of some rules in the operator "divide".

Am I missing something ?

Edit :

This is the UaEqn :

UaEqn =
(
(scalar(1) + Cvm*rhob*beta/rhoa)*
(
fvm::ddt(Ua)
+ fvm::div(phia,Ua,"div(phia,Ua)")
- fvm::Sp(fvc::div(phia),Ua)
)
- fvm::laplacian(nuEffa,Ua)
+fvm::Sp(fvc::div(phiRa),Ua)
+ (fvc::grad(alpha)/(fvc::average(alpha)+scalar(0.001)) & Rca)
==
-fvm::Sp(beta/rhoa*dragCoef,Ua)
- beta/rhoa*(liftCoef - Cvm*rhob*DDtUb)
-rhob*k*fvc::grad(alpha)/(alpha*rhoa) //this is the turbulent dispersion force
);

This is the error message :

#0 Foam::error:rintStack(Foam::Ostream&) in "......../libOpenFOAM.so"
#1 Foam::sigFpe::sigHandler..../libOpenFOAM.so
#2 ?? in "lib64/libc.so.6"
#3 void Foam::divide<Foam::Vector<double>,Foam::fvPatchFie ld,Foam::volMesh>
(........)

[alberto]

Hi Aurelien,

if your term is

Mtd,a = C * rhob * k * grad(alpha)

you have to add

C*rhob/rhoa*k*grad(alpha)/(alpha + 0.0001) to UaEqn

and

C*k*grad(alpha)/(beta + 0.0001)

since alpha and beta can be zero.

If this term is large, it's worth moving it to the flux term.

Best,

[Aurelien Thinat]

Thank you Alberto for your answer.

It's correct, alpha = 0 was the problem. For information, I'm testing right now this implementation :

UaEqn : C*rhob/rhoa*k*fvc::grad(alpha)/(fvc::average(alpha) + scalar(0.0001))
UbEqn : C*k*fvc::grad(alpha)/(fvc::average(beta) + scalar(0.0001))

Next step is the Simonin's model and working on a degassing boundary condition.

Thank you again.