[Pflaume1891]

Hello, I have a question aboute the pressure correction using the SIMPLE Algorithm.

Velocity correction is defined as:
v' =Df * grad (p')

Df = Volume / aP

aP is the pivot matrix coefficient of the momentum equation.

Imagine a hexaeder with equal edge length of 0.1, density 1.0 and viscosity of 1.0e-5

In the first iteration, massflow for each face is 0.0.

The coefficiert of each neighbour element in my matrix line is:

aE =(1.0/Edgelength) * Area * Viscosity = 1.0e-6.

As the hexaeder has six sides, the pivot coefficient is:

aP = 6 * aE = 6.0e-6

That means for Df:

Df = Volume / aP = 0.001 / 6.0e-6 = 166!

Multipling Df with a pressure correction gradient of 1.0, I get a velocity correction value of -166m/s

I have this problem in a cartesian channel with length of 1 m. After the first iteration I get a linear profile from the inlet to the outlet with -1.0 Pa/m for the pressure correction and here we go...

The pressure value makes sense, because i have a poisson eqation with the value 0.0 at the outlet, zero flux at inlet and wall.
At the first row of cells behind the inlet, i have a source of 0.01 with a inlet velocity of 1.0 m/s.

Can any body tell me, where I may be wrong?

Thank you in advance.

[Pflaume1891]

Hello.

I have some plots, to show my problem.

I want to simulate a elbow with my incompressible FVM Solver using SIMPLE Pressure Correction and Rhie Chow interpolation.

With this pictures, I want to illustrate my problem with the velocity correction. All this picture where taken after the first iteration.

The first picture shows the central matrix coefficient of the momentum equation. Sorry, "Insert Image" function didn't work.

https://goo.gl/photos/9au32vBMANxGLx2CA


This is shows the pressure correction after the first iteration:
https://goo.gl/photos/SUEoK5vrs5ntpdFn7


And the gradient of the pressure correction:
https://goo.gl/photos/u4nPVW9yvbV1MrH36


The cell volume:
https://goo.gl/photos/1XxD62XpzirVJCvn6


And DC = Volume / Matrix Coefficient
https://goo.gl/photos/TrQF9byJY4qkgFue9


And this is the problem, because as long as there is no massflow, the matrix coefficient only contains the diffusive fluxes which are multiplied with a viscosity of 1.0e-005 very small.

Finally I get this extrem velocity correction:
https://goo.gl/photos/6QSjtvrmUhLWxDEc7


What could be wrong?

[Pflaume1891]

Mayby, somebody don't want to click at any link, so I attached the pictures of the previous post in this message. As I can only attach five files, the other files are posted in the next message.

The pictures are in the same order as linked in the previous post.

The matrix coefficient at the inlet is much higher than displayed on the legend, because of the convective flux at the inlet. This has not reached any of the other cells in the first iteration.

[Pflaume1891]

And here the missing velocity correction, which is much to high for 1 m/s inlet velocity.

[arjun]

Originally Posted by Pflaume1891

And here the missing velocity correction, which is much to high for 1 m/s inlet velocity.

After multiplying with pressure under relaxation factor of say 0.1 or so correction drops by a magnitude.

Now if pressure is corrected by p = p_old + urf * p_correction, should velocity be also corrected from under relaxed pressure correction gradient. urf x grad (p correction).

22 in stead of 220 does not sound bad for 1m/sec inlet. Note that in the start what you are getting is indeed the situation. Just that velocity correction shall be corresponding to the pressure correction made to pressure and not according to pressure correction you got out of linear solver.

Hope it makes sense.

[Pflaume1891]

Thank you, but I'm still underrelaxing my pressure correction. Usually, I'm using 0.2 for the update. But the velocity correction is based on the pressure correction, not on the corrected pressure, isn't it?

[arjun]

Originally Posted by Pflaume1891

Thank you, but I'm still underrelaxing my pressure correction. Usually, I'm using 0.2 for the update. But the velocity correction is based on the pressure correction, not on the corrected pressure, isn't it?

See the point is that pressure is update as p new = p old + urf x p correction

So the pressure is corrected by (urf x p correction) and not ( p correction) so the corresponding change in velocity shall be based on (urf x p correction) and not based on
( p correction) to be consistent (otherwise you may face convergence issues).

What you are doing now is that velocity correction = -vol/ap x grad(p corr) , which should be the case is you accepted all the pressure correction that came out of linear system.

Another way to think is this, if your pressure urf was 0, then there shall be no change in velocities due to pressure correction. Isn't it.

[Pflaume1891]

Ok, that means

v' = DC * grad(p')

is correct, but more stable is

v' = DC * uRelax * grad(p').

Did I get this right?

[arjun]

Originally Posted by Pflaume1891

Ok, that means

v' = DC * grad(p')

is correct, but more stable is

v' = DC * uRelax * grad(p').

Did I get this right?

It is not issue of more stable.


The equation
v' = DC * grad(p')

Assuming that everything of p correction is accepted. Which is not the case.
What is accepted is

v' = DC * grad( urf . p')

Think over it.

[Pflaume1891]

