[Jost K]

Dear Foamers,

I am having some trouble setting up a large scale oceanic flow case in buoyantPimpleFoam.
Since I cannot set a (uniform) fixed pressure for any of my boundaries I need to specify the reference pressure through pRefCell and pRefValue.
However, if I do so, the simulation immediately crashes and it looks like I have a "black hole" inside my pRefCell.



This problem can easily be reproduced using the hotRoom tutorial for buoyantPimpleFoam:
just set the gravitational acceleration in constant/g to (0 -98.1 0) (This has the same effect as increasing the scale of the case and changing the medium to water.)
The simulation should crash after 10s

If it does, change the writeInterval in system/controlDict to 1 and run again (to get some output you can view in paraview).


The result should look like this:
Selection_010.png



Does anyone have an idea how to avoid this behavior?

Are we seeing a general bug here?

I would really appreciate any suggestions!


Cheers,
Jost






P.S.: Please note that I am not using g=-98.1 [m/s^2] in my real case but this is just a simple way of reproducing my problem without having to increase the size of the hotRoom tutorial and changing the medium to water, which would represent my actual case.

[Jost K]

Hi All,


I realized for the example case (hotRoom) the problem can be solved by choosing a smaller time step size.
(I was already about to be very embarrassed.)

But I found that this does not help for my actual case, which is a flow of water in a domain with a depth of 150m.
It even seems like the solver crashes after less time steps if i choose a smaller time step size.
Also, I realized that the solution looks good after the first time step and then the "black hole" develops from there.

I have almost tried everything now:

• Smaller time steps
• Different fvShemes
• Relaxation of p_rgh
• Tolerance of p_rgh

Any more Ideas?


Cheers,
Jost

[piotr.mecht]

Originally Posted by Jost K

It even seems like the solver crashes after less time steps if i choose a smaller time step size.

I had similiar problem with running buoyant water simulation (same quoted symptom). I based on the hotRoom tutorial and i didn't know the equations simulated back then. Simulations for much less compressible water don't need solving the pressure-work term in the energy equation.
https://caefn.com/openfoam/solvers-buoyantpimplefoam


you can disable this by adding non-default entry to your thermophysicalProperties dictionary:

dpdt off;

My simulation went immediately stable. Maybe this is also your problem?

[Jost K]

Hi Piotr,


Thanks a lot for your your reply!
Switching dpdt off definitely sounds like a good idea for my case.
However, my problem still persists.
It looks like GAMG needs a lot of iterations in my second (and final) pressure loop.
Maybe it's because of my thermo model, I am using:

    
thermoType
    {
        type                    heRhoThermo;
        mixture               pureMixture;
        specie                specie;
         equationOfState icoPolynomial;
         thermo               hPolynomial;
         transport            polynomial;      
         energy               sensibleInternalEnergy;
    }

If I use the thermophysical properties from the hotRoom tutorial it works fine (but that is an air flow).

[piotr.mecht]

Hi
You are using the sensibleInternalEnergy formulation for the energy equation, which is equivalent for sensibleEnthalpy (it is shown in the article in the link from my previous post). If you'll look at the equations (4) and (5) you'll see, that the pressure change in the form of dp/dt is expressed only in the sensibleEnthalpy equation. You'll not disable pressure-work term by disabling dpdt in the sensibleInternalEnergy formulation.


My thermophysicalProperties looked something like:

thermoType
{
    type            heRhoThermo;
    mixture         pureMixture;
    transport       polynomial;
    thermo          hPolynomial;
    equationOfState icoPolynomial;
    specie          specie;
    energy          sensibleEnthalpy;
}

pRef            400000;
dpdt off;

