hi,

I am trying to implement projection 2 algorithm of Gresho using finite element method. I would like to discuss some points with anybody who already implemented this scheme. I would be very happy if you can contact me.

Do you mean the incremental scheme? What are your questions?

I am talking about the projection 2 that gresho presents in his paper " on the theory of semi-implicit projection methods for viscous incompressible flow and its implementation via a finite element method that also introduces a nearly consistent mass matrix Part 2: implementation "

I try to implement the scheme using a Q2Q1 element (quadratic shape functions for velocities and linear shape functions for pressure). I am trying to debug it using a test problem. It is the classical channel flow with parabolic velocity profile. but i don't get logical results for the velocities at the first time step. the pressures look quite logical because pressure is linearly decreasing starting from inlet to outlet as it is in the analytical results. but I wonder why it doesn't give logical results for the velocities. my first question is do you any boundary condition for the solution of laplace equation (CT*ML-1*C)*PHI=CT*U~n+1-gn+1 (in discrete form as it appears in his paper) second question is about the term (CT*ML-1*C). How does it look like for one element, or if you already implemented this scheme is it possible to verify it for one element from the results of your code. In fact I still have some tiny questions in my mind but i keep them for my self.

I will try to answer your questions, though I have not personally used Q2Q1.

(CT*ML-1*C) is the so-called consistent Laplace matrix with M usually being properly diagnosed or lumped, since otherwise it is computationally unfeasible. Even using lumped M, the entries of this matrix are larger than the traditional Laplace matrix. I am not sure that it can be assembled in the element level. For example, for a given node, its support for the traditional Laplace is the intermediate neighbor nodes though element connectivity as we all know. But for the consistent matrix, the support involves one more level of the elements. This probably answers your second question.

The implied boundary condition for the pressure is dPHI/dn=0 where n the unit vector normal to the boundary. Since PHI is the increment of pressure (probably scaled by dt), this condition means dP^n+1/dn = dP^n/dn.

Hope this helps.

Tony

thanks for your suggestions and directions... what i have done up to now is to consider no boundary condition for the pressure poisson equation (CT*ML-1*C)*PHI = CT*U~n+1 - g_n+1 because in gresho's and chorin's paper and also in gresho book, it says that the boundary conditions are automatically built in this equation, therefore i don't implement any boundary condition for this equation. and the laplace operator i compute it globally as it is mentioned in the appendix of gresho's book not in element level...

another thing that i wonder about is the g term. it say that g term is the complete gradient matrix tranposed times the u_n+1 g_n+1 = CT * u_n+1. i would like to know why it appears in all of the gresho's papers but not in others. and what does it physically represent.. and did you also implement this term in your code..

and lastly i would like to know which discretization did you implemeted in your code, i guess Q1Q0 ??

thanks a lot

ilker

Ilker,

My experience with this issue was using Taylor-Hood type elements (P2P1). The consistent matrix was tested as the pressure preconditioner in Uzawa iteration. It worked fine.

I bet the ¡°g¡± term is related to some BC¡¯s, and you may ignore it in your case.

As you mentioned earlier that the pressure seems correct but the velocity does not. My guess is that something might be wrong with the gradient operator in the projection step.

Also do not forget to check if the prescribed velocity (no-slip?) BC is enforced.

Good luck.

Probably, i found the major error... up to now i was calculating the A matrix (CT*ML-1*C) after reducing the dirichlet boundary conditions from the C and lumped mass matrix (the contributions of the boundary nodes to gradient operator and lumped mass matrix are not taken into account ). But this should be done just before reducing the dirichlet boundary conditions from the equations.. at the end we get the same dimension for A either i calculate without reducing the C or after reducing the C. but the final matrices for A in both cases are totaly different. now for the first time step the velocity values look better than before ..

thanks for your suggestions

ilker basol

hi Tony,

Now i have some improvements in the results but i still have problems and questions. My question is about the K (viscous term). I checked Mark A.Christon's paper named "Dealing with pressure:FEM solution strategies for the pressure in the time-dependent Navier-Stokes equations" and he says that he uses K which is augmented by balancing tensor diffusivity (BTD) and uses trapezoidal rule for this term as below; (M+dt/2*K)U~n+1 = (M-dt/2*K)Un + dt*(.......) did you also augmented the K term and also implemented the trapezoidal rule. and the last question is whether you can solve with any reynolds number. Did you ever try reynolds number of one for example?

best regards, ilker

Yes, I did try very low Reynolds numbers. It is very stable as long as you treat the viscous term implicitly, like the Crank Nicolson as you mentioned above.

I believe the trapezoidal rule is good enough for integrating the viscous term.

Good luck.