I just made the change. Now, I'm calculating p' = urelax * p'. This way, grad (urelax * p') is calculated. Well, the velocity overshots are less, but still velocity correction is about 71m/s. Without velocity correction, I have no overshots and get a maximum velocity of about 2 m/s. At the outlet, I'm calculating the gradient of p' using zero at the outlet face and there, I always have grad (p')>>0 so the velocity there is exploding after some iterations

[Pflaume1891]

I'm not really sure about my pressure. My solver needs 1 s to solve three times the momentum equation, but 10 s to solve the pressure correction equation.

I also calculated this case in ANSYS Fluent, incompressible, laminar. The attached pictures are showing the results after the first iteration.

The pressure field is only visible at the inlet. And the velocities are looking like my central matrix coefficients.

The pressure coefficient gradient is looking like the final results of the flow - fast at the inner radius and slow at the outer.

[arjun]

Originally Posted by Pflaume1891

I just made the change. Now, I'm calculating p' = urelax * p'. This way, grad (urelax * p') is calculated. Well, the velocity overshots are less, but still velocity correction is about 71m/s. Without velocity correction, I have no overshots and get a maximum velocity of about 2 m/s. At the outlet, I'm calculating the gradient of p' using zero at the outlet face and there, I always have grad (p')>>0 so the velocity there is exploding after some iterations

See you are in start of the calculation and it is pretty much the scenario. From what you have written here, you dont seem to have bug.

Solver should recover from here in next iterations, that is the corrections become smaller and smaller.

[Pflaume1891]

If I want to use static pressure for outlet condition, do I need to ensure global mass balance in the pressure equation? Actually, I'm calculating outlet massflow using rhie chow for boundary faces with extrapolated values and correct it. I thing, watching the fluent pictures, that my pressure corrrction is much to large. In fluent, a underrelaxation factor of 0.3 was used.

[arjun]

Originally Posted by Pflaume1891

If I want to use static pressure for outlet condition, do I need to ensure global mass balance in the pressure equation? Actually, I'm calculating outlet massflow using rhie chow for boundary faces with extrapolated values and correct it. I thing, watching the fluent pictures, that my pressure corrrction is much to large. In fluent, a underrelaxation factor of 0.3 was used.

Fluent does lot of other stuff to keep pressure correction "nice" and velocity correction under control. This is why places where your code may diverge their would work. Dont expect that you will be writing exact same thing as you and me dont know what they exactly do.

I thought you wanted to understand if there was a bug in your code, it seems there is no bug.

[Pflaume1891]

Well, that is not really what happens. The correction is increasing.

Attached, a youtube link to a video, what happens:
https://www.youtube.com/watch?v=AHug2wF-xFM
Finally, the velocity correction is 1.4e004 m/s!!!
Velocity is updated using: v' = DC * grad(urelax * p')

If i disable pressure correction equation, means pressure correction is alwaly zero, I get a well looking convective transport:
https://www.youtube.com/watch?v=QenBt_KbzzI

Sorry, I have no numbers at the legend yet. The first video is ranged from min to max - you can see the values in the left, lower corner.
In the second video, the range is from 0 m/s to 1 m/s as i used a inlet velocity of 1 m/s

[arjun]

You should try with something more simpler set up that avoids skew as much as possible.

It is very difficult to find the issue from here, but velocity shooting up in start is quite normal. Carefully inspect, sometimes what you think you are doing and what actually done is different.

[Pflaume1891]

Originally Posted by arjun

You should try with something more simpler set up that avoids skew as much as possible.

It is very difficult to find the issue from here, but velocity shooting up in start is quite normal. Carefully inspect, sometimes what you think you are doing and what actually done is different.

I can exclude some error sources, because I'm only using functions, which are debugged. I never calculate the velocity gradient matrix at a face or cell center - I always use a function GetFaceCenterGradientMatric(TCell) or something like that.
The diffusive coefficients have to be correct, because the wall distance is calculated correctly.

My program has a mathematica interface where I checked all the parameter used to calculate the fluxes.

I uploaded a cartesian channel case. It's the same.
https://youtu.be/TPVbhE8_hp8

The eye-catching thing is, that the velocity at the outlet is increasing extremely. Here I update the massflow by:
m = m* + m' = m* + Density * FluxCoefficient * (p'_Outletface - p'_Cell)

[Pflaume1891]

The point is, that there are some things I don't know if I do correctly.

How to calculate the massflow at the outlet for MOM and Pressure equation?
Rhie Chow, Velocity extrapolation or something else?

Can i underrelax the pressure correction equation it self or just the correction?

Is my pressure field correct, because mass imbalace doen't decrease in the outer loops.

If there are more unknown settings, it's hard to find out...

Ps.: If i increase viscosity to 1.0 Pa s - the solution is swining into its final solution in a complette other way,
because of the higher matrix central coefficient.

[arthursoprano]

You could try to use the under-relaxation factor for the Momentum equations, like the following



where is the under-relaxation factor and see if it helps.

[arjun]

Originally Posted by Pflaume1891

The point is, that there are some things I don't know if I do correctly.

How to calculate the massflow at the outlet for MOM and Pressure equation?
Rhie Chow, Velocity extrapolation or something else?

Can i underrelax the pressure correction equation it self or just the correction?

Is my pressure field correct, because mass imbalace doen't decrease in the outer loops.

If there are more unknown settings, it's hard to find out...

Ps.: If i increase viscosity to 1.0 Pa s - the solution is swining into its final solution in a complette other way,
because of the higher matrix central coefficient.

It is very hard to tell from here what is wrong in your code as it very hard to tell even when the code is in front of you.

I once wrote LES code and upon testing for flow around sphere every thing looked very good. except one small thing Cp at front was 2 and not 1 as it shall be.
I wrote the whole code again from scratch because I thought that is faster thing to do than to find for this bug.

Carefully go step by step.

1. At out flow even if you left Rhie and chow it would not hurt. I have run simulations where I completely ignored rhie and chow flux so leaving it on outlet is not a big deal.

2. You dont under relax pressure equation by scaling diagonal term. You update pressure using urf as I mentioned earlier.

3. You pressure field is definitely wrong as your continuity is diverging.

This is why I said try a simpler case first. At the moment your case has skew and all and its hard to tell if it is diverging due to bug or it diverged due to some other reason. (Like its pointed out that you may not be under relaxing momentum).