
[Sponsors] 
Convergence Issues with Poisson equation in Fractional Step Code 

LinkBack  Thread Tools  Search this Thread  Display Modes 
July 13, 2016, 06:58 
Convergence Issues with Poisson equation in Fractional Step Code

#1 
New Member
Jack
Join Date: Jul 2016
Posts: 8
Rep Power: 9 
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 cylindricalpolar case, and am experiencing problems with the Poisson equation for the pseudopressure 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 GaussSeidel 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. nonzero 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. 

July 13, 2016, 07:07 

#2 
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,675
Rep Power: 71 
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 divergencefree constraint in each cell satisfied only up to the local truncation error. 

July 13, 2016, 07:55 

#3 
New Member
Jack
Join Date: Jul 2016
Posts: 8
Rep Power: 9 
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 cellcentered 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. 

July 13, 2016, 08:05 

#4  
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,675
Rep Power: 71 
Quote:
Have you set n.v*=0 on the boundary? 

July 13, 2016, 08:22 

#5 
New Member
Jack
Join Date: Jul 2016
Posts: 8
Rep Power: 9 
The no slipcondition (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 GaussSeidel 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. 

July 13, 2016, 09:01 

#6 
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,675
Rep Power: 71 
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)i1/2]/dx = u*i+1/2u*i1/2 if i1/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 

Tags 
convergence issues, fractional step method, poisson equation, poisson solver 
Thread Tools  Search this Thread 
Display Modes  


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  STARCCM+  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 