Incremental projection method

moshe · January 24, 2018, 12:39

So I've been working on trying to get a non-incremental projection method working, but to no avail. I've narrowed it do the pressure update step which I don't seem to quite understand how to implement.

In the predictor step we have $\nabla P^{n}$ . We then solve a PPE for $\phi^{n+1}$ . With this in hand we do the FINAL pressure step update as $P^{n+1} = P^{n} + \phi^{n+1}$

Shouldnt it be as simple as

Code:

while (t < tend)
% solve us, vs
% solve PPE for phi_new
% project onto divergence free velocity field
% time up date velocities

%update pressure
Pnew(1:nx,1:ny) = P(1:nx,1:ny) + phi_new(1:nx,1:ny)
P(1:nx,1:ny) = Pnew(1:nx,1:ny)
t = t + dt
end

I feel foolish for asking

. My velocity field is free of oscillations, but there are many oscillations in my pressure field. Given both methods (pressure-free and incremental) are first-order accurate they should (according to literature) have agreeable results.....

FMDenaro · January 24, 2018, 12:46

Are you still working on the lid-driven cavity? Are you using a compact stencil for the pressure? Could you show the pressure field? Add a further plot, the (div v) field obtained after the first time step.

moshe · January 24, 2018, 13:01

1. Yes, I am still working on the lid driven cavity for simplicity.
2. I am using a compact stencil (Rhie-Chow interpolation for coupling)
3. I have included some photos of the divergence and pressure

Some extra information to help you help me

. $\nabla \cdot \mathbf{u}$ is quite small 10^-4 - 10^-3 which would be consistent with the fact I am dealing with an approximate projection method. My code is quite small and very readable if that might help (not that I am asking for debugging help, but I am certain everything except the PPE works)

FMDenaro · January 24, 2018, 13:10

The pressure plot seems quite regular apart some points...what about the grid? However, the div v plot reaches peak values of almost O(10^-1) that is larger than the magnitude of the local truncation error, what about the time-step value?

I suggest to do a first test excluding totally the convective term (very small Re number), that will produce a symmetric field so we can better test the pressure

moshe · January 24, 2018, 13:41

My time step size is 1e-3 with a 32 x 32 grid. Here is a plot of my pressure gradient without the convection term and viscosity of 0.1. As you can see there are still oscillations. Note after a few time steps it gets more unstable. In my pressure-free PM I do not have this effect.

FMDenaro · January 24, 2018, 14:00

I can suspect your dt is high and drives to a numerical instability. Try to use dt=10^5. Then, what about the level of convergence of the pressure equation? are you using a low threshold on the residual?

moshe · January 24, 2018, 17:35

Even if I set my dt = 1e-4 I still oscillations in my solution. In terms of my tolerance, My tolerance constraint is 10^-6. I have also tried 10^-9 and even 10^-12 and I still get the oscillations. I have experimented with both the L1 and L2 norm and theres not too much of a difference except the L1 needs less iterations to convergence than L2.

FMDenaro · January 24, 2018, 18:14

Do oscillations appear just after the first time step? What about if you refine the mesh? I would think some bugs in the interpolation or in the bc.s for the pressure. Could you write your pressure equation in a node close to a wall or a corner?

moshe · January 24, 2018, 21:48

I feel bad for not providing more details. I did perform a mesh refinement. It did not help. I always start wit ha 32 x 32, nu = 0.01, etc as my base benchmark. In terms of my PPE at the boundary. I have for the left boundary and right boundary

Code:

%Left BC
R(2,j) = (rho/dt)*((1/dx)*(ufs(2,j) - un(1,j)) + (1/dy)*(vfs(2,j) - vfs(2,j-1)));
P(2,j) = (1-omega)*P(2,j) + omega/(1/dx^2 + 2/dy^2)*((1/dx^2)*P(3,j) + (1/dy^2)*(P(2,j+1) + P(2,j-1)) - R(2,j));

