convezione nei fluidi incomprimibili: analisi ai volumi finiti

CONVEZIONE NEI FLUIDI INCOMPRIMIBILI:
ANALISI AI VOLUMI FINITI
Enrico Nobile
Dipartimento di Ingegneria Navale, del Mare e per l’Ambiente – DINMA
Università degli Studi Trieste, 34127 TRIESTE
15 marzo 2010
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
1 / 143
15 marzo 2010
2 / 143
OUTLINE
Parte I
Fondamenti
1
INTRODUZIONE
2
L’EQUAZIONE DI TRASPORTO
3
L’IDEA DI BASE
La griglia di calcolo
Griglie Cartesiane per il metodo dei volumi finiti
4
DISCRETIZZAZIONE SPAZIALE
Integrali di superficie
Integrali di volume
Integrazione del termine sorgente
Tecniche di interpolazione
Correzione differita
Equazione algebrica finale
Condizioni al contorno
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
Metodo dei volumi finiti per fluidi incomprimibili
Griglie Cartesiane in domini bidimensionali:
I
I
I
I
Generica equazione di conservazione
Problema termofluidodinamico completo
Modifiche necessarie per griglie Cartesiane in domini tridimensionali
Caso monodimensionale: derivazione elementare
Griglie non strutturate:
I
I
I
I
Generica equazione di conservazione
Quantità geometriche e disposizione delle variabili sulla griglia
Ricostruzione del gradiente
Griglie ibride
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
3 / 143
Metodo dei volumi finiti – FVM (Finite Volume Method)
Equazioni algebriche ottenute dal bilancio su ciascun volume di controllo;
Metodo popolare in termofluidodinamica computazionale;
Utilizzato in numerosi pacchetti CFD (Computational Fluid Dynamics)
commerciali;
Semplicità e corrispondenza fisica fra il FVM ed i principi di conservazione;
Spesso adottato nell’introduzione all’uso di tecniche numeriche in trasmissione
del calore e termofluidodinamica.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
4 / 143
Tutte le equazioni di conservazione hanno una struttura simile, a parte l’equazione di
conservazione della massa:
Equazione dell’energia in forma conservativa:
∂
(ρcp t) + ∇ · (ρwcp t) = ∇ · (λ∇t) + q̇
∂ϑ
Equazione di conservazione della quantità di moto, sempre in forma
conservativa:
∂
(ρw) + ∇ · (ρww) = −∇p + ∇ · (µ∇w) − ρ β (t − t0 ) g
∂ϑ
Forma generale dell’equazione di trasporto:
∂
(ρφ) + ∇ · (ρwφ) = ∇ · (Γ∇φ) + s
∂ϑ
φ è un generico scalare, e Γ la proprietà di trasporto molecolare.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
5 / 143
In altre parole:
(Accumulo) + (Convezione) = (Diffusione) + (Sorgente)
Equazione dell’energia: per fluidi incomprimibili con c costante, φ è dato
dall’entalpia specifica (c t), Γ è dato da λ/c, ed il termine sorgente, s, include
l’eventuale generazione interna di calore.
Equazione della quantità di moto: φ è la componente di velocità, Γ è la viscosità
dinamica µ ed il termine sorgente include, oltre al gradiente di pressione, anche
le forze di galleggiamento, se presenti, ed altri campi di forze.
Sarà sufficiente considerare dapprima la soluzione della generica equazione di
trasporto, per passare poi alla soluzione delle equazioni di Navier-Stokes e di
continuità:
I
I
Si assumeranno noti il campo delle velocità w, le proprietà termofisiche ρ e Γ, ed il
termine sorgente s.
Nel caso in cui le velocità siano incognite, come per le equazioni di Navier-Stokes,
la procedura descritta mantiene la sua validità.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
6 / 143
Utilizzo della formulazione integrale, o finita, dell’equazione di conservazione, scritta
per un generico volume di controllo V:
Z ∂
(ρφ) + ∇ · (ρwφ) − ∇ · (Γ∇φ) − s dV = 0
V ∂ϑ
Applicazione del teorema di Gauss, con A superficie del volume V:
Z
Z
Z
Z
∂
(ρφ) dV + ρφw · n dA = Γ∇φ · n dA + s dV
V ∂ϑ
A
A
V
ed in forma più compatta:
Z
∂
(ρφ) dV +
∂ϑ
V
()
Z
00
Z
J · n dA =
A
s dV
V
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
7 / 143
00
Con J indichiamo il vettore flusso specifico della variabile φ, comprensivo del flusso
00
00
00
00
convettivo, Jc , e diffusivo, Jd , e quindi J = J · n rappresenta la componente del
flusso in direzione normale alla superficie:
00
J
J
()
00
00
00
=
Jc − Jd = ρφw − Γ∇φ
=
J · n = ρφw · n − Γ∇φ · n
00
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
8 / 143
Griglie strutturate Cartesiane
Costituite da famiglie, mutuamente ortogonali, di rette parallele.
Linee di ciascuna famiglia, o le celle definite da tali linee, individuate con un set
di due indici (i, j) in due dimensioni, o di tre indici (i, j, k) in tre dimensioni.
Limitata flessibilità geometrica.
Semplicità ed efficienza dei metodi di calcolo basati su tali griglie.
Esempi di griglie strutturate Cartesiane:
Monoblocco
()
Multiblocco
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
9 / 143
Griglie strutturate curvilinee
Costituite da famiglie di linee curvilinee, nelle quali ciascuna linea di una
famiglia non interseca mai una linea della stessa famiglia, ed interseca una sola
volta le linee delle altre famiglie.
Identiche alle griglie Cartesiane dal punto di vista logico.
Maggiore flessibilità geometrica rispetto alle griglie Cartesiane.
Variante in cui le famiglie di linee sono mutuamente ortogonali.
Esempi di griglie strutturate curvilinee:
Monoblocco
()
Multiblocco
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
10 / 143
Griglie non-strutturate
Le più adatte a trattare geometrie complesse di interesse industriale.
Il dominio è suddiviso in elementi, o celle, di forma arbitraria, tipicamente
triangoli o quadrilateri in 2D, e tetraedri ed esaedri in 3D.
Possibilità di addensare la griglia nelle zone di interesse, anche in modo
automatico (griglie adattive).
Maggiore complessità e onere computazionale.
Esempi di griglie non-strutturate:
A triangoli
()
Ibrida
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
11 / 143
15 marzo 2010
12 / 143
Un esempio
Meshatura a esaedri e ibrida (elica navale, M. Morgut, 2009).
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
Griglie Cartesiane per il metodo dei volumi finiti
Il dominio viene suddiviso in un numero finito di volumi di controllo (VC),
come illustrato in figura nel caso bidimensionale.
La griglia, diversamente dalle Differenze Finite, definisce le facce dei VC, e non
i nodi sui quali sono definite le grandezze incognite.
L’approccio più comune è quello di posizionare i nodi, sui quali sono collocate le
variabili incognite, nei centri dei VC.
I nodi al centro dei VC sono indicati con ’◦’, mentre i nodi sul contorno sono
rappresentati dal simbolo ’•’.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
13 / 143
Conservatività globale
L’equazione di conservazione, scritta in forma integrale, può venire applicata a
qualunque volume di controllo arbitrario, ed in particolare ad ogni VC:
Sommando le equazioni ottenute per tutti i VC che costituiscono l’intero
dominio di calcolo, si ottiene nuovamente l’equazione di conservazione globale
per tutto il dominio, poiché i flussi su ciascuna faccia dei VC si cancellano, dato
che vengono valutati con segno diverso nei due VC adiacenti.
Tale proprietà si verifica solo se la procedura di valutazione dei flussi sulle facce
è univoca.
Sebbene tale proprietà (proprietà telescopica) sembri ovvia, nel caso di griglie
non strutturate arbitrarie co-locate, è necessario adottare alcuni stratagemmi per
rispettare tale importante requisito, che costituisce la base, ed il vantaggio
principale, del metodo dei volumi finiti.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
14 / 143
Caso stazionario
Approssimazione numerica degli integrali di superficie e di volume che appaiono
nell’equazione di conservazione scritta in forma integrale.
Adozione di tecniche di interpolazione, per esprimere il valore di alcune
grandezze in punti diversi da quelli in cui sono collocate.
Equazione di conservazione per il caso stazionario
Z
Z
ρφw · n dA =
A
()
Z
Γ∇φ · n dA +
A
s dV
V
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
15 / 143
Notazione
Tipici VC ottenuti da una griglia Cartesiana, assieme alla notazione adottata:
Bidimensionale
Tridimensionale
Il caso 2D può essere visto come un particolare caso 3D nel quale le variabili
non dipendono da z.
Attenzione rivolta al caso 2D: semplice estensione, per griglie Cartesiane, al
caso 3D.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
16 / 143
Flusso netto
Flusso netto attraverso il contorno del VC
Z
00
J · n dA =
A
XZ
k
=
XZ
k
Ak
J00c
J00 · n dA
Ak
· n dA −
XZ
k
J00d · n dA
Ak
Per garantire la conservatività a livello di VC (locale) e dell’intero dominio
(globale), non vi deve essere sovrapposizione fra i VC:
I
I
Ogni faccia appartiene ad un VC, se giace sul contorno, o a due VC se si trova
all’interno del dominio.
Tale proprietà va garantita anche nei casi più complessi, ad esempio griglie mobili
(sliding grids).
Nel seguito sarà sufficiente considerare solo una faccia, quella indicata con e in
figura; espressioni analoghe si ricavano per le altre con opportune sostituzioni
degli indici.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
17 / 143
Approssimazioni
Per calcolare l’integrale di superficie sulla faccia e in modo esatto, sarebbe
necessario conoscere il valore della funzione integranda (J00 · n) su ogni punto
della faccia.
Ciò non è possibile, poiché la variabile φ, e quindi i flussi ad essa associati, sono
noti solo nei nodi (centri dei VC).
È necessario introdurre due approssimazioni:
1
2
L’integrale è espresso in funzione di uno o più valori della variabile sulla faccia del
VC;
I valori della variabile sulla faccia del VC sono approssimati per mezzo dei valori
nodali della variabile.
Rimandando al seguito la determinazione approssimata delle variabili sulle facce
dei volumi di controllo in funzione dei valori nodali, vediamo come possano
venire valutati gli integrali.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
18 / 143
Calcolo integrali di superficie
Regola del punto medio
Z
J00 · n dA = Je00 Ae ≈ Je00 Ae
Fe =
Ae
Regola dei trapezi
Z
J00 · n dA ≈
Fe =
Ae
Ae 00
00
(Jne + Jse
)
2
Approssimazioni di ordine più elevato: formula di Simpson
Z
Fe =
J00 · n dA ≈
Ae
()
Ae 00
00
(Jne + 4Je00 + Jse
)
6
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
19 / 143
15 marzo 2010
20 / 143
Calcolo integrali di volume
Regola del punto medio (2o ordine)
Z
s dV = s ∆V ≈ sP ∆V
SP =
V
che nel caso bidimensionale si esprime come:
SP ≈ sP ∆xi ∆yj
Approssimazione di ordine elevato (4o ordine)
∆x ∆y
SP =
s dV ≈
16 sP + 4 ss + 4 sn
36
V
+4 sw + 4 se + sne + sse + snw + ssw
Z
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
Integrazione del termine sorgente
Significato del termine sorgente s:
I
Nel termine sorgente vanno inseriti tutti quei contributi che non possono venire
inseriti negli altri termini.
Se sP è noto, e non dipende da φ, non ci sono difficoltà.
Se, come spesso accade, esso dipende dalla variable stessa, φ, è necessario
linearizzarlo:
I
I
È opportuno - e conveniente - tener conto di tale dipendenza nella formulazione
dell’equazione discretizzata;
Nell’esplicitare tale dipendenza, essa può essere al più lineare, poichè il risultato
finale verrà inglobato in un sistema di equazioni lineari.
Linearizzazione del termine sorgente
lhs
sP = srhs
P + sP φP
slhs
P
()
con
< 0 e srhs
P > 0.
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
21 / 143
Linearizzazione del termine sorgente
Se sP è funzione non-lineare di φ, è necessario linearizzarla, cioè specificare i
rhs
valori di slhs
P e sP , che possono ambedue essere funzione di φ:
I
I
rhs
Ad ogni iterazione, o passo di tempo, slhs
P e sP andranno ricalcolati sulla base dei
nuovi valori (in store) di φP , che qua indichiamo per brevità con φ∗P (φkP o φnP );
La linearizzazione scelta per sP dovrebbe rappresentare in modo adeguato la
relazione sP = sP (φP ).
La regola vista slhs
P < 0 dovrebbe venir sempre rispettata:
I
I
Tale requisito non sarebbe necessario se, per la soluzione del sistema di equazioni
lineari ad ogni iterazione o passo di tempo, venisse utilizzato un metodo diretto;
Viceversa, questa regola è molto importante se, come accade sempre in CFD, la
soluzione del sistema di equazioni avviene tramite metodi iterativi.
Vediamo nel seguito alcuni esempi.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
22 / 143
Esempi (1)
Esempio 1
Dato:
s = 6 − 3φ
vediamo le possibili linearizzazioni:
1
2
3
lhs
srhs
P = 6; sP = −3. Ovvia e raccomandata.
∗ lhs
srhs
P = 6 − 3φP ; sP = 0. Poco efficiente, ma spesso è l’unico approccio quando
l’espressione del termine sorgente S è molto complicato (es. Modelli di
Turbolenza).
∗ lhs
srhs
P = 6 + 5φP ; sP = −8. La relazione sP = sP (φP ), proposta in tale
formulazione, è più pendente di quella reale (ricordiamo che φ∗P va considerato a
tutti gli effetti una costante). Corrisponde ad un sottorilassamento, con aumento
di peso della diagonale principale della matrice dei coefficienti, e perciò riduce la
velocità di convergenza. Opportuna in presenza di altre, pesanti non-linearità.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
23 / 143
Esempi (2)
Esempio 2
Dato:
s = 4 + 11φ
vediamo le possibili linearizzazioni:
1
2
3
lhs
lhs
srhs
P = 4; sP = 11. Non accettabile, poichè sP > 0.
∗ lhs
srhs
P = 4 + 11φP ; sP = 0. Corretto, poichè non è possibile individuare una
formulazione naturale che dia slhs
P < 0.
∗ lhs
lhs
srhs
P = 4 + 15φP ; sP = −4. Creazione artificiale di sP < 0. Può rallentare la
convergenza, ma talvolta può portare a benefici (robustezza, minori rischi di
divergenza della simulazione).
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
24 / 143
Esempi (3)
Esempio 3
s = 3 − 5φ3
1
2
3
∗3 lhs
srhs
P = 3 − 5φP ; sP = 0. Poco efficiente, non sfrutta la dipendenza di S da φ.
lhs
∗2
srhs
P = 3; sP = −5φP . In apparenza corretto, ma la relazione sP = sP (φP ) è più
pendente di quella proposta in tale formulazione.
Metodo raccomandato:
∗2
∗
s = s∗ + (ds/dφ)∗ (φP − φ∗P ) = 3 − 5φ∗3
P − 15φP (φP − φP )
∗3
lhs
∗2
∗
quindi: srhs
P = 3 + 10φP ; sP = −15φP : tangente, in φP , alla curva
sP = sP (φP ).
4
∗3 lhs
∗2
srhs
P = 3 + 20φP ; sP = −25φP . Più pendente della sP = sP (φP ), potrebbe
rallentare la convergenza.
Le linearizzazioni ora viste sono illustrate nella slide successiva.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
25 / 143
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
26 / 143
Esempi (3)
()
Interpolazioni
La valutazione degli integrali richiede di conoscere il valore delle variabili in
posizioni diverse da quelle in cui esse sono definite (centri dei VC).
La funzione integranda è il prodotto di diverse variabili, e/o gradienti delle stesse:
Jc00 = ρφw · n per il flusso convettivo, e Jd00 = Γ∇φ · n per il flusso diffusivo.
Per calcolare il valore dei flussi è necessario individuare il valore di φ e la
componente del gradiente normale alla faccia del VC, su uno o più punti di
quest’ultima.
Distinzione fra componente normale del gradiente - flusso diffusivo - e valore
della variabile - flusso convettivo.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
27 / 143
15 marzo 2010
28 / 143
Schemi per il flusso diffusivo
CDS (Central Difference Scheme) 2o ordine
φe = φE λe,PE + φP (1 − λe,PE )
λe,PE =
xe − xP
xE − xP
Componente del gradiente normale alla faccia:
φE − φP
∂φ
≈
∂x e
xE − xP
CDS 4o ordine
φ (x) = a + b x + c x2 + d x3
Per griglie uniformi:
()
∂φ
∂x
≈
e
27 φE − 27 φP + φW − φEE
24 ∆x
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
Considerazioni sugli schemi per il flusso diffusivo
Uso, per i flussi diffusivi, di schemi del 2o ordine;
Schemi di ordine elevato aumentano l’onere computazionale - maggiore
dimensione della molecola di calcolo;
Minore robustezza degli schemi di ordine elevato;
Utilizzo della correzione differita;
Approssimazioni di ordine più elevato non garantiscono una soluzione più
accurata:
I
I
La griglia dev’essere sufficientemente fine.
Ricorso al raffinamento della griglia.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
29 / 143
Flusso diffusivo
Diffusività sulla faccia
Valore della diffusività Γe sulla faccia del VC in presenza di variazioni spaziali di tale
proprietà:
−1
λe,PE
(1 − λe,PE )
+
Γe =
ΓP
ΓE
Per griglia uniforme media armonica di ΓP e ΓE :
Γe =
2ΓE ΓP
ΓE + Γ P
Flusso diffusivo specifico
00
Jd,e
()
= (Γ ∇φ · n)e = Γe
∂φ
∂x
e
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
30 / 143
Schemi base per il flusso convettivo
Flusso convettivo specifico:
00
Jc,e
= (ρ φ w · n)e = ρφe ue
ue è la componente di velocità normale alla faccia e, ed è assunta nota.
CDS - Central Difference Scheme
φe = φE λe,PE + φP (1 − λe,PE )
UDS - Upwind Difference Scheme
(
φe =
()
se (w · n)e > 0
se (w · n)e < 0
φP
φE
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
31 / 143
Diffusività numerica
Lo schema UDS garantisce assenza di oscillazioni, ma introduce una rilevante
diffusione numerica.
Nell’ipotesi (w · n)e > 0:
φe = φP + (xe − xP )
∂φ
∂x
(xe − xP )
+
2
P
2
∂2φ
∂x2
+ ...
P
Lo schema UDS contiene solo il primo termine a destra, e quindi ha
un’accuratezza del primo ordine, mentre l’errore di troncamento (secondo
termine a destra) ricorda proprio l’espressione del flusso diffusivo.
Lo schema UDS introduce quindi un falso flusso diffusivo dato dalla:
∂φ
∆xP
00
JNum
= ΓNum
con
ΓNum = ρ ue
∂x P
2
(1)
La diffusione artificiale aumenta inoltre nei casi in cui le linee di corrente non
siano allineate alla griglia.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
32 / 143
Uso dello schema UDS
Come mai lo schema UDS (o il suo cugino HY - Hybrid, nel quale si utilizza lo
schema CDS nei casi in cui il flusso convettivo sia modesto, passando all’UDS negli
altri casi) è ancora diffuso nei prodotti CFD commerciali?
Semplicità ed assenza di oscillazioni: maggiori garanzie di ottenere una prima
soluzione (approssimata) del problema;
Avviamento della simulazione nei casi difficili, continuando poi con schemi più
accurati;
Utilizzo con schemi di ordine più elevato in modalità correzione differita.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
33 / 143
Schemi di ordine elevato per il flusso convettivo
SOUDS - Second Order Upwind Difference Scheme
(
φe =
φP + [λe,PW (φW − φP )]
φE + [λe,E EE (φEE − φE )]
λe,PW =
Nel caso (w · n)e > 0:
xe − xP
;
xW − xP
φe = φP +
se (w · n)e > 0
se (w · n)e < 0
λe,E EE =
xe − xE
xEE − xE
φW − φP
(xe − xP )
xW − xP
QUICK - Quadratic Upstream Interpolation for Convective
Kinematics
(
φe =
φP + [γ1e φW − (γ1e + γ2e ) φP + γ2e φE ] se (w · n)e > 0
φE + [γ3e φP − (γ3e + γ4e ) φE + γ4e φEE ] se (w · n)e < 0
γ1e = λe,EW λe,PW
γ3e = λe,EP λe,EE P
()
γ2e = λe,PE λe,WE
γ4e = λe,P EE λe,E EE
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
34 / 143
Sugli schemi di ordine elevato per il flusso diffusivo
Sebbene più robusti dello schema CDS, e nel contempo più accurati dello
schema UDS, gli schemi di ordine elevato possono comunque dar luogo ad
oscillazioni non-fisiche della soluzione.
Per rendere più stabili, e privi di oscillazioni, gli schemi ora visti ed altri di
ordine ancora più elevato, i sistemi più diffusi possono venire raggruppati in due
classi, flux blending methods e composite flux limiter methods.
Flux blending methods:
I
I
Un flusso anti-diffusivo viene sommato ad uno schema upwind del 1o ordine (es.
FCT - Flux Corrected Transport).
Una certa quantità di diffusività artificiale viene sommata - in modo selettivo - ad
uno schema di ordine elevato per smorzare le oscillazioni (es. FRAM - Filtering
Remedy and Methodology).
Composite flux limiter methods.
Il flusso numerico sulla faccia della cella viene modificato attraverso un criterio
che imponga di limitarne il valore (boundedness):
I
I
Total Variational Diminishing (TVD) flux limiters.
Normalized Variable Formulation (NVF), o Normalized Variable and Space
Formulation (NVSF).
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
35 / 143
Cenni ad altri schemi di ordine elevato
Su tali basi sono stati sviluppati alcuni schemi di ordine elevato quali:
I
I
I
I
I
MSOU - Monotonic Second Order Upwind scheme;
HLPA - Hybrid Linear/Parabolic Approximation;
SHARP - Simple High-Accuracy Resolution Program);
SMART - Sharp and Monotonic Algorithm for Realistic Transport;
NIRVANA - Nonoscillatory, Integrally Reconstructed Volume-Averaged Numerical
Advection.
Svantaggi di tali schemi:
I
I
Notevole dimensione della molecola di calcolo, a meno di non utilizzare la
correzione differita;
Difficile estensione di tali schemi alle griglie non strutturate.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
36 / 143
Esempio - descrizione
Comportamento di alcuni schemi convettivi sul test del trasporto bidimensionale
puramente convettivo di uno scalare:
Campo di moto assegnato e costante, inclinato di un angolo α = π/4 rispetto
alla griglia.
Condizioni al contorno indicate in figura, e l’equazione di trasporto considerata è
∇ · (ρwφ) = 0, che si traduce nella:
u
∂φ
∂φ
+v
=0
∂x
∂y
Mancando la diffusione, il gradino presente nel profilo di φ all’ingresso,
dovrebbe venire trasportato, senza alcuna diffusione, indefinitamente ed
inalterato.
La presenza di diffusione artificiale, ed eventuali oscillazioni, alterano il risultato.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
37 / 143
Esempio - risultati
I profili, ottenuti sulla sezione orizzontale posta a metà altezza, sono confrontati con il
profilo esatto, cioè quello più accurato ottenibile con la griglia utilizzata.
Lo schema UDS non da luogo a oscillazioni, ma altera completamente il profilo,
a seguito della diffusività numerica che lo caratterizza.
Lo schema CDS preserva in modo accettabile il profilo (la griglia è
particolarmente rada), ma introduce delle oscillazioni.
Lo schema QUICK ha un comportamento intermedio fra UDS e CDS.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
38 / 143
Correzione differita
Difficoltà nell’uso di schemi di ordine elevato.
Insufficiente accuratezza degli schemi semplici.
Correzione differita (Deferred Correction)
Indicando con k l’iterazione k-esima si ha:
oe, k−1
ob, k−1
k
+
φ
−
φ
φke = φob,
e
e
} | e
| {z
{z
}
lhs
rhs
ob indica uno schema di ordine basso (es. UDS), e oe uno di ordine elevato (es.
QUICK).
Formulazione delle espressioni per lo schema SOUDS e per lo schema QUICK.
Basso costo computazionale, migliore convergenza, ridotta possibilità di
oscillazioni non fisiche della soluzione.
Necessità di una procedura iterativa, già presente nella metodologia di soluzione
dei problemi termofluidodinamici (procedura di tipo iterativo o, in modo
equivalente, di integrazione temporale).
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
39 / 143
Equazione algebrica finale
Problemi bidimensionali
AP φP + AE φE + AW φW + AN φN + AS φS = SP
Problemi tridimensionali
AP φP + AE φE + AW φW + AN φN + AS φS + AT φT + AB φB = SP
Forma compatta
A P φP +
X
Anb φnb = SP
nb
con nb (da neighbor) si intende che la sommatoria va estesa ai VC adiacenti.
Le espressioni da utilizzarsi per i coefficienti AP , AE , AW etc. dipendono dagli schemi
adottati per i termini convettivo e diffusivo, e dalle relazioni usate per la valutazione
degli integrali.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
40 / 143
Esempio di calcolo dei coefficienti
Usando lo schema CDS del 2◦ ordine per i flussi convettivo e diffusivo, e
dell’integrazione secondo la regola del punto medio:
ṁe = ρ ue ∆yj
AcE = ṁe λe,PE
Γe ∆yj
AdE = −
xE − xP
AE = AcE + AdE
AP = − (AE + AW+ AN + AS ) − slhs
P ∆xi ∆yj
SP = srhs
P ∆xi ∆yj
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
41 / 143
Condizioni al contorno
Le condizioni al contorno vanno applicate prima di assemblare, e
successivamente risolvere, il sistema di equazioni.
Per le griglie Cartesiane, o più in generale strutturate, sono possibili due diversi
approcci:
b
Metodo 1: nodo sul
contorno
()
Metodo 2: cella fantasma
(ghost cell)
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
42 / 143
Condizione al contorno convettiva - metodo 1
n · Γ ∇φ + α (φBC − φ∞ ) = 0
In altra forma:
−Γ
∂φ
+ α φBC = α φ∞
∂x
ed in termini discreti:
φP − φW
+ α φW = α φ∞
∆x1 /2
#
"
2Γ
α ∆x1
φW = φP
+
φ∞
2 Γ + α ∆x1
2 Γ + α ∆x1
|
{z
} |
{z
}
Clhs
Crhs
−Γ
Equazione algebrica per la cella P:
(AP + AW Clhs ) φP + AE φE + AN φN + AS φS = (SP − AW Crhs )
|
{z
}
|
{z
}
AP
SP
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
43 / 143
Condizione al contorno di Dirichlet - metodo 2
È la condizione al contorno più semplice, e corrisponde, nel caso dell’equazione
dell’energia, al valore imposto della temperatura. Detto φBC il valore assegnato, si ha:
φBC =
φP + φW
2
da cui:
φW = 2 φBC − φP
Equazione algebrica per la cella P:
(AP − AW ) φP + AE φE + AN φN + AS φS = (SP − 2 AW φBC )
| {z }
|
{z
}
AP
SP
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
44 / 143
OUTLINE
Parte II
Integrazione temporale
5
INTEGRAZIONE TEMPORALE
Metodi base
Metodi multilivello e metodi predictor-corrector
Applicazione alla generica equazione di trasporto
Metodi particolari
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
45 / 143
Integrazione temporale
Equazione di trasporto nella formulazione generale
tempo-variante:
∂
(ρφ) + ∇ · (ρwφ) = ∇ · (Γ∇φ) + s
∂ϑ
Integrazione nel tempo alle differenze finite: l’intervallo temporale ϑ è suddiviso
in un certo numero di intervalli temporali ∆ϑ.
Integrazione temporale marciando nel tempo, nota la condizione iniziale al
tempo ϑ = ϑ0 , e le condizioni al contorno in ciascun istante successivo
ϑ = l ∆ϑ, con l = 1, . . . , n, n + 1, . . .
L’equazione può essere scritta nella forma compatta:
∂
(ρφ) = F (ϑ, φ (ϑ))
∂ϑ
F = −∇ · (ρwφ) + ∇ · (Γ∇φ) + s
Questa ricorda la forma delle equazioni differenziali ordinarie ai valori iniziali.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
46 / 143
ODE - Ordinary Differential Equations
Equazione differenziale ordinaria (ODE) del primo ordine, assieme alla condizione
iniziale al tempo ϑ0 :
d φ (ϑ)
= f (ϑ, φ (ϑ))
dϑ
con φ (ϑ0 ) = φ0
Si desidera determinare φ dopo un intervallo ∆ϑ dall’istante iniziale, ottenendo così
φ1 all’istante ϑ1 = ϑ0 + ∆ϑ; successivamente si determinerà φ2 all’istante
ϑ2 = ϑ1 + ∆ϑ e così via. Nel modo più semplice, integrando la precedente fra ϑn e
ϑn+1 = ϑn + ∆ϑ:
ϑ
Zn+1
d φ (ϑ)
dϑ = φn+1 − φn =
dϑ
ϑ
Zn+1
f (ϑ, φ (ϑ)) dϑ
ϑn
ϑn
dove φn+1 = φ (ϑn+1 ). L’espressione scritta è esatta, tuttavia il termine a destra non
può venire valutato senza conoscere la soluzione al nuovo istante, ed è quindi
necessario approssimarlo. Vediamo in breve i metodi più comuni.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
47 / 143
15 marzo 2010
48 / 143
(1) - Metodo di Eulero esplicito
φn+1 = φn + f (ϑn , φn ) ∆ϑ
(2) - Metodo di Eulero implicito
φn+1 = φn + f ϑn+1 , φn+1 ∆ϑ
(3) - Metodo del punto medio
φ
n+1
n
n+1/2
= φ + f ϑn+1/2 , φ
∆ϑ
(4) - Metodo di Crank-Nicolson
φn+1 = φn +
()
1
f (ϑn , φn ) + f ϑn+1 , φn+1 ∆ϑ
2
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
Stabilità ed accuratezza dei metodi base
I metodi visti, a parte (3), sono detti a due livelli.
A parte il metodo (1), tutti gli altri metodi richiedono la funzione incognita ad
istanti diversi da ϑn , per il quale è nota, e quindi è necessario procedere ad
approssimazioni e/o iterazioni.
Il metodo (1) è detto esplicito, mentre gli altri sono metodi impliciti.
Il metodo (1), come tutti i metodi espliciti, è più economico e semplice, ma non
può essere utilizzato per valori di ∆ϑ superiori al limite di stabilità:
∂
f
(ϑ,
φ)
<2
∆ϑ
(2)
∂φ I metodi impliciti sono più costosi, ma sono incondizionatamente stabili:
I
I
metodi impliciti – scala temporale maggiore del limite di stabilità metodi espliciti;
metodi espliciti – scale temporali dello stesso ordine del limite di stabilità.
I metodi (1) e (2) hanno un’accuratezza temporale del prim’ordine (∆ϑ), mentre
2
i metodi (3) e (4) hanno un’accuratezza del second’ordine, ((∆ϑ) ).
Con metodi a due livelli, l’accuratezza può essere, al più, del second’ordine.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
49 / 143
Metodi multilivello e metodi predictor-corrector
Metodi di ordine più elevato:
Metodi multilivello – scelti fra gli istanti di tempo già considerati, e per i quali la
soluzione è disponibile: metodi espliciti di Adams-Bashfort:
1
3 f (ϑn , φn ) − f ϑn−1 , φn−1 ∆ϑ
φn+1 = φn +
2
1 φn+1 = φn +
23 f (ϑn , φn ) − 16 f ϑn−1 , φn−1 + 5 f ϑn−2 , φn−2 ∆ϑ
12
Metodi di Runge-Kutta – compresi fra ϑn e ϑn+1 : metodo del quarto ordine:
φ∗n+1/2
φ∗∗
n+1/2
φ∗n+1
φn+1
()
∆ϑ
f (ϑn , φn )
2
∆ϑ ∗
n
f ϑn+1/2 , φn+1/2
= φ +
2 = φn +
= φn + ∆ϑ f ϑn+1/2 , φ∗∗
n+1/2
∆ϑ h
n
n
∗
= φ +
f (ϑn , φ ) + 2 f ϑn+1/2 , φn+1/2 +
6
i
∗∗
∗
2 f ϑn+1/2 , φn+1/2 + f ϑn+1 , φn+1
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
50 / 143
Applicazione alla generica equazione di trasporto
Caso bidimensionale, metodo implicito di Eulero del prim’ordine e schema CDS :
n+1
AP φn+1
+ AE φn+1
+ AW φn+1
+ AS φn+1
= SP
P
E
W + A N φN
S
ṁe
=
ρ ue ∆yj
AcE
=
AdE
=
AE
=
AP
=
SP
=
ṁe λe,PE
Γe ∆yj
−
xE − xP
c
AE + AdE
ρ
∆xi ∆yj − (AE + AW + AN + AS ) − slhs
P ∆xi ∆yj
∆ϑ
n
ρ
φ
srhs
∆xi ∆yj
P +
∆ϑ
Il termine non-stazionario aumenta il peso del coefficiente AP .
Utilizzo della forma non-stazionaria delle equazioni per problemi stazionari.
Gli schemi espliciti non richiedono la soluzione di sistemi di equazioni, ma il
passo di integrazione temporale è limitato.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
51 / 143
Metodi particolari
Schemi indirizzati ad applicazioni particolari o caratterizzati da particolari vantaggi.
Schema di Eulero implicito del secondo ordine, o schema Gear:
dφ
3 φn+1 − 4 φn + φn−1
≈
dϑ n+1
2 ∆ϑ
4
1
2
φn+1 = φn − φn−1 + f ϑn+1 , φn+1 ∆ϑ
3
3
3
Metodi misti esplicito-implicito, rispettivamente per il flusso convettivo e per il flusso
diffusivo, usati in simulazioni LES e DNS su geometrie semplici.
Schema di Adams-Bashfort del second’ordine per il termine convettivo, e
schema Gear per il contributo diffusivo:
3φn+1 − 4φn + φn−1
2 ∆ϑ
1
3 ∇ · (ρwφn ) − ∇ · ρwφn−1 +
2
n+1
(∇ · (Γ∇φ))
+ sn+1
= −
Schema Runge-Kutta per i termini convettivi, e Crank-Nicolson o Gear per i
contributi diffusivi.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
52 / 143
OUTLINE
Parte III
Sistemi di equazioni
6
METODI DI SOLUZIONE DEI SISTEMI DI EQUAZIONI LINEARI
Metodi diretti
Metodi iterativi
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
53 / 143
Soluzione dei sistemi di equazioni - 1
Sistema di equazioni:
AΦ = S
dove A è la matrice dei coefficienti, Φ è il vettore delle incognite, e S il vettore dei
termini noti:






A11 A12 . . . A1N
Φ
S


1 
1 






 A21 A22 . . . A2N 
 Φ2 

 S2 



A= .
Φ=
S=
..
.. 
..
..
..



 ..
.
.
. 
.
. 












AN1 AN2 . . . ANN
ΦN
SN
La matrice A è sparsa.
La disposizione degli elementi non nulli, all’interno di A, dipende dalla modalità
di numerazione delle incognite.
L’approccio più comune, per griglie strutturate, è quello di utilizzare la
numerazione lessicografica:
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
54 / 143
Soluzione dei sistemi di equazioni - 2
Numerazione lessicografica:
Per un problema 2D con Nj = 3 ed Ni = 4:
AP1
 AN2



 AW4




A=









AS1
AP2
AN3
AE1
AS2
AP3
AE2
AE3
AP4
AN5
AW5

AW6
AS4
AP5
AN6
AE4
AS5
AP6
AW7
AE5
AE6
AP7
AN8
AW8
AW9
AS7
AP8
AN9
AE7
AS8
AP9
AW10
AE9
AP10
AN11
A11
AW12
()
AE8
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
AS10
AP11
AN12
AS11
AP12


















15 marzo 2010
55 / 143
Soluzione dei sistemi di equazioni - 3
La matrice A possiede una struttura polidiagonale, detta anche diagonale a blocchi o
tridiagonale a blocchi.
Non è necessario memorizzare l’intera matrice.
Per tali matrici esistono efficienti algoritmi di soluzione di tipo iterativo e, per
casi particolari, anche diretto.
La struttura della matrice non cambia – alcuni termini sono trattati in modalità
correzione differita - lavorando con griglie curvilinee.
La struttura della matrice si conserva, all’interno di ciascun blocco, anche per
griglie strutturate multiblocco.
La struttura della matrice, sempre sparsa, cambia utilizzando griglie non
strutturate.
Distinzione fra metodi diretti e metodi iterativi:
I
I
Nei primi la soluzione è ottenuta in un numero predeterminato di operazioni,
funzione del numero di incognite.
Nei metodi iterativi la soluzione è ottenuta attraverso iterazioni successive partendo
da un valore di tentativo, ed il numero di operazioni non è predeterminabile, ma
dipende da diversi fattori.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
56 / 143
Eliminazione di Gauss - 1
Fra i metodi diretti il più semplice è il metodo di Gauss. N è il numero totale delle
incognite, pari a Ni × Nj per problemi 2D e Ni × Nj × Nk in 3D.
Eliminazione di A21 dalla matrice A:
I
I
Si moltiplica la prima equazione (prima riga della matrice) per A21 /A11 e la si
sottrae dalla seconda equazione.
Anche tutti gli altri elementi della seconda riga risulteranno modificati, così come il
secondo elemento S2 del termine noto.
Poi si moltiplica la prima equazione per A31 /A11 e la si sottrae dalla terza.
Si procede così per tutti gli altri elementi della prima colonna della matrice.
Ai1
A0ij = Aij −
A1j
A11
i = 2, . . . N, j = 1, . . . N
A
i1
Si0 = Si −
S1
A11
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
57 / 143
Eliminazione di Gauss - 2
Dopo questo primo passaggio la matrice risulta così modificata:


A11 A12 . . . A1N
 0
A022 . . . A02N 


A0 =  .
.
.
.
..
..
.. 
 ..

0
A0N2
...
A0NN
Nessuna delle equazioni 2, 3, . . . , N contiene più la la variabile φ1 .
Si passa quindi ad eliminare la variabile φ2 da questo sistema ridotto, nello
stesso modo, e si continua quindi per le colonne 3, 4 . . . , N − 1.
In termini generali:
A0ik
Akj
A0kk
A0ij
= A0ij −
Si0
A0ik 0
0
= Si − 0 Sk
Akk
()
k = 1, . . . N − 1, i = k + 1, . . . N, j = k . . . , N
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
58 / 143
Eliminazione di Gauss - 3
I fattori A0ik /A0kk sono detti moltiplicatori, o coefficienti moltiplicativi, e la
matrice risultante U è in forma triangolare superiore:




A11 A12 . . . A1N
S

1 



 0 A022 . . . A02N 
 S20 



0
U= .
S =
..
.. 
..
..

 ..
.
.
. 
. 



 0 

0
0
0 . . . ANN
SN
Il sistema così modificato è facilmente risolto. Osservando che:
φN =
SN0
A0NN
è sufficiente procedere all’indietro:
Si0
φi =
()
−
N
X
A0ij φj
j=i+1
A0ii
i = N − 1, . . . 1
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
59 / 143
Eliminazione di Gauss - 4
La sequenza di operazioni che, partendo dalla matrice triangolare superiore,
fornisce i valori delle incognite, è detta sostituzione all’indietro
(backsubstitution).
Per N sufficientemente grande il numero totale di operazioni richieste dal
metodo di Gauss è proporzionale a N 3 /3, sebbene la fase di sostituzione
all’indietro richiede N 2 /2 operazioni, ed è quindi molto meno costosa.
Il costo elevato, in termini di numero di operazioni per N grande, dell’algoritmo
di Gauss, giustifica la ricerca di metodi più economici per la soluzione dei
sistemi di equazioni.
Il metodo di Gauss trova scarso utilizzo, in questa forma, nei problemi di
termofluidodinamica computazionale.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
60 / 143
Decomposizione LU
L’idea è di fattorizzare la matrice A:
A = LU
Per rendere la fattorizzazione unica, si impone che gli elementi diagonali della
matrice L, Lii , siano unitari (o, in modo equivalente, siano unitari gli elementi
diagonali Uii della matrice U). La matrice U è quella ottenuta nella fase di
eliminazione in avanti del metodo di Gauss, mentre gli elementi della matrice L sono
i coefficienti moltiplicativi, Aji /Aii , utilizzati in tale fase.
UΦ = Y
LY = S
La soluzione, in sequenza, di tali due sistemi consente di risolvere il problema.
L’importanza del metodo di decomposizione LU deriva da:
1
la fattorizzazione non richiede la conoscenza del vettore S dei termini noti;
2
la decomposizione LU costituisce la base di alcuni dei migliori metodi iterativi.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
61 / 143
Metodi particolari
Casi speciali di sistemi lineari:
limitati a geometrie semplici;
efficienti per soluzioni ripetute del medesimo sistema con differenti termini noti.
Algoritmi di maggiore interesse nella termofluidodinamica computazionale:
Algoritmi per sistemi tridiagonali e tridiagonali ciclici: Thomas, Temperton, etc..
Il numero di operazioni scala linearmente con N, anzichè con N 3 , come
nell’eliminazione di Gauss.
Algoritmi per griglie Cartesiane (ed anche cilindriche e sferiche) 2D e 3D
uniformi: basati su metodi FFT (Fast Fourier Transform) e di riduzione ciclica. Il
numero di operazioni, escludendo la fase di fattorizzazione, è dell’ordine di
N × (ln2 N).
Algoritmi per griglie Cartesiane (cilindriche e sferiche) 2D e 3D non uniformi:
decomposizione matriciale, per il quale il numero di operazioni, per la sola fase
di soluzione, varia con N 4/3 .
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
62 / 143
Metodi iterativi - 1
L’utilizzo di metodi iterativi è giustificato dalle seguenti osservazioni:
Con il metodo LU, sebbene A sia sparsa, le matrici L ed U sono piene.
Gli errori di discretizzazione sono usualmente ben superiori agli errori di
arrotondamento accumulati con i metodi diretti.
Nei problemi applicativi è necessario risolvere sistemi di rilevante dimensione.
L’idea alla base dei metodi iterativi consiste nel partire da una soluzione di tentativo, e
migliorarla sino al livello desiderato: dopo n iterazioni avremo una soluzione
approssimata Φn , che non soddisfa l’equazione in modo esatto:
AΦn = S − rn
Sottraendo questa dal sistema di equazioni di partenza si ottiene una relazione fra
l’errore di convergenza n , definito dalla:
n = Φ − Φn
n
ed il residuo r :
An = rn
L’obiettivo di un metodo iterativo è quello di diminuire, sino a zero, il residuo, ed in
tal caso anche l’errore di convergenza va a zero.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
63 / 143
Metodi iterativi - 2
Schema iterativo:
MΦn+1 = NΦn + B
La soluzione, a convergenza avvenuta, deve soddisfare l’equazione di partenza. A
convergenza, inoltre, Φn+1 = Φ, da cui:
A=M−N
e
B=S
PA = M − N
e
B = PS
o, più in generale:
dove P è una matrice (non singolare) di precondizionamento.
Versione alternativa dello schema iterativo:
M Φn+1 − Φn = B − (M − N) Φn
o
Mδ n = rn
dove δ n = Φn+1 − Φn è la correzione, o aggiornamento, e rappresenta
un’approssimazione dell’errore di convergenza.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
64 / 143
Metodi iterativi - 3
Un metodo iterativo è tanto più efficace quanto più economica è la soluzione del
sistema di equazioni modificato, e quanto più rapida è la convergenza, cioè
minore il numero di iterazioni necessarie.
La velocità di convergenza è funzione del raggio spettrale λ1 , definito come il
modulo dell’autovalore massimo della matrice M−1 N.
Un’espressione approssimata del numero di iterazioni necessario è dato dalla:
ln aδ1
n≈
ln λ1
dove a1 è una costante e δ è la tolleranza richiesta per l’errore di convergenza.
Da questa relazione approssimata per n si può notare che, per λ1 prossimo
all’unità, la convergenza può risultare molto lenta.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
65 / 143
Metodi di base - 1
Metodo di Jacobi
M è una matrice diagonale, i cui elementi sono quelli della diagonale principale di A:
Φn+1
=
P
SP − AE ΦnE − AW ΦnW − AN Φn − AS ΦnS
AP
Il numero di iterazioni richiesto è proporzionale al quadrato del numero di VC in una
direzione. È quindi un metodo piuttosto lento, e perciò di scarsa utilità.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
66 / 143
Metodi di base - 2
Metodo di sovrarilassamento o SOR (Successive
Over-Relaxation)
La matrice M è costituita dalla porzione triangolare inferiore della matrice A.
Seguendo l’ordine lessicografico:
Φn+1
P
= (1 − ω)
ΦnP
n+1
SP − AE ΦnE − AW Φn+1
− AS ΦnS
W − AN Φ
+ω
AP
dove ω è il fattore di sovrarilassamento, (1.5 ≤ ω ≤ 1.85), necessario per accelerare
la convergenza.
Il metodo SOR è più efficiente del metodo di Jacobi poichè con il SOR si utilizzano i
valori nuovi delle variabili non appena questi si rendono disponibili.
Metodo di Gauss-Seidel
Corrisponde allo schema SOR con ω = 1:
Φn+1
P
()
=
ΦnP
n+1
SP − AE ΦnE − AW Φn+1
− AS ΦnS
W − AN Φ
+
AP
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
67 / 143
Metodi di base - 3
Metodo di sovrarilassamento per linee o SLOR (Successive Line
Over-Relaxation)
Aggiornamento delle variabili di un’intera linea – riga o colonna – di VC.
M è costituita da un sottoinsieme della A, relativo ai soli coefficienti delle variabili
sulla linea in esame.
Procedendo per colonne, e per valori di i crescenti:
i
h
n+1
n+1
n+1
n+1
n
AP Φi,j + AN Φi,j+1 + AS Φi,j−1 = SP − AE Φi+1,j − AW Φi−1,j
Si ottiene un sistema tridiagonale.
Possibile utilizzo di un coefficiente di sovrarilassamento ω.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
68 / 143
Cenni sui metodi basati sul gradiente
Trovano fondamento nel fatto che è possibile minimizzare una funzione, rispetto
a più direzioni, lavorando su una direzione alla volta.
Possono venire classificati in funzione del tipo di sistemi (matrici) lineari a cui si
rivolgono:
1
2
Sistemi simmetrici. CG (Conjugate Gradient) e ICCG (Incomplete Cholesky
Conjugate Gradient).
Sistemi non simmetrici. BCG (BiConjugate Gradient), CGS (Conjugate Gradient
Squared), CGSTAB (CGS Stabilized) e GMRES (Generalized Minimal Residual).
I metodi basati sul gradiente fanno uso di opportuni precondizionatori al fine
migliorarne le prestazioni.
Si ricorda infatti che la velocità di convergenza di tali metodi diminuisce
all’aumentare dell’indice di condizionamento κ della matrice A, definito dalla:
κ=
λmax
λmin
L’idea alla base delle tecniche di precondizionamento è quello di sostituire il
problema in esame con un altro caratterizzato da un indice di condizionamento
inferiore.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
69 / 143
Metodi multigriglia - 1
Basati sull’utilizzo di più griglie: la griglia effettiva di calcolo (livello l), ed una
sequenza di griglie via via più rade (l − 1, l − 2, . . . 1), dove ciascuna cella è ottenuta
sommando (agglomerando) più celle della griglia immediatamente più fine.
L’idea alla base di tali metodi, limitando la descrizione al caso di due sole griglie,
nasce dalle seguenti osservazioni:
Tutti i metodi iterativi sono efficaci nella riduzione degli errori (residui) di
piccola lunghezza d’onda, cioè degli errori di elevata frequenza. Si tratta di
errori la cui lunghezza d’onda è comparabile con la dimensione della griglia.
La riduzione degli errori di piccola lunghezza d’onda avviene nelle prime
iterazioni.
Viceversa, la riduzione degli errori di lunghezza d’onda via via più elevata (bassa
frequenza) avviene più lentamente, ed in effetti, dopo le prime iterazioni, la
velocità di convergenza si riduce.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
70 / 143
Metodi multigriglia - 2
Riduzione degli errori di bassa frequenza, in modo molto più economico, su una
griglia più rada (interpolazione o somma dei residui).
Ottenuta la soluzione sulla griglia rada, questa viene iniettata o sommata sulla
griglia fine.
Le iterazioni sulla griglia più rada sono molto più economiche.
Maggiore velocità di convergenza, poichè sulla griglia rada lo stesso errore ha
frequenza più elevata.
Il vantaggio aumenta utilizzando più griglie, o livelli.
Idealmente, il numero di iterazioni equivalenti, e quindi il tempo di calcolo, varia
linearmente con il numero N di variabili.
Metodi multigriglia di maggiore interesse in CFD:
1
Multigriglia di tipo Additivo - ACM (Additive Correction Multigrid).
2
Multigriglia di tipo Algebrico - AGM (AlGebraic Multigrid).
Come risolutore iterativo nelle varie griglie, detto smoother, può essere utilizzato uno
qualsiasi dei metodi iterativi visti.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
71 / 143
OUTLINE
Parte IV
Problemi termofluidodinamici - approccio FVM
7
SOLUZIONE DEI PROBLEMI TERMOFLUIDODINAMICI
Disposizione delle variabili sulla griglia
8
PROCEDURA AI VOLUMI FINITI
Metodi segregati
Metodi accoppiati
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
72 / 143
PROBLEMI TERMOFLUIDODINAMICI - 1
FVM: approccio comune per la risoluzione numerica di problemi
termofluidodinamici.
Ampia diffusione di codici commerciali basati su tale tecnica.
Le equazioni algebriche ottenute dalla discretizzazione delle equazioni di
Navier-Stokes, ed altre grandezze, potrebbero venire risolte simultaneamente,
oppure in modo sequenziale: la soluzione procede, per ogni variabile,
congelando le altre al passo di tempo o iterazione precedente.
I primi codici FVM utilizzavano l’approccio sequenziale, o segregato, e sole
griglie strutturate di tipo Cartesiano.
Estensione alle griglie strutturate curvilinee, prima monoblocco, poi multiblocco.
FVM per griglie non strutturate:
I
I
I
I
generazione completamente automatica della griglia per geometrie complesse;
importazione di geometrie CAD;
griglie adattive;
metodologie di soluzione di tipo segregato ed accoppiato.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
73 / 143
PROBLEMI TERMOFLUIDODINAMICI - 2
Con il FVM, si utilizzano procedure di calcolo di tipo iterativo per problemi
stazionari o, per problemi non-stazionari, si fa ricorso a tecniche di integrazione
di tipo implicito.
Metodo SIMPLE (Semi-Implicit Method for Pressure-Linked Equations), e
derivati da questo, come SIMPLER, SIMPLEC, SIMPLEST etc.:
I
I
Non differiscono sostanzialmente dai procedimenti di calcolo di tipo Projection;
Per ragioni di omogeneità faremo riferimento ai metodi di proiezione.
Metodi di Proiezione: determinazione di una equazione differenziale per la
pressione a partire dal vincolo di conservazione della massa.
Proiezione: il campo di velocità approssimato, trovato nella prima fase del
calcolo, viene successivamente proiettato in un campo a divergenza nulla
tenendo conto delle correzioni di pressione.
Differenze fra le equazioni di Navier-Stokes e la generica equazione di trasporto:
conseguenze delle possibili distribuzioni delle variabili velocità e pressione sulla
griglia.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
74 / 143
Disposizione delle variabili sulla griglia - 1
Le equazioni di Navier-Stokes e della continuità, in forma conservativa in coordinate
Cartesiane bidimensionali:
∂
∂
∂
∂p
(ρu) +
(ρuu) +
(ρvu) = −
∂ϑ
∂x
∂y
∂x
∂
∂u
∂
∂u
+
µ
+
µ
− ρ0 β (t − t0 ) gx
∂x
∂x
∂y
∂y
∂
∂
∂p
∂
(ρv) +
(ρuv) +
(ρvv) = −
∂ϑ
∂x
∂y
∂y
∂
∂v
∂
∂v
+
µ
+
µ
− ρ0 β (t − t0 ) gy
∂x
∂x
∂y
∂y
∂u ∂v
+
=0
∂x ∂y
Termine sorgente - gradiente di pressione - richiede particolare attenzione.
Il termine di galleggiamento non da luogo a particolari problemi.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
75 / 143
Disposizione delle variabili sulla griglia - 2
Consideriamo la prima delle due equazioni di Navier-Stokes.
Con il FVM il termine relativo al gradiente di pressione viene usualmente
espresso (approccio conservativo) come forza di superficie:
Z
Z
Z
∂p
dV = −
p dA −
p dA ≈ ∆yj (pw − pe )
−
∂x
Ae
Aw
V
Usando l’interpolazione lineare (differenze centrali), avremo:
pW + pP
pP + pE
pW − pE
∆yj (pw − pe ) ≈ ∆yj
−
= ∆yj
2
2
2
Risulta quindi che il termine relativo alla pressione viene valutato utilizzando
una griglia due volte più grossolana.
Debole accoppiamento fra il campo delle velocità e quello delle pressioni:
campo di pressione a scacchiera (checkerboard pressure field).
È quindi necessario eliminare questo inconveniente.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
76 / 143
Disposizione delle variabili sulla griglia - 3
Utilizzo di griglie diverse, per le componenti di velocità, rispetto a quella,
principale.
Griglie ottenute traslando (sfalsando) i VC della griglia principale, per ciascuna
componente di velocità e nella direzione corrispondente, di una quantità pari a
metà del lato del VC.
Con questa disposizione delle variabili (staggered grid), l’integrazione del
termine relativo al gradiente di pressione per ue fornisce:
Z
Z
Z
∂p
−
dV = −
p dA −
p dA ≈ ∆yj (pP − pE )
Ve ∂x
AE
AP
Assenza di campi di pressione checkerboarded: il gradiente viene valutato con
riferimento alla pressione in due nodi adiacenti.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
77 / 143
Disposizione delle variabili sulla griglia - 4
Le griglie sfalsate nel FVM corrispondono all’utilizzo di funzioni di forma di
ordine inferiore per la pressione (unequal order interpolation) nel FEM.
Ulteriori vantaggi offerti da tale tipo di griglia:
I
I
Alcuni termini delle equazioni, che richiederebbero altrimenti delle interpolazioni,
possono venire valutati, con un’accuratezza del second’ordine, in modo semplice.
Viene garantita la conservazione dell’energia cinetica.
Difficoltà di utilizzo delle griglie staggered di tipo curvilineo, ed ancor più non
strutturate.
Utilizzo della distribuzione co-locata (collocated) delle variabili:
La proprietà di conservazione dell’energia cinetica delle griglie sfalsate, le rende
particolarmente indicate in LES e DNS.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
78 / 143
Metodi segregati - 1
Deflussi bidimensionali di fluidi con proprietà costanti.
Termine di galleggiamento: ipotesi di Oberbeck-Boussinesq.
Calcolo delle componenti della velocità di tentativo (o di stima) u∗ e v∗ ,
assumendo la pressione di stima p∗ pari a pn :
u∗ − un
+ ∇ · (ρwn γu∗ ) + ∇ · (ρwn (1 − γ) un )
ρ
∆ϑ
∂p∗
∗
n
= ∇ · (µ∇γu ) + ∇ · (µ∇ (1 − γ) u ) −
− ρβ (tn − t0 ) gx
∂x
ρ
v∗ − vn
+ ∇ · (ρwn γv∗ ) + ∇ · (ρwn (1 − γ) vn )
∆ϑ
∂p∗
= ∇ · (µ∇γv∗ ) + ∇ · (µ∇ (1 − γ) vn ) −
− ρβ (tn − t0 ) gy
∂y
dove i termini convettivi sono linearizzati secondo il metodo di Picard:
ρwγw ≈ (ρwn ) γw
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
79 / 143
Metodi segregati - 2
Portando a destra dell’ uguale i termini che non contribuiscono alla matrice dei
coefficienti, ed integrando sui rispettivi volumi di controllo si ottiene:
Z
Z
Z
ρ
∗
∗ n
u dV + ργ
u w · n dA − µγ
∇u∗ · n dA
∆ϑ Vu
A
Au
Zu
Z
n n
= −ρ (1 − γ)
u w · n dA + µ (1 − γ)
∇un · n dA
Au
Au
"Z
#
Z
Z
Z
ρ
−
p∗ dA −
p∗ dA − ρβgx
(tn − t0 ) dV +
un dV
∆ϑ Vu
AE,x
AP,x
Vu
ρ
∆ϑ
Z
v∗ dV + ργ
Z
∇v∗ · n dA
Vv
A
Av
Zv
Z
= −ρ (1 − γ)
vn wn · n dA + µ (1 − γ)
∇vn · n dA
Av
Av
"Z
#
Z
Z
∗
−
∗
p dA − ρβgy
p dA −
AN,y
()
v∗ wn · n dA − µγ
Z
AP,y
ρ
(t − t0 ) dV +
∆ϑ
Vv
n
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
Z
vn dV
Vv
15 marzo 2010
80 / 143
Metodi
segregati - 3
La correzione di pressione viene calcolata risolvendo l’equazione di Poisson
corrispondente che, integrata sul volume di controllo della pressione fornisce:
Z
Z
ρ
w∗ · n dA
∇p0 · n dA =
∆ϑ A
A
che in forma discreta, utilizzando la regola del punto medio, diventa:
Z
Z
Z
Z
∂p0 ∂p0 ∂p0 ∂p0 dy −
dy +
dx −
dx
∂x
∂x
∂y
∂y
Ae
Aw
An
As
e
w
n
s
ρ
= [u∗e − u∗w ] ∆yj + [v∗n − v∗s ] ∆xi
∆ϑ
I gradienti delle pressioni di correzione sulle facce vengono valutati secondo uno
schema alle differenze centrali, ottenendo:
0
0
p0P − p0W
pN − p0P
p0P − p0S
pE − p0P
−
∆yj +
−
∆xi
xE − xP
xP − xW
yN − yP
yP − yS
∗
ρ
=
ui,j − u∗i−1,j ∆yj + v∗i,j − v∗i,j−1 ∆xi
∆ϑ
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
81 / 143
Metodi segregati - 4
Le componenti u0 e v0 delle correzioni di velocità vengono calcolate in funzione
del gradiente della correzione di pressione:
"Z
#
Z
Z
ρ
u0 dV = −
p0 dA −
p0 dA
∆ϑ Vu
AE,x
AP,x
"Z
#
Z
Z
ρ
v0 dV = −
p0 dA −
p0 dA
∆ϑ Vv
AN,y
AP,y
In forma discreta, ed utilizzando le più semplici espressioni del second’ordine
per l’integrazione e l’interpolazione, fornisce:
()
u0e
=
v0n
=
∆ϑ
ρ
∆ϑ
−
ρ
−
(p0P − p0E )
(xP − xE )
(p0P − p0N )
(yP − yN )
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
82 / 143
Metodi segregati - 5
I valori della pressione e delle componenti di velocità alla fine del passo
temporale sono dati dalle:
pn+1
P
= p∗P + p0P
un+1
e
= u∗e + u0e
vn+1
n
= v∗n + v0n
Per ciò che riguarda la temperatura, questa viene ottenuta integrando sul VC
l’equazione dell’energia, utilizzando lo stesso algoritmo di integrazione
temporale utilizzato per le componenti di velocità:
Z
Z
Z
ρ cp
n+1
n+1 n
t
dV + ρ cp γ t w · n dA − λγ ∇tn+1 · n dA
∆ϑ V
A
Z A
Z
n n
= −ρ cp (1 − γ) t w · n dA + λ (1 − γ) ∇tn · n dA
A
A
Z
ρ cp
+
tn dV
∆ϑ V
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
83 / 143
Metodi segregati - 6
Osservazioni
In analogia con il FEM, le proprietà termofisiche sono state ritenute costanti.
L’algoritmo di integrazione temporale tipo projection ora visto presenta molte
varianti possibili, quale ad esempio quella data dalla combinazione di uno
schema Adams-Bashfort esplicito per il termine convettivo, ed uno schema Gear
(Eulero second’ordine) implicito per il flusso diffusivo.
La metodologia illustrata può servire, viste le notevoli somiglianze, come
riferimento per altri algoritmi di soluzione, come il SIMPLE e derivati.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
84 / 143
Metodi accoppiati - 1
I metodi segregati possono dar luogo a problemi di convergenza quando
l’accoppiamento fra le variabili è molto pronunciato.
Necessità di dover adottare ridotti passi di integrazione temporale o valori piccoli
dei coefficienti di sottorilassamento.
Per queste situazioni può risultare conveniente utilizzare una strategia di
soluzione che preservi, in modo implicito, l’accoppiamento fra le variabili.
Le variabili risolte simultaneamente appartengono tipicamente ad opportuni
sotto-domini, che vengono visitati in sequenza. I sotto-domini sono usualmente
costituiti da:
1
2
Un singolo VC. Metodo SCGS (Simmetrically Coupled Gauss-Seidel). Può essere
utilizzato anche per griglie non strutturate.
Una linea, riga o colonna, di VC. Ad esempio CELS (Coupled Equation Line
Solver). Limitato a griglie strutturate.
(SCGS)
()
(CELS)
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
85 / 143
Metodi accoppiati - 2
Metodo SCGS, sole variabili pressione e velocità nel caso 2D
stazionario
AuPi,j ui,j + AuEi,j ui+1,j + AuWi,j ui−1,j + AuNi,j ui,j+1
u
+AuSi,j ui,j−1 + (pi+1,j − pi,j ) ∆yj = Si,j
AuPi−1,j ui−1,j + AuEi−1,j ui,j + AuWi−1,j ui−2,j + AuNi−1,j ui−1,j+1
u
+AuSi−1,j ui−1,j−1 + (pi,j − pi−1,j ) ∆yj = Si−1,j
AvPi,j vi,j + AvEi,j vi+1,j + AvWi,j vi−1,j + AvNi,j vi,j+1
v
+AvSi,j vi,j−1 + (pi,j+1 − pi,j ) ∆xi = Si,j
AvPi,j−1 vi,j−1 + AvEi,j−1 vi+1,j−1 + AvWi,j−1 vi−1,j−1 + AvNi,j−1 vi,j
v
+AvSi,j−1 vi,j−2 + (pi,j − pi,j−1 ) ∆xi = Si,j−1
(ui,j − ui−1,j ) ∆yj + (vi,j − vi,j−1 ) ∆xi = 0
Il termine relativo al gradiente di pressione è stato scritto in modo esplicito.
Sistema lineare di cinque equazioni nelle incognite ui−1,j , ui,j , vi,j−1 , vi,j e pi,j .
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
86 / 143
Metodi accoppiati - 3
Per poter semplificare il sistema, al fine di ottenere economicamente la soluzione per
via analitica, si adotta la seguente strategia:
1
Per le equazioni di quantità di moto, si portano nel termine noto tutti i contributi,
con l’esclusione del termine diagonale AP e del contributo della pressione pi,j .
2
L’equazione di continuità viene mantenuta in forma completa.
Così facendo si ottiene il seguente sistema:






AuPi−1,j
0
0
0
−∆yj
0
AuPi,j
0
0
∆yj
0
0
AvPi,j−1
0
−∆xi
()
0
0
0
AvPi,j
∆xi
∆yj
−∆yj
∆xi
−∆xi
0
 
ui−1,j



 
  ui,j

  vi,j−1
 
vi,j



pi,j
 u
bi−1,j




 bui,j
bvi,j−1
=






bvi,j




0






Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile











15 marzo 2010
87 / 143
Metodi accoppiati - 4
I termini noti contengono i contributi che non appaiono esplicitamente in forma
matriciale, ad esempio:
u
bui,j = −AuEi,j u∗i+1,j − AuWi,j u∗i−1,j − AuNi,j u∗i,j+1 − AuSi,j u∗i,j−1 − p∗i+1,j ∆yj + Si,j
Soluzione del sistema di equazioni:
r1 = −∆yj /AuPi−1,j ;
r2 = ∆yj /AuPi,j ;
r3 = −∆xi /AvPi,j−1 ;
r4 = ∆xi /AvPi,j
DEN = r1 ∆yj − r2 ∆yj + r3 ∆xi − r4 ∆xi
pi,j = r1 bui−1,j + r2 bui,j + r3 bvi,j−1 + r4 bvi,j /DEN
ui−1,j = bui−1,j − ∆yj pi,j /AuPi−1,j ; ui,j = bui,j + ∆yj pi,j /AuPi,j
vi,j−1 = bvi,j−1 − ∆xi p,j /AvPi,j−1 ; vi,j = bvi,j + ∆xi pi,j /AvPi,j
SCGS applicabile anche alle griglie non strutturate, e altre variabili.
La convenienza dei metodi accoppiati deriva, soprattutto, dal loro utilizzo
contestuale ad un metodo multigriglia.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
88 / 143
OUTLINE
Parte V
Geometrie complesse: griglie non strutturate
9
GEOMETRIE COMPLESSE: GRIGLIE NON STRUTTURATE
Scelta delle componenti di velocità
Disposizione delle variabili sulla griglia e quantità geometriche
Distribuzione spaziale delle variabili ed integrazione
Calcolo del gradiente
Termine transitorio e termine sorgente
Flusso convettivo e flusso diffusivo
Condizioni iniziali ed al contorno
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
89 / 143
Geometrie complesse
La procedura illustrata per griglie strutturate Cartesiane è semplice ed accurata, e
di agevole implementazione, ed infatti ha trovato applicazione in ricerche di base
e nei primi codici commerciali.
Difficoltà di utilizzo per geometrie complesse:
Necessità di modificare il FVM in modo da affrontare anche tali problemi.
Inizialmente ricorrendo a griglie strutturate curvilinee (approccio ancora
utilizzato in qualche programma commerciale), passando poi a griglie non
strutturate.
Tendenza attuale è quella di utilizzare griglie non strutturate: l’uso di griglie
strutturate, o Cartesiane, rimane garantito, sebbene con minore efficienza
computazionale.
Vista la maggiore generalità delle griglie non strutturate, in quanto segue
considereremo solo queste.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
90 / 143
Classificazione FV
Svariate modalità di costruzione dei Volumi Finiti (celle) a partire da griglie non
strutturate.
Le più comuni sono:
(a) Cell-centered
(b) Vertex-centered (o node-centered)
(c) Circum-centered
Nel seguito, faremo riferimento al solo caso cell-centered (a).
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
91 / 143
Generalità e notazione
Dominio discretizzato con una griglia non strutturata, costituita da VC di forma
poliedrica arbitraria con un numero di facce n ≥ 4:
Accuratezza maggiore con celle a forma di esaedri regolari.
Generatori di griglia che producono mesh a poliedri:
I
I
L’uso di poliedri, in generale, può garantire una maggiore accuratezza, a causa del
più elevato grado di ortogonalità della griglia (vedi in seguito).
Per griglie node-centered i VC sono già dei poliedri centrati sui nodi.
Operando con FVM, non è necessario effettuare trasformazioni di coordinate.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
92 / 143
Centroidi
Centroide C0 (baricentro) del VC
Z
(r − rC0 ) dV = 0
V
Centroide cj della faccia Aj compresa fra C0 e Cj
Z
(r − rj ) dA = 0
Aj
Aj è il vettore area della faccia:
Aj = A j n
In generale dj e Aj non sono paralleli, e l’angolo da essi formato aumenta per
griglie molto distorte.
Il punto di applicazione di Aj , cioè cj , è in generale diverso dal punto ottenuto
dall’intersezione di dj con la faccia stessa.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
93 / 143
Equazione di trasporto
È importante operare con griglie nelle quali i due vettori siano il più possibile
paralleli e disposti lungo la stessa retta, e cioè con griglie il più possibile
ortogonali.
Nel FVM l’ortogonalità si misura fra congiungente i centroidi dj , e vettore area
Aj della faccia interposta.
Generica equazione di trasporto integrata sul VC:
Z
V
∂
(ρφ) dV +
∂ϑ
Z
Z
ρφw · n dA =
A
Z
Γ ∇φ · n dA +
A
s dV
V
In analogia a quanto visto per le griglie Cartesiane, considereremo singolarmente i
vari termini dell’equazione di trasporto, dopo aver prima definito alcuni aspetti della
procedura di calcolo.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
94 / 143
Scelta delle componenti di velocità
Nel metodo FV, le equazioni della quantità di moto sono espresse in forma
conservativa:
I
I
I
Tutti i termini delle equazioni possono venire rappresentati come divergenza di un
vettore o di un tensore.
L’utilizzo della forma conservativa delle equazioni con il FVM, garantisce la
conservazione globale della quantità di moto anche nel calcolo discreto.
Tuttavia, la forma conservativa delle equazioni richiede di esprimere le componenti
della quantità di moto in un sistema di riferimento fisso, pena l’introduzione di
campi di forza non-conservativi, secondo la definizione vista.
Per tale ragione, è opportuno utilizzare componenti di tale tipo, ed in particolare
Cartesiane, che sono le più semplici ed intuitive.
Nel seguito, faremo riferimento a grandezze (vettori e tensori) espresse in
termini di componenti Cartesiane.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
95 / 143
Disposizione delle variabili sulla griglia - 1
La disposizione sfalsata (staggered) delle variabili sulla griglia garantisce un
buon accoppiamento fra componenti di velocità e pressione.
Ciò vale solo per griglie strutturate Cartesiane, per le quali le componenti della
velocità risultano normali alle facce dei VC.
Per griglie non strutturate, viceversa, sarebbe necessario memorizzare tutte le
componenti della velocità su ciascuna faccia.
L’utilizzo di componenti di velocità allineate localmente alla griglia (componenti
contravarianti), risolve tale problema, ma non sarebbe garantita la conservatività
globale della quantità di moto.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
96 / 143
Disposizione delle variabili sulla griglia - 2
Pratica diffusa co-locare tutte le variabili nel centroide dei VC:
Agevole gestione della struttura dati e programmazione (informazioni
geometriche di un solo set di volumi di controllo).
Utilizzo di opportune procedure di interpolazione per ovviare al debole
accoppiamento fra il campo di velocità ed il campo di pressione.
Nel FVM su griglie non strutturate, l’utilizzo di componenti Cartesiane della
velocità, e la disposizione co-locata delle variabili, rappresenta un buon
compromesso fra semplicità di implementazione, accuratezza e rispetto della
conservatività (applicazione in numerosi codici, commerciali e open source).
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
97 / 143
Quantità geometriche
Il dominio di calcolo viene suddiviso in un numero finito di VC contigui
attraverso una griglia.
La griglia, in analogia con il metodo degli Elementi Finiti, definisce i vertici
(vertices) o nodi, che collegati due a due costituiscono a loro volta gli spigoli
(edges) dei VC.
Gli spigoli definiscono le facce (faces) dei VC, che non sono necessariamente
piane; tuttavia le loro proiezioni sui piani Cartesiani risultano indipendenti dalla
forma effettiva della faccia.
Tali proiezioni sono le componenti Cartesiane del vettore area A, univocamente
definito sulla base delle coordinate dei vertici della faccia.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
98 / 143
Richiami di calcolo vettoriale - 1
Prodotto scalare
Il prodotto scalare di due vettori a e b (dot
product) è usualmente scritto a · b, ed è
definito dalla quantità scalare:
a · b = |a| |b| cos θ
Alcune proprietà del prodotto scalare:
a · (b + c) = a · b + a · c
2
a · a = |a|
a · b = (ax i + ay j + az k) · (bx i + by j + bz k)
= ax bx + ay by + az bz
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
99 / 143
Richiami di calcolo vettoriale - 2
Prodotto vettoriale
Il prodotto vettoriale di due vettori a e b
(cross product) è usualmente scritto a × b
(o a ∧ b), ed è definito dalla quantità
vettoriale:
a × b = |a| |b| sin θ n
Alcune proprietà del prodotto vettoriale:
b × a = −a × b
a × (b + c) = a × b + a × c
a × b = (ax i + ay j + az k) × (bx i + by j + bz k)
i
a×b = ax
bx
j
ay
by
k
az
bz
= (ay bz − az by ) i + (az bx − ax bz ) j + (ax by − ay bx ) k
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
100 / 143
Richiami di calcolo vettoriale - 3
Vettore area A
A = a × b = |a| |b| sin θ n
Volume di un parallelepipedo
V = c·(a × b) = a·(b × c) = b·(c × a)
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
101 / 143
Area di una faccia
Qualunque faccia congiungente NV vertici Vl , con l = 1 . . . NV può essere
decomposta in NV − 2 triangoli con un vertice comune.
Vettore area Aj della faccia
Ottenuto dalla somma dei vettori area dei triangoli che la compongono:
NV
1X
Aj =
rVl−1 − rV1 × (rVl − rV1 )
2
l=3
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
102 / 143
Centroide di una faccia
Per il triangolo definito dai vertici V1 , V2 e V3 , la posizione del centroide è data
dalla:
1
r123 = (rV1 + rV2 + rV3 )
3
Centroide rj della faccia
Ottenuto dalla media, pesata per la rispettiva area, dei centroidi di tutti i triangoli
utilizzati per la sua suddivisione:
NV
1X
rV1 + rVl−1 + rVl rVl−1 − rV1 × (rVl − rV1 )
3
rj =
l=3
NV
X
rVl−1 − rV1 × (rVl − rV1 )
l=3
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
103 / 143
15 marzo 2010
104 / 143
Volume e centroide di un tetraedro
Volume di un tetraedro
V1234 =
1
(rV3 − rV1 ) · [(rV2 − rV1 ) × (rV4 − rV1 )]
6
Centroide di un tetraedro
r1234 =
()
1
(rV1 + rV2 + rV3 + rV4 )
4
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
Volume di un poliedro
Volume di un poliedro
Scomposizione in tetraedri e sommando il contributo di ciascuno di questi:
Nf
NV,i
1 XX
V=
rVi,1 − rVref · rVi,l−1 − rVref × rVi,l − rVref
6
i=1 l=3
con Nf numero facce del poliedro, NV,i numero vertici della faccia i-esima, e Vref un
vertice (o qualunque altro punto) di riferimento arbitrario.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
105 / 143
Centroide di un poliedro
Centroide di un poliedro
Dalla media pesata, con il rispettivo volume, dei vettori posizione dei centroidi dei
tetraedri costituenti il poliedro:
rC0
()
Nf NV,i
1 1 XX
=
rVref + rVi,l + rVi,l−1 + rVi,l
24 V
i=1 l=3
rVi,1 − rVref · rVi,l−1 − rVref × rVi,l − rVref
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
106 / 143
Relazioni semplificate
Sfruttando l’identità:
∇r = ∇ (xi + yj + zk) ≡ 3
il volume della cella si può calcolare sfruttando il teorema di Gauss:
Nf
Z
Z
1
1
1X
V=
∇r dV =
r · ndA =
ri · ni Ai
3 V
3 A
3
i=1
Con considerazioni analoghe, attraverso una derivazione un pò più complessa, il
centroide può valutarsi con la relazione:
Nf
X
3
(ri · ni ) ri Ai
rC0 =
i=1
Nf
X
4 ri · ni Ai
i=1
Le espressioni viste sono in generale approssimate, ma risultano esatte per celle
aventi facce triangolari o quadrilatere piane.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
107 / 143
Tipologia di celle
Vista la possibilità di suddividere il dominio in VC caratterizzati da numero
arbitrario di facce, è comune l’adozione, nell’implementazione, di una struttura
dati basata sulle facce, anziché basata sui VC o sui nodi:
I
Agevole considerare celle e combinazioni di celle quali quelle rappresentate in
figura:
I
Per ciascuna faccia è sufficiente conoscere le celle adiacenti C1 e C2 , ed il vettore
d’area A, per poter valutare i flussi convettivo e diffusivo.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
108 / 143
Forme comuni delle celle
Tipologie di celle più comuni - (a) tetraedro; (b) piramide; (c) prisma; (d)
esaedro:
(a)
(b)
(c)
(d)
La valutazione numerica degli integrali, di superficie e di volume, richiede la
conoscenza delle coordinate dei centroidi dei VC, Cj , e delle facce, cj .
È inoltre necessario conoscere il vettore d’area di ogni faccia, Aj , e il volume V
del VC.
Questi dati vengono in genere forniti dai generatori di griglia, altrimenti si può
procedere alla loro valutazione, attraverso le relazioni viste, note le coordinate
dei vertici Vj del CV.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
109 / 143
Distribuzione spaziale delle variabili - 1
Le variabili dipendenti, ed i valori delle proprietà termofisiche, sono definite nei
centroidi dei VC.
Tuttavia, per calcolare gli integrali di superficie, è necessario disporre del loro
valore nei centroidi delle facce.
Si suppone una variazione lineare della variabile φ all’interno della cella:
φ = φC0 + (∇φ)C0 · (r − rC0 )
dove (∇φ)C0 indica il gradiente valutato in C0 , costante su tutto il volume di
controllo.
L’espressione vista fornisce un valore diverso di φ sulla faccia, a seconda di
quale dei due VC questa si consideri appartenente.
Formulazione più complessa che, in analogia al caso Cartesiano, garantisce
perfetta simmetria.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
110 / 143
Distribuzione spaziale delle variabili - 2
Formulazione simmetrica - CDS
1
1
φj =
φC 0 + φC j +
(∇φ)C0 · (rj − rC0 ) + (∇φ)Cj · rj − rCj
2
2
dove il pedice j indica la faccia compresa fra i VC C0 e Cj , adiacenti a tale faccia.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
111 / 143
Integrali di superficie e di volume
L’equazione di trasporto dev’essere integrata su ogni VC.
Adozione di schemi del second’ordine, adeguati per le applicazioni industriali
della CFD.
Utilizzo della formula del punto medio.
Integrazione
Z
φ dV ≈ φC0 ∆V
V
Z
φ dA ≈ φj Aj
Aj
Z
∇φ · dA ≈ (∇φ)j · Aj
Aj
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
112 / 143
Calcolo del gradiente
Per poter calcolare i valori della variabile a centro-faccia, e per valutare termine
diffusivo, è necessario conoscere il gradiente di φ nel centroide della faccia.
Anche il calcolo del termine convettivo necessita della conoscenza del gradiente.
Vi sono diverse strategie, adottate nei programmi commerciali, per ottenere
(ricostruire) tale grandezza sulla base dei soli valori della variabile sui centroidi.
Quelle indicati nel seguito rappresentano due delle alternative più diffuse.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
113 / 143
Metodo di Gauss-Green
φx0 , φy0 e φz0 sono le componenti del gradiente di φ in C0 , e φx0 può essere
considerato la divergenza del vettore φi.
Applicando il teorema di Gauss-Green:
Z
Z
Z
X
φx0 ∆V ≈
φx dV =
∇ · (φ i) dV = φ i · dA ≈
φj Aj,x
V
V
A
j
Procedendo in maniera analoga per le altre componenti del gradiente si ottiene:
φx0 i + φy0 j + φz0 k
P
P
P
φ
A
φ
A
j
j,x
j
j,y
j
j
j φj Aj,z
≈
i+
j+
k
∆V
∆V
∆V
Ricostruzione del gradiente con il metodo di Gauss-Green
P
(∇φ)C0 ≈
()
j
φ j Aj
∆V
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
114 / 143
Metodo di Gauss-Green - esempio - 1
La relazione trovata fornisce un valore del gradiente, in C0 , che non dipende da φC0 .
Esempio 2D
Si ottiene
1
φx0 ≈
∆V
(φC0 + φC1 )
(φC0 + φC2 )
(φC0 + φC3 )
A1,x
A2,x +
A3,x
2
2
2
1
φy0 ≈
∆V
(φC0 + φC1 )
(φC0 + φC2 )
(φC0 + φC3 )
A1,y
A2,y +
A3,y
2
2
2
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
115 / 143
15 marzo 2010
116 / 143
Metodo di Gauss-Green - esempio - 2
Esempio 2D (continua)
Semplificando le espressioni trovate:
()
φx 0 ≈
1
[φC1 A1,x + φC2 A2,x + φC3 A3,x ]
2 ∆V
φy0 ≈
1
[φC1 A1,y + φC2 A2,y + φC3 A3,y ]
2 ∆V
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
Metodo di Gauss-Green - esempio - 3
Adozione di un dominio d’integrazione esteso: riduce gli errori associati alla
valutazione dei valori della variabile sulle facce tramite interpolazione.
Esempio 2D - Dominio d’integrazione esteso
1
[φC1 (yC2 − yC3 ) + φC2 (yC3 − yC1 ) + φC3 (yC1 − yC2 )]
2 ∆V
1
φy0 ≈
[φC1 (xC3 − xC2 ) + φC2 (xC1 − xC3 ) + φC3 (xC2 − xC1 )]
2 ∆V
φx0 ≈
Anche in questo caso, però, le componenti del gradiente in C0 non dipendono da φC0 .
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
117 / 143
Metodo di Gauss-Green - commenti
Il metodo di Gauss-Green, per la valutazione dei gradienti, è accurato solo per
griglie Cartesiane, o comunque strutturate ortogonali.
Piuttosto impreciso per griglie non-uniformi e griglie non strutturate costituite da
celle di tipo diverso, quale è il caso delle griglie ibride.
Migliore accuratezza, con un onere computazionale leggermente maggiore,
ottenuta tramite un’approssimazione ai minimi quadrati.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
118 / 143
Metodo ai minimi quadrati - 1
La relazione vista (non simmetrica) che fornisce la legge di variazione lineare della
generica variabile φ all’interno della cella, può essere scritta in forma esplicita:
φ = φC0 + φx0 (x − xC0 ) + φy0 (y − yC0 ) + φz0 (z − zC0 )
Imponendo che questa relazione sia valida non solo sulla cella in esame, ma descriva
l’andamento spaziale della variabile sino ai centroidi delle celle adiacenti, si ha:
φCj = φC0 + φx0 xCj − xC0 + φy0 yCj − yC0 + φz0 zCj − zC0 − rj
dove rj è il residuo: sistema sovradeterminato.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
119 / 143
15 marzo 2010
120 / 143
Metodo ai minimi quadrati - 2
In forma compatta:
rj = φx0 dxj + φy0 dyj + φz0 dzj − ∆j
dove:
∆j
dxj
dyj
dzj
=
=
=
=
φC j − φC 0
xCj − xC0
yCj − yC0
zCj − zC0
Il quadrato del residuo si traduce quindi nella:
2
2
2
rj2 = ∆2j + (φx0 dxj ) + (φy0 dyj ) + (φz0 dzj )
−2 (∆j φx0 dxj + ∆j φy0 dyj + ∆j φz0 dzj )
+2 (φx0 φy0 dxj dyj + φx0 φz0 dxj dzj + φy0 φz0 dyj dzj )
che sommato per tutte le nb (neighbours) celle adiacenti fornisce:
R=
nb
X
rj2
j=1
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
Metodo ai minimi quadrati - 3
Minimizzando ai minimi quadrati:
∂R
∂φx0
∂R
∂φy0
∂R
∂φz0
= 0
= 0
= 0
Sviluppando, raccogliendo ed esprimendo in forma matriciale:
 nb
 nb

nb
nb
X
X
X
X


2

(dxj )
dxj dyj
dxj dzj 
∆j dxj







j=1
j=1
j=1

 
 j=1

 X


nb
nb
nb
nb
X
X

  φx0   X
2

dxj dyj
(dyj )
dyj dzj 
∆j dyj

  φy0  = 
 j=1


j=1
j=1
j=1
φ z0

 nb



nb
nb
nb
 X


X
X
X

2 



dxj dzj
dyj dzj
(dzj )
∆j dzj


j=1
j=1
()
j=1
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
j=1























15 marzo 2010
121 / 143
Metodo ai minimi quadrati - 4
Con simbologia compatta
D0 ∇φC0 = f0
Ricostruzione del gradiente con il metodo ai minimi quadrati
∇φC0 = D−1
0 f0
D0 è una matrice simmetrica 3 × 3, i cui elementi dipendono solo dalla griglia.
Pertanto è necessario calcolare la matrice inversa D−1
0 , con griglia indeformabile,
solo all’inizio del calcolo.
La ricostruzione del gradiente col metodo dei minimi quadrati richiede quindi di
memorizzare, ed invertire, la matrice D0 per ogni cella, per un totale di 6 × N
numeri reali, con N numero totale di celle.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
122 / 143
Metodo ai minimi quadrati - commenti
Per griglie molto distorte la matrice D0 può risultare mal condizionata:
I
I
Adozione di una metodologia di soluzione basata sulla fattorizzazione QR
(fattorizzazione della matrice di partenza nel prodotto di una matrice ortogonale Q e
una matrice triangolare superiore R tramite il processo di Gram-Schmidt).
Introduzione di pesi inversamente proporzionali alla distanza delle celle vicine.
Il metodo è caratterizzato da un’accuratezza del primo ordine per griglie
costituite da celle di forma arbitraria, ed è inoltre consistente, cioè il gradiente di
una funzione lineare è calcolato esattamente per qualunque tipo di cella. Esso è
perciò particolarmente adatto a trattare griglie ibride.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
123 / 143
Termine transitorio
La procedura, ed il risultato, dipende dallo schema di integrazione temporale
adottato.
Per problemi stazionari o con lente variazioni temporali, è sufficiente l’adozione
dello schema di Eulero implicito del primo ordine.
Per flussi con rapide variazioni temporali, risulta conveniente lo schema di
Crank-Nicolson del second’ordine.
Integrazione del termine transitorio
Z
V
()
∂
∂
∆V n+1
(ρφ) dV ≈
(ρC0 φC0 ∆V) ≈ ρC0
φC0 − φnC0
∂ϑ
∂ϑ
∆ϑ
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
124 / 143
Termine sorgente
Integrazione del termine sorgente
Z
sdV ≈ sC0 ∆V
V
Qualora il termine sorgente dipenda dalla variabile stessa, si procede alla sua
linearizzazione secondo quanto visto nel caso delle griglie Cartesiane.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
125 / 143
Flusso convettivo - 1
Il termine convettivo viene calcolato considerando il contributo di tutte le facce
costituenti il VC:
Z
XZ
ρφw · dA =
ρφw · dA
A
j
Aj
e, sulla base della formula del punto medio:
Flusso convettivo
XZ
j
ρφw · dA ≈
Aj
X
ṁj φj
j
dove ṁj rappresenta la portata di massa attraverso la faccia Aj :
Z
ṁj =
ρw · dA ≈ ρj (wj · Aj )
Aj
Linearizzazione del termine convettivo attraverso il metodo di Picard.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
126 / 143
Flusso convettivo - 2
Valutazione di φj :
Schema CDS
1
1
φC 0 + φC j +
(∇φ)C0 · (rj − rC0 ) + (∇φ)Cj · rj − rCj
φj =
2
2
Instabilità ed oscillazioni (wiggles), e difficoltà di convergenza.
Cura: raffinamento selettivo della griglia.
Schema UDS
(
φj =
φC 0
φCj
se (wj · Aj ) > 0
se (wj · Aj ) < 0
Introduzione di diffusività artificiale: usare, se proprio necessario, con cautela.
Accettabile per equazioni di trasporto source-dominated, es. modelli di
turbolenza.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
127 / 143
15 marzo 2010
128 / 143
Flusso convettivo - 3
Rappresentazione grafica dello schema UDS:
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
Flusso convettivo - 4
Schema upwind di ordine più elevato
(
φj =
φC0 + (∇φ)C0 · (rj − rC0 )
φCj + (∇φ)Cj · rj − rCj
se (wj · Aj ) > 0
se (wj · Aj ) < 0
Ricorso a schemi misti, attraverso un fattore di blending γ compreso fra 0 e 1:
Schema misto
ob
ob
oe
ob
φj = γφoe
+
(1
−
γ)
φ
=
φ
+
γ
φ
−
φ
j
j
j
j
j
γ = 1 corrisponde allo schema CDS, e γ = 0 corrisponde allo schema UDS.
ob indica uno schema di ordine basso (UDS), e oe indica uno schema di ordine
elevato(CDS).
Utilizzo in modalità correzione differita.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
129 / 143
15 marzo 2010
130 / 143
Flusso convettivo - 5
Rappresentazione grafica dello schema upwind di ordine elevato:
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
Flusso diffusivo - 1
Il termine diffusivo viene calcolato considerando il contributo di tutte le facce
costituenti il VC:
Z
XZ
Γ ∇φ · dA =
Γ ∇φ · dA
A
Aj
j
e, sulla base della formula del punto medio::
XZ
X
Γ ∇φ · dA ≈
Γj (∇φ)j · Aj
j
Aj
j
Il valore di (∇φ)j viene approssimato tramite la:
i
1h
(∇φ)j ≈
(∇φ)C0 + (∇φ)Cj
2
che sostituito nella precedente risulta:
i
XZ
1X h
Γ ∇φ · dA ≈
Γj (∇φ)C0 + (∇φ)Cj · Aj
2
A
j
j
j
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
131 / 143
Flusso diffusivo - 2
P
P
Raccogliendo, ed osservando che j Γj Aj ≈ 0, ed in particolare Γ j Aj = 0 nel caso
di proprietà termofisiche costanti, otteniamo:
Z
1X
Γj (∇φ)Cj · Aj
Γ ∇φ · dA ≈
2 j
A
L’espressione ora ricavata presenta due inconvenienti:
1
2
Abbiamo incluso VC non adiacenti (second neighbors) nella discretizzazione.
Sebbene i contributi dei VC adiacenti (first neighbors) siano formalmente
presenti, essi si cancellano a vicenda: soluzioni con oscillazioni non fisiche
(wiggles o checkerboarding).
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
132 / 143
Flusso diffusivo - 3
Discretizzazione del termine diffusivo per la cella C0 : solo i VC ombreggiati
contribuiscono con coefficienti non nulli.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
133 / 143
Flusso diffusivo - 4
Utilizzo di un’espressione modificata per il gradiente, che includa l’effetto delle celle
adiacenti nella discretizzazione, ad esempio:
!
(∇φ)
·
d
φ
−
φ
j Aj
Cj
C0 Aj
j
] = (∇φ) +
(∇φ)
−
j
j
|dj |
|Aj |
|dj |
|Aj |
Con tali modifiche, l’espressione discreta del flusso diffusivo, scritta utilizzando la
correzione differita, diventa:
Flusso diffusivo
(
Z
Γ ∇φ · dA ≈
A
X
j
Γj
φCj − φC0
|dj |
"
#
)
(∇φ)j · dj Aj
|Aj | + (∇φ)j −
· Aj
|dj |
|Aj |
Il primo termine a destra viene trattato in modo implicito, mentre l’espressione in
parentesi quadra, nota con il nome di termine cross-diffusivo, è trattata in modo
esplicito.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
134 / 143
Flusso diffusivo - osservazioni
La procedura vista per la discretizzazione del termine diffusivo costituisce la base di
alcuni codici commerciali, ma presenta alcuni difetti:
L’accuratezza del metodo, e le proprietà di convergenza, si riducono per griglie
lontane dall’ortogonalità: aumenta il contributo del termine cross-diffusivo.
Per griglie molto distorte il termine cross-diffusivo può addirittura diventare
negativo, con conseguenze catastrofiche per variabili definite positive, ad
esempio l’energia cinetica turbolenta.
Le griglie a tetraedri sono meno accurate, con la discretizzazione vista, rispetto a
griglie costituite da soli esaedri. Le prime, tuttavia, sono spesso preferite per
ragioni legate alla facilità di generazione automatica della griglia, ed eventuale
utilizzo di griglie adattive.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
135 / 143
Griglie ibride - 1
L’approccio FVM per griglie non strutturate è molto diffuso ed è adottato, pur
con qualche variante, in numerosi pacchetti CFD.
Esso presenta un inconveniente qualora venga utilizzato con griglie a tetraedri, in
3D, o a triangoli in 2D.
Nel calcolo di flussi turbolenti con modelli RANS o LES, il centroide delle celle
di calcolo disposte sulle pareti solide deve trovarsi ad una distanza y+ , da queste
ultime, compresa in un range opportuno.
Ciò si traduce nella necessità di schiacciare le prime celle in prossimità delle
pareti:
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
136 / 143
Griglie ibride - 2
Con il metodo FEM non vi sarebbero problemi di stabilità poiché tutti i nodi,
indicati nel particolare, darebbero un contributo alla riga della matrice
corrispondente al nodo in esame, indicato con simbolo pieno.
Vi sarebbero, eventualmente, problemi legati alla ridotta accuratezza di elementi
molto allungati).
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
137 / 143
Griglie ibride - 3
Con il metodo FVM descritto la molecola di calcolo risulta molto deformata:
alcune rette congiungenti i centri delle celle formano un angolo rilevante con le
normali alle facce e, nel caso limite di celle molto allungate, possono risultare
quasi parallele alle facce stesse:
I
I
I
Aumenta il contributo del termine diffusivo esplicito, e si riduce il termine implicito:
difficoltà di convergenza.
Aumenta il contributo dei termini cross-diffusivi che possono talvolta diventare
negativi: ridotta accuratezza.
Per le equazioni di quantità di moto, diminuisce l’accuratezza nel calcolo del
gradiente di pressione in prossimità delle pareti solide.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
138 / 143
Griglie ibride - 4
Soluzione: utilizzo di griglie ibride, cioè esaedri e/o prismi (rettangoli in 2D) per
i primi strati di celle, e tetraedri (triangoli in 2D) al centro del dominio.
La molecola di calcolo è molto più regolare, a causa dell’ortogonalità fra le
normali alle facce e le congiungenti i centroidi.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
139 / 143
Griglie ibride - 5
Opportunità (necessità) di utilizzare griglie ibride, anziché griglie a soli tetraedri,
come confermato dalla gran parte dei prodotti commerciali, che ne prevedono la
generazione automatica quale opzione.
Adozione, qualora possibile, di griglie non strutturate o strutturate a soli esaedri:
possibilità offerta da numerosi pacchetti CFD, e prodotti specifici per la
generazione di griglie.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
140 / 143
Condizioni iniziali
Per problemi non stazionari, è necessario fissare i valori delle variabili all’istante
iniziale.
Se interessa solo il comportamento periodico a regime (es. motori a combustione
interna, organi rotanti di macchine ecc.), tali valori possono essere arbitrari, ma è
necessario integrare per un tempo adeguato, in modo da eliminare gli effetti
dovuti alle condizioni iniziali.
Per innescare alcuni fenomeni, quali ad esempio separazioni non stazionarie
degli strati limite, è conveniente perturbare la soluzione.
Per problemi stazionari, l’utilizzo di condizioni iniziali il più possibile prossime
alla soluzione cercata rende più economico il calcolo, e riduce il rischio di
fallimento (divergenza).
Tali condizioni iniziali possono ottenersi risolvendo, ad esempio:
I
I
Problemi semplificati (es. flussi e scambi termici puramente diffusivi)
Interpolazione/estrapolazione della soluzione di problemi più economici (es. griglia
rada, problema bidimensionale).
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
141 / 143
Condizioni al contorno
Le componenti Cartesiane della velocità non sono generalmente allineate con le
facce delle celle sul contorno.
Non si fa uso di celle fantasma, ed in generale gli integrali sulle facce dei VC sul
contorno vanno espressi, come nel caso Cartesiano, in funzione delle condizioni
al contorno e dei valori delle variabili nei centri dei medesimi VC:
I
I
Con condizioni al contorno di Dirichlet, il flusso convettivo è immediatamente
disponibile, mentre il flusso diffusivo richiede di approssimare il gradiente sul
contorno.
Con condizioni al contorno di Neumann, è immediatamente disponibile il valore del
flusso diffusivo, mentre il valore della variabile, necessario nel caso di condizione al
contorno convettiva, viene ricavato in base all’approssimazione (discreta) del
gradiente.
Le condizioni al contorno sulle pareti solide, nei problemi di flusso turbolento,
dipendono dal modello di turbolenza adottato, e dalle modalità di trattamento
della zona di parete.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
142 / 143
Equazione algebrica finale
Dopo aver assemblato, per ogni VC, tutti i termini che appaiono nella equazione
di trasporto integrata sul VC stesso, si ottiene un’equazione algebrica lineare:
Equazione algebrica finale
A0 φ0 +
X
Anb φnb = S0
nb
dove con nb si intende ancora che la sommatoria va estesa ai VC adiacenti, il cui
numero però, a differenza del caso di griglie strutturate, varia da cella a cella.
I coefficienti che appaiono nell’equazione dipendono dalle approssimazioni
usate, e contengono i contributi dei flussi convettivi e diffusivi, del termine non
stazionario, e del termine sorgente.
La matrice risultante, ottenuta scrivendo l’equazione per ogni VC, è sparsa, e
viene di solito risolta con metodi iterativi.
()
Termofluidodinamica Computazionale - Volumi Finiti, E. Nobile
15 marzo 2010
143 / 143