CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > General Forums > Main CFD Forum

Convergence Issues with Poisson equation in Fractional Step Code

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   July 13, 2016, 06:58
Default Convergence Issues with Poisson equation in Fractional Step Code
  #1
New Member
 
Jack
Join Date: Jul 2016
Posts: 8
Rep Power: 6
FluidFox is on a distinguished road
Hi There,

I have a fractional step code that was written following the description given by Hirsch in 'Numerical Computation of Internal & External Flows', but modified for an axisymmetric cylindrical-polar case, and am experiencing problems with the Poisson equation for the pseudo-pressure that is used to correct the velocity to satisfy continuity.

I am using the boundary conditions dp/dn = 0 normal derivative to the wall being zero as discussed in Peyret and Taylor and by Romper and others. The method I use to solve is either TDMA algorithm or Gauss-Seidel until a convergence limit is reached based on the maximum change between iterations.

When I isolate the Poisson equation solver and run for an idealised test case all is fine, and the solution converges to an error of around 10^-16. However, when I run it as part of the overall code solving for pressure with the divergence of u_star as the source term, the best it can manage is 10^-6, and it never converges, instead the change between iterations reaches a constant value such that the solution everywhere keeps growing by this constant each iteration.
If I use the solution obtained at the 10^-6 limit to correct the velocity, the maximum divergence drops from around 0.5 to 0.05, so the method is working, but this still seems quite large to say that mathematically incompressibility should be satisfied by this method...

I know that mathematically for zero gradient boundary conditions all the way around the boundary to be satisfied, the integral of the source function across the domain must be zero, and that is not the case for my divergence of u_star.


So my questions are:
Should this be the case, or should the divergence of the intermediate velocity across the domain integrate to zero to satisfy the boundary conditions and therefore might there be something wrong elsewhere in my code?

I know there is some discussion of how physical these boundary conditions are... Is there a better option that would give me improved convergence (e.g. non-zero gradient)?

Am I using the best algorithm for very simple, unsteady axisymmetric flow (rectangular mesh, one inlet, one outlet) or is there a better option? I did some research and this seemed the best pressure correction scheme, though I believe transient SIMPLE would also have been an option... If anyone has any suggestions and a source where I can learn in detail how the algorithm works I would be grateful.


Apologies if this is long, I have not used this service directly before.
Many thanks.
FluidFox is offline   Reply With Quote

Old   July 13, 2016, 07:07
Default
  #2
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 5,503
Rep Power: 58
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
The main comment is that the surface integral of the source term MUST

Int[S] n.v* dS =0

this is a mandatory condition that stems out from the decomposition

v* = vn+1 + Grad phi

In conclusion you have to check the BC.s you are using in computing Divv*.

Is the grid arrangement staggered? If you are working with colocated method are you working with the Exact or Approximate Projection Method? Be careful that the APM will lead to a divergence-free constraint in each cell satisfied only up to the local truncation error.
FMDenaro is online now   Reply With Quote

Old   July 13, 2016, 07:55
Default
  #3
New Member
 
Jack
Join Date: Jul 2016
Posts: 8
Rep Power: 6
FluidFox is on a distinguished road
Hi, thank you for the quick response.
The integral around the boundary is zero (enforced by boundary conditions), but it is the integral over the surface of the domain (via Gauss' Theorem) that does not appear to be zero (I am currently investigating this).

I am using a staggered grid with velocities calculated at the edges of cells and pressure in the centre (The cell-centered finite volume layout using a staggered grid approach as described in Hirsch section 12.4.3). The boundary conditions are imposed by the use of a layer of 'ghost' cells outside the domain.

I am using the method by Chorin described in Hirsch, Peyret and Taylor, and Biringen. All the code is written to double precision, so although I know there will be some numerical error, the values I am seeing at present seem quite large. I will check through the divergence and boundary calculations carefully again to see if there are issues here.
FluidFox is offline   Reply With Quote

Old   July 13, 2016, 08:05
Default
  #4
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 5,503
Rep Power: 58
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
Quote:
Originally Posted by FluidFox View Post
Hi, thank you for the quick response.
The integral around the boundary is zero (enforced by boundary conditions), but it is the integral over the surface of the domain (via Gauss' Theorem) that does not appear to be zero (I am currently investigating this).

I am using a staggered grid with velocities calculated at the edges of cells and pressure in the centre (The cell-centered finite volume layout using a staggered grid approach as described in Hirsch section 12.4.3). The boundary conditions are imposed by the use of a layer of 'ghost' cells outside the domain.

I am using the method by Chorin described in Hirsch, Peyret and Taylor, and Biringen. All the code is written to double precision, so although I know there will be some numerical error, the values I am seeing at present seem quite large. I will check through the divergence and boundary calculations carefully again to see if there are issues here.
no, you do not have to use ghost points...and you do not have to write the Laplacian operator but the Div Grad one. This way, on the boundary you directly substitute the condition n.Grad phi =0 into the matrix without any further discretization.
Have you set n.v*=0 on the boundary?
FMDenaro is online now   Reply With Quote

Old   July 13, 2016, 08:22
Default
  #5
New Member
 
Jack
Join Date: Jul 2016
Posts: 8
Rep Power: 6
FluidFox is on a distinguished road
The no slip-condition (including zero flux through) is imposed at all walls for both the actual velocities as well as the predictor velocities u*. You are correct that I do not actually use the ghost cells in the TDMA or Gauss-Seidel algorithms (instead modifying the ap coefficient for the edge cells).

How do you mean replacer Laplacian with Div(Grad)? I thought that was the definition of the laplacian... sorry for not understanding.
FluidFox is offline   Reply With Quote

Old   July 13, 2016, 09:01
Default
  #6
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 5,503
Rep Power: 58
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
Since the Grad phi enters into the momentum equation, you apply the Div operator to all terms and enforces the condition Div vn+1=0.
That means you get the Poisson equation

Div Grad phi = Div v*

that must be discretized using the same discrete DIv and Grad operation you adopt in the continuity and momentum. On the boundary you apply dphi/dn directly in it, for example in x direction:

[(d phi/dx)|i+1/2 - (d phi/dx)|i-1/2]/dx = u*|i+1/2-u*|i-1/2

if i-1/2 is for example a wall you get directly

[(d phi/dx)|i+1/2 ]/dx = u*|i+1/2

with

(d phi/dx)|i+1/2= (phi(i+1)-phi(i))/dx
FMDenaro is online now   Reply With Quote

Reply

Tags
convergence issues, fractional step method, poisson equation, poisson solver

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Help for the small implementation in turbulence model shipman OpenFOAM Programming & Development 25 March 19, 2014 11:08
Micro Scale Pore, icoFoam gooya_kabir OpenFOAM Running, Solving & CFD 2 November 2, 2013 14:58
User code, How can I get information from last time step onlyacan STAR-CCM+ 1 July 24, 2012 14:52
TVD RK and fractional step method HaKu Main CFD Forum 0 October 19, 2009 03:56
Design Integration with CFD? John C. Chien Main CFD Forum 19 May 17, 2001 16:56


All times are GMT -4. The time now is 14:15.