Hey all,

I am just curious as to how one would develop an implicit solver for the incompressible Navier-Stokes equations, with finite differences to handle the spatial derivatives.

I am quite familiar with methods for ODEs (like backward Euler, Adams-Bashford/Adams-Moulton, BDFs, RKs, etc.). I have also seen the usual explicit Eulerian approach to solving NS, which differs from how one would solve an ODE. I am talking about the explicit method presented by Griebel, Dornseifer, Neunhoeffer in their book "Numerical Simulation in Fluid Dynamics". They way they describe it, one would take a forward-Euler step, and then solve the Poisson equation for pressure to make the system divergence-free.

I can't seem to find any decent explanations of how NS is solved with implicit methods, or with higher-order methods. I can think of two ways, but haven't been able to find info on either one of them. Either there is a way to reduce NS to an ODE, or there are special methods (like the explicit one by Griebel et al.) that only work with NS.

So, am I missing something? Can someone point me to some info?

Thanks.

Dave

You need to look into solving differential-algebraic equations. The constraint equation for strictly incompressible fluids makes the DAE system index-2. If you are not strictly incompressible then you need to determine the proper constraint manifold. Once you've done this, go to the book

http://www.amazon.com/Solving-Ordina...e=UTF8&s=books

Always ask yourself what is the physical basis for the constraint manifold and am I enforcing it. The derivation of the constraint manifold for non-isothermal flows is a mess.

Thanks. I had seen that the incompressible Navier-Stokes equations could be represented as a DAE. However, I was under the impression that it was really just a non-linear PDE. Certainly, all of the textbooks on CFD that I have seen spend a chapter or two on ODE solvers, but not DAE solvers. However, I haven't seen a form of NS that I can just feed to an ODE solver.

So, if the *incompresible* version of the Navier-Stokes equations are an index-2 DAE, then (as far as I know) the only way to solve them numerically would be to reduce it to an index-1 DAE or to a pure ODE. So before I buy your recommended text and bone up on my DAE solvers, I'd like to know which is it? In either case, I presume there is literature that does the reduction for me. Could you point me to such literature?

