CFD Online Discussion Forums

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   Main CFD Forum (https://www.cfd-online.com/Forums/main/)
-   -   Incremental projection method (https://www.cfd-online.com/Forums/main/197744-incremental-projection-method.html)

moshe January 15, 2018 12:18

Incremental projection method
 
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

Practical implementation of projection method
 
1 Attachment(s)
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 (Post 678182)
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

Projection Method
 
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 (Post 678189)
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

Boundary condition for Pnew
 
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

Projection Method
 
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 (Post 678201)
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

Pressure free projection method vs incremental projection method
 
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 (Post 678351)
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

Pressure-free PM
 
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 (Post 678357)
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

Pressure profile
 
1 Attachment(s)
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 n+1/2 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

Projection Methods
 
Thank you for letting me pick your brain :D. 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 (Post 678593)
Thank you for letting me pick your brain :D. 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


All times are GMT -4. The time now is 03:55.