[Daniel_Khazaei]

Dear foamers,

I'm working on a problem dealing with evaporation and condensation of water inside a closed domain. The fixedFluxPressure boundary condition is applied on all the boundaries which necessitate the need for setting a reference pressure before solving the pdEqn.

The evaporation-condensation model adds a mass transfer source in pdEqn as follow:

        fvScalarMatrix pdEqn
        (
            fvc::div(phiHbyA)
          - fvm::laplacian(rAUf, pd_)
          - (vDotc - vDotv)
        );

The problem is the reference pressure location has an unexpected effect on the simulation which seems totally unrealistic. If I remove the source term from the equation, then everything becomes completely normal. This finally cause the temperature field to become unbounded and drops way below physical temperature range allowed inside the domain. The contours of velocity and temperature shows the problem:





If it helps further investigating the issue, here is the full pEqn function:

void interPhaseChangeFluid::pEqn
(
    const fvMesh& mesh, const fvVectorMatrix& UEqn
)
{
    rAU_ = 1.0/UEqn.A();
    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU_));

    volVectorField HbyA("HbyA", U_);
    HbyA = rAU_*UEqn.H();

    surfaceScalarField phiHbyA
    (
        "phiHbyA",
        (fvc::interpolate(HbyA) & mesh.Sf())
      + fvc::ddtPhiCorr(rAU_, rho_, U_, phi_)
    );

    #include "updateRobinFsiInterface.H"

    adjustPhi(phiHbyA, U_, pd_);

    //- calculating surface tension flux
    //- + gravitational flux
    surfaceScalarField phig
    (
        interface_.capillaryFlux(rAUf)
      - ghf_*fvc::snGrad(rho_)*rAUf*mesh.magSf()
    );

    //- adding the resulting flux to the main flux
    phiHbyA += phig;

    Pair<tmp<volScalarField> > vDot = phaseChangeModel_->vDot();
    const volScalarField& vDotc = vDot[0]();
    const volScalarField& vDotv = vDot[1]();

    while (pimple_.correctNonOrthogonal())
    {
        fvScalarMatrix pdEqn
        (
            fvc::div(phiHbyA)
          - fvm::laplacian(rAUf, pd_)
          - (vDotc - vDotv)
        );

        pdEqn.setReference(pRefCell_, getRefCellValue(pd_, pRefCell_));

        pdEqn.solve
        (
            mesh.solutionDict().solver(pd_.select(pimple_.finalInnerIter()))
        );

        if (pimple_.finalNonOrthogonalIter())
        {
            phi_ = phiHbyA + pdEqn.flux();

            U_ = HbyA + rAU_*fvc::reconstruct((phig + pdEqn.flux())/rAUf);
            U_.correctBoundaryConditions();
            gradU_ = fvc::grad(U_);
        }
    }

    //- re-evaluate pressure using additional terms
    p_ = pd_ + rho_*gh_;
    p_.correctBoundaryConditions();
    gradp_ = fvc::grad(p_);

    if (pd_.needReference())
    {
        p_ += dimensionedScalar
        (
           "p",
            p_.dimensions(),
            pRefValue_ - getRefCellValue(p_, pRefCell_)
        );
        pd_ = p_ - rho_*gh_;
    }

    //- Make the fluxes relative to the mesh motion
    fvc::makeRelative(phi_, U_);
}

Regards,
D.Khazaei

[t.teschner]

So just to clarify, the source term that you added is

 (vDotc - vDotv)

right? It is actually funny, I am dealing with a similar issue (adding a source term to the equations) which then also gives me an unreasonable pressure distribution. So my thinking is that there is probably a better way (a more "openfoam" way) to incorporate the source term into the equations but I had not had time to return to the problem. If you are interested, here is my issue:Dual time stepping in OpenFoam

Anyways, so what I know so far is that solving the equation explicitly or implicitly does not seem to make a change (no influence of the matrix coefficients due to the implicit solver), I also thought that the source term may be multiplied by the volume but that doesn't change anything (also it seems strange, other threads on adding source terms to openfoam do not do this, but I tried anyway, I guess to calm my paranoia ...).