could you please give some information about concept "index" of PDE? (pls don't refer me to any paper or book, i only need overal veiw)

i see it only in methematical literature.

I believe the claim that Mr. Runge_Kutta made was that Navier-Stokes is a DAE, rather than a PDE (or at least, it's both). If it is a DAE, then index-2 makes sense. Though, as far as I know, there is no concept of an index for a PDE. For more info on what the index of a DAE is, have a look at section 2 of this guy:

http://se.wtb.tue.nl/~vanbeek/pub/mcmds01.pdf

I know that you don't want a paper, but that section is quite short an explains things better than I could in a raw text posting.

The strictly incompressible (isothermal) NSEs are what are called partial differential algebraic equations (PDAEs). Once you discretize space with finite {difference,volume,element} methods, you then have a system of DAEs. First off, this stuff is very complicated and very mathematical.

In DAEs, there are several different definitions for the index. One is the differential index and other is the perturbation index. They are measures of the sensitivity of the solution to perturbations introduced into the differential and algebraic equations. Higher index values imply a more difficult problem to solve. The isothermal NSEs are index-2. You can differentiate the constraint equation w.r.t. time and derive a new constraint equation. For the isothermal NSE, differentiate Div(u)=0 with the material derivative. You get the dilatation transport equation. This is related to the index-1 version of the NSE. Substitute in the index-2 constraint, Div(u)=0, and you get the pressure Poisson equation. This is called index reduction and it can lead to some pathological behavior because you are enforcing the derivative of the constraint but not the constraint.

I haven't mentioned the nastiest problem. When you have variable temperature, species, and chemical reactions, what is the exact form of the constraint manifold. This manifold is the one that suppresses acoustics. If you do not know the correct form of this manifold then what are you doing?? If you aren't suppressing the "fast time scale," why solve DAEs?? Also, the DAE system does nothing towards the diffusive and reactive terms.

http://scholar.google.com/scholar?q=...r=&btnG=Search

The Navier-Stokes equations are a parabolic PDE which is all you need to really know. The best place to find out how to solve the Navier-Stokes equations using finite differences is probably the book by Ferziger & Peric - there's quite a detailed discussion and some example codes which can be downloaded from Springer's website.

Basically there are a number of ways to solve the incompressible Navier-Stokes equations. One possible method is the "Jacobian Free Krylov Method". This method entails discretizing the equations in space and time and iterating at each timestep to obtain a solution of the nonlinear algebraic system. This method is relatively easy to understand - you have a system of nonlinear algebraic equations obtained by discretizing the equations which you solve iteratively (no pressure projection).

More common methods seggragate the equations to form an iterative scheme; e.g.

(1) Estimate an increment to get the velocity at timelevel n+1 using the momentum equation (minus the pressure term).

(2) Project (Hodge decomposition) the result to obtain an increment to the velocity which ensures that the new velocity field is divergence free.

(3) repeat until converged.

Have a look at these lecture notes

http://www.math.rug.nl/~veldman/Coll...ictaat2006.pdf

>The isothermal NSEs are index-2

i think that, say that NS eq. is index-2 is not correct!

- NS is PDAE, index with respect to whose variable?

- index of PDAE is function of initial and boundary conditions, so withouth especification of them index is not known.

- NS is non-linear and classic analysis are (based on my knowledge) is limitted to linear/quasi-linear systems.

is it?

- if u look for spatially high oreder implicit, it is very simple, refer to compact high order methods for convection-diffusion and poisson eq., and use Projection or unsteady-simple, unsteady-simpler, PISO method for pressure-vel coupling (PISO has higher order time accuracy too)

- but high order (arbitrary) implicit temporal direction is difficult,

i don't knowe why this thread streched to index analysis (math not CFD). Index analysis of DAE is difficul, for PDAE it is more difficult, for non-linear PDAE is very difficult and usually analysis are limited to linear and quasi-linear systems,

while NS is non-linear PDAE. Also i don't think after such complex and tedious analysis you find desired high order implicit method (es. arbitrary order as for explicit RK).

- also i think that the pressent methods such as streamline-vorticity, projection methods are nothing else results of index reduction.

so my advice: don't put yours in such ambiguos trade (i have simpler cure). But if you are interested (i face such analysis in derivation of continous adjoint of NS eq.) i recommend L. Petzold's book and Barton's papers (SIAM J SC, give u good look, i have PDF, if u need notify)

But, why high order implicit in time?

i see some physical barier to this (due to temporal-parabolic nature) that give this guide that if you derive such method it can be very non-linear, stiff and ill-posed.

Probably aim of high order implicit in time is beter temporal accuracy (ok), so as cure:

1) decrease temporal step size and use avalible tools.

2) use present method with posteriori error estimation and adaptation, i guide you to seris of work in "Heidelberg University" and their opensource tool (Deal II library, very powerfull), they also publish a book, they present methods for high order solution of PDE and estimation of general functional. Also another work that i knowe is supra-convergent method with adjoint error recovery due to Giles and Pierce (SIAM Rev).

Maybe i'm not so addicted to math and CFD as other people here but i think that the simpler way (pedagogical) to obtain an implicit, 2° order in time, scheme for the N.S. equations is the fractional step method.

That is, the intermediate velocity field(V*) is decomposed in a divergence-free part (V(n+1)) and a curl-free part (Nabla(Phi)). Obviously the curl-free part is NOT the gradient of the pressure field (this is true only for first-order in time).

Once made such a decomposition, discretize the conv.term explicitly, the diff. term implicitly (it is the first source of instabilities) and integrate with Adam-Bashfort. In this way you don't need to handle non-linearities implicitly, you have a 2° order scheme and a much more wide region of stability than in a 2° order fully explicit scheme.

I'm sorry, of course i mean the projection method, not the fractional step.