Tag: finite element methods

Reproducing a transport instability in convection-diffusion equation

Drawing from Larson-Bengzon FEM book, I wanted to experiment with transport instabilities. It looks there might be an instability in my ocean-ice model but before being able to address that, I wanted to wrap my head around the 1D simplest example one could find. And that is the Convection-Diffusion Equation:

(1)   \begin{equation*} \begin{split} - \epsilon \Delta u + b \cdot \nabla u &= f \ \ in \ \Omega \\ u &= 0 \ \ on \ \partial \Omega \end{split} \end{equation*}

The first term -\epsilon \Delta u is responsible of the smearing of the velocity field u proportionally to \epsilon, and is thus called the diffusion term. Intuitively, it controls how much the neighbors of a given point x are influenced by the behavior of x; how much velocities (or temperatures: think of whatever you want!) diffuse in the domain.

The second term b \cdot \nabla u controls how much the velocities u are transported in the direction of the vector field b, and is thus called the convective term. A requirement for this problem to be well-posed is that \nabla \cdot b = 0 — otherwise it would mean that we allow velocities to vanish or to come into existence.

The instability is particularly likely to happen close to a boundary where a Dirichlet boundary condition is enforced. The problem is not localized only close to the boundary though, as fluctuations can travel throughout the whole domain and lead the whole solution astray.

Transport instability in 1D

The simplest convection-diffusion equation in 1D has the following form:

(2)   \begin{equation*} \epsilon u_{xx} + u_x = 1 \ \ in \ (0,1), \ \ u(0) = u(1) = 0 \end{equation*}

whose solution, for small \epsilon is approximately just u = x. This goes well with the boundary condition at 0, but not with the condition at 1, where the solution needs to go down to 0 quite abruptly to satisfy the boundary condition.

It’s easy to simulate the scenario with FEniCS and get this result (with \epsilon = 0.01 and the unit interval divided in 10 nodes):

transport-instability

in which we can infer two different trends: one with odd points and one with even points! In fact, if we discretize the 1D equation with finite elements we obtain:

(3)   \begin{equation*} \epsilon \frac{u_{i+1}-2u_i-u_{i-1}}{h^2} + \frac{u_{i+1}-u_{i-1}}{2h} = 1$ \end{equation*}

Continue reading “Reproducing a transport instability in convection-diffusion equation”

How do Dirichlet and Neumann boundary conditions affect Finite Element Methods variational formulations?

To solve a classical second-order differential problem

    \begin{equation*} -(au')' + bu' + cu = f \ in \ \Omega \end{equation*}

with FEM, we first need to derive its weak formulation. This is achieved by multiplying the equation by a test function \phi and then integrating by parts to get rid of second order derivatives:

(1)   \begin{equation*} \begin{split} 0 &= \int_\Omega ((-a u')' + b u' + c u - f) \phi dx \\ &= \underbrace{\int_\Omega ((a u' \phi' + b u' \phi + c u \phi) dx}_{a(u, \phi)} \underbrace{- \int_\Omega f \phi dx - (a u' \phi)|_{\partial\Omega}}_{L(\phi)} \end{split} \end{equation*}

A typical FEM problem then reads like:

    \begin{equation*} \begin{split} \text{Find } u \in H_0'(\Omega) \ s.t. \ a(u, \phi) + L(\phi) = 0 \ \ \forall \phi \in H_0'(\Omega), \\ \text{where } H_0'(\Omega) = \{ v: \Omega \rightarrow \mathbb{R} : \int_0^1v^2(x) + v'(x)^2 dx < \infty \}. \end{split} \end{equation*}

What is the difference between imposing Dirichlet boundary conditions (ex. u(\partial \Omega) = k) and Neumann ones (u'(\partial \Omega) = k(x)) from a math perspective? Dirichlet conditions go into the definition of the space H_0', while Neumann conditions do not. Neumann conditions only affect the variational problem formulation straight away.

Continue reading “How do Dirichlet and Neumann boundary conditions affect Finite Element Methods variational formulations?”

What is the difference between Finite Differences and Finite Element Methods?

With Finite Differences, we discretize space (i.e. we put  a grid on it) and we seek the values of the solution function at the mesh points. We still solve a discretized differential problem.

With Finite Elements, we approximate the solution as a (finite) sum of functions defined on the discretized space. These functions make up a basis of the space, and the most commonly used are the hat functions. We end up with a linear system whose unknowns are the weights associated with each of the basis functions: i.e., how much does each basis function count for out particular solution to our particular problem?

Brutally, it is finding the value of the solution function at grid points (finite differences) vs the weight of the linear combinations of the hat functions (finite elements).