[FluidFox]

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.

[FMDenaro]

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.

[FluidFox]

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.

[FMDenaro]

Originally Posted by FluidFox

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?

[FluidFox]

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.

[FMDenaro]

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