Equations Used in CFD


Contents

NS equations are a set of non-linear partial derivative equations of conservation of mass, momentum, and energy:

ρt+(ρv)=0(ρv)t+(ρvv+pI)=τ+fm(ρE)t+(ρvH)=(kT)+(τv)+fmv\begin{aligned} \frac{\partial \rho}{\partial t}+\nabla\cdot(\rho\mathbf v)&=0\\ \frac{\partial(\rho\mathbf v)}{\partial t}+\nabla\cdot \left( \rho\mathbf v\mathbf v + p\mathbf I \right)&=\nabla\cdot\mathbf\tau+\mathbf f_m\\ \frac{\partial(\rho E)}{\partial t}+\nabla\cdot(\rho\mathbf vH)&=\nabla\cdot(k\mathbf\nabla T)+\nabla\cdot(\mathbf\tau\cdot\mathbf v)+\mathbf f_m\cdot\mathbf v \end{aligned}

Equations of state:

  • p=ρRgTp=\rho R_gT
  • e=cvTe=c_vT
    • Rg=cpcvR_g=c_p-c_v
  • E=e+V2/2E=e+V^2/2
  • H=E+p/ρH=E+p/\rho

These equations are even more complicated in CFD as one needs to add artificial viscosity and equations for modeling turbulences.

Some notes:

  • One should use conservative variables to solve discontinuities correctly
  • Convective terms are more completed than viscous terms

Linearized Euler (1D)

Approximation of ideal gas, no viscosity and heat transfer (τ=q=0\mathbf\tau'=\mathbf q=0, with τ=pI+τ\mathbf\tau=-p\mathbf I+\mathbf\tau'):

t[ρρuρE]+x[ρuρu2+pρuH]=[000]\frac{\partial}{\partial t} \begin{bmatrix} \rho\\ \rho u\\ \rho E \end{bmatrix}+\frac{\partial}{\partial x} \begin{bmatrix} \rho u\\ \rho u^2+p\\ \rho u H \end{bmatrix}= \begin{bmatrix} 0\\0\\0 \end{bmatrix}

Defining conservative variables:

  • u1=ρu_1=\rho
  • u2=ρuu_2=\rho u
  • u3=ρEu_3=\rho E

And with the vector of conservative variable: U=(u1,u2,u3)T\mathbf U=(u_1,u_2,u_3)^T

Ut+Fx=0\frac{\partial\mathbf U}{\partial t}+\frac{\partial\mathbf F}{\partial x}=\mathbf 0

Where F\mathbf F is the vector of conservative flux:

f1=u2f2=u22u1+(γ1)u312(γ1)u22u1f3=u2(γu3u112(γ1)u22u12)\begin{aligned} f_1&=u_2\\ f_2&=\frac{u_2^2}{u_1}+(\gamma-1)u_3-\frac{1}{2}(\gamma-1)\frac{u_2^2}{u_1}\\ f_3&=u_2\left( \gamma\frac{u_3}{u_1}-\frac{1}{2}(\gamma-1)\frac{u_2^2}{u_1^2} \right) \end{aligned}

In quasi-linear form:

Ut+FUUx=0\boxed{\frac{\partial\mathbf U}{\partial t}+\frac{\partial\mathbf F}{\partial\mathbf U}\frac{\partial \mathbf U}{\partial x}=\mathbf 0}

And by defining A=F/UA=\partial\mathbf F/\partial\mathbf U:

A=fiuj=[01012(γ3)(u2u1)2(3γ)u2u1γ1γu2u3u12+(γ1)(u2u3)3γu3u132(γ1)(u2u1)2γu2u1]A=\frac{\partial f_i}{\partial u_j}= \begin{bmatrix} 0 & 1 & 0 \\ \dfrac{1}{2}(\gamma-3)\left( \dfrac{u_2}{u_1} \right)^2 & (3-\gamma)\dfrac{u_2}{u_1} & \gamma-1 \\ -\gamma\dfrac{u_2u_3}{u_1^2}+(\gamma-1)\left( \dfrac{u_2}{u_3} \right)^3 & \gamma\dfrac{u_3}{u_1}-\dfrac{3}{2}(\gamma-1)\left( \dfrac{u_2}{u_1} \right)^2 & \gamma\dfrac{u_2}{u_1} \end{bmatrix}

Using conservative variables and a=γRgTa=\sqrt{\gamma R_g T}:

A=[01012(γ3)u2(3γ)uγ112(γ2)u3a2uγ132γ2u2+a2γ1γu]A= \begin{bmatrix} 0 & 1 & 0 \\ \dfrac{1}{2}(\gamma-3)u^2 & (3-\gamma)u & \gamma-1 \\ \dfrac{1}{2}(\gamma-2)u^3-\dfrac{a^2u}{\gamma-1} & \dfrac{3-2\gamma}{2}u^2+\dfrac{a^2}{\gamma-1} & \gamma u \end{bmatrix}

Perturbations

Supposing U(x,t)=U0+u(x,t)\mathbf U(x,t)=\mathbf U_0+\mathbf u'(x,t) with uU0||\mathbf u'||\ll||\mathbf U_0|| where:

  • U0\mathbf U_0: Uniform and stationary solution
  • u(x,t)\mathbf u'(x,t): General perturbation

Using Taylor expansion with the quasi-linear form, retaining only 1st order components:

ut+A(U0)ux=0\frac{\partial\mathbf u'}{\partial t}+A(\mathbf U_0)\frac{\partial\mathbf u'}{\partial x}=\mathbf 0

This equation is a hyperbolic system if the eigenvalues for A(U0)A(\mathbf U_0) are reals: A(U0)uiR=λiuiRA(\mathbf U_0)\mathbf u_i^R=\lambda_i\mathbf u_i^R

  • λ1=u0\lambda_1=u_0
  • λ2=u0+a0\lambda_2=u_0+a_0
  • λ2=u0a0\lambda_2=u_0-a_0

With a0=γRgT0=γp0/ρ0a_0=\sqrt{\gamma R_gT_0}=\sqrt{\gamma p_0/\rho_0}, and right eigenvectors:

u1R=[1u0u022];u2R=[1u0+a0u022+a02γ1+a0u0];u3R=[1u0a0u022+a02γ1a0u0]\mathbf u_1^R= \begin{bmatrix} 1\\u_0\\\dfrac{u_0^2}{2} \end{bmatrix}\quad;\quad \mathbf u_2^R= \begin{bmatrix} 1\\ u_0+a_0\\ \dfrac{u_0^2}{2}+\dfrac{a_0^2}{\gamma-1}+a_0u_0 \end{bmatrix}\quad;\quad \mathbf u_3^R= \begin{bmatrix} 1\\ u_0-a_0\\ \dfrac{u_0^2}{2}+\dfrac{a_0^2}{\gamma-1}-a_0u_0 \end{bmatrix}

Defining T0=(u1R,u2R,u3R)T_0=(\mathbf u_1^R,\mathbf u_2^R,\mathbf u_3^R) and w\mathbf w:

T01A0T0=Λ0=[λ1000λ2000λ3];w=T01uT_0^{-1}A_0T_0=\Lambda_0= \begin{bmatrix} \lambda_1&0&0\\ 0&\lambda_2&0\\ 0&0&\lambda_3 \end{bmatrix}\quad;\quad \mathbf w=T_0^{-1}\mathbf u'

Adding to the linear equation:

wt+Λ0wx=0\frac{\partial\mathbf w}{\partial t}+\Lambda_0\frac{\partial\mathbf w}{\partial x}=\mathbf 0

where:

  • w1=pa02ρw_1=p'-a_0^2\rho'
  • w2=p+ρ0a0uw_2=p'+\rho_0a_0u'
  • w3=pρ0a0uw_3=p'-\rho_0a_0u'

Boundary conditions:

  • M0<1M_0<1:
    • 2 conditions in x=0x=0 for downstream waves (entropy and acoustic)
    • 1 condition in x=Lx=L for upstream wave
  • M0>1M_0>1:
    • 3 conditions in x=0x=0 for all three waves

Linearized Navier-Stokes

t[ρρuρE]+x[ρuρu2+pρuH]=x[0μxμuux+κTx]\frac{\partial}{\partial t} \begin{bmatrix} \rho\\ \rho u\\ \rho E \end{bmatrix}+\frac{\partial}{\partial x} \begin{bmatrix} \rho u\\ \rho u^2+p\\ \rho u H \end{bmatrix}= \frac{\partial}{\partial x} \begin{bmatrix} 0\\\mu_x\\\mu uu_x+\kappa T_x \end{bmatrix}

Compact form:

Ut+f(U)x=fv(U)x\frac{\partial\mathbf U}{\partial t}+\frac{\partial\mathbf f(\mathbf U)}{\partial x}=\frac{\partial\mathbf f_v(\mathbf U)}{\partial x}
  • f\mathbf f: Convective flux
  • fv\mathbf f_v: Viscous flux

For analytical study, better user primitive variables: V=(ρ,u,p)T\mathbf V=(\rho,u,p)^T. Therefore:

t[ρup]+[uρ00u1/ρ0γpu]x[ρup]=[0(μux)x(γ1)((kTx)x+μux2)]\frac{\partial}{\partial t} \begin{bmatrix} \rho\\u\\p \end{bmatrix}+ \begin{bmatrix} u&\rho&0\\ 0&u&1/\rho\\ 0&\gamma p&u \end{bmatrix} \frac{\partial}{\partial x} \begin{bmatrix} \rho\\u\\p \end{bmatrix}= \begin{bmatrix} 0\\ (\mu u_x)_x\\ (\gamma-1)((kT_x)_x+\mu u_x^2) \end{bmatrix}

By linearizing the system:

t[ρup]+[u0ρ000u01/ρ00γp0u0]A0x[ρup]=[0μuxx(γ1)k0ρ0Rg(pxxp0ρ0ρxx)]\frac{\partial}{\partial t} \begin{bmatrix} \rho'\\u'\\p' \end{bmatrix}+ \underbrace{ \begin{bmatrix} u_0&\rho_0&0\\ 0&u_0&1/\rho_0\\ 0&\gamma p_0&u_0 \end{bmatrix}}_{A_0} \frac{\partial}{\partial x} \begin{bmatrix} \rho'\\u'\\p' \end{bmatrix}= \begin{bmatrix} 0\\ \mu u_{xx}'\\ \dfrac{(\gamma-1)k_0}{\rho_0R_g} \left( p'_{xx}-\dfrac{p_0}{\rho_0}\rho'_{xx} \right) \end{bmatrix}

Eigenvalues for A0A_0 are the same, but the eigenvectors are now the columns of Q0Q_0:

Q0=[1ρ02a0ρ02a0012120ρ0a02ρ0a0a];Q01=[101a02011ρ0a0011ρ0a0]Q_0= \begin{bmatrix} 1&\dfrac{\rho_0}{2 a_0}&-\dfrac{\rho_0}{2a_0}\\ 0&\dfrac{1}{2}&\dfrac{1}{2}\\ 0&\dfrac{\rho_0a_0}{2}&-\dfrac{\rho_0a_0}{a} \end{bmatrix}\quad;\quad Q_0^{-1}= \begin{bmatrix} 1&0&-\dfrac{1}{a_0^2}\\ 0&1&\dfrac{1}{\rho_0a_0}\\ 0&1&-\dfrac{1}{\rho_0a_0} \end{bmatrix}

By multiplying Q01Q_0^{-1} to the linearized equation:

t[w1w2w3]+[u0000u0+a0000u0a0]x[w1w2w3]=B2x2[w1w2w3]\frac{\partial}{\partial t} \begin{bmatrix} w_1\\w_2\\w_3 \end{bmatrix}+ \begin{bmatrix} u_0&0&0\\ 0&u_0+a_0&0\\ 0&0&u_0-a_0 \end{bmatrix} \frac{\partial}{\partial x} \begin{bmatrix} w_1\\w_2\\w_3 \end{bmatrix}= B\frac{\partial^2}{\partial x^2} \begin{bmatrix} w_1\\w_2\\w_3 \end{bmatrix}

Stripping dimensions with the base field values:

  • u=u/a0\overline u=u'/a_0
  • p=p/(ρ0a02)\overline p=p'/(\rho_0a_0^2)
  • t=t/(Lc/a0)\overline t=t/(L_c/a_0)
  • x=x/Lc\overline x=x/L_c

The equation is therefore:

t[w1w2w3]+[u0000u0+1000u01]x[w1w2w3]=B2x2[w1w2w3]B=1Re[00001/21/201/21/2]+1RePr[1γ12γ121γ2γ221γ22γ2]\begin{gathered} \frac{\partial}{\partial t} \begin{bmatrix} \overline w_1\\\overline w_2\\\overline w_3 \end{bmatrix}+ \begin{bmatrix} \overline u_0&0&0\\ 0&\overline u_0+1&0\\ 0&0&\overline u_0-1 \end{bmatrix} \frac{\partial}{\partial x} \begin{bmatrix} \overline w_1\\\overline w_2\\\overline w_3 \end{bmatrix}= \overline B\frac{\partial^2}{\partial x^2} \begin{bmatrix} \overline w'_1\\\overline w'_2\\\overline w'_3 \end{bmatrix}\\ B=\frac{1}{\mathrm{Re}} \begin{bmatrix} 0&0&0\\ 0&1/2&-1/2&\\ 0&-1/2&1/2 \end{bmatrix}+\frac{1}{\mathrm{Re\,Pr}} \begin{bmatrix} 1&\dfrac{\gamma-1}{2}&\dfrac{\gamma-1}{2}\\ 1&\dfrac{\gamma}{2}&\dfrac{\gamma-2}{2}\\ 1&\dfrac{\gamma-2}{2}&\dfrac{\gamma}{2} \end{bmatrix} \end{gathered}
  • Re=ρ0u0Lc/μ0\mathrm{Re}=\rho_0u_0L_c/\mu_0
  • Pr=cpμ0/κ0\mathrm{Pr}=c_p\mu_0/\kappa_0
  • w1=pρ\overline w_1=\overline p'-\overline\rho'
  • w1=p+u\overline w_1=\overline p'+\overline u'
  • w1=pu\overline w_1=\overline p'-\overline u'

Some conclusions:

  • If κ0=0\kappa_0=0: Entropy equation decouples, and ss' is conserved along the path lines
  • If Re1\mathrm{Re}\gg1: Waves decouple. Perturbed characteristic variables are conserved along the path lines, slightly damped due to the viscosity

Entropy equation:

ρT ⁣Ds ⁣Dt=(κT)+Φv+Q\rho T\frac{\DD s}{\DD t}=\mathbf\nabla\cdot(\kappa\mathbf\nabla T)+\mathbf\Phi_v+Q

Linearized entropy equation:

ρ0T0(ts+u0xs)=κ0x2T+Φv\rho_0T_0(\partial_ts'+u_0\partial_xs')=\kappa_0\partial_x^2T'+\cancel{\mathbf\Phi'_v}

Vorticity equation (ω=×u\mathbf\omega=\mathbf\nabla\times\mathbf u) for isentropic and ideal flow, without forces from masses:

 ⁣D(ω/ρ) ⁣Dt=ωρV\frac{\DD (\mathbf\omega/\rho)}{\DD t}=\frac{\mathbf\omega}{\rho}\cdot\mathbf\nabla\mathbf V

Acoustics Equations

Using the following decomposition in the Euler’s equation:

u=uρ=ρ0+ρp=p0+p\begin{aligned} u&=u'\\ \rho&=\rho_0+\rho'\\ p&=p_0+p' \end{aligned}

We can obtain:

0=t(ρ0+ρ)+x[(ρ0+ρ)u]0=t[(ρ0+ρ)u]+x[(ρ0+ρ)uu+p0+p]\begin{aligned} &0=\partial_t(\rho_0+\rho')+\partial_x[(\rho_0+\rho')u']\\ &0=\partial_t[(\rho_0+\rho')u']+\partial_x[(\rho_0+\rho')u'u'+p_0+p'] \end{aligned}

As ρ0\rho_0 and p0p_0 (with u0=0u_0=0) also verifies Euler’s equation:

tρ0=0;xp0=0\partial_t\rho_0=0\quad;\quad\partial_xp_0=0

Using these values we can simplify the equation above:

0=tρ+x[(ρ0+ρ)u]0=ρ0tu+t(ρu)+x[(ρ0+ρ)uu]+xp\begin{aligned} &0=\partial_t\rho'+\partial_x[(\rho_0+\rho')u']\\ &0=\rho_0\partial_tu'+\partial_t(\rho'u')+\partial_x[(\rho_0+\rho')u'u']+\partial_xp' \end{aligned}

Looking for linear solutions, neglecting quadratic terms of the perturbations:

0=tρ+x(ρ0u)0=ρ0tu+xρ\begin{aligned} 0&=\partial_t\rho'+\partial_x(\rho_0u')\\ 0&=\rho_0\partial_tu'+\partial_x\rho' \end{aligned}

The relationship between pp' and ρ\rho' is:

p=p0+pp(ρ)p(ρ0)+pρρ0ρ}p=pρs0,ρ0ρ=a02ρ\begin{rcases} p&=p_0+p'\\ p(\rho)&\approx p(\rho_0)+ \left. \dfrac{\partial p}{\partial\rho} \right|_{\rho_0}\rho' \end{rcases}p'=\left. \frac{\partial p}{\partial \rho} \right|_{s_0,\rho_0}\rho'=a_0^2\rho'

