Skip to content

Equation Types

PDE's can be classified as elliptic, parabolic, or hyperbolic. The type of equation will suggest the choice of an appropriate algorithm.

The equation type can change for different flow parameters, or even for different regions within the same domain.

20.4.1 Scalar Equations

All of the scalar equations in fluid dynamics can be written in canonical form as the following 2nd order PDE.

aϕxx+bϕxy+cϕyy+dϕx+eϕy+fϕ=g a \phi_{xx} + b \phi_{xy} + c \phi_{yy} + d \phi_x + e \phi_y + f \phi = g

where x,y x, y are any two independent variables. I.e. we could let x→t x \rightarrow t . In this form, we can identify the equation type by transforming into a quadratic form (derivation not shown)

a(d yd x)2+b(d yd x)+c=0 a \left( \dv{y}{x} \right) ^2 + b \left( \dv{y}{x} \right) + c = 0

The equation types are

  • Elliptic: b2<4ac b^2 < 4ac
  • Parabolic: b2=4ac b^2 = 4ac
  • Hyperbolic: b2>4ac b^2 > 4ac

Apparently, the equation type is going to have something to do with the types of roots of the scalar PDE.

Some examples:

  1. Laplace's Eq.

    ϕxx+ϕyy=0(0)2<4(1)(1)→elliptic \phi_{xx} + \phi_{yy} = 0 \qquad (0) ^2 < 4(1)(1) \rightarrow \text{elliptic}

    We know that typical solutions to Laplace's equation are entirely determined by the value at the boundary. The solutions are the smoothest possible functions that satisfy the boundary condition. There are no waves.
  2. Heat conduction equation

    θt=αθxx(0)2=4(α)(0)→parabolic \theta _t = \alpha \theta_{xx} \qquad (0)^2 = 4 (\alpha)(0) \rightarrow \text{parabolic}

    We know that the typical solutions have diffusion-like behavior.
  3. Linear advection equation

    ut+aux=0 u_t + a u_x = 0

    We're not already in canonical form, but we can get there with some algebraic manipulation

    0=utt+auxy=utt+a(ut)x=utt−a2uxx 0 = u_{tt} + a u_{xy} = u_{tt} + a (u_t)_x = u_{tt} - a^2 u_{xx}

    →0>4(1)(−a2)→hyperbolic \rightarrow 0 > 4(1)(-a^2) \rightarrow \text{hyperbolic}

    Typical solutions to the linear advection equation are waves which propagate at a set speed.

Let's look back at the linear potential equation

(1−M∞2)ϕxx+ϕyy=0 (1 - M_{\infty}^2) \phi_{xx} + \phi_{yy} = 0

The equation type depends on the flow speed (Mach number) of the bulk M∞ M_{\infty} . For M∞2<1 M_{\infty}^2 < 1 it is elliptic. For M∞2>1 M_{\infty}^2 > 1 it is hyperbolic. We don't get parabolic solutions because near M∞=1 M_\infty = 1 we are going to get shocks and we are obviously not in linear flow.

It is also possible for the equation type to change within the domain. In particular, look at the TSD (transonic disturbance) equation. The selection parameter is

