Incremental projection method

moshe · January 15, 2018, 12:18

I am curious about the practical implementation of an incremental projection method. For pressure we have 3 time levels, $P^{n+1}, \phi^{n+1}, P^{n}$ . In the prediction step we have

$\frac{\mathbf{u}^{*} - u^{n}}{\Delta t} = - \mathbf{u}^{n} \cdot \nabla \mathbf{u}^{n} - \frac{1}{\rho} \nabla P^{n} - \nu \nabla^{2} \mathbf{u}^{n}$

We solve for pressure $\phi^{n+1}$

$\nabla^{2} \phi^{n+1} = \frac{\rho}{\Delta t} \nabla \cdot \mathbf{u}^{*}$

We then project onto

$\frac{\mathbf{u}^{n+1} - \mathbf{u}^{n}}{\Delta t} = - \nabla \phi^{n+1}$

Finally we update pressure

$P^{n+1} = P^{n} + \phi^{n+1}$

In literature $P^{n}$ is referred to as an approximate to $P^{n+1}$ . My question is in practice how do I initialize $P^{n}$ and then update pressure. In my code I have, after updating the velocities

Code:

Pnew(i,j) = P(i,j) + phi(i,j)
phi(i,j) = P(i,j)

I have tried setting P = Pnew, but my simulation blows up.

FMDenaro · January 15, 2018, 12:27

In your case, the discrete momentum equation is explicit and you have no particular problem in setting an initial condition for Pn. Be careful to the BC.s.
Have a look to this paper http://129.81.170.14/~cortez/Prints/accproj.pdf

moshe · January 15, 2018, 12:38

I have recently read this very paper. For pnew, phi, and p, I am setting them to be initially 0. When I set pnew = p + phi I am not applying anything boundary condition to pnew or p. I only apply a BC for phi for the pressure poisson equation. Given p^n is the previous time level of p^n+1 shouldnt I set p^n = p^n+1 for the pressure update?

In terms of the paper I'm following I am studying this paper: http://www.math.tamu.edu/~guermond/P...o_jcp_2009.pdf

FMDenaro · January 15, 2018, 12:43

Quote:

Originally Posted by moshe

I have recently read this very paper. For pnew, phi, and p, I am setting them to be initially 0. When I set pnew = p + phi I am not applying anything boundary condition to pnew or p. I only apply a BC for phi for the pressure poisson equation. Given p^n is the previous time level of p^n+1 shouldnt I set p^n = p^n+1 for the pressure update?

In terms of the paper I'm following I am studying this paper: http://www.math.tamu.edu/~guermond/P...o_jcp_2009.pdf

What do you mean for "I am not applying anything boundary condition to pnew or p"? Depending on the type of collocation, you might need to prescribe the values at the boundaries to compute the correct gradients to correct the velocity.
I suggest to do simple test and running the code in debugging only for the first time step. Check that the divergence-free constraint is ensured at the end of the time step cell-by-cell.

moshe · January 15, 2018, 13:05

What I meant was I only apply $\nabla \phi \cdot n = \frac{1}{\Delta t} (\mathbf{u}^{*} - \mathbf{u}^{n+1})$ to $\phi$ . Not Pnew or Pn. I am trying a collocated arrangement and I am not getting a divergence free velocity field (up to machine precision.) I believe what I am getting is results to that of an approximate projection method. My velocity profile looks correct, but my pressure profile does not seem to be right.

FMDenaro · January 15, 2018, 13:21

Quote:

Originally Posted by moshe

What I meant was I only apply $\nabla \phi \cdot n = \frac{1}{\Delta t} (\mathbf{u}^{*} - \mathbf{u}^{n+1})$ to $\phi$ . Not Pnew or Pn. I am trying a collocated arrangement and I am not getting a divergence free velocity field (up to machine precision.) I believe what I am getting is results to that of an approximate projection method. My velocity profile looks correct, but my pressure profile does not seem to be right.

Yes, using colocated arrangement and a compact stencil for the pressure leads to the approximate projection, therefore you should have the error in the divergence-free constraint of the magnitude order of the local truncation error. However, when you correct the velocities close to the boundaries you need to prescribe a condition for Pnew to compute correctly the pressure gradients.