Are you sure that the issue is due to the reference pressure? There is also a method called constrainPressure(), maybe that would help in your case?

[Daniel_Khazaei]

Originally Posted by t.teschner

So just to clarify, the source term that you added is

 (vDotc - vDotv)

right?

Yes, it is the net mass transfer rate due to the evaporation/condensation.

Originally Posted by t.teschner

It is actually funny, I am dealing with a similar issue (adding a source term to the equations) which then also gives me an unreasonable pressure distribution.

Well we are both dealing with adding source terms, But the difference is that I'm not getting a unreasonable pressure field. As you can see both the velocity and temperature fields are acting strangely. The pressure field is OK!
I think it has something to do with velocity prediction after solving for pd where the velocity is corrected using pdEqn.flux().

Originally Posted by t.teschner

So my thinking is that there is probably a better way (a more "openfoam" way) to incorporate the source term into the equations but I had not had time to return to the problem.

Unfortunately I have tried both implicit and explicit approaches with slight differences in the results but the problem is there for both of them:

explicit implementation: if you look into interPhaseChangeFoam and interCondensingEvaporatingFoam solvers it can be seen that this is a common approach.

    fvScalarMatrix pdEqn
    (
          fvc::div(phiHbyA)
        - fvm::laplacian(rAUf, pd_)
        - (vDotc - vDotv)
    );

implicit implementation:

    fvScalarMatrix pdEqn
    (
          fvc::div(phiHbyA)
        - fvm::laplacian(rAUf, pd_)
        - fvm::Sp((vDotc - vDotv) / (pd_ + pd0), pd_)
    );

let OpenFOAM decide implicit or explicit incorporation based on the sign of (vDotc - vDotv):

    fvScalarMatrix pdEqn
    (
          fvc::div(phiHbyA)
        - fvm::laplacian(rAUf, pd_)
        - fvm::SuSp((vDotc - vDotv) / (pd_ + pd0), pd_)
    );

non of the above methods work correctly when

pd_.needReference()

returns true and subsequently

pdEqn.setReference(pRefCell, pdRefValue)

is used to enforce pressure level. Otherwise, everything works perfectly if one of the pressure boundaries fixes a pressure value.

Originally Posted by t.teschner

If you are interested, here is my issue:Dual time stepping in OpenFoam

I will take a look at your thread in details, thank you.

Originally Posted by t.teschner

Anyways, so what I know so far is that solving the equation explicitly or implicitly does not seem to make a change (no influence of the matrix coefficients due to the implicit solver), I also thought that the source term may be multiplied by the volume but that doesn't change anything (also it seems strange, other threads on adding source terms to openfoam do not do this, but I tried anyway, I guess to calm my paranoia ...).

Yes, as I had said earlier implicit or explicit doesn't change anything.
However, I don't have any problem with the source term if I don't try to model a closed-domain. So my problem seems to be related to setting a reference pressure when having a source term in pEqn.

Originally Posted by t.teschner

Are you sure that the issue is due to the reference pressure? There is also a method called constrainPressure(), maybe that would help in your case?

I'm not sure but the signs point me to that direction as the problematic location changes when pRefCell is moved. It follows the location of the reference pressure location, but the pressure field seems unaffected. Only the velocity and temperature fields are wrong which are related by phi.
I don't think foam-extend has such a library called constrainPressure(), I will look into it. However, both adjustPhi and constrainPressure modify boundaries of the domain. My problem is not limited to boundaries, if I put the reference location inside the domain, then I will see unrealistic velocity and Temperature there.


Regards

[t.teschner]

my bad, could you then indicate what the correct velocity and temperature field should look like? (I presume it is just a homogeneous field of velocity and temperature, i.e. no spikes as seen in your contour plots then?)

It seems you have already put more effort into understanding the issue than I could, so I have some new ideas to try for my issue (thanks for that), but now that I understand your problem a bit better, I agree that our problems are probably fundamentally different.

