Lecture 06 - Finite volumes

Examples: linear advection

Simplest hyperbolic PDE:

\[ \partial_t u + \mathbf{a} \cdot \nabla u = 0, \]

i.e. \(\mathbf{u} = u\), \(\mathbf{f}(\mathbf{u}) = \mathbf{a} u\), \(\mathbf{s}(\mathbf{u}) = 0\).

Examples: inviscid Burgers’ equation

Non-linear equation that can exhibit formation of shocks:

\[ \partial_t u + \partial_x \left( \frac{u^2}{2} \right) = 0, \] i.e. \(\mathbf{u} = u\), \(\mathbf{f}(\mathbf{u}) = \frac{u^2}{2}\), \(\mathbf{s}(\mathbf{u}) = 0\).

Examples: Compressible Euler equations in 2D

\[ \partial_t \rho + \nabla \cdot (\rho \mathbf{v}) = 0, \\ \partial_t (\rho \mathbf{v}) + \nabla \cdot (\rho \mathbf{v} \otimes \mathbf{v}) + \nabla p = 0, \\ \partial_t E + \nabla \cdot \left( (E + p) \mathbf{v} \right) = 0, \]

i.e. \(\mathbf{u} = (\rho, \rho \mathbf{v}, E)^T\), \(\mathbf{f}(\mathbf{u}) = (f_x, f_y)\), \(\mathbf{s}(\mathbf{u}) = 0\), flux components \(f_x, f_y\) are given by

\[ f_x = (\rho v_x, \rho v_x^2 + p, \rho v_x v_y, (E + p) v_x)^T, \\ f_y = (\rho v_y, \rho v_x v_y, \rho v_y^2 + p, (E + p) v_y)^T \]

Closure:

\[ E = \rho \varepsilon_{int} + \rho \frac{v^2}{2},\quad p=p(\rho, \varepsilon_{int}) \]

Why methods for unsteady equations?

In steady-state \(\partial_t \equiv 0\), the type of the problem depends on the regime:

• subsonic: elliptic
• supersonic: hyperbolic

Unsteady approaches more general.

Finite volumes

Compared to finite differences, FVM allows for

• complex geometries and meshes
• easier handling of discontinuities
• more robust conservation properties

Integral conservation law

Integrate PDE over control volume \(V_i\):

\[ \int_{V_i} \partial_t \mathbf{u} \, \mathrm{d}V + \int_{V_i} \nabla \cdot \mathbf{f}(\mathbf{u}) = \int_{V_i} \mathbf{s}(\mathbf{u}) \, \mathrm{d}V. \]

Fluxes at cell interfaces

Split the cell boundary \(\partial V_i\) into cell faces: \(\partial V_i = \cup_j f_{ij},\, j\in \mathcal{N}(i)\), \(\mathcal{N}(i)\) is the set of neighbors of cell \(i\), \(f_{ij} = \partial V_i \cap \partial V_j\).

Then we have \[ \int_{\partial V_i} \mathbf{f}(\mathbf{u}) \cdot \mathbf{n} = \sum_{j\in \mathcal{N}(i)} \int_{f_{ij}} \mathbf{f}(\mathbf{u})\cdot \mathbf{n} \, \mathrm{d}S. \]

Fluxes at cell interfaces

We can approximate the flux integral over a cell face \(f_{ij}\) as \[ \int_{f_{ij}} \mathbf{f}(\mathbf{u})\cdot \mathbf{n} \, \mathrm{d}S \approx g_{ij}(\mathbf{u}_{i},\mathbf{u}_{j}), \]

\(g_{ij}\) is a numerical approximation of the flux from cell \(i\) into cell \(j\).

Final semi-discrete form

Define

\[ \mathbf{u}_i = \frac{1}{V_i} \int_{V_i} \mathbf{u} \, \mathrm{d}V, \]

then we have the semi-discrete form of the conservation law:

\[ \partial_t \mathbf{u}_i + \frac{1}{V_i} \sum_{j \in \mathcal{N}(i)} g_{ij}(\mathbf{u}_i, \mathbf{u}_j) = \mathbf{s}_i. \]

Cell-centered vs node/vertex-centered

Assume we have a grid we would like to compute our solution on. What does that mean exactly?

There are two ways to take a grid and define a solution on it: cell-centered and vertex-centered.

Cell-centered

In a cell-centered approach, we define the solution at the center of the grid elements.

Riemann problem

For a 1-D hyperbolic conservation law

\[ \partial_t \mathbf{u} + \partial_x \mathbf{f}(\mathbf{u}) = 0 \]

define initial conditions

\[ \mathbf{u}(x, t=0) = \begin{cases} \mathbf{u}_L & x < 0 \\ \mathbf{u}_R & x > 0 \end{cases} \]

This is called the Riemann problem, we denote the solution as \(\mathbf{u}^{\mathcal{RP}}(x,t)\).

Example: linear advection

Simplest textbook example:

\[ \partial_t u + a \partial_x u = 0. \]

Solution to Riemann problem is given by

\[ u^{\mathcal{RP}}(x,t) = \begin{cases} u_L & x/t < a, \\ u_R & x/t > a. \end{cases} \]

In fact, any system of linear transport equations

\[ \partial_t \mathbf{u} + \mathbf{A} \partial_x \mathbf{u} = 0. \]

has a simple linear solution related to the eigenvalues and eigenvectors of \(\mathbf{A}\).

Godunov’s method

Take our conservation law (in 1D for simplicity)

\[ \partial_t \int_{x_{i-1/2}}^{x_{i+1/2}} \mathbf{u} \, \mathrm{d}x + \int_{x_{i-1/2}}^{x_{i+1/2}} \partial_x \mathbf{f}(\mathbf{u}) \mathrm{d}x = 0, \]

