Hello everybody!

I am working on a quite simple CFD-application (FEM). This should solve the simple (steady, incompressible, const. viscosity, const. density) stokes-eqn. in 2D, so my unknowns are u,v and the pressure p, the extenal forces (=gravity) usually set to zero. I have tried discretizing several approaches (with and without lambda), even eliminated the pressure with penalty-method. BUT: the solver fails unless i set all pressure to zero (homogeneus b.c.), and even then (with velocity at inflow or gravity) the mass conservation is not fulfilled! (speed sum at inflow is much bigger than at outflow, or at symmetric problems it is lower in the middle) Or when i set a pressure distibution, which should give me a u and v result, the velocities remain zero.

Same picture with penalty-solution: though the directions of the speeds seem correct and "natural", the speed decreases continuously, so mass conservation is not fulfilled.

Please can anybody give me a hint whats the problem with my thoughts or a helpful electronic source? (i have no access to libraries at the moment, books are really expensive and i can only buy them when they contain THE SOLUTION to my problem)

I think if somebody knows it, i would need a step-by-step-explanation (for steady viscous stokes) from the diff-eqn. up to the compilation of the stiffness matrix. Guess that would solve it all for me.

And before somebody brings it up: yes, my solver works correct (i solved mech. problems and they were correct up to the 8th number after the komma), and i checked the matrix compilation by hand (well, some hours and tons of paper went in there)

