Appunti Simulation and modelling of turbulent flows

Space discretization

Wavenumber: k_m = \frac{2 \pi k}{\Delta} is a measure of the spatial frequency of a wave.

Wavelength: \lambda_m = \frac{2 \pi}{k_m}

In Fourier space a convolution integral becomes a simple multiplication:

(f*g)_\Delta=\int f(t)g(x-t)dt=\frac{1}{\Delta}(f^* g)_k = \frac{2 \pi}{\Delta} \hat f_k \hat g_k

Parseval's theorem:

\int | f(x) |^2 dx = \frac{2 \pi}{\Delta} \sum_{k=-\infty}^{+\infty} | \hat{f}_k |^2

Derivative property:

• The Fourier coefficients of g = df/dx is:

\frac{df}{dx} = \frac{d}{dx} \sum_{k=-\infty}^{\infty} f_ke^{ikx} = \sum_{k=-\infty}^{\infty} ik f_ke^{ikx}

\hat{g}_k = \left( \frac{df}{dx} \right)_k = jk \hat{f}_k

• The Fourier coefficients of g = d^2 f/dx^2 is:

\frac{d^2f}{dx^2}=\frac{d}{dx}\sum_{k=-\infty}^{\infty}ikf_ke^{ikx}=\sum_{k=-\infty}^{\infty}(ik)^2f_ke^{ikx}

\hat{g}_k=\left(\frac{d^2f}{dx^2}\right)_k=-k^2\hat{f}_k

We now consider a finite number of harmonics: the truncated Fourier series is defined as:

S_N(f(x))=\frac{1}{\Delta}\sum_{k=-N}^{N}\hat{f}_ke^{ikx}

where the number of harmonics is truncated to N+1. This truncation will introduce truncation error as a function of N: err_1=\left| f(x)-S_N(f(x))\right|

There is a second source of error: the aliasing, where the integral in the definition of Fourier coefficient:

f_{m}=\frac{1}{\Delta}\int f(x)e^{-ijx_k}dx

is computed by using a finite number of points thus introducing the discrete Fourier transform. When using N discretization points it is impossible to distinguish an harmonic with wavenumber k with another harmonic with a higher wavenumber (s). In fraction flows the largest waves number of the scale \left(k_m\right)_{max}=\frac{2\pi}{\Delta}

It's important to find how many sampling points M let us avoid aliasing error when analyzing a potential with harmonics up to e^{-i2\pi}

We need Fourier coefficients for all harmonics below N/2 or M=2N. We discovered that we need at least: M=N.

Finite difference method

We now consider 1D space N+1 points where the velocity u(x) in the discrete 1D-space assumes:

u(x_i) where x_i=(i-1)\Delta x, i=0, ... , N+1 and \Delta x=\frac{L}{N}

We want algebraic version of differential operators: the continuum derivative is defined as:

\frac{d\varphi}{dx}=\lim_{h \to 0}\frac{U(x+h)-U(x)}{h}

d_x = lim_h→0 ᵤ(x+h)−ᵤ(x) / h

by neglecting the limit we achieve,

(dᵤ)x ≃ ᵤ(x+Δx) − ᵤ(x) / Δx = ᵤ_i+1 − ᵤ_i / Δx

It is important now to characterize the quality of the approximation by introducing the Taylor Series in the expression that reproduces the exact differential operator while x→0. That is why we consider the Taylor series:

ᵤ_i+1 = ᵤ_i+(dᵤ)_xΔx + (d^2ᵤ / dx²)_x (Δx² / 2) + o(Δx³) (1)

and truncate it at the first order to obtain:

(dᵤ)_x = ᵤ_i+1−ᵤ_i / Δx + o(Δx) . . . . first order accurate

We can do the same thing with backward Taylor:

(dᵤ)_x = ᵤ_i−ᵤ_i−1 / Δx + o(Δx) . . . . first order accurate

These are FFD and BFD, we can subtract them (2)-(1) to achieve CFD:

ᵤ_i+1 − ᵤ_i−1 = 2(dᵤ)_dxΔx + (d^3ᵤ / 3) (Δx³) + o(Δx⁴)

(d^oᵤ / dx¹)_x = ᵤ_i+1−ᵤ_i−1 / 2Δx + o(Δx²) . . . . second order accurate

We use the same approach to find the second order derivative:

ᵤ_i+1 + ᵤ_i-1 − 2ᵤ = d^2ᵤ / dx² (Δx²) + o(Δx⁴)

(d^2ᵤ / dx²)_x = ᵤ_i+1 − 2ᵤ_i + ᵤ_i-1 / Δx² + o(Δx²) . . . . second order accurate

Modified WNB

An important property of numerical techniques in CFD is the spectral accuracy. A way to address it is given by the analysis of the modified wnb.

The spectral accuracy of a method represented as a function of the length scale.

We can use the exact expression of the roundoff coeffcient definitions as terms of reference to evaluate the spectral accuracy of numerical methods.

We know:

dᵤ / dx = -jkᵤ_k