But if I understand that correctly, you do not have the issue when you use at least on one boundary a dirichlet condition for the pressure? I.e. the problem only arises when you have fully Neumann boundary conditions? I remember a similar problem when I was writing a simplified code for the lid driven cavity problem, there the pressure also has fully Neumann BCs and when prescribing a reference cell and value that approach was always giving wrong results for the cell where I put the reference pressure (similar to what you are showing). In the end I removed the reference pressure approach and just calculated the average pressure and removed that from the simulation (did the job, I guess you could say it is more of a global approach, whereas prescribing the reference cell could be seen as a local approach). Anyways, in openfoam that may not work as easily (and may also not be generalised for different cases), so I am not sure what to do about that.

One thing I am a bit more surprised about is that the pressure field seems to be OK while the velocity field is not. They should be coupled and even if the issue is only within the velocity field, at the next update it should then also affect the pressure. Thats a bit odd.

And the last thing, you said that you expect the problem to be with "pdEqn.flux()", but in-fact that method is always called (?!), while you say the issue only arises when if (pd_.needReference()) evaluates to true? Then I would expect the issue to be there, what makes you think it is the flux method then? (just trying to understand your thinking, maybe I missunderstood something along the way)

[Daniel_Khazaei]

Originally Posted by t.teschner

my bad, could you then indicate what the correct velocity and temperature field should look like? (I presume it is just a homogeneous field of velocity and temperature, i.e. no spikes as seen in your contour plots then?)

Thank you for looking into my thread.

Yes, exactly. I have also installed ANSYS Fluent student version and performed the simulation there. There should never be such a spike in velocity and temperature fields and those fields should start changing smoothly as the simulation goes on.

Originally Posted by t.teschner

It seems you have already put more effort into understanding the issue than I could, so I have some new ideas to try for my issue (thanks for that), but now that I understand your problem a bit better, I agree that our problems are probably fundamentally different.

Well yes, the problem I have only occurs when pd_ needs a reference level to be set.

Originally Posted by t.teschner

But if I understand that correctly, you do not have the issue when you use at least on one boundary a dirichlet condition for the pressure? I.e. the problem only arises when you have fully Neumann boundary conditions?

exactly.

Originally Posted by t.teschner

I remember a similar problem when I was writing a simplified code for the lid driven cavity problem, there the pressure also has fully Neumann BCs and when prescribing a reference cell and value that approach was always giving wrong results for the cell where I put the reference pressure (similar to what you are showing). In the end I removed the reference pressure approach and just calculated the average pressure and removed that from the simulation (did the job, I guess you could say it is more of a global approach, whereas prescribing the reference cell could be seen as a local approach). Anyways, in openfoam that may not work as easily (and may also not be generalised for different cases), so I am not sure what to do about that.

Can you describe the approach in more details?
When I was analyzing the same case with ANSYS Fluent, it says that to see the actual reference location you need to run the following TUI command which is only available when the simulation doesn't have any pressure related boundary conditions:

define -> operating-conditions -> used-ref-pressure-location

When I execute this command, it returns nothing which I think it means that Fluent is not using a reference pressure for my simulation.
I have tried to remove pdEqn.setReference and that just cause immediate divergence.

Originally Posted by t.teschner

One thing I am a bit more surprised about is that the pressure field seems to be OK while the velocity field is not. They should be coupled and even if the issue is only within the velocity field, at the next update it should then also affect the pressure. Thats a bit odd.

Yes, that is also the same issue that keep me thinking how on earth the solution of pd is causing such a problem on momentum equation and pd itself is not affected

Originally Posted by t.teschner

And the last thing, you said that you expect the problem to be with "pdEqn.flux()", but in-fact that method is always called (?!), while you say the issue only arises when if (pd_.needReference()) evaluates to true? Then I would expect the issue to be there, what makes you think it is the flux method then? (just trying to understand your thinking, maybe I missunderstood something along the way)

I think I didn't put my word as clear as possible here, what I mean is that when (pd_.needReference()) evaluates to true, then something goes wrong in the solution of pd which subsequently results in wrong calculation of the flux, e.g. pdEqn.flux(). As the whole domain is at rest before the start of the simulation, then a slight spike in phi causes nonphysical velocity field.

