By ignoring viscosity and heat conduction, which are less important in numerical computations, the resulting equations are the Euler equations which are based in Navier-Stokes equations in an ideal flow.
As numerical schemes should behave correctly with small perturbations, instead of working with a generic flow U \mathbf U U u ′ \mathbf u' u ′ U 0 \mathbf U_0 U 0
U = U 0 + u ′ \mathbf U=\mathbf U_0+\mathbf u' U = U 0 + u ′ Then, if the scheme works for u ′ \mathbf u' u ′ U \mathbf U U
The 1D linearized Euler equations can be de split into 3 equations of characteristic variables with the following form:
∂ u ∂ t + c ∂ u ∂ x = 0 \frac{\partial u}{\partial t}+c\frac{\partial u}{\partial x}=0 ∂ t ∂ u + c ∂ x ∂ u = 0 Contents
Analysis
Linear wave equation with constant c > 0 c>0 c > 0
∂ u ∂ t + c ∂ u ∂ x = 0 \frac{\partial u}{\partial t}+c\frac{\partial u}{\partial x}=0 ∂ t ∂ u + c ∂ x ∂ u = 0 Boundary conditions:
Non-periodic (general): u 0 = u 0 ( t ) u_0=u_0(t) u 0 = u 0 ( t )
Periodic: u ( 0 ) = u ( L ) u(0)=u(L) u ( 0 ) = u ( L )
For periodic boundary conditions, solutions are plain waves.1
u = ∑ k = − ∞ ∞ u ^ k ( t ) e i k x u=\sum_{k=-\infty}^\infty\hat u_k(t)e^{\iu kx} u = k = − ∞ ∑ ∞ u ^ k ( t ) e i k x
i = − 1 \iu=\sqrt{-1} i = − 1 k = 2 π / l k=2\pi/l k = 2 π / l l l l
Hypothesis:
Variable separation: Solution can be split into the product of a spatial and a temporal part
Linearity: Superposition, u = u ^ k ( t ) e i k x + conditions u=\hat u_k(t)e^{\iu kx}+\text{conditions} u = u ^ k ( t ) e i k x + conditions
Note that u ^ k ( t ) e i k x \hat u_k(t)e^{\iu kx} u ^ k ( t ) e i k x u ^ − k ( t ) e − i k x \hat u_{-k}(t)e^{-\iu kx} u ^ − k ( t ) e − i k x
e i k x d u ^ k ( t ) d t = − i c k u ^ k ( t ) e i k x \cancel{e^{\iu kx}}\frac{\dd\hat u_k(t)}{\dd t}=-\iu ck\hat u_k(t)\cancel{e^{\iu kx}} e i k x d t d u ^ k ( t ) = − i c k u ^ k ( t ) e i k x Which results into an initial value problem:
d u ^ k ( t ) d t = λ u ^ k ( t ) ; λ = − i c k ; \boxed{\frac{\dd \hat u_k(t)}{\dd t}=\lambda\hat u_k(t)\quad;\quad\lambda=-\iu ck\quad;\quad} d t d u ^ k ( t ) = λ u ^ k ( t ) ; λ = − i c k ;
λ \lambda λ u ^ k \hat u_k u ^ k u ^ k ( t = 0 ) = u ^ k , 0 \hat u_k(t=0)=\hat u_{k,0} u ^ k ( t = 0 ) = u ^ k , 0
Integrations
The equation can be easily integrated, obtaining:
u ^ k ( t ) = u ^ k ( t = 0 ) e λ t = u ^ k ( t = 0 ) e − i c k t \hat u_k(t)=\hat u_k(t=0)e^{\lambda t}=\hat u_k(t=0)e^{-\iu ckt} u ^ k ( t ) = u ^ k ( t = 0 ) e λ t = u ^ k ( t = 0 ) e − i c k t The analytical amplification factor (gain) G G G
G = u ( t + Δ t ) u ( t ) = ∣ u ( t + Δ t ) u ( t ) ∣ e i Φ a \boxed{G=\frac{u(t+\Delta t)}{u(t)}= \left| \frac{u(t+\Delta t)}{u(t)} \right|e^{\iu\Phi_a}\\} G = u ( t ) u ( t + Δ t ) = u ( t ) u ( t + Δ t ) e i Φ a
u = u ^ k u=\hat u_k u = u ^ k Φ a \Phi_a Φ a
When using numerical integration for the time:
G = u n + 1 u n = ∣ u n + 1 u n ∣ e i Φ n \boxed{G=\frac{u^{n+1}}{u_n}=\left| \frac{u^{n+1}}{u^n} \right|e^{\iu\Phi_n}} G = u n u n + 1 = u n u n + 1 e i Φ n Using solutions like u ^ k ( t ) = u ^ k ( t = 0 ) e λ t \hat u_k(t)=\hat u_k(t=0)e^{\lambda t} u ^ k ( t ) = u ^ k ( t = 0 ) e λ t analytical integration :
G = e λ ( t + Δ t ) e λ t = e λ Δ t \boxed{G=\frac{e^{\lambda(t+\Delta t)}}{e^{\lambda t}}=e^{\lambda\Delta t}} G = e λ t e λ ( t + Δ t ) = e λ Δ t
In general, λ = λ R + i λ I \lambda=\lambda_R+\iu\lambda_I λ = λ R + i λ I G = e λ R Δ t e i λ I δ t G=e^{\lambda_R\Delta t}e^{\iu\lambda_I\delta t} G = e λ R Δ t e i λ I δ t
The modulus is therefore ∣ G ∣ = e λ R Δ t |G|=e^{\lambda_R\Delta t} ∣ G ∣ = e λ R Δ t
Integration is stable when ∣ G ∣ ≤ 1 |G|\le1 ∣ G ∣ ≤ 1
Computational Variables
k ^ = k Δ x \hat k=k\Delta x k ^ = k Δ x
k ^ = 2 π Δ x / l = 2 π / N points per wave lengh \hat k=2\pi\Delta x/l=2\pi/N_{\text{points per wave lengh}} k ^ = 2 π Δ x / l = 2 π / N points per wave lengh k ^ → 0 \hat k\rightarrow0 k ^ → 0 k ^ → π \hat k\rightarrow\pi k ^ → π G > 1 G>1 G > 1 Refining effect: Keeping k k k Δ x ∝ k ^ \Delta x\propto\hat k Δ x ∝ k ^
C F L = c Δ t / Δ x \mathrm{CFL}=c\Delta t/\Delta x CFL = c Δ t /Δ x
λ Δ t = − i c k Δ t = − i C F L k ^ \boxed{\lambda\Delta t=-\iu ck\Delta t=-\iu\mathrm{CFL}\hat k} λ Δ t = − i c k Δ t = − i CFL k ^ Dispersion Relation
In some problems, perturbations depend on time harmonically : u ^ k ( t ) = u ^ k 0 e i ω t \hat u_k(t)=\hat u_k^0e^{\iu\omega t} u ^ k ( t ) = u ^ k 0 e i ω t
ω = 2 π / T \omega=2\pi/T ω = 2 π / T T T T
Introducing this solution into PVI:
d u ^ k ( t ) d t = λ u ^ k ( t ) → i ω u k 0 e i ω t = λ u k 0 e i ω t → i ω − λ = 0 \frac{\dd\hat u_k(t)}{\dd t}=\lambda\hat u_k(t)\quad \rightarrow \quad \iu\omega\cancel{u_k^0e^{\iu\omega t}}=\lambda\cancel{u_k^0e^{\iu\omega t}} \quad\rightarrow\quad \iu\omega-\lambda=0 d t d u ^ k ( t ) = λ u ^ k ( t ) → i ω u k 0 e i ω t = λ u k 0 e i ω t → i ω − λ = 0 With λ = − i c k \lambda=-\iu ck λ = − i c k dispersion relation is:
i ω + i c k = 0 → ω + c k = 0 \iu\omega+\iu ck=0\quad \rightarrow \quad\omega+ck=0 i ω + i c k = 0 → ω + c k = 0
For a harmonic perturbation ω \omega ω c c c k = − ω / c k=-\omega/c k = − ω / c
Linear Wave Equation with Dissipation
The equation:
∂ u ∂ t + c ∂ u ∂ x = ν ∂ 2 u ∂ x 2 \frac{\partial u}{\partial t}+c\frac{\partial u}{\partial x}=\nu\frac{\partial^2u}{\partial x^2} ∂ t ∂ u + c ∂ x ∂ u = ν ∂ x 2 ∂ 2 u
Initial conditions: u ( x , 0 ) = u I ( x ) u(x,0)=u_I(x) u ( x , 0 ) = u I ( x )
Contour conditions:
u ( 0 , t ) = u 0 ( t ) u(0,t)=u_0(t) u ( 0 , t ) = u 0 ( t ) x x x u ( L , t ) = u L ( t ) u(L,t)=u_L(t) u ( L , t ) = u L ( t ) t t t
Supposing solutions like u ( x , t ) = u ^ k ( t ) e i k x u(x,t)=\hat u_k(t)e^{\iu kx} u ( x , t ) = u ^ k ( t ) e i k x
d u ^ k d t = λ u ^ k ( t ) \frac{\dd\hat u_k}{\dd t}=\lambda\hat u_k(t) d t d u ^ k = λ u ^ k ( t ) Where λ = − i c k − ν k 2 \lambda=-\iu ck-\nu k^2 λ = − i c k − ν k 2 λ R = − ν k 2 \lambda_R=-\nu k^2 λ R = − ν k 2 ∣ G ∣ < 1 |G|<1 ∣ G ∣ < 1
Integrating the equation gives:
u ^ k ( t ) = u ^ k ( t = 0 ) e λ t = u ^ k ( t = 0 ) e ( − i c k − ν k 2 ) t \hat u_k(t)=\hat u_k(t=0)e^{\lambda t}=\hat u_k(t=0)e^{(-\iu ck-\nu k^2)t} u ^ k ( t ) = u ^ k ( t = 0 ) e λ t = u ^ k ( t = 0 ) e ( − i c k − ν k 2 ) t As u ( x , t ) = u ^ k ( t ) e i k x u(x,t)=\hat u_k(t)e^{\iu kx} u ( x , t ) = u ^ k ( t ) e i k x u ^ k ( t ) = u ^ k ( t = 0 ) e λ t \hat u_k(t)=\hat u_k(t=0)e^{\lambda t} u ^ k ( t ) = u ^ k ( t = 0 ) e λ t
u ( x , t ) = u ^ k ( t = 0 ) e i k ( x − i λ t / k ) u(x,t)=\hat u_k(t=0)e^{\iu k(x-\iu\lambda t/k)} u ( x , t ) = u ^ k ( t = 0 ) e i k ( x − i λ t / k ) Separating by parts:
u ( x , t ) = u ^ k ( t = 0 ) e i k ( x + λ I t / k ) e λ R t u(x,t)=\hat u_k(t=0)e^{\iu k(x+\lambda_I t/k)}e^{\lambda_Rt} u ( x , t ) = u ^ k ( t = 0 ) e i k ( x + λ I t / k ) e λ R t In this case, the dispersion relation is:
i ω + i c k + ν k 2 = 0 \iu\omega+\iu ck+\nu k^2=0 i ω + i c k + ν k 2 = 0 Numerical Resolution
Method of lines for a generic internal node j j j
d u j d t = − c ∂ u ∂ x ∣ j \frac{\dd u_j}{\dd t}=-c \left. \frac{\partial u}{\partial x} \right|_j d t d u j = − c ∂ x ∂ u j Central Discretization without Viscosity
Given the central discretization:
d u j d t = − c 2 Δ x ( u j + 1 − u j − 1 ) \frac{\dd u_j}{\dd t}=-\frac{c}{2\Delta x}(u_{j+1}-u_{j-1}) d t d u j = − 2Δ x c ( u j + 1 − u j − 1 ) Using Fourier:
u j = u ^ k ( t ) e i k x j u_j=\hat u_k(t)e^{\iu kx_j} u j = u ^ k ( t ) e i k x j u j − 1 = u ^ k ( t ) e i k x j − 1 = u ^ k ( t ) e i k x j − i k Δ x = u ^ k ( t ) e i k x j e − i k ^ u_{j-1}=\hat u_k(t)e^{\iu kx_{j-1}}=\hat u_k(t)e^{\iu kx_j-\iu k\Delta x}=\hat u_k(t)e^{\iu kx_j}e^{-\iu\hat k} u j − 1 = u ^ k ( t ) e i k x j − 1 = u ^ k ( t ) e i k x j − i k Δ x = u ^ k ( t ) e i k x j e − i k ^ u j + 1 = u ^ k ( t ) e i k x j + 1 = u ^ k ( t ) e i k x j + i k Δ x = u ^ k ( t ) e i k x j e i k ^ u_{j+1}=\hat u_k(t)e^{\iu kx_{j+1}}=\hat u_k(t)e^{\iu kx_j+\iu k\Delta x}=\hat u_k(t)e^{\iu kx_j}e^{\iu\hat k} u j + 1 = u ^ k ( t ) e i k x j + 1 = u ^ k ( t ) e i k x j + i k Δ x = u ^ k ( t ) e i k x j e i k ^
Injecting into the previous equation and cancel out e i k x j e^{\iu kx_j} e i k x j
d u ^ k d t = − c 2 Δ x ( e i k ^ − e − i k ^ ) u ^ k \frac{\dd \hat u_k}{\dd t}=-\frac{c}{2\Delta x} \left( e^{\iu\hat k}-e^{-\iu\hat k} \right)\hat u_k d t d u ^ k = − 2Δ x c ( e i k ^ − e − i k ^ ) u ^ k As e i k ^ − e i k ^ = 2 i sin ( k ^ ) e^{\iu\hat k}-e^{\iu\hat k}=2\iu\sin(\hat k) e i k ^ − e i k ^ = 2 i sin ( k ^ )
d u ^ k d t = − i c Δ x sin ( k ^ ) u ^ k = λ u ^ k \frac{\dd\hat u_k}{\dd t}=-\iu\frac{c}{\Delta x}\sin(\hat k)\hat u_k=\lambda\hat u_k d t d u ^ k = − i Δ x c sin ( k ^ ) u ^ k = λ u ^ k The new eigenvalue λ = − i c sin ( k ^ ) / Δ x \lambda=-\iu c\sin(\hat k)/\Delta x λ = − i c sin ( k ^ ) /Δ x k ^ ≪ 1 \hat k\ll1 k ^ ≪ 1 k ^ = π \hat k=\pi k ^ = π
The phase velocity v f = − λ I / k v_f=-\lambda_I/k v f = − λ I / k
v f = c sin ( k ^ ) k ^ { c if k ^ ≪ 1 0 if k ^ = π v_f=\frac{c\sin(\hat k)}{\hat k}
\begin{cases}
c&\text{if }\hat k\ll1\\
0&\text{if }\hat k=\pi
\end{cases} v f = k ^ c sin ( k ^ ) { c 0 if k ^ ≪ 1 if k ^ = π High frequency waves are therefore, wrongly computed.
Central Discretization with Artificial Viscosity
We will be using artificial viscosity to damp high frequency waves, keeping the low frequency waves intact:
∂ u ∂ t + c ∂ u ∂ x = μ x ∣ c ∣ Δ x ∂ 2 u ∂ x 2 \frac{\partial u}{\partial t}+c\frac{\partial u}{\partial x}=\mu_x|c|\Delta x\frac{\partial^2 u}{\partial x^2} ∂ t ∂ u + c ∂ x ∂ u = μ x ∣ c ∣Δ x ∂ x 2 ∂ 2 u The above equations results in, after using method of lines:
d u ^ k d t = [ − i c Δ x sin ( k ^ ) − 2 μ 2 c Δ x ( 1 − cos ( k ^ ) ) ] u ^ k = λ u ^ k \frac{\dd\hat u_k}{\dd t}=
\left[ -\iu\frac{c}{\Delta x}\sin(\hat k)-\frac{2\mu_2c}{\Delta x}(1-\cos(\hat k)) \right]\hat u_k=\lambda\hat u_k d t d u ^ k = [ − i Δ x c sin ( k ^ ) − Δ x 2 μ 2 c ( 1 − cos ( k ^ )) ] u ^ k = λ u ^ k That is:
Δ t λ = − i C F L sin ( k ^ ) − 2 μ 2 C F L ( 1 − cos ( k ^ ) ) \boxed{\Delta t\lambda=-\iu\mathrm{CFL}\sin(\hat k)-2\mu_2\mathrm{CFL}(1-\cos(\hat k))} Δ t λ = − i CFL sin ( k ^ ) − 2 μ 2 CFL ( 1 − cos ( k ^ )) Where the maximum damping is at k ^ = π \hat k=\pi k ^ = π