%Right BC
R(nx-1,j) = (rho/dt)*((1/dx)*(un(nx,j) - ufs(nx-2,j)) + (1/dy)*(vfs(nx-1,j) - vfs(nx-1,j-1)));
P(nx-1,j) = (1-omega)*P(nx-1,j) + omega/(1/dx^2 + 2/dy^2)*((1/dx^2)*P(nx-2,j) + (1/dy^2)*(P(nx-1,j+1) + P(nx-1,j-1)) - R(nx-1,j));

The R arrays is the RHS of the PPE eand P is pressure. In terms of my linear solver I am starting off simple and using SOR. In terms of the mathematics I am doing

$R_{nx-1,j} = \frac{\rho}{dt} \left(\frac{u^{n+1}_{nx,j} - uf^{*}_{nx-2,j}}{dx} + \frac{vf^{*}_{nx-1,j}-vf_{nx-1,j-1}^{*}}{dy}\right)$

$P_{nx-1,j} = (1-\omega)P_{nx-1,j} + \frac{\omega}{1/dx^{2} + 1/dy^{2}} (\frac{1}{dx^{2}}P_{nx-2,j} + \frac{1}{dy^{2}} \left(P_{nx-1,j+1} + P_{nx-1,j-1} \right)$ $- R_{nx-1,j})$

I apply this condition 4 boundaries and 4 corners. I've gone through my code again and I don't see any issues. Maybe another set of eyes might help

FMDenaro · January 25, 2018, 04:43

I am not sure about your equation and the fact that the pressure is colocated at the grid nodes. I assume you have i=1,j=1 and i=N,j=M as walls index. Now for example I write the x-component in i=2 using a compact stencil

(dp/dx|i+1/2 - dp/dx|i-1/2)/dx = (rho/dt)*(u*|i+1/2-u*|i-1/2)/dx

while the Neumann BC at i=1 is

dp/dx|1=(rho/dt)*(u_wall - u*|1)

W............X.............X -> i direction
1.............2.............3

Is that your pressure colocation? If yes, how do you use the BC into dp/dx|i-1/2?

FMDenaro · January 25, 2018, 05:45

And why the y derivative of vfs is first order non-centered?

moshe · January 25, 2018, 10:30

To add some more details, here is how I do it in full. First off $uf_{i,j},vf_{i,j}$ represent the face velocities in which they are defined as $uf_{i,j} = 0.5(u_{i+1,j} + u_{i,j}), vf_{i,j} = 0.5(v_{i,j+1}+v_{i,j})$

The way I define the RHS of the PPE is

$R_{i,j} \frac{uf_{i,j}^{*} - uf^{*}_{i-1,j}}{dx} + \frac{vf^{*}_{i,j} - vf^{*}_{i,j-1}}{dy}$

My Pressure BC is x (which is defined in a similar fashion for y)

$\frac{P_{i+1,j} - P_{i,j}}{dx} = \frac{uf^{n+1}_{i,j} - uf^{*}_{i,j}}{dt}$

In short I modify the RHS to take into account u^{n+1} and then remove $\frac{\partial P}{\partial x}$ from $\frac{\partial^{2} P}{\partial x^{2}}$
In terms of your question on why the y direction is not centered, it is because vf is the face so doing the algebra we have

$\frac{vf_{i,j} - vf_{i,j-1}}{dy} = \frac{v_{i,j+1} - v_{i,j-1}}{2dy}$

moshe · January 25, 2018, 10:40

My complete PPE routine is

Code:

