UNIVERSITÀ DEGLI STUDI DI TRIESTE
FACOLTÀ DI INGEGNERIA
Corso di Laurea in Ingegneria Meccanica
Dipartimento di Ingegneria Navale, del Mare e per l’Ambiente
Sezione di Fisica Tecnica
Tesi di laurea in
FISICA TECNICA
PARAMETRIZZAZIONE GEOMETRICA ED
OTTIMIZZAZIONE DI FORMA IN PROBLEMI DI
SCAMBIO TERMICO CONVETTIVO
Laureando:
Gino RIZZETTO
Relatore:
Chiar.mo Prof. Enrico NOBILE
Correlatore:
Dott. Ing. Francesco PINTO
Anno Accademico 2003-2004
A Coloro
che con dolce e ferma pazienza
mi sopportano
e a Nuvola Rapida.
Chiunque insegue un sogno non desidera,
in realtà la sua realizzazione, ma vuole solo
poter continuare a sognare.
All’orizzonte di quell’oceano ci sarebbe stata sempre un’altra isola,
per ripararsi durante un tifone,
o per riposarsi e amare.
Quell’orizzonte aperto sarebbe sempre stato lı́,
un invito ad andare.
Hugo Pratt
Corto Maltese
Una ballata del mare salato
”Non affliggerti troppo nel tuo cuore.
Nessuno riuscirà a uccidermi se non lo vorrà il destino;
e se il destino lo vorrà, allora pensa che al destino nessun uomo,
una volta che è nato, può sfuggire.
Vile o coraggioso che sia.
Nessuno. [...]”
Alessandro Baricco
Omero,
Illiade
Indice
Introduzione
1
2
XI
Utilizzo di NURBS per la rappresentazione geometrica di profili
1.1 Rappresentazione parametrica della geometria in FEMLAB . .
1.2 NURBS 2D . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2.1 Curve 2D non periodiche . . . . . . . . . . . . . . . .
1.2.2 Curve 2D periodiche . . . . . . . . . . . . . . . . . .
1.2.3 Trasformazione NURBS-Bézier 2D . . . . . . . . . .
1.3 Superfici 3D . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3.1 Trasformazione NURBS-Bézier 3D . . . . . . . . . .
Definizione del problema nel caso bidimensionale
2.1 Definizione della geometria . . . . . . . . . . . . . . . . . .
2.1.1 Parametrizzazione del profilo inferiore del canale . .
2.1.2 Generazione del profilo superiore . . . . . . . . . .
2.2 Definizione e risoluzione del problema termofluidodinamico
2.2.1 Grandezze adimensionali e grandezze di riferimento
2.2.2 Equazioni adimensionali . . . . . . . . . . . . . . .
2.2.3 Impostazione delle equazioni per il campo fluido . .
2.2.4 Fattore di attrito adimensionale . . . . . . . . . . .
2.2.5 Tecniche di soluzione per il campo fluidodinamico .
2.2.6 Primo approccio . . . . . . . . . . . . . . . . . . .
2.2.7 Secondo approccio . . . . . . . . . . . . . . . . . .
2.2.8 Terzo approccio . . . . . . . . . . . . . . . . . . . .
2.2.9 Impostazione delle equazioni per il campo termico .
2.2.10 Tecniche di soluzione per il campo termico . . . . .
VII
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
1
2
4
6
8
10
12
13
.
.
.
.
.
.
.
.
.
.
.
.
.
.
17
17
20
23
23
24
26
28
31
32
33
35
37
40
42
VIII
3
4
5
6
INDICE
Caratteristiche della procedura risolutiva in FEMLAB
3.1 Sezione geometrica . . . . . . . . . . . . . . . . . .
3.1.1 Controllo di unicità del dominio . . . . . . .
3.1.2 Controllo di stampabilità del profilo . . . . .
3.2 Sezione numerica . . . . . . . . . . . . . . . . . . .
3.2.1 Parte fluidodinamica. . . . . . . . . . . . . .
3.2.2 Parte termica. . . . . . . . . . . . . . . . . .
3.2.3 Discretizzazione del dominio. . . . . . . . .
.
.
.
.
.
.
.
49
49
49
50
53
54
54
55
Ottimizzazione 2D
4.1 Il processo di ottimizzazione . . . . . . . . . . . . . . . . . . . . .
4.2 Prima fase: ottimizzazione MOGA multi-obiettivo . . . . . . . . .
4.2.1 Primo lancio . . . . . . . . . . . . . . . . . . . . . . . . .
4.2.2 Secondo lancio . . . . . . . . . . . . . . . . . . . . . . . .
4.2.3 Terzo lancio . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3 Seconda fase: eliminazione di un obiettivo . . . . . . . . . . . . . .
4.3.1 MOGA-II: massimizzazione di scambio termico . . . . . .
4.3.2 MOGA-II: minimizzazione del fattore d’attrito . . . . . . .
4.4 Terza fase: ottimizzazione mono-obiettivo mediante utility function .
4.4.1 Utility function con MOGA-II e SIMPLEX . . . . . . . . .
4.4.2 MOGA-II con differente Utility function . . . . . . . . . . .
4.5 Analisi dei Risultati . . . . . . . . . . . . . . . . . . . . . . . . . .
4.6 Confronto geometrie smooth e spigolose . . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
57
57
58
58
62
65
66
66
68
71
72
77
83
89
Definizione geometrica e termofluidodinamica del caso 3D
5.1 Canale 3D estruso . . . . . . . . . . . . . . . . . . . . .
5.2 Canale 3D realizzato con superfici NURBS . . . . . . .
5.3 Definizione del problema termofluidodinamico . . . . .
5.3.1 Condizioni al contorno . . . . . . . . . . . . . .
5.3.2 Cenni alla procedura risolutiva . . . . . . . . . .
5.4 Meshing & Solving . . . . . . . . . . . . . . . . . . . .
93
. 94
. 95
. 99
. 100
. 101
. 102
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
Risultati dell’analisi tridimensionale
103
6.1 Prima strategia di analisi . . . . . . . . . . . . . . . . . . . . . . . . 103
6.2 Seconda strategia di analisi . . . . . . . . . . . . . . . . . . . . . . . 107
Conclusioni
A I Software Utilizzati.
A.1 modeFRONTIER . . . . . . . . . . . . . . .
A.1.1 Caratteristiche essenziali . . . . . . .
A.1.2 Algoritmo Genetico: MOGA-II . . .
A.1.3 Algoritmo mono-obiettivo: SIMPLEX
113
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
117
117
118
120
121
IX
INDICE
A.1.4 Definizione di soluzione migliore
A.1.5 M.C.D.M. . . . . . . . . . . . . .
A.2 FEMLAB . . . . . . . . . . . . . . . . .
A.2.1 Caratteristiche essenziali . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
121
122
123
123
Introduzione
Gli scambiatori di calore, ed in particolare gli scambiatori di calore compatti liquidogas, sono dispositivi molto diffusi in numerosi settori applicativi. Si citano, a titolo
di esempio, il settore automobilistico e dei trasporti ed i settori energetico e chimico.
Nella pratica industriale la configurazione geometrica delle superfici alettate e dei canali degli scambiatori di calore compatti presenta le forme più disparate. Nonostante la
loro semplicità costruttiva e di manutenzione, gli scambiatori a canale piano-parallelo
non sono caratterizzati da prestazioni elevate di scambio termico. Per tale ragione, numerose aziende e diversi autori hanno affrontato lo studio di forme più complesse atte
a incrementare le performance.
Il presente lavoro di tesi rientra in questo filone di analisi, focalizzandone l’aspetto metodologico. L’aspetto immediatamente applicativo è posto in secondo piano,
semplificando la fisica del sistema, senza per questo nulla togliere alla generalità del
metodo. Il problema dello scambio termico in convezione forzata in scambiatori di
calore compatti è sviluppato numericamente in ottica CAE (Computer Aided Engineering) con l’ausilio di strumenti CFD (Computational Fluid Dynamics). In particolare, si adottano un solutore agli elementi finiti (FEM) con griglia non-strutturata e un
applicativo che dispone di diversi algoritmi di ottimizzazione. Seguendo un sentiero
già tracciato da altri autori, si intraprende un processo di ottimizzazione di forma al
fine di migliorare le prestazioni di un canale convettivo periodico, che costituisce il
building block fondamentale degli scambiatori di calore compatti. Come indicatori di
prestazione si individuano lo scambio termico, che va aumentato, le perdite di carico,
che vanno ridotte e alcuni vincoli costruttivi.
L’intero lavoro può essere suddiviso in tre fasi:
1. parametrizzazione geometrica del canale
2. definizione di una metodologia risolutiva per i campi fluidodinamico e termico
all’interno dello scambiatore
3. implementazione e interfacciamento di tale metodologia con sistema di ottimizzazione capace di legare forme e prestazioni
XI
XII
Introduzione
Il lavoro ha inizio con la parametrizzazione della geometria del canale. Un primo elemento di novità è l’utilizzo di curve NURBS (Non-Uniform Rational B-Spline)
per la descrizione del profilo del canale, sia per il caso bidimensionale sia per il caso tridimensionale. Le NURBS consentono una agevole definizione del concetto di
continuità della geometria, ed i loro punti di controllo rappresentano, assieme ad altri
parametri, le variabili di progetto del problema.
Successivamente si descrive a livello teorico la metodologia di risoluzione. Sfruttando la periodicità spaziale del problema, l’attenzione si concentra su di un singolo
modulo ripetitivo del canale, con grande vantaggio per il costo computazionale. Lo
studio è condotto in forma adimensionale. Il fluido utilizzato, in condizioni di incompressibilità, è caratterizzato da un valore del numero di Prandtl pari a 0.7, rappresentativo dell’aria. Si assume fissato il numero di Reynolds, in regime laminare,
stazionario, quindi idrodinamicamente e termicamente sviluppato. Per quel che riguarda il campo termico si introducono condizioni di temperatura imposta alle pareti,
caso poco trattato in letteratura. Lo schema di calcolo prodotto, del tutto generale, è
applicabile e al caso bidimensionale, canali caratterizzati dalla presenza di due sole
componenti di velocità, e al più complesso ed interessante caso tridimensionale.
Dopo la descrizione della messa a punto di tale metodo, si riporta la fase di ottimizzazione bidimensionale, dove gli input sono i parametri geometrici, mentre gli
output sono le prestazioni. Poiché si è in presenza di un problema multi-obiettivo la
definizione di ottimo introduce al concetto di fronte di Pareto e dominanza di Pareto.
Al termine della fase di ottimizzazione si procede a un confronto fra geometrie con
parametrizzazione NURBS, smooth profile, e geometrie parametrizzate mediante tratti rettilinei, linear-piecewise profile. L’adozione di profili smooth permette, a parità di
scambio termico, una notevole riduzione del coefficiente d’attrito.
Quanto sviluppato finora è propedeutico per l’ottimizzazione di forma in canali tridimensionali. L’obiettivo diviene la verifica della possibilità di incremento di scambio
termico rispetto al caso bidimensionale grazie alla presenza di eventuali moti secondari. Il costo computazionale cresce notevolmente, soprattutto dal punto di vista della
soluzione numerica, oltre che per il maggior numero di gradi di libertà geometrici. La
bontà della metodologia è inizialmente saggiata mediante un’analisi parametrica, che
risulta fedele alle aspettative. Si analizza cioè un numero ristretto di canali tridimensionali, al variare di un unico parametro, sensibile alla generazione dei moti secondari.
Si affronta infine il problema di ottimizzazione tridimensionale vera e propria, che, sia
pur a livello qualitativo, rappresenta un aggiuntivo fattore di innovazione.
Volutamente i software utilizzati sono di tipo commerciale e non dedicato. Questo
per garantire generalità alla procedura. La parametrizzazione geometrica e l’analisi
CFD sono condotte mediante FEMLAB, pacchetto agli elementi finiti che può essere utilizzato anche come toolbox MATLAB. Il sistema di ottimizzazione presente in
questa tesi è modeFRONTIER, il quale raggruppa in forma organica e accessibile una
larga serie di algoritmi e strumenti per l’ottimizzazione.
Quanto realizzato si dimostra un procedimento robusto e sufficientemente effi-
XIII
ciente. I risultati sono incoraggianti ed aprono la strada a successive investigazioni di
sistemi con geometria e fisica più complesse, i.e. moto non stazionario o turbolento.
Gli unici limiti riguardano ancora l’onere di calcolo e l’accuratezza degli eventuali
modelli che si dovrebbero introdurre.
1
Utilizzo di NURBS per la
rappresentazione geometrica di
profili
Buona parte degli elementi geometrici trattati in questo lavoro appartengono alla categoria delle cosiddette forme non-definibili. Sono forme complesse, capaci di rispondere alle più svariate richieste progettuali, come continuità, smoothness e definizione
agevole [1]. A queste si contrappongono le forme definibili, ovvero quelle descritte
dalla geometria classica, incapaci di ottemperare in modo semplice alle odierne necessità progettuali. Si immagini ad esempio di dover disegnare il profilo di un’autovettura
utilizzando unicamente coniche.
La rappresentazione matematica di una curva (per una superficie si tratta di una
semplice estensione) può essere:
1. non-parametrica, cioè y = f (x), forma esplicita, oppure f (x, y) = 0, forma implicita, la quale consente la descrizione di curve chiuse;
2. parametrica, in cui la curva C è espressa in forma vettoriale al variare del
parametro t come in (1.1)
C(t) = [x(t), y(t)]
(1.1)
con t ∈ [a, b], dove a e b sono numeri reali non-negativi arbitrari. Il parametro t
può essere visto come la coordinata curvilinea che percorre la curva.
La rappresentazione parametrica e polinomiale presenta diverse proprietà:
- la forma della curva risulta intrinsecamente indipendente dal sistema di coordinate prescelto: soltanto le relazioni che intercorrono tra una particolare serie di
punti determina la forma della curva;
1
2
Utilizzo di NURBS per la rappresentazione geometrica di profili
- le forme non-parametriche non consentono di descrivere curve che presentano pendenze infinite rispetto un qualsiasi sistema di riferimento, mentre per le
forme parametriche vale
dy/ dt
dy
=
dx
dx/ dt
per cui quando il termine a sinistra è infinito, nel termine a destra è semplicemente dx/ dt = 0;
- gli elementi geometrici sono facilmente trattabili sotto forma di vettori e matrici,
il che rende più naturale la loro rappresentazione al computer;
- in genere il numero di gradi di libertà concessi per il controllo della curva è
maggiore rispetto al caso non-parametrico a parità di ordine della curva, ovvero
l’ordine del polinomio che la descrive.
Si è appena fatto riferimento all’ordine della curva polinomiale. Il grado della curva,
detto p l’ordine, è pari a p −1 e si identifica con il grado del polinomio che la descrive.
Il grado di un polinomio e’ uguale al grado del suo monomio di grado più alto.
1.1
Rappresentazione parametrica della geometria in
FEMLAB
I problemi di ottimizzazione svolti sono incentrati sull’utilizzo di FEMLAB. La definizione geometrica del dominio dei problemi trattati passa attraverso il sistema CAD
interno al software. La versione in uso di FEMLAB permette la creazione di geometrie complesse in 2D e in 3D con l’utilizzo di curve e superfici di Bézier del terzo e
quarto ordine.
Le rappresentazioni parametriche secondo curve e superfici di Bézier consentono
di approssimare [1] un insieme di punti di controllo senza il vincolo del passaggio
- interpolazione - per i punti considerati, tipica delle curve spline, di largo impiego
in campo navale. Con la tecnica dell’approssimazione una eventuale modificazione
della geometria è legata alla variazione di pochi parametri e la stessa definizione della
geometria richiede un numero limitato di punti. Inoltre è agevole garantire la continuità della curva. Facendo riferimento al caso 2D, una curva di Bézier gode di alcune
principali proprietà [1, 2]:
1. la curva è definita da un poligono di controllo, il cui numero di vertici è pari
all’ordine della curva;
2. la curva ha gli estremi in comune con il poligono e i vettori ivi tangenti sono
paralleli al primo e all’ultimo segmento del poligono;
3. la curva è contenuta nell’involucro convesso delimitato dal poligono di controllo, di cui segue l’andamento spaziale.
3
1.1 Rappresentazione parametrica della geometria in FEMLAB
La rappresentazione alla Bézier di una curva è una delle più semplici all’interno della
famiglia delle funzioni parametriche e polinomiali. La generica curva di ordine p (il
grado della curva è p − 1) C = C(t) viene descritta come in (1.2)
n
t ∈ [0, 1]
C(t) = ∑ Bi,p (t)Pi
(1.2)
i=1
dove Pi sono i punti di controllo, la cui unione mediante spezzate dà origine al poligono di controllo tratteggiato in figura 1.1(a). Bi,p (t) sono le funzioni di base, figura
1.1(b).
1
0.5
P2
0.4
0.9
0.8
0.3
0.7
0.2
0.6
0.1
0.5
y
0
P4
P0
−0.1
0.4
0.3
−0.2
0.2
−0.3
P3
0.1
−0.4
0
0.2
0.4
x
0.6
0.8
(a) Esempio di curva di Bézier.
1
0
0
0.2
0.4
t
0.6
0.8
1
(b) Esempio di funzioni di base IV ordine.
Figura 1.1 Curva di Bézier, poligono di controllo e funzioni di miscelamento. L’ordine delle
polinomiali è 4.
Secondo la trattazione di Bézier le funzioni di base - o di miscelamento - che
garantiscano le proprietà sopra citate [1, 2] sono le cosiddette polinomiali di Bernstein.
La generica funzione di base per una curva di Bézier di ordine p, riferita al punto i è
Bi,p (t) =
p!
ui (1 − u) p−i
i!(p − i)!
u ∈ [0; 1]
(1.3)
La modellazione geometrica alla Bézier non consente tuttavia di controllare localmente la forma delle curve realizzate: una modifica locale (i.e. spostamento di un
punto di controllo) causa una perturbazione che si manifesta in modo globale sull’intera entità geometrica. Questo effetto è dovuto alla tipologia delle particolari funzioni di
miscelamento utilizzate dalle Bézier, non nulle lungo l’intero intervallo di definizione
del parametro, figura 1.1(b). L’aumento del numero di punti di controllo è accompagnato inoltre da un pari aumento dell’ordine polinomiale della curva, cosa che rende
più pesante la trattazione numerica.
4
Utilizzo di NURBS per la rappresentazione geometrica di profili
Questi limiti sono superati con la modellazione mediante NURBS, qui utilizzate
nella forma di B-Spline.
1.2
NURBS 2D
Le NURBS (Non-Uniform Rational B-Spline) [1, 3] sono una estensione delle BSpline. Queste ultime sono curve polinomiali a tratti, ossia non viene utilizzato un unico polinomio per la rappresentazione di forme complesse che necessitino un numero
elevato di punti di controllo.
Le funzioni di miscelamento delle B-Spline, a differenza di quelle delle Bézier,
sono non-nulle solo in un tratto limitato del dominio del parametro, per cui non accendono i punti di controllo o danno loro peso lungo l’intero dominio: una eventuale
modifica a un punto di controllo ha effetto solo nell’intorno del punto stesso, local
support. Per questa e ulteriori proprietà si vedano [1, 2, 3]. Le B-Spline non consentono però la definizione di coniche, per le quali sono necessarie forme razionali;
entrano quindi in gioco le NURBS, ossia curve o superfici B-Spline razionali con
distribuzione nodale non-uniforme.
La generica formulazione di una curva NURBS è analoga a quella di una curva di
Bézier. L’unica differenza appare nella tipologia delle funzioni di miscelamento. Per
una NURBS di ordine p, grado p − 1 e un numero n di punti di controllo
n
C(t) = ∑ Ri,p (t)Pi
(1.4)
i=1
dove
- il parametro t ∈ [0, 1];
- C(t) rappresenta i punti della curva al variare del parametro t;
- Pi sono i punti di controllo;
- Ri,p (t) sono le funzioni base che miscelano i punti di controllo.
Le funzioni base derivano dalle funzioni base polinomiali utilizzate per le B-Spline,
per cui detti w il vettore dei pesi e Ni,p (t) la funzione base B-Spline
Ri,p (t) =
Ni,p (t)wi
n
(1.5)
∑ N j,p (t)w j
j=1
Per le curve NURBS, cosı̀ come per le B-Spline, l’ordine p delle funzioni polinomiali
non ha quel forte legame con il numero dei punti di controllo com’è per le Bézier. Non
è infatti più necessario che il numero dei punti di controllo eguagli l’ordine (n = p),
bensı̀ compare un legame fra numero di punti, ordine del polinomio e numero dei nodi.
5
1.2 NURBS 2D
I nodi, knots, sono dei valori, normalizzati a 1 come il parametro t, i quali definiscono il tratto in cui sono non-nulle le funzioni di base; essi entrano dunque in
relazione con il parametro t e sono caratterizzati da una distribuzione non uniforme,
da cui il termine Non-Uniform. Nella definizione di una NURBS il vettore dei nodi U
è costruito in modo tale che:
U = [0, ..., 0, u p+1 , ..., um−p−1 , 1, ..., 1]
| {z }
| {z }
p
p
Il legame che intercorre fra numero di punti n, ordine della curva p e numero di knots
m è dato dalla (1.6).
m = n+ p
(1.6)
É dunque possibile giocare, fissato che sia l’ordine, tra punti e nodi in modo da ottenere la geometria che meglio soddisfi l’utente. Tra le proprietà notevoli [1, 3] si ricorda
che:
- se n = p e U = [0, ..., 0, 1, ..., 1] la NURBS degenera in una Bézier;
| {z } | {z }
p
p
- è comunque garantito il passaggio per i punti estremi del poligono di controllo;
- i pesi avvicinano o meno la curva al punto di controllo cui sono associati;
- local support, cioè modifiche al punto di controllo Pi o al peso wi hanno influenza solo per t ∈ [ui , ui+p [;
- la continuità della curva e delle sue derivate è garantita all’interno dei sottointervalli definiti dai nodi, poiché le funzioni di base sono polinomi (funzioni
C∞ ), mentre agli estremi dell’intervallo la NURBS è C p−k , dove k è la molteplicità del nodo.
Dall’ultima si ottiene che accrescere la molteplicità del nodo, nel rispetto della relazione sopra ricordata, porta alla riduzione del grado di continuità della curva stessa.
Le funzioni di miscelamento B-Spline Ni,p sono ricavate secondo l’algoritmo di
Cox-De Boor (1.7), generalizzazione sempre in forma ricorsiva dell’algoritmo di de
Casteljau [1].
1 se ui ≤ t < ui+1
Ni,1 =
(1.7a)
0 altrimenti
Ni,k =
(t − ui )Ni,k−1 (t) (ui+k − t)Ni+1,k−1 (t)
+
ti+k−1 − ti
ti+k − ti+1
(1.7b)
6
1.2.1
Utilizzo di NURBS per la rappresentazione geometrica di profili
Curve 2D non periodiche
Il processo di definizione della geometria è partito dal caso bidimensionale, utilizzando gli esempi nei testi citati in bibliografia come test-case. Il software con cui si sono
sviluppati gli algoritmi è MATLAB, cui FEMLAB è interfacciato.
Il primo passo ha visto la generazione delle curve di miscelamento B-Spline secondo l’algoritmo di Cox-De Boor. In questo modo il risultato ottenuto sarebbe stato
di validità generale per la creazione di qualsiasi tipo di curva parametrica, ovvero curve di Bézier, B-Spline, NURBS. L’algoritmo in letteratura viene presentato in forma
ricorsiva, mentre in questo lavoro si è utilizzata la forma iterativa, partendo cioè dalle
funzioni costanti e risalendo fino alle funzioni dell’ordine voluto; è stato inoltre necessario modificare gli indici in (1.7a) e (1.7b) poiché MATLAB permette di utilizzare
come contatori per array solo numeri interi positivi per cui:
Designazione
i
k
Ruolo
contatore punto
contatore livello
Campo di variazione
1:numero di punti
1:ordine della curva
Essendo nullo il valore delle funzioni di miscelamento B-Spline al di fuori dell’intervallo di variazione del parametro, per l’ultimo punto viene eliminato il secondo
addendo della (1.7b).
Nell’intero lavoro si è cercato di fare uso prevalentemente di funzioni MATLAB
(function), in modo da lasciare il più possibile libera la memoria, estraendo come
variabili di output soltanto quelle strettamente necessarie.
Il risultato finale di questa fase è una funzione che, noti i parametri di definizione
di una curva, calcoli le coordinate del punto della curva corrispondente ad un certo
valore del parametro t, con t ∈ [0; 1]. La successiva tabella 1.1 raccoglie gli elementi di
definizione della generica NURBS, B-Spline. I valori sono raccolti nello struct-array
data.
data.point
data.sizeP
data.ord
data.knots
data.sizeK
data.pesi
=
=
=
=
=
=
punti di controllo
numero di punti di controllo
ordine della curva
vettore dei nodi
numero di nodi
vettore dei pesi
Tabella 1.1 Elementi di definizione di una curva.
In questo lavoro si sono utilizzati pesi unitari, per cui le NURBS utilizzate degenerano in B-Spline. La presenza del vettore dei pesi consente alla trattazione un carattere
più generale.
7
1.2 NURBS 2D
La figura 1.2 presenta un esempio di curva NURBS non-periodica con distribuzione dei nodi uniforme e pesi unitari. L’ordine della curva è 3, il numero di punti di
controllo è pari a 7, il numero dei nodi è dunque 10. I dati caratterizzanti la curva sono
riportati in tabella 1.2. In figura 1.3 sono riportate le funzioni di miscelamento che
Punti di controllo
X
Y
0
1.5
3.5
3
6
6.5
7
P1
P2
P3
P4
P5
P6
P7
Nodi
0
0
0
0.2
0.4
0.6
0.8
1
1
1
u1
u2
u3
u4
u5
u6
u7
u8
u9
u10
0
2.5
4
1.5
-1.5
2
1
Tabella 1.2 Parametri di definizione della curva di figura 1.2.
P3
4
3
P2
P6
y
2
P4
1
P7
P1
0
−1
−2
P5
0
1
2
3
x
4
5
6
7
Figura 1.2 NURBS di ordine 3 creata con 7 punti di controllo. Tratteggiato il poligono di
controllo.
8
Utilizzo di NURBS per la rappresentazione geometrica di profili
caratterizzano una qualsiasi curva NURBS di ordine 3, a 7 punti di controllo e distribuzione nodale uniforme. Dalle figura 1.2 si nota la caratteristica interpolazione degli
1
R
R1,3
0.9
7,3
R3,3
0.8
R2,3
funzioni base
0.7
R
4,3
R5,3
R
6,3
0.6
0.5
0.4
0.3
0.2
0.1
0
0
0.2
0.4
0.6
parametro
0.8
1
Figura 1.3 Funzioni base per una NURBS con vettore dei pesi unitari, 7 punti di controllo e
ordine 3.
estremi, ottenuta mediante la ripetizione dei nodi. Questo porta alla condensazione
delle funzioni di miscelamento, come la R1,3 , la quale è non-nulla in un tratto limitato
rispetto alla R3,3 (figura 1.3).
Dalla figura 1.4 è interessante notare come muovendo l’i-esimo punto di controllo,
la curva venga modificata solo nell’intervallo [ui ; ui+p [. Ciascun punto di controllo
influenza una porzione limitata di curva. L’intervallo [0; 1] di variazione del parametro
t descrive l’intera curva, mentre il vettore dei nodi è una suddivisione in tratti di questo
intervallo: ciascuna funzione di miscelamento è non-nulla su di un numero di tratti non
superiore all’ordine della curva stessa. Ad esempio, la funzione R4,3 è non-nulla per
t ∈ [0.2; 0.8].
1.2.2
Curve 2D periodiche
Come poi sarà sottolineato in sezione 2.1, in questo lavoro si è reso necessario l’utilizzo di curve periodiche, in modo da poter agevolmente traslare il profilo superiore
rispetto a quello inferiore del canale. La descrizione mediante NURBS consente l’utilizzo di curve periodiche, in questo caso a distribuzione nodale uniforme e vettore dei
pesi unitario, ovvero B-Spline. Rispetto al caso non-periodico il vettore degli knots
non presenta la ripetizione dei valori estremi per un numero di volte pari all’ordine,
9
1.2 NURBS 2D
P3
4
3
P2
P6
y
2
P4
1
P7
P1
0
P’5
−1
−2
P5
0
1
2
3
x
4
5
6
7
Figura 1.4 Esempio di local support: una modifica al punto di controllo P5 ha influenzato un
tratto limitato della curva di figura 1.2.
bensı̀:
U = [0, ..., 1]
dove la distribuzione tra 0 e 1 è uniforme. Questo porta alla impossibilità di interpolare
i punti di controllo estremi. Le funzioni di base sono tutte uguali fra loro e non vi è
condensazione alle estremità, ovvero non vi sono intervalli degeneri di ampiezza nulla.
Il dominio di variazione del parametro t è quindi ristretto:
t ∈ [u p−1 ; un+1 ]
dove
- u p−1 è il (p − 1) − esimo nodo, con p ordine della curva;
- un+1 è l’(n + 1) − esimo nodo, con n numero di punti di controllo della curva.
In figura 1.5 è proposta una curva periodica con il medesimo poligono di controllo di
figura 1.2. Il vettore dei nodi è uniformemente spaziato, tutti i nodi sono a molteplicità
singola e il poligono di controllo si ripete periodicamente. Confrontando figura 1.3
e figura 1.5(a) si nota l’effetto della mancata condensazione delle funzioni base alle
estremità di curve periodiche. In figura 1.5(b) si osserva la mancanza di interpolazione
agli estremi del poligono di controllo. Questo aspetto è più marcato sul lato destro
della curva, cioè per valori prossimi a 1 del parametro t, a causa della modalità di
implementazione dell’algoritmo (1.7).
10
Utilizzo di NURBS per la rappresentazione geometrica di profili
1
10
0.9
8
0.8
6
0.7
4
0.6
y
0.5
0
0.4
−2
0.3
−4
0.2
−6
0.1
0
2
−8
0.35
0.4
0.45
0.5
t
0.55
0.6
0
0.65
(a) Funzioni di miscelamento: selezione. Le curve
curve sono tutte simili fra loro.
5
(b) Curva periodica:
elementare.
10
x
15
20
è evidenziato il modulo
Figura 1.5 Esempio di curva periodica con vettore dei nodi a distribuzione uniforme.
1.2.3
Trasformazione NURBS-Bézier 2D
Il software FEMLAB tratta le curve sotto forma di Bézier del terzo o del quarto ordine. Una efficiente procedura di ottimizzazione richiede modifiche alla geometria
di tipo locale, consentite da curve polinomiali a tratti quali B-Spline e Non-UniformRational B-Spline (NURBS) [3]. Le NURBS, essendo forme parametriche razionali,
permettono la definizione corretta delle coniche, tuttavia in questo lavoro i pesi che
compaiono nella (1.5) sono assunti unitari e la distribuzione dei nodi è uniforme, per
cui le NURBS ottenute sono in realtà B-Spline(1) .
Per non rinunciare ai vantaggi offerti dall’uso di curve NURBS, si è creata una
procedura che determina i punti di controllo di un opportuno numero di curve di Bézier
dello stesso ordine, tali da riprodurne fedelmente l’andamento. In questo lavoro si è
deciso di utilizzare NURBS quarto ordine e di ridisegnarle con Bézier di pari ordine.
Detta C(u) per la generica NURBS vale la (1.4), mentre per una generica Bézier è
4
CB (t) =
∑ Bi,4 (t)Qi
(1.8)
k=1
dove
- Bi,4 è la funzione di miscelamento delle Bézier, il secondo indice è relativo
all’ordine;
(1) Per
la cosiddetta partition of unity [3], il denominatore della (1.5) assume valore unitario.
11
1.2 NURBS 2D
- Qi è uno dei 4 punti di controllo necessari.
Il numero di curve di Bézier necessarie dipende dalla distribuzione dei nodi della
NURBS; le NURBS sono infatti ottenute con tante funzioni polinomiali a tratti quanti sono gli intervalli definiti dai nodi, mentre le curve di Bézier sono ottenute con
un’unica funzione polinomiale, in questo particolare caso di grado 3. Sono dunque
necessarie tante curve di Bézier quante sono le polinomiali a tratti, detto numint il
numero di intervalli sono richiesti (numint · 4) punti di controllo.
Una curva di Bézier IV ordine è univocamente determinata una volta definiti i suoi
punti di controllo, il calcolo di questi punti si può cosı̀ riassumere:
1. si rigenera il vettore dei nodi della NURBS in modo che tutti i valori abbiano
molteplicità singola;
2. si fanno corrispondere ai valori del parametro u, con u ∈ [Uk−1 ,Uk ], i valori del
parametro t con t ∈ [0, 1], detto k il nuovo contatore dei nodi;
3. i punti di controllo agli estremi della curva di Bézier sono Q(0) = C(Uk−1 ) e
Q(1) = C(Uk );
4. si fissano due valori del parametro u distinti e tali che u ∈]Uk−1 ,Uk [ in cui
calcolare due punti della NURBS: siano essi u1 e u2 cui corrispondono t1 e
t2 ;
5. si risolve il seguente sistema lineare in cui si impone il passaggio della Bézier
per i 4 punti scelti della NURBS:

1
 B1,4 (t1 )

 B1,4 (t2 )
0
0
B2,4 (t1 )
B2,4 (t2 )
0
0
B3,4 (t1 )
B3,4 (t2 )
0

Q1
0
 Q2
B4,4 (t1 ) 

B4,4 (t2 )   Q3
Q4
1

C(Uk−1 )

  C(u1 )

=

  C(u2 )
C(Uk )


6. si avanza di 1 il contatore k e si ritorna a 2. finché non si esauriscono gli
intervalli.
Il risultato di questa operazione su una curva non-periodica a 7 punti di controllo e
distribuzione uniforme dei nodi è presentato in figura 1.6.
12
Utilizzo di NURBS per la rappresentazione geometrica di profili
8
7
6
y
5
4
3
2
1
1
2
3
4
x
5
6
7
8
Figura 1.6 Esempio di trasformazione NURBS-Bézier.
Seguendo il verso indicato dalla freccia, la curva tratteggiata va a sovrapporsi alla
curva di partenza Osservando figura la 1.6, è evidente il moltiplicarsi dei punti di
controllo necessari alla descrizione della curva di partenza mediante curve di Bézier.
1.3
Superfici 3D
Le superfici NURBS nascono come estensione nello spazio tridimensionale delle curve NURBS 2D. Compare in questo caso un poliedro caratteristico che la superficie
va ad approssimare, con qualità differente a seconda dell’ordine dei polinomi che
definiscono le funzioni di miscelamento [1].
Si osserva rispetto alle curve bidimensionali la duplicazione degli elementi caratteristici:
- si definiscono due insiemi di valori per i nodi;
- due sono i parametri u e v, normalizzati a 1, che percorrono la superficie;
- i punti di controllo vanno a formare un array bidimensionale di dimensioni
n x m, dove m è il numero di punti di controllo associati al parametro u e n è
il numero di punti di controllo associati al parametro w;
- per ciascuno dei due insiemi di nodi vale la relazione (1.6) che li lega all’ordine
e al numero di punti di controllo.
13
1.3 Superfici 3D
Si riportano le equazioni che definiscono i punti di una superficie NURBS al variare
dei due parametri (u, v) e le funzioni di base polinomiali razionali
n
S(u, v) =
m
∑ ∑ Rk,h (u, v)Pk,h
(1.9)
k=1 h=1
Rk,h (u, v) =
Nk,p (u)Nh,q (v)wk,h
n
∑
m
(1.10)
∑ Nk,p (u)Nh,q (v)wk,h
k=1 h=1
dove:
- i pedici p e q rappresentano l’ordine delle funzioni di miscelamento Nk,p e Nh,q ;
- w è la matrice dei pesi.
Come per il caso bidimensionale, in questo lavoro si sono utilizzate superfici NURBS
con matrice dei pesi unitaria e distribuzione nodale uniforme, per cui le superfici
degenerano in superfici B-Spline.
Si vuole far notare come le dimensioni dei vettori dei nodi, cosı̀ come dei punti
di controllo possono essere tra loro differenti, cosicché la superficie risultante può
avere ordine differente lungo il senso di percorrenza dei 2 parametri. Le superfici di
Bézier sono ancora un caso particolare delle superfici NURBS, quando, oltre alle già
citate matrice dei pesi unitaria e nodi uniformemente distribuiti, la matrice dei punti di
controllo ha dimensioni pari rispettivamente all’ordine della superficie lungo ciascuna
direzione.
È usuale rappresentare una superficie proiettando su di un piano i suoi punti di
controllo. In questo modo si può identificare la matrice di punti con questa proiezione
o patch. In figura 1.7 si dà un esempio di NURBS 3D definita da una matrice di punti
di controllo 4 x 4, ordine 3 lungo la direzione u (in figura lungo l’asse y) e ordine 2
lungo la direzione w (in figura lungo l’asse x). I vettori degli knots non presentano una
distribuzione uniforme, mentre la matrice dei pesi è unitaria. Come lecito attendersi,
la superficie mostra un andamento lineare a tratti lungo la direzione w (direzione asse
x), mentre è quadratica lungo la direzione u (direzione asse y). Sul piano z = −1
è riportata, come rappresentazione grafica della matrice dei punti, la proiezione in
pianta del poligono di controllo.
Per le superfici periodiche vale l’estensione al caso 3D di quanto detto riguardo le
curve. Come si possono presentare ordini differenti lungo le due direzioni dei parametri, cosı̀ è possibile imporre la periodicità soltanto lungo una direzione di uno dei
due parametri. Ne è un esempio la figura 1.8: lungo u (direzione asse y) la superficie
è periodica e non interpola gli estremi.
1.3.1
Trasformazione NURBS-Bézier 3D
Il procedimento è analogo al caso delle curve in due dimensioni. Si tratta infatti di
suddividere la superficie NURBS in più patch di Bézier, in numero pari al prodotto
14
Utilizzo di NURBS per la rappresentazione geometrica di profili
tra gli intervalli con cui i vettori dei nodi suddividono i due domini di variazione dei
parametri u e w. Si ricorda che le curve di Bézier non sono polinomiali a tratti e che
loro funzioni di miscelamento sono non nulle lungo l’intero intervallo di definizione.
Poichè l’ordine di una NURBS determina il numero massimo di funzioni di base afferenti ad un intervallo nodale, l’utilizzo di curve di Bézier di pari ordine deve essere
limitato a ciascuno di quegli intervalli.
Tuttavia, anche se concettualmente semplice, questo procedimento per superfici
ragionevolmente complesse dimostra di essere oneroso dal punto di vista computazionale. Il numero di patch necessarie, rispetto ad un analogo caso bidimensionale è
decisamente maggiore. Per caso analogo si intende un canale tridimensionale ottenuto per estrusione di un generico canale bidimensionale. Un altro aspetto rilevante è
garantire che due superfici affiancate presentino un bordo in comune; questo fatto può
inoltre essere sfruttato per limitare il numero di punti di controllo da calcolarsi.
Nelle figure 1.9 e 1.10 si riporta in due differenti angolazioni un esempio di trasformazione di superficie NURBS in più superfici di Bézier. La NURBS di partenza
è di tipo doppio-periodico con matrice dei punti di controllo 10 x 4 e ordine 4. Due
delle superfici di Bézier che la ricoprono sono evidenziate con wireframe giallo.
15
1.3 Superfici 3D
Figura 1.7 Esempio di superficie NURBS.
Figura 1.8 Esempio di superficie NURBS periodica lungo la direzione y.
16
Utilizzo di NURBS per la rappresentazione geometrica di profili
Figura 1.9 Esempio di superficie doppio-periodica su cui si è operata la trasformazione
NURBS-Bézier.
Figura 1.10 In questa vista si osservano due delle superfici di Bézier che coprono la NURBS
(wireframe giallo).
2
Definizione del problema nel
caso bidimensionale
Nel presente capitolo si illustrano le tecniche numeriche che si intende utilizzare nell’ottimizzazione di un canale bidimensionale periodico, rappresentativo del modulo
fondamentale di un’ampia classe di scambiatori di calore. Alcuni esempi reali di tali dispositivi sono presentati in figura 2.1. L’obiettivo è l’incremento dello scambio
termico alle pareti del canale, in riferimento al caso del canale piano, cercando di contenere il più possibile l’aumento delle perdite di carico, indotto da forme geometriche
complesse. Questo progetto è strutturato in due fasi:
a. definizione della geometria del canale;
b. definizione della procedura iterativa di soluzione per il campo fluidodinamico e
per il campo termico, con i quali si descrive il regime di convezione forzata.
2.1
Definizione della geometria
Si assume che lo scambiatore di calore sia di profondità infinita, si possono cosı̀ in
pratica trascurare gli effetti di bordo e forzare la presenza di due sole componenti di
velocità: lo studio viene effettuato in forma bidimensionale. Si vuole inoltre ottenere la soluzione termofluidodinamica del canale in condizioni di moto completamente
sviluppato, ragion per cui le geometrie prese in considerazione sono di tipo periodico.
In questo modo, a tutto vantaggio dell’onere computazionale, risulta possibile concentrare l’ottimizzazione su un singolo modulo dell’intera geometria, a patto di porre
adeguate condizioni al contorno di periodicità alle sezioni di ingresso e di uscita del
flusso [4]. Un vincolo sulle geometrie generate riguarda la possibilità di poter impilare i singoli moduli, cosı̀ da simulare al meglio le caratteristiche di uno scambiatore
alettato o a piastre [5], realizzato dall’impaccamento di più piastre metalliche, separate da un’intercapedine entro cui scorre il fluido, o i fluidi, di processo. Per avere
un controllo locale e prevedibile della forma del profilo del canale, la superficie del
17
18
Definizione del problema nel caso bidimensionale
(a) Complessivo di scambiatore di calore a
piastre.
(b) Particolare di scambiatore di calore a
piastre.
(c) Esempi di superfici di scambiatori di calore.
Figura 2.1 Esempi di scambiatori di calore.
2.1 Definizione della geometria
19
Figura 2.2 Esempio di canale 2D.
mantello è modellata mediante una B-Spline, ovvero una NURBS con vettore di pesi
unitari(1) . Infatti in queste particolari espressioni di curve parametriche la modifica di
un punto di controllo perturba una porzione limitata dell’intera curva.
La generazione del canale è ottenuta con i seguenti passi:
1. generazione di una NURBS periodica per il profilo inferiore di lunghezza pari
a 3 volte il periodo L;
2. trasformazione della forma NURBS in più forme Bézier equivalenti; questo
passaggio è necessario poiché FEMLAB tratta curve in forma di Bézier del III
o del IV ordine;
3. traslazione dei punti di controllo lungo l’asse x di una quantità trasl ∈] − L, +L[,
dove L è il periodo della curva;
4. traslazione dei punti di controllo appena trovati lungo l’asse y di una quantità
pari a Hav : si è ottenuto il profilo superiore con un’altezza media del canale pari
a Hav ;
5. il canale cosı̀ ottenuto viene tagliato dalle sezioni di ingresso e uscita.
Il risultato è visibile in figura 2.2.
(1) D’ora in avanti si utilizzeranno in modo ambivalente il termine B-Spline e NURBS, intendendo che
queste ultime abbiano un vettore dei pesi unitari.
20
Definizione del problema nel caso bidimensionale
Per rendere le modifiche alla geometria lungo il percorso di ottimizzazione intuitive e facilmente controllabili si è data una rappresentazione parte cartesiana e parte
polare dei punti di controllo del profilo base del canale.
Si è deciso per un numero di punti di controllo pari a 9, come compromesso fra una
buona risoluzione e trattabilità della curva: un numero eccessivo di punti di controllo
equivarrebbe ad un numero elevato di variabili indipendenti, cosa che appesantisce la
fase di ottimizzazione in FRONTIER. Questi punti vengono, come al passo 1., traslati
in direzione positiva e negativa del periodo L, eliminando i punti doppi nelle giunzioni
figura 2.3: in questo modo il numero di punti di controllo è pari a 25. Avendo scelto
Figura 2.3 Eliminazione punti di controllo doppi.
curve di ordine 4 (grado 3), il vettore dei nodi presenta 31 elementi, uniformemente
distribuiti nell’intervallo [0, 1]. A tal proposito si veda la sezione 1.2.2. A questo punto
è applicabile lo schema di generazione precedente.
Come in [5] si è scelto un modulo elementare costituito da 2 parti rettilinee che
comprendono una parte curvilinea. La scelta di ottenere il profilo superiore del canale
per semplice traslazione consente di garantire la sovrapponibilità dei moduli. In questo
lavoro tuttavia si è potuto trattare i profili, rispettivamente inferiore e superiore, come
un’unica curva NURBS, andando a controllarne la generazione in modo efficace come
di seguito.
2.1.1
Parametrizzazione del profilo inferiore del canale
Come sopra riportato la generazione del canale ha inizio definendo i punti di controllo di una NURBS periodica del IV ordine da utilizzare per il profilo inferiore.
21
2.1 Definizione della geometria
Coerentemente alla numerazione dei punti di controllo, si sono utilizzate le relazioni
Figura 2.4 Relazioni fra i punti di controllo.
di tabella 2.1 e tabella 2.2 per parametrizzare il profilo inferiore. Il punto 1 coincide
Punto
Ascissa x
Ordinata y
Numero variabili
1
0
0
0
2
x1 + ∆x2
0
0
3
x2 + ∆x3
0
1
4
x5 + ρ4 cos ϑ4
y5 + ρ4 sin ϑ4
2
5
x1 + ∆x5
y1 + ∆y5
2
6
x5 + ρ6 cos ϑ6
y5 + ρ6 sin ϑ6
2
7
x8 − ∆x7
0
1
8
x2 + Lwave
0
1
9
L
0
1
Tabella 2.1 Punti di controllo del modulo elementare lungo il profilo inferiore.
con l’origine degli assi e può essere quindi assunto come riferimento. I punti 3-7 sono
definiti in modo relativo rispetto a 2-8. I punti 4-6 sono definiti in forma polare rispetto al punto 5 assunto come polo, le anomalie sono riferite all’orizzontale. Il punto 2 è
automaticamente determinato, una volta fissate la lunghezza del canale L e la lunghezza proiettata del tratto curvo Lwave , in modo che il tratto curvo sia sempre posizionato
al centro del profilo inferiore. Togliendo un grado di libertà sul posizionamento si evi-
22
Definizione del problema nel caso bidimensionale
ta che, per taluni valori di traslazione del profilo superiore, si ottengano canali dalle
geometrie linearmente dipendenti. ∆x2 ha il ruolo di funzione di trasferimento ed è:
∆x2 =
L − Lwave
2
Il punto 9 è posizionato sulla medesima ordinata del punto 1 e dista da esso la lunghezza del canale.
I campi di validità delle variabili introdotte sono definiti rispetto al diametro idraulico (2.3), unico riferimento geometrico fisso; essendo l’altezza media pari a 0.5, esso
assume valore unitario e se ne può sottintendere la presenza. Il valore dei parametri
Parametro
L
Lwave
∆x3
ρ4
ϑ4
∆x5
∆y5
ρ6
ϑ6
∆7
Descrizione
lunghezza lungo x
del modulo elementare
lunghezza lungo x
del profilo curvilineo
distanza lungo x
del punto 3 rispetto al punto 2
distanza polare
del punto 3 rispetto al punto 2
anomalia polare
del punto 3 rispetto al punto 2
distanza lungo x
del punto 5 rispetto al punto 2
distanza lungo y
del punto 5 rispetto al punto 2
distanza polare
del punto 6 rispetto al punto 5
anomalia polare
del punto 6 rispetto al punto 5
distanza lungo x
del punto 8 rispetto al punto 7
Campo
ainf ; asup
binf L; bsup L
cinf L; csup L
dinf ; dsup
einf ; esup
finf L; fsup L
ginf ; gsup
hinf ; hsup
iinf ; isup
linf L; lsup L
Tabella 2.2 Dominio di definizione delle variabili.
che definiscono i range delle variabili verranno riportati nella presentazione dei risultati. Il processo di ottimizzazione prevede infatti un graduale affinamento di tali valori
in modo da indirizzare l’ottimizzatore nella corretta strada. Finora si contano 10 va-
2.2 Definizione e risoluzione del problema termofluidodinamico
23
riabili, ne verrà introdotta una ulteriore a definire la traslazione rispetto all’origine (o
al punto 1 equivalentemente).
2.1.2
Generazione del profilo superiore
Come già specificato il profilo superiore è ottenuto traslando lungo la direzione delle
ordinate di una quantità Hav = 0.5 i punti di controllo del profilo inferiore e traslando
gli stessi lungo l’asse delle ascisse di una quantità trasl definita come in tabella 2.3.
Vale per i parametri minf e msup quanto detto alla fine della precedente sezione.
Parametro
trasl
Descrizione
entità di traslazione lungo x dei punti
di controllo del profilo superiore
Campo
minf L; msup L
Tabella 2.3 Entità di traslazione del profilo superiore rispetto quello inferiore.
Il campo di variazione dei parametri di generazione geometrica consente la presenza di combinazioni tali da ottenere domini non coerenti con la fisica del problema,
oppure su cui non è possibile costruire mesh valide, oppure ancora che non possono
essere processate dal solver. Per ovviare a questi inconvenienti che portano al blocco
delle operazioni di ottimizzazione, si sono previsti nello script MATLAB alcuni controlli che consentono il salvataggio dell’eventuale errore in un file di testo e l’uscita
dal programma in esecuzione. Nella sezione 3.1 si darà una più amplia descrizione
del problema. Le geometrie eliminate dai controlli di cui sopra non sono realizzabili
o comunque il loro utilizzo non è conveniente e quindi, dopo una eventuale ulteriore
analisi da parte dell’operatore, vengono scartate.
2.2
Definizione e risoluzione del problema termofluidodinamico
Il problema termofluidodinamico preso in considerazione è lo scambio termico in regime di convezione forzata, per cui la riduzione delle perdite di carico nel canale si
traduce in una riduzione della potenza di pompaggio richiesta. Questo consente di risolvere in un primo tempo il campo fluidodinamico, separatamente da quello termico,
dal quale non è influenzato. Una volta noto il campo di moto, si affronta la risoluzione
della parte termica influenzata dalla particolare distribuzione di velocità: l’accoppiamento dei due campi è dunque a senso unico, cosa che consente una più veloce ricerca
della soluzione numerica rispetto al caso di convezione naturale dove il campo di moto
è influenzato dal campo di temperatura.
Il flusso nel canale viene studiato assumendo il regime laminare e stazionario. Si
individuano quindi le caratteristiche di periodicità delle grandezze che governano i
24
Definizione del problema nel caso bidimensionale
fenomeni fisici occorrenti, ossia velocità, pressione e temperatura. La scelta del regime laminare stazionario è legata ad una serie di motivazioni di diverso carattere. Dal
punto di vista numerico, visto il carattere metodologico di questo lavoro, si è deciso
di contenere l’onere computazionale. Inoltre, non essendo richiesto alcun modello di
turbolenza, è possibile in questo modo garantire per il caso bidimensionale che i risultati siano grid-independent, i.e. con bassi errori di discretizzazione e convergenza [6].
Dal punto di vista fisico invece, in svariati casi pratici e in particolare nei dispositivi
di scambio termico a scala ridotta, il regime è necessariamente laminare a seguito del
basso numero di Reynolds [7, 8]. Non va poi trascurato che l’aumento dello scambio
termico si ottiene favorendo il mixing del flusso. Il regime turbolento è caratterizzato da un forte rimescolamento, tuttavia l’incremento del fattore d’attrito richiede alti
costi di pompaggio. Potrebbe quindi risultare vantaggioso sviluppare il fenomeno del
mixing con particolari geometrie in regime non-turbolento [9].
Al fine di svincolarsi dalla necessità di definire le proprietà termofisiche del fluido
in opera e di ottenere soluzioni generali che prescindano dalle dimensioni geometriche, si è scelto di affrontare il problema in forma adimensionale: in questo modo [10, 11] è infatti possibile disporre di una rapida confrontabilità fra due sistemi fra di loro in rapporto di similitudine. La similitudine tra due sistemi è definita
nell’uguaglianza a livello puntuale dei parametri geometrici, cinematici e dinamici.
2.2.1
Grandezze adimensionali e grandezze di riferimento
Il processo di adimensionalizzazione consiste nel rapporto fra due grandezze dimensionalmente omogenee ed esprime l’importanza dell’una rispetto all’altra assunta come riferimento, fisso o variabile. Per procedere alla adimensionalizzazione delle equazioni necessarie alla risoluzione del problema di convezione forzata si fissano tre grandezze di riferimento, di conseguenza si impongono tre scale per le variabili in gioco.
In questo lavoro si è scelto di fissare una scala per le lunghezze, una per le velocità
e una per la temperatura. Sono possibili ovviamente differenti modalità di scalaggio, operazione del tutto convenzionale, tuttavia in questo modo si ha una valutazione
coerente di tutte le grandezze ricavate, si veda paragrafo 2.2.4.
Lunghezze Tutte le lunghezze sono rapportate ad una dimensione di riferimento:
nel nostro caso il diametro idraulico medio del canale Dh . Questa definizione consente una facile confrontabilità tra risultati ottenuti e quelli presenti in
letteratura, nonché tra casi 2D e 3D.
Velocità Le due componenti di velocità u e v sono adimensionalizzate rispetto la velocità media nel canale Uav nella direzione x, ossia la sola componente di velocità
che contribuisce alla portata.
Temperatura Detta Tw la temperatura della parete del canale e Tb, in la temperatura
2.2 Definizione e risoluzione del problema termofluidodinamico
25
media entalpica (detta di bulk) nella sezione di ingresso si ha che
θ=
T (x, y) − Tw
Tb, in − Tw
(2.1)
Pressione In osservanza del teorema di Buckingham, il termine della pressione è
2 , dove ρ è la densità del
implicitamente adimensionalizzato dal gruppo ρ · Uav
fluido.
Il diametro idraulico è definito come
Dh =
4Ac
C
(2.2)
dove Ac è la sezione di efflusso e C è il perimetro bagnato del canale. Nel caso di un
canale 2D con altezza media Hav vale che
Dh = lim
z→∞
4 · z · Hav
= 2Hav
2(z + Hav )
(2.3)
La velocità media viene definita attraverso la portata volumica V̇ , costante per
l’incompressibilità del fluido, come
Uav =
V̇
Hav
(2.4)
La temperatura media entalpica, detta di bulk, o mixing cup temperature [12], è
definita da un bilancio entalpico sulla generica sezione del canale normale all’asse x
come
Z
1
Tb,x =
u T dS
(2.5)
Ac uav Ac
dove (Ac uav ) è la portata volumetrica, uav è la velocità media nella sezione Ac , mentre
u è la componente di velocità normale alla sezione generica Ac . Il pedice x indica il
valore locale in corrispondenza della generica sezione.
Il flusso e lo scambio termico bidimensionali vengono descritti da tre equazioni di
conservazione che in forma adimensionale presentano i seguenti raggruppamenti.
Reynolds Questo gruppo rappresenta il rapporto fra le forze di inerzia e le forze
viscose in un canale di specificata geometria:
Re =
ρ Uav Dh
µ
dove µ è la viscosità dinamica del fluido. Si fa notare come nel caso bidimensionale per flusso incomprimibile il prodotto Uav Hav è la portata volumica per unità
di profondità ed è costante; segue che anche il prodotto Uav Dh = Uav 2Hav = cost
26
Definizione del problema nel caso bidimensionale
Prandtl Questo numero esprime il rapporto fra la diffusione di quantità di moto e la
diffusione termica:
µ cp
ν
Pr =
=
k
α
Péclet Il prodotto Re Pr è detto numero di Péclet, il quale rappresenta il rapporto fra
l’energia termica fornita al fluido per convezione e l’energia termica trasmessa
assialmente per conduzione lungo il flusso.
La tabella 2.4 riassume quanto esposto, le variabili adimensionali sono indicate con
l’asterisco in apice. Con Tw si intende la temperatura in corrispondenza delle pareti
Grandezza
Espressione
x∗
x/Dh
y∗
y/Dh
u∗
u/Uav
v∗
v/Uav
p∗
2)
p/(ρ Uav
θ
(T (x, y) − Tw )/(Tb, in − Tw )
Re
(ρ Uav Dh )/µ
Pr
(µ cp )/k
Pe
Re Pr
Tabella 2.4 Elenco delle grandezze adimensionali per il caso 2D.
e in genere il pedice w indica wall. Con i pedici in e out si intende fare riferimento
rispettivamente alla sezioni di ingresso e di uscita del canale, per cui xin = 0 e xout = L,
se L è il periodo del canale.
2.2.2
Equazioni adimensionali
Si riportano di seguito le equazioni di conservazione in forma adimensionale. D’ora
in avanti, se non specificatamente espresso, si sottintende l’asterisco per indicare le
grandezze adimensionali. Il regime del flusso all’interno del canale è assunto:
- laminare;
- stazionario;
- completamente sviluppato dinamicamente e termicamente;
2.2 Definizione e risoluzione del problema termofluidodinamico
27
- incompressibile.
Conservazione della massa
∂u ∂v
+
=0
∂x ∂y
(2.6)
Conservazione della quantità di moto
u
∂u
∂u
1
+v
=
∂x
∂y
Re
∂v
∂v
1
+v
=
∂x
∂y
Re
u
∂ 2u ∂ 2u
∂p
+
−
∂ x2 ∂ y2
∂x
(2.7a)
∂ 2v ∂ 2v
∂p
+
−
∂ x2 ∂ y2
∂y
(2.7b)
Conservazione dell’energia
u
∂θ
1
∂θ
+v
=
∂x
∂y
Pe
∂ 2θ ∂ 2θ
+ 2
∂ x2
∂y
(2.8)
La simulazione numerica viene condotta per un flusso con le seguenti caratteristiche:
Re
200
Pr
0.7
Tabella 2.5 Caratteristiche del flusso
Il numero di Prandtl scelto è rappresentativo dell’aria in condizioni standard.
Il canale creato viene studiato in condizioni di flusso completamente sviluppato.
Dal punto di vista fluidodinamico questo comporta che i profili di velocità varino
periodicamente, ritrovando la medesima configurazione ad ogni modulo [13]. Questo
non accade invece nella regione di ingresso o nel caso di flusso non stazionario. La
regione di ingresso in un canale piano, con moto a falde parallele, è caratterizzata dallo
sviluppo dello strato limite nel senso del moto, finché questo non raggiunge un’altezza
pari alla semi-altezza del canale [12]. Il generico canale trattato in questo lavoro non
è tuttavia assialsimmetrico, quindi lo strato limite non si sviluppa in modo simmetrico
nel condotto.
La condizione di completo sviluppo è ottenuta dal punto di vista fluidodinamico imponendo [4] condizioni di periodicità sul campo di velocità e pressione come
riportato nel paragrafo 2.2.3.
Dal punto di vista termico lungo un condotto caratterizzato da un fluido in moto
che scambia calore con l’esterno, detta s la coordinata curvilinea e n la coordinata in
28
Definizione del problema nel caso bidimensionale
direzione ortogonale a s si ha che:
dTb
,0
ds
cosı̀ come
∂ T (s, n)
,0
∂s
mentre considerando il gruppo adimensionale
θp =
T (x, y) − Tw
Tb,x − Tw
(2.9)
la condizione di completo sviluppo per una canale a sezione costante [4] è espressa
dalla
∂ θp
=0
(2.10)
∂s
Si può inoltre dimostrare [12] che in un deflusso interno nelle condizioni citate e con
un fluido di processo a proprietà costanti, il coefficiente convettivo, cosı̀ come il numero di Nusselt, sono costanti e indipendente da s. Le condizioni per l’esistenza di una
zona a completo sviluppo termico in canali circolari in regime di convezione forzata
sono riportate in [14]. Per canali a sezione variabile con periodicità L la precedente
relazione (2.10) vale solo in termini integrali (2.11), ovvero il profilo di temperatura
adimensionale θp si ripete uguale nelle diverse stazioni distanti L fra loro lungo la
direzione del flusso [4].
θp (x, y) = θp (x + L, y)
(2.11)
È interessante notare che un canale a sezione costante può essere visto come un canale
periodico di periodo infinitesimo.
L’aver fissato le proprietà del fluido e le caratteristiche del problema in forma adimensionale con Pr e Re a questo livello di analisi non è rilevante. Infatti il numero di
Reynolds determina la lunghezza del tratto di tubazione detta di sviluppo idrodinamico, mentre il numero di Prandtl la lunghezza di tubazione in cui si verifica lo sviluppo
termico: al variare di questi due parametri si osservano variazioni della soluzione nella
sola regione di ingresso (o eventualmente in condizioni di non stazionarietà).
Si è temporaneamente introdotto un sistema di coordinate curvilinee, il discorso
non cambia però nel piano (x, y) di coordinate cartesiane, scelta peraltro obbligata per
l’utilizzo di FEMLAB.
2.2.3
Impostazione delle equazioni per il campo fluido
Si consideri un flusso all’interno di un canale bidimensionale sviluppato idrodinamicamente con l’asse delle ascisse positive nel senso di avanzamento del fluido come in
figura 2.5. La periodicità si osserva nella direzione x, ossia il canale viene idealmen-
2.2 Definizione e risoluzione del problema termofluidodinamico
29
Figura 2.5 Soluzione del campo fluido: rappresentazione del campo con vettori di velocità e
linee di corrente.
te montato unendo lungo tale direzione un certo numero di moduli elementari, come
quello in figura 2.5.
Come in [4] ed altri autori in articoli successivi, detta L la lunghezza del canale e
Y = Y (x) la generica ordinata del profilo vale che
Y (x) = Y (x + n · L)
(2.12)
con n ∈ Z.
Condizioni al contorno per il campo fluidodinamico.
Per quel che riguarda le condizioni al contorno le componenti di velocità devono
assumere lo stesso valore in ingresso e in uscita al canale, mentre saranno nulle in
corrispondenza del mantello laterale.
La pressione merita invece un discorso a parte. L’effetto del campo di pressione
è quello di consentire un efflusso che soddisfi la (2.6) contro la presenza dell’attrito
dovuto alla viscosità del fluido di processo: si assiste quindi ad una perdita di carico
lungo la direzione positiva dell’asse delle ascisse. Questo non significa però che non
si possa dedurre una condizione di periodicità per la pressione. Il termine di pressione
può infatti essere pensato come somma di 2 contributi [4, 5]:
1. un termine periodico al pari delle componenti di velocità;
2. un termine a-periodico che tiene conto della continua caduta di pressione lungo
il canale.
L’andamento della pressione media in ogni sezione in funzione di x può essere espresso da una curva come in figura 2.6. Se L è la lunghezza del modulo del canale, si
osserva che
p(x, y) − p(x + L, y) = β · L
(2.13)
dove β > 0 è il gradiente pressorio ed esprime la caduta di pressione per unità di
lunghezza necessaria a garantire la portata desiderata, mentre P(x, y) rappresenta la
30
px
Definizione del problema nel caso bidimensionale
andamento globale
caduta di pressione dovuta al gradiente pressorio
componente periodica
x
Figura 2.6 Andamento caratteristico della pressione media in ciascuna sezione del canale e
sue componenti.
componente periodica della pressione, responsabile della distribuzione del campo di
moto nel modulo. La pressione all’interno del campo fluido può quindi essere espressa
come
p(x, y) = −β · x + P(x, y)
(2.14)
Verranno ora riportate le condizioni al contorno e le equazioni per la risoluzione
nel dominio di interesse, ossia quello per cui è 0 ≤ x ≤ L .
Pareti Sulle pareti sono imposte condizioni di no-slip, ovvero componenti di velocità
nulle per effetto della viscosità.
u=0 v=0
(2.15)
Inlet e Outlet L’ingresso e l’uscita del canale necessitano di condizioni di periodicità
per le componenti di velocità e per la pressione, in questo modo è infatti possibile studiare un solo modulo dell’intero canale e non ci si deve preoccupare
della regione di ingresso in cui si assiste alla condizione di hydrodynamically
developing flow.
u(xin , y) = u(xout , y) v(xin , y) = v(xout , y) P(xin , y) = P(xout , y)
(2.16)
Infine è necessario fissare un riferimento per la pressione, imponendone un valore
arbitrario in un punto: i.e. P = 0 nel punto in basso a destra del canale. Senza quest’ultima infatti, le condizioni al contorno fornite definirebbero la pressione a meno di
31
2.2 Definizione e risoluzione del problema termofluidodinamico
Figura 2.7 Imposizione di un punto a pressione nulla come valore di riferimento.
una costante additiva: la soluzione trovata non sarebbe univoca. Per quanto riguarda
le equazioni da risolvere per ottenere il campo fluidodinamico la (2.6) rimane immutata, mentre la (2.7a) e la (2.7b) presentano il termine di pressione periodica assieme
al gradiente gradiente pressorio β nella forma di forzante
∂P
∂u
∂u
1 ∂ 2u ∂ 2u
+
+β −
u +v
=
(2.17a)
∂x
∂y
Re ∂ x2 ∂ y2
∂x
∂v
1
∂v
=
u +v
∂x
∂y
Re
2.2.4
∂ 2v ∂ 2v
∂P
+
−
∂ x 2 ∂ y2
∂y
(2.17b)
Fattore di attrito adimensionale
Come indice della potenza spesa dalla pompa per movimentare il fluido all’interno del
canale, si è scelto un parametro adimensionale, quale il fattore d’attrito f , in modo da
poter subito confrontare diverse geometrie. Autori diversi presentano differenti definizioni di tale valore (Fanning friction factor, Darcy friction factor [12]), tuttavia non
sempre c’è coerenza fra le varie grandezze correlate a tale fattore. In questo lavoro si
è data una definizione autonoma, comunque conforme alla letteratura corrente, in modo da avere coerenza tra il fattore d’attrito ottenuto dalla (2.22), ovvero dall’integrale
degli sforzi a parete, e lo stesso ricavato dalla (2.23).
In questa sezione le forme adimensionali sono indicate con l’asterisco. Utilizzan2 per adimensionalizzare tensioni e pressioni, si definisce il fattore
do il termine ρUav
d’attrito locale
τnx
fx =
(2.18)
2
ρUav
dove τnx rappresenta lo sforzo viscoso a parete in direzione x, wall shear stress, il
quale può essere espresso in funzione della pressione e delle componenti di velocità
32
Definizione del problema nel caso bidimensionale
[11]
∂u ∂v
∂u
nx + µ
ny
τnx = −p + µ 2
+
∂x
∂y ∂x
(2.19)
In forma adimensionale, seguendo quanto già fatto per l’equazione della conservazione della quantità di moto
2 ∂ u∗
1 ∂ u∗ ∂ v∗
∗
+
fx = τnx
= −p∗ +
n
+
ny
(2.20)
x
Re ∂ x∗
Re ∂ y∗ ∂ x∗
Noto il fattore d’attrito locale fx , quello globale si ricava con la media sul canale
f=
1
L
Z
wall
fx ds
Sostituendo l’espressione di fx con la (2.20) e ricordando la (2.14)
Z
1
2 ∂ u∗
1 ∂ u∗ ∂ v∗
f=
−P∗ + β ∗ x∗ +
+
n
+
ny ds
x
L wall
Re ∂ x∗
Re ∂ y∗ ∂ x∗
(2.21)
(2.22)
Il coefficiente f altro non è se non la forma adimensionale della forza d’attrito media
sul canale. Per garantire la condizione di completo sviluppo, dalla conservazione
della quantità di moto fra ingresso e uscita, come semplice bilancio di forze mediate
nel canale
Z
Fshear stress =
τnx ds = ∆p Hav
(2.23)
wall
ed elaborando la (2.23) assieme alle (2.18) e (2.21) si ricava infine
f =β
β∗
Hav
=
2
ρU
2
(2.24)
La (2.24) è la forma utilizzata per determinare il parametro d’attrito medio del canale,
il quale coincide con la metà del gradiente pressorio in forma adimensionale. Procedendo in questo senso, si può verificare che le (2.22) e (2.24) portano allo stesso
risultato.
2.2.5
Tecniche di soluzione per il campo fluidodinamico
Le equazioni (2.6), (2.17a), (2.17b) costituiscono un sistema di equazioni alle derivate
parziali [4] sulle 3 incognite:
u(x, y)
v(x, y)
P(x, y)
2.2 Definizione e risoluzione del problema termofluidodinamico
33
Si osservi che però il sistema presenta, oltre a Re, che è un parametro, un’ulteriore
incognita, ovvero β , legato come noto alla portata che scorre all’interno del canale.
Per affrontare la soluzione si è scelto di fissare un valore per la portata V̇ = 0.5,
che unitamente alla scelta per il valore di Hav , altezza del canale, consenta di utilizzare
scale di lunghezza e velocità normalizzate sul valore unitario. In questo modo si ha:
Re ≡
1
ν
Per ovviare alla singolarità del sistema di cui sopra, il processo di soluzione è di tipo
iterativo: si fissa per β un valore di tentativo, si procede alla soluzione del sistema e
lo si corregge finché la velocità media ottenuta non è unitaria, come nello schema di
figura 2.8.
Il gradiente pressorio è dunque il responsabile della velocità media del flusso nel
canale, ad una sua graduale diminuzione si osserva una riduzione di portata. Si sono
tentati 3 differenti approcci per risolvere il campo fluido con il valore di portata desiderato, cui corrisponde un valore unitario per la velocità media, Uav = 1. Di seguito si
riportano le modalità di correzione di cui in figura 2.8.
2.2.6
Primo approccio
Una prima soluzione consiste nel generare una correzione per il gradiente pressorio di
primo tentativo a partire dall’errore sulla portata media nel canale [5, 15]. Il controllo
che viene applicato su β può essere una combinazione di
proporzionale P, integrativo I, derivativo D.
Nel lavoro [5] viene applicato un controllo P, in cui la correzione su β è proporzionale all’errore percentuale sulla portata. Per agevolare la convergenza il coefficiente
di proporzionalità è di tipo adattativo a seconda dell’evoluzione della funzione errore.
Nel lavoro [15] il controllo è di tipo PI (proporzionale-integrativo), il che consente
un certo smorzamento per evitare eccessivi overshoots della funzione di correzione.
Entrambe le strade precedenti sono tuttavia risultate poco efficienti e la dipendenza
dal valore del gradiente pressorio di primo tentativo si è dimostrata notevole. Anche
con un controllo PID [16], con termine integrativo che influenza un numero limitato
di errori, si è reso necessario avere un controllo in itinere dei vari coefficienti e la
ricerca di valori costanti e di compromesso mediante modeFRONTIER su diverse
geometrie e con diversi β iniziali non ha fornito risultati soddisfacenti. Si ritiene che
questo approccio alla soluzione sia efficiente nel caso di problemi specifici, in cui la
geometria non abbia i gradi di libertà necessari in un problema di ottimizzazione. In
questo lavoro è richiesta una forte flessibilità e una metodologia troppo conservativa
(i.e. correzioni contenute ad ogni iterazione) rischia di portare a tempi di calcolo non
sostenibili.
34
Definizione del problema nel caso bidimensionale
Imposizione del numero di Reynolds Re
e del gradiente pressorio β di tentativo
Risoluzione del problema fluidodinamico
Calcolo della portata
e della velocità
media
Nuovo valore per
il gradiente
pressorio
Valutazione dell'errore
rispetto la velocità
di riferimento
SI
L'errore supera la
tolleranza richiesta
NO
Il campo di velocità e
il valore di β ottenuti
sono quelli ricercati
Figura 2.8 Diagramma di flusso per la risoluzione del campo fluidodinamico.
2.2 Definizione e risoluzione del problema termofluidodinamico
2.2.7
35
Secondo approccio
2 , si
Guardando alla forma dimensionale per maggior semplicità, vale che β ∝ fUav
veda equazione (2.24). La caduta di pressione all’interno del modulo è infatti assunta, attraverso il fattore d’attrito, proporzionale all’energia cinetica del flusso. Si può
dunque procedere in forma iterativa nel seguente modo:
1. si sceglie un valore di tentativo per il gradiente pressorio βn
2. si risolve il campo fluidodinamico e si calcola il valore della portata volumica e
della velocità media Uav,n
3. la coppia (Uav,n , βn ) ottenuta individua il coefficiente d’attrito fn che viene
calcolato con la (2.24), si veda figura 2.9
4. si determina un nuovo valore del gradiente pressorio βn+1 ponendo Uav = 1,
ovvero il valore obiettivo per la velocità media, nella (2.24)
5. si risolve nuovamente il campo di moto, estraendo la portata volumica e la
velocità media Uav,n+1
6. se l’errore tra il nuovo valore di velocità media e il valore unitario obiettivo è
inferiore ad una fissata tolleranza, si va al punto 7., altrimenti si riparte dal punto
1. con βn+1 appena trovato come nuovo valore di tentativo
7. la coppia (Uav,n+1 , βn+1 ) ≡ Uobj , βobj e rispetta il valore di portata richiesto ed
è coerente con il coefficiente d’attrito calcolato fn+1 ≡ fobj
Il pedice n rappresenta un contatore di iterazioni.
36
Definizione del problema nel caso bidimensionale
Tentativo numero 1
1
0.9
gradiente pressorio
0.8
0.7
1
1
(Uav,β )
0.6
0.5
0.4
0.3
(1,β2)
0.2
0.1
0
0
0.5
1
U
1.5
2
av
1 ,β1
Figura 2.9 Gradiente pressorio di primo tentativo e prima correzione: la coppia Uav
1 > 1 e il nuovo valore di β 2 viene calcolato
individua un coefficiente f 1 < fobj in quanto Uav
con Uav = 1.
Tentativo numero 2
gradiente pressorio
1.5
1
(1,β3)
0.5
2
2
(Uav,β )
0
0
0.5
1
Uav
1.5
2
2 < 1, quindi questa coppia individua f 2 >
Figura 2.10 Seconda correzione: β 2 < βobj e Uav
fobj .
37
2.2 Definizione e risoluzione del problema termofluidodinamico
Tentativo numero 3
2
1.8
gradiente pressorio
1.6
1.4
1.2
1
0.8
(U3av,β3)
0.6
(1,β4)
0.4
0.2
0
0
0.5
1
U
1.5
2
av
Figura 2.11 Terza correzione: come 2.9 con però un errore minore sulla portata/velocità
media.
Tale metodo di correzione è stato tuttavia abbandonato perché la velocità di convergenza è piuttosto lenta a causa della dipendenza quadratica del gradiente pressorio
dalla velocità media nel canale. Dalle figure riportate si può notare l’oscillazione della soluzione attorno al presunto valore obiettivo: questo problema è assimilabile ad
un sistema del secondo ordine con un forte overshooting, causato dalla quasi totale
assenza dello smorzamento. Rispetto al caso precedente dunque non si è osservato
un miglioramento dell’efficienza, tuttavia il meccanismo presenta una maggiore robustezza, avendo il valore di primo tentativo di β una minor influenza sulla possibilità di
convergenza.
2.2.8
Terzo approccio
In regime laminare [11, 12, 17], per un flusso in un canale cilindrico o in un canale
piano parallelo, il fattore d’attrito f , sia dal punto di vista analitico sia dal punto di
vista sperimentale, è inversamente proporzionale al numero di Reynolds: si utilizza
questo aspetto fisico del problema per ottenere una efficiente convergenza numerica.
Sempre in forma dimensionale
K
f=
(2.25)
Re
dalla definizione di Re si evince che
f=
K0
Uav
(2.26)
38
Definizione del problema nel caso bidimensionale
Riproponendo la dipendenza del gradiente pressorio dalla velocità media (2.24)
2
β ∝ fUav
(2.27)
1 2
U = K 00Uav
Uav av
(2.28)
Inserendo nella precedente la (2.26)
β = K 00
Ovviamente lo stesso discorso si ripropone con la portata volumetrica in quanto V̇ =
Uav Hav ∝ Uav . Analogamente al caso di cui sopra, ora è K 00 incognito e può essere
determinato nota la generica coppia di tentativo (β ,Uav ). La relazione (2.25) non vale
in modo esatto per il flusso in un canale con sezione di efflusso variabile, anzi si può
osservare, anche in virtù della modalità precendente, che
m
β = KUav
(2.29)
con 1 6 m 6 2. La condizione di regime laminare, lontano dalla condizione di transizione, assicura però che il coefficiente m sia più vicino al valore unitario. A questo
aspetto è dovuta l’eventuale non perfetta linearità della relazione (2.28), rappresentata
in figura 2.12 dalla curva a trattini. Questo ed errori di round off fanno sı̀ che per
convergere alla portata obiettivo sia necessario un numero di iterazioni superiore a 1.
Tuttavia si fa notare come nel canale piano parallelo la soluzione corretta è ottenuta
già alla prima iterazione.
Si può procedere quindi nel modo seguente:
1. si sceglie un valore del gradiente pressorio βn di tentativo
2. si risolve il campo di moto calcola il valore della portata e della velocità media
Uav,n
3. la coppia (Uav,n , βn ) individua la costante di proporzionalità K 00 si veda figura
2.12
4. si determina un nuovo valore del gradiente pressorio βn+1 =
βn
con cui
Uav,n
ricalcolare il campo fluidodinamico
5. si calcola nuovamente il valore della velocità media Uav,n+1
6. se l’errore fra la velocità media appena calcolata e la velocità obiettivo è inferiore ad un fissato valore di tolleranza, si prosegue al punto 7., altrimenti si ritorna
al punto 1. con βn+1 come nuovo valore di tentativo
7. la coppia (Uav,n+1 , βn+1 ) ≡ Uobj , βobj e rispetta il valore di portata richiesto
2.2 Definizione e risoluzione del problema termofluidodinamico
39
In realtà nello script la correzione è ottenuta con il rapporto fra la portata calcolata
nella sezione di ingresso e la portata obiettivo: a convergenza avvenuta tale rapporto è unitario, cosı̀ come il rapporto fra le velocità medie, qualora calcolate lungo la
medesima sezione.
Tentativo numero 1
3
1
gradiente pressorio
1
(Uav,β )
2.5
2
2
1.5
(1,β )
1
0.5
0
0
0.5
1
U
1.5
2
av
Figura 2.12 Gradiente pressorio di primo tentativo e prima correzione.
Tentativo numero 2
3
gradiente pressorio
2.5
2
(U2av,β2)
1.5
3
(1,β )
1
0.5
0
0
0.5
1
Uav
1.5
2
Figura 2.13 Gradiente pressorio di secondo tentativo e seconda correzione.
40
Definizione del problema nel caso bidimensionale
Tentativo numero 3
3
gradiente pressorio
2.5
2
1.5
(U3 ,β3)
av
(1,β4)
1
0.5
0
0
0.5
1
U
1.5
2
av
Figura 2.14 Gradiente pressorio di terzo tentativo e terza correzione: ormai si è molto vicini
alla soluzione obiettivo e l’inviluppo delle rette dà l’effettiva relazione fra velocità/portata e
gradiente pressorio.
L’evoluzione delle figure 2.12, 2.13, 2.14 mette in evidenza come la curva che
esprime la (2.28) in forma non lineare sia l’inviluppo delle (2.28) a K 00 = const; la
non linearità della (2.28) è giustificata in presenza di un moto laminare su canali con
profili più o meno distanti dal canale piano parallelo, per il quale la (2.25) è valida con
K = const.
La convergenza è di tipo monotono, ovvero se il gradiente pressorio di primo tentativo βn è maggiore di βobj le successive correzioni consentiranno di giungere a tale
valore rimanendo alla sua destra rispetto al sistema di coordinate in figura. L’opposto se βn < βobj . La motivazione di questa peculiarità risiede nella monotonia della
derivata prima della (2.28).
2.2.9
Impostazione delle equazioni per il campo termico
Tra le numerose condizioni al contorno che si possono imporre per il campo termico in
regime di convezione forzata [12], si è optato per temperatura uniforme alle pareti del
condotto. Questa, oltre ad essere rappresentativa di alcune situazioni reali, i.e. elevata
conduttività termica delle pareti, è più complessa da trattare numericamente rispetto
al caso di flusso termico imposto alle pareti. Questa scelta comporta infatti per il
flusso termicamente sviluppato una formulazione della periodicità diversa rispetto a
quella del problema fluidodinamico [4], mentre per flusso termico imposto alle pareti
si procede similmente.
2.2 Definizione e risoluzione del problema termofluidodinamico
41
Come già evidenziato, la condizione di periodicità (2.11) per il profilo di temperatura adimensionale riguarda θp e non la forma adimensionale θ scelta per la temperatura in questo lavoro e attraverso la quale è stata riscritta l’equazione di conservazione
dell’energia (2.8).
Seguendo il processo di adimensionalizzazione di [4], compare un termine, sotto
forma di forzante di volume, il quale tiene conto del fatto che
∂ Tb,x
,0
∂x
(2.30)
Il vantaggio di risolvere l’equazione dell’energia nella variabile primitiva θp , automaticamente periodica, comporta quindi la necessità di determinare lungo il dominio la
distribuzione della variazione della temperatura di bulk in forma dimensionale (2.30).
Questo significa introdurre una ulteriore equazione [4], piuttosto complessa da implementare sul tipo di griglia non strutturata che FEMLAB utilizza. Attraverso invece
una adimensionalizzazione per la temperatura θ a riferimento fisso, consistente nella
temperatura di bulk in corrispondenza della sezione di ingresso al canale, si aggira
questo problema e l’equazione dell’energia assume la forma di (2.8). In questo modo
però la variabile indipendente non è più periodica fra inlet e outlet, poiché la periodicità si ha sul profilo di temperatura del tipo (2.9). Si osserva tuttavia che, nelle
condizioni di completo sviluppo termico, valgono le seguenti relazioni
ovvero
θp (xin , y) = θp (xout , y)
(2.31)
T (L, y) − Tw
T (0, y) − Tw
=
Tb,in − Tw
Tb,out − Tw
(2.32)
Per riportarsi alla forma adimensionale a punto fisso, si moltiplica e si divide il secondo membro della 2.32 per il termine (Tb,in − Tw )
T (xout , y) − Tw Tb,in − Tw
T (xin , y) − Tw
=
Tb,in − Tw
Tb,out − Tw Tb,in − Tw
Detto
σ=
Tb,in − Tw
Tb,out − Tw
(2.33)
(2.34)
dalle (2.33) e (2.34) si ottiene che
θ (xin , y) = θ (xout , y)σ
(2.35)
Si è quindi sostituita una funzione di x [4], con una costante che rispetti la successiva
condizione (2.42): ovvero un salto di temperatura coerente con il flusso termico al
mantello ed il salto entalpico tra inlet e outlet del canale. Riguardo la determinazione
del parametro σ si discuterà più avanti.
42
Definizione del problema nel caso bidimensionale
Condizioni al contorno per il campo termico
Viene seguita la strada di [5], riproposta in forma adimensionale e considerando, senza
togliere generalità al metodo, che il fluido subisca una fase di raffreddamento.
Pareti Sulle pareti solide la condizione è di temperatura imposta e costante, per cui
dalla definizione di θ si ottiene
θ =0
(2.36)
Inlet e Outlet Per garantire la condizione di completo sviluppo termico all’interno
del canale si impone, come già espresso, la (2.35) unita alla (2.34).
Tuttavia come già osservato in [5], con un tale set di condizioni al contorno FEMLAB
fornisce unicamente la soluzione banale, ossia
θ (x, y) = θw = 0
Infatti una soluzione di questo tipo è coerente con le (2.34) (2.35) e (2.36), al pari
della soluzione non-banale ricercata e, almeno qualitativamente, si può dire che il
solutore è portato a questa soluzione poiché gli estremi della sezione di ingresso e di
uscita hanno la medesima condizione al contorno del mantello. Per risolvere questo
problema si impone in ingresso una ulteriore condizione di tipo integrale, tenendo
conto della (2.5), insieme alla incomprimibilità del fluido e alla definizione di θ ,
R
θ u ds
= θb,in = 1
in u ds
Rin
(2.37)
Cosı̀ facendo, si costringe il solutore a trovare quella particolare distribuzione di
temperatura lungo la sezione di ingresso che rispetti le condizioni agli estremi e al
contempo non sia ovunque nulla.
In figura 2.15 è proposto un esempio di canale bidimensionale con le principali
condizioni al contorno imposte e per il campo di moto e per il campo termico.
2.2.10
Tecniche di soluzione per il campo termico
Analogamente al campo fluidodinamico, la (2.8) presenta come unica variabile primitiva la temperatura adimensionale θ , tuttavia nelle condizioni di periodicità compare
il parametro incognito σ , la cui determinazione viene ora affrontata.
L’andamento di Tb,out in funzione di σ (2.35), figura 2.16, è di tipo iperbolico ed
è interessante notare come il parametro σ assuma valori maggiori dell’unità sia nel
caso di canale con flusso in fase di raffreddamento sia di canale con flusso in fase di
riscaldamento. Infatti, dall’applicazione del secondo principio della termodinamica a
questa tipologia di problema, la temperatura di bulk in uscita è sempre intermedia tra la
la temperatura di bulk in ingresso e la temperatura di parete. Se (|σ | = 1) dalla (2.34)
si ottiene che non è avvenuto alcuno scambio termico fra ingresso e uscita del canale.
43
2.2 Definizione e risoluzione del problema termofluidodinamico
Figura 2.15 Principali condizioni al contorno imposte.
La condizione (|σ | < 1) è sufficiente, ma non necessaria, per garantire la presenza di
pozzi o sorgenti di calore interne al canale (i.e.: reazioni chimiche). La presenza del
modulo è necessaria poiché valori negativi di σ non hanno significato fisico, anche se
numericamente nel processo di soluzione può presentarsi la condizione di σ < 0.
20
Tb,out=f(σ)
Tb,out=f(|σ|)
15
10
Tb,out
5
0
−5
−10
−15
−20
−5
−4
−3
−2
−1
0
σ
1
2
3
4
5
Figura 2.16 Andamento della temperatura dimensionale media all’uscita in funzione del
parametro σ : la zona con |σ | > 1 rappresenta il caso di canale in fase di raffreddamento.
Il processo di risoluzione del campo termico è iterativo. Una volta risolto il campo
fluidodinamico, si utilizza il campo di velocità trovato per ottenere la distribuzione
44
Definizione del problema nel caso bidimensionale
della temperatura adimensionale nel caso di raffreddamento: Tb,in > Tb,out e σ > 1. I
passi del processo risolutivo prevedono di:
1. fissare un valore per σ di tentativo, con il solo vincolo che sia σ > 0
2. risolvere il campo termico con le condizioni al contorno (si veda sezione 2.2.9)
3. imporre l’uguaglianza in forma adimensionale fra il salto entalpico fra ingresso
e uscita del canale e il calore ceduto attraverso il mantello
4. dalla precedente espressione ricavare un nuovo valore per il parametro σnew
5. se l’errore relativo fra σ e σnew supera una certa tolleranza tornare a 1. con
σ = σnew , altrimenti andare a 6.
6. il valore σ determinato è tale per cui il flusso di calore al mantello bilancia il
salto entalpico fra ingresso e uscita del canale
Lo schema è proposto in figura 2.17. L’assunzione di un valore non corretto per σ
causa la nascita di un termine di generazione in corrispondenza della sezione di uscita
del canale. Questo serve a compensare il bilancio di primo principio, ma non ha alcun
significato fisico per il problema trattato. L’uguaglianza in forma adimensionale fra
salto entalpico e flusso termico al mantello si può esprimere, partendo al solito dalla
forma dimensionale, come
ṁcp (Tb,in − Tb,out ) =
Z
w
−k
∂T
ds
∂n
(2.38)
dove ṁ è la portata massica, cp è il calore specifico a pressione costante e k è la
conducibilità termica molecolare del fluido. Dalla definizione di Pr e Re e per la
condizione di temperatura uniforme al mantello
Tb,in − Tb,out =
2
Re Pr
Z
−
w
∂θ
(Tb,in − Tw ) ds
∂n
(2.39)
Sommando e sottraendo la temperatura di parete al primo membro della (2.39) si
ottiene
Z
Tb,in − Tw + Tw − Tb,out
∂θ
1
−
ds
(2.40)
=
Tb,in − Tw
Re Pr w ∂ n
Ricordando la definizione data a σ (2.34)
1−
2
1
=
σ
Re Pr
Z
w
−
∂θ
ds
∂n
(2.41)
Per evitare la determinazione durante il processo iterativo di valori negativi di σ , la
(2.41) può essere riscritta
Z
2
∂ θ 1 = 1+
ds
(2.42)
σ Re Pr w ∂ n 2.2 Definizione e risoluzione del problema termofluidodinamico
Imposizione del valore di θb,in e
del rapporto fra le differenze di
temperatura in ingresso e uscita σ
di tentativo
Risoluzione del problema termico
in convezione forzata
Nuovo valore
per σ
Imposizione dell'uguaglianza in forma adimensionale
tra il salto entalpico fra ingresso e uscita del canale
e il calore ceduto attraverso il mantello
Dall'espressione precedente valutazione
di un nuovo valore per σ
L'errore fra il nuovo
e il vecchio σ supera
la tolleranza richiesta
SI
NO
Il campo termico e il valore
di σ sono fra loro coerenti
Figura 2.17 Diagramma di flusso per la risoluzione del termico.
45
46
Definizione del problema nel caso bidimensionale
Se il coefficiente di attrito valuta la prestanza del canale in esame dal punto di vista
fluidodinamico, le prestazioni di scambio termico sono caratterizzate dal numero di
Nusselt Nu. Come per il coefficiente f , si parla ci numero di Nusselt locale Nux e
mediato lungo una parete Nu.
Nux =
Nu =
hx Dh
k
1
L
(2.43a)
Z L
0
Nux ds
(2.43b)
Poiché un canale è delimitato da 2 pareti la (2.43b), diventa
Nu =
1
2L
Z
wall
Nux ds
(2.44)
in questo modo si è mediato lo scambio termico adimensionale su ciascuna parete
con la lunghezza del modulo. Si sono introdotti il coefficiente di scambio termico
convettivo locale hx e globale h. Esiste infine un legame fra σ e Nu, che ne consente
un facile calcolo, una volta noto il rapporto fra la temperatura di bulk in ingresso e in
uscita. Riscrivendo la (2.38) in forma differenziale
ṁcp dTb,x = −hx (Tb,x − Tw ) ds
(2.45)
Questa è una equazione differenziale del primo ordine, integrabile separando le variabili
dTb,x
hx
=−
ds
(2.46)
Tb,x − Tw
ṁcp
Moltiplicando e dividendo il secondo membro per (µkDh ) ed esprimendo la portata
come ṁ = ρUav Hav
dTb,x
hx Dh 1 ds
=−
(2.47)
Tb,x − Tw
k Pe Hav
Nella (2.47) compare il numero di Nusselt locale, passando alla relazione integrale
ln
Tb,out − Tw
1
=−
Tb,in − Tw
PeHav
Z
w
Nux ds
(2.48)
Con la (2.44) ed esplicitando Nu si ottiene quanto cercato
Nu = ln σ
Hav
Dh
Pe = ln σ
Pe
2L
4L
(2.49)
Questa modalità di soluzione è stata testata affrontando la soluzione del campo
termico nel caso di canale piano parallelo. Secondo la soluzione analitica per questa
2.2 Definizione e risoluzione del problema termofluidodinamico
47
geometria, il numero di Nusselt per flusso pienamente sviluppato non dipende dal
numero di Reynolds e per temperature alle pareti costanti e uguali vale [12]
NuT = 7.54070087
(2.50)
Operando in forma dimensionale, come già notato, il coefficiente hx è definito come
hx =
q00x
Tb,x − Twall
(2.51)
dove q00x è il flusso termico specifico [ W/ m2 ].
In un canale generico si ha che
∃ (y1 , y2 ) | ∀x T (x, y1 ) = T (x, y2 ) = Tb,x
e in assenza di assialsimmetria ciascuna delle due pareti (di dimensione S1 e S2 rispettivamente) avrà un differente valore di Nu. Conseguentemente la potenza termica q
è
q = h(Tb,x − Tw )S = h1 (Tb,x − Tw )S1 + h2 (Tb,x − Tw )S2
(2.52)
dove S = S1 + S2 . Il coefficiente di scambio termico medio è quindi somma pesata dei
coefficienti medi a parete e lo stesso dicasi per il numero di Nusselt.
Figura 2.18 Visualizzazione del campo termico per un generico canale. La zona all’interno
delle due curve nere è caratterizzata da T (x, y) > Tb,x .
In figura 2.18 si riporta graficamente quanto detto per un generico canale. Sovrapposte al campo di temperatura, sono tracciate due curve che delimitano la zona, all’interno della quale la temperatura è maggiore di quella di bulk. Si osservi la periodicità
del profilo di temperatura θp : θp,in = θp,out .
3
Caratteristiche della procedura
risolutiva in FEMLAB
Nel capitolo precedente, sezione 2.2, si è data una descrizione dei passi necessari alla
risoluzione del problema di convezione forzata in un canale periodico bidimensionale.
In questa capitolo sono riportate alcune problematiche affrontate nell’implementazione di tale metodologia in FEMLAB, e dal punto di vista della geometria e dal punto
di vista numerico. FEMLAB è un solutore agli elementi finiti (FEM, Finite Element
Method) su griglia non-strutturata, si vedano a proposito appendice A.2 e [18]. Lo
scripting utilizzato da FEMLAB ne permette un diretto accoppiamento a MATLAB,
per cui è possibile interfacciare le funzionalità di entrambi i software. In questo modo
FEMLAB si rivela uno strumento molto potente e flessibile, molto più che nella sola
veste d’uso grafica GUI (Graphical User Interface).
3.1
Sezione geometrica
Durante la fase di test il programma principale ha subito diversi miglioramenti in
termini di controllo della geometria generata, in particolare riguardo la presenza di
più di un dominio e la caratteristica di stampabilità del profilo del canale.
3.1.1
Controllo di unicità del dominio
Come anche si vedrà in seguito, i.e. capitolo 4 sull’ottimizzazione, l’adozione di
NURBS per la generazione del canale ha come conseguenza un’enorme variabilità di
forme. Spesso però queste risultano patologicamente non processabili, perché nonsingolarmente connesse, caratterizzate cioè da un numero di domini maggiore ad uno
o prive di una adeguata sezione di inlet e outlet, perché FEMLAB non è stato in grado
di tagliarne i bordi. Due chiari esempi sono riportati nelle figure 3.1(a) e 3.1(b). Il
CAD interno a FEMLAB permette di riconoscere tali disfunzioni nelle geometrie,
per cui queste vengono scartate. Tuttavia una prima applicazione di tale controllo
portava a eliminare anche un canale come quello rappresentato in figura 3.2. Si tratta
49
50
Caratteristiche della procedura risolutiva in FEMLAB
(a) Geometria multi-connessa.
(b) Inlet e outlet definiti non correttamente.
Figura 3.1 Esempi di canali non processabili. Tratteggiata la NURBS periodica con punti di
controllo.
però questa di una geometria singolarmente connessa, ottenuta con una traslazione
del profilo superiore, trasl, in direzione positiva. Si può osservare che, traslando i due
Figura 3.2 Esempio di canale 2D non processabile.
profili in verso opposto a quanto espresso da trasl lungo la direzione x, il dominio
diviene unico, la geometria processabile ed equivalente alla precedente dal punto di
vista termofluidodinamico.
3.1.2
Controllo di stampabilità del profilo
Già dopo i primi lanci di prova, l’ottimizzatore trovava in geometrie simili alla 3.3
caratteristiche performanti e per scambio termico e per attrito, tuttavia si tratta di
un risultato puramente accademico, vista l’antieconomicità, o meglio l’impossibilità,
3.1 Sezione geometrica
51
Figura 3.3 Il canale precedente è ora processabile dopo la traslazione dei due profili.
di realizzazione del canale. Infatti appare evidente come risulti difficile, se non impossibile, ottenere una delle due superfici tramite stampaggio. Questo procedimento
consiste nell’obbligare una lamiera metallica ad occupare una cavità ricavata fra due
stampi contrapposti fra loro; in realtà sarebbe forse più corretto parlare di imbutitura, che per definizione consiste nel costringere una lamiera piana ad assumere una
certa concavità. Ad ogni modo la conca non può allargarsi dalla base all’apice,pena
l’irrealizzabilità del pezzo.
Visto il forte impatto della stampabilità sulla geometria, lo script è stato quindi
addizionato di un controllo di realizzabilità, inteso tuttavia come solo esempio di vincolo, (constraint), di fabbricabilità. In questo lavoro, dal carattere accademico, non si
pretende di assicurare canali realizzabili, quanto più di fornire un approccio metodologico a questo evidente problema. Nei casi reali, anche in funzione delle particolari
tecnologie produttive, ulteriori vincoli andrebbero inseriti.
Si è proceduto come segue. Una volta generato il profilo inferiore, questo viene
interamente spazzolato lungo la direzione x da rette via via di diversa inclinazione,
corrispondente alla eventuale direzione di stampaggio/imbutitura. Se esiste almeno
una inclinazione per cui tutte le rette intersecano il profilo in un unico punto, allora
e solo allora il profilo è stampabile. Le direzioni di stampaggio scelte sono descritte
da un angolo rispetto la direzione y compreso tra −45 ◦ e +45 ◦ . In figura 3.4 sono
rappresentate le direzioni ψ = −25 ◦ e ψ = 35 ◦ e in entrambi i casi esistono posizioni in cui le rette intersecano per più di una volta il profilo. Si è deciso di estrarre
comunque il coefficiente di attrito e il numero di Nusselt in canali di questo genere,
dichiarandoli poi geometrie unfeasible, ovvero irrealizzabili, ma utili nel successivo
processo di ottimizzazione. Un altro motivo, che ha spinto a conservare questi canali,
è la possibilità di renderli realizzabili con una lieve modifica alla geometria, tale cioè
da consentirne lo stampaggio e non influire sulle caratteristiche termofluidodinamiche
52
Caratteristiche della procedura risolutiva in FEMLAB
1.5
1
y
0.5
0
−0.5
−1
−1.5
−1
−0.5
0
0.5
x
1
1.5
2
Figura 3.4 Esempio di canale 2D non ottenibile tramite stampaggio.
del canale. Un esempio è il profilo riportato in figura 3.5. É sufficiente correggere
il tratto di curva indicato dalla freccia per ottenere un canale stampabile, senza per
questo alterare significativamente le prestazioni del canale. La verifica di quanto detto
è riportata in figura 3.6, dove sono riportati per ciascun canale il modulo di velocità
per il campo di velocità e le curve a temperatura costante per il campo termico.
Correzione
Figura 3.5 Esempio di profilo ottenibile tramite stampaggio con una lieve modifica alla
geometria: profilo non-marcato (non stampabile), profilo marcato (stampabile).
3.2 Sezione numerica
(a) Esempio di geometria non stampabile.
53
(b) Con una lieve correzione si ottiene un profilo stampabile senza penalizzare sensibilmente le
prestazioni.
Figura 3.6 Confronto fra le prestazioni dei due canali, il cui profilo è riportato in figura 3.5.
3.2
Sezione numerica
FEMLAB è un solutore agli elementi finiti (FEM). Il generico problema [19], descritto
da un sistema di equazioni differenziali alle derivate parziali (PDE), non viene risolto
calcolando le variabili negli infiniti punti del suo dominio, bensı̀ in un numero finito
di punti, nodi, del dominio. L’approssimazione del continuo in un numero finito di
parametri incogniti, o gradi di libertà (DOFs), è detta discretizzazione. Le modalità
di variazione dei parametri fra i vari nodi sono descritte dalle cosiddette funzioni di
forma, che esprimono la possibile forma della soluzione approssimata. L’applicazione
di questo metodo numerico richiede la seguente sequenza di operazioni:
1. il punto di partenza è la generazione di una mesh del dominio, ovvero la discretizzazione della geometria attraverso un certo numero di piccole unità geometriche di forma semplice;
2. ottenuta la griglia di calcolo, si introducono le approssimazioni delle variabili,
ovvero le funzioni di forma;
3. le PDE vengono scritte per ciascun elemento in forma integrale, ottenendo una
equazione matriciale;
4. le diverse equazioni create vengono assemblate fino a formare un sistema di
equazioni lineare, che con diversi metodi (diretti o iterativi) viene risolto fornendo i valori nodali delle incognite.
In questa tesi il modello è affrontato con dominio e mesh unici. In figura 3.7 è
presentato un esempio di mesh triangolare non-strutturata su un canale 2D.
54
Caratteristiche della procedura risolutiva in FEMLAB
Figura 3.7 Esempio di mesh.
3.2.1
Parte fluidodinamica.
Si è imposto che le funzioni di forma, utilizzate per descrivere l’andamento di una
variabile all’interno del singolo elemento, siano di tipo parabolico per la velocità e
per la temperatura, mentre per la pressione di tipo lineare. È infatti opportuno che la
pressione sia descritta all’interno dell’elemento da una funzione continua di classe C1 ,
ed in generale di ordine inferiore alla velocità. Si evita cosı̀ il pericolo che, a causa del
debole accoppiamento fra il campo delle velocità e quello della pressione, oscillazioni
sulla pressione non abbiano influenza sulla velocità [19]. Questo modo di procedere
è detto unequal order interpolation e corrisponde all’adozione di una staggered grid
nell’analisi ai volumi finiti.
È poi evidente come il numero di Reynolds, a cui viene studiato il canale, dia un
considerevole peso alla parte non-lineare delle equazioni di conservazione della quantità di moto, (2.7a) e (2.7b). Partendo con un tale valore di Re, il solver non-lineare
rischia di non andare in convergenza. La parte fluidodinamica è quindi stata suddivisa
in due fasi. Per ciascuna geometria nella prima fase il valore del numero di Reynolds
è portato da 10 a 200 con un conveniente valore del gradiente pressorio, a mezzo di
successivi restart, ossia di volta in volta la soluzione del campo fluido diventa punto di
partenza per la risoluzione successiva. Nella seconda l’ultimo campo di velocità trovato diventa punto di partenza per la determinazione del valore del gradiente pressorio
che dia la portata desiderata.
3.2.2
Parte termica.
Il problema termico di convezione forzata, al contrario di quello di moto, è tipicamente
lineare, essendo noto il campo di velocità all’interno del canale. Tuttavia l’utilizzo del
solutore lineare ha dato luogo a sporadiche oscillazioni, del tutto risolte con l’adozione
di un solutore non-lineare. Non deve stupire l’utilizzo di un algoritmo non-lineare per
55
3.2 Sezione numerica
problemi lineari. Infatti si può pensare il problema lineare, come un non-lineare in cui
il termine diffusivo (vedi (2.8)) è elevato e al limite infinito [20].
3.2.3
Discretizzazione del dominio.
La griglia utilizzata per il caso 2D è non-strutturata, del tipo triangolare a 3 nodi.
Questa è l’unica soluzione offerta dalla versione di FEMLAB utilizzata. Una griglia
di questo tipo è meno accurata di una ad elementi quadrangolari e a parità di nodi
e dominio richiede un maggior numero di elementi, tuttavia riesce a seguire meglio
l’andamento di geometrie complesse. Rispetto a lavori precedenti [5, 10] non si è
potuto utilizzare il solutore adattativo, ovvero un solutore di tipo iterativo capace di
comandare il raffinamento della griglia in zone in cui l’errore locale supera un certo
valore prescelto dall’utente. Infatti questa opzione non avrebbe permesso di effettuare
i restart. Si è quindi cercato di fissare i parametri di generazione della mesh in modo
da ottemperare mediamente alle richieste di qualità del maggior numero di domini,
senza però eccedere con le risorse di calcolo. FEMLAB permette di impostare un
certo numero di parametri globalmente e localmente su domini (se ve n’è più di uno!),
bordi e punti. Si riassume di seguito le scelte effettuate, tabella 3.1.
Parametro
Dimensione massima globale dell’elemento
Tasso di incremento di dimensione
Fattore di curvatura
Fattore di troncamento della curvatura
Dimensione massima dell’elemento su inlet e outlet
Dimensione massima dell’elemento su wall
Tasso di incremento di dimensione su inlet e outlet
Valore
1/16 Hav
1.08
0.4
0.04
1/20 Hav
1/18 Hav
1.04
Tabella 3.1 Imposizione parametri per la generazione della mesh.
Il fattore di curvatura della mesh fornisce una indicazione circa la dimensione
dell’elemento in presenza di raggi di curvatura; il fattore di troncamento invece evita
che in presenza di spigoli si addensi un elevato numero di elementi. Tutto ciò che
non viene dichiarato è posto da FEMLAB al valore di default e per maggiori dettagli
si rimanda a [18]. La fase di generazione della mesh si è rivelata critica qualora le
dimensioni dei bordi non risultassero dimensionalmente simili. In queste zone si è
riscontrato un infittimento del numero degli elementi non necessario e che talvolta ha
portato mesh non coerenti e geometrie non processabili. Una stima del numero medio
di elementi cosı̀ ottenuti si aggira attorno a 4500, con un tempo medio di risoluzione
per canale pari a 6 minuti di CPU time su processore INTEL 2GHz.
4
Ottimizzazione 2D
L’insieme delle procedure precedentemente descritte diventa ora il motore per la risoluzione del problema di ottimizzazione geometrica. L’ottimizzazione di un sistema è
la scelta, o determinazione, di opportuni valori per una serie di parametri, con cui il
sistema è descritto, in modo tale che siano soddisfatte specifiche richieste, obiettivi.
Le strategie utilizzabili per ottimizzare sono svariate [21]. In questo lavoro è richiesta
una procedura il più possibile generale, vista l’estrema variabilità lasciata alla forma
geometrica del canale.
4.1
Il processo di ottimizzazione
Il software modeFRONTIER, si vedano appendice A.1 per una breve descrizione e
[22], racchiude in forma organica diversi algoritmi di ottimizzazione e ne consente
una agevole applicabilità al problema in questione, i.e. ottimizzazione di forma di
canali di scambiatori di calore.
Le grandezze che caratterizzano la performance del canale sono il coefficiente
d’attrito adimensionale e il numero di Nusselt. Entrambe sono in funzione dei parametri geometrici, ma la relazione non è nota. Detto X l’insieme delle variabili
geometriche, generate da modeFRONTIER, e φ l’insieme degli obiettivi da massimizzare/minimizzare, l’ottimizzatore valuta la seguente espressione [10]
φ = g(X)
(4.1)
[ f , Nu] = g (L, Lwave , ∆x3 , ρ4 , ϑ4 , ∆x5 , ∆y5 , ρ6 , ϑ6 , ∆x7 ,trasl)
(4.2)
e in particolare la
attraverso quella che è detta application. L’applicazione in uso è FEMLAB che, accoppiato a MATLAB, esternamente carica i parametri forniti da modeFRONTIER e
ritorna ad esso le caratteristiche prestazionali del canale.
Nella pagine che seguono, sono riportati i risultati dei diversi lanci dell’ottimizzatore. Si parte con un algoritmo genetico, MOGA-II, molto robusto, per sondare
57
58
Ottimizzazione 2D
l’intero campo delle soluzioni, sfruttando anche la sua peculiarità di essere multiobiettivo. In seguito con la definizione della utility function, il problema è trasformato
in mono-obiettivo per provare un algoritmo più accurato ed efficiente, quale SIMPLEX. Infine verrà presentato un confronto fra l’ottimizzazione dei canali descritti
da profili NURBS e canali descritti da profili rettilinei (linear-piecewise). La maggior complessità geometrica, apportata dalla parametrizzazione NURBS, è ripagata
da prestazioni evidentemente migliori.
4.2
Prima fase: ottimizzazione MOGA multi-obiettivo
Tra gli algoritmi disponibili, quello che meglio si presta ad una prima indagine della relazione (4.2) è l’algoritmo genetico (GA), MOGA-II (Multi Objective Genetic
Algorithm) in modeFRONTIER. Come meglio riportato in appendice A.1 questo algoritmo, oltre ad essere multi-obiettivo, è molto robusto. Pur pagando una minore
efficienza quindi, non rischia di adagiarsi su massimi o minimi locali, avvallamenti
da cui scheduler più efficienti e accurati non riuscirebbero a discostarsi. L’idea alla
base di questo algoritmo è l’evoluzione genetica di una popolazione iniziale, D.O.E.
(Design of Experiments), attraverso le operazioni di cross-over, mutazione, selezione,
appendice A.1.
Per ottenere una buona uniformità all’interno dello spazio delle soluzioni per la
generazione del D.O.E., affidata a modeFRONTIER, si è utilizzato l’algoritmo Sobol. Come in [22] per ottenere risultati soddisfacenti il numero di individui della
generazione di partenza NDOE deve soddisfare la seguente relazione empirica
NDOE > 2 NDOF Nobj
(4.3)
dove NDOF è il numero delle variabili libere (nel nostro caso 11), mentre Nobj è il
numero degli obiettivi (nel nostro caso 2).
4.2.1
Primo lancio
Per mantenere un certo margine di sicurezza si è generata una popolazione iniziale di
50 individui. Il numero di generazioni testate è pari a 15, per un totale di 750 individui
processati. Infine si sono impostati dei vincoli, constraints, per meglio guidare l’ottimizzatore. Siano Nu0 il numero di Nusselt e f0 il coefficiente d’attrito per il canale
piano, dove f0 = 24/ Re [12]: per Re = 200 f0 = 0.12. Durante questo primo lancio
si è imposto che:
59
4.2 Prima fase: ottimizzazione MOGA multi-obiettivo
Vincolo
Nu/ Nu0
f / f0
Relazione
>
<
Limite
1
10
Tabella 4.1 Constraints per il primo lancio MOGA.
In aggiunta a ciò il canale deve essere realizzabile tramite stampaggio o imbutitura.
Se uno di questi vincoli viene violato, il particolare design è dichiarato unfeasible
e viene penalizzato dall’ottimizzatore, come fosse un individuo a basse prestazioni.
Avrà quindi scarse possibilità di trasmettere i suoi geni alle generazioni future.
In tabella 4.2 sono riportati i campi di variazione delle variabili libere, assieme alle
relative basi, ovvero il numero di passi discreti con cui si suddivide detto intervallo.
Per maggiore chiarezza si riporta quanto indicato nelle tabelle 2.2 e 2.3 completando
i valori degli estremi degli intervalli. Si osservi che gli angoli sono riportati in gradi e
Parametro
Campo
Base
L
[1.0; 2.5]
1001
Lwave
[0.15 L; 0.85 L]
1001
∆x3
[0.01 L; 0.20 L]
1001
ρ4
[0.05; 1.20]
1001
ϑ4
[50 ◦ ; 80 ◦ ]
501
∆x5
[0.2 L; 0.7 L]
1001
∆y5
[0.0; 0.5]
1001
ρ6
[0.05; 1.10]
1001
ϑ6
[−100 ◦ ; 140 ◦ ]
501
∆7
[0.01 L; 0.75 L]
1001
trasl
[−0.499 L; 0.500 L]
101
Tabella 4.2 Dominio di definizione delle variabili per il primo lancio MOGA.
sono misurati rispetto la direzione positiva dell’asse delle ascisse. L’angolo è positivo
se la rotazione è antioraria. Seguendo di pari passo il modo in cui modeFRONTIER
importa gli estremi degli intervalli, questi ultimi sono tutti chiusi e con il valore 0.01
si vuole intendere un valore prossimo a 0.
Si può notare come rispetto al lavoro in [5], grazie alla maggior potenza di calcolo
disponibile, sia stato possibile incrementare i gradi di libertà del modello, cha passano
da 10 a 11. In particolare, viene lasciata libera di variare la lunghezza del modulo
L. Tale parametro riveste un ruolo importante nella pesantezza del calcolo. Infatti,
60
Ottimizzazione 2D
poiché la densità di mesh è definita in funzione dell’altezza media Hav (tabella 3.1)
che è costante, al crescere della lunghezza del modulo, ne aumenta proporzionalmente
la superficie e quindi il numero di elementi della griglia.
Considerando le usuali tecniche, l’incremento dello scambio termico in canali di
scambiatori di calore è ottenibile (a parità di altre condizioni)
- aumentando la superficie di scambio fluido-parete (i.e. mediante corrugazioni)
- aumentando il coefficiente di scambio termico attraverso promotori di vorticità
La prima delle due possibilità porta a chiari aumenti del coefficiente d’attrito, dovuti
ad una maggiore superficie esposta al cosiddetto wall shear stress. Con la seconda
invece si può ottenere un maggior scambio termico grazie al mescolamento di flussi appartenenti alla zona di bulk (i.e. più caldi) con flussi appartenenti alla zona di
wall (i.e. più freddi), limitando i fenomeni dissipativi comunque presenti in un flusso
rotazionale. Per ottenere vortici di questa specie è richiesta una opportuna configurazione geometrica. Si è quindi pensato che lasciare ampia variabilità alla lunghezza
del canale L potesse favorire il conseguimento di configurazioni tali da promuovere lo
sviluppo dei vortici e soprattutto il fenomeno del miscelamento. Si vuole qui ripetere
che il flusso è laminare, si parla quindi di vorticità e non di turbolenza. La turbolenza
è infatti caratterizzata da un mescolamento diffusivo dei filetti fluidi che vengono a
perdere la loro identità; la vorticità in regime laminare è caratterizzata dal moto di
particelle che continuano a seguire ciascuna la propria linea di corrente [11].
Il risultato di questo primo lancio per i soli individui feasible, ovvero tali da rispettare i vincoli, è presentato in figura 4.1. Come già evidenziato i due obiettivi sono
Figura 4.1 Grafico del coefficiente adimensionale d’attrito f in funzione del coefficiente di
scambio termico adimensionale Nu per il primo lancio MOGA.
contrastanti per cui il fronte di Pareto, ovvero l’insieme delle soluzioni migliori (si
veda A.1.4), assume l’andamento di figura 4.2. Osservando le figure 4.1 e 4.2, appare evidente come i design si addensino sul valore minimo del Nu imposto, salendo
con Nu gli individui invece si diradano. Per validare i risultati ottenuti si è scelto di
eseguire una verifica di grid independence su una dei migliori individui del fronte. Il
4.2 Prima fase: ottimizzazione MOGA multi-obiettivo
61
Figura 4.2 Grafico di f in funzione di Nu per i soli individui dominanti del primo lancio
MOGA.
design scelto è il 739, si veda figura 4.3. Questo test consiste nel risolvere il campo
(a) ID 739
(b) mesh su ID 739
Figura 4.3 Individuo 739 su cui si è eseguita la grid independence.
termofluidodinamico all’interno del canale con una mesh via via migliore: se il coefficiente di attrito e il numero di Nusselt trovati tendono asintoticamente ad un valore
costante, i risultati risentono solo minimamente dell’errore dovuto alla discretizzazione. In tabella 4.3 sono riportati in ordine crescente il numero di elementi di cui la
mesh è costituita e i risultati estratti: i valori evidenziati in grassetto sono ottenuti con
la griglia di cui in tabella 3.1 e rappresentata in figura 4.3(b). Nelle figure 4.4(a) e
4.4(b) si presenta graficamente la tabella 4.3.
62
Ottimizzazione 2D
Numero di elementi
2015
2899
4280
4888
5734
8610
9897
12498
17120
Nu
13.427
13.296
13.273
13.180
13.052
13.314
13.227
13.202
13.205
f
1.206
1.210
1.199
1.198
1.208
1.204
1.204
1.202
1.202
13.7
1.33
13.6
1.29
13.5
fattore d’attrito
numero di Nusselt
Tabella 4.3 Grid Independence Test su ID 739.
13.4
13.3
13.2
1.21
1.17
1.13
13.1
1.09
13
12.9
2000
1.25
4000
6000
8000 10000 12000 14000 16000 18000
numero di elementi
(a) Grid Independence Test per Nu
1.05
2000
4000
6000
8000 10000 12000 14000 16000 18000
numero di elementi
(b) Grid Independence Test per f
Figura 4.4 Grid Independence Test per l’individuo 739.
Dall’analisi delle due figure precedenti, si può ritenere sufficientemente attendibile
il risultato ottenuto con la griglia generata con i parametri di tabella 3.1.
4.2.2
Secondo lancio
Il primo lancio dell’ottimizzatore è destinato a valutare la bontà dei campi di variazione scelti per le variabili libere. Dall’analisi delle history chart [22] di ciascuna
variabile è quindi possibile restringere i vari domini, in modo da consentire un’analisi più accurata. Le basi sono state modificate in modo da lasciare inalterati gli step
di ciascuna variabile, con l’unica eccezione per la variabile trasl, il cui passo è stato
ridotto di un fattore 10, aumentando la base da 101 a 1001; in tabella 4.4 vengono
4.2 Prima fase: ottimizzazione MOGA multi-obiettivo
63
presentati solo gli estremi inferiore e superiore. Pensando che, per un aumento del
Parametro
Campo
L
[1.0; 2.65]
Lwave
[0.29 L; 0.78 L]
∆x3
[0.012 L; 0.185 L]
ρ4
[0.021; 1]
ϑ4
[100 ◦ ; 300 ◦ ]
∆x5
[0.2 L; 0.7 L]
∆y5
[0.0; 0.4]
ρ6
[0.05; 0.9]
ϑ6
[−120 ◦ ; 120 ◦ ]
∆7
[0.01 L; 0.65 L]
trasl
[−0.499 L; 0.500 L]
Tabella 4.4 Dominio di definizione delle variabili per il secondo lancio MOGA.
vincolo sul numero di Nusselt, si dovesse accettare una relativa maggiorazione delle
perdite di carico, si sono ridefiniti i constraints, tabella 4.5. Il D.O.E. è ancora compoVincolo
Nu/ Nu0
f / f0
Relazione
>
<
Limite
1.13
15
Tabella 4.5 Constraints per il secondo lancio MOGA.
sto da 50 individui, di cui 34 appartenenti al fronte di Pareto del precedente lancio, 1
unfeasible e i restanti 15 ottenuti con l’algoritmo Sobol sui nuovi campi di variazione.
In questo modo l’algoritmo di ottimizzazione è indirizzato già in partenza su individui
prestazionali; si è scelto di includere un individuo unfeasible perché comunque caratterizzato da un buon valore di scambio termico. Il numero di generazioni è pari a 15,
per un totale di 750 individui.
Rispetto al primo lancio il numero di geometrie che FEMLAB non riesce a caricare si è notevolmente ridotto, segno che la riduzione dei campi di variazione aiuta la
generazione di individui processabili. A questo si deve aggiungere l’opera di selezione a carico dello scheduler, che fa in modo di non permettere la riproduzione di error
design.
64
Ottimizzazione 2D
Come per il lancio precedente si riportano i grafici che consentono una agevole,
sia pur qualitativa, lettura dei risultati ottenuti. In figura 4.5 sono riportati insieme ai
design dell’attuale ottimizzazione (evidenziati), anche i risultati del precedente lancio; il fronte di Pareto e in generale l’ideale curva, che interpola i risultati, si sono
spostati verso destra e in basso, ovvero verso valori maggiori per quanto concerne lo
scambio termico e inferiori per l’attrito. Il fronte di Pareto è piuttosto scarno. Ciò
Figura 4.5 Grafico del coefficiente adimensionale d’attrito f in funzione del coefficiente di
scambio termico adimensionale Nu per il secondo lancio MOGA. In evidenza i design attuali.
si ritiene dovuto al non avere concesso allo scheduler un maggior numero di generazioni: l’andamento è in figura 4.6. Vi è un certo numero di design isolati rispetto
Figura 4.6 Grafico di f in funzione di Nu per i soli individui dominanti del secondo lancio
MOGA.
alla tendenza generale, sia in riferimento al coefficiente d’attrito, sia in riferimento al
numero di Nusselt. Tuttavia mentre canali ad attrito elevato, giustificano un considerevole scambio termico, canali con Nu superiore a 13 e f nella norma hanno destato
un certo sospetto. Analizzando tali individui separatamente si è notata, in linea con le
aspettative, una mancata convergenza della sezione termica. È stato infatti sufficiente
65
4.2 Prima fase: ottimizzazione MOGA multi-obiettivo
modificare leggermente la griglia per riportare il numero di Nusselt a valori più credibili. Come precedentemente spiegato al paragrafo 2.2.10, il corretto rapporto fra le
temperature di bulk tra ingresso e uscita del canale è ottenuto tramite la valutazione
del flusso termico adimensionale sulle pareti. Questo processo è molto sensibile alla qualità della mesh utilizzata, poiché entra in gioco il gradiente della temperatura
adimensionale, il cui calcolo necessita oltre che del valore della funzione nel nodo
i − esimo, anche dello stesso nei nodi adiacenti. Una configurazione distorta porta a
prevedivili errori numerici. Altre volte, anche se per fortuna più raramente, si è osservata una netta differenza sempre nel calcolo di Nu per il medesimo canale tra un
lancio su macchina LINUX e lo stesso su macchina WINDOWS. In questo frangente
concorrono il diverso trattamento degli errori nei due sistemi operativi ed eventuali
disuniformità tra gli stessi applicativi nelle due differenti versioni.
4.2.3
Terzo lancio
Il miglioramento nella popolazione apportato dal secondo lancio MOGA non è stato
molto rilevante, probabilmente a causa dell’incremento del vincolo sul coefficiente
d’attrito. Questa volta si opera in direzione opposta: il vincolo sullo scambio termico
viene ridotto di un punto percentuale, mentre quello su f diviene quasi la metà, tabella
4.6. Il D.O.E. è composto sempre da 50 individui, di cui 38 appartenenti al fronte di
Vincolo
Relazione
Limite
Nu/ Nu0
>
1.12
f / f0
<
8
Tabella 4.6 Constraints per il terzo lancio MOGA.
Pareto del precedente lancio e i restanti 12 ottenuti con l’algoritmo Sobol. I nuovi
intervalli di definizione per le variabili libere sono presentati in tabella 4.7. Il numero
di generazioni è pari a 10, per un totale di 500 nuovi individui. Per maggior chiarezza
si riportano i nuovi risultati assieme a quelli del solo lancio precedente: la verifica
dell’eventuale miglioramento ottenuto è immediata. I design dell’attuale lancio sono
posti in evidenza. Il rinnovato vincolo sul coefficiente d’attrito ha dato un certo
effetto, osservabile nella traslazione del gruppo verso il basso. Per quanto riguarda
Nu invece non si ha una sostanziale modifica. Non ci si deve infatti lasciar fuorviare
dalla presenza di individui all’estremità inferiore destra di figura 4.7, caratterizzati da
una mancata convergenza della parte termica dell’algoritmo risolutivo.
66
Ottimizzazione 2D
Parametro
Campo
L
[0.91; 2.65]
Lwave
[0.35 L; 0.78 L]
∆x3
[0.020 L; 0.180 L]
ρ4
[0.020; 0.80]
ϑ4
[110 ◦ ; 310 ◦ ]
∆x5
[0.10 L; 0.65 L]
∆y5
[0.0; 0.38]
ρ6
[0.05; 0.91]
ϑ6
[−151 ◦ ; 87 ◦ ]
∆7
[0.01 L; 0.51 L]
trasl
[−0.499 L; 0.500 L]
Tabella 4.7 Dominio di definizione delle variabili per il terzo lancio MOGA.
Figura 4.7 Grafico del coefficiente adimensionale d’attrito f in funzione del coefficiente di
scambio termico adimensionale Nu per il terzo lancio MOGA. In evidenza i design attuali.
4.3
Seconda fase: eliminazione di un obiettivo
Prima di passare all’ottimizzazione mono-obiettivo con l’uso delle utility function, si
è valutato il comportamento di MOGA-II con l’unica richiesta prima di massimizzare
Nu e poi di minimizzare f .
4.3.1
MOGA-II: massimizzazione di scambio termico
In questo lancio si passa dalla configurazione bi-obiettivo, alla mono-obiettivo. L’obiettivo di minimizzare il fattore d’attrito medio viene sostituito da un vincolo sul-
67
4.3 Seconda fase: eliminazione di un obiettivo
Figura 4.8 Grafico di f in funzione di Nu per i soli individui dominanti del terzo lancio
MOGA.
lo stesso (tabella 4.8). L’unico obiettivo rimane dunque rendere massimo lo scambio termico adimensionale e non risulta più necessario porre un vincolo sul rapporto
Nu/ Nu0 minimo, essendo Nu l’unico obiettivo. È mantenuta come ulteriore constraint la stampabilità del profilo. Gli estremi di variazione delle variabili libere sono
Vincolo
f / f0
Relazione
<
Limite
6
Tabella 4.8 Constraints per il quarto lancio MOGA mono-obiettivo.
riportati in tabella 4.9.
Figura 4.9 Grafico del coefficiente adimensionale d’attrito f in funzione del coefficiente di
scambio termico adimensionale Nu per il quarto lancio MOGA mono-obiettivo. In evidenza i
design attuali.
68
Ottimizzazione 2D
Parametro
L
Lwave
∆x3
ρ4
ϑ4
∆x5
∆y5
ρ6
ϑ6
∆7
trasl
Campo
[0.91; 2.65]
[0.40 L; 0.78 L]
[0.020 L; 0.180 L]
[0.020; 0.65]
[118 ◦ ; 320 ◦ ]
[0.20 L; 0.65 L]
[0.0; 0.38]
[0.1; 0.88]
[−151 ◦ ; 75 ◦ ]
[0.01 L; 0.41 L]
[−0.499 L; 0.440 L]
Tabella 4.9 Dominio di definizione delle variabili per il quarto lancio MOGA mono-obiettivo.
Figura 4.10 Andamento di Nu per i soli individui feasible del quarto lancio MOGA monoobiettivo. In evidenza i design attuali.
Non si osserva alcun sensibile miglioramento rispetto la configurazione precedente; i design con buone prestazioni sono quindi tutti derivati con piccole modifiche
dalle generazioni precedenti. Le figure 4.10 e 4.11 rappresentano le history charts
rispettivamente per il numero di Nusselt e per il coefficiente d’attrito.
4.3.2
MOGA-II: minimizzazione del fattore d’attrito
In questo lancio si rovescia la configurazione precedente, ovvero si assume come unico obiettivo la minimizzazione del fattore d’attrito, unito ad un vincolo sul valore
69
4.3 Seconda fase: eliminazione di un obiettivo
Figura 4.11 Andamento di f per i soli individui feasible del quarto lancio MOGA monoobiettivo. In evidenza i design attuali.
minimo del numero di Nusselt (tabella 4.10). Come nel caso precedente non risulta
più necessario porre un vincolo sul rapporto f / f0 massimo, essendo f l’unico obiettivo. È ancora mantenuta come ulteriore constraint la stampabilità del profilo. Gli
Vincolo
Nu/ Nu0
Relazione
>
Limite
1.55
Tabella 4.10 Constraints per il quinto lancio MOGA mono-obiettivo.
estremi di variazione delle variabili libere sono riportati in tabella 4.11.
70
Ottimizzazione 2D
Parametro
L
Lwave
∆x3
ρ4
ϑ4
∆x5
∆y5
ρ6
ϑ6
∆7
trasl
Campo
[1.0; 2.65]
[0.29 L; 0.78 L]
[0.020 L; 0.180 L]
[0.020; 0.80]
[101 ◦ ; 320 ◦ ]
[0.20 L; 0.71 L]
[0.0; 0.38]
[0.1; 0.88]
[−151 ◦ ; 75 ◦ ]
[0.01 L; 0.41 L]
[−0.499 L; 0.500 L]
Tabella 4.11 Dominio di definizione delle variabili per il quinto lancio MOGA
mono-obiettivo.
Rispetto al lancio precedente, per quel che riguarda la lunghezza del canale L e
la traslazione del profilo superiore trasl si è tornati ai campi del terzo lancio (tabella
4.7).
Figura 4.12 Grafico del coefficiente adimensionale d’attrito f in funzione del coefficiente di
scambio termico adimensionale Nu per il quinto lancio MOGA mono-obiettivo. In evidenza i
design attuali.
4.4 Terza fase: ottimizzazione mono-obiettivo mediante utility function
71
Figura 4.13 Andamento di Nu per i soli individui feasible del quinto lancio MOGA monoobiettivo. In evidenza i design attuali.
Figura 4.14 Andamento di f per i soli individui feasible del quinto lancio MOGA monoobiettivo. In evidenza i design attuali.
I grafici delle history chart, figure 4.13 e 4.14, sono piuttosto radi per il gran numero di individui dichiarati unfeasible dal vincolo di stampabilità e dal valore minimo
del numero di Nusselt.
4.4
Terza fase: ottimizzazione mono-obiettivo mediante utility function
In questa sezione si mantiene la forma mono-obiettivo, finalizzata questa volta all’utilizzo di un algoritmo di ottimizzazione quale SIMPLEX, più efficiente di MOGAII, anche se meno robusto. Un ottimizzatore di questo tipo tratta però soltanto problemi mono-obiettivo e modeFRONTIER offre svariati algoritmi che permettono di
trasformare problemi multi-obiettivo, in problemi mono-obiettivo.
L’algoritmo in questione è detto M.C.D.M. (Multi Criteria Decision Making). Nei
72
Ottimizzazione 2D
problemi multi-obiettivo, non esiste un’unica soluzione eccellente, bensı̀ una serie
di casi che hanno caratteristiche tali da prevalere su altri: questi design si dicono
appartenere al fronte di Pareto. La funzione obiettivo, di tipo monotono, combina
sotto forma di somma pesata i diversi obiettivi di partenza(1) . Detti fi gli n obiettivi
multipli, la utility function Futility è
n
Futility = ∑ pi ( fi ) fi
(4.4)
i=1
dove con pi ( fi ) si intendono i pesi associati a ciascun obiettivo e variabili da obiettivo a obiettivo; i pesi sono caratterizzati dall’avere somma unitaria. L’algoritmo
M.C.D.M. è in grado di fornire i valori dei pesi. L’utente sceglie un certo numero
di individui ed esprime su di essi un giudizio, mettendoli a confronto su particolari
attributi che li caratterizzano o sui loro dati di input. Dei design scelti viene creato un
ranking, ovvero una scala di merito ottenuta dalla somma pesata dei vecchi obiettivi,
da cui ricavare i pesi. A questo punto allo scheduler è richiesto di rendere estrema la
funzione obiettivo.
4.4.1
Utility function con MOGA-II e SIMPLEX
Prima di effettuare un lancio con un algoritmo basato sul gradiente A.1, viene riproposto per la sua intrinseca robustezza l’algoritmo MOGA-II, in modo da verificare
l’opportunità nella scelta delle relazioni fra gli individui. La funzione obiettivo è stata
formata con gli individui e le relazioni riportati in tabella 4.12. La scelta delle relazioID
Relazione
ID
896
è equivalente a
2526
2861
è preferito a
2526
3027
è preferito a
3023
3027
è equivalente a
3029
Tabella 4.12 Relazioni fra gli individui per la generazione della utility function.
ni è legata al valore delle due funzioni obiettivo. In questa fase si è cercato di sondare
lo spazio delle soluzioni che privilegia lo scambio termico. Gli individui, con i quali
sono state create le relazioni di tabella 4.12, sono riportati in figura 4.15. Le geometrie sono ordinate in ordine crescente rispetto al nuovo valore assunto dalla funzione
obiettivo, si veda anche figura 4.16. A ciascun canale sono allegati il numero di design, il valore del coefficiente d’attrito, del numero di Nusselt e il valore normalizzato
a 1 della nuova funzione obiettivo (4.5). La sezione M.C.D.M. consente di utilizzare
(1) Lo
stesso procedimento si può realizzare sulle variabili di input.
4.4 Terza fase: ottimizzazione mono-obiettivo mediante utility function
(a) Rank value 0.283
73
(b) Rank value 0.289
(c) Rank value 0.556
(d) Rank value 0.717
(e) Rank value 0.966
(f) Rank value 0.984
Figura 4.15 Geometrie degli individui utilizzati per la definizione della funzione obiettivo.
74
Ottimizzazione 2D
diversi sotto-algoritmi per generare le utility function, tra cui Linear MCDM e Genetic
Algorithm MCDM. Il primo di questi è da utilizzarsi quando il numero delle variabili
decisionali è inferiore a 4. In questo caso i due metodi hanno fornito i medesimi risultati. Le figure 4.16 e 4.17 rappresentano rispettivamente lo scalaggio degli individui
secondo le relazioni fornite e i pesi forniti dall’algoritmo. La funzione obiettivo è
Figura 4.16 Relazioni di preferenza e ranking dei design scelti.
Figura 4.17 Pesi attribuiti a ciascun obiettivo di partenza.
quindi
Futility = 0.717 Nuobj + 0.283 fobj
(4.5)
Sono state effettuate alcune modifiche ai campi di variazione delle variabili, soprattutto per quanto riguarda l’entità della traslazione del profilo superiore, tabella 4.13. Il
D.O.E. è composto da 28 individui, di cui 18 appartenenti al processo di ottimizzazione precedente. I restanti 10 sono stati generati con Sobol. Il numero di generazioni è
pari a 15 per un totale di 420 individui.
4.4 Terza fase: ottimizzazione mono-obiettivo mediante utility function
Parametro
Campo
L
[1.0; 2.65]
Lwave
[0.28 L; 0.78 L]
∆x3
[0.020 L; 0.180 L]
ρ4
[0.020; 0.80]
ϑ4
[101 ◦ ; 329 ◦ ]
∆x5
[0.20 L; 0.71 L]
∆y5
[0.0; 0.38]
ρ6
[0.1; 0.88]
ϑ6
[−171 ◦ ; 75 ◦ ]
∆7
[0.01 L; 0.41 L]
trasl
[−0.499 L; 0.300 L]
75
Tabella 4.13 Dominio di definizione delle variabili per il lancio con utility function.
Sono stati eliminati i vincoli su Nu e su f , mentre rimane il vincolo sulla stampabilità. Dall’analisi dei canali ottenuti in precedenza ci si è accorti che lo scheduler
manteneva in vita geometrie performanti caratterizzate però da curve spigolose. Si è
quindi introdotto un ulteriore vincolo sulla curvatura delle pareti: se questa è inferiore
a 1/230, la geometria è dichiarata unfeasible.
Dal fronte di Pareto ottenuto si è creata la popolazione iniziale per una ottimizzazione tramite SIMPLEX, utilizzando la stessa funzione di peso. Il numero di design
iniziali è pari al numero di variabili più una, ovvero 12 [22]. I risultati trovati sono
presentati per i soli design appartenenti al lancio attuale.
Figura 4.18 Grafico del coefficiente adimensionale d’attrito f in funzione del coefficiente di
scambio termico adimensionale Nu a seguito del lancio SIMPLEX.
76
Ottimizzazione 2D
Figura 4.19 Grafico di f in funzione di Nu per i soli individui dominanti a seguito del lancio
SIMPLEX.
Figura 4.20 Andamento di Nu per i soli individui feasible a seguito del lancio SIMPLEX.
Figura 4.21 Andamento di f per i soli individui feasible a seguito del lancio SIMPLEX.
Dalle figure 4.20, 4.21 e 4.22 si osserva come dopo appena 174 interazioni l’algoritmo abbia ottenuto la posizione estrema per la utility function, la quale tende asinto-
4.4 Terza fase: ottimizzazione mono-obiettivo mediante utility function
77
Figura 4.22 Andamento della utility function per i soli individui feasible a seguito del lancio
SIMPLEX.
ticamente al valore massimo normalizzato all’unità. Una volta raggiunta questa conca
l’algoritmo non vi si solleva più, da cui il rischio di trovare estremi locali (mancanza
di robustezza).
4.4.2
MOGA-II con differente Utility function
I risultati appena trovati con il precedente procedimento non si prestano ad un agevole confronto con [5], poiché in genere caratterizzati da valori di Nu e assieme f
maggiori. Gli ultimi due lanci effettuati utilizzano funzioni mono-obiettivo con un
peso maggiore sulla riduzione del fattore d’attrito. In questo modo si è cercato di
individuare geometrie in linea con i risultati di [5].
UT1
La nuova utility function è stata calcolata con gli individui e le relazioni di tabella 4.14
Gli individui, con i quali sono state create le relazioni di tabella 4.14, sono riportati
ID
2526
3023
2595
3023
Relazione
è preferito a
è preferito a
è preferito a
è equivalente a
ID
2529
2526
3526
2829
Tabella 4.14 Relazioni fra gli individui per la generazione della utility function UT1.
in figura 4.23. Le geometrie sono ordinate in ordine crescente rispetto il nuovo valore
assunto dalla funzione obiettivo. La funzione obiettivo creata sempre dall’algoritmo
78
Ottimizzazione 2D
(a) Rank value 0.333
(b) Rank value 0.666
(c) Rank value 0.667
(d) Rank value 0.683
(e) Rank value 1
Figura 4.23 Geometrie degli individui utilizzati per la definizione della funzione obiettivo
UT1.
79
4.4 Terza fase: ottimizzazione mono-obiettivo mediante utility function
lineare è
Futility = 0.667 Nuobj + 0.333 fobj
(4.6)
Ci si aspetta un fronte di Pareto traslato verso il basso e a sinistra.
Figura 4.24 Grafico del coefficiente adimensionale d’attrito f in funzione del coefficiente di
scambio termico adimensionale Nu con UT1.
Figura 4.25 Grafico di f in funzione di Nu per i soli individui dominanti con UT1.
80
Ottimizzazione 2D
Figura 4.26 Andamento di Nu per i soli individui feasible.
Figura 4.27 Andamento di f per i soli individui feasible con UT1.
UT2
Si è provato a spostare ancora verso bassi valori del coefficiente d’attrito il fronte di
Pareto, definendo una nuova utility function. Finora le funzioni obiettivo sono state
create con pochi individui nel ranking. In questo lancio il numero di individui selezionati per un possibile confronto è maggiore. L’algoritmo di generazione è questa volta
di tipo genetico. Gli individui, con i quali sono state create le relazioni di tabella 4.15,
sono riportati in figura 4.28. Le geometrie sono ordinate in ordine crescente rispetto il
nuovo valore assunto dalla funzione obiettivo, si veda anche figura 4.29. La funzione
obiettivo, creata mediante l’algoritmo genetico GA, è
Futility = 0.336 Nuobj + 0.664 fobj
(4.7)
Il ranking 4.29 degli individui prescelti è decisamente più cospicuo dei precedenti, si
veda ad esempio 4.16. Questo dovrebbe aiutare il processo di ottimizzazione basato
sulla funzione di peso.
4.4 Terza fase: ottimizzazione mono-obiettivo mediante utility function
(a) Rank value 0.830
(b) Rank value 0.880
(c) Rank value 0.880
(d) Rank value 0.950
(e) Rank value 0.882
81
(f) Rank value 0.941
(g) Rank value 0.954
Figura 4.28 Geometrie degli individui utilizzati per la definizione della funzione obiettivo con
UT2.
82
Ottimizzazione 2D
ID
Relazione
ID
3054
è preferito a
1391
3197
è preferito a
1452
1391
è preferito a
3729
2571
è preferito a
1553
Tabella 4.15 Relazioni fra gli individui per la generazione della utility function UT2.
Figura 4.29 Ranking dei design scelti per UT2.
Ancora ci si aspetta un fronte di Pareto traslato verso il basso e a sinistra, poiché
si indirizza lo scheduler verso quella zona dello spazio delle soluzioni.
Figura 4.30 Grafico del coefficiente adimensionale d’attrito f in funzione del coefficiente di
scambio termico adimensionale Nu con UT2.
4.5 Analisi dei Risultati
83
Figura 4.31 Grafico di f in funzione di Nu per i soli individui dominanti con UT2.
4.5
Analisi dei Risultati
Come più volte ripetuto, si sono trovate diverse forme dominanti, a seconda di quale
dei due aspetti, i.e. f o Nu, si sia voluto valorizzare. Riassumendo, quattro sono le
tipologie di canale che si è ritenuto scegliere, e in particolare una, identificata con la
lettera c, simile in forma e prestazioni a quella determinata nel lavoro di [5]. La figura 4.32(a) riporta il fronte di Pareto risultante dall’intero processo di ottimizzazione.
Quattro canali, rappresentanti le suddette tipologie, sono evidenziati. La figura 4.32(b)
presenta gli stessi valori scalati in riferimento ai risultati per un canale piano parallelo
a pari Re con pieno sviluppo idrodinamico e termico: f0 = 24/ Re e Nu0 = 7.5407
[12]. Questi quattro canali sono presentati, figure 4.33, 4.34, 4.35, 4.36, visualizzando
la geometria e le prestazioni, il campo fluidodinamico e il campo termico rispettivamente. Dall’analisi della sequenza di canali è evidente che l’aumento di scambio
termico sia favorito dalle zone di accelerazione e decelerazione del flusso sotto forma
di vortici in corrispondenza delle anse. Sono le geometrie a profilo curvo a indurre rimescolamento della corrente fluida e di conseguenza consentire distribuzioni di
temperatura non-uniformi [9]. Semplificando, nel problema trattato si riconosce l’influenza principalmente di due parametri, ovvero la lunghezza del modulo e l’altezza
dell’asperità o corrugazione. L’altezza media del canale, a cui le prime due grandezze
sono relazionate (tabella 2.2), è infatti costante.
84
Ottimizzazione 2D
0.9
0.8
0.7
c
f
0.6
d
0.5
0.4
0.3
a
0.2
b
0.1
8
9
10
Nu
11
12
(a) Fronte di Pareto 2D.
7
6
c
f/f0
5
d
4
3
a
2
b
1
1
1.1
1.2
1.3
1.4
Nu/Nu0
1.5
1.6
1.7
(b) Fronte di Pareto 2D scalato sul canale piano.
Figura 4.32 Fronte di Pareto, risultato dell’ottimizzazione 2D: cerchiati i quattro individui
selezionati.
85
4.5 Analisi dei Risultati
(a) Geometria e prestazioni.
(b) Campo di velocità: in scala il modulo di velocità.
(c) Campo termico: in scala θ .
Figura 4.33 Canale a: f = 0.151 Nu = 7.864.
86
Ottimizzazione 2D
(a) Geometria e prestazioni.
(b) Campo di velocità: in scala il modulo di velocità.
(c) Campo termico: in scala θ .
Figura 4.34 Canale b: f = 0.285 Nu = 8.721.
87
4.5 Analisi dei Risultati
(a) Geometria e prestazioni.
(b) Campo di velocità: in scala il modulo di velocità.
(c) Campo termico: in scala θ .
Figura 4.35 Canale c: f = 0.577 Nu = 10.921.
88
Ottimizzazione 2D
(a) Geometria e prestazioni.
(b) Campo di velocità:
velocità.
in scala il modulo di
(c) Campo termico: in scala θ .
Figura 4.36 Canale d: f = 0.691 Nu = 11.614.
4.6 Confronto geometrie smooth e spigolose
4.6
89
Confronto geometrie smooth e spigolose
In un secondo momento, l’elemento di novità di questo lavoro, ovvero l’utilizzo di
curve NURBS nella definizione della geometria, è stato messo a confronto con canali realizzati mediante spezzate (linear-piecewise), quindi in assenza di smoothness.
Questi ultimi sono sicuramente più vicini alla pratica industriale.
Rispetto alle geometrie finora trattate, questi nuovi canali sono caratterizzati da un
profilo di più semplice definizione con un minor numero di variabili. In figura 4.37
sono riportati i nuovi gradi di libertà, a cui va aggiunta, coerentemente con quanto
fatto per il profilo NURBS, l’entità della traslazione del profilo superiore.
Figura 4.37 Esempio di costruzione del profilo a spezzate .
Le variabili sono quindi riportate in tabella 4.16.
Variabile
L
h
φin
φout
trasl
Descrizione
lunghezza del modulo
altezza della corrugazione
inclinazione spigolo ingresso
inclinazione spigolo uscita
traslazione del profilo superiore
Tabella 4.16 Variabili per il canale a spigolo.
In figura 4.38 si può osservare un esempio di siffatto canale. I campi di variazione delle variabili libere consentono di ottenere geometrie coerenti, per dimensioni,
a quelle con profilo smooth. In figura 4.39 sono rappresentati i due fronti di Pareto,
rispettivamente dei canali con profilo a spigolo (∗) e con profilo NURBS (). Del
processo di ottimizzazione di questa più rudimentale tipologia di canali, si accenna
soltanto al fatto che il minor numero di variabili in gioco ha consentito di ottenere
una più rapida convergenza rispetto al caso precedentemente trattato. Si è unicamente
90
Ottimizzazione 2D
Figura 4.38 Esempio di canale con profilo a spezzate.
eseguita una ottimizzazione MOGA-II di tipo multi-obiettivo. Il fronte di Pareto del
canale a spigolo si può con buona ragione considerare come limite per questo tipo di
geometria: in più ottimizzazioni successive non si è notato alcun sensibile miglioramento delle prestazioni. Come lecito attendersi, i canali con profilo NURBS, nono-
Figura 4.39 Confronto tra individui dominanti delle due differenti tipologie di canali per Re =
200.
stante la maggior complessità geometrica e quindi costruttiva, consentono di ottenere
un pari scambio termico con un coefficiente d’attrito che in media è del 25% inferiore
a quello del profilo spigoloso. Ovvero si riesce ad ottenere un pari mescolamento del
4.6 Confronto geometrie smooth e spigolose
91
flusso, in modo più dolce ed accompagnato. In figura 4.40 sono rappresentati il campo
fluidodinamico e termico allo stesso modo che per il profilo NURBS. La geometria
in questione è quella di figura 4.38: il coefficiente d’attrito è f = 0.296, il numero di
Nusselt è Nu = 8.236.
(a) Campo di velocità: in scala il modulo di
velocità.
(b) Campo termico: in scala θ .
Figura 4.40 Esempio di soluzione per un canale con profilo spigoloso: f = 0.296 Nu = 8.236.
5
Definizione geometrica e
termofluidodinamica del caso
3D
Fin’ora ci si è occupati di canali infinitamente lunghi in direzione z, fatto che consente
di studiarne il campo termofluidodinamico in forma bidimensionale. In questo modo
tuttavia si viene a perdere la possibilità di incrementare lo scambio termico attraverso
il rimescolamento ottenuto con vorticosità aventi come asse la direzione principale di
moto. La nuova componente non-nulla di velocità promuove lo scambio termico con
perturbazioni periodiche delle zone di ricircolazione in corrispondenza delle anse del
canale. L’estensione al caso tridimensionale, con l’introduzione della terza componente di moto, dei problemi visti, vuole studiare se sia possibile avere un aumento del
numero di Nusselt, senza una eccessiva penalizzazione del coefficiente d’attrito, con
l’instaurarsi di vorticosità sul piano ortogonale alla direzione di moto [9, 20, 23].
La complessità del problema aumenta non tanto dal punto di vista metodologico,
è infatti sufficiente l’addizione di due piani di simmetria, quanto più dal punto di vista
computazionale. Crescono infatti il numero parametri di definizione della geometria,
oltre che, patologicamente, il numero di gradi di libertà che il solver deve affrontare
nella risoluzione del problema termofluidodinamico.
Per questo motivo si sono tentate, in processo di progressiva semplificazione del
modello, diverse modalità di costruzione geometrica del canale tridimensionale. Concettualmente sono tutte valide, ma non tutte scendono a patti con la realtà.
Questa parte del lavoro, presentata nel corrente capitolo e nel successivo, vuole
essere una fase preliminare di risoluzione ed ottimizzazione di canali tridimensionali.
In realtà non si dovrebbe nemmeno parlare di ottimizzazione. Si vuole qui anticipare
che per le non sufficienti risorse temporali e computazionali, si è preferito eseguire,
come parte rilevante di analisi, una user sequence (ovvero una sequenza basata su di
un database creato dall’utente [22]) di alcune geometrie. In questo modo si è voluta
verificare l’applicabilità del metodo e l’eventuale tendenza di miglioramento nel caso tridimensionale. Per analoghe problematiche, si è inoltre stati costretti a ridurre il
93
94
Definizione geometrica e termofluidodinamica del caso 3D
numero di Reynolds da 200 a 100, per garantire un adeguato valore di convergenza
numerica. Un confronto fra i risultati bidimensionali e tridimensionali è stato approntato, ricalcolando gli individui del fronte di Pareto ottenuto nel caso 2D a Re = 100.
Soltanto in un secondo momento si è fatto un tentativo di ottimizzazione.
Le diverse tipologie geometriche saranno presentate in ordine crescente di complessità. In seguito si farà poi cenno a cosa si è dovuto aggiornare nella procedura di
calcolo per poter applicare la metodologia descritta in sezione 2.2 al caso tridimensionale. I risultati di questa analisi sono riportati nel capitolo successivo 6, coerentemente
a quanto fatto per il caso bidimensionale.
5.1
Canale 3D estruso
Questa tipologia di canale è ottenuta per estrusione lungo la direzione z delle due differenti geometrie trattate nella fase 2D, si vedano figure 4.38 e 2.2. In entrambi i casi
le variabili aggiunte sono la profondità di estrusione, Lz , e la direzione di estrusione
rispetto l’asse z, ψ. La modalità di costruzione è intuibile nelle figure 5.1 e 5.2. Nel
Figura 5.1 Esempio di canale con profilo linear-piecewise.
caso di estrusione da canale con profilo linear-piecewise, figura 5.1, i gradi di libertà
sono 7, si veda tabella 5.1. Ovvero le stesse del caso bidimensionale, si vedano figura
4.37 e tabella 4.16, più due.
5.2 Canale 3D realizzato con superfici NURBS
95
Figura 5.2 Esempio di canale con profilo NURBS.
Variabile
L
h
φin
φout
trasl
Lz
ψ
Descrizione
lunghezza del modulo
altezza della corrugazione
inclinazione spigolo ingresso
inclinazione spigolo uscita
traslazione del profilo superiore
profondità lungo z del modulo elementare
angolo di estrusione
Tabella 5.1 Variabili per il canale 3D con profilo linear-piecewise.
Mentre nel caso di estrusione da canale con profilo NURBS, figura 5.2, le variabili
salgono a 13, ovvero quante nel caso bidimensionale più 2, Lz e ψ. Per convenienza
i gradi di libertà per la NURBS 2D sono riportati nella prima parte successiva tabella
5.2.
5.2
Canale 3D realizzato con superfici NURBS
Inizialmente si intendeva realizzare il canale periodico tridimensionale, modellandone i profili con due superfici NURBS due volte periodiche. Ovvero lungo la direzione principale del flusso si mantiene la periodicità già presente nel caso puramen-
96
Definizione geometrica e termofluidodinamica del caso 3D
te bidimensionale, mentre si introduce un ulteriore effetto periodico lungo la terza
dimensione z.
La costruzione del profilo mediante la superficie NURBS ha come origine ancora
il profilo bidimensionale a 9 punti di controllo ed è schematizzata come segue e in
figura 5.3:
1. il profilo bidimensionale di periodo Lx , si veda sezione 2.1, viene caricato nel
piano xy attraverso i suoi 25 punti di controllo;
2. si procede con l’estensione lungo la terza dimensione, z, dei punti di controllo
che descrivono il canale 2D. I nuovi punti di controllo vengono determinati
per semplice interpolazione lineare lungo la direzione inclinata dell’angolo ψ
rispetto all’asse z. La profondità di interpolazione è pari a Lz . Si individuano
cosı̀ quattro file di punti di controllo, equamente spaziate, la cui proiezione sul
piano xz appare in figura 5.3;
3. la periodicità viene riproposta, oltre che in direzione x, anche in direzione z mediante due operazioni di specchio sui punti di controllo (simmetria pari rispetto
a piani di normale z) rispettivamente a sinistra e a destra(1) del modulo centrale;
4. l’insieme dei punti di controllo costituisce la patch, 25 x 10 punti, che definisce
la NURBS 2 volte periodica.
Per ulteriori dettagli sulle superfici periodiche si rimanda alla sezione 1.3. Il singolo
Figura 5.3 Esempio di patch per superfici NURBS due volte periodiche. In evidenza il singolo
modulo periodico, ovvero la patch iniziale.
modulo è riportato per 3 periodi lungo la direzione x e lungo la direzione z.
(1) Rispetto
la direzione z.
5.2 Canale 3D realizzato con superfici NURBS
97
Una volta realizzata la generica superficie inferiore, la costruzione del canale ricalca quanto fatto per il caso bidimensionale. Il procedimento di taglio mediante due
rette del caso 2D, viene qui effettuato mediante quattro piani come nelle figure 5.4 e
5.5, dove lo stesso canale è visualizzato su due differenti angolazioni. Si procede con:
1. generazione di una superficie NURBS doppiamente periodica per il profilo inferiore di lunghezza pari a 3 volte il periodo Lx (lungo la direzione x) e larghezza
pari a 3 volte il periodo Lz (lungo la direzione z);
2. trasformazione della forma NURBS in più patch di Bézier equivalenti; questo
passaggio è necessario poiché FEMLAB tratta superfici in forma di Bézier del
III o del IV ordine;
3. traslazione dei punti di controllo della superficie inferiore lungo l’asse x di
una quantità traslx ∈] − Lx , +Lx [, dove Lx è il periodo della superficie lungo
la direzione x;
4. traslazione dei punti di controllo della nuova superficie lungo l’asse z di una
quantità traslz ∈] − Lz , +Lz [, dove Lz è il periodo della della superficie lungo la
direzione z;
5. traslazione dei punti di controllo appena trovati lungo l’asse y di una quantità
pari a Hav : si è ottenuto il profilo superiore con un’altezza media del canale pari
a Hav ;
6. il canale cosı̀ ottenuto viene tagliato dalle sezioni di ingresso e uscita e da due
ulteriori piani di normale z per le due superfici laterali di simmetria.
98
Definizione geometrica e termofluidodinamica del caso 3D
Figura 5.4 Esempio di costruzione di canale 3D con superfici NURBS.
Figura 5.5 Altra angolazione del canale.
Le variabili geometriche che definiscono questa parametrizzazione sono riportate in
tabella 5.2. Le prime 11 sono le medesime del caso bidimensionale, si veda anche
figura 2.4 per individuare i punti. Si introducono 3 ulteriori gradi di libertà per un
totale di 14. Oltre all’accresciuto numero di variabili, si sono riscontrare numerose
5.3 Definizione del problema termofluidodinamico
Parametro
Lx
Lwave
Descrizione
lunghezza lungo x del modulo elementare
lunghezza lungo x del profilo curvilineo
∆x3
distanza lungo x del punto 3 rispetto al punto 2
ρ4
distanza polare del punto 3 rispetto al punto 2
ϑ4
anomalia polare del punto 3 rispetto al punto 2
∆x5
distanza lungo x del punto 5 rispetto al punto 2
∆y5
distanza lungo y del punto 5 rispetto al punto 2
ρ6
distanza polare del punto 6 rispetto al punto 5
ϑ6
anomalia polare del punto 6 rispetto al punto 5
∆7
distanza lungo x del punto 8 rispetto al punto 7
traslx
99
entità di traslazione lungo x dei punti
di controllo del profilo superiore
Lz
ψ
profondità lungo z del modulo elementare
angolo rispetto la direzione positiva dell’asse z
per definire la direzione di interpolazione
traslz
entità di traslazione lungo z dei punti
di controllo del profilo superiore
Tabella 5.2 Variabili per parametrizzazione 3D con superfici NURBS. Nella prima parte le
variabili del caso 2D.
difficoltà nel caricare l’elevato numero di patch di Bézier nel CAD di FEMLAB,
nonché nelle operazioni di taglio, in cui errori puramente numerici erano sempre in
agguato. Anche per questi motivi l’utilizzo di queste geometrie è stato abbandonato,
con l’auspicio di poterle riproporre in futuro con un software geometrico adeguato e
con elaboratori più potenti. Purtroppo FEMLAB non si è rivelato lo strumento adatto
a trattare geometrie tridimensionali complesse.
5.3
Definizione del problema termofluidodinamico
L’estensione al caso tridimensionale delle equazioni non-dimensionali caratteristiche
del problema di scambio termico in regime di convezione forzata non pone particolari
difficoltà. La loro forma rimane sostanzialmente la stessa, si tratta infatti di ridefinire
100
Definizione geometrica e termofluidodinamica del caso 3D
il vettore velocità aggiungendo la terza componente w:
u = (u, v, w)
Brevemente per il campo di moto si ha in forma vettoriale
1 2
∇ u − ∇P + F
Re
(5.1a)
∇·u = 0
(5.1b)
u · ∇u =
dove F = (β , 0, 0). Per il termine β si veda il paragrafo successivo e in generale la
sezione 2.2. Lo stesso dicasi per gli altri simboli: nulla cambia rispetto alle definizioni
già date. Per il campo termico è
u · ∇θ =
1 2
∇ θ
Pe
(5.2)
Inoltre tutti gli integrali che nel caso bidimensionale sono definiti su una linea, diventano ora integrali di superficie, e analogamente quelli di superficie vengono definiti
sul volume.
5.3.1
Condizioni al contorno
Si faccia riferimento a figura 5.6.
Pareti Sulle pareti wall sono imposte rispettivamente una condizione di no-slip per il
campo fluidodinamico e di temperatura assegnata per il campo termico:
u=0 θ =0
Inlet e Outlet In ingresso e uscita del canale si prescrivono condizioni di periodicità
su velocità e pressione per il campo fluidodinamico e temperatura per quello
termico allo stesso modo che nel caso bidimensionale, si vedano a proposito le
sezioni 2.2.3 e 2.2.9. Per evitare la soluzione banale del campo termico, i.e.
θ = θw = 0 ∀(x, y, z), si è imposto all’ingresso che
R
θ u · n dS
= θb,in = 1
in u · n dS
Rin
(5.3)
Pareti laterali L’elemento di novità rispetto al caso 2D riguarda la presenza delle
due pareti laterali, sulle quali per il campo di moto si impongono condizioni di
simmetria, i.e. symmetry. Tale condizione si traduce in componenti di velocità
normali alla superficie nulle:
u·n = 0
5.3 Definizione del problema termofluidodinamico
101
dove n è il versore normale alla superficie, e assenza di componenti tangenziali
delle forze viscose:
τ ·t = 0
dove τ rappresenta il generico sforzo viscoso e t il generico versore tangente
alla superficie.
Figura 5.6 Condizioni al contorno per il canale 3D.
La condizione di flusso nullo è assunta per la temperatura adimensionale:
∇θ · n = 0
5.3.2
Cenni alla procedura risolutiva
Analogamente al caso bidimensionale, per affrontare la soluzione in forma adimensionale si ricerca quel gradiente pressorio adimensionale β che dia una velocità media
Uav unitaria. Cosı̀ facendo il numero di Reynolds è:
Re =
La definizione di Uav è ora
R
1
ν
u dV
Ω dV
Uav = RΩ
dove Ω è il dominio tridimensionale e dV = dx dy dz è l’infinitesimo di volume.
(5.4)
102
Definizione geometrica e termofluidodinamica del caso 3D
Nulla cambia nella determinazione del coefficiente d’attrito. Il diametro idraulico
del canale è ancora Dh = 2Hav = 1, per cui
β
2
Per la determinazione delle caratteristiche di scambio termico l’equazione 2.42,
subisce una lieve modifica. Il coefficiente σ , ovvero il rapporto delle temperature di
bulk tra inlet e outlet, si ricava in forma iterativa come
Z
2 1
∂ θ 1 (5.5)
= 1+
dS
σ Lz Pe w ∂ n f=
dove dS è l’infinitesimo di superficie. Una volta ottenuto per σ un valore che garantisca il bilancio fra il salto entalpico fra ingresso e uscita del canale e il flusso termico a
parete, si ha:
Hav
Dh
Nu = ln σ
Pe = ln σ
Pe
2L
4L
Per quanto ancora sia stato sottinteso, si faccia sempre riferimento al caso bidimensionale.
5.4
Meshing & Solving
FEMLAB consente di generare per domini tridimensionali griglie non strutturate tetraedriche. Analogamente ai domini bidimensionali è possibile controllare un certo
numero di parametri, i.e. la dimensione massima di un elemento, o il rapporto fra raggio di curvatura e la dimensione del tetraedro su di un bordo. In questa fase la scelta
di tali parametri si è sbilanciata sul mantenimento di un ridotto numero di elementi
per ottenere tempi di calcolo sostenibili. Questo a scapito purtroppo dell’accuratezza
della soluzione. Per avere costi computazionali accettabili il numero di elementi non
deve superare l’intorno di 10000, sicuramente negando una opportuna valutazione dei
gradienti a parete. Ad ogni modo un elevato numero di geometrie è risultato non processabile a causa di mesh non congruenti. Tale disfunzione si è verificata soprattutto
in presenza di un rilevante numero di bordi dalle dimensioni non omogenee a seguito
della trasformazione NURBS-Bézier. Questo testimonia i limiti del meshatore interno
a FEMLAB. Si peraltro riscontrato che l’implementazione delle condizioni di periodicità è piuttosto critica nei problemi tridimensionali. Infatti a parità di geometria e
gradi di libertà, la presenza di tali condizioni alle sezioni di ingresso e uscita aumenta
considerevolmente i tempi di calcolo del problema fluidodinamico. I solver utilizzati
sono e per il campo di moto e per il campo termico di tipo non-lineare e diretto. Il
time-consuming per ciascun canale è in media di 2 ore di CPU su processore AMD 64
3.6GHz.
Come test-case dello script realizzato si sono valutati il coefficiente d’attrito f0 e
il numero di Nusselt Nu0 per il canale piano parallelo. Si ritiene di puntualizzare che,
mentre nulla cambi sul valore di Nu0 , per Re = 100, f0 = 0.24.
6
Risultati dell’analisi
tridimensionale
Come già anticipato nel precedente capitolo, in questa sezione sono riportati i risultati
dell’analisi di ottimizzazione effettuata su canali tridimensionali in regime di convezione forzata. La tipologia geometrica trattata è quella di canali ottenuti per estrusione
di domini bidimensionali. In questa ultima parte del lavoro si possono riconoscere due
differenti strategie di analisi. Con la prima, tenendo in debito conto le limitate risorse
di calcolo, si è cercato di verificare l’applicabilità del metodo risolutivo al caso tridimensionale. modeFRONTIER permette infatti all’utente di definire un database di
variabili, ovvero un certo numero di individui che successivamente vengono processati. Si definisce una cosiddetta user sequence, cioè un’analisi parametrica [22]. In
questo modo, selezionando ad hoc i canali, si è sondata la possibilità di incremento
di scambio termico a parità di coefficiente d’attrito nel passaggio da un canale 2D
ad uno 3D. Con la seconda, pensata come suggerimento a sviluppi futuri, si è tentata una ottimizzazione mediante algoritmo genetico MOGA-II su geometrie realizzate
per estrusione di canali con profilo linear-piecewise, si veda 5.1. Per agevolare la
procedura di calcolo il numero di Reynolds è stato portato da 200 a 100.
6.1
Prima strategia di analisi
Per valutare i miglioramenti nelle prestazioni nel caso tridimensionale, le geometrie
ottenute sono state messe a confronto con quelle appartenenti ai fronti di Pareto bidimensionali da cui hanno origine. Tale confronto è stato portato avanti parallelamente
per profili linear-piecewise e profili NURBS. Per poter effettuare questo accostamento
i fronti di Pareto per le geometrie linear-piecewise 2D e NURBS 2D sono stati ricalcolati a Re = 100. Per ciascuno dei due casi il data-base è stato creato integrando le
design variables bidimensionali con i parametri Lz e ψ. Per semplicità la profondità
di estrusione Lz è stata mantenuta costante al valore di Hav , mentre si sono studiati
4 angoli di estrusione ψ differenti, come riportato in tabella 6.1. Si hanno dunque
quattro prove per ciascun canale. La prima, con ψ = 0 ◦ , saggia la bontà del calcolo,
103
104
Risultati dell’analisi tridimensionale
Variabile
Lz
Descrizione
profondità di estrusione
ψ
angolo di estrusione
Valore
0.5
0◦
20 ◦
30 ◦
40 ◦
Tabella 6.1 Definizione dei parametri per il canale 3D estruso.
infatti si dovrebbe ottenere il medesimo risultato del caso bidimensionale. Per facilitare la convergenza della procedura di calcolo 3D, i valori di tentativo per il gradiente
pressorio, β , e per il rapporto fra le temperature adimensionali di bulk fra ingresso e
uscita del canale, σ , sono quelli corretti già ottenuti nelle prove 2D. Per convenienza
sia per il caso di profilo linear-piecewise, sia per il caso NURBS si è scelto un numero ragionevolmente limitato di design da estrudere, tale però da rappresentare diverse
aree dello spazio delle soluzioni.
Nelle figure 6.1 e 6.2 sono riportati i risultati del confronto. Si ricorda ancora che,
mentre il fronte di Pareto bidimensionale è ripresentato per intero, dei canali tridimensionali si è fatta una selezione. Va poi detto che molte geometrie, soprattutto del più
complesso tipo NURBS, figura 6.2, non sono risultate processabili a causa di mesh
incongruenti. Il confronto è eseguito sulle funzioni obiettivo Nu/Nu0 e f / f0 , dove
24
22
canali 2D
canale 3D (a)
canale 3D (b)
canale 3D (c)
canale 3D (d)
20
18
16
f/f0
14
12
10
8
6
4
1.2
1.4
1.6
Nu/Nu0
1.8
2
2.2
Figura 6.1 Confronto 2D-3D per geometrie con profilo linear-piecewise.
Nu0 = 7.5407 e f0 = 24/ Re. Dalle figure è evidente come, per entrambe le tipologie
105
6.1 Prima strategia di analisi
24
22
canali 2D
canale 3D (a)
canale 3D (b)
canale 3D (c)
canale 3D (d)
20
18
f/f0
16
14
12
10
8
6
4
1.2
1.4
1.6
Nu/Nu0
1.8
2
2.2
Figura 6.2 Confronto 2D-3D per geometrie con profilo NURBS.
di canale 3D estruso, i.e. linear-piecewise and smooth, il numero di Nusselt tenda
a crescere senza essere accompagnato da un sensibile aumento del fattore d’attrito.
Questo effetto si verifica quando il campo di moto è tipicamente tridimensionale, ovvero non riconducibile alla forma 2D. In tale situazione la presenza di moti secondari
[9], che spezzano le zone di ricircolazione bidimensionali, favorisce il trasferimento di massa e di calore (mixing). Si potrebbe obiettare che l’incremento del numero
di Nusselt sia principalmente dovuto a una maggior superficie esposta allo scambio
termico. Da un’analisi attenta si può invece notare che l’aumento di superficie alle
pareti, a seguito di una estrusione con ψ , 0, ha una sensibilità sull’incremento di Nu
inferiore rispetto all’induzione dei moti secondari, ottenuta con la deformazione della
parete stessa. A conferma di quanto detto, si riportano in tabella 6.2 l’incremento di
area delle pareti (wall) e l’incremento del numero di Nusselt per i diversi angoli di
estrusione ψ del canale (a) di figura 6.1. Con A2D = 1.863 e Nu2D = 9.832 si vogliono indicare i valori rispettivamente di area e di Nusselt per il canale estruso con
angolo ψ = 0 ◦ , caratterizzato cioè da un campo termofluidodinamico bidimensionale.
Ancora da tabella 6.2 si può ragionevolmente evincere come al crescere dell’angolo di
estrusione, migliori lo scambio termico. Infatti quanto più ci si allontana dalla configurazione bidimensionale, tanto più considerevole è la presenza dei vortici secondari
tridimensionali.
D’altro canto il citato aumento di superficie per estrusione con ψ , 0, rispetto
ψ = 0, si manifesta soltanto per quei tratti di parete la cui normale abbia componente
z non nulla, cioè proiezione sul piano xy non degenerante su una linea. Si veda a tal
proposito l’esempio riportato in figura 6.3, dove uno stesso profilo linear-piecewise è
106
Risultati dell’analisi tridimensionale
A/A2D
Nu/ Nu2D
ψ
1.000
1.000
0◦
1.015
1.013
20 ◦
1.037
1.060
30 ◦
1.076
1.160
40 ◦
Tabella 6.2 Incremento del numero di Nusselt in funzione dell’aumento di superficie, wall,
relativi per diversi angoli di estrusione.
stato estruso con due differenti valori di angolo ψ. In figura 6.3(a) è riportata una vista
(a) Canale con ψ = 0 ◦
(b) Canale con ψ = 30 ◦
Figura 6.3 Esempio di estrusione con aumento di superficie per la parete di un canale.
della superficie ottenuta con angolo di estrusione nullo. La freccia indica la proiezione
di tale superficie su di un piano parallelo a yz. Facendo riferimento al caso con angolo
107
6.2 Seconda strategia di analisi
di estrusione ψ = 30 ◦ , figura 6.3(b), si può osservare come tale superficie abbia, oltre
che proiezione sul piano yz pari alla precedente, proiezione non nulla sul piano xy.
La proiezione sul piano xz è la medesima nei due casi in termini di area. In questo
semplice esempio l’aumento di superficie, indicato dalla freccia nella zona ∆, è pari a:
1
−1
Lz h
cos ψ
Lo stesso dicasi per geometrie più complesse con profili smooth NURBS. É possibile
infatti approssimare l’andamento della parete curvilinea mediante un’opportuna serie
di profili a spigolo.
In tabella 6.3 sono esposti i risultati per un canale NURBS estruso a diversi angoli
ψ. Ancora una volta è manifesto per lo scambio termico il positivo effetto combinato
di moti secondari e dei vortici longitudinali stazionari: il numero di Nusselt cresce
all’aumentare dell’angolo di estrusione. A scopo puramente illustrativo, la forma di
tali vortici secondari è rappresentata in figura 6.4 per il medesimo canale con due
differenti estrusioni, ψ = 20 ◦ figura 6.4(a) e ψ = 40 ◦ figura 6.4(b).
f/f0
Nu/Nu0
ψ
3.807
1.307
0◦
3.932
1.366
20 ◦
3.937
1.508
30 ◦
3.909
1.530
40 ◦
Tabella 6.3 Variazione del fattore d’attrito e del numero di Nusselt in funzione dell’angolo di
estrusione ψ per un canale smooth.
Si vuole infine rivolgere l’attenzione sul fattore d’attrito f . Come già evidente
dalla sezione 2.2.4, vale che:
f ∝ (τx , S)
dove τx è lo sforzo a parete medio in direzione x, mentre S è la superficie bagnata.
É evidente che un aumento della superficie giochi a sfavore delle perdite di carico,
tuttavia la tridimensionalità e l’adozione di profili smooth possono garantire un flusso
maggiormente guidato con la conseguente riduzione di τx .
6.2
Seconda strategia di analisi
Data la natura pionieristica di quest’ultima parte del lavoro, si è cercato di mantenere il più possibile ridotto l’onere computazionale. La tipologia geometrica su cui si
è tentata l’ottimizzazione è la linear-piecewise, che rispetto la NURBS presenta un
108
Risultati dell’analisi tridimensionale
(a) ψ = 20 ◦
(b) ψ = 40 ◦
Figura 6.4 Esempio di flussi secondari 3D.
numero di variabili libere decisamente inferiore. Delle 7 individuate in tabella 5.1, si
blocca la profondità di estrusione Lz al valore 0.5. Ricordando quanto espresso dalle
4.1 e 4.2 nel caso di ottimizzazione del canale bidimensionale, si ha ora che
[ f , Nu] = g (L, h, φin , φout ,trasl, ψ)
(6.1)
L’algoritmo scelto in questa fase non può che essere MOGA-II viste le sue caratteristiche di robustezza. Il numero di gradi di libertà è pari a 6, ragion per cui il numero
di individui costituenti la popolazione iniziale, D.O.E., è fissato a 25. Il numero di
generazioni è pari a 10, per un totale di 250 individui processati. Non si è rinunciato
a porre dei vincoli e sul numero di Nusselt e sul fattore d’attrito, si veda tabella 6.4.
Siano Nu0 = 7.5407 il numero di Nusselt e f0 = 0.24 il coefficiente d’attrito per il
canale piano per Re = 100 [12]. La generazione del D.O.E. è al solito affidata alVincolo
Nu/ Nu0
f / f0
Relazione
>
<
Limite
1.1
10
Tabella 6.4 Constraints per il lancio MOGA-II.
l’algoritmo Sobol. I campi di variabilità dei parametri liberi sono riportati in tabella
6.5. Nonostante l’amplificazione notevole dei tempi di calcolo, i risultati estratti sono
incoraggianti. Di seguito si riportano due tipologie di confronti. La prima, figura 6.5,
109
6.2 Seconda strategia di analisi
Parametro
Campo
Base
L
[0.800; 2.000]
501
h
[0.050 L; 0.400 L]
501
φin
[0 ◦ ; 60 ◦ ]
501
φout
[0 ◦ ; 60 ◦ ]
501
trasl
[−0.499 L; 0.499 L]
501
ψ
[5 ◦ ; 50 ◦ ]
51
Tabella 6.5 Dominio di definizione delle variabili.
rappresenta il fronte di Pareto bidimensionale a confronto con quello tridimensionale sempre per profili linear-piecewise. Anche se il fronte 3D è molto più rado (solo
26 design), per il limitato numero di individui processati, è evidente il miglioramento
delle prestazioni rispetto al caso bidimensionale sempre per profili linear-piecewise.
Si fa notare che i risultati 2D sono intesi per Re = 100, al pari del 3D. Anche quei
20
2D linear−piecewise
3D linear−piecewise
18
16
14
f/f0
12
10
8
6
4
2
0
1
1.2
1.4
1.6
Nu/Nu0
1.8
2
2.2
Figura 6.5 Fronti di Pareto 2D-3D per geometrie con profilo linear-piecewise.
canali 3D, cui corrispondono valori minimi del numero di Nusselt, zona in basso a sinistra di figura 6.5, sono comunque caratterizzati da un coefficiente d’attrito inferiore
al caso 2D, fatto che testimonia ulteriormente come i flussi secondari promuovano lo
scambio termico. In figura 6.6 si riporta la distribuzione dei valori della variabile ψ. A
110
Risultati dell’analisi tridimensionale
0.31
0.27
frequenza
0.23
0.19
0.15
0.12
0.08
0.04
0
32
34
36
38
ψ
40
42
44
46
Figura 6.6 Istogramma di frequenza per la variabile ψ.
questo primo livello di analisi sembra che l’ottimizzatore concentri i migliori risultati
in termini di scambio termico su un valore medio dell’angolo di estrusione pari a 39 ◦ .
La seconda tipologia, figura 6.7, introduce nel grafico di figura 6.5 anche il fronte di
Pareto bidimensionale per canali smooth. Anche qui il progresso non manca, tutta25
2D linear−piecewise
2D NURBS
3D linear−piecewise
20
f/f0
15
10
5
0
1
1.2
1.4
1.6
1.8
Nu/Nu0
2
2.2
2.4
Figura 6.7 Fronti di Pareto a confronto: profilo linear-piecewise 2D-3D e profilo smooth 2D.
6.2 Seconda strategia di analisi
111
via è meno palese, com’era lecito attendersi. I profili NURBS infatti sono in genere
caratterizzati da perdite di carico contenute rispetto a profili a spigolo. Il grafico di figura 6.7 lascia intravedere future prospettive di sviluppo con l’ottimizzazione di profili
NURBS-derivati. A livello attuale tuttavia i costi in termini di CPU time sono proibitivi. Si dovrebbe infatti passare dai 6 gradi di libertà attuali (7 se si lasciasse libera
anche la profondità di estrusione) ai 13 in caso di estrusione da profili NURBS 2D.
Conclusioni
In questo lavoro di tesi si è descritta una procedura di ottimizzazione di forma di tipo
multi-obiettivo, applicata a canali corrugati periodici, nei quali vi è scambio termico
convettivo fra le loro pareti ed il fluido che li attraversa. Questi sistemi rappresentano
l’elemento ripetitivo fondamentale di numerosi scambiatori di calore, in particolare
quelli liquido-gas e gas-gas. Il fine ultimo è quello di migliorarne le prestazioni di
scambio termico e fluidodinamiche.
Come emerso in lavori simili, la parametrizzazione della geometria riveste un ruolo fondamentale nel processo di ottimizzazione. Nel capitolo 1 si è dato spazio alla
documentazione delle forme NURBS (Non-Uniform Rational B-Spline) utilizzate per
la descrizione del generico canale in due e in tre dimensioni. Il loro utilizzo ha consentito di definire agevolmente il concetto di continuità e smoothness del profilo dei
canali e di tracciare forme complesse semplicemente giocando su ordine, numero di
punti di controllo e nodi (knots). Ci si è inoltre soffermati sulla problematica della trasformazione di un’entità NURBS in più forme semplici di Bézier e sono state
introdotte le NURBS in forma periodica.
Nel capitolo 2 si è sviluppata in modo generale la metodologia di risoluzione per
il campo termofluidodinamico nel problema di convezione forzata, tipico degli scambiatori di calore. Allo scopo di ridurre le richieste computazionali si è preso atto della
periodicità spaziale del dispositivo e conseguentemente si è potuto studiarne il campo
termofluidodinamico di un unico modulo ripetitivo. Il flusso preso in considerazione
è caratterizzato da un numero di Prandtl pari a 0.7 ed è stato assunto in condizioni
laminari, stazionarie e in regime di pieno sviluppo fluidodinamico e termico. Rispetto
ad un precedente lavoro di tesi, di cui questo vorrebbe essere l’ideale continuazione, si è individuata una tecnica di risoluzione più efficiente in forma completamente
adimensionale.
La fase di implementazione in FEMLAB è stata descritta nel capitolo 3, dove è
stata riportata la traduzione in forma CAD/FEMLAB dei due vincoli costruttivi imposti durante la successiva fase di ottimizzazione, ovvero la stampabilità e la curvatura
non eccessiva dei profili.
Quanto realizzato finora è divenuto al capitolo 4 il motore risolutivo, legato al software modeFRONTIER, per l’ottimizzazione bidimensionale, descritta in tutta la sua
113
114
Conclusioni
fase evolutiva. Le variabili geometriche del canale costituivano i gradi di libertà del
problema, input. Gli obiettivi, output, sono stati individuati nel numero di Nusselt,
indicatore di scambio termico, e nel coefficiente d’attrito, indicatore delle perdite di
carico all’interno dello scambiatore di calore. La notevole libertà nelle forme concessa
dall’uso di NURBS, ha fatto emergere i limiti di FEMLAB nella procedura di discretizzazione. La trasformazione delle curve NURBS in più curve di Bézier, richiesta
da FEMLAB, ha causato il moltiplicarsi del numero di bordi del canale, spesso dalle dimensioni tra loro non omogenee. Questo aspetto ha messo in crisi il meshatore,
causando la creazione di griglie non coerenti. Nella successiva analisi tridimensionale questo problema si è rivelato in forma esasperata, poiché il numero di individui
con raggi di curvatura molto piccoli, che di conseguenza si è dovuto scartare, si è
dimostrato crescere notevolmente.
Come accennato, si è quindi passati a descrivere il medesimo problema in forma
tridimensionale, attendendosi un miglioramento delle prestazioni grazie alla presenza
dei moti secondari. Tali moti sono costituiti da vortici che giacciono su piani normali
alla direzione principale del moto. Nel capitolo 5 si è voluto segnalare le diverse tappe
che hanno portato alla scelta di quali tipologie geometriche descrivere. La richiesta
computazionale per la risoluzione numerica di questi problemi tridimensionali di convezione forzata si è rivelata notevole, tant’è che si è dovuto dimezzare il numero di
Reynolds della simulazione. In aggiunta il CAD FEMLAB si è dimostrato poco adatto
a trattare canali limitati da superfici NURBS complesse.
Come espresso nel capitolo 6 lo studio tridimensionale è stato suddiviso in due fasi. Nella prima è stata affrontata un’analisi parametrica di alcuni canali, realizzati per
semplice estrusione di geometrie bidimensionali, al variare del solo angolo di estrusione. Si è cosı̀ provata l’estensibilità della metodologia proposta al capitolo 2 al caso
tridimensionale e si è verificato l’incremento di scambio termico auspicato per canali smooth profile e per canali linear-piecewise profile. Si è successivamente tentata
una pura ottimizzazione sulle sole geometrie linear-piecewise profile, ovviamente più
semplici e caratterizzate da un minor numero di variabili progettuali rispetto al caso
NURBS profile. Nonostante i problemi di mesh e l’impossibilità di garantire risultati
indipendenti dalla griglia, si sono ottenuti risultati molto incoraggianti. Il fronte di
Pareto di quest’ultima ottimizzazione tridimensionale è stato posto a confronto con
quelli ottenuti nell’ottimizzazione bidimensionale per il caso NURBS profile e per il
caso linear-piecewise profile, ed ha fornito un’indicazione del vantaggio di utilizzare canali corrugati caratterizzati da campi termofluidodinamici tridimensionali. Per
ragioni di tempo non si è potuto procedere nell’ottimizzazione di canali tridimensionali NURBS profile, anche se la lettura dei risultati lasciava sperare che l’utilizzo di
tali forme avrebbe portato verso un ulteriore miglioramento delle prestazioni, almeno per quel che riguardava la riduzione delle perdite di carico. La parametrizzazione
NURBS ha infatti permesso di ottenere profili dolci, che bene hanno dimostrato accompagnare il flusso. L’ottimismo che nell’ultimo capitolo si è lasciato trasparire non
deve nascondere i limiti di questo studio. Non è stato possibile processare individui
115
con mesh superiori a 11000 elementi circa, per cui i risultati sono stati ottenuti con un
livello di accuratezza minimo, sia pur accettabile per un primo confronto con il caso
bidimensionale.
Ulteriori indagini, magari accompagnate da un’evoluzione di FEMLAB ed affiancate anche dall’utilizzo di altri solutori CFD, sono raccomandate. In entrambe le fasi
di ottimizzazione si è fatto uso del solo solutore diretto, ossia non-iterativo. Si è
sentita la mancanza di un solver efficiente, che non entrasse in conflitto con l’imposizione delle condizioni al contorno di periodicità, nonché della possibilità di un rapido
utilizzo di un solver segregato, sicuramente più indicato nell’analisi tridimensionale.
In tutta la tesi i problemi sono stati affrontati seguendo un indirizzo prevalentemente metodologico. L’imposizione dei vincoli costruttivi si è fermata ad un livello
accademico e ulteriori condizioni sarebbero state necessarie perché le forme più prestanti emerse in questo lavoro potessero passare rapidamente e in modo indolore alla
fase applicativa.
Come ultima indicazione si noti che, nell’ottimizzazione delle prestazioni di sistemi come quelli cui si fa riferimento in questo lavoro, spesso viene tralasciata una
caratteristica fondamentale, ossia la stabilità della performance al variare dei parametri di funzionamento. Nella realtà pratico-industriale devono invece essere sviluppati
prodotti che, entro un certo range, abbiano scarsa sensibilità alla alterazione di tali
parametri. Una analisi di questo tipo è concessa da modeFRONTIER, i.e. Robust
Design, ed è fondata non più su una distribuzione deterministica delle variabili, come
in questo lavoro, bensı̀ su una distribuzione probabilistica con l’introduzione del concetto di media e deviazione standard delle variabili di input e conseguentemente della
soluzione ottenuta. Il modello di una tale indagine è dunque statistico e richiede un
campionamento di dati che appesantisce notevolmente il processo di ottimizzazione,
pur tuttavia coprendo un importante aspetto dal punto di vista applicativo.
Appendice A
I Software Utilizzati.
In questa sezione si darà una breve descrizione dei pacchetti utilizzati nel presente
lavoro di tesi, senza peraltro volere sovrapporsi alla documentazione di cui ciascun
programma è corredato [18], [22] e facendo riferimento a quanto già ampiamente
riportato in precedenti lavori [5], [10].
A.1
modeFRONTIER
modeFRONTIER è per questo lavoro l’organizzatore delle operazioni e del flusso di
dati, avente come scopo ultimo l’ottimizzazione. Ottimizzare significa ottenere la
migliore configurazione di un sistema, caratterizzato da un certo numero di gradi
di libertà, detti anche variabili decisionali, e sottoposto a specifici vincoli, o constraints. Si tratta quindi scegliere la migliore tra un certo numero di possibilità, al
limite infinite! Il processo di ottimizzazione vede quindi la necessità di definire:
- cosa si vuole migliorare, ovvero il modello, semplificazione del sistema che si
intende analizzare;
- quali aspetti o caratteristiche del modello si desidera migliorare;
- con quale metodo e procedimento effettuare la scelta.
Come visto, nei problemi trattati in questa tesi il risultato da raggiungere consta di
due obiettivi, peraltro in contrasto fra loro. Questi due obiettivi(1) sono legati, oltre
che alle caratteristiche del fluido di processo, ad un certo numero di parametri che
vanno a definire la forma geometrica del canale. In generale detto X l’insieme delle
variabili, che il progettista può, più o meno liberamente, sondare e φ l’insieme degli
obiettivi da massimizzare/minimizzare, detti funzioni obiettivo, il suddetto legame si
esprime come:




 φ1 
 X1 
...
...
= g(X) = f
φ=




Xn
φm
(1) Minimizzazione del coefficiente di attrito adimensionale f e massimizzazione dello scambio termico
adimensionale Nu.
117
118
I Software Utilizzati.
Mancando una espressione analitica per g, il processo di ottimizzazione deve basarsi
sull’esperienza, secondo uno schema trial & error. In sostanza è quanto comunemente
si applica nelle scelte quotidiane, circa quali abiti indossare a seconda del clima, o
quale percorso intraprendere a seconda del traffico etc.
modeFRONTIER consente di affrontare queste operazioni di scelta in maniera
sistematica su problemi tecnico-ingegneristici ed economici, attraverso svariati algoritmi, i.e. basati sul gradiente di una funzione, evoluzione genetica, teoria dei giochi. Notevole è la flessibilità che modeFRONTIER dimostra nell’affrontare diverse
tipologie di problemi su altrettanto differenti modelli, non necessitando di algoritmi
dedicati.
A.1.1
Caratteristiche essenziali
In questo lavoro si è utilizzata la versione del software 3.1.0. Attraverso la GUI il
programma permette di:
- descrivere il problema;
- fissare gli obiettivi da raggiungere;
- impostare ed eseguire uno tra i diversi algoritmi di ottimizzazione presenti;
- analizzare i dati ottenuti e far interagire le scelte-macchina con le
scelte-progettista.
modeFrontier consente una ottimizzazione multi-obiettivo e multi-disciplinare, attraverso un agevole accoppiamento con i più svariati programmi di calcolo. Gli algoritmi
di ottimizzazione presenti si distinguono fra multi-obiettivo e mono-obiettivo, ciascuno con differenti caratteristiche in termini di accuratezza, robustezza, efficienza [5].
L’interfaccia grafica del programma presenta 3 ambienti che definiscono le 3 fasi, in
cui viene suddiviso il processo di ottimizzazione.
Work Flow In questa finestra si stabiliscono:
• il modello, riportando le variabili che lo descrivono;
• gli obiettivi che si vogliono perseguire;
• l’application, ovvero il programma di calcolo interfacciato a modeFRONTIER;
• eventuali vincoli, a cui l’individuo ottimo deve attenersi;
• l’algoritmo di ottimizzazione, cioè lo scheduler che si intende utilizzare;
• il set di individui che si intende analizzare in partenza.
119
A.1 modeFRONTIER
In figura A.1 è riportato un esempio di work-flow, con cui si è definito il problema di
ottimizzazione bidimensionale.
Figura A.1 Esempio di work flow
Si possono osservare due flussi. Il flusso logico va dalla generazione del D.O.E.
(Design of Experiments), ovvero il set di individui iniziali, all’application (denominata operaio 2D). Il flusso di dati parte invece dalla definizione dei nodi delle input variables e delle eventuali transfer variables per giungere come input per l’application.
La generazione del D.O.E. è affidata a diversi algoritmi, tra i quali si è scelto per
i primi lanci Sobol, vista la capacità di coprire al meglio e omogeneamente l’intero
dominio delle variabili decisionali.
A questo punto all’utente è richiesto di definire la strategia di ottimizzazione, ovvero di scegliere il tipo di scheduler. modeFRONTIER offre al livello base MOGA-II
(evoluzione del Multi Objective Genetic Algorithm), SIMPLEX, B-BFGS Bounded
BFGS basato sul metodo del gradiente, SA (Simulated Annealing). Vi è poi una scelta
di scheduler avanzati come MOSA, MACK, FMOGA, FSIMPLEX [22]. In questa tesi si è scelto per i primi lanci MOGA-II, Multi Objective Genetic Algorithm, data sua
robustezza e la capacità di affrontare problemi multi-obiettivo. Per i successivi lanci di
120
I Software Utilizzati.
affinamento, si è invece utilizzato un algoritmo mono-obiettivo, previa trasformazione
del problema mediante l’inserimento di una utility function.
Run Logs Sotto questa tendina durante la fase di running il software dialoga con
l’utente, fornendo un resoconto sul numero di design in corso, sui tempi di calcolo, su
eventuali errori etc.
Design Space È questa la sezione dedicata al post-processing. I dati forniti dalla
fase computazionale vengono visualizzati a schermo mediante opportuni grafici. In
questo modo è possibile valutare la bontà del processo appena compiuto e quale sia la
direzione di ottimizzazione che il software vuole fornire al disegnatore e infine compiere una analisi statistica avanzata dei risultati. Sempre in questa sezione sono messi a disposizione dell’utente procedimenti passo-passo per costruire utility function,
attraverso il M.C.D.M. (Multi Criteria Decision Making).
A.1.2
Algoritmo Genetico: MOGA-II
L’algoritmo MOGA-II si è rivelato una delle risorse principali di questo lavoro di tesi.
Esso interpreta l’ottimizzazione al pari del processo evolutivo con cui si modificano
nel tempo gli esseri viventi. Ad ogni individuo, indicato con una stringa binaria in luogo del DNA, viene associata una fitness, ovvero un punteggio in base alla particolare
posizione nella scala degli obiettivi. Il processo evolutivo segue le operazioni di:
cross over ovvero una ricombinazione genetica tra due individui, le cui stringhe vengono spezzate casualmente e ricombinate. In questo modo ciascuno dei figli
eredita alcuni geni da entrambi i genitori;
mutazione ovvero la casuale variazione di uno o più parametri di una stringa. L’utilizzo di tale metodo aumenta la probabilità di esplorare un numero maggiore di
configurazioni;
selezione ovvero viene consentita la sopravvivenza dei soli individui migliori nel
passaggio da una generazione alla successiva.
In tale modo gli individui del D.O.E. iniziale si evolvono in più generazioni proponendo nuove stringe di D.N.A, tali da generare fitness che più si avvicinino agli obiettivi
prefissati. Ciascuno degli eventi con cui avviene l’evoluzione della stringa ha una
certa probabilità di avverarsi, definibile dall’utente.
MOGA-II è una particolare versione di MOGA, algoritmo genetico presente nelle
versioni precedenti di modeFRONTIER; rispetto a MOGA c’è un miglioramento dell’elitismo, ovvero la capacità di scegliere di generazione in generazione gli individui
migliori senza però avere replicanti (memoria storica) e l’introduzione del cross-over
direzionale, ovvero una particolare combinazione fra due genitori che si può vedere
come una somma vettoriale pesata. Rispetto agli algoritmi mono-obiettivo, MOGA-II
121
A.1 modeFRONTIER
è caratterizzato da una maggiore robustezza: assicura cioè che la soluzione trovata
sia in posizione di minimo/massimo assoluto e non locale. Per contro questo algoritmo è meno accurato ed efficiente di un mono-obiettivo, ovvero assicura una minor
precisione e una minore velocità di convergenza.
A.1.3
Algoritmo mono-obiettivo: SIMPLEX
Una volta colta la posizione stabile di estremo assoluto con l’algoritmo genetico, con
ulteriori lanci mediante un algoritmo mono-obiettivo quale SIMPLEX si può ottenere
in modo efficiente una soluzione più accurata. SIMPLEX è un algoritmo per l’ottimizzazione di problemi non lineari. Detto n il numero dei gradi di libertà del modello,
il simplesso [5] è una poliedro formato da n + 1 vertici in uno spazio n-dimensionale.
In modeFRONTIER i vertici sono n+1 individui linearmente indipendenti del D.O.E..
Lo scheduler valuta i design ai vertici e deforma in modo iterativo il poliedro iniziale
fino a raggiungere la soluzione ottimale. I movimenti consentiti sono:
riflessione ovvero nello spazio n + 1-dimensionale viene eseguita una riflessione del
vertice, cui corrisponde il valore peggiore dell’obiettivo, rispetto allo spazio ndimensionale formato dai restanti punti. Se il nuovo vertice assume un miglior
valore dell’obiettivo, questo va a sostituire il vertice di partenza;
espansione ovvero se la riflessione ha prodotto un vertice con l’obiettivo migliore, è lecito attendersi che proseguendo in quel verso e direzione si ottengano
individui migliori del precedente;
contrazione ovvero, se al contrario, la riflessione ha portato ad un punto peggiore di
tutti gli altri meno che del vertice di partenza, è possibile che l’ottimo di trovi in
un punto interno al poliedro per cui dal punto riflesso ci si addentra nella figura.
Per la sua architettura SIMPLEX è un algoritmo mono-obiettivo.
A.1.4
Definizione di soluzione migliore
Dalle design chart presentate, i.e. figure 4.1 e 4.2, appare evidente come, in un problema multi-obiettivo con obiettivi fra loro in contrasto, esistano più soluzioni, comunque equivalentemente valide. È questo il caso del presente lavoro e non è possibile
trovare un unico ottimo. Ciononostante è possibile individuare un set di individui che
prevalgono gli altri, il cosiddetto fronte di Pareto, vedi figure A.2 e A.3. Si parla
di dominanza di un gruppo di individui e formalmente, detta fi la i-esima funzione
obiettivo, il design x domina il design y (si scriva x > y) se e solo se:
(∀i fi (x) > fi (y)) ∩ (∃ j | f j (x) > f j (y))
Per cui x domina y se almeno per un aspetto lo supera. Se nobj è il numero degli
obiettivi, il fronte di Pareto è un ente geometrico di dimensione nobj − 1.
122
I Software Utilizzati.
Figura A.2 Esempio di fronte di Pareto (individui evidenziati): gli individui del fronte
dominano i restanti.
Figura A.3 Ulteriore esempio di fronte di Pareto: entrambi gli obiettivi vanno massimizzati.
A.1.5
M.C.D.M.
Spesso i problemi di ottimizzazione sono di tipo multi-obiettivo. Come accennato
però, un algoritmo efficiente ed accurato quale SIMPLEX è di tipo mono-obiettivo.
Inoltre la definizione di un’unica soluzione migliore di tutte le altre richiede la massimizzazione di un’unica funzione, in cui un obiettivo venga privilegiato rispetto ad
altri(2) . Può risultare quindi necessario trasformare un problema multi-obiettivo in
(2) Il
che equivale ad esplorare una specifica zona dello spazio delle soluzioni.
A.2 FEMLAB
123
un problema mono-obiettivo attraverso una utility function, ovvero una somma pesata
degli obiettivi originari:
Fobj = ∑ pi ϕi
dove Fobj è il nuovo obiettivo da rendere estremo, pi è il peso dell’obiettivo i-esimo
e ϕi è l’obiettivo i-esimo della multi-obiettivo. I pesi sono caratterizzati dall’avere
somma unitaria e lo strumento M.C.D.M. (Multi Criteria Decision Making) interno
a modeFRONTIER aiuta la generazione di una utility function con funzioni di peso
lineari o non-lineari. Per una breve descrizione del funzionamento di questa opzione
si veda sezione 4.4.
A.2
FEMLAB
FEMLAB è il software utilizzato per la soluzione numerica del problema termicoconvettivo affrontato in questa tesi. FEMLAB è un pacchetto agli elementi finiti
(FEM) con generatore di mesh non-strutturata, per la risoluzione di problemi riconducibili a sistemi di equazioni differenziali del secondo ordine alle derivate parziali
(PDE). Rispetto le versioni precedenti, quella in uso, FEMLAB 3.0, non è più solo un
toolbox di MATLAB, pur mantenendo inalterata la piena integrazione con questo applicativo. Analogamente ai lavori [5, 10] si è utilizzato FEMLAB per la possibilità di
accedere direttamente alle equazioni risolutive, adattandole entro certi limiti alla proprie esigenze, e di interfacciarsi con un linguaggio di programmazione agevole quale
MATLAB.
FEMLAB annovera inoltre una serie di modelli, application mode, di svariati campi fisici, che consentono la definizione di uno specifico problema anche senza entrare
nella modalità di scrittura delle PDE. Tra questi campi si possono citare: modulo per
meccanica strutturale, elettromagnetismo, conduzione del calore, fluidodinamica e
quant’altro [18].
Infine FEMLAB permette la definizione di problemi in multiphysics mode, i.e.
multi-fisica, ovvero è possibile accoppiare nello stesso modello più sotto-modelli dalla
fisica differente, i.e. interazione fluido-struttura in turbomacchine.
A.2.1
Caratteristiche essenziali
I passi per la risoluzione di un modello sono i seguenti:
1. definizione geometrica del dominio;
2. definizione delle equazioni differenziali che governano il modello;
3. imposizione delle condizioni al contorno (boundary conditions);
4. definizione dei parametri di griglia e generazione della mesh;
124
I Software Utilizzati.
5. scelta del risolutore, Solver, (lineare, non-lineare, non-stazionario, agli autovalori, parametrico lineare o non-lineare, adattativo);
6. risoluzione del problema;
7. eventuale fase di post-processing per l’estrapolazione dei risultati.
Per maggiori dettagli si veda [18].
Definizione della geometria
FEMLAB consente di generare modelli geometrici anche complessi 2D e 3D, oppure
di importare strutture da CAD esterni. Il CAD interno offre la possibilità di effettuare
modellazione solida o per superfici a livello base e un certo numero di informazioni
sulla geometria possono essere agevolmente estratte.
PDE general form
In questo lavoro per la scrittura delle equazioni fluidodinamiche e termiche, si è
utilizzata la formulazione diretta delle PDE nella cosiddetta general form (A.1).

∂ ui


in Ω

 da i ∂t + ∇ · Γi = Fi
∂Rj
(A.1)
−n · Γi = Gi +
µ j su ∂ Ω


∂
u

i

0 = Rj
su ∂ Ω
dove
- ui è l’i-esima variabile, u = (u1 , ..., uN );
- da i è l’i-esimo coefficiente di massa;
- Γi è il cosiddetto flux vector;
- Fi è il termine forzante;
- n è il versore normale al boundary;
- Gi è un coefficiente per l’imposizione delle BC;
- R j è un coefficiente per l’imposizione del j-esimo vincolo, R = (R1 , ..., RM );
- µi è il j-esimo moltiplicatore di Lagrange, i moltiplicatori sono in numero pari
ai vincoli.
La prima di A.1 è l’equazione differenziale alle derivate parziali, le seconda e la terza
rappresentano rispettivamente le condizioni al contorno di Neumann e Dirichlet.
A.2 FEMLAB
125
Mesh
La versione in uso consente di generare una griglia non-strutturata che si traduce nel
caso bi-dimensionale in una mesh triangolare, mentre nel caso tri-dimensionale in una
mesh a tetraedri. FEMLAB consente di avere un controllo sulla qualità degli elementi
a livello di dominio e sotto-domini, quali bordi, facce, punti.
Solver
In questo lavoro si è utilizzato unicamente il solver non-lineare stazionario (femnlin),
capace di risolvere anche problemi lineari riconoscendo un peso elevato alla parte
diffusiva della prima equazione di A.1. Si ricorda che un modello non-lineare è tale
se descritto da una o più equazioni differenziali, in cui uno, o più coefficienti, delle
derivate è non-lineare. Il solver, che può essere diretto o iterativo [18, 19], provvede a
discretizzare le equazioni come
L − NT Λ = 0
(A.2)
M=0
dove L, N e M sono funzioni del vettore soluzione U. Il sistema A.2 viene risolto
con una variante del metodo di Newton. Da una soluzione iniziale U0 si calcola una
correzione δ U attraverso (A.3)
Kδ U + N T Λ = L
(A.3)
Nδ U = M
con K, N, L e M calcolate per U = U0 . Il nuovo vettore soluzione è ora U = U0 +λ δ U,
dove λ ∈ [0; 1] è un coefficiente di smorzamento (damping factor). Il procedimento è
iterato fino al livello di convergenza desiderato.
Per le altre tipologie di solver si faccia riferimento a [18].
Post-processing
FEMLAB fornisce molteplici strumenti per la fase di post-processing e visualizzazione risultati. Consente di creare grafici quali:
- surface plots
- isosurface plots
- contour plots
- streamline plots
nonché combinazioni di questi.
Bibliografia
[1] Michael E. Mortenson. Modelli geometrici in computer graphics. McGraw-Hill,
Milano, 1989.
[2] Gerald Farin. Curves and Surfaces in Computer-Aided Geometric Design: A
Practical Guide. ACADEMIC PRESS, Milano, fourth edition, 1997.
[3] Andrea Rimondi. Generazione di configurazioni aerodinamiche mediante
NURBS (Non-Uniform Rational B-Spline). Tesi di Laurea, Università degli
Studi di Pisa, Facoltà di Ingegneria, Dipartimento di Ingegneria Aerospaziale,
2002/2003.
[4] S.V.Patankar C.H.Liu E.M.Sparrow. Fully developed flow and heat transfer in ducts having streamwise-periodic variations of cross-sectional area.
Transactions of the ASME, 99, May 1977.
[5] Luca Fogar. Ottimizzazione termofluidodinamica di sistemi per lo scambio
termico convettivo. Tesi di Laurea, Università degli Studi di Trieste, Facoltà
di Ingegneria, Dipartimento di Ingegneria Navale, del Mare e per l’Ambiente,
2000/2001.
[6] J.H.Ferziger M.Peric. Computational Methods for Fluid Dynamics. Springer,
2002.
[7] N.G.Hagjicostantinou O.Simek. Constant-wall-temperature nusselt number in
micro and nano-channels. J. of Heat Transfer, 124:356–364, 2002.
[8] R.W.Barber D.R.Emerson. The influence of knudsen number on the hydrodynamic development length within parallel plate micro-channels. Advances in Fluid
Mechanics IV, 32:207–216, 2002.
[9] D.R.Sawyers M.Sen H.C.Chang.
Heat transfer enhancement in threedimensional corrugated channel flow. Int. J. Heat Mass Transfer, 41:3559–3573,
1998.
127
128
BIBLIOGRAFIA
[10] Francesco Pinto. Ottimizzazione di forma in problemi inversi di trasmissione del
calore. Tesi di Laurea, Università degli Studi di Trieste, Facoltà di Ingegneria,
Dipartimento di Ingegneria Navale, del Mare e per l’Ambiente, 2002/2003.
[11] Carlo Vecile. Appunti delle lezioni di fluidodinamica. LINT, 1998.
[12] R.K.Shah A.L.London. Laminar Flow Forced Convection in Ducts, A source
book for compact Heat Exchanger Analytical Data. ACADEMIC PRESS, New
York, San Francisco, London, 1978.
[13] B.Niceno E.Nobile. Numerical analysis of fluid flow and heat trasfer in periodic wavy channels. International Journal of HEAT AND FLUID FLOW, 23
December 2000.
[14] A.Barletta E.Zanchini. The existence of an asymptotic thermally developed region for laminar forced convection in a circular duct. Int. J. Heat Mass Transfer,
39(13):2735–2744, 1996.
[15] C.Nonino G.Comini. Finite-element analysis of convection problems in spatially periodic domains. Numerical Heat Tranfer, Taylor & Francis, Part
B(34):361:378, 1998.
[16] E.Funaioli A.Maggiore U.Meneghetti. Lezioni di meccanica applicata alle
macchine, volume 2. Patron Editore, Bologna, 1988.
[17] G.Wang S.P.Vanka. Convective heat transfer in periodic wavy passages. International Journal of Heat and Mass Tranfer, 38(17):3219–3230, 19 January
1995.
[18] COMSOL FEMLAB version 3 Documentation.
http://www.comsol.com.
Si veda anche URL
[19] M.Ciofalo G.Comini G.Croce E.Nobile C.Nonino.
Fondamenti di TERMOFLUIDODINAMICA COMPUTAZIONALE.
SGE editoriali, Padova,
2001.
[20] M.Greiner P.F.Fisher H.M.Tufo R.A.Wirtz. Three dimensional simulations of
enhanced heat transfer in a flat passage downstream from a grooved channel.
Proceedings of NHTC’00, 34th National Heat Transfer Conference, Pittsburgh,
Pennsylvania, August 2000.
[21] Singiresu Rao. Engineering Optimization, Theory and Practice.
Interscience Publication. John Wiley & Sons, Inc., third edition, 1996.
[22] modeFRONTIER version 3 Documentation.
http://www.esteco.com.
Wiley-
Si veda anche URL
BIBLIOGRAFIA
129
[23] Marco Manzan. Utilizzo di codici commerciali per lo studio dello scambio termico in geometrie periodiche. 18◦ Congresso Nazionale sulla Trasmissione del
Calore -U.I.T.- Cernobbio, I:371–382, 2000.

Scarica

universit `a degli studi di trieste parametrizzazione geometrica ed

universit `a degli studi di trieste parametrizzazione geometrica ed

PowerPoint

AI CITTADINI DI CANALE MONTERANO In questi giorni si sono

Gianluca Ferranti, Reitek

Unità di lavoro 1 - Elementi strutturali della comunicazione

gigi

Lezione n. 5

Bollettino tecnico ADS-R8 - ASA-RT

COMUNICATO STAMPA Canale ARTE e la Divina Marchesa

COMUNICATO STAMPA Canale ARTE cattura la bellezza nelle