mixture
{
    specie
    {
        molWeight       18;
    }
    equationOfState
    {
        rhoCoeffs<8> (117.045534 7.913913 -0.02254739 1.97987E-05 0 0 0 0);
    }
    thermodynamics
    {
        Hf              0;
        Sf                0;
        CpCoeffs<8> (4200 0 0 0 0 0 0 0);
        
    }
    transport
    {
        muCoeffs<8> (0.1368992 -0.0012040579 3.5597468E-06 -3.526231E-09 0 0 0 0);
        kappaCoeffs<8> (-2.4283206 0.0233779 -5.9953665E-05 5.2602452E-08 0 0 0 0);
    }
}

Make sure, that your fvSolution has the solver for energy equation specified.
(h for sensibleEnthalpy and e for sensibleInternalEnergy).
For example:

"(U|h|e|k|omega|R)"
    {
        solver          PBiCG; // PBiCGStab;
        preconditioner  DILU;
        tolerance       1e-6;
        relTol          0.1;
    }

    "(U|h|e|k|omega|R)Final"
    {
        $U;
        relTol          0;
    }

[Jost K]

Originally Posted by piotr.mecht

Hi
You are using the sensibleInternalEnergy formulation for the energy equation, which is equivalent for sensibleEnthalpy (it is shown in the article in the link from my previous post). If you'll look at the equations (4) and (5) you'll see, that the pressure change in the form of dp/dt is expressed only in the sensibleEnthalpy equation. You'll not disable pressure-work term by disabling dpdt in the sensibleInternalEnergy formulation.

Infinite thanks to you Poitr!!!
You just made my month!

(yes this has been bugging me for quite a while now)



Using sensibleEnthalpy with dpdt off works perfectly.

To be honest, I was not really expecting that this problem would come from the energy equation. Well, yet another lesson on how important it is not only to know the equations but to actually understand them.


Thanks again,
Jost

[Jost K]

Hi All,
Sorry, I have more bad news:
Even though after removing the pressure work term from the energy equation ( dpdt off ) my case does not crash any more, my little "black hole" inside the pRefCell still persists.
I set up a simple example case to show the Problem.
It's a closed volume without any flow.
The dimensions are 20x20x150m (Length x Width x depth)
All the walls are slip walls.
The pressure reference point (pRefPoint) is in the middle of the top patch.
The Mesh looks like this:
Mesh.png
If you run the case for 1sec the pressure generally looks okay
p_after_1sec.png
However, there is a strange velocity at the top patch
U_after_1sec.png
Taking a closer look at the pressure on the top patch reveals the problem
p_onTopPatch_after_1sec.png
The test case is attched to this post. It runs with buoyantPimpleFoam on OpenFOAM-v1812.
I have tried the following things to fix the issue (without success):

• Different grad, snGrad and interpolation schemes
• Different tolerances for p_rgh (in fvSoulution)
• Different solvers for p_rgh
• Different positions of pRefCell and pRefValue
• Some tweaks on the pEqn
• Using buoyantSimpleFoam instead of buoyantPimpleFoam

This looks like a rather general problem.
Has anybody ever seen this before?
Your help would be highly appreciated!
Cheers,
Jost

[Jost K]

Hi All,

I just found out that switching off the momentum predictor step in fvSolution helps a lot.
For my case this might actually be an okay solution.
I am just trying to figure out what effect this has on the solution.

[piotr.mecht]

Originally Posted by Jost K

However, there is a strange velocity at the top patch

Well: your boundary conditions for p_rgh are that there is no gradient of pressure across your boundaries, nothing else, nothing more. This doesn't show any direction for equilibrium state, thats why any pressure reference point is needed. Your initial state for pressure is the uniform field, but when it calculates static pressure distribution it assumes that some unsteady flow caused this initial pressure difference, so there is your flow. It exist for a short time, and disappears after few iterations (it disappears a bit faster with smaller timestep).


Look at your p_rgh values. It is ~0.3 Pa difference, and your flow patch values are in magnitude of 1e-6 m/s, which are rather extremely small values if you'll imagine 150m water column.


You can get rid of this artifact by initializing your pressure field according to static distribution (non-uniform value instead of 7.5E+5 for the whole domain). I did this whith funkySetFields (tool from swak4Foam). My funkySetFieldsDict:

FoamFile
{
    version         2.0;
    format          ascii;
    class           dictionary;
    object          funkySetFieldsDict;
}
// usage: 
// funkySetFields -time 0
expressions
(
    pressureWater
    {
     field p; //field to initialize
     // 1022.601443 - water density according to your polynomial value
     // 1504758.0233745 - static pressure at the bottom according to your polynomial value
     expression "9.81*1022.601443*(0-pos().z)+1504758.0233745";
     condition  "(pos().z<=0) && (pos().z>=-150)";
         keepPatches 1; 
    }
    
);

With pressure initialization there's still a spot after first 0.1s, but the pressure difference in the whole domain is in 1e-30 order of magnitude (Attached picture). it disappears soon after. Looks strange on the default scale range, it's because there is no other flow in the domain.
I didn't have to disable momentum predictor. I'd leave it on;

[Jost K]

Hi Piotr,

thanks again for your reply.

You are definitely right in saying that both the pressure and velocity artifacts we are seeing here are very small (actually below our usual range of accuracy).

But in my actual case they are of the order of 30 Pa and the Velocity is 0.1 m/s. This is still not huge but has a significant influence on my solution. And the problem does not disappear after 5000s (even without momentum predictor).

My actual case is cold water from a depth of 75m being pumped up to the ocean surface region.
It looks like this:
myCase.png
Here is a slice that cuts trough my pRefPoint showing velocity after 5000s simulated time.
U_5000s.png
The effect becomes even clearer when looking at the temperature.
T_5000s.png
I initialize my hydrostatic pressure using the hydrostatic background pressure calculation from the fireFoam solver, which I adopted to buoyantPimpleFoam.
This gives me an initial pressure which includes the effects of my nonuniform initial density field (because of the temperature stratification).

My case setup is comparable to the one in the example case.

I still do not understand why I am seeing this artifact.
I understand that establishing the hyrostatic pressure may produce some initial flow. But why does that flow come from my pRefCell?

[Jost K]

Hi,

I just found out that changing the algorithm arrangement of buoyantPimpleFoam helps to completely get rid of the the described problem.

Solving the energy equation at the end of the pimple loop seems to work fine.
The solver then looks like this:

// --- Pressure-velocity PIMPLE corrector loop
        while (pimple.loop())
        {
            #include "UEqn.H"

            // --- Pressure corrector loop
            while (pimple.correct())
            {
                #include "pEqn.H"
            }

            if (pimple.turbCorr())
            {
                turbulence->correct();
            }


            #include "EEqn.H"
        }

The issue that, solving the energy equation based on a velocity (from the velocity predictor) that does not obey continuity, might lead to problems has been raised before in this forum (see this thread) and is also mentioned in Oliveira and Issa's paper on "AN IMPROVED PISO ALGORITHM FOR THE COMPUTATION OF BUOYANCY-DRIVEN FLOWS" (see here).

However I am not sure if that is really the problem here, since if it was, something like this:
(solving the energy predictor before the velocity predictor)

// --- Pressure-velocity PIMPLE corrector loop
        while (pimple.loop())
        {
            #include "EEqn.H"

            #include "UEqn.H"

            // --- Pressure corrector loop
            while (pimple.correct())
            {
                #include "pEqn.H"
            }

            if (pimple.turbCorr())
            {
                turbulence->correct();
            }
        }