FMDenaro · January 15, 2018, 13:25

Finally, consider that in the non-staggered grids some wiggles can be produced, for example see https://www.researchgate.net/publica...y-driven_flows

moshe · January 15, 2018, 13:41

Yes, this is one of the frustrations I am finding in using a collocated arrangement and a compact stencil for pressure. May I ask you what I should prescribe to Pnew at the boundaries to get the correct pressure gradient profile?

FMDenaro · January 15, 2018, 13:48

When you write the increment [LaTeX Error: Syntax error][/COLOR]
at the boundary you should consider that the normal derivative on phi was prescribed to be zero therefore you can use the relation d Pn+1/dn = d Pn/dn.
However, I think that this is the reason why you get a numerical boundary layer in the pressure field as reported by several papers.

Unfortunately I use the pressure-free projection, not the incremental, therefore withour a direct experience I should better think about the use of the Hodge decomposition to avoid such problem.

moshe · January 15, 2018, 14:12

I appreciate your help professor Denaro! I was curious: For DNS or LES do people still use staggered grids in light of the problems of collocated grids?

FMDenaro · January 15, 2018, 14:18

Quote:

Originally Posted by moshe

I appreciate your help professor Denaro! I was curious: For DNS or LES do people still use staggered grids in light of the problems of collocated grids?

Hystorically the staggered grid with second order central discretization was the standard in DNS/LES but now the non-staggered grid is largely used.

moshe · January 16, 2018, 11:28

I have been doing some readings on your papers with regards to pressure-free projection methods. It has led me to reading some papers by Kim and Moin. With regards to their statement of the "actual pressure"

$P^{n+1} = \phi^{n+1} - \frac{\Delta t \nu}{2} \nabla^{2} \phi^{n+1}$

do I apply $\nabla P^{n+1} \cdot n = 0$ to P^{n+1} on the boundary or can I leave this term without any boundary treatment. Their paper is not very clear. Thank you!

FMDenaro · January 16, 2018, 11:35

Quote:

Originally Posted by moshe

I have been doing some readings on your papers with regards to pressure-free projection methods. It has led me to reading some papers by Kim and Moin. With regards to their statement of the "actual pressure"

$P^{n+1} = \phi^{n+1} - \frac{\Delta t \nu}{2} \nabla^{2} \phi^{n+1}$

do I apply $\nabla P^{n+1} \cdot n = 0$ to P^{n+1} on the boundary or can I leave this term without any boundary treatment. Their paper is not very clear. Thank you!

You have to be careful of the fact that such relation appears only when you use the implicit Crank-Nicolson time integration. In your method this relation does not apply.
However, rigorously speaking, the Neuman BC.s are just the projection of the Hodge decomposition along the normal direction. This is valid for any type of projection method.

moshe · January 16, 2018, 11:55

I agree, in my case P^{n+1} = \phi^{n+1} as I am using FE integration. I quickly changed my solver to a pressure-free method and I update my pressure as follows

Code:

for j=1:ny
for i=1:nx
Pnew(i,j) = phi(i,j)
end
end

If I were to use AB2-CN I would have

Code:

for j=2:ny-1
for i=2:nx-1
Pnew(i,j) = phi(i,j) - (dt*nu)/2.0*Laplace(P)
end
end

Since Pnew recovers the physical pressure, I do not need to impose any BC to Pnew?

FMDenaro · January 16, 2018, 12:46

Quote:

Originally Posted by moshe

I agree, in my case P^{n+1} = \phi^{n+1} as I am using FE integration. I quickly changed my solver to a pressure-free method and I update my pressure as follows

Code:

for j=1:ny
for i=1:nx
Pnew(i,j) = phi(i,j)
end
end

If I were to use AB2-CN I would have

Code:

for j=2:ny-1
for i=2:nx-1
Pnew(i,j) = phi(i,j) - (dt*nu)/2.0*Laplace(P)
end
end

Since Pnew recovers the physical pressure, I do not need to impose any BC to Pnew?