%=====Construct face velocites====
    for j=1:ny
        for i=1:nx-1
            ufs(i,j) = 0.5*(us(i+1,j) + us(i,j));
        end
    end
      
    for j=1:ny-1
        for i=1:nx
            vfs(i,j) = 0.5*(vs(i,j+1) + vs(i,j));
        end
    end
    %==================================
    
    %Right-hand side of PPE
    for j=3:ny-2
        for i=3:nx-2
            R(i,j) = (rho/dt)*((1/dx)*(ufs(i,j) - ufs(i-1,j)) + (1/dy)*(vfs(i,j) - vfs(i,j-1)));
        end
    end
      
    Iter = 0;
      
    while (Iter <= IterMax)       
       %Iterate interior of P       
        for j=3:ny-2
            for i=3:nx-2
                P(i,j) = (1-omega)*P(i,j) + omega/(2/dx^2 + 2/dy^2)*((P(i+1,j) + P(i-1,j))/dx^2 + (P(i,j+1) + P(i,j-1))/dy^2 - R(i,j));
            end
        end
         
       %Left and right BC
        for j=3:ny-2
            R(2,j) = (rho/dt)*((1/dx)*(ufs(2,j) - un(1,j)) + (1/dy)*(vfs(2,j) - vfs(2,j-1)));
            P(2,j) = (1-omega)*P(2,j) + omega/(1/dx^2 + 2/dy^2)*((1/dx^2)*P(3,j) + (1/dy^2)*(P(2,j+1) + P(2,j-1)) - R(2,j));
      
            R(nx-1,j) = (rho/dt)*((1/dx)*(un(nx,j) - ufs(nx-2,j)) + (1/dy)*(vfs(nx-1,j) - vfs(nx-1,j-1)));
            P(nx-1,j) = (1-omega)*P(nx-1,j) + omega/(1/dx^2 + 2/dy^2)*((1/dx^2)*P(nx-2,j) + (1/dy^2)*(P(nx-1,j+1) + P(nx-1,j-1)) - R(nx-1,j));
        end
   
        %Top and bottom BC
        for i=3:nx-2
            R(i,2) = (rho/dt)*((1/dx)*(ufs(i,2) - ufs(i-1,2)) + (1/dy)*(vfs(i,2) - v(i,1)));
            P(i,2) = (1-omega)*P(i,2) + omega/(2/dx^2 + 1/dy^2)*((1/dx^2)*(P(i+1,2) + P(i-1,2)) + (1/dy^2)*P(i,3) - R(i,2));
     
            R(i,ny-1) = (rho/dt)*((1/dx)*(ufs(i,ny-1) - ufs(i-1,ny-1)) + (1/dy)*(v(i,ny) - vfs(i,ny-2)));
            P(i,ny-1) = (1-omega)*P(i,ny-1) + omega/(2/dx^2 + 1/dy^2)*((1/dx^2)*(P(i+1,ny-1) + P(i-1,ny-1)) + (1/dy^2)*P(i,ny-2) - R(i,ny-1));
        end
        
       %Corner BCs
        R(2,ny-1) = (rho/dt)*((1/dx)*(ufs(2,ny-1) - un(1,ny-1)) + (1/dy)*(vn(2,ny) - vfs(2,ny-1)));
        P(2,ny-1) = (1-omega)*P(2,ny-1) + omega/(1/dx^2 + 1/dy^2)*((1/dx^2)*P(3,ny-1) + (1/dy^2)*P(2,ny-2) - R(2,ny-1));
            
        R(nx-1,2) = (rho/dt)*((1/dx)*(u(nx,2) - ufs(nx-2,2)) + (1/dy)*(vfs(nx-1,2) - vn(nx-1,1)));
        P(nx-1,2) = (1-omega)*P(nx-1,2) + omega/(1/dx^2 + 1/dy^2)*((1/dx^2)*P(nx-2,2) + (1/dy^2)*P(nx-1,3) - R(nx-1,2));
            
        R(nx-1,ny-1) = (rho/dt)*((1/dx)*(un(nx,ny-1) - ufs(nx-2,ny-1)) + (1/dy)*(vn(nx-1,ny) - vfs(nx-1,ny-2)));
        P(nx-1,ny-1) = (1-omega)*P(nx-1,ny-1) + omega/(1/dx^2 + 1/dy^2)*((1/dx^2)*P(nx-2,ny-1) + (1/dy^2)*P(nx-1,ny-2) - R(nx-1,ny-1));
         
        R(2,2) = (rho/dt)*((1/dx)*(ufs(2,2) - un(1,2)) + (1/dy)*(vfs(2,2) - vn(2,1)));
        P(2,2) = (1-omega)*P(2,2) + omega/(1/dx^2 + 1/dy^2)*((1/dx^2)*P(3,2) + (1/dy^2)*P(2,3) - R(2,2));

        %Define residual: Residual = B - Ax        
        Res = 0.0;
         
        for j=3:ny-2
            for i=3:nx-2
                Res(i,j) = R(i,j) - ((1/dx^2)*(P(i-1,j) - 2*P(i,j) + P(i+1,j)) + (1/dy^2)*(P(i,j-1) - 2*P(i,j) + P(i,j+1)));
            end
        end
          
        %Define norm.  In this case we use L1 norm
        normE = dx*dy*abs(sum(sum(Res(3:nx-2,3:ny-2))));
        error = normE;
          
        if (error >= tol)
            Iter = Iter + 1;
        else
            fprintf("PPE has Converged")
            break;
        end
    end 
     
    %After PPE has converged we set the BCs for i=1,nx, and j=1,ny
    P(1,:) = P(2,:) + (dx/dt)*(us(1,:) - u(1,:));
    P(nx,:) = P(nx-1,:) - (dx/dt)*(us(nx,:) - u(nx,:));
    P(:,1) = P(:,2) + (dy/dt)*(vs(:,1) - v(:,1));
    P(:,ny) = P(:,ny-1) - (dy/dt)*(vs(:,ny) - v(:,ny));

