[Matt_B]

Hi all,
I'm running a test on buoyantPimpleFoam solver using the default "hotRoom" example, with OpenFoam 2.1.1.
This is a solver which involves a compressible fluid, solving a momentum, pressure and enthalpy equations; the domain is a 3D box and there is also a dependence on time; the convergence method is "pimple".
This is a brand new solver for me, so I'm trying to understand how does it work and which are the potentials to be exploited.
I figured out that, before starting to customize this solver, I should verify if it works well either with light conditions and with heavy conditions.
Thus, I started with the lightest condition possible: zero energy transfer.
For time "0", I defined a homogeneous temperature inside the box (300K) with all the walls at the same, fixed temperature.
I also defined a zero velocity at all the walls and of course inside the domain; for the pressure, the boundary conditions remained as default definition (p_rgh as "buoyantPressure" and p as "calculated" at all walls).
In brief, I would have expected that all the parameters involved would have kept unchanged values because the box is an isolated system with no interaction with the external and no energy fluxes into it.
But, with great surprise, I found out that the pressure increases: p_rgh passes from 100000 Pa to 100028 Pa just after 0.1 sec, with a quite constant temperature. There is also a very small velocity flux from the bottom to the top, depending on p_rgh gradient.
Could somebody explain to me what is happening with this? Is this a known issue? Why does pressure increase even if the boundary conditions for all parameters involved are set with the aim not to transfer energy inside the domain?
I believe the problem shall be in the pressure calculation; I first tried to under-relax p_rgh but the issue persisted.
Decreasing the time step leads to same results and the pressure is always higher than the original 100000 Pa.
I read in this Forum about pressure in a closed domain, that when there is no explicit definition of it on a single wall at least, it could bring to wrong results....but I also read that this shall be an issue for steady-state solvers and not for transient solvers as it is.

Every help would be appreciated.

Thank you.

Matteo

[Matt_B]

None experienced such an issue with this solver? This is driving me crazy, because I didn't change a thing in the solver so I cannot even imagine where the problem is.

[karlli]

Dear Matt,
Please have a look at the log from your simulation. You will find a line like this:

time step continuity errors : sum local = 3.86557e-06, global = -1.66908e-06, cumulative = 0.00127436

Running the hotRoom tutorial with the standard settings (only one pressure corrector loop) will cause some continuity errors since the continuity (i.e. p_rgh) equation is only solved approximately in each time step. In my case the accumulated error is of the order 10^-3. Since mass is proportional to pressure, we can expect a similar error in the pressure. Your error of 28 Pa therefore looks quite reasonable, just slightly higher than 10^-4.

Hope this helps.

Best regards,
Karl

[Matt_B]

Originally Posted by karlli

Dear Matt,
Please have a look at the log from your simulation. You will find a line like this:

time step continuity errors : sum local = 3.86557e-06, global = -1.66908e-06, cumulative = 0.00127436

Running the hotRoom tutorial with the standard settings (only one pressure corrector loop) will cause some continuity errors since the continuity (i.e. p_rgh) equation is only solved approximately in each time step. In my case the accumulated error is of the order 10^-3. Since mass is proportional to pressure, we can expect a similar error in the pressure. Your error of 28 Pa therefore looks quite reasonable, just slightly higher than 10^-4.

Hope this helps.

Best regards,
Karl

Thank you Karl for your answer.
I know about continuity errors and here following you can see what I found after one step at 0.1 sec:

time step continuity errors : sum local = 8.528246093e-12, global = 2.864539688e-12, cumulative = -3.093536523e-06

Hence, the cumulative error is very small. By the way, I have to say that not all the equations involved converge below a satisfied tolerance: with the 4000 cells original mesh and the initial conditions I described in the first post, after 0.1 sec I only see the velocity Uz and p_rgh converging below a tolerance of 1e-08; Ux and Uy remain above 1e-03; I just uncommented out the energy equation from the solver for this test case, just to further simplify it.
I have also to highlight that the z axis is parallel to the gravity vector g; the original mesh got the gravity vector parallel to y axis.
I still don't understand why Ux and Uy don't converge, because there is no momentum transportation inside the domain and the velocity is everywhere set to (0, 0, 0)... actually, I see that velocity is generated by a kind of p_rgh gradient that appears during the calculation but, in my opinion this is not physical but very odd!