But you don’t need to recover the pressure, the method uses only phi. And still a Neumann condition is required.

FMDenaro · January 16, 2018, 15:15

I forgot to mention that the method in the paper of Kim and Moin uses a staggered grid

moshe · January 17, 2018, 11:58

So after reading the paper by Brown. I noticed that $P^{n+1} = \phi^{n+1}$ gives reduced accuracy as opposed to update pressure at the

time level. As such I have $P^{n+1/2} = \phi^{n+1}$ (i.e. I pressure update before I perform the corrector step) As such, my pressure profile is shown below. Since the paper by Moin used a staggered arrangement, I am under the impression I will get slightly less accurate results than that of the staggered mesh.

I have come across that there exists an improvement of the Rhie-Chow interpolation, titled "Improved Rhie-Chow Interpolation for Unsteady Flow Computations." Unfortunately, I do not have access to this paper, would this possibly increase the accuracy of my code? In the original paper by Rhie and Chow, it is used for a steady state turbulence problem.

FMDenaro · January 17, 2018, 12:17

Your questions implies several issues... And for each one the answer is not simple to explain here.

1) Despite of what stated in the paper of Brown et al., the accuracy of the pressure-free projection method (be careful that here one refers as only to the time accuracy of the splitting method) can be high as you want. You need to consider the proper BC.s for the intermediate velocity. Then pressure can be recovered from phi.

2) The staggered grid in the paper of Kim and Moin is associated to a classical second order discretization that allows to get an exact projection method. Non-staggered grids with second order discretization have the problem of non-communicating grids for the pressure equation (large stencil) if you perform an exact projection method. Therefore, generally the approximate projection method is used on a compact stencil. However, the spatial accuracy is formally the same.

3) The R-C interpolation does not increase the order of accuracy. It is used to get a strong coupling. You can see also the discussion in the book of Peric and Ferziger.

However, I just want to alert you abouyt the fact that the deeper and often hidden details of the projections methods have taken me engaged for some years...

moshe · January 18, 2018, 11:06

Thank you for letting me pick your brain

. I did not mean that the Rhie-Chow interpolation increases the accuracy of the method. What I was trying to get at is the improved Rhie-Chow method will provide a more stable solution for unsteady flow computations.

I just want to pick your brain on one other thing. So the BCs for the intermediate velocity step is a source of ambiguity for pressure-free PMs. To get second-order accuracy Moin proposed $u^{*} \cdot n = u^{n+1} \cdot n+ \Delta t \nabla \phi^{n} \cdot n$ . In my current solver I have implemented CN-AB2. I only have one phi variable and that is the n+1 time level. Is that sufficient to just use n+1 in place of the n level for the BC? In a practical implementation I'm not sure how one would have $\phi^{n}$ using an SOR type iterator as it only uses one variable.

FMDenaro · January 18, 2018, 11:46

Quote:

Originally Posted by moshe

Thank you for letting me pick your brain

. I did not mean that the Rhie-Chow interpolation increases the accuracy of the method. What I was trying to get at is the improved Rhie-Chow method will provide a more stable solution for unsteady flow computations.

I just want to pick your brain on one other thing. So the BCs for the intermediate velocity step is a source of ambiguity for pressure-free PMs. To get second-order accuracy Moin proposed $u^{*} \cdot n = u^{n+1} \cdot n+ \Delta t \nabla \phi^{n} \cdot n$ . In my current solver I have implemented CN-AB2. I only have one phi variable and that is the n+1 time level. Is that sufficient to just use n+1 in place of the n level for the BC? In a practical implementation I'm not sure how one would have $\phi^{n}$ using an SOR type iterator as it only uses one variable.

Well, first the problem of the BC.s for the intermediate velocity exists only when using an implicit time integration. Then, the formulation proposed by Kim and Moin is actually only first order accurate therefore decreases the accuracy order close to a boundary.
Now, the evaluation of the pressure field (or phi) in time has to be carefully considered. Consider for example the momentum equation (convection and diffusion is denoted by R) dv/dt + grad p = R, let integrate exactly in time
vn+1-vn + Int [tn,tn+1] grad p dt = Int [tn,tn+1] R dt (1)