d^2ᵤ / dx² = -k²ᵤ_k

1^o order derivative CFD

(d^oᵤ / dx)_x = ᵤ_i+1 - ᵤ_i-1 / 2Δx

Now substitute with Fourier series:

(d^o / dx)^c = 1 / 2Δx ( ∑_k=−n/2^n/2 ᵤ_k e^jkxi+1 − ∑ ᵤ_k e^jk(xi-Δx) )

= 1 / 2Δx ( ∑ ᵤ_ke^jᴋx_i )

This is an additional amplitude error that could be:

• ARTIFICAL DISSIPATION
• ARTIFICAL AMPLIFICATION

It's important to avoid both of them because:

• _{Im ω} < 0: as t → ∞, û _k → 0
• _{Im ω} > 0: as t → ∞, û _k → ∞

Art. dis. causes unphysical damping of high wn interonicsArt. amp. causes numerical instabilities because increases amplitude

We also know that Im part of FFD and BFD have opposite sign.To avoid artificial amplification we use FFD or BFD depending on convective velocity sign:

• c : > 0 ⇝ BFD
• c : < 0 ⇝ FFD

This is the UPWIND SCHEME

DIFFUSION EQUATION

∂x∂t = v ∂²xu(x,t) = u(λn,t)u(x,0) = u(x)

by applying Fourier :û_(x,t) = Σ_k e^jxw_k-νk²te^-iwt

because:

dû_k^t_dt=-vk^²iû_k û_k(t) - û_ke^-v²k²t

The diffusion process will smooth all irregularities of the initial condition: the duamping increases by increasing the wvn kn like a quadratic function:

In physical space becomes:u(x,t) = Σ_k e^{jxw_k-v²k²t}e^-ixw_k

Let's now solve it by applying 1° order CFD with M=N:d^tû_k = -ν(k²n²)u_kn̂

for k=-n/2..n/2

The solution of the system is:ûk(t) = ûe^-v²kn²dt_t

for k=n+/2...n/2

The damping effect is underestimated at high wvn wavenumbers:If we use two times height, the solution is drastically worse.u_x=1, on which we only to dec top, while it should be _im, the harmonic is in the first body.As a consequence as t→∞, all other harmonics will be damped with except for k=0 (with u_m)

The linearized equation is:

ξ_i+1^m+1 = ξ_i+1^m - Ñ^mΔt + σ^m+1 -

ξ_i+1^m = Û^m - Û^m-1 - Û^m - Ñ^mΔt + σ^m - Ñ^mΔt =

= Ũ^m Û^m+1 + { Ñ^' - Ñ^m }Δt + σ^m**

We do the same trick as before but with Ñ^m:

Ñ^m = Ñ̆^m + ( dÑ^m / dṪ ) ( Û̆^m - Û̆^m-1 ) + O( ( Û̆^m - Û̆^m-1 ) )²

Ñ^m - Ñ̆^m = ( dÑ^m / dṪ ) ξ^m Δt + O( ( ξ^m-1 ξ ) )²

by substituting:

ξ^m+1 = ξ^m + ( dÑ^m / dṪ ) ξ^m Δt + O ( ( ξ^m-1 ξ ) )² Δt + σ^m

= [1 + ( dÑ^m / dṪ ) Δt ] ξ^m + O( ( ξ^m-1 ξ ) )² Δt + σ^m

amplification factor = α

Contrary to linear problems this α cannot be known a priori, so it needs to be computed at each iteration to continously adapt Δt to verify | α | ≤ 1.

CN FOR NONLINEAR PROBLEMS

The algorithm is:

Û^m+1 = Û̆^m + ( Ñ^m - Ñ̂^m / 2 ) Δt / 2

Ñ˘^m is unknown and it’s our problem, meaning that explicit schemes are no use for nonlinear problems. A method to solve it then is the:

PREDICTOR CORRECTOR METHOD

The idea is to approximate Ñ˘^m using an explicit method like Euler so we can write the updated solution as:

Û^m+1 = Û̆^m + Ñ^m Δt

that can be used to estimate Ñ˘^m:

Ñ˘^m = N[ Û̆^m - Ñ˘[ Û̆^m-1 + Ñ˘ Δt / 2

where as means that is the Euler approximated solution not the CN one. This allows us to compute Ñ˘^m with only available information like so:

Û^m+1 = Û̆^m + ( Ñ˘^m Û̆^m ) Δt / 2 ^*

However CN is second order accurate while Euler is just first. We need to verify the accuracy.

by using Taylor on Ñ˘^m as:

Ñ ˘^m = Ñ˘^m + ( dÑ˘^m / dṪ ) ( Û̆^m - Û̆˘^m ) + O ( ( Û̆˘^m - Û̆^m+1 ) )

Now taking into account the truncation error of Euler and CN:

Û^m+1 = Û̆^m + O(Δt²) Ê_error approx solution = exact + trunc error

Û̆^m = Û˘^m + O(Δt³) CNtheir difference is: Û^m+1 - Û̆^m+1 = O(Δt²)