CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Main CFD Forum

Pressure correction in fractional step method

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

Reply
 
LinkBack Thread Tools Display Modes
Old   May 22, 2012, 03:41
Default
  #21
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 1,651
Rep Power: 23
FMDenaro will become famous soon enough
Quote:
Originally Posted by new_at_this View Post
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 ...
FMDenaro is offline   Reply With Quote

Old   May 22, 2012, 14:44
Default
  #22
Member
 
Join Date: Jul 2011
Posts: 59
Rep Power: 6
rmh26 is on a distinguished road
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?
rmh26 is offline   Reply With Quote

Old   May 23, 2012, 11:39
Default
  #23
Member
 
Peter
Join Date: Oct 2011
Posts: 52
Rep Power: 5
new_at_this is on a distinguished road
Quote:
Originally Posted by FMDenaro View Post
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 View Post
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.
new_at_this is offline   Reply With Quote

Old   May 23, 2012, 12:01
Default
  #24
Member
 
Join Date: Jul 2011
Posts: 59
Rep Power: 6
rmh26 is on a distinguished road
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.
rmh26 is offline   Reply With Quote

Old   May 23, 2012, 13:34
Default
  #25
Member
 
Peter
Join Date: Oct 2011
Posts: 52
Rep Power: 5
new_at_this is on a distinguished road
Quote:
Originally Posted by rmh26 View Post
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
new_at_this is offline   Reply With Quote

Old   May 23, 2012, 13:52
Default
  #26
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 1,651
Rep Power: 23
FMDenaro will become famous soon enough
do you modify congruently near the boundaries the source term in the Poisson equation?
FMDenaro is offline   Reply With Quote

Old   May 23, 2012, 14:20
Default
  #27
Member
 
Join Date: Jul 2011
Posts: 59
Rep Power: 6
rmh26 is on a distinguished road
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.
rmh26 is offline   Reply With Quote

Old   May 23, 2012, 16:13
Default
  #28
Member
 
Peter
Join Date: Oct 2011
Posts: 52
Rep Power: 5
new_at_this is on a distinguished road
Quote:
Originally Posted by rmh26 View Post
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 View Post
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.
new_at_this is offline   Reply With Quote

Reply

Thread Tools
Display Modes

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 Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Comparing between the Fractional step method and the SIMPLE method ghlee Main CFD Forum 1 April 10, 2012 16:59
Time Accuracy of the NITA - FRACTIONAL STEP method Paolo Lampitella FLUENT 1 March 19, 2008 06:45
IcoFoam parallel woes msrinath80 OpenFOAM Running, Solving & CFD 9 July 22, 2007 02:58
free surface bc using pressure correction method Sal Main CFD Forum 0 March 10, 2007 16:50
Gas pressure question Dan Moskal Main CFD Forum 0 October 24, 2002 22:02


All times are GMT -4. The time now is 22:56.