or this:
(which would be somewhat along the lines of Oliveira and Issa's algorithm)

// --- Pressure-velocity PIMPLE corrector loop
        while (pimple.loop())
        {
            #include "UEqn.H"

            // --- Pressure corrector loop
            while (pimple.correct())
            {
                 #include "pEqn.H"
                 #include "EEqn.H"
             }

            if (pimple.turbCorr())
            {
                turbulence->correct();
            }
        }

should also work fine, which it doesn't (with regard to my problem).
Maybe someone else can make sense of all these findings.

Cheers,
Jost

[Jost K]

Hi All,
Still working on this problem.

I really seems to come down to a continuity problem.
Apparently, buoyantPimleFoam becomes very sensitive to continuity errors when cases have large spatial dimensions and high density.
Thus, everything that minimizes continuity errors also helps to get rid of the pressure problems described in this thread.
I have tried several things, which all had a positive effect:

1. Modify the Order of the PIMPLE algorithm as described in my previous post.
2. Get rid of any cyclicAMI boundary conditions (they introduce continuity errors as I found out here. Setting a high lowWeightCorrection e.g. 0.999 for the AMI boundariey also helps slightly.)
3. Set one velocity boundary condition to fluxCorrectedVelocity rather than keeping the volume closed completely.
   (Note that I set all the pressure boundary conditions to fixedFluxPressure, since setting them to fixedValue or totalPressure is generally a bad idea when dealing with nonuniform density)
4. Finally, I made some modifications to buoyantPimpleFoams pEqn. The most important one seems to be adding

   adjustPhi

   to adjust the boundary fluxes so that continuity is maintained. (This is done in buoyantSimpleFoam. I don't know why it is missing in buoyantPimpleFoam.)

My pEqn now looks like this:

    dimensionedScalar compressibility = fvc::domainIntegrate(psi);
    bool compressible = (compressibility.value() > SMALL);

    rho = thermo.rho();

    // Thermodynamic density needs to be updated by psi*d(p) after the
    // pressure solution
    const volScalarField psip0(psi*p);

    volScalarField rAU(1.0/UEqn.A());
    surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));

    surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());

    surfaceScalarField phiHbyA
    (
        "phiHbyA",
        fvc::flux(rho*HbyA)
        + MRF.zeroFilter(rhorAUf*fvc::ddtCorr(rho, U, phi))
    );

    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);

    phiHbyA += phig;
    bool closedVolume = adjustPhi(phiHbyA, U, p_rgh);

    // Update the pressure BCs to ensure flux consistency
    constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);

    fvScalarMatrix p_rghDDtEqn
    (
        fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
      + fvc::div(phiHbyA)
      ==
        fvOptions(psi, p_rgh, rho.name())
    );

    while (pimple.correctNonOrthogonal())
    {
        fvScalarMatrix p_rghEqn
        (
            p_rghDDtEqn
          - fvm::laplacian(rhorAUf, p_rgh)
        );

        p_rghEqn.setReference
        (
            pRefCell,
    //         compressible ? getRefCellValue(p_rgh, pRefCell) : pRefValue
            getRefCellValue(p_rgh, pRefCell)
        );

        p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));

        if (pimple.finalNonOrthogonalIter())
        {
            // Calculate the conservative fluxes
            phi = phiHbyA + p_rghEqn.flux();

            // Explicitly relax pressure for momentum corrector
            p_rgh.relax();

            // Correct the momentum source with the pressure gradient flux
            // calculated from the relaxed pressure
            U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
            U.correctBoundaryConditions();
            fvOptions.correct(U);
            K = 0.5*magSqr(U);
        }
    }

    p = p_rgh + rho*gh;

    #include "rhoEqn.H"
    #include "compressibleContinuityErrs.H"

    if (closedVolume)
    {
        if (!compressible)
        {
            p += dimensionedScalar
            (
                "p",
                p.dimensions(),
                pRefValue - getRefCellValue(p, pRefCell)
            );
        }
        else
        {
            p += (initialMass - fvc::domainIntegrate(psi*p))
                /compressibility;
            thermo.correctRho(psi*p - psip0);
            rho = thermo.rho();
    //         p_rgh = p - rho*gh;  //JK: moved down
        }
        p_rgh = p - rho*gh;     //JK: make ph_rgh consistent with changed p in any case!
    }
    else
    {
        thermo.correctRho(psi*p - psip0);
    }

    rho = thermo.rho();

    if (thermo.dpdt())
    {
        dpdt = fvc::ddt(p);
    }

