|
[Sponsors] |
Large pressure gradients from 3D pressure poisson equation |
|
LinkBack | Thread Tools | Search this Thread | Display Modes |
August 1, 2017, 14:50 |
Large pressure gradients from 3D pressure poisson equation
|
#1 |
Senior Member
Selig
Join Date: Jul 2016
Posts: 213
Rep Power: 10 |
I have recently starting making my numerical solver 3D from 2D and I have ran into a problem when solving my PPE. For my 2D code, I have a convergent solution and pressure field never exceeds values of 1.2. For my 3D solver after 1 iteration I get pressure values of 12000.
I am almost certain it has to do with my pressure solver as my velocity components are not reaching values even near those of the pressure. The method in which I am solving the PPE is SOR. My relaxation factor I am using is 1.93 and my tolerance is 1e-6. Also my max iterations is 30,000. Code:
for (int j=2; j < ny-2; j++) { for (int k=2; k < nz-2; k++) { R[ID(1,j,k)] = (rho/dt)*(idx*(ufs[ID(1,j,k)] - u[ID(0,j,k)]) + idy*(vfs[ID(1,j,k)] - vfs[ID(1,j-1,k)]) + idz*(wfs[ID(1,j,k)] - wfs[ID(1,j,k-1)])); P[ID(1,j,k)] = (1.0-omega)*P[ID(1,j,k)] + Aw*(idx2*P[ID(2,j,k)] + idy2*(P[ID(1,j+1,k)] + P[ID(1,j-1,k)]) + idz2*(P[ID(1,j,k+1)] + P[ID(1,j,k-1)]) - R[ID(1,j,k)]); R[ID(nx-2,j,k)] = (rho/dt)*(idx*(u[ID(nx-1,j,k)] - ufs[ID(nx-3,j,k)]) + idy*(vfs[ID(nx-2,j,k)] - vfs[ID(nx-2,j-1,k)]) + idz*(wfs[ID(nx-2,j,k)] - wfs[ID(nx-2,j,k-1)])); P[ID(nx-2,j,k)] = (1.0-omega)*P[ID(nx-2,j,k)] + Ae*(idx2*P[ID(nx-3,j,k)] + idy2*(P[ID(nx-2,j+1,k)] + P[ID(nx-2,j,j-1)]) + idz2*(P[ID(nx-2,j,k+1)] + P[ID(nx-2,j,k-1)]) - R[ID(nx-2,j,k)]); } } for (int i=2; i < nx-2; i++) { for (int k=2; k < nz-2; k++) { R[ID(i,1,k)] = (rho/dt)*(idx*(ufs[ID(i,1,k)] - ufs[ID(i-1,1,k)]) + idy*(vfs[ID(i,1,k)] - v[ID(i,0,k)]) + idz*(wfs[ID(i,1,k)] - wfs[ID(i,1,k-1)])); P[ID(i,1,k)] = (1.0-omega)*P[ID(i,1,k)] + An*(idx2*(P[ID(i+1,1,k)] + P[ID(i-1,1,k)]) + idy2*P[ID(i,2,k)] + idz2*(P[ID(i,1,k+1)] + P[ID(i,1,k-1)]) - R[ID(i,1,k)]); R[ID(i,ny-2,k)] = (rho/dt)*(idx*(ufs[ID(i,ny-2,k)] - ufs[ID(i-1,ny-2,k)]) + idy*(v[ID(i,ny-1,k)] - vfs[ID(i,ny-3,k)]) + idz*(wfs[ID(i,ny-2,k)] - wfs[ID(i,ny-2,k-1)])); P[ID(i,ny-2,k)] = (1.0-omega)*P[ID(i,ny-2,k)] + As*(idx2*(P[ID(i+1,ny-2,k)] + P[ID(i-1,ny-2,k)]) + idy2*P[ID(i,ny-3,k)] + idz2*(P[ID(i,ny-2,k+1)] + P[ID(i,ny-2,k-1)]) - R[ID(i,ny-2,k)]); } } for (int i=2; i < nx-2; i++) { for (int j=2; j < ny-2; j++) { R[ID(i,j,1)] = (rho/dt)*(idx*(ufs[ID(i,j,1)] - ufs[ID(i-1,j,1)]) + idy*(vfs[ID(i,j,1)] - vfs[ID(i,j-1,1)]) + idz*(wfs[ID(i,j,1)] - w[ID(i,j,0)])); P[ID(i,j,1)] = (1.0-omega)*P[ID(i,j,1)] + Ab*(idx2*(P[ID(i+1,j,1)] + P[ID(i-1,j,1)]) + idy2*(P[ID(i,j+1,1)] + P[ID(i,j-1,1)]) + idz2*P[ID(i,j,2)] - R[ID(i,j,1)]); R[ID(i,j,nz-2)] = (rho/dt)*(idx*(ufs[ID(i,j,nz-2)] - ufs[ID(i-1,j,nz-2)]) + idy*(vfs[ID(i,j,nz-2)] - vfs[ID(i,j-1,nz-2)]) + idz*(w[ID(i,j,nz-1)] - wfs[ID(i,j,nz-3)])); P[ID(i,j,nz-2)] = (1.0-omega)*P[ID(i,j,nz-2)] + At*(idx2*(P[ID(i+1,j,nz-2)] + P[ID(i-1,j,nz-2)]) + idy2*(P[ID(i,j+1,nz-2)] + P[ID(i,j-1,nz-2)]) + idz2*P[ID(i,j,nz-3)] - R[ID(i,j,nz-2)]); } } |
|
August 1, 2017, 15:02 |
|
#2 |
Senior Member
Michael Prinkey
Join Date: Mar 2009
Location: Pittsburgh PA
Posts: 363
Rep Power: 25 |
None of your updates touch the inner nodes. These all seem to loop over oddly chosen 2D boundary planes. The first block loops over j,k. The second over i,k, and the last over i,j. Somewhere, you need to loop over i,j,k all together. I don't know exactly your notation here, but you should have a loop that runs:
Code:
for (int i =...) { for (int j = ...) { for (int k = ...} { R[ID(i,j,k)] = ... } } } |
|
August 1, 2017, 15:02 |
|
#3 |
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,777
Rep Power: 71 |
for incompressible flow, the absolute value of the pressure has no relevance...check if your pressure gradients are correct, that is if they ensure a divergence-free velocity field
|
|
August 1, 2017, 15:09 |
Pressure
|
#4 |
Senior Member
Selig
Join Date: Jul 2016
Posts: 213
Rep Power: 10 |
mprinkey: My indexing is perfectly fine. What I am doing is writing a 1D array and then extracting the X,Y,Z dimensions from the 1D array, this is a very common method as having a 3D array as a 1D array is more efficient than P[nx][ny][nz] or vector<vector<vector>>>. In terms of accessing the data of my pressure array, I do indeed loop over the interior P and PPE RHS (R). I just did not include it in the code snippet. With that said here is what i am doing.
Note: #define ID (i, j, k) (int) i + ny*j + ny*nz*k. Code:
void PressureSolver (dVector& us, dVector& vs, dVector& ws, dVector& u, dVector& v, dVector& w, dVector& P) { dVector ufs, vfs, wfs; dVector Residual, R; const double omega = 1.53; const double tol = 1e-6; const int max_it = 65000; double L2_Error = 0.0; int it = 0.0; const double Ad = omega/(2.0*idx2 + 2.0*idy2 + 2.0*idz2); const double Aw = omega/(idx2 + 2.0*idy2 + 2.0*idz2); const double Ae = omega/(idx2 + 2.0*idy2 + 2.0*idz2); const double An = omega/(2.0*idx2 + idy2 + 2.0*idz2); const double As = omega/(2.0*idx2 + idy2 + 2.0*idz2); const double At = omega/(2.0*idx2 + 2.0*idy2 + idz2); const double Ab = omega/(2.0*idx2 + 2.0*idy2 + idz2); ufs.resize(nxyz); vfs.resize(nxyz); wfs.resize(nxyz); R.resize(nxyz); Residual.resize(nxyz); FaceVelocities(us, vs, ws, ufs, vfs, wfs); for (int i=2; i < nx-2; i++) { for (int j=2; j < ny-2; j++) { for (int k=2; k < nz-2; k++) { R[ID(i,j,k)] = (rho/dt)*(idx*(ufs[ID(i,j,k)] - ufs[ID(i-1,j,k)]) + idy*(vfs[ID(i,j,k)] - vfs[ID(i,j-1,k)]) + idz*(wfs[ID(i,j,k)] - wfs[ID(i,j,k-1)])); } } } it = 0; while (it <= max_it) { for (int i=2; i < nx-2; i++) { for (int j=2; j < ny-2; j++) { for (int k=2; k < nz-2; k++) { P[ID(i,j,k)] = (1.0-omega)*P[ID(i,j,k)] + Ad*(idx2*(P[ID(i+1,j,k)] + P[ID(i-1,j,k)]) + idy2*(P[ID(i,j+1,k)] + P[ID(i,j-1,k)]) + idz2*(P[ID(i,j,k+1)] + P[ID(i,j,k-1)]) - R[ID(i,j,k)]); } } } for (int j=2; j < ny-2; j++) { for (int k=2; k < nz-2; k++) { R[ID(1,j,k)] = (rho/dt)*(idx*(ufs[ID(1,j,k)] - u[ID(0,j,k)]) + idy*(vfs[ID(1,j,k)] - vfs[ID(1,j-1,k)]) + idz*(wfs[ID(1,j,k)] - wfs[ID(1,j,k-1)])); P[ID(1,j,k)] = (1.0-omega)*P[ID(1,j,k)] + Aw*(idx2*P[ID(2,j,k)] + idy2*(P[ID(1,j+1,k)] + P[ID(1,j-1,k)]) + idz2*(P[ID(1,j,k+1)] + P[ID(1,j,k-1)]) - R[ID(1,j,k)]); R[ID(nx-2,j,k)] = (rho/dt)*(idx*(u[ID(nx-1,j,k)] - ufs[ID(nx-3,j,k)]) + idy*(vfs[ID(nx-2,j,k)] - vfs[ID(nx-2,j-1,k)]) + idz*(wfs[ID(nx-2,j,k)] - wfs[ID(nx-2,j,k-1)])); P[ID(nx-2,j,k)] = (1.0-omega)*P[ID(nx-2,j,k)] + Ae*(idx2*P[ID(nx-3,j,k)] + idy2*(P[ID(nx-2,j+1,k)] + P[ID(nx-2,j,j-1)]) + idz2*(P[ID(nx-2,j,k+1)] + P[ID(nx-2,j,k-1)]) - R[ID(nx-2,j,k)]); } } for (int i=2; i < nx-2; i++) { for (int k=2; k < nz-2; k++) { R[ID(i,1,k)] = (rho/dt)*(idx*(ufs[ID(i,1,k)] - ufs[ID(i-1,1,k)]) + idy*(vfs[ID(i,1,k)] - v[ID(i,0,k)]) + idz*(wfs[ID(i,1,k)] - wfs[ID(i,1,k-1)])); P[ID(i,1,k)] = (1.0-omega)*P[ID(i,1,k)] + An*(idx2*(P[ID(i+1,1,k)] + P[ID(i-1,1,k)]) + idy2*P[ID(i,2,k)] + idz2*(P[ID(i,1,k+1)] + P[ID(i,1,k-1)]) - R[ID(i,1,k)]); R[ID(i,ny-2,k)] = (rho/dt)*(idx*(ufs[ID(i,ny-2,k)] - ufs[ID(i-1,ny-2,k)]) + idy*(v[ID(i,ny-1,k)] - vfs[ID(i,ny-3,k)]) + idz*(wfs[ID(i,ny-2,k)] - wfs[ID(i,ny-2,k-1)])); P[ID(i,ny-2,k)] = (1.0-omega)*P[ID(i,ny-2,k)] + As*(idx2*(P[ID(i+1,ny-2,k)] + P[ID(i-1,ny-2,k)]) + idy2*P[ID(i,ny-3,k)] + idz2*(P[ID(i,ny-2,k+1)] + P[ID(i,ny-2,k-1)]) - R[ID(i,ny-2,k)]); } } for (int i=2; i < nx-2; i++) { for (int j=2; j < ny-2; j++) { R[ID(i,j,1)] = (rho/dt)*(idx*(ufs[ID(i,j,1)] - ufs[ID(i-1,j,1)]) + idy*(vfs[ID(i,j,1)] - vfs[ID(i,j-1,1)]) + idz*(wfs[ID(i,j,1)] - w[ID(i,j,0)])); P[ID(i,j,1)] = (1.0-omega)*P[ID(i,j,1)] + Ab*(idx2*(P[ID(i+1,j,1)] + P[ID(i-1,j,1)]) + idy2*(P[ID(i,j+1,1)] + P[ID(i,j-1,1)]) + idz2*P[ID(i,j,2)] - R[ID(i,j,1)]); R[ID(i,j,nz-2)] = (rho/dt)*(idx*(ufs[ID(i,j,nz-2)] - ufs[ID(i-1,j,nz-2)]) + idy*(vfs[ID(i,j,nz-2)] - vfs[ID(i,j-1,nz-2)]) + idz*(w[ID(i,j,nz-1)] - wfs[ID(i,j,nz-3)])); P[ID(i,j,nz-2)] = (1.0-omega)*P[ID(i,j,nz-2)] + At*(idx2*(P[ID(i+1,j,nz-2)] + P[ID(i-1,j,nz-2)]) + idy2*(P[ID(i,j+1,nz-2)] + P[ID(i,j-1,nz-2)]) + idz2*P[ID(i,j,nz-3)] - R[ID(i,j,nz-2)]); } } for (int i=2; i < nx-2; i++) { for (int j=2; j < ny-2; j++) { for (int k=2; k < nz-2; k++) { Residual[ID(i,j,k)] = Laplace(P, i, j, k) - R[ID(i,j,k)]; } } } L2_Error = 0.0; for (int i=2; i < nx-2; i++) { for (int j=2; j < ny-2; j++) { for (int k=2; k < nz-2; k++) { L2_Error += abs(dx*dy*dz*Residual[ID(i,j,k)]); } } } if (L2_Error < tol) { break; } else { it += 1; } } for (int j=0; j < ny; j++) { for (int k=0; k < nz; k++) { P[ID(0,j,k)] = P[ID(1,j,k)] + (dx/dt)*(us[ID(0,j,k)] - u[ID(0,j,k)]); P[ID(nx-1,j,k)] = P[ID(nx-2,j,k)] - (dx/dt)*(us[ID(nx-1,j,k)] - u[ID(nx-1,j,k)]); } } for (int i=0; i < nx; i++) { for (int k=0; k < nz; k++) { P[ID(i,0,k)] = P[ID(i,1,k)] + (dy/dt)*(vs[ID(i,0,k)] - v[ID(i,0,k)]); P[ID(i,ny-1,k)] = P[ID(i,ny-2,k)] - (dy/dt)*(vs[ID(i,ny-1,k)] - v[ID(i,ny-1,k)]); } } for (int i=0; i < nx; i++) { for (int j=0; j < ny; j++) { P[ID(i,j,0)] = P[ID(i,j,1)] + (dz/dt)*(ws[ID(i,j,0)] - w[ID(i,j,0)]); P[ID(i,j,nz-1)] = P[ID(i,j,nz-2)] - (dz/dt)*(ws[ID(i,j,nz-1)] - w[ID(i,j,nz-1)]); } } } |
|
August 1, 2017, 15:33 |
|
#5 |
Senior Member
Michael Prinkey
Join Date: Mar 2009
Location: Pittsburgh PA
Posts: 363
Rep Power: 25 |
Your first code snippet didn't include the i,j,k relaxation loop.
EDIT. Does the pressure system converge or does the pressure solver exit because it hit the maximum number of iterations? Have you monitored convergence? Have you tried a smaller value for omega? It is 1.53 in the code, but you mention 1.93 in your first post. Last edited by mprinkey; August 1, 2017 at 21:15. |
|
August 3, 2017, 16:17 |
Pressure
|
#6 |
Senior Member
Selig
Join Date: Jul 2016
Posts: 213
Rep Power: 10 |
After doing some testing. I deduced that for my 3D solver, after 1 time step it diverges and then starts "converging" after the second step. Checking continuity (using upwind FDM) I see that some values are large, i.e. -20.0896, etc. Does anyone have suggestions? I can post code.
|
|
August 3, 2017, 17:08 |
|
#7 | |
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,777
Rep Power: 71 |
Quote:
Furthermore, what does it mean that you check the continuity using upwind?? |
||
August 4, 2017, 10:10 |
Clarification
|
#8 |
Senior Member
Selig
Join Date: Jul 2016
Posts: 213
Rep Power: 10 |
Sorry for the unclear message I wrote the post in a rush, apologies. What I meant to say is to calculate the continuity equation I used 1st-order upwind method for the derivatives. Subsequently, looking at the values of the array I'm getting rather large numbers, which I know to be incorrect as they are have values 2.09, 20.8983, etc.
In terms of the pressure, I'm finding that with a tolerance of 1e-6 and an maximum iterations of 25000, after the first iteration of the PPE (it converges after 435 time steps), however the second step never converges. It runs through the maximum iterations. I am 95% sure my problems are due to the PPE, but I am not sure why with a tolerance of 1e-6, it is not working. If I pick a tolerance of 1e-5, it properly solves the PPE. |
|
August 4, 2017, 11:23 |
|
#9 |
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,777
Rep Power: 71 |
I suspect a bug in the code that affects the convergence....
However, you cannot use an upwind discretization for the Div v term. You must use exactly the same discrete operators used in the pressure equation, otherwise the check for the discrete continuity error has no meaning. |
|
August 4, 2017, 12:52 |
Ppe
|
#10 |
Senior Member
Selig
Join Date: Jul 2016
Posts: 213
Rep Power: 10 |
Hi,
I managed to "fix" the bug with my solver, so it accepts a tolerance of 1e-6. Thank you for pointing that out, I was not aware. I am doing the discretization of the face velocities at the future time step n+1, Upon fixing my tolerance bug, I am also seeing smaller pressure values as I saw in my convergent 2D solver. In terms of the conservation of mass, I am seeing much smaller numbers, but after the first time step it will sum to -29.23. After each time step it will get more negative, i.e. after 100 time steps it will get to -10,043.2, which is not right as it should be 10^-4 or in that ballpark. |
|
August 4, 2017, 14:48 |
|
#11 |
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,777
Rep Power: 71 |
This appears a staggered notation, how have you colocated the velocity components with respect to the pressure node?
|
|
August 5, 2017, 10:58 |
|
#12 |
Senior Member
Join Date: Aug 2011
Posts: 272
Rep Power: 15 |
Hi Selig,
I do agree with mprinkey what you do is whatever you want but certainly not SOR. You have to loop on the 3 indices i, j and k Leflix |
|
August 5, 2017, 12:16 |
|
#13 |
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,777
Rep Power: 71 |
I do not understand if the posted code is complete or not...
however, I suggest to test always the new coded subroutines in a controlled test case, not in the NS solver. For example, prescribing as source term a sum of the elements in the row to check for a solution constant equal to 1. |
|
August 7, 2017, 11:32 |
Fixed problem
|
#14 |
Senior Member
Selig
Join Date: Jul 2016
Posts: 213
Rep Power: 10 |
Hi everyone,
thank you for the replies I believe I was able to fix my problem. So one reason why I was not getting "proper" conservation of mass is I was not using the right terms in the continuity, i.e. I was not using the RC interpolation properly. When I use This gives me a more appropriate conservation of mass continuity, i.e. 1e-8. I'm not sure why I forgot to do this |
|
August 7, 2017, 11:46 |
|
#15 |
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,777
Rep Power: 71 |
What you wrote is not exactly the RC interpolation but only the Hodge decomposition in discrete form. If you substitute the terms in the discrete continuity equation that you wrote at bottom, you simply obtain the pressure equation you solve.
|
|
August 7, 2017, 11:51 |
Continuity
|
#16 |
Senior Member
Selig
Join Date: Jul 2016
Posts: 213
Rep Power: 10 |
I understand that the the face velocities will form the PPE, but is what I'm doing for the continuity not right? I did some reading and found that for proper continuity, one should construct the continuity equation from the face velocities (n+1)..
|
|
August 7, 2017, 11:57 |
|
#17 | |
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,777
Rep Power: 71 |
Quote:
Yes it's right, and the correct process is to write the continuity equation in the form you wrote and then substitute the face velocities at n+1 using the Hodge decomposition. That explain also how to set the proper BC.s. Of course, after the pressure equation is solved, you can check the continuity equation that is satisfied in each cell up to the tolerance of the pressure equation. |
||
August 7, 2017, 13:37 |
Ppe
|
#18 |
Senior Member
Selig
Join Date: Jul 2016
Posts: 213
Rep Power: 10 |
OK, good to know. When I look at my continuity cells, they are very small, which I take as a good sign.
|
|
August 7, 2017, 14:27 |
|
#19 |
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,777
Rep Power: 71 |
Of course, if you check the continuity constraint by using the velocities in the cell centre, they do not satisfy the divergence-free constraint up to the tolerance of the pressure equation. The error will be of the order of magnitude of the local truncation error which is a typical issue in the approximate projection method.
|
|
August 8, 2017, 09:48 |
Last question
|
#20 |
Senior Member
Selig
Join Date: Jul 2016
Posts: 213
Rep Power: 10 |
I'm not sure if this deserves a new thread, but a quick question:
I am using Runge-Kutta for both advection and diffusion terms, and as I am using a pressure-free projection method according(roughly) to Kim and Moins paper (1985) they do their final pressure update with It is mentioned it defined like that to accomodate for the Crank-Nicholson for the viscous terms. Given I am not doing implicit integration, do I still need that type of pressure update? In terms of my intermediate velocity BC, I am using Kim and Moin's formulation properly fulfill the condition to second-order. |
|
|
|
Similar Threads | ||||
Thread | Thread Starter | Forum | Replies | Last Post |
Calculation of the Governing Equations | Mihail | CFX | 7 | September 7, 2014 06:27 |
Problems with pressure equation in SIMPLER algorithm | Corby | Main CFD Forum | 0 | June 17, 2013 09:24 |
The correction on pressure equation of SIMPLE algorithm in MRFSimpleFOAM solver | renyun0511 | OpenFOAM Running, Solving & CFD | 0 | November 10, 2010 01:47 |
Neumann pressure BC and velocity field | Antech | Main CFD Forum | 0 | April 25, 2006 02:15 |
FEM pressure poisson equation (Implicit) | cfd-beginner | Main CFD Forum | 0 | August 9, 2005 13:32 |