Hi all,

I have a question about Navier Stokes in FEM. The dissipative (viscous) term in Navier Stokes is something like

mu * d^2 u / dx^2

If I have quadratic basis functions, I can use it as it is. If I make linear basis functions, then I can not because the second derivative is null. Thus the weak form will be something with

mu * integral{ du/dx * n } ds (I don't mention the weighting function)

Integrated on the border of a cell.

The problem is: du/dx is not well defined on a cell border, because the derivative is constant on each cell, and at a cell border, the derivative is not continuous/derivable.

How do you solve it ? Thanx in advance, LBRS

A d^2 T / d x^2 term comes up in heat conduction. Luckily for you, almost all books on the finite element method that have ever been written describe the use of linear basis functions for the heat equation!

So, I suggest you read about how to apply the finite element method to the d^2 T / d x^2 term in the heat equation, then use the same ideas that you read about to deal with the viscous dissipation term in the NS equations.

There are other complexities in using the finite element method for flows (upwinding, 'locking', etc.) so reading a wide range of papers and books will be necessary to create a full NS solver.

Regards, Andy

"almost all books on the finite element method that have ever been written describe the use of linear basis functions for the heat equation! "

Yes, exactly. And in all of these books, they use the Green identity to get an integral on the domain.

But they seem to forget that actually, the green identity has a second term, an integral of (grad u)*n on the border of the domain.

My question is: where is that term gone ? Many articles does not answer that question clearly.

Not sure I fully understand what bothers you, since some details are missing. In the following I assume your mu is constant, and that you apply Galerkin method, so that the weight and basis functions are identical (let's call them N(I) for node I).

When you go through the usual FEM procedure, i.e., multiply by the weight function, integrate, use the basis function to approximate u, apply the Gauss divergence theorem... you end up with two terms:

(1) mu*u(I)*integral over the volume of { N(I),k * N(J),k }

(2) mu*u(I)*integral over the surface of { nk* N(I),k * N(J) }

where ,k is partial derivative wrt coordinate k (covariant if non Cartesian coordinates are used) and nk is the k-th component of the unit normal on the surface.

So, first of all, there should not be any problem with linear (or bi- or tri- linear) basis functions, as the function and its derivatives are well posed, both in the interior domain and on the element boundaries.

The other question you ask, where this second term disappears - it does not. However, when you have adiabatic BC, it is what is known as "natural BC", and you may note that in this case this term is identically 0. In the other common case of prescribed boundary temperature, this term is known, and may be eliminated from the solved system to reduce its size. In more general cases, however, e.g., a prescribed non-zero flux or even more general - you have to include this term. However, I do no see any special difficulty here.

I hope I helped, Rami

"However, when you have adiabatic BC, it is what is known as "natural BC", and you may note that in this case this term is identically 0"

That is exactly the problem ! The adiabatic BC applies on boundaries. If you take a cell in the middle of your domain, which does not touch any boundary condition, then I don't understand how this term can disappear.

In some FEM books, they just take the weak FEM formulation from the global approach (galerkin) and use it for local approach (FEM). But you cannot apply a BC to a cell in the middle of your domain, no ?

> If you take a cell in the middle of your domain, which does not touch any boundary condition, then I don't understand how this term can disappear.

This is because you are using a continuous space so the terms on either side of an interior surface exactly cancel. This does not happen, for instance, with the discontinuous Galerkin method where flux terms show up on all the interior boundaries.

"the terms on either side of an interior surface exactly cancel"

I don't agree. If you integrate on a cell A, you don't add a term coming from the next cell B on the other side of the surface.

Jed is right. When you assemble neighboring elements (assuming continuous Galerkin) the internal surfaces are cancelling each other. Check my former post, and pay attention to the second term. Everything on both sides of a surface are identical, except to the opposite signs on the normal nk.

"When you assemble neighboring elements " Do you mean when you assemble the big (global) matrix from the small matrices (from local finite element) ?

>"When you assemble neighboring elements " Do you mean when you assemble the big (global) matrix from the small matrices (from local finite element) ?

yes...

Ok, then it helps, thanx