Now we introduce an implicit discretization in R and disregard the pressure term

v*-vn = Int_dis [tn,tn+1] R* dt (2)

and, after the Poisson equation allows to solve for the field phi we can write

vn+1=v* - dt*grad phi (3)

If you substitute in this relation the equation (2) for the intermediate field

vn+1-vn + dt*grad phi = Int_dis [tn,tn+1] R* dt (4)

Now, comparing (1) and (2), you see that grad phi is actually correlated to a the time averaging in the interval tn,tn+1, is not a function evaluated at a certain time.

If you are interested in the details of the formulation I analysed some years ago, you will find useful these papers:

https://www.researchgate.net/publica...ary_conditions
https://www.researchgate.net/publica...ection_methods
https://www.researchgate.net/publica...blicationTitle

January 15, 2018, 12:18	Incremental projection method	#1
moshe Member Moshe De Leon Join Date: Nov 2017 Location: Portugal Posts: 31 Rep Power: 8	I am curious about the practical implementation of an incremental projection method. For pressure we have 3 time levels, $P^{n+1}, \phi^{n+1}, P^{n}$ . In the prediction step we have $\frac{\mathbf{u}^{} - u^{n}}{\Delta t} = - \mathbf{u}^{n} \cdot \nabla \mathbf{u}^{n} - \frac{1}{\rho} \nabla P^{n} - \nu \nabla^{2} \mathbf{u}^{n}$ We solve for pressure $\phi^{n+1}$ $\nabla^{2} \phi^{n+1} = \frac{\rho}{\Delta t} \nabla \cdot \mathbf{u}^{}$ We then project onto $\frac{\mathbf{u}^{n+1} - \mathbf{u}^{n}}{\Delta t} = - \nabla \phi^{n+1}$ Finally we update pressure $P^{n+1} = P^{n} + \phi^{n+1}$ In literature $P^{n}$ is referred to as an approximate to $P^{n+1}$ . My question is in practice how do I initialize $P^{n}$ and then update pressure. In my code I have, after updating the velocities Code: Pnew(i,j) = P(i,j) + phi(i,j) phi(i,j) = P(i,j) I have tried setting P = Pnew, but my simulation blows up.

January 15, 2018, 13:05	Projection Method	#5
moshe Member Moshe De Leon Join Date: Nov 2017 Location: Portugal Posts: 31 Rep Power: 8	What I meant was I only apply $\nabla \phi \cdot n = \frac{1}{\Delta t} (\mathbf{u}^{*} - \mathbf{u}^{n+1})$ to $\phi$ . Not Pnew or Pn. I am trying a collocated arrangement and I am not getting a divergence free velocity field (up to machine precision.) I believe what I am getting is results to that of an approximate projection method. My velocity profile looks correct, but my pressure profile does not seem to be right.

January 15, 2018, 13:41	Boundary condition for Pnew	#8
moshe Member Moshe De Leon Join Date: Nov 2017 Location: Portugal Posts: 31 Rep Power: 8	Yes, this is one of the frustrations I am finding in using a collocated arrangement and a compact stencil for pressure. May I ask you what I should prescribe to Pnew at the boundaries to get the correct pressure gradient profile?

January 15, 2018, 13:48		#9
FMDenaro Senior Member Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,777 Rep Power: 71	When you write the increment [LaTeX Error: Syntax error][/COLOR] at the boundary you should consider that the normal derivative on phi was prescribed to be zero therefore you can use the relation d Pn+1/dn = d Pn/dn. However, I think that this is the reason why you get a numerical boundary layer in the pressure field as reported by several papers. Unfortunately I use the pressure-free projection, not the incremental, therefore withour a direct experience I should better think about the use of the Hodge decomposition to avoid such problem. moshe likes this.

January 15, 2018, 14:12	Projection Method	#10
moshe Member Moshe De Leon Join Date: Nov 2017 Location: Portugal Posts: 31 Rep Power: 8	I appreciate your help professor Denaro! I was curious: For DNS or LES do people still use staggered grids in light of the problems of collocated grids?