FMDenaro · January 25, 2018, 11:47

I tried to figure it out how do you use the BC after you have discretized but I see problems...for example, what is the goal to express the ghost point for the pressure when you have a compact stencil that does not involve it at all?? The pressure equation is written from i,j=2,2 to nx-1,ny-1, right? Therefore the pressure values involved are at i=1 and nx and j=1 and ny.
I don't think that the problem is in the coding but in the formulation you are adopting close to the walls.
I will try to work on and write some extended lines that I will post later.

FMDenaro · January 25, 2018, 12:26

I just wrote an idea how to work, hope it does ...

moshe · January 26, 2018, 10:36

So I spent some time yesterday rewriting my solver with a staggered grid and found some interesting things!

First off, I was able to get the incremental projection method to work on the staggered grid!

. For my collocated based solver, I was using central interpolation to form the gradient of p^{n} in the intermediate step. In the staggered solver I used one sided interpolation, i.e.
$\nabla_{x} q^{n} \approx \frac{q_{E} - q_{P}}{\Delta x}$
$\nabla_{y} q^{n} \approx \frac{q_{N} - q_{P}}{\Delta y}$

When I used a central approximation I would get severe oscillations (numerical instabilities), but this seemed to fix the problem, at least for the staggered grid. I know using a one sided approximation is invalid for the collocated formulation as u,v,w,P are at the same location.

In literature, it is not mentioned what type of approximation they use for the q^{n} pressure time level, but it would make sense why a central approximation gave bad results for the staggered formulation......

FMDenaro · January 26, 2018, 10:46

Quote:

Originally Posted by moshe

So I spent some time yesterday rewriting my solver with a staggered grid and found some interesting things!

First off, I was able to get the incremental projection method to work on the staggered grid!

. For my collocated based solver, I was using central interpolation to form the gradient of p^{n} in the intermediate step. In the staggered solver I used one sided interpolation, i.e.
$\nabla_{x} q^{n} \approx \frac{q_{E} - q_{P}}{\Delta x}$
$\nabla_{y} q^{n} \approx \frac{q_{N} - q_{P}}{\Delta y}$

When I used a central approximation I would get severe oscillations (numerical instabilities), but this seemed to fix the problem, at least for the staggered grid. I know using a one sided approximation is invalid for the collocated formulation as u,v,w,P are at the same location.

In literature, it is not mentioned what type of approximation they use for the q^{n} pressure time level, but it would make sense why a central approximation gave bad results for the staggered formulation......

I am not sure to understand... q is now the pressure? But if you write the pressure equation on staggered grid, you have no difference from the pressure stencil when compared to the compact one used in colocated grid. The pressure gradient are centered at step h at half distance between pressure nodes. The only difference appears in the source term.

moshe · January 26, 2018, 10:52

Sorry, I was a bit sloppy in my variables let me clarify. So in terms of the pressure time levels we have P^{n+1}, P^{n}, and q^{n+1}. In the intermediate step we have P^{n}, we solve the PPE for q^{n+1} and we do the pressure update with P^{n+1}. Again, sorry for the confusion.