Maybe the solution doesn't converge properly when I insert the source and the pdEqn.setReference() is active simultaneously. Although that doesn't answer why pd_ field doesn't show any anomaly if there is a convergence problem.

        if (pimple_.finalNonOrthogonalIter())
        {
            phi_ = phiHbyA + pdEqn.flux();

            U_ = HbyA + rAU_*fvc::reconstruct((phig + pdEqn.flux())/rAUf);
            U_.correctBoundaryConditions();
            gradU_ = fvc::grad(U_);
        }

[t.teschner]

I see, thanks for the additional explanation. This might not work (or be pointless), but what happens if you call pd_.correctBoundaryConditions(); ? And could you also explain what the pd equation actually solves?

Regarding the pressure boundary implementation, I lterally just looped over all cells, added the pressure to a counter and then diveded the counter by the number of elements. That was then my average pressure and I removed that from my pressure field. Or in pseudo code

averagePressure = 0
for element in totalElements {
  averagePressure += pressure[element]
}

averagePressure /= totalElements
for element in totalElements {
  pressure -= averagePressure
}

as I said, that may not work for general cases but for those fully Neumann problems I looked at it was an alternative to setting the reference pressure.


I am actually not sure how fluent handles the reference pressure. You can specify your operational pressure which, for incompressible flows, should be anyways irrelevant. During post processing that value is used, however, to calculate your total pressure but I am not sure if the value is used in the simulation itself.

[Daniel_Khazaei]

Originally Posted by t.teschner

I see, thanks for the additional explanation. This might not work (or be pointless), but what happens if you call pd_.correctBoundaryConditions(); ? And could you also explain what the pd equation actually solves?

Unfortunately, that doesn't change anything.
Well it's just the same as p_rgh in OpenFOAM, with different naming in foam-extend.
It solves for the pressure subtracted by rho*gh.

Originally Posted by t.teschner

Regarding the pressure boundary implementation, I lterally just looped over all cells, added the pressure to a counter and then diveded the counter by the number of elements. That was then my average pressure and I removed that from my pressure field. Or in pseudo code

averagePressure = 0
for element in totalElements {
  averagePressure += pressure[element]
}

averagePressure /= totalElements
for element in totalElements {
  pressure -= averagePressure
}

as I said, that may not work for general cases but for those fully Neumann problems I looked at it was an alternative to setting the reference pressure.

That sounds like an interesting replacement for local referencing approach, I will look into it.

Originally Posted by t.teschner

I am actually not sure how fluent handles the reference pressure. You can specify your operational pressure which, for incompressible flows, should be anyways irrelevant. During post processing that value is used, however, to calculate your total pressure but I am not sure if the value is used in the simulation itself.

Well then the question will be how fluent handle closed-domain simulations when non of the boundaries fixes a value on the pressure and the working fluids are in-compressible.

Regards

[Daniel_Khazaei]

today I have noticed something that I can't understand. I have modified the code to output pd_ values at each step:

        fvScalarMatrix pdEqn
        (
            fvc::div(phiHbyA)
          - fvm::laplacian(rAUf, pd_)
          - (vDotc - vDotv)
        );

        Info << "1- pd before setting a reference = " << pd_[pRefCell_] << endl;

        pdEqn.setReference(pRefCell_, getRefCellValue(pd_, pRefCell_));

        Info << "2- pd after setting a reference = " << pd_[pRefCell_] << endl;

        pdEqn.solve
        (
            mesh.solutionDict().solver(pd_.select(pimple_.finalInnerIter()))
        );

        Info << "3- pd after solve = " << pd_[pRefCell_] << endl;

When the source term is active in pdEqn, pd_ values at step 1 and 2 are the same. However the pd_ value at step 3 (after solving for pd_) differs significantly. Is this the expected behavior when we are dealing with source terms in pdEqn?

I though that pdEqn.setReference is supposed to actually keep pd_ value at pRefCell constant. Am I wrong here? I think that may reveal why prediction of velocity using pdEqn.flux() causes non-physical values.

If I remove the source term from pdEqn, then pd_ value remains constant at all the steps and I don't get non-physical velocity or temperature field anymore.

[t.teschner]

ahhh, I see, in-fact, the following line doesn't make any sense (at least to me)

