Negative absolute pressure in rhoPorousMRFSimpleFoam
 

Hi,

I am trying to solve a simple problem of compressible flow through a 2D rectangular duct with a porous zone in between.

The domain is described in the attached pdf

I ran 150 iterations and I am getting negative absolute pressure in the domain, which is non-physical. I am monitoring pressure at inlet. The 0 constant and system folders are in the attached tar file.

run the script "./runIt". Observe the values of inlet pressure in log file.

How can absolute pressure go negative at all? What is going wrong in such a simple problem?

Thank you!

Hello,

It seems the solver is having problem with the boundary condition. In your case, why don't you put the BC for the pressure as;
internalField uniform 0;
outlet
{
type fixedValue;
value uniform 0;
}
Since, at the inlet boundary of your problem case you are defining the Velocity(U).

Also, if you want to use the same BC which you have mentioned in your case, then try changing the BC for velocity(U) too.

Thanks, Tushar. The solution diverges in the first iteration itself when I put zero on outlet and internalField. It gives a floating point exception.

Hello friend,

I rechecked your case, Have you changed the "controlDict" file??

Since, I am not an expert but, it seems the function that you might have added in it, could be responsible for the divergence of pressure.

Do, run the solver without the function and recheck the solution.

Yes I have changed the controlDict file. But as per my understanding, the additional piece of code I have written is only for post processing purpose.

Removing that piece of code did not affect the solution.

Hello friend,
 

How you are defining "porosityWall" in your case? I think the problem lies somewhere with "porosityWall" you might have mistaken some parameters.
You need to check those with the available case models.
I wish this will solve your problem.

-----
Best regards,
Tushar

you have a massflow as inlet bc and you are running a 2D case,
so my question is then if the value of the massflow is specified in kg/s/m or in kg/s
and if the width of the channel is 1 meter or just 'thin' because that will matter.

remember that the area of the inlet will be the actual area used when calculating the inlet velocity

Thanks, Niklas. You are absolutely right. I wanted to simulate 1 meter depth, but the depth in blockMesh I gave was 10mm. I was not aware of the fact that the actual area of the inlet will be used by OpenFOAM for a 2D case.

Now I have increased the depth to 1 meter and the mass flow rate I have entered is in kg/s.

Actually I am trying to simulate a canonical problem for a bigger 3D problem that I am trying to solve. My original problem is same (compressible+porous+high flow rates but subsonic flow with Ma < 0.4), but it is 3D. Unfortunately I am not allowed to share the data of that problem.

I have increased the mass flow rate to 0.5 kg/s. Upto 0.5 kg/s, there is no negative absolute pressure. But as I increase the flow rate, the negative pressure problem starts becoming more and more prominent.

For 0.5 kg/s the solution stabilizes, and I get a convergence. But for higher flow rates there is instability, and in some cases it causes the solution to blow up.

I am not able to understand this negative absolute pressure, and the instability at the start of the solution. I simulated the same problem in Fluent, and the solution is smooth. The OpenFOAM results and the Fluent results are with 1% difference.

I am attaching the updated files for simulation. Script to run is "runIt.sh". Also, I am attaching plots of the convergence history.

how does these settings work

Thanks for the settings.

The damping for pressure is amazing. Attached jpeg file shows residuals for first few iterations.

I will try these setting for my main problem.

I would like to know the reason behind the change in settings.

Also, the problem of negative pressure at the start of iterations still exists (though eventually it disappears). Though this does not have effect on final solution since it is a steady state problem, I am not able to understand how the equation of state is getting satisfied with negative absolute pressure. The domain does not show negative absolute temperature anywhere.

Hi Niklas,

I have an important question about the boundary condition at outlet. For pressure, you have suggested totalPressure boundary condition at outlet with the value of p0 =101325 Pa. But when I post-process, it shows static pressure = 101325 Pa at the outlet. So is OpenFOAM really fixing the static pressure at outlet instead of total pressure? Why is this happening despite of the value specified for total pressure?

So I looked at the totalPressureFvPatchScalarField.C and the code is:

operator==(p0p - 0.5*(1.0 - pos(phip))*magSqr(Up));