Maybe someone with a better Idea of the theory can explain why buoyantPimpleFoam gets so sensitive for continuity errors.
I am pretty sure there is a better way of dealing with this issue that I haven't found yet.

Cheers,
Jost

[mxylondon]

Hi dear foamers,

Unlucky that I just met this issue in my project. The solver is buoyantSimpleFoam, which is a steady-state version.

In my analysis, it was observed that the predicted pressure and velocity are very strange at the location of reference cell (pRefCell=0 and pRefValue=101325). Then the analysis crashed after some iterations.

The mesh quality is not poor, since the geometry is simple in this test case.

If I changed reference cell to another location, then the strange result also moves to another location, which is still at the reference cell.

Can anyone kindly help? Really appreciate. Great thanks!

XM

[jherb]

If you specify the pressure at the outlet, there is no need to give pRefCell. What happens, if you remove it?

Originally Posted by mxylondon

Hi dear foamers,

Unlucky that I just met this issue in my project. The solver is buoyantSimpleFoam, which is a steady-state version.

In my analysis, it was observed that the predicted pressure and velocity are very strange at the location of reference cell (pRefCell=0 and pRefValue=101325). Then the analysis crashed after some iterations.

The mesh quality is not poor, since the geometry is simple in this test case.

If I changed reference cell to another location, then the strange result also moves to another location, which is still at the reference cell.

Can anyone kindly help? Really appreciate. Great thanks!

XM

[Jost K]

Hi XM,

did you solve your problem?
In my case the modifications I stated earlier still work fine.
Also a good initialisation of the hydrostatic pressure seems to be important.
If you can specify the pressure at any of the boundaries that would probably be the easiest way to solve your problem (as jherb stated). But be careful, with nonuniform density it is quite tricky to specify the pressure anywhere since you can not just assume the hydrostatic pressure to be rho*g*h (see buoyantSimpleFoam and watertank).

Cheeres,
Jost

[mxylondon]

Hi Joachim,

Many thanks for your attention. My case is a cavity-like problem, although it looks like a pipe. That's the reason why I used a reference pressure instead of a pressure at a patch, such as outlet in a normal flow case.

Cheers

[mxylondon]

Hi Jost,

Thank you very much for your response. I think your solution of additional work in energy equation does make a lot of sense. If my memory was correct I also did similar things around 2013 to fix a very similar problem for my PhD project.

However unfortunately I am working for a company role now and it is not convenient any more to compile a solver...

Currently I am still working on this issue and trying to find a solution without adding corrections to the solver itself.

Definitely I will let you know if any progress could be made before I am laid off because of this, haha

XM

[Jost K]

Hi XM,

I like the approach of not getting your hands dirty on the solver itself!

You have probably tried this already, but I recently found that sometimes silly amounts of nonOrthoCorrector steps are needed to get p_rgh converged.
Also, my guess is that this problem (at least partially) arises from the fact that the temperature and thus the density is computed on the basis of a velocity field that does not obey continuity (i.e. the predicted velocity field). Thus, if you do not want to rearrange the SIMPLE algorithm itself, my guess would be that you have to crank up the relaxation for U (and maybe h/e).
But you are probably already way past this by now...

Cheers,
Jost

[mxylondon]

Morning Jost and dear all,

I think I have fixed my problem yesterday. Here is to share my experience:

1. I moved the reference cell from boundary wall to an internal point

2. I realised that I don't really need to create cyclic interfaces between my porous zones and outer normal flow domain. Passing data via these cyclic boundaries will definitely induce potential issues (compared to common interpolations across numerical cells). So I deleted all those cyclic boundaries between porous zones and outer flow domain

3. I reset the boundary conditions for k and epsilon according to the general velocity of natural flow in the domain

Running the analysis again, it is working well so far and the sum local error goes back to very small values again, at e-6 level.

Thank you so much for your attention. I am getting so much power from here, since 2012. Now I am ready for the next annoying problem. Into the unknown, ah...ah.....