[doubtsincfd]

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]

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]

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]

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]

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]

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]

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]

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]

how does these settings work

[doubtsincfd]

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]

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]

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]

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]

and/or your convergence criteria for pressure is too large

[doubtsincfd]

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]

Originally Posted by doubtsincfd

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]

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

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]

Originally Posted by doubtsincfd

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]

Originally Posted by doubtsincfd

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]

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.