[FMDenaro]

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 ...

[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?

[new_at_this]

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.

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.

[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.

[new_at_this]

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

[FMDenaro]

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

[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.

[new_at_this]

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?

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.