CFD Online Discussion Forums

CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   OpenFOAM (http://www.cfd-online.com/Forums/openfoam/)
-   -   Negative absolute pressure in rhoPorousMRFSimpleFoam (http://www.cfd-online.com/Forums/openfoam/111290-negative-absolute-pressure-rhoporousmrfsimplefoam.html)

doubtsincfd January 2, 2013 14:55

Negative absolute pressure in rhoPorousMRFSimpleFoam
 
2 Attachment(s)
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!

Tushar@cfd January 4, 2013 07:14

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.

doubtsincfd January 4, 2013 09:16

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

Tushar@cfd January 4, 2013 23:58

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.

doubtsincfd January 7, 2013 14:42

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.

Tushar@cfd January 9, 2013 02:38

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

niklas January 9, 2013 07:51

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

doubtsincfd January 9, 2013 14:36

2 Attachment(s)
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.

niklas January 9, 2013 15:23

3 Attachment(s)
how does these settings work

doubtsincfd January 9, 2013 15:53

1 Attachment(s)
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.

doubtsincfd January 10, 2013 10:29

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?

doubtsincfd January 10, 2013 10:46

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?

niklas January 10, 2013 10:49

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?

niklas January 10, 2013 10:52

and/or your convergence criteria for pressure is too large

doubtsincfd January 10, 2013 12:32

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

niklas January 11, 2013 04:09

Quote:

Originally Posted by doubtsincfd (Post 401240)
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.

doubtsincfd January 11, 2013 10:47

Quote:

Yes, thats what I think.
I cant explain why the static pressure is 101325, other than that you havent converged the pressure
enough.
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

Quote:

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 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

niklas January 13, 2013 13:59

Quote:

Originally Posted by doubtsincfd (Post 401386)
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?

maHein January 13, 2013 17:01

Quote:

Originally Posted by doubtsincfd (Post 401213)
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?

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.

niklas January 14, 2013 09:10

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.


All times are GMT -4. The time now is 06:25.