and integrate from time \(t^n\) to \(t^{n+1}\):

\[ \int_{t^n}^{t^{n+1}} \partial_t \int_{x_{i-1/2}}^{x_{i+1/2}} \mathbf{u} \, \mathrm{d}x \mathrm{d}t + \int_{t^n}^{t^{n+1}} \int_{x_{i-1/2}}^{x_{i+1/2}} \partial_x \mathbf{f}(\mathbf{u}) \mathrm{d}x \mathrm{d}t = 0 \implies \]

\[ \frac{\mathbf{u}_{i}^{n+1} - \mathbf{u}_{i}^{n}}{\Delta t} + \frac{1}{\Delta x} \int_{t^n}^{t^{n+1}} (\mathbf{f}(\mathbf{u}_{i+1/2}) - \mathbf{f}(\mathbf{u}_{i-1/2}))\mathrm{d}t = 0. \]

Godunov’s method

If we take \(\mathbf{u} \equiv const\) in each cell, then we have

1. Riemann problem at each interface
2. Solution of Riemann problem at \(x/t = 0\) stays constant \(\implies\) flux of Riemann problem solution also stays constant in time

So we can write \[ \int_{t^n}^{t^{n+1}} (\mathbf{f}(\mathbf{u}_{i+1/2}) - \mathbf{f}(\mathbf{u}_{i-1/2}))\mathrm{d}t = \Delta t (\mathbf{f}(u^{\mathcal{RP}}_{i+1/2}) - \mathbf{f}(u^{\mathcal{RP}}_{i-1/2})) \]

The birth of modern(-er) CFD

Roe’s method

The key insight is: better to solve Riemann problem approximately and gain solution accuracy elsewhere. Which approximation? Local linearization (at cell interface):

\[ \partial_t \mathbf{u} + \nabla \mathbf{f}(\mathbf{u}) = 0 \implies \partial_t \mathbf{u} + \mathbf{A}(\mathbf{u}_L, \mathbf{u}_R) \nabla \mathbf{u} = 0. \]

Conditions on \(\mathbf{A}\) (Roe matrix):

1. \(\mathbf{A}\) is diagonalizable with real eigenvalues (remain hyperbolic)
2. \(\mathbf{A}\) is consistent: \(\mathbf{A}(\mathbf{u}, \mathbf{u}) = \mathbf{A}(\mathbf{u})\)
3. \(\mathbf{f}(\mathbf{u}_R) - \mathbf{f}(\mathbf{u}_L) = \mathbf{A}(\mathbf{u}_L, \mathbf{u}_R)(\mathbf{u}_R - \mathbf{u}_L)\)

One way to define this matrix is as \(\mathbf{A} = \frac{\partial \mathbf{f}}{\partial \mathbf{u}} \left( \mathbf{u}_{avg} \right)\), where \(\mathbf{u}_{avg}\) is some average of the states \(\mathbf{u}_L\) and \(\mathbf{u}_R\).

Roe’s method

Recipe:

1. Compute \(\mathbf{u}_{avg}\), \(\mathbf{A}\)
2. \(\mathbf{A}\) is diagonalizable, i.e. \(\mathbf{A} = \mathbf{R} \mathbf{\Lambda} \mathbf{R}^{-1}\)
3. Compute \(|\mathbf{A}| := \mathbf{R} |\mathbf{\Lambda}| \mathbf{R}^{-1}\)

For Euler (1D):

• \(\rho_{Roe} = \sqrt{\rho_L \rho_R}\)
• \(u_{Roe} = \dfrac{\sqrt{\rho_L} u_L + \sqrt{\rho_R} u_R}{\sqrt{\rho_L} + \sqrt{\rho_R}}\)
• \(h_{Roe} = \dfrac{\sqrt{\rho_L} h_L + \sqrt{\rho_R} h_R}{\sqrt{\rho_L} + \sqrt{\rho_R}}\) (\(h = e + p/\rho\))

Roe’s method: summary

• Computationally expensive
• Not always stable (hypersonic flows):
  • if \(f'(u_L) < 0 < f'(u_R)\) (rarefaction) leads to unphysical solution
  • “Carbuncle phenomena”

Other Riemann solvers

• HLL (Harten, Lax, Leer): assumes we know two fastest wave speeds of Riemann problem, resolves only them: need wavespeed estimates, i.e. min/max eigenvalues of \(\mathbf{A}\)
  • cannot resolve contact waves
• HLLC
  • add contact wave (need speed estimate, etc.)
• Hybrid methods (Roe + HLL)
• AUSM, AUSM+
  • high-speed flows, splitting of transport due to pressure and convection
  • efficient (no eigendecomposition required)

TVD and Godunov’s theorem

MUSCL

• Reconstruct piecewise linear solution from neighbouring cell averages
• Use values at cell interfaces to compute fluxes

\[ \mathbf{u}_L = \mathbf{u}_i + \nabla \mathbf{u}_i \cdot (\mathbf{r}_m - \mathbf{r}_i),\quad \mathbf{u}_R = \mathbf{u}_j + \nabla \mathbf{u}_j \cdot (\mathbf{r}_m - \mathbf{r}_j) \]

• Gradient computed either based on linear interpolation or least-squares

Slope limiting

Replace \(\mathbf{u}_i + \nabla \mathbf{u}_i \cdot (\mathbf{r}_m - \mathbf{r}_i)\) with \(\mathbf{u}_i + \Psi \nabla \mathbf{u}_i \cdot (\mathbf{r}_m - \mathbf{r}_i)\).

\(\Psi = 0\) recovers 1st order scheme, \(\Psi\) - slope limiter. Restricts local solution slope to avoid oscillations.