Thanks for your help in advance! (And I'm sorry for my bad scientific english)

Juergen,

it is a little difficult to directly "hit" your problem without further knowledge about your solution method.

- Coupling of velocity and pressure: "coupled" or "segregated" approach?

- If "segregated", how is the pressure computed? (the pressure "matrix" ?)

- What kind of elements are used? Are the elements LBB-stable?

- Does your approach work for Stoke's flow ? (Re-number -> 0)

Regards

Chris

Thank you Chris for taking care of my problem! But after some days of reading papers and considering the hardware available to me (and the computation time) I found the "penalty-method" approach suits my needs best. I got rid of my mass-conservation-problem, which was caused by my wrong interpretation of the "penalty-parameter" and some minor mistakes in my formulas.

I went through some papers during the last few days and found some "errors" of mine which I want post here, so others can avoid them. These errors were quite basic and therefore SERIOUS. If someone finds incorrectness in the following please do correct it, so others can take advantage of the collected knowledge.

I tried a coupled approach at first, and I found out you have to take a higher order approximation for the velocity than for the pressure. So if you want linear pressure you should take quadratic velocities. These elements may not be LBB-compliant but do work.

I expected the stiffness matrix to be nice "symmetric, positive definite" as in linear mech problems, but this is very difficult (or impossible?) to obtain. -> The more general the equations become, the less symmetric the matrix will be.

IF you manage to set up your linear equation set and solve it (which I couldn't), you will find the instability of pressure, so you need a stabilizing procedure. Several ways for this are possible, but I didn't go for that approach any further.

At this point I almost started crying, but then I remebered my first try: the penalty method. If you can make so many assumptions as I can, then this is the nicest solution, I think. It reduces your unknowns (no pressure!), which means less memory and less computation required. The pressure is not lost, but can be computed in a post-processing-step. There is also the "lambda"-value (in some papers called epsilon, which is the inverse value: lambda=1/epsilon). Setting lambda to almost infinite (or epsilon to almost zero) gives the incompressible solution. My calculations showed a 100% match to the analytic results!!

One last thing is there still open: some of the papers say, the lambda-value is the strain-viscosity (in opposite to the usual shear viscosity). For my applications (polymer flow) there are measured strain-viscosity-values available. But oher papers do just call this the "penalty-parameter", without any physical meaning.

Please can you tell me what is correct?

Thanks for your help!

Juergen

Just to describe the basics of the penalty method for incompressible fluids: the continuity equation "div(v)=0" is somewhat relaxed and taken as

eps * p + div(v) = 0; with esp << 1.

The penalty parameter "eps" in combination with the pressure "p" allows a small compressibility. Often the term is justified by introducing the speed of sound - but you can as well see the term as some kind of "trick" to simplify the equations. Any further interpretation may lead to an interesting story - and of course sounds better than "i took the numerical parameter eps as ...". Anyway, in my experience "eps" is just a numerical parameter and it can depend on the type of flow, the flow geometry and it for sure depends on the Reynolds number.

Well, step by step i'm getting into that math-trick-idea. Obviously, some of my reference papers have a symbol mix: they use the lambda-symbol for both the strain viscosity AND the inverse penalty-parameter, i guess thats where i thought they were the same, but they surely aren't. (just tried a calculation with water-values, where lambda(=strain vicosity) is 3x of the shear visosity, and the results are far from correct.

As I see you are a really pro in this fluid-fem-stuff, so may I beg another hint?

As I already mentioned, my velocities do match the analytic results exactly (at least up to the 8th number beyond the komma). So I tried the pressure postprocessing, what I thought to be an easy task. In my understanding, it is just multiplying the velocity derivates with the lambda value (or dividing the deriates by epsion). I read: p=(1/eps)*(u i,i) or in my formulation p=lambda*(u i,i) where (u i,i) is the sum of derivates: du/dx+dv/dy+dw/dz so I calculated the derivates per element (while my velocities are linear, 3-node-triangles). The sums of derivates show an almost alternating + - pattern, and this multiplied by the penalty-param lambda (~1e7) gives just crap. So there IS an error in my understanding of the post-processing, not in my code. I checked the results by hand, which gave me the same + - alternation.

Please can you tell me how I can get my pressure? And is it possible to calcutate it by node, or do I have to go to a lower order than the velocities are? I found a (too) brief description of linear-v-constant-p (what I tried) and quadratic-v-linear-p (what I want at least at the end of my work), but there was also mentioned a linear-v-linear-p way, but not described in any word.

Thanks!

Juergen

Uhh - thanks for the compliment; actually it is a very long time ago since I coded FEM solvers.

Anyway and without looking at any old documents: As far as I remember most people using the penalty approach used a single pressure node within the element or used pressure nodes inside the element (not on the element boundary). Thus the pressure approximation was of lower order than the velocity approximation and is was simple to compute the pressure "matrix". And most of them used a "reduced integration" for the pressure term (less gauss nodes than for the velocity matrices). Have you tried to change the parameter "eps"? And have you tried to change the boundary conditions? If you have dirichlet boundary conditions everywhere - did you enforce the pressure value at a single node? Did you try to compute a simple channel flow with neumann boundary conditions at the outlet? What happended?

So far so good: just based on your last reply I can tell you what I did and what I did NOT

I my actual testing code I am working with 2D, linear 3-node triangles. Because they are simple and I can check results by hand. My "mesh" consists of a 10x10 node rectangular area, regulary connected by 162 elements.

Velocities work well. For my testing calculations I usually set homogeneous b.c. for velocities on the upper and lower line of nodes. For "nice" flow I just impose gravity in flow direction (from left to right)

----------

**********

**********

**********

**********

**********

**********

**********

**********

----------

This gives me a nice parabolic velocity distribution, that matches the analytic results exactly.

Now to the part that does not work: I thought, I could calculate the pressure by simple multiplying. So I calculate the velcity derivates in the middle of each element (NOT per node). This should be enough lower order, and is is quite easy. (And correct, as I checked) Then, I calculated the sum for each element: du/dx+dv/dy

(this sum shows already the + - alternation)

And this sum I multiplied with (1/eps).

So, I didn't do ANY integration for pressure. And I think, this is my error! I don't know where to integrate, should I assume another (virtual) mesh with "middle nodes"? This would get me some new form functions, which I actually do not have for my virtual pressure node.

I know, the result for the above described flow should be an almost zero pressure distribution. But, as I mentioned in my last post, I get huge (~1e4) + - pressure values.

As I put the pices together, I missed out any integration. And as I understood the "pressure post-processing" as a simple multiplication, I didn't set up a matrix, or solved anything, I just worked value by value, of course they were stored in a (n x 1) matrix.

So to put it in a single sentence: what should I integrate and what should I solve?

Thanks for spending your time with this, as my brainpower is obviously not enough for getting this on my own.

Jürgen,

I am afraid you still have some of the basics wrong (forgive me if I am wrong in stating this):

Velocity and pressure are COUPLED in incompressible flows. That means you have to solve BOTH the momentun and the continuity equation. You have 3 unknowns (u,v,p) and 3 equations (x- and y-momentum and continuity). As far as I can see you just tried to decouple the equations ...

OK, the idea of the penalty method is to replace the pressure in the weak formulation of the momentum equations by the divergence of the velocity and the penalty parameter "eps". You then have only two equations to solve and you only get (u,v). In order to do so you need velocity shape functions.

After that you can calculate the pressure by "solving" the continuity equation. In order to so you need pressure shape functions within the elements. Typically you end up with an equation simular to "eps * Mp * P = trans(C) * V" where "Mp" is the pressure "mass matrix", "trans(C)" is the divergence matrix and "V" and "P" are vectors containing the node values of (u,v,p).

The matrix "Mp" is generated by integration (here you have it) using the pressure shape functions; reduced integration recommended.

There are variants of the penalty method, one of them the "discrete penalty method" where the matrix formulation of the continuity equation is inserted in the matrix formulation of the momentum equation.

I strongly suggest to use a test case where the solution for the pressure is different from zero. Try a simple channel flow with Re=1 using a very high viscosity (thus you will get pressure values much higher than zero).

I also strongly suggest you either do a WWW-research (google with "incompressible, finite, element, penalty") or try to find a textbook.

Regards, Chris

Hi Chris,

I fear you are right, I really didn't get some of the basics. I was just happy that I got the velcities so quick (from zero to complete FEM-program in 2 days) so I thought I needn't to go deeper into the background.

I want to thank you for your spent time, and especially for that last hint. I think that line with the formula and the symbol explanation gave me more "aha!" than the few dozen pdfs I read on the topic.

I will go to the beginning and try to finally get the pressure, based on your formulation. I will post some beer for you when I'm finished!

Thanks for your help!

Juergen

I honor those with power in their brain!