Pressure correction in fractional step method
 

Can someone explain to me how the pressure correction step works in cylindrical coordinates? I have been working on a code for a while and my pressure correction step is not producing a divergence free field.

did you check if your method is divergence-free in Cartesian coordinates? fractional step is a method to split pressure and velocity in the time integration, the system of coordinate is not involved

How does the divergence change before and after the correction step?
Is the divergence localized in certain regions after the correction step?
How accurately is your pressure equation being solved?
Are you using an axisymmetric cylindrical grid, is weird stuff happening at the axis of symmetry?

I do have a working code in cartesian coordinates that has been verified against ghia for the lid driven cavity. I'm a bit confused about what you are saying... Are you saying I should use the exact same pressure correction step for cylindrical coordinates even though the divergence is different?

I am using an axisymmetric grid to solve for the internal flow in a cylinder. Before the pressure correction, the maximum divergence in the field is on the order of 10^-1. After the correction step, it is on the order of 10^-2. After the pressure correction step, the majority of the field has a divergence on the order of 10^-3 and the max divergence values are located in the upper left and right corners near the cylinder wall. In my 2D code I get divergence values orders of magnitude smaller.

I initially start with a velocity profile of 1 everywhere in the field. When integrated it produces a value of 1. I expect that at any time, if I integrate the velocity profile, I should also get a value of 1. However as my flow progresses to steady state, I get values up to 1.2

Are you using an exact or approximate projection method?
You must be careful in the inner-outer product of discretized divergence-gradient operators

I believe that i am using an exact projection method. In each step of my code I update the u_z and u_r field by the advection and diffusion terms. Then I solve the following equation



Could you clarify more about the inner outer divergence-gradient operators?

The pressure equation must be discretized from the Div Grad (*) operators, you should not use the Laplace operator. This allows you to substitute the correct BC in terms of the normal component of the pressure gradient

So instead of solving for the pressure in one step, could I split it into two steps? the first step would be multiplying the divergence of the intermediate velocity (not divergent free) by the inverse of the discretized divergence matrix. Then multiply that result by the inverse of the discretized gradient matrix.

Also are there any books or papers that I can look up that talk about how to do this?

for those interested, I found this paper which essentially talks about my problem. Unfortunately I still have not found the bug in my code. I expect my code to produce a poiseuielle like profile that peaks at 1.5 however my code peaks at ~1.8 and is not divergence free.

If anyone has seen similar problems, please let me know. thanks

http://www.sciencedirect.com/science...45793004001148

How do you arrive at this equation?? If it is a pressure correction equation then it should have a 'delta-t' term somewhere!!!

maybe the time step is simply inserted in the P variable ...

the equation is just the cylindrical expansion for:


I didn't put in the time because it cancels once applied in the navier stokes equations.

I have checked and I don't think that is the problem.

is not a problem if you then compute

Vn+1= V* - Grad P

In any case, I am quite sure you have problems in setting pressure BCs ... you must write Div (Grad P) and substitute n.Grad P at boundaries

I think I have accounted for that. Since my code is on a staggered grid, I'm only calculating pressure at the interior points. At the edges, I assume that the ghost point outside the boundary is the same as the first interior point. For example the second derivative in the R direction should be . Substituting for neumann bcs I get at the axis and at the wall.

I do not have time right now to figure this out. But I will try to imbibe what you wrote.

I am assuming the P in your equation is pressure correction and not the pressure itself.

I asked that because I never remember deriving pressure correction in fractional step method that is independent of delta T.

You can not apply zero gradient at the axis by ghost points. Axis has pressure gradients.


Yes this is the case, he is calculating for ( deltT * Pressure-correction).

There normal component of the gradient of the pressure should be zero at the axis, if we are assuming axisymmetric flow and the fact that Ur must be zero at the axis.

What are you using for the inlet/oulet conditions?

So i'm working on a sliding cylinder problem, which is axisymmetric. Basically its a cylindrical cavity where the end caps are stationary and the wall is moving with a given velocity in the z direction. So expect that the pressure gradient at the axis should be zero. I've defined U_r to be 0 at all the boundaries. U_z has a velocity of 1 at the outer wall and symmetry condition at the axis.