My reading of the equation:
pos(1) function returns 1, since flow is leaving the domain. This make the term - 0.5*(1.0 - pos(phip))*magSqr(Up) go to zero and so static pressure value is equal to total pressure value.

But, I am not able to understand what p0p does and what the = operator does. Can someone please check my explanation and help me understand p0p and = operator?

Im going to make a wild guess here and assume that your massflow is too low to really make a difference.

What happens if you increase the massflow?

and/or your convergence criteria for pressure is too large

Thank you, Niklas.

If I have understood you correctly, do you mean to say that for higher mass flow rates and/or tighter convergence criterion for pressure, the static pressure might deviate from the desired value?

For the current flow conditions, inlet mass flow rate = 2kg/s
Velocity at outlet = 160 m/s
Static pressure at outlet = 101325 Pa
Total pressure at outlet (cal from ptot utility) = 106703 Pa
Maximum mach number = 0.285

For higher flow rates (> 4kg/s) and with current settings, the solution diverges after 2 iterations

Yes, thats what I think.
I cant explain why the static pressure is 101325, other than that you havent converged the pressure
enough.

About the divergence, try to reduce the under-relaxation for rho even more and observe...
What that low under-relaxation I would run more than 2000 iteration just to see if there is
something changing slowly after that.

For the same simulation (2kg/s), I ran the simulation for 6000 iterations, and here are my observations:

1. Residuals < 1e-06
2. Inlet pressure value remains constant for >2000 iterations
3. Outlet pressure value rock solid at 101325 for all iterations suggesting that the solver is not solving for static pressure at outlet

For higher flow rate [> 5kg/s], I reduced the k-epsilon relaxation to 0.2 and the simulation ran without divergence. URFs for p and rho were the same

ok, btw how do you look at the pressure at the outlet?

Your explanation is correct! If the flow enters the domain, p0 is interpreted as total pressure (for this reason the dynamic pressure is subtracted from p0 to get the static pressure). Is the flow leaving the domain, p0 is used as static pressure. That's why you will always get a static pressure of 101325 Pa at the outlet.

The =-operator just defines what happens when the boundary conditions are updated. So in this case, the static pressure at the outlet is calculated using the above equation.

Hmm, I've never thought about this but its correct.

Using total pressure as outlet bc is most likely unstable considering that you can
reduce the pressure and increase the velocity, while maintaining the value of p0.

use fixedValue instead.

Thank you, Niklas and maHein.

I monitored the outlet pressure using this code:
pOutlet
{
type patchAverage;
functionObjectLibs ("libsimpleFunctionObjects.so");
verbose true;
patches (outlet);
log true;
fields
(
p
);
}

then its incorrect to set totalPressure at outlet?
I have 2 totalPressure BC at inlet and outlet.at inlet it behaves somewhat good(except supersonic velocity at 3-4 begining time steps that should be subsonic in real but it modifies at later time steps and becomes subsonic) but when outlet opens pressure goes too little and velocity increases a lot and problem diverges.
could anyone argue about the formula I have written at other thread also?
thanks.

dear Omkar I have a foolish question!:D
how do you plot residuals and othe variables like you did above?
I want to plot residuals by iteration numbers like you did but pyFoam draws residuals( and not variables at all) by time step in x-axis.

can anyone guide me through two questions i propounded at two posts above?

You can plot residuals using gnuplot as mentioned in this post:

http://www.cfd-online.com/Forums/ope...tml#post410537

I am not sure about using total pressure at both ends. When you specify total pressure at both ends, it is like fixing the loss across the system rather than letting the system decide. This is what I think, I am not sure about it though.

Plus, as per the discussions above, if the flow is leaving the domain, OF will considered the input value as static pressure because the term 1.0 - pos(phip) will go to zero.

I was looking at the files posted by you and those modified by Niklas to understand the reason of the improvement. When you say "The damping for pressure is amazing" are you referring to the decrease in tolerance from 1e-6; to in 1e-9; the fvSolution?
...
p
{
solver GAMG;
tolerance 1e-9;
relTol 0.1;

No. What I meant was that the pressure values stabilized quickly.

how can avoid high velocities at first iterations?maybe limiting of velocity to mach velocity?i tried but without useful result.
Efforts are thanked

has anyone done a progress on the issue of high velocities at first few iterations in total pressure model of BC?