[elBombastico]

Dear Foamers,

I am modelling the simultaneous transport of several passive scalars over a phase field. Because of this, the advection term of the transport equation reads as follows:

,


where alpha is the phase field (scalar), U is the velocity field (vector) and Psi is the quantity being transported (scalar). Since both alpha and Psi vary on the space, it seems to me that this term should be written as:

.


Since the flow is incompressible, I have implemented it in openFoam the following way:

 (Psi * U & fvc::grad(alpha)) 
  + alpha * Psi * fvc::div(phi)
   + (alpha * U & fvc::grad(Psi))

For the term that requires the divergence of the velocity field, I have used fvc::div(phi), where phi is the product of the velocity times the area at the cell faces, since, as far as I know:




where u_i would be the velocity through each face and s_i, the surface of the face.

So far, it compiles and runs without complaining. Does anybody know if this would be correct and/or if there is a more efficient way of doing it?

Regards and thanks in advance

[Tobermory]

The problem that I see with you expanding out the convection term into your three terms is that the three terms are not in a conservative form. This means that when you integrate them over a cell volume, as is required for the finite volume approach, you don't get simple expressions based on the average cell value or the average face value ... By contrast, with the single convection term, you can convert the divergence into a sum of fluxes on the cell face, with face averaged values of psi, U and alpha ... and this is second order accurate if you assume a linear variation of values across a face/cell (which is the underpinning of the FV approach).

Indeed, if you expand out the Taylor series for your approach, you'll probably find that your version includes a bunch of additional cross terms and is not equivalent to the "correct" implementation using the single term.

So in short - whilst you can do it this way (in the sense that it compiles), I wouldn't recommend it. I hope that helps.

[elBombastico]

Hi,


Thanks for your response. As you correctly guessed, the implementation with the expansion does not work - it blows up after a few time steps. After some tests, the reason is most likely the advective term, so I dropped the expansion.
First, I calculate the face fluxes of U times alpha the following way:

surfaceScalarField uAlpha
(
    IOobject
    (
        "uAlpha",
        runTime.timeName(),
        mesh,
        IOobject::READ_IF_PRESENT,
        IOobject::AUTO_WRITE
    ),
    linearInterpolate(alpha*U) & mesh.Sf()
);

and then, the divergence term as:

fvm::div(uAlpha,Psi)

But it blows up as well....
Do you happen to have any other suggestion?

Thanks in advance

[Tobermory]

I am not very conversant with VoF type solvers (for example - should there not be a density term in the divergence as well? probably my ignorance ...), but could you not simply multiply alpha and psi since they are both volScalarFields, and then just write a conventional divergence as div(phi,alphaPsi) or am I missing something?

[elBombastico]

Hi,


I have now tried multiple approaches, with different, yet disappointing results...

1. With the approach from my previous post, the divergence is calculated as

    uAlpha = linearInterpolate(alpha*U) & mesh.Sf()

   , and then

   fvm::div(uAlpha,Psi)

   . The simulation blows up, unless I use a bounded discretization scheme for the advective term. With "bounded Gauss upwind", the simulation converges, but the scalar leaks through the interface, which is not supposed to. I am now trying to determine whether this is an issue caused by numerical diffusion of the upwind scheme, a mesh that is too coarse, or by the velocity or phase fields not being appropriately resolved. Any other discretization scheme causes the simulation to blow out, or the scalar values to become unbounded beyond acceptable.
2. Your suggestion, which consists on creating a new volScalarField alphaPsi = alpha * Psi requires to specify the boundary conditions, as well as the initial values of the new volScalarField alphaPsi, or otherwise openFoam complains about the calculated BC not being applicable to that variable. I tried to do that, and resulted in the simulation blowing out as well .

It appears I am out of luck for today. I will try with a finer mesh and see what happens.
Thanks anyway for your suggestion. It was worth a shot