1−M∞2M∞2>(γ+1)u′u∞ \frac{1 - M_{\infty}^2}{M_{\infty}^2} > (\gamma + 1) \frac{u'}{u_{\infty}}

So for u∞+u′<a u_{\infty} + u' < a we have an elliptic equation, but for u∞+u′>a u_{\infty} + u' > a we have a hyperbolic equation. These are locally evaluated, so we can change types within the same domain. We can imagine the flow around a wing for which u∞>a∞ u_{\infty} > a_{\infty} , then in most of the domain u>a u > a and the problem is hyperbolic, but near the front of the wing where u<a u < a the problem becomes elliptic.

20.4.2 Systems of Equations

Moving away from scalar equations to the world of our Euler equations,

∂Q∂t+∂F∂x=0 \pdv{\vec Q}{t} + \pdv{\vec F}{x} = 0

Q=[ρρue]=[ρme] Q = \begin{bmatrix} \rho \\ \rho u \\ e \end{bmatrix} = \begin{bmatrix} \rho \\ m \\ e \end{bmatrix}

Where we've written ρu \rho u as a separate independent variable m m (we can do that because ρ \rho and u u are independent).

∂Q∂t+A⋅∂F∂x=0A≡∂F∂Q(Flux Jacobian) \pdv{\vec{Q}}{t} + \vec A \cdot \pdv{\vec F}{x} = 0 \qquad \vec A \equiv \pdv{\vec F}{\vec Q} \quad \text{(Flux Jacobian)}

The equation type is given by the characteristics of the flux Jacobian:

  1. Elliptic: The eigenvalues of A \vec A are imaginary.
  2. Parabolic: The eigenvalues are repeated and the eigenvectors are not unique.
  3. Hyperbolic: The eigenvalues are real and the eigenvectors are unique.

20.4.3 Flux Jacobian for Euler Equations

Recall the definition

F‾=[mm2ρ+pmρ(e+p)] \overline{F}= \begin{bmatrix} m \\ \frac{m^2}{\rho} + p \\ \frac{m}{\rho} (e + p) \end{bmatrix}

∂F‾∂Q‾=[∂F1∂Q1∣Q2Q3∂F1∂Q2∣Q1Q3∂F1∂Q1∣Q2Q3∂F2∂Q1∣Q2Q3……] \pdv{\overline{F}}{\overline{Q}} = \begin{bmatrix} \left. \pdv{F_1}{Q_1} \right|_{Q_2 Q_3} & \left. \pdv{F_1}{Q_2} \right|_{Q_1 Q_3} & \left. \pdv{F_1}{Q_1} \right|_{Q_2 Q_3} \\ \left. \pdv{F_2}{Q_1} \right|_{Q_2 Q_3} & \ldots & \ldots \end{bmatrix}

Starting at the first element, and recalling that we take the partial derivative holding Q2 Q_2 and Q3 Q_3 constant,

∂F1∂Q1=∂m∂ρ=0 \pdv{F_1}{Q_1} = \pdv{m}{\rho} = 0

Skipping to the end result,

A‾=[010−m2ρ2+pρ2mρ+pmpe−mρ2(e+p)+mρpρe+pρ+mρpmmρ(1+pe)] \overline{A} = \begin{bmatrix} 0 & 1 & 0 \\ - \frac{m^2}{\rho^2} + p_{\rho} & \frac{2 m}{\rho} + p_m & p_e \\ - \frac{m}{\rho^2} (e + p) + \frac{m}{\rho} p_{\rho} & \frac{e + p}{\rho} + \frac{m}{\rho} p_m & \frac{m}{\rho} (1 + p_e) \end{bmatrix}

Let's talk about how we deal with these partial derivatives of pressure with respect to our conserved quantities

pρ=∂p∂ρ∣me p_{\rho} = \left. \pdv{p}{\rho} \right| _{m e}

From the definition,

e=pγ−1+ρv22 e = \frac{p}{\gamma - 1} + \frac{\rho v^2}{2}

Solve for pressure,

p=(γ−1)(e−ρv22) p = ( \gamma - 1) \left( e - \frac{\rho v^2}{2} \right)

=(γ−1)(e−m22ρ) = ( \gamma - 1 ) \left( e - \frac{m ^2}{2 \rho} \right)

Now we can get our partial derivatives easily

pρ=∂p∂ρ∣me=(γ−1)m2ρ2=12(γ−1)u2 p_{\rho} = \left. \pdv{p}{\rho} \right| _{m e} = (\gamma - 1) \frac{m^2}{\rho ^2} = \frac{1}{2} ( \gamma - 1) u^2

We also define the sound speed

cS2=γpρ c_S ^2 = \frac{\gamma p}{\rho}

Writing our the flux Jacobian in these terms,

A‾=[010γ−32u2−(γ−3)uγ−1−ucs2γ−1+γ−22u3cs2γ−1−(γ−3/2)u2γu] \overline{A} = \begin{bmatrix} 0 & 1 & 0 \\ \frac{\gamma - 3}{2} u^2 & - (\gamma - 3) u & \gamma - 1 \\ - \frac{u c_s ^2}{\gamma - 1} + \frac{\gamma - 2}{2} u^3 & \frac{c_s ^2}{\gamma - 1} - ( \gamma - 3/2) u^2 & \gamma u \end{bmatrix}

If we perform an eigen decomposition of A‾ \overline{A} we find that it has eigenvalues

Λ=(uu−csu+cs) \Lambda = \begin{pmatrix} u & u-c_s & u + c_s \end{pmatrix}

And a complete set of eigenvectors are

X‾=[111uu−csu+csu22u22−csu+cs2γ−1u22+csu+cs2γ−1] \overline{X} = \begin{bmatrix} 1 & 1 & 1 \\ u & u - c_s & u + c_s \\ \frac{u^2}{2} & \frac{u^2}{2} - c_s u + \frac{c_s ^2}{\gamma - 1} & \frac{u^2}{2} + c_s u + \frac{c_s ^2}{\gamma - 1} \end{bmatrix}

This set of eigenvectors are linearly independent and the eigenvalues are all real, which means the Euler equations are a hyperbolic system of equations. In the same way that the scalar advection equation gives wave solutions, the Euler equations will give solutions in which information propagates along characteristics (the eigenvalues) with a wave-like behavior. The eigenvalues we got here are the same characteristics you would get via the method of characteristics in linear fluid dynamics:

  • u+cs u + c_s : Shock wave
  • u u : Contact discontinuity
  • u−cs u - c_s : Rarefaction (expansion) waves

Mathematically, by doing the eigen decomposition we have solved A‾X‾=X‾Λ‾ \overline{A} \overline{X} = \overline{X} \overline{\Lambda} . This solves the complete Euler equations

∂Q∂t+A∂Q∂x=0 \pdv{Q}{t} + A \pdv{Q}{x} = 0

We can right-multiply by the inverse of X \vec X

AXX−1=XΛX−1 \vec A \vec X \vec X^{-1} = \vec X \vec{\Lambda} \vec X ^{-1}

→A=XΛX−1 \rightarrow A = \vec X \vec{\Lambda} \vec X ^{-1}

→∂Q∂t+XΛX−1=0 \rightarrow \pdv{Q}{t} + \vec X \vec \Lambda \vec X^{-1} = 0

Now we make an approximation: Assume that the eigenvector X \vec X is locally constant (in both space in time) such that

∂X∂t=0∂X∂x=0 \pdv{\vec X}{t} = 0 \qquad \pdv{\vec X}{x} = 0

We can re-arrange our X \vec X 's again by pre-multiplying by X−1 \vec X ^{-1}

X−1∂Q∂t+(X−1X)ΛX−1∂Q∂x=0 \vec X^{-1} \pdv{Q}{t} + ( \vec X^{-1} \vec X) \vec \Lambda \vec X ^{-1} \pdv{Q}{x} = 0

→∂∂t(X−1Q)+Λ∂∂x(X−1Q)=0 \rightarrow \pdv{}{t} (\vec X ^{-1} Q) + \vec \Lambda \pdv{}{x} ( \vec X^{-1} \vec Q) = 0

We can define a new variable

W≡X−1Q \vec W \equiv \vec X^{-1} \vec Q

So we have

∂W∂t+Λ∂W∂x=0 \pdv{\vec W}{t} + \vec \Lambda \pdv{\vec W}{x} = 0

Since Λ \Lambda is simply a diagonal 3x3 matrix, we have now de-coupled our system into three independent equations

∂Wi∂t+λi∂Wi∂x=0i=1,2,3 \pdv{W_i}{t} + \lambda_i \pdv{W_i}{x} = 0 \qquad i = 1, 2, 3

While this is not equivalent to the Euler equations, it is a pretty close approximation and it maintains all of the non-linear properties of the Euler equations. Indeed, because of the values λi \lambda_i the system has the same characteristics (u,u+cs,u−cs) (u, u+c_s, u-c_s) as the Euler equations. Now we have a set of 3 advection equations. Therefore, a convenient model for the Euler equations is the linear wave equation

∂u∂t+a∂u∂x=0 \pdv{u}{t} + a \pdv{u}{x} = 0

because it will capture the same wave-like behavior of the Euler equations, but most importantly because we have an analytic solution of the linear advection equation which allows us to check both the stability and the accuracy of our solution.

u(x−at)=const. u ( x - at) = \text{const.}

With periodic boundary conditions, we expect the solution to wrap back around to its initial conditions, allowing us to check for errors we may have introduced in our algorithm.

Another model equation we can use for the Euler equations is Burger's equation, which is a nonlinear advection equation

∂u∂t+u∂u∂x=0 \pdv{u}{t} + u \pdv{u}{x} = 0

Other model equations we can use to model the Navier-Stokes equations is an advection-diffusion equation, which is simply an addition of a diffusion term to the linear wave equation

∂u∂t+a∂u∂x=b∂2u∂x2 \pdv{u}{t} + a \pdv{u}{x} = b \pdv{^2 u }{x^2}

A viscous version of Burger's equation also works as a good model for Navier-Stokes

∂u∂t+u∂u∂x=b∂2u∂x2 \pdv{u}{t} + u \pdv{u}{x} = b \pdv{^2 u}{x^2}

This presents a development path for CFD algorithms. We take an algorithm, apply it to the advection equations and check stability/accuracy, then apply Burger's equation, then add viscous effects, then we may be able to solve the full Euler's equations.