2 Computational continuum fluid dynamics

2.1 Finite Element Method

2.1.1 Stationary Convection-Diffusion Equation

2.1.1.1 Strong & weak forms

Strong form of the stationary convection-diffusion equation

The strong form of the stationary convection-diffusion equation is: \[ \nabla\cdot(\vec{u}\Phi)-\nabla\cdot(\epsilon\nabla\Phi)=\vec{f} \]

Peclet Number: $\epsilon=\frac{1}{\text{Pe}}$

For boundaries we have: \[ \begin{align} \Phi(x) &= g_1(x) &, x \in \Gamma_1 &,\text{Dirichlet BC}\\ \nabla\Phi\cdot\vec{n} &= g_2(x) &, x\in \Gamma_2 &,\text{Neumann BC}\\ \end{align} \]

Here, $\Omega=\Gamma_1\cap\Gamma_2$ and $\partial\Omega=\Gamma_1\cup\Gamma_2$.

From the strong form, we can then derive the weak form by integrating the strong form over the full domain after multiplying with a test function $w$: \[ \int_\Omega w\left[\nabla\cdot(\vec{u}\Phi) - \nabla\cdot(\epsilon\nabla\Phi)\ dx\right] = \int_\Omega wf\ dx \]

We then apply the divergence theorem to the higher order terms, i.e. the term $\nabla\cdot(\epsilon\nabla\Phi)$. After rewriting using the divergence theorem, we arrive at the weak form.

Weak form of the stationary convection-diffusion equation

Given $f\in F$, $\Phi\in V$ and $w\in W$, the strong form of the equation can be derived through integrating both sides of the strong form equation across the domain and then removing all higher-order terms through the divergence theorem. The resulting weak form equation for the stationary convection-diffusion equation is then: \[ \int_\Omega w\nabla\cdot(\vec{u}\Phi) + \nabla w\cdot(\epsilon\nabla\Phi)\ dx = \int_\Omega wf\ dx + \int_{\Gamma_2} w\epsilon g_2\ dx \]

We designate the left-hand side (which is of bi-linear form) as $B(w,\Phi)$ and the right-hand side as $F(w)$.

Function spaces :\[ \begin{align} F&=\left\{f: \Omega \rightarrow \mathbb{R}, \int_\Omega f^2\ dx < \infty\right\} &, f\in\mathcal{L}^2(\Omega)\\ V&=\left\{v: \Omega\rightarrow\mathbb{R}, \int_\Omega v^2\ dx < \infty \right\} &, v\in\mathcal{L}^2(\Omega)\\ & \left\{\int_\Omega ||\nabla v||^2\ dx < \infty\right\} &, \nabla v\in\mathcal{L}^2(\Omega) \end{align} \]

The problem now becomes to find $\Phi\in V$ such that the weak form equation holds.

2.1.1.2 Galerkin approximation

Consider $V^h$ V$, $W^h \subset W$ as the discrete subspace, which are a subset of their infinite-dimensional counterparts. We call $V^h$ the trial space and $W^h$ the test space. The Galerkin approximation states that the trial functions and the test functions are the same (i.e. come from the same space), except at the essential boundary.

\[ \Phi^h = w^h + v_0\ , v_0\big|_{\Gamma_1}=g_1 \]

Given $f\in F$, we want to find $\Phi^h \in V^h$ such that: \[ B(w_h, \Phi_n)=F(w_h)\ \forall w_h \in W^h \] This is called the “Galerkin problem”.

2.1.1.3 Finite Elements

To solve the Galerkin problem, we follow a process:

Discretize the domain into elements with $N$ unknowns
Choose a basis for $W^h$, i.e. choose ${w_i(\vec{x})}^N_{i=1}$ such that they span $W^h$. Find coefficients $c_j$ such that $\Phi(\vec{x})\approx\Phi_n(\vec{x})=\sum_{i=1}^N c_jw_j(\vec{x}) + v_0(\vec{x})$
Solve for coefficients $c_j$ using $B\left(w_i, \sum_{i=1}^N c_jw_j(\vec{x}) + v_0(\vec{x})\right) = F(w_i)$

In practice it can be useful to choose basis functions that are only non-zero in the proximity of the element, as this leads to a sparse matrix system instead of a dense one.

We find several properties of the FE solution:

For $f=0$, the scheme is equivalent to using central difference approximation
For $u=0$ (i.e. pure diffusion), the scheme is nodally exact
For the “mesh Peclet number” $\text{Pe}_h\geq1$, the method produces unphysical oscillations. This comes from articifical negative diffusion (workaround: upwinding, or use Petrov-Galerkin instead of Bubnov-Galerkin)