[moomba]

Hi

During my training period, I have to develop a 2D DNS solver for incompressible flow (dP/dx=0 and dRho/dx=0). The system is supposed to be air open, so P=P_atm and dP/dt=0.


I am trying to compute a simple vortex in a convecting flow with periodical boundaries in X and Y. The program works well for a specified duration (the vortex convect and disappear from right to reappear on the left), but after a fixed number of computation time, it explode.
I have found that even with the poisson correction for momentum equation, the term Div(Ui) is not null. (which is not correct according to the conservation law for incompressible flow : Div(Ui)=0 )
This term is increasing, and when it is sufficient to weight, it make the system explode.


I think it is due to the fact that my poisson solver is second order accurate (FishPack90) when my discretization schema is fourth order accurate (FD4). I first tried with a 6th order PADE, but it was too much instable.


Does anyone know a poisson solver that is made for fourth order finite difference ?


Thanks in advance.



I do apologize for my very bad English

[harishg]

If the system is unstable for second order pressure poisson equation, I doubt it would be stable for a higher order system. How do you setup your system and how do you implement boundary conditions ? Be careful to use a smaller stencil at the domain boundaries while using high order schemes to maintain stability.

[moomba]

Hi

Thank you for your reply.

In fact, my domain is fully periodic, so I don't have implemented the boundaries yet. I wait that it work with periodic conditions to code the wall and the outflows I need. So my derivatives are FD4 on the whole system, even at boundaries.

I have plot the Div(U) value depending on iterations :




As you can see, there is a problem.

I am doing this :

Density (Rho) =1.0d0
Pressure = 1.0d0
Temperature = 1.0d0

Then, I calculate Hi=dTij/dxi - dRhoUiUj/dxi
Then the advanced velocity in time without respect do continuity : RhoUi*=RhoUi^n + Dt*Hi

From here, P means the dynamique pressure.

My poisson equation is then : d/dx1 dP/dx1 + d/dx2 dP/dx2 = 1/Dt * ( dRhoU1*/dx1 + dRhoU2*/dx2 )

To finish, I find the real velocity advanced in time (that respect continuity) :
RhoUi^n+1=RhoUi*-Dt*dP/dxi

But it doesn't respect continuity, and I don't understand why...

I hope these explanations will help you to find the problem.

I thank you in advance.

[harishg]

How do you implement your convection term ? When you perform DNS, with higher order schemes, it is recommended that you use the skew symmetric form of the convection term to avoid aliasing errors. I guess one of the possible reasons for the code blowing up might be the aliasing error.

Do you use staggered or collocated grids ? From the equations it seems that you are trying to satisfy the continuity constraint using fractional step method. How do you implement your boundary conditions for pressure ? As a final thought you can use the spectral channel flow code www.channelflow.org to perform your simulations, if it is not mandatory to develop your own code.

[moomba]

Hi. Sorry to be late.

My convection term is calculated as follow. Q refer to conservative variables (Rho*Ui), and U to simple Variable (Ui). There are NP points, and
Ndim dimensions (here, Ndim=2 because it is a 2D calculation).
So I calculate for the axis 'Ivar' the value RhoUi(Ivar)*U1. It is derivated and store in H(Ivar). Then I calculate RhoUi(Ivar)*U2 and it is added to H(Ivar).

    
do Ivar=1,Ndim    ! axis Ivar

    J=1
    do IP = 1,NP
        WK1(IP) = Q(IP,Ivar)*U(IP,J)    ! WK1 = RhoUi*UJ = RhoUi*U1
    enddo
    CALL DERIVATIVE(J,WK1,H(:,Ivar))        ! H=dRhoUiUj/dxj = dRhoUiU1/dx1

    do J = 2,Ndim
        do IP = 1,NP
            WK1(IP) = Q(IP,Ivar)*U(IP,J)    ! WK1 = RhoUi*U2 et 3
        enddo
        CALL DERIVATIVE(J,WK1,WK2)    ! H=dRhoUiUj/dxj 
        do IP = 1,NP
            H(IP,Ivar) = H(IP,Ivar)+WK2(IP)    ! H=Som(dRhoUiUj/dxj )
        enddo
    enddo
end do

Concerning your recommandation :

it is recommended that you use the skew symmetric form of the convection term to avoid aliasing errors

I don't know this way to calculate convection term. Do you have some reference where I could read about it ? I found many things on internet, and I am not sure about it.

Do you use staggered or collocated grids ?

The code is Collocated for now, and will be staggered when finished. I thought it would be easier to construct it collocated and then transform it staggered.

How do you implement your boundary conditions for pressure ?

In fact, there are no boundary conditions, the flow is fully periodic for the moment. I was thinking that it could be the source of the problem, because there is nowhere for the code to see the 'outside'.