Matteo

[pupo]

I've been dealing with this problem (and trying to fix it) for quite some time.

The reason for this is quite simple. buoyantPimpleFoam does the following step at the cretefield.H

p_rgh = p - rho*gh;

this re-writes p_rgh to somthing that is NOT constant throughout the domain, and, therefore, the solver starts with a "bad" solution. The initial parameters with which the solver starts are

rho max/avg/min : 1.15862 1.15862 1.15862
T   max/avg/min : 300 300 300
P   max/avg/min : 100000 100000 100000
Prg max/avg/min : 100057 100028 100000

The correct solution would be p_rgh = constant.

I've been trying to fix this by replacing the above mentioned code of the createfield.h by

thermo.rho() = p_rgh / (  RperfectGas * thermo.T() - gh );
thermo.p()   = p_rgh + thermo.rho()*gh;

this, generate a starting solution with the following values:

rho max/avg/min : 1.15928 1.15895 1.15862
T   max/avg/min : 300 300 300
P   max/avg/min : 100000 99971.6 99943.1
Prg max/avg/min : 100000 100000 100000

Which are, according to the ideal gas law, the stable "correct" solution of the system.

However, this initial solution behaves rather strangely: it sometimes work, it sometimes doesn't. In your particular case it does NOT work.

Using my edited solver i get at Time = 0.0892992

rho max/avg/min : 1.15895 1.15862 1.15829
T   max/avg/min : 300.004 300 299.996
P   max/avg/min : 100028 100000 99971.7
Prg max/avg/min : 100028 100028 100028

And using the default Pimple Sover at Time = 0.0892992

rho max/avg/min : 1.15895 1.15862 1.1583
T   max/avg/min : 300.004 300 299.996
P   max/avg/min : 100028 100000 99971.8
Prg max/avg/min : 100029 100029 100028

Implying my change is having no practical effect. I will be trying to figure this out tomorrow.

[pupo]

I've "fixed the issue".

Currently, at Time = 0.0892992 i have:

rho max/avg/min : 1.15862 1.15829 1.15796
T   max/avg/min : 300 300 300
P   max/avg/min : 100000 99971.6 99943.2
Prg max/avg/min : 100000 100000 100000

where the initial soluiton was:

rho max/avg/min : 1.15862 1.15829 1.15796
T   max/avg/min : 300 300 300
P   max/avg/min : 100000 99971.6 99943.2
Prg max/avg/min : 100000 100000 100000

As in, exactly the same. The maximum velocity vector inside the domain has a module of the order of -8. This is hardly significant and can probably be traced to computational errors.

the code I've used to achieve this is rather simple . In the CreateFields.H you have to change the P_rgh calculation for a calculation of the pressure and density according to the ideal Gas Law.

// definition of the Ideal Gas constant. 
specie theSpecie( thermo.subDict("mixture") );
dimensionedScalar RperfectGas ("RperfectGas",dimensionSet(0,2,-2,-1,0,0,0),scalar(theSpecie.R())); 

// p_rgh = p - rho*gh; // comment this out,  it's a less than ideal initial solution
thermo.p()   = p_rgh * exp( gh / (  RperfectGas * thermo.T() )); // New Pressure
thermo.rho() = thermo.p() / ( RperfectGas * thermo.T() );          // New Rho
rho = thermo.rho();    // Necessary for the first iteration of the solver to consider the new density!

you will have to do an include of "species.h" in the main solver code, as well as the libraries in the Make folder.

Additonaly, if you want the a short report of the Density, pressure etc as i've been showing add the following code:

Info<< "rho max/avg/min : " << max(thermo.rho()).value() << " " << average(thermo.rho()).value() << " " << min(thermo.rho()).value() << endl;
Info<< "T   max/avg/min : " << max(thermo.T()).value()   << " " << average(thermo.T()).value()   << " " << min(thermo.T()).value()   << endl;
Info<< "P   max/avg/min : " << max(thermo.p()).value()   << " " << average(thermo.p()).value()   << " " << min(thermo.p()).value()   << endl;
Info<< "Prg max/avg/min : " << max(p_rgh).value()        << " " << average(p_rgh).value()        << " " << min(p_rgh).value()        << endl;
Info<< " " << endl;

