Hi all!

I am extending an FEM code for incompressible viscous fluid flows. I am just stucking in the pressure poisson equation.

The FEM discretization looks like: BT*ML-1*B p = BT*u

1. Assuming the above equation is assambled for one element, the pressure vector must have 3 components to be consistent with the RHS (3-components vector)... But as far as I know, the pressure is an scalar element value ???

2. I use linear shape functions for the velocity. After research I think that a constant shape function for the pressure is better. Isn't it?

Big thanks in advance for any comments.

boni

In the previous post: B is the gradient matrix. ML is the lumped mass matrix. BT is the divergence matrix.

BT*ML-1*B is the Laplace operator.

Thanks again for any suggestions after reading.

pressure is a scalar and so is your poisson equation. The problem is with your RHS; it should be a scalar. Start with the Navier Stokes equation and apply the divergence operator; or check out any book on basic fluid mechanics (and/or CFD book)

Adrin Gharakhani

Thank you for your help!

The RHS is for sure a vector. And the laplace operator (LHS) is also a matrix. So, it's not consistent treating the poisson equation as a scalar equation. More ever, that equation should be assambled in a Ax=b system for whole system...

thanks again!!

The RHS the way you have it is for sure a vector, but it doesn't mean it's correct! The pressure poisson equation is obtained by applying the divergence to NS.

Adrin Gharakhani

I got you. You are right, the divergence of the velocity vector is a scalar. But now, if the whole PPE is scalar, why people do not solve if directly? why the assamble it?

Thank you for posting!!

What do you mean "solving it directly" vs. "assembling"? These are two separate and unrelated sub-processes of the whole solution ...

Adrin Gharakhani

well, I am trying to use a fractional-step method. First, the mom. eq. is calculated neglecting pressure terms. Second, the pressure (poisson) eq. is calculated by using intermediate values from step 1. A later step updates the velocity.

Step one is done explicit: a = (f_external - f_internal)/m. This is done for each node...

The second step however must be done implicit (poisson eq. is a elliptic eq.). Here, I undestand that PPE must be assambled for solution...Ax=b

Did I understand the whole process wrong???

I strongly recommend you read the basics of FEM (patiently) to understand the underlying processes. Matrix assembly, to put it in the most primitive way, is the discrete version of the weighted integral form of the poisson equation. Each finite element is just one sub-volume that satisfies the equation _only_ at that sub-volume. But you need to satisfy the equation for the whole volume, so you need each finite element to "see" the influence of its neighbors (as a local boundary condition, to put it simply) - and as a result you end up "assembling" the matrix (i.e., sum the influence of all elements on all other elements, since, as you put it, this is an elliptic equation). Then you solve this matrix. How you solve it is up to you. It can be a direct approach, or any other whole host of iterative schemes.

Again, I urge you to read a basic book to understand the concept before you jump into the much more complex fluid mechanics problem.

Adrin Gharakhani

wow, that hurts me. That beginner I am not.

For your information, an explicit dynamic FEM does not need to assamble any matrices asumming the mass is lumped. An implicit dynamic FEM must. And I am new with implicit FEM coding.

Once again, my doubt initially posted here was about the dimension of matrices and vectors of the FEM-discretized pressure Poisson equation, which is solved implicit (again, here i am new). For solving it I will use a conjugate gradient interative solver.

I am sure the answer here is easier than coming back to continuum mechanics basics where virtual work and variational issues makes one mad.

No insult or ill-will was intended, I assure you.

Adrin Gharakhani

Hi!

I use the method of weighted residual FEM in solving laminar convection using the primitive variables. I know about the pressure Poisson equation superficially.

From your first post, I have few questions.

How can you get the continuous pressure solutions across element edges or let them be discontinuous? (I use out-of-date linear pressure, quadratic velocity and temperature elements.)

How do you impose the physical boundary conditions into the pressure Poisson equation?

Thanks

Hi, thanks for posting!

For velocity I use linear shape functions and for pressure I try now with constant ones.

The boundaries conditions are indeed absolutely important, but I still have another doubt. And here is it:

I use also weighted residual FEM and the PPE posted here results in a scalar equation for the element pressure.

Further, I use a fractional step method, which calculates the intermidiate velocity with explicit time integration. And for the pressure, I must solve it implicitly. Means, I need to calculate the "global" pressure. I guess, I need to assamble the PPE in a matrix system like Ax=b. But how?? if the PPE is a scalar equation...