Using this property results in:

0=1a02tp+ρ0xu0=ρ0tu+xp\begin{aligned} 0&=\frac{1}{a_0^2}\partial_tp'+\rho_0\partial_xu'\\ 0&=\rho_0\partial_t u'+\partial_xp' \end{aligned}

Which is a system of two linear equations with two unknowns.

Deriving the first equation with respect to time and the second with respect to xx:

1a02ttp+ρ0tx(u)=0ρ0xtu+xxp=0}1a02ttp=xxp\begin{rcases} \dfrac{1}{a_0^2}\partial_{tt}p'+\rho_0\partial_{tx}(u')&=0\\ \rho_0\partial_{xt}u'+\partial_{xx}p'&=0 \end{rcases} \frac{1}{a_0^2}\partial_{tt}p'=\partial_{xx}p'

Which is the equation of the perturbation in the pressure field. In 2D:

1a02ttp=2p;2=xx+yy\frac{1}{a_0^2}\partial_{tt}p'=\nabla^2p'\quad;\quad\nabla^2=\partial_{xx}+\partial_{yy}

Generalized Equation

Considering the acoustic equation with base flux speed u00u_0\ne0, the linearized Euler’s equations in 1D:

0=1a02pt+ρ0ux+u0a02px0=ut+1ρ0px+uxux\begin{aligned} 0&=\frac{1}{a_0^2}\frac{\partial p'}{\partial t}+\rho_0\frac{\partial u'}{\partial x}+\frac{u_0}{a_0^2}\frac{\partial p'}{\partial x}\\ 0&=\frac{\partial u'}{\partial t}+\frac{1}{\rho_0}\frac{\partial p'}{\partial x}+u_x\frac{\partial u'}{\partial x} \end{aligned}

Combining both:

1a02 ⁣D02p ⁣Dt2=2p; ⁣D0 ⁣Dt=t+u01a02(2pt2+2u02ptx+u022px2)=2px2\begin{gathered} \frac{1}{a_0^2}\frac{\DD_0^2p'}{\DD t^2}=\nabla^2p'\quad;\quad\frac{\DD_0}{\DD t}=\frac{\partial}{\partial t}+u_0\nabla\\ \frac{1}{a_0^2}\left( \frac{\partial^2p'}{\partial t^2}+2u_0\frac{\partial^2p'}{\partial t\partial x}+u_0^2\frac{\partial^2p'}{\partial x^2} \right)=\frac{\partial^2p'}{\partial x^2} \end{gathered}

which transforms into an equation with a 2nd temporal derivative with respect to time into two of 1st order:

t(p1p2)=(01(a02u02)(2/x2)2u0(/x))(p1p2)\frac{\partial}{\partial t} \begin{pmatrix} p_1\\p_2 \end{pmatrix}= \begin{pmatrix} 0&1\\ (a_0^2-u_0^2)(\partial^2/\partial x^2)&-2u_0(\partial/\partial x) \end{pmatrix} \begin{pmatrix} p_1\\p_2 \end{pmatrix}

where p1=pp_1=p' and p2=p˙p_2=\dot p'

Other Equations

  • Linear waves: tu+axu=0\partial_tu+a\partial_xu=0
  • Viscous lineal waves: tu+axu=νx2u\partial_tu+a\partial_xu=\nu\partial_x^2u
  • Burgers: tu+uxu=0\partial_tu+u\partial_xu=0
  • Viscous burgers (Convection-Diffusion): tu+uxu=νx2u\partial_tu+u\partial_xu=\nu\partial_x^2u
  • Heat: tu=αx2u\partial_tu=\alpha\partial_x^2u
  • Stationary heat (Poisson): x2u=0\partial_x^2u=0
  • Linear wave with source: tu+axu=αu\partial_tu+a\partial_xu=\alpha u
  • 2D linear waves: tu+Axu+Byu=0\partial_tu+A\partial_xu+B\partial_yu=0