In terms of central vs one-sided approximations, I think the pressure in the intermediate step has to be consistent with the corrector-step, though I do not have a rigorous justification yet.

FMDenaro · January 26, 2018, 11:05

Considering the node (i,j) for the q field (colocated at the same grid locations of p^n+1 and p^n), you write

d/dx(dq/dx)|i,j = (dq/dx|i+1/2,j - dq/dx|i-1/2,j)/dx
d/dy(dq/dy)|i,j = (dq/dy|i,j+1/2 - dq/dy|i,j-1/2)/dy

Therefore I don't se any one-sided discretization....You see the effect of the staggering when computing the divergence of the v* field

moshe · January 26, 2018, 11:11

I agree with what you say. Though, I have found a paper that does an incremental PM method (with a staggered grid) and they( use one sided approximations for P^{n} in the intermediate velocity field.

For a collocated grid, the one-sided does not work (I tried it for the sake of curiosity), and it would make sense.

* Calculating three-dimensional flows around structures and over rough terrain, C.W. Hirt, J.L Cook, JCP, 1972

January 24, 2018, 12:39	Results of incremental projection method	#21
moshe Member Moshe De Leon Join Date: Nov 2017 Location: Portugal Posts: 31 Rep Power: 8	So I've been working on trying to get a non-incremental projection method working, but to no avail. I've narrowed it do the pressure update step which I don't seem to quite understand how to implement. In the predictor step we have $\nabla P^{n}$ . We then solve a PPE for $\phi^{n+1}$ . With this in hand we do the FINAL pressure step update as $P^{n+1} = P^{n} + \phi^{n+1}$ Shouldnt it be as simple as Code: while (t < tend) % solve us, vs % solve PPE for phi_new % project onto divergence free velocity field % time up date velocities %update pressure Pnew(1:nx,1:ny) = P(1:nx,1:ny) + phi_new(1:nx,1:ny) P(1:nx,1:ny) = Pnew(1:nx,1:ny) t = t + dt end I feel foolish for asking . My velocity field is free of oscillations, but there are many oscillations in my pressure field. Given both methods (pressure-free and incremental) are first-order accurate they should (according to literature) have agreeable results.....

January 24, 2018, 17:35	Stability	#27
moshe Member Moshe De Leon Join Date: Nov 2017 Location: Portugal Posts: 31 Rep Power: 8	Even if I set my dt = 1e-4 I still oscillations in my solution. In terms of my tolerance, My tolerance constraint is 10^-6. I have also tried 10^-9 and even 10^-12 and I still get the oscillations. I have experimented with both the L1 and L2 norm and theres not too much of a difference except the L1 needs less iterations to convergence than L2.

January 24, 2018, 21:48	Incremental PM	#29
moshe Member Moshe De Leon Join Date: Nov 2017 Location: Portugal Posts: 31 Rep Power: 8	I feel bad for not providing more details. I did perform a mesh refinement. It did not help. I always start wit ha 32 x 32, nu = 0.01, etc as my base benchmark. In terms of my PPE at the boundary. I have for the left boundary and right boundary Code: %Left BC R(2,j) = (rho/dt)((1/dx)(ufs(2,j) - un(1,j)) + (1/dy)(vfs(2,j) - vfs(2,j-1))); P(2,j) = (1-omega)P(2,j) + omega/(1/dx^2 + 2/dy^2)((1/dx^2)P(3,j) + (1/dy^2)(P(2,j+1) + P(2,j-1)) - R(2,j)); %Right BC R(nx-1,j) = (rho/dt)((1/dx)(un(nx,j) - ufs(nx-2,j)) + (1/dy)(vfs(nx-1,j) - vfs(nx-1,j-1))); P(nx-1,j) = (1-omega)P(nx-1,j) + omega/(1/dx^2 + 2/dy^2)((1/dx^2)P(nx-2,j) + (1/dy^2)(P(nx-1,j+1) + P(nx-1,j-1)) - R(nx-1,j)); The R arrays is the RHS of the PPE eand P is pressure. In terms of my linear solver I am starting off simple and using SOR. In terms of the mathematics I am doing $R_{nx-1,j} = \frac{\rho}{dt} \left(\frac{u^{n+1}_{nx,j} - uf^{}_{nx-2,j}}{dx} + \frac{vf^{}_{nx-1,j}-vf_{nx-1,j-1}^{}}{dy}\right)$ $P_{nx-1,j} = (1-\omega)P_{nx-1,j} + \frac{\omega}{1/dx^{2} + 1/dy^{2}} (\frac{1}{dx^{2}}P_{nx-2,j} + \frac{1}{dy^{2}} \left(P_{nx-1,j+1} + P_{nx-1,j-1} \right)$ $- R_{nx-1,j})$ I apply this condition 4 boundaries and 4 corners. I've gone through my code again and I don't see any issues. Maybe another set of eyes might help Last edited by moshe; January 25, 2018 at 01:43.*

January 25, 2018, 10:30	Ppe bc	#32
moshe Member Moshe De Leon Join Date: Nov 2017 Location: Portugal Posts: 31 Rep Power: 8	To add some more details, here is how I do it in full. First off $uf_{i,j},vf_{i,j}$ represent the face velocities in which they are defined as $uf_{i,j} = 0.5(u_{i+1,j} + u_{i,j}), vf_{i,j} = 0.5(v_{i,j+1}+v_{i,j})$ The way I define the RHS of the PPE is $R_{i,j} \frac{uf_{i,j}^{} - uf^{}_{i-1,j}}{dx} + \frac{vf^{}_{i,j} - vf^{}_{i,j-1}}{dy}$ My Pressure BC is x (which is defined in a similar fashion for y) $\frac{P_{i+1,j} - P_{i,j}}{dx} = \frac{uf^{n+1}_{i,j} - uf^{*}_{i,j}}{dt}$ In short I modify the RHS to take into account u^{n+1} and then remove $\frac{\partial P}{\partial x}$ from $\frac{\partial^{2} P}{\partial x^{2}}$ In terms of your question on why the y direction is not centered, it is because vf is the face so doing the algebra we have $\frac{vf_{i,j} - vf_{i,j-1}}{dy} = \frac{v_{i,j+1} - v_{i,j-1}}{2dy}$

January 26, 2018, 10:36	Some findings and progress	#36
moshe Member Moshe De Leon Join Date: Nov 2017 Location: Portugal Posts: 31 Rep Power: 8	So I spent some time yesterday rewriting my solver with a staggered grid and found some interesting things! First off, I was able to get the incremental projection method to work on the staggered grid! . For my collocated based solver, I was using central interpolation to form the gradient of p^{n} in the intermediate step. In the staggered solver I used one sided interpolation, i.e. $\nabla_{x} q^{n} \approx \frac{q_{E} - q_{P}}{\Delta x}$ $\nabla_{y} q^{n} \approx \frac{q_{N} - q_{P}}{\Delta y}$ When I used a central approximation I would get severe oscillations (numerical instabilities), but this seemed to fix the problem, at least for the staggered grid. I know using a one sided approximation is invalid for the collocated formulation as u,v,w,P are at the same location. In literature, it is not mentioned what type of approximation they use for the q^{n} pressure time level, but it would make sense why a central approximation gave bad results for the staggered formulation......

January 24, 2018, 12:46		#22
FMDenaro Senior Member Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,769 Rep Power: 71	Are you still working on the lid-driven cavity? Are you using a compact stencil for the pressure? Could you show the pressure field? Add a further plot, the (div v) field obtained after the first time step.

January 24, 2018, 13:10		#24
FMDenaro Senior Member Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,769 Rep Power: 71	The pressure plot seems quite regular apart some points...what about the grid? However, the div v plot reaches peak values of almost O(10^-1) that is larger than the magnitude of the local truncation error, what about the time-step value? I suggest to do a first test excluding totally the convective term (very small Re number), that will produce a symmetric field so we can better test the pressure

January 24, 2018, 14:00		#26
FMDenaro Senior Member Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,769 Rep Power: 71	I can suspect your dt is high and drives to a numerical instability. Try to use dt=10^5. Then, what about the level of convergence of the pressure equation? are you using a low threshold on the residual?

January 24, 2018, 18:14		#28
FMDenaro Senior Member Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,769 Rep Power: 71	Do oscillations appear just after the first time step? What about if you refine the mesh? I would think some bugs in the interpolation or in the bc.s for the pressure. Could you write your pressure equation in a node close to a wall or a corner?

January 25, 2018, 04:43		#30
FMDenaro Senior Member Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,769 Rep Power: 71	I am not sure about your equation and the fact that the pressure is colocated at the grid nodes. I assume you have i=1,j=1 and i=N,j=M as walls index. Now for example I write the x-component in i=2 using a compact stencil (dp/dx\|i+1/2 - dp/dx\|i-1/2)/dx = (rho/dt)(u\|i+1/2-u\|i-1/2)/dx while the Neumann BC at i=1 is dp/dx\|1=(rho/dt)(u_wall - u*\|1) W............X.............X -> i direction 1.............2.............3 Is that your pressure colocation? If yes, how do you use the BC into dp/dx\|i-1/2?

January 25, 2018, 05:45		#31
FMDenaro Senior Member Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,769 Rep Power: 71	And why the y derivative of vfs is first order non-centered?

January 25, 2018, 11:47		#34
FMDenaro Senior Member Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,769 Rep Power: 71	I tried to figure it out how do you use the BC after you have discretized but I see problems...for example, what is the goal to express the ghost point for the pressure when you have a compact stencil that does not involve it at all?? The pressure equation is written from i,j=2,2 to nx-1,ny-1, right? Therefore the pressure values involved are at i=1 and nx and j=1 and ny. I don't think that the problem is in the coding but in the formulation you are adopting close to the walls. I will try to work on and write some extended lines that I will post later.

January 26, 2018, 10:52	Incremental PM	#38
moshe Member Moshe De Leon Join Date: Nov 2017 Location: Portugal Posts: 31 Rep Power: 8	Sorry, I was a bit sloppy in my variables let me clarify. So in terms of the pressure time levels we have P^{n+1}, P^{n}, and q^{n+1}. In the intermediate step we have P^{n}, we solve the PPE for q^{n+1} and we do the pressure update with P^{n+1}. Again, sorry for the confusion. In terms of central vs one-sided approximations, I think the pressure in the intermediate step has to be consistent with the corrector-step, though I do not have a rigorous justification yet.

January 26, 2018, 11:05		#39
FMDenaro Senior Member Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,769 Rep Power: 71	Considering the node (i,j) for the q field (colocated at the same grid locations of p^n+1 and p^n), you write d/dx(dq/dx)\|i,j = (dq/dx\|i+1/2,j - dq/dx\|i-1/2,j)/dx d/dy(dq/dy)\|i,j = (dq/dy\|i,j+1/2 - dq/dy\|i,j-1/2)/dy Therefore I don't se any one-sided discretization....You see the effect of the staggering when computing the divergence of the v* field moshe likes this.

January 26, 2018, 11:11	Incremental PM	#40
moshe Member Moshe De Leon Join Date: Nov 2017 Location: Portugal Posts: 31 Rep Power: 8	I agree with what you say. Though, I have found a paper that does an incremental PM method (with a staggered grid) and they( use one sided approximations for P^{n} in the intermediate velocity field. For a collocated grid, the one-sided does not work (I tried it for the sake of curiosity), and it would make sense. * Calculating three-dimensional flows around structures and over rough terrain, C.W. Hirt, J.L Cook, JCP, 1972

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
On the alpha Eqn of VOF method when using Immersed boundary method in OpenFOAM	keepfit	OpenFOAM	4	January 31, 2014 14:32
Attempt to implement the Chorin Projection Method	McFly	OpenFOAM Programming & Development	0	October 26, 2013 06:46
Turbulence inflow generation - recycling method	panda60	OpenFOAM Running, Solving & CFD	15	April 25, 2013 01:34
Fluent 6.3.26 vs 12.1 and partition method	Anorky	FLUENT	0	April 27, 2010 10:55
Comparison: Finite Volume Method vs. Analytic Method	m-fry	Main CFD Forum	1	April 20, 2010 14:40