How do solve your PPE?

After this I must of course touch BC issues.

Thanks for any help!

Hi! cfd-newbie

I solve the set of mass, momentum, and energy equations by weighted residual FEM directly and only steady state problems. So, I have never solve ppe and avoid to modify the governing equations.

Now, I have an idea (different from fractional step) to solve transient problems. Firstly, I use it to solve 1-D transient heat conduction efficiently. And I will try to extend it to solve 1-D phase change and transient 2-D laminar convection, respectively.

In my opinion, your fractional step method may be work for transient problems. But I expect something as follow.

In any flow field, pressure must be varied continuously. If you use constant pressure elements, you cannot represent the continuous pressure distributions. If you still need to use these elements, assembling element matrices maybe unnecessary (the discontinuity bring about special term in assembling procedure) but the pressure distributions on each element should be directly computed from velocity.

I want to recapitulate that this is my opinion only.

Hi worasit, thanks again!

yes, I want to use the code for transient problems. I have a explicit dynamic solver which I want to extend for incompressible viscous fluid flows.

I meant, constant shape fuctions for pressure, which means that the pressure is same along the element. Of course pressure may vary along the whole system since it is directly calculated from the velocity.

The FEM discretization of PPE looks like:CT*ML-1*C*p = CT*v The vector and matrix dimensions for a 4-node quadrilateral element are then: (1x8)(8x8)(8x1)(p=1x1) = (1x8)(8x1) Please correct me is dimensions are wrong.

However, the PPE is a elliptic equation which forces one to calculate it globaly. Also boundary values must be known for solution. As read in papers, here people solve PPE implicit by using interative solvers...Exactly this is my doubt, how they do it technically?

1000 thanks for your postings!!! Boni

From your previous post, you mentioned that 'people solve PPE implicitly by using iterative solvers'.

I think you have some paper about solving PPE. You must find out how they do.

If not, I suggest you to use linear pressure elements and start to derive element equations carefully.

I don't know how can you derive the 8x8 matrices by constant pressure elements. (They must come from two velocity components at each node.)

But the equation you expressed is not the weighted residual equation. So, you don't know the problem about discontinuity pressure at element edges in assembling procedure.

I'm glad to discuss with you!

Yes, I think I got to derive the equations carefully again. Actually I did it many times... I was hoping somebody did the same mistake (if any) so that I come faster forward. BTW the 8x8 matrix is the lamped mass matrix and , you are right, 4 nodes, 2 components.

I also can use linear shape functions for pressure. I just read in papers that it's recommended to use shape functions with one order lower for pressure than for velocity.

Did I express the PPE in Galerkin FEM? My pardon for the definition mistake.

OK pal, thanks a lot for your help!

As you known, a paper in 1975 suggested to use pressure shape functions one order lower than velocities. In my work, I use linear pressure, quadratic velocity and temperature element. I have 21 unknowns per element. And surely, I must construct 21 equations for them.

But your approach is not. You neglect pressure to find intermediate velocity explicitly. Then, you compute pressure distributions by implicit method.

I don't know whether your weighted residual is the same version with me. Because I must performing by-parts integration and use Gauss' theorem. I expect to see the order of pressure derivetives to be reduced by one.

I guess (I've never solved PPE) that your equation is 'underived' Galerkin FEM because I didn't see the term that I expected and the vector of external loads (physical BC's) from Gauss' theorem.

If your Galerkin FEM differ from mine, I can't tell you anymore.

Good Luck

Excuse me!

I forget some information you gived. Your equation is derived because you have divergence and gradient matrices in your equation.

But the term at boundary still missing. I think RHS of your equation does not contain it.

The solution of the problem posted here is as following.

Ignoring BC for instance, the 2D FEM discretized PPE looks globally like this:

[CxT*Mx-1*Cx + CyT*My-1*Cy] p = CxT*u + CyT*v

Dimensions: ========== - Inverse of Mx and My [numnodes x numnodes] - CxT and CyT [numelems x numnodes] - p [numelems] - velocity components u and v [numnodes]

numnodes=total number of nodes numelems=total number of elements

The system to solve is then Ax=b A [numelems x numelems] x [numelems] b [numelens]

So that is for now. I am still fighting to get an ICCG solver running but hope above info can help a future PPE-programmer.