January 15, 2018, 12:27		#2
FMDenaro Senior Member Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,777 Rep Power: 71	In your case, the discrete momentum equation is explicit and you have no particular problem in setting an initial condition for Pn. Be careful to the BC.s. Have a look to this paper http://129.81.170.14/~cortez/Prints/accproj.pdf

January 15, 2018, 13:25		#7
FMDenaro Senior Member Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,777 Rep Power: 71	Finally, consider that in the non-staggered grids some wiggles can be produced, for example see https://www.researchgate.net/publica...y-driven_flows

January 16, 2018, 11:28	Pressure free projection method vs incremental projection method	#12
moshe Member Moshe De Leon Join Date: Nov 2017 Location: Portugal Posts: 31 Rep Power: 8	I have been doing some readings on your papers with regards to pressure-free projection methods. It has led me to reading some papers by Kim and Moin. With regards to their statement of the "actual pressure" $P^{n+1} = \phi^{n+1} - \frac{\Delta t \nu}{2} \nabla^{2} \phi^{n+1}$ do I apply $\nabla P^{n+1} \cdot n = 0$ to P^{n+1} on the boundary or can I leave this term without any boundary treatment. Their paper is not very clear. Thank you!

January 16, 2018, 11:55	Pressure-free PM	#14
moshe Member Moshe De Leon Join Date: Nov 2017 Location: Portugal Posts: 31 Rep Power: 8	I agree, in my case P^{n+1} = \phi^{n+1} as I am using FE integration. I quickly changed my solver to a pressure-free method and I update my pressure as follows Code: for j=1:ny for i=1:nx Pnew(i,j) = phi(i,j) end end If I were to use AB2-CN I would have Code: for j=2:ny-1 for i=2:nx-1 Pnew(i,j) = phi(i,j) - (dtnu)/2.0Laplace(P) end end Since Pnew recovers the physical pressure, I do not need to impose any BC to Pnew?

January 16, 2018, 15:15		#16
FMDenaro Senior Member Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,777 Rep Power: 71	I forgot to mention that the method in the paper of Kim and Moin uses a staggered grid moshe likes this.

January 17, 2018, 12:17		#18
FMDenaro Senior Member Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,777 Rep Power: 71	Your questions implies several issues... And for each one the answer is not simple to explain here. 1) Despite of what stated in the paper of Brown et al., the accuracy of the pressure-free projection method (be careful that here one refers as only to the time accuracy of the splitting method) can be high as you want. You need to consider the proper BC.s for the intermediate velocity. Then pressure can be recovered from phi. 2) The staggered grid in the paper of Kim and Moin is associated to a classical second order discretization that allows to get an exact projection method. Non-staggered grids with second order discretization have the problem of non-communicating grids for the pressure equation (large stencil) if you perform an exact projection method. Therefore, generally the approximate projection method is used on a compact stencil. However, the spatial accuracy is formally the same. 3) The R-C interpolation does not increase the order of accuracy. It is used to get a strong coupling. You can see also the discussion in the book of Peric and Ferziger. However, I just want to alert you abouyt the fact that the deeper and often hidden details of the projections methods have taken me engaged for some years...

January 18, 2018, 11:06	Projection Methods	#19
moshe Member Moshe De Leon Join Date: Nov 2017 Location: Portugal Posts: 31 Rep Power: 8	Thank you for letting me pick your brain . I did not mean that the Rhie-Chow interpolation increases the accuracy of the method. What I was trying to get at is the improved Rhie-Chow method will provide a more stable solution for unsteady flow computations. I just want to pick your brain on one other thing. So the BCs for the intermediate velocity step is a source of ambiguity for pressure-free PMs. To get second-order accuracy Moin proposed $u^{*} \cdot n = u^{n+1} \cdot n+ \Delta t \nabla \phi^{n} \cdot n$ . In my current solver I have implemented CN-AB2. I only have one phi variable and that is the n+1 time level. Is that sufficient to just use n+1 in place of the n level for the BC? In a practical implementation I'm not sure how one would have $\phi^{n}$ using an SOR type iterator as it only uses one variable.

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