
[Sponsors] 
May 18, 2009, 11:21 
DNS Incompressible : poisson order and Div(Ui)

#1 
New Member
Join Date: May 2009
Posts: 9
Rep Power: 9 
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 

May 19, 2009, 23:08 

#2 
Senior Member
N/A
Join Date: Mar 2009
Posts: 189
Rep Power: 9 
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.


May 20, 2009, 05:04 

#3 
New Member
Join Date: May 2009
Posts: 9
Rep Power: 9 
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. 

May 20, 2009, 10:30 

#4 
Senior Member
N/A
Join Date: Mar 2009
Posts: 189
Rep Power: 9 
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. 

May 25, 2009, 04:46 

#5  
New Member
Join Date: May 2009
Posts: 9
Rep Power: 9 
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). Code:
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 Quote:
Quote:
Quote:
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 

May 25, 2009, 06:34 

#6 
Senior Member
andy
Join Date: May 2009
Posts: 129
Rep Power: 9 
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.  nonperiodic 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. 

May 25, 2009, 10:05 

#7  
New Member
Join Date: May 2009
Posts: 9
Rep Power: 9 
Thank you for this answer.
Quote:
Quote:
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 

May 25, 2009, 11:35 

#8 
Senior Member
andy
Join Date: May 2009
Posts: 129
Rep Power: 9 
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 nonlinear 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. 

May 26, 2009, 08:17 

#9  
New Member
Join Date: May 2009
Posts: 9
Rep Power: 9 
Sorry to be late
I am full of work today... Thank you very much for all this help Quote:
In a staggered grid, there are the following grid :
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 

Thread Tools  
Display Modes  