If anyone is interested or if anyone thinks seeing the code might help, then i can post my matlab code. I'm pretty sure that this poisson equation solver is the last bug in my code.

I think I figured out where the error is. When I formulate the finite difference matrix for the laplacian, I essentially solve Ax = b. However A is nearly singular and so my solution ends up being innacurate. How can I increase the accuracy?

be careful ... for Neumann conditions, the matrix is singular but you have to fulfill the compatibility relation (that involves BCs) such that the system admits infinite solutions...this is an issue very different from the accuracy in solving the system ...

Yeah the nuemann condition part is can be tricky. I think a lot of people set one of the boundary values to zero instead of using a nuemann condition. You mentioned that you solved this problem in rectangular coordinates. Did you need any special treatment for the neumann condition then?

Can you elaborate or point me towards a reference that talks about this more. I have looked into a couple of cfd references and they don't really seem to talk about this issue and in most cases only talk about cartesian cases.

In the cartesian code I don't think I did anything special. For my Laplacian finite difference matrix, I replaced the ghost pressure point outside the boundary with the first interior point due to nuemann BCs, as i described earlier in the thread. To solve the singular matrix i used the gmres function in matlab that makes use of the GMRES method and solves the near singular linear system.

In cylindrical coordinates, I have made the Laplacian matrix, except in contrast to the cartesian version it isn't symmetric but is still near singular. I have done the same thing as before and replaced pressure outside the boundary with the first interior point for all the boundary finite differences. Then I try solving the system Ax = b where A is the near singular discrete laplacian matrix, and b is the divergence of the intermediate velocity field. I use the same GMRES method that I used in the cartesian case but I still produce a field that has a divergence on the order of 1E-2 to 1E-3 when it should be much lower.

Well first off if you are solving a poisson type problem with nueman conditions it is 100% singular not just nearly. Try solving a 1-D poisson problem with nuemann conditions using gaussian elimination and you'll see what i mean.

Suggestions
1. Examine the residual r = Ax-b. See if it is maximum where your divergence errors are max.
2. Look at the pressure field itself to see if you really are achieving a nuemann condition on the pressure.
3. Post your MATLAB code, meshplots of pressure, divergence etc.

1. So I checked the residual and yes the residual is maxed out in the same regions as the maximum divergence values.

2. I don't ever calculate the pressure points on the boundary or outside because I am computing on a staggered grid. However based on my formulation of the finite difference matrix I believe that this is enforced.

3. My code can be found in the .zip folder in the link below. If you run the script cartesian_sliding_cavity.m, it will compute the cartesian equivalent to my problem. If you run axisymmetric_sliding_cylinder.m, it will run the code I am currently trying to debug. If anything is unclear, please ask since it will force me to justify my actions and perhaps help me better understand what I am doing wrong. Thanks for all the help

https://www.dropbox.com/s/1tiuq4aobu...ric%20Code.zip

do you modify congruently near the boundaries the source term in the Poisson equation?

Program is surprisingly short, I guess all the built in MATLAB functions help but its also makes it very tough to read/understand.

Anyways
% pressure correction
rhs = reshape(diff([uW;U;uE])/hx+diff([vS' V vN']')'/hy+avg([vS' V vN']')'./rvp,[],1);
[p,flag(k)] = gmres(-Lp,rhs,[],[],length(Lp));
.
.
.
div = div_cyl(nx,ny,hx,hy,Ue,Ve,X);
maxdiv(k) = max(max(abs(div(2:end-1,2:end-1))));


It looks like you are using two different methods of calculating the divergence. One way is using diff and that is fed into the poisson solver then you use div_cyl to check the divergence. If you check the diveregence after with the diff method it apperas to be ok.

You are right. I didn't consider that. I guess what confuses me right now is that I expect that the area under the velocity curve should be constant in time. However as time progresses the area under my velocity curve increases. I thought that it was a divergence problem, but could it be due to something else?

I'm sorry, but I'm not exactly sure what you mean.