# Pressure correction in fractional step method

 Register Blogs Members List Search Today's Posts Mark Forums Read

May 22, 2012, 03:41
#21
Senior Member

Filippo Maria Denaro
Join Date: Jul 2010
Posts: 3,360
Rep Power: 38
Quote:
 Originally Posted by new_at_this 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 ...

 May 22, 2012, 14:44 #22 Member   Join Date: Jul 2011 Posts: 59 Rep Power: 8 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?

May 23, 2012, 11:39
#23
Member

Peter
Join Date: Oct 2011
Posts: 52
Rep Power: 7
Quote:
 Originally Posted by FMDenaro 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 ...
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.

Quote:
 Originally Posted by rmh26 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?
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.

 May 23, 2012, 12:01 #24 Member   Join Date: Jul 2011 Posts: 59 Rep Power: 8 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.

May 23, 2012, 13:34
#25
Member

Peter
Join Date: Oct 2011
Posts: 52
Rep Power: 7
Quote:
 Originally Posted by rmh26 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

 May 23, 2012, 13:52 #26 Senior Member   Filippo Maria Denaro Join Date: Jul 2010 Posts: 3,360 Rep Power: 38 do you modify congruently near the boundaries the source term in the Poisson equation?

 May 23, 2012, 14:20 #27 Member   Join Date: Jul 2011 Posts: 59 Rep Power: 8 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.

May 23, 2012, 16:13
#28
Member

Peter
Join Date: Oct 2011
Posts: 52
Rep Power: 7
Quote:
 Originally Posted by rmh26 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?

Quote:
 Originally Posted by FMDenaro do you modify congruently near the boundaries the source term in the Poisson equation?
I'm sorry, but I'm not exactly sure what you mean.

 Thread Tools Display Modes Linear Mode

 Posting Rules You may not post new threads You may not post replies You may not post attachments You may not edit your posts BB code is On Smilies are On [IMG] code is On HTML code is OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post ghlee Main CFD Forum 1 April 10, 2012 16:59 Paolo Lampitella FLUENT 1 March 19, 2008 06:45 msrinath80 OpenFOAM Running, Solving & CFD 9 July 22, 2007 02:58 Sal Main CFD Forum 0 March 10, 2007 16:50 Dan Moskal Main CFD Forum 0 October 24, 2002 22:02