Presently, I am coding a FD4 Poisson solver using Conjugate Gradient algorithm, to make it fully compatible with the FD4 of the code. (If I can't find the problem before the end of my training period, I will be then able to prove I have done everything I could to make it work. )


Thank you for your software recomandation. In fact, constructing the code is part of the training perdioc, so I don't have the possibility to use an other one. But with this software, I will have something to make a comparison with my resulys.


I thank you in advance

[andy_]

I think you need to read a bit about incompressible solvers. A few points to get you going:
- a staggered Cartesian grid will be easier to sort out than a collocated one because you will not need to consider how to handle pressure smoothing.
- the pressure equation with the usual Neumann boundary conditions is singular and therefore needs to be solved with care. This largely boils down to ensuring the source terms always sum exactly to zero and preventing the average value from wandering too far from zero.
- when you perform your continuity correction step it must also account for any residual continuity error left from the previous step. If you assume the continuity error from the previous step is exactly zero (perhaps the obvious thing to do when deriving a solution procedure) it tends to accumulate with time which might be what you are seeing.
- non-periodic boundary conditions for time varying two step schemes which advance the "velocity" field using some approximate pressure field and then correct to impose continuity require careful consideration if they are to maintain formal accuracy and not generate strong gradients in the intermediate variables next to boundaries.

[moomba]

Thank you for this answer.

- when you perform your continuity correction step it must also account for any residual continuity error left from the previous step. If you assume the continuity error from the previous step is exactly zero (perhaps the obvious thing to do when deriving a solution procedure) it tends to accumulate with time which might be what you are seeing.

I thought about it, but I can't figure out how to implement it. I follow the alogorithme I described earlier :

Density (Rho) =1.0d0
Pressure = 1.0d0
Temperature = 1.0d0

Then, I calculate Hi=dTij/dxi - dRhoUiUj/dxi
Then the advanced velocity in time without respect do continuity : RhoUi*=RhoUi^n + Dt*Hi

From here, P means the dynamique pressure.

My poisson equation is then : d/dx1 dP/dx1 + d/dx2 dP/dx2 = 1/Dt * ( dRhoU1*/dx1 + dRhoU2*/dx2 )

To finish, I find the real velocity advanced in time (that respect continuity) :
RhoUi^n+1=RhoUi*-Dt*dP/dxi

In fact, it is not d/dx1 dP/dx1 + d/dx2 dP/dx2 = 1/Dt * ( dRhoU1*/dx1 + dRhoU2*/dx2 ) that should be solved for the poisson equation, but :
d/dx1 dP/dx1 + d/dx2 dP/dx2 = 1/Dt * ( dRhoU1*/dx1 + dRhoU2*/dx2 - dRhoU1^n+1/dx1 - dRhoU2^n+1/dx1)
that is to say :

Delta(P)=1/Dt * (Div(RhoUi*) - Div(RhoUi^n+1))

But Div(RhoUi^n+1) should be null because Rho is constant. Then, where in the calcul could I "account for any residual continuity error left from the previous step" ? It must be obvious, but not for me

[andy_]

After you have solved your pressure correction equation, the pressure gradient between the pressure nodes can be used directly to correct the velocity at the cell faces to satisfy continuity exactly. But you are not doing this. You are averaging the pressure gradients at the cell walls and applying the correction to the cell centre velocity. The continuity error in each individual cell will be large even though the error will sum to zero throughout the region.

The velocity term in the momentum equation is non-linear and you have a local continuity error that your scheme cannot "see" to remove. If you do nothing, over time the gradients in your fields will get steeper, the local continuity error larger and the scheme will eventually blow up.

Staggering the velocity components directly couples the velocity gradient and its pressure gradient removing all these problems. Alternatively, you can use pressure smoothing which is a weaker mechanism but in widespread use.

[moomba]

Sorry to be late
I am full of work today...

Thank you very much for all this help

After you have solved your pressure correction equation, the pressure gradient between the pressure nodes can be used directly to correct the velocity at the cell faces to satisfy continuity exactly. But you are not doing this. You are averaging the pressure gradients at the cell walls and applying the correction to the cell centre velocity. The continuity error in each individual cell will be large even though the error will sum to zero throughout the region.

I don't understand, because I am using a centered derivative (FD4 centered) so the value dP/dxi should be at the same position than the velocity points, am I wrong ? If not, then I will consider modifing my code to make it use staggered grid.

In a staggered grid, there are the following grid :

• Scalar variable grid (contain Rho, Temperature, Pressure, Species Fraction)
• Velocity Grids. One Grid per velocity (Grid U1 and grid U2 for me)
• Derivative Grid. Contain all derivatives values.
• The structure grid, that is to say the reference grid.

Is all of this correct ?

I am sorry to ask so much question. In fact, I am alone to construct this code and I don't have a lot of time to finish it, so it's difficult