[randolph]

Hi everyone,

I am interested in simulating a scalar transport within the primary phase (i.e., water) in the interFoam.

I follow the approach of phaseScalarTransport (https://github.com/OpenFOAM/OpenFOAM...larTransport.H) and successfully added a scalar transport to the primary phase with MULES compressed flux (i.e., alphaPhiUn).

    fvScalarMatrix YEqn
    (
          fvm::ddt(Y)
        + fvm::div(alphaPhiUn, Y)
        - fvm::laplacian(Deff, Y)
        ==
        alpha1*fvOptions(Y)
    );

    YEqn.relax();
    fvOptions.constrain(YEqn);
    YEqn.solve();

The code runs well and as expected the scalar is transported only within the water phase. When the water surface arises, the scalar will arise too. However, when the water surface falls, the part of the scalar will "stuck" above the surface (see attachment). I would like to have scalar fall as the water falls. Could you please shed some light on this issue?

Thanks,
Rdf

[randolph]

Okay. I think I figured it out. I will post my findings here for future reference. There are several post concerns about this issue and provide valuable information.

InterFoam: Add an equation to only one phase
interFoam with transport: add a transport equation to only one phase
Numerical problem in species Transport in InterFoam at the interface
interFoam + a scalar transport equation
New phaseScalarTransport function object for a single phase in interFoam



In summary, the simplest way (at this moment) to simulate a scalar transport in a single phase (e.g., water) with interFoam is using the "phaseScalarTransport" utility from the OpenFOAM 7 or 8 (org version). The field alphaS.water = alpha.water * S.water should be used as the results, not the scalar itself (S.water). This is because the continuity equation is calculated for alphaS.water, as shown in the following code. If one looks at the S.water field, the scalar itself will be dispersed into the other phase. This is not a bug.

I use the dambreak as a testing case and found the mass of alphaS.water is generally conserved with an error of less than 1% while the mass of the scalar itself (S.water) will have an error range from 100% to 400%.

 
             fvScalarMatrix fieldEqn
             (
                 fvm::ddt(alpha, s_)
               + fvm::div(alphaPhi, s_, divScheme)
               - fvm::laplacian
                 (
                     fvc::interpolate(alpha)*fvc::interpolate(D),
                     s_,
                     laplacianScheme
                 )
              ==
                 fvOptions_(alpha, s_)
               - fvm::ddt(residualAlpha_, s_)
               + fvc::ddt(residualAlpha_, s_)
             );
 
             fieldEqn.relax(relaxCoeff);
             fvOptions_.constrain(fieldEqn);
             fieldEqn.solve(schemesField_);

In the ESI-OpenFOAM, although there is a phase option in the "scalarTransport" utility and they even provide tutorials (e.g, angledDuct, waterChannel), it does not seems that phase constrained transport is implemented at all. In these tutorials, the scalars only diffusive but not convective, as observed in this post (Foam::error::printStack(Foam::Ostream&) with interFoam -parallel). This is likely because the default option for phasePhiCompressed, alphaPhiUn is zero. Although alphaPhiUn is created in the interFoam, the interFoam solver does not pass a value to alphaPhiUn. This seems quite odd to me and please correct me if I am wrong.

There is a method suggested for convecting the scalar in the "scalarTransport" utility of ESI OpenFOAM. Post (Foam::error::printStack(Foam::Ostream&) with interFoam -parallel) suggest that defined the "alphaPhi0.water" for the "phasePhiCompressed". However, after some tests, this method will let some scalar stuck in the air. The "scalarTransport" utility of ESI OpenFOAM tries to conserve the mass for s.water instead of alphaS.water, as shown in the following code.

 
        tmp<surfaceScalarField> tTPhiUD;
        for (label i = 0; i <= nCorr_; i++)
        {
            fvScalarMatrix sEqn
            (
                fvm::ddt(s)
              + fvm::div(limitedPhiAlpha, s, divScheme)
              - fvm::laplacian(D, s, laplacianScheme)
              ==
                alpha*fvOptions_(s)
            );

            sEqn.relax(relaxCoeff);
            fvOptions_.constrain(sEqn);
            sEqn.solve(mesh_.solverDict(schemesField_));

            tTPhiUD = sEqn.flux();
        }

I hope this information helps,

Rdf

[Sud_fluid]

Hi,


I'm also trying to simulate a case (similar to sand in water + air) however I do not understand how to implement the phaseTransportFoam (in particular could you please explain what you meant by interfoam + phaseTransportFoam).

I tried to implement this in the dambreak case. In the controlDict, file I included the following lines before the end.

phaseScalarTransport1
{
type phaseScalarTransport;
libs ("libsolverFunctionObjects.so");
field s.water;
p p_rgh;
}

phaseScalarTransport2
{
type phaseScalarTransport;
libs ("libsolverFunctionObjects.so");
field alpha.water * s.water;
p p_rgh;
}
Then I run the command "interfoam" which ends well but when I enter the command "phaseTransportFoam" it is not found. Could you please tell the steps clearly on how to implement phaseTransportFoam.


Thanks.
Sudarshan.

Originally Posted by randolph

Okay. I think I figured it out. I will post my findings here for future reference. There are several post concerns about this issue and provide valuable information.

InterFoam: Add an equation to only one phase
interFoam with transport: add a transport equation to only one phase
Numerical problem in species Transport in InterFoam at the interface
interFoam + a scalar transport equation
New phaseScalarTransport function object for a single phase in interFoam



In summary, the simplest way (at this moment) to simulate a scalar transport in a single phase (e.g., water) with interFoam is using the "phaseScalarTransport" utility from the OpenFOAM 7 or 8 (org version). The field alphaS.water = alpha.water * S.water should be used as the results, not the scalar itself (S.water). This is because the continuity equation is calculated for alphaS.water, as shown in the following code. If one looks at the S.water field, the scalar itself will be dispersed into the other phase. This is not a bug.

I use the dambreak as a testing case and found the mass of alphaS.water is generally conserved with an error of less than 1% while the mass of the scalar itself (S.water) will have an error range from 100% to 400%.

 
             fvScalarMatrix fieldEqn
             (
                 fvm::ddt(alpha, s_)
               + fvm::div(alphaPhi, s_, divScheme)
               - fvm::laplacian
                 (
                     fvc::interpolate(alpha)*fvc::interpolate(D),
                     s_,
                     laplacianScheme
                 )
              ==
                 fvOptions_(alpha, s_)
               - fvm::ddt(residualAlpha_, s_)
               + fvc::ddt(residualAlpha_, s_)
             );
 
             fieldEqn.relax(relaxCoeff);
             fvOptions_.constrain(fieldEqn);
             fieldEqn.solve(schemesField_);

In the ESI-OpenFOAM, although there is a phase option in the "scalarTransport" utility and they even provide tutorials (e.g, angledDuct, waterChannel), it does not seems that phase constrained transport is implemented at all. In these tutorials, the scalars only diffusive but not convective, as observed in this post (Foam::error::printStack(Foam::Ostream&) with interFoam -parallel). This is likely because the default option for phasePhiCompressed, alphaPhiUn is zero. Although alphaPhiUn is created in the interFoam, the interFoam solver does not pass a value to alphaPhiUn. This seems quite odd to me and please correct me if I am wrong.

There is a method suggested for convecting the scalar in the "scalarTransport" utility of ESI OpenFOAM. Post (Foam::error::printStack(Foam::Ostream&) with interFoam -parallel) suggest that defined the "alphaPhi0.water" for the "phasePhiCompressed". However, after some tests, this method will let some scalar stuck in the air. The "scalarTransport" utility of ESI OpenFOAM tries to conserve the mass for s.water instead of alphaS.water, as shown in the following code.

 
        tmp<surfaceScalarField> tTPhiUD;
        for (label i = 0; i <= nCorr_; i++)
        {
            fvScalarMatrix sEqn
            (
                fvm::ddt(s)
              + fvm::div(limitedPhiAlpha, s, divScheme)
              - fvm::laplacian(D, s, laplacianScheme)
              ==
                alpha*fvOptions_(s)
            );

            sEqn.relax(relaxCoeff);
            fvOptions_.constrain(sEqn);
            sEqn.solve(mesh_.solverDict(schemesField_));

            tTPhiUD = sEqn.flux();
        }

I hope this information helps,

Rdf

[37269]

It looks like they are working on a fix for this in the ESI-version:

https://develop.openfoam.com/Develop.../-/issues/2053

[randolph]

Hi,

In the case of fox may encounter this issue. Here is one of my juggling by integrating the "phaseScalarTransport" (org version) with ESI. https://github.com/Rdfing/interAdsFoam

Thanks,
Rdf

[hhu_lulu]

Hi, randolph.
Your comment is of great help in imposing phaseScalarTransport.However, I'm still confused about its usage.
For instance, what do you mean by 'The field alphaS.water = alpha.water * S.water should be used'? Should I specify a new variable named alphaS.water?
I post part of my functions in controlDict file below(interFoam), could you help checking its correctness? Thank you in advance!

functions
{
phaseScalarTransport1
{
type phaseScalarTransport;
libs ("libsolverFunctionObjects.so");
alpha alpha.water;
field alpha.water * s.water;
p p_rgh;
writeAlphaField true;
nCorr 1;
resetOnStartUp false;
writeControl outputTime; // write scalar field results
writeInterval 20;
// D 0.001; //difussion coefficient
// nut nut; //name of field to use as diffusivity, default = 'none'
// alphaDt 2;
// schemesField U; //name of field to use when looking up schemes from fvSchemes
bounded01 true;

[randolph]

Originally Posted by hhu_lulu

Hi, randolph.
Your comment is of great help in imposing phaseScalarTransport.However, I'm still confused about its usage.
For instance, what do you mean by 'The field alphaS.water = alpha.water * S.water should be used'? Should I specify a new variable named alphaS.water?
I post part of my functions in controlDict file below(interFoam), could you help checking its correctness? Thank you in advance!

functions
{
phaseScalarTransport1
{
type phaseScalarTransport;
libs ("libsolverFunctionObjects.so");
alpha alpha.water;
field alpha.water * s.water;
p p_rgh;
writeAlphaField true;
nCorr 1;
resetOnStartUp false;
writeControl outputTime; // write scalar field results
writeInterval 20;
// D 0.001; //difussion coefficient
// nut nut; //name of field to use as diffusivity, default = 'none'
// alphaDt 2;
// schemesField U; //name of field to use when looking up schemes from fvSchemes
bounded01 true;

Shengjie,

You can just use phaseScalarTransport as explained in the header file (https://cpp.openfoam.org/v7/classFoa...Transport.html). alphaS.water will be the field for postprocessing.

Thanks,
Rdf

[foamF]

Thanks Randolph for explaining phaseScalarTransport.

I'm new to this function object. It looks your conclusion is alphaS.water is the scalar retained in water phase, and should be used for processing.

I wonder why S.water is present? Is it involved in the scalar transport equation? What's it's physical meaning and under what situation it should be used in processing?

Looking forward to your reply!

[randolph]

Originally Posted by foamF

Thanks Randolph for explaining phaseScalarTransport.

I'm new to this function object. It looks your conclusion is alphaS.water is the scalar retained in water phase, and should be used for processing.

Hi,

This is not the exact reason why alphaS.water should be used. The alphaS.water should be used is because (i) "the continuity equation is calculated for alphaS.water" in the phaseScalarTransport as illustrated in my previous post and (ii) ad-hoc flux correction method implemented in phaseScalarTransport (OpenFOAM 7 or 8 [org version]).

In the ESI-OpenFOAM (e.g., v1912) formulation (i.e., phase option in the "scalarTransport" utility), the mass of the scalar will not be conserved even if alphaS.water is used for the processing.

The way that phaseScalarTransport formulates the phase-constrained scalar transport for interFoam is very similar to how the scalar transport handled in a typical dispersed multiphase flow simulation. The only added effort is on flux correction routine, which is also very similar to the (momentum) flux correction in a typical N-S solver (if my memory serves).

Thanks,
Rdf

[Swift]

Hi Randolph,

I have read your posts with interest because I am attempting to model the scalar transport between phases, whereas you are interested in the scalar transport limited to a single phase.

I thought I would share my method with you and am interested in your comments. I am sure my approach would work in your case as well.

In my approach, I search through the mesh and compare the values of alpha. If alpha is between 0.05 and 0.95, then I set the diffusivity to a high value and if outside these bounds I set it to 0. This basically meant I have a diffusivity field for which I can solve the scalar transport equation. I actually took the scalar equation directly from scalarTransportFoam.

To make this work for your case, all that would be required would be to set a high diffusivity for the phase that you want and set diffusivity everywhere else to zero.

I created a test case based on the dam break case. I was particularly looking for the problem you identified in your first post, i.e. that the scalar was trapped in the one phase. For my approach, I do not seem to have this problem.

I was then interested in your evaluation of the continuity of the scalar. How did you do this so I can verify for myself with my approach? Anyway, seeing as for interFoam, both phases share the velocity field and the scalar equation comes straight from scalarTransportFoam, I can't imagine there would be continuity issues. Maybe I'm being niave though... Let me know your thoughts.

I've put the code as it is so far on github as well as a test case for you to evaluate. I compiled in on the 2106 ESI version. Here is the link to it:
It is a work in progress, so it may not be as well documented as it could be. Let me know if you have any questions.


My intention is to extend the code and make the diffusivity at the interface dependent on the turbulence as well as alpha and finally to have a maximum limit for the scalar in the phases. In this case, even if there is a difference in value between two cells, but one cell has reached its limit, no more scalar will be diffused across the faces. This is analogous to oxygen transfer of air to water. Air can have 23% O2, whereas the maximum concentration of O2 in water is about 8.5%

I'd appreciate your (or anyone else's) thoughts and comments on my approach.

Thanks,
Thomas

[JVNR]

Hi everyone! I'm on the brink of a nervous breakdown and about to torch my pc.
Hopefully some of you can give me a hand.

I'm trying to determine the mixing time of a fermenter (stirred vessel with water and bubbles)

To do so I was hoping to use phaseScalarTransport on twoPhaseEulerFoam. Inject some known quatity of scalar S (topoSet) and probe some points in paraview to do the post processing.

The simulation runs normally but no s.water file is created in each of the timesteps. When looking at the log file I see that no calculations are done for the scalar.


Here what I consider the relevant parts of my code:

Regarding boundary conditions the scalar has zerogrdient everywhere since it's being generated by the injectionRateSuSp

In controlDict I put the following (based on https://cpp.openfoam.org/v7/classFoa...ransport.html)

functions
{
 s
  {
    type            phaseScalarTransport;
    libs            ("libsolverFunctionObjects.so");
    field           s.water;
    alphaPhi        alphaRhoPhi.water;
    rho             thermo:rho.water;
    resetOnStartUp  no;
    fvOptions
        {
            s
            {
                type            scalarSemiImplicitSource;

                timeStart       0;
                duration        1;
                active          true;
                selectionMode   cellSet;
                cellSet         tracer;
                volumeMode      absolute;
                injectionRateSuSp
                {
                    s.water     (10 0); // kg/s
                }
            }
    }
  };
};

Additionally I added my fvScheme and my fvSolution

fvSchemes
    div(phi,s.water)     Gauss linearUpwind grad(s.water); 

fvSolution
    "s.*"
    {
        solver          PBiCGStab;
        preconditioner  DILU;
        tolerance       1e-08;
        relTol          0;
        minIter	     1;
    }

Can anyone help me? I'm really struggling here.