Appunti Fluidodinamica computazionale

CFD

2016 - 2017

APPUNTI DI LUCA VERZEROLI

Consideriamo per facilità una griglia di rette

parallele e cartesiana. In esso il valore della

funzione u al punto (x_i,y_i) all'istante di

tempo tⁿ può essere scritto come:

uⁿ_i,j ≈ U(x_i, y_j, tⁿ)

Vediamo come procedere con la discretizzazione

della derivata con il metodo delle differenze

finite. Per definizione la derivata parziale

∂u/∂x

calcolata nel punto x_i, y_j è definita

come

lim_{Δx → 0}

(u(x₀+Δx_i,y₀)-U(x₀,y₀))/Δx

Essendo a livello computazionale impossibile

assegnare per un dominio continuo, la

funzione U(x,y) sarò noto solo in alcuni punti

(x_i, y_j) che sono anche quelli della griglia. Per

cui il limite precedente verrà così scritto:

∂u/∂x_1,3 ≈ (u_i+1,j - u_i,j)/Δx

con Δx sufficientemente piccolo. Vediamo

come determinare l’errore di approssimazione.

Scriviamo l’espansione in serie di Taylor della

funzione U(x,y) nel punto (x_i, y_j) :

u_i+1,j = u_i,j +

∂u/∂x_i,j * Δx + (∂²u/∂x²)_i,j * Δx²/2! +

O((∂³u/∂x³)_i,j * Δx³/3!)

dove ξ è un punto nell’intervallo

[x_i, x_i+1]

U_i+1,j - U_i,j = ^∂u⁄_∂xi,j Δx + ^∂2u⁄_∂x²i,j ^Δx2⁄_2! + ^∂3u⁄_∂x³i,j ^Δx3⁄_3! + o(Δx³)

Sommando due espansioni di Taylor:

U_i+1,j + U_i-1,j - 2U_i,j = ^∂2u⁄_∂x²i,j (Δx²) + O[(Δx)⁴]

^∂2u⁄_∂x²i,j = ^{U_i+1,j - 2U_i,j + U_i-1,j}⁄_(Δx)² + Ο (Δx)²

Da ricordare anche le formule dei 3 punti in avanti e indietro:

^∂2u⁄_∂x²i,j = ^{u_i,j - 2u_i+1,j + U_i+2,j}⁄_(Δx)² + O (Δx)²

^∂2u⁄_∂x²i,j = ^{u_i,j - 2u_i-1,j + U_i-2,j}⁄_(Δx)² + Ο (Δx)²

Per le derivate miste quello che si fa è discretizzare prima la derivata in una direzione e poi procedere nella stessa modo con la rimanente. Vediamo l'esempio per

^∂2u⁄_∂x∂yi,j.

^{(∂u)⁄(∂y)}_i+1,j - ^{(∂u)⁄(∂y)}_i,j ⁄ Δx + Ο (Δx)

^∂u⁄_∂yi+1,j = ^{U_i+1,j - U_i,j-1}⁄_Δy + O(Δy)

Si ottiene così:

ν = _cΔt⁄^Δx

Analizziamo ora il RHS

- _cΔx⁄² (1-λ) U_xx: rappresenta il processo di diffusione o dispersione attenuando i gradienti della soluzione. Ciò introduce dissipazione numerica e _cΔx⁄ ² (1-λ) viene chiamata viscosità numerica o artificiale.

- _(Δx)²⁄⁶ (r³ - 3r² + 1) U_xxx può essere identificato come un effetto artificiale pseudo-fisico. Esso porta a dispersione che è un fenomeno di dipendenza della velocità dell’onda dalla lunghezza d’onda. La soluzione esatta può essere rappresentata come la

Come trovare la soluzione in tutti i punti del dominio e che evolve nel tempo? Quello che si deve fare è scrivere l'equazione delle differenze finite per ogni punto intorno al dominio. Ciò che si ottiene è un sistema di equazioni che deve essere risolto per trovare il valore delle incognite nei punti della griglia. Le equazioni sono accoppiate a causa dell'utilizzo di più punti per la descrizione della funzione nel singolo punto della griglia. Così il sistema risulta irriducibile e deve essere risolto considerando tutte le equazioni contemporaneamente.

ESEMPIO 1 (Schema esplicito)

Schema:

U_i+1ⁿ⁺¹ - U_iⁿ⁺¹ U_iⁿ + 2q U_i-1^n-1 + U₁^n-1 + f₁ⁿ

Δt q² (Δx)²

Sottriamo le due precedenti equazioni n, otteniamo:

ε_i+1ⁿ - ε_iⁿ - ∑εⁿ₁ ε_n-1

Δt (Δx)²

(Questa equazione vale solo per l'esempio)

L'ultima equazione scritta può essere

vista come la rappresentazione delle

differenze finite dell'equazione:

∂ε             a² ∂²ε

∂t           ∂x²

che può essere risolta analiticamente con

il metodo di separazione delle variabili.

Avremo quindi:

Ε(x,t)                    ∑ b_n(t) ε^ikmx

                                          +CC

dove CC è il complesso coniugato

Avermo così scritto Ε(x,t) come

può essere visto come una funzione propria di x

di valore di x sempre come numero

k_m = 2πm (m=1,2,3)

sistema di Fourier però può essere fatto altro

come funzioni continue quello che effetto

Non sempre la stima analitica di un parametro può essere fatta, per cui si potrebbe procedere per via numerica.

Passi per l’analisi di stabilità:

1. Sostituire u_numⁿ = u_esⁿ + εⁿ nello schema delle IF che si vuole utilizzare.
2. Usare la decomposizione in fattore per descrivere ε(X,t).
3. Riscrivere l’equazione PDE per l’errore e trovare la soluzione per ogni modo E_m(k,T), trovare l’espressione del fattore di amplificazione (G).
4. Imponendo G ≤ 1 si trova un criterio di stabilità risolvendo la disequazione.

NB Il metodo di Neumann non può essere sempre se applicato a PDE complesse contenenti non linearità in quanto non esiste la soluzione esponenziale {e^m(k,T) = eⁿ e^-k•im}

Per uno schema implicito G(β) ≤ 1 è solamente soddisfatto automaticamente. Molti di questi schemi risolvono con errore incredibilmente nobile.

5

Un approccio comune è quello di porre questi punti al centro della cella, così da avere un valore della quantità rappresentativo per l'intera cella. Altre possibilità sono: posizionare i punti cosa che i bordi del volume siano a metà strada tra i punti di due celle vicine; oppure posizionare i punti negli estremi del volume. Descriviamo ora la proprietà di conservazione globale dei volumi finiti. Vediamo un esempio. 1-D \(0 \leq x \leq l\) C.C: \(q_0, \, x = 0\) \(q_l, \, x = l\) con \(q\) il flusso Scriviamo la conservazione per il generico settore sull'intero dominio: \(\frac{d}{dt} \int_0^l Q dx - \int Q dx + q_l - q_0\) Per ogni cella avremo invece:

\(\frac{d}{dt} \int_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}} Q dx = - \int_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}} Q dx + \Phi_{i+\frac{1}{2}} - \Phi{i-\frac{1}{2}}\)