pdEqn.setReference(pRefCell_, getRefCellValue(pd_, pRefCell_));

it looks like you are first asking the solution field for what the reference pressure at the reference cell location is with the second argument (i.e. getRefCellValue(pd_, pRefCell_)) and then you set that value, that you just got, at the same location where you have read if from (first argument, pRefCell_). So in effect, this call does not do anything (well, in the sense that you are not changing the reference pressure), can you prescribe some reference pressure here? you could try the following and see if that works:

pdEqn.setReference(pRefCell_, pRefValue);

where you have to define pRefValue (scalar pRefValue = 0.0; ). But it is probably already defined in you createFields.H file?

[Daniel_Khazaei]

Originally Posted by t.teschner

ahhh, I see, in-fact, the following line doesn't make any sense (at least to me)

pdEqn.setReference(pRefCell_, getRefCellValue(pd_, pRefCell_));

it looks like you are first asking the solution field for what the reference pressure at the reference cell location is with the second argument (i.e. getRefCellValue(pd_, pRefCell_)) and then you set that value, that you just got, at the same location where you have read if from (first argument, pRefCell_). So in effect, this call does not do anything (well, in the sense that you are not changing the reference pressure), can you prescribe some reference pressure here? you could try the following and see if that works:

pdEqn.setReference(pRefCell_, pRefValue);

where you have to define pRefValue (scalar pRefValue = 0.0; ). But it is probably already defined in you createFields.H file?

well this line actually in combination with the following lines makes sense to me

    //- re-evaluate pressure using additional terms
    p_ = pd_ + rho_*gh_;
    p_.correctBoundaryConditions();

    if (pd_.needReference())
    {
        p_ += dimensionedScalar
        (
           "p",
            p_.dimensions(),
            pRefValue_ - getRefCellValue(p_, pRefCell_)
        );
        pd_ = p_ - rho_*gh_;
    }

What we need is keeping the gauge pressure equal to zero at reference pressure location. The gauge pressure is calculated from the following line:

    p_ = pd_ + rho_*gh_;
    p_.correctBoundaryConditions();

Then in order to keep the gauge pressure zero at ref. location, we subtract the value of p_ at that location from the p_ field by executing the following line and finally re-evaluating pd_ field based on the new p_ field:

    if (pd_.needReference())
    {
        p_ += dimensionedScalar
        (
           "p",
            p_.dimensions(),
            pRefValue_ - getRefCellValue(p_, pRefCell_)
        );
        pd_ = p_ - rho_*gh_;
    }

This implies that pd_ value at the reference pressure location should be equal to - rho_*gh_. Then this value is used to set a reference for pdEqn equation using the following line:

pdEqn.setReference(pRefCell_, getRefCellValue(pd_, pRefCell_));

Of course I have already tried using your suggestion which just shift the pressure values but the actual difference is still there which I think is reasonable as the working fluid is in-compressible.

pdEqn.setReference(pRefCell_, pRefValue);

I think the problem is why pdEqn.setReference can't keep the pd_ value constant after solving the pdEqn when the pdEqn contains any mass transfer source term.

The interesting think is that I tried using TEqn.setReference to enforce constant temperature at ref. cell but the exact same problem is there. If TEqn contains any source term then setReference can't keep T constant at the specified value.

[Daniel_Khazaei]

just looking for any reason behind this issue, I have just noticed that there is a bug report about it:

https://develop.openfoam.com/Develop...lus/issues/465

As I'm using the exact same mass transfer model in my solver, I think the same limitation is inherited:

"The mass can not be conserved when any mass transfer model is involved."

So, no closed-doamin for now until I have enough time to come up with a solution.
Probably ANSYS found a solution for this, because it handles these kind of simulation without any problem.

[t.teschner]

I see, well, half satisfying but at least we know its not your code then

[mxylondon]

Unfortunately, I got the same issue. Having no idea why and how to fix.

As you can see from my screenshot, the pressure spike is at the reference cell (pRefCell=0). This leads to a crash later.

My solver for this case is buoyantSimpleFoam. Mesh quality is quite good, since the geometry is simple.

Help needed and really appreciate in advance!!!