I have this code block at the end of the createFields.H and at the end of every iteration.

Good luck!

[jherb]

Just one remark about your last code block: The methods min,average,max are execute on each processor and only the results of processor0 are printed in case of a parallel run. There are functions/macros gMax, gMin, gAverage which are working across processors.

See e. g. https://github.com/OpenFOAM/OpenFOAM...nctions.H#L218

[pupo]

Originally Posted by jherb

See e. g. https://github.com/OpenFOAM/OpenFOAM...nctions.H#L218

Thanks for the sugestion, this works just fine ( you have to remove .value() )

Info<< "rho max/avg/min : " << gMax(thermo.rho()) << " " << gAverage(thermo.rho()) << " " << gMin(thermo.rho()) << endl;
Info<< "T   max/avg/min : " << gMax(thermo.T())   << " " << gAverage(thermo.T())   << " " << gMin(thermo.T())   << endl;
Info<< "P   max/avg/min : " << gMax(thermo.p())   << " " << gAverage(thermo.p())   << " " << gMin(thermo.p())   << endl;
Info<< "Prg max/avg/min : " << gMax(p_rgh)        << " " << gAverage(p_rgh)        << " " << gMin(p_rgh)        << endl;

[Matt_B]

Originally Posted by pupo

I've "fixed the issue".

Currently, at Time = 0.0892992 i have:

rho max/avg/min : 1.15862 1.15829 1.15796
T   max/avg/min : 300 300 300
P   max/avg/min : 100000 99971.6 99943.2
Prg max/avg/min : 100000 100000 100000

where the initial soluiton was:

rho max/avg/min : 1.15862 1.15829 1.15796
T   max/avg/min : 300 300 300
P   max/avg/min : 100000 99971.6 99943.2
Prg max/avg/min : 100000 100000 100000

As in, exactly the same. The maximum velocity vector inside the domain has a module of the order of -8. This is hardly significant and can probably be traced to computational errors.

the code I've used to achieve this is rather simple . In the CreateFields.H you have to change the P_rgh calculation for a calculation of the pressure and density according to the ideal Gas Law.

// definition of the Ideal Gas constant. 
specie theSpecie( thermo.subDict("mixture") );
dimensionedScalar RperfectGas ("RperfectGas",dimensionSet(0,2,-2,-1,0,0,0),scalar(theSpecie.R())); 

// p_rgh = p - rho*gh; // comment this out,  it's a less than ideal initial solution
thermo.p()   = p_rgh * exp( gh / (  RperfectGas * thermo.T() )); // New Pressure
thermo.rho() = thermo.p() / ( RperfectGas * thermo.T() );          // New Rho
rho = thermo.rho();    // Necessary for the first iteration of the solver to consider the new density!

you will have to do an include of "species.h" in the main solver code, as well as the libraries in the Make folder.

Additonaly, if you want the a short report of the Density, pressure etc as i've been showing add the following code:

Info<< "rho max/avg/min : " << max(thermo.rho()).value() << " " << average(thermo.rho()).value() << " " << min(thermo.rho()).value() << endl;
Info<< "T   max/avg/min : " << max(thermo.T()).value()   << " " << average(thermo.T()).value()   << " " << min(thermo.T()).value()   << endl;
Info<< "P   max/avg/min : " << max(thermo.p()).value()   << " " << average(thermo.p()).value()   << " " << min(thermo.p()).value()   << endl;
Info<< "Prg max/avg/min : " << max(p_rgh).value()        << " " << average(p_rgh).value()        << " " << min(p_rgh).value()        << endl;
Info<< " " << endl;

I have this code block at the end of the createFields.H and at the end of every iteration.

Good luck!

Thanks Pedro for your help.
I changed as you here indicated and I got much better results for what is related to p_rgh: now I also have values stuck at nearly 100000 Pa, consistent with the initial conditions:

rho max/avg/min : 1.15859 1.15829 1.158
T max/avg/min : 300 300 300
P max/avg/min : 99997.2 99971.6 99946
Prg max/avg/min : 100000 100000 100000

The problem is that the errors for the velocity components don't allow to get below the set U equation error residuals, PIMPLE conditions that I set and are necessary to pass further to the next 0.1 sec step; I set all equations error below 10^-4. How to resolve this issue?
Below you find the plot of equation residuals:



Ux and Uy it's true they have absolute values very small, below 10^-11, but the equation errors are still very high, close to 1, as you can see.
Any ideas?

Thank you.

Matteo

[pupo]

Hi Matt,

Unfortunately this is something I'm also struggling with, and I have no solution (nor i fully understand the problem)

[Matt_B]

Hi you all,
from last post I made one modification to the code and now the situation is better I think, but not already completely resolved I'm afraid.
The modification is relevant to the header pEqn.H at the first line:

//    rho = thermo.rho();  //<---- original code uncommented out
    thermo.rho() = rho;    // modification

I noticed that with the previous expression, probably the code wasn't taking into account the contribution of the continuity equation which is stated with the expression #include "rhoEqn.H" at the end of pEqn.H.

In this way, now the solver converges very fast to the solution and p_rgh keeps his initial value of 100000 Pa (still considering the modification to be taken in createFields.H, suggested previously by Pedro).
For the simulation I took a time step of 0.1 sec and no thermal contribution, with a constant T=300 K into the domain.



There are still some problems though; for example the velocity, even without any buoyant force in the domain, is once again generated as can be checked in the below image, even if with a smaller value than before:



I checked that, in that point where velocity is generated, the density is higher than the surroundings and has a similar texture distribution in the domain.
An other strange behaviour is that if I take a smaller time step, for example DT=0.0001 sec, it does not converge anymore.

A help on this subject is welcome.

Thank you.

Matteo

[pupo]

hi matt, I've made some progresses in this issue. So, first let me clear up some mistakes i made:

- P_rgh will not be constant in the domain! if you look at the momentum balance equation, and since we expect the velocity to be 0 we are left with an equation like [LaTeX Error: Syntax error]. This means that a really small difference between the P_rgh at the top and at the bottom of the box is expected. This difference is actually small enough that this line of code:

Info<< "Prg max/avg/min : " << max(p_rgh).value()        << " " << average(p_rgh).value()        << " " << min(p_rgh).value()        << endl;

actually truncates the difference and shows Prg max/avg/min : 10000 10000 10000 all the time.

This forced me to make some changes to my createfields.h

specie theSpecie( thermo.subDict("mixture") ); // unchanged
dimensionedScalar RperfectGas ("RperfectGas",dimensionSet(0,2,-2,-1,0,0,0),scalar(theSpecie.R())); // unchanged


thermo.p()  = thermo.p() * Foam::exp( gh / RperfectGas / thermo.T() ); // this will now use the initial p value,
thermo.correct(); // this will update thermo.rho() and all the other thermophysical propertis to the new pressure

volScalarField& p = thermo.p(); // says that the "p" used in the equations is equal to thermo.p()

//....

p_rgh = p - rho*gh; // we are now back to the original openFOAM formulation!

You should remove your change:

//    rho = thermo.rho();  //<---- original code uncommented out
    thermo.rho() = rho;    // modification

since thermo.rho() is updated with thermo.correct(); your changes did help in this particular case where there is no heat, but as soon as we considered heat again the solver would almost for sure start acting weird. It would just keep the density from the first iteration and never update it.

This changes are working for me in an open Case problem. Although i did have to program a new outlet boundary condition for P_rgh. I tried this in your case (but using Simple method) and the velocity also converged to 0 after some time. It is however not ideal and i don't understand why. The solution should not converge to 0 velocity, it should start there and not progress from there...

Hope this helps, and I hope my mistakes earlier did not cause too much of an hassle

PS.: something i noticed, you should disable the radiation models in the no heat scenario, if you are not doing it already. They misbehave badly when there is nothing to heat by radiation and generate spurious temperature peaks.