Large pressure gradients from 3D pressure poisson equation

selig5576 · August 1, 2017, 14:50

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)]);
            }
        }

mprinkey · August 1, 2017, 15:02

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)] = ...
    }
  }
}

Without such a structure, the core portions of your field are never updated.

FMDenaro · August 1, 2017, 15:02

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

selig5576 · August 1, 2017, 15:09

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)]);
        }
    }
}

mprinkey · August 1, 2017, 15:33

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.

selig5576 · August 3, 2017, 16:17

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.

FMDenaro · August 3, 2017, 17:08

Quote:

Originally Posted by selig5576

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.

I don't understand what you are doing... diverges and the converges?? It means nothing for me...
Furthermore, what does it mean that you check the continuity using upwind??

selig5576 · August 4, 2017, 10: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.

FMDenaro · August 4, 2017, 11:23

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.

selig5576 · August 4, 2017, 12:52

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,

$\frac{uf^{n+1}_{i,j,k} - uf^{n+1}_{i-1,j,k}}{dx} + \frac{vf^{n+1}_{i,j,k} - vf^{n+1}_{i,j-1}}{dy} + \frac{wf^{n+1}_{i,j,k} - wf^{n+1}_{i,j,k-1}}{dz} = 0.0$

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.

FMDenaro · August 4, 2017, 14:48

This appears a staggered notation, how have you colocated the velocity components with respect to the pressure node?

leflix · August 5, 2017, 10:58

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

FMDenaro · August 5, 2017, 12:16

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.

selig5576 · August 7, 2017, 11:32

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

$uf^{n+1}_{i,j,k} = 0.5*(u^{*}_{i+1,j,k} + u^{*}_{i,j,k}) - \frac{dt}{\rho}\left(\frac{P_{i+1,j,k} - P_{i,j,k}}{dx}\right)$
$vf^{n+1}_{i,j,k} = 0.5*(v^{*}_{i,j+1,k} + v^{*}_{i,j,k}) - \frac{dt}{\rho}\left(\frac{P_{i,j+1,k} - P_{i,j,k}}{dy}\right)$
$wf^{n+1}_{i,j,k} = 0.5*(w^{*}_{i,j,k+1} + w^{*}_{i,j,k}) - \frac{dt}{\rho} \left(\frac{P_{i,j,k+1} - P_{i,j,k}}{dz}\right)$
$\frac{uf^{n+1}_{i,j,k} - uf^{n+1}_{i-1,j,k}}{dx} + \frac{vf^{n+1}_{i,j,k} - vf^{n+1}_{i,j-1,k}}{dy} + \frac{wf^{n+1}_{i,j,k} - wf^{n+1}_{i,j,k-1}}{dz} = 0$

This gives me a more appropriate conservation of mass continuity, i.e. 1e-8. I'm not sure why I forgot to do this

FMDenaro · August 7, 2017, 11:46

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.

selig5576 · August 7, 2017, 11:51

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

FMDenaro · August 7, 2017, 11:57

Quote:

Originally Posted by selig5576

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

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.

selig5576 · August 7, 2017, 13:37

OK, good to know. When I look at my continuity cells, they are very small, which I take as a good sign.

FMDenaro · August 7, 2017, 14:27

Quote:

Originally Posted by selig5576

OK, good to know. When I look at my continuity cells, they are very small, which I take as a good sign.

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.

selig5576 · August 8, 2017, 09:48

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

$P^{n+1} = \phi^{n+1} + (Re*dt/2)*\nabla^{2} \phi^{n+1}$

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.

August 1, 2017, 15:02		#2
mprinkey 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)] = ... } } } Without such a structure, the core portions of your field are never updated.

August 1, 2017, 15:33		#5
mprinkey 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
selig5576 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 4, 2017, 10:10	Clarification	#8
selig5576 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, 12:52	Ppe	#10
selig5576 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, $\frac{uf^{n+1}_{i,j,k} - uf^{n+1}_{i-1,j,k}}{dx} + \frac{vf^{n+1}_{i,j,k} - vf^{n+1}_{i,j-1}}{dy} + \frac{wf^{n+1}_{i,j,k} - wf^{n+1}_{i,j,k-1}}{dz} = 0.0$ 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 1, 2017, 15:02		#3
FMDenaro Senior Member Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,768 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 4, 2017, 11:23		#9
FMDenaro Senior Member Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,768 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, 14:48		#11
FMDenaro Senior Member Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,768 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
leflix 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
FMDenaro Senior Member Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,768 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
selig5576 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 $uf^{n+1}_{i,j,k} = 0.5(u^{}_{i+1,j,k} + u^{}_{i,j,k}) - \frac{dt}{\rho}\left(\frac{P_{i+1,j,k} - P_{i,j,k}}{dx}\right)$ $vf^{n+1}_{i,j,k} = 0.5(v^{}_{i,j+1,k} + v^{}_{i,j,k}) - \frac{dt}{\rho}\left(\frac{P_{i,j+1,k} - P_{i,j,k}}{dy}\right)$ $wf^{n+1}_{i,j,k} = 0.5(w^{}_{i,j,k+1} + w^{*}_{i,j,k}) - \frac{dt}{\rho} \left(\frac{P_{i,j,k+1} - P_{i,j,k}}{dz}\right)$ $\frac{uf^{n+1}_{i,j,k} - uf^{n+1}_{i-1,j,k}}{dx} + \frac{vf^{n+1}_{i,j,k} - vf^{n+1}_{i,j-1,k}}{dy} + \frac{wf^{n+1}_{i,j,k} - wf^{n+1}_{i,j,k-1}}{dz} = 0$ 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
FMDenaro Senior Member Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,768 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. leflix likes this.

August 7, 2017, 11:51	Continuity	#16
selig5576 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, 13:37	Ppe	#18
selig5576 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 8, 2017, 09:48	Last question	#20
selig5576 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 $P^{n+1} = \phi^{n+1} + (Redt/2)\nabla^{2} \phi^{n+1}$ 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