Sensitivita aerodinamica nelle autovetture` da competizione: uno … · 2017-07-13 · corpi con...

Politecnico di Milano

FACOLTA DI INGEGNERIA INDUSTRIALECorso di Laurea Magistrale in Ingegneria Aeronautica

Sensitivita aerodinamica nelle autovettureda competizione: uno studio basato su

codice open source e operatore aggiunto

Relatore:Prof. Maurizio Quadrio

Correlatori:Ing. Fabrizio TessiciniIng. Silvio Lorenzo

Tesi di Laurea di:Simone Battaini

Matr. 765130

Anno Accademico 2012/2013

SOMMAR IO

In questa tesi ci si propone di sfruttare il concetto dell’operatore ag-giunto nell’ambito dell’ottimizzazione di forma in flussi esterni, al finedi ricavare l’andamento della sensitivita di un corpo rispetto alla suadeformazione e relativamente ad una funzione obiettivo, sia su casibidimensionali che tridimensionali. Per ottenere questo scopo, vieneutilizzato il modello delle equazioni RANS, nonche il programma opensource OpenFOAM. Il calcolo dell’aggiunto e stato applicato dapprimaad un caso 2D, in cui la validita del procedimento e stata dimostratatramite la scrittura di un algoritmo di deformazione per il profilo NA-CA preso in esame, e sia in regime di moto laminare che turbolento.L’analisi si e poi spostata al 3D, dove il metodo e stato applicato acorpi con forma aerodinamica di varia natura: anche qui la sensitivitye stata estratta con successo, ed ha permesso sia di individuare lezone del corpo piu sensibili ad una modifica, che la direzione in cuieseguire la deformazione stessa, al fine di ottenere un miglioramentodella prestazione desiderata.

Parole chiave: operatore aggiunto, RANS, OpenFOAM

i

ABSTRACT

The purpose of this thesis is to exploit the concept of the adjoint opera-tor in the field of the shape optimization in external flows, in order toobtain the sensitivity of the deformation of a body with respect to anobjective function, both with two-dimensional and three-dimensionalcases. To achieve this purpose, RANS model has been used, togetherwith the open source program OpenFOAM. The calculation of theadjoint has been first applied on a 2D case, where the validity ofthe procedure has been demonstrated by developing a deformationscript for the considered NACA airfoil, both in laminar and turbulentregime. The analysis then moved to the 3D, where the method hasbeen applied to various aerodynamic shape bodies: the sensitivityhas been successfully extracted also in this case, and allowed either tohighlight the zones of the body more suitable to a modification, andthe direction on which perform the deformation, in order to obtain animprovement for the desired performance.

Keywords: adjoint operator, RANS, OpenFOAM

iii

R INGRAZ IAMENT I

Vorrei innanzitutto ringraziare la Scuderia Toro Rosso per avermipermesso di vivere l’esperienza unica di far parte di un team di F1,e per avermi messo nelle condizioni di lavorare nel miglior modopossibile. La mia gratitudine va soprattutto a Fabrizio, i cui consigli eindicazioni si sono rivelati sempre corretti e propositivi.

Desidero allo stesso modo esprimere riconoscenza al Prof. MaurizioQuadrio che, seppur con tutte le difficolta che nascono da una “comu-nicazione a distanza”, e sempre riuscito a dedicarmi il suo tempo ead indirizzarmi nella retta via, e non di meno per avermi scelto tra ipossibili candidati di questo stage.

Un grazie particolare va ai nuovi amici: innanzitutto a Jacques, concui e nata una sincera amicizia, e poi a Marco1 e a Matteo, senza il cuiaiuto alcuni problemi sarebbero ancora probabilmente irrisolti. Tantariconoscenza va anche a Marcello, per l’aiuto che mi ha dato primaancora di conoscermi, e grazie anche a chiunque abbia scambiato dueparole o una risata con me2.

Come poter dimenticare voi, vecchi compagni con cui ho condivisoquesta avventura che e il Polimi, e con cui ho passato alcuni anni trai piu duri ma anche tra i piu divertenti e importanti della mia vita:Giulio, Marco e Mega c©, del nostro gruppo iniziale siamo “rimasti”soltanto noi, e ora che siamo quasi alla fine del nostro percorso distudi non posso fare altro che augurarvi ogni bene. Boga, Foranz,Mitch, Pino e Rondine, ci siamo conosciuti a meta percorso, e quantomai non prima! Ci saremmo sicuramente divertiti molto piu di quantonon abbiamo fatto3!

... and last but not least i ringraziamenti piu importanti: grazie a tuttala mia famiglia, in particolare a mamma Rosy e a papa Fausto, chesono da sempre due pilastri nella mia vita: se oggi posso camminare atesta alta lo devo in gran parte a voi; grazie anche per avermi sempresUpportato (piu che sOpportato) e per non avermi mai tagliato leali, accettando ogni mia scelta. Stefano, ti auguro con tutto me stes-so di riuscire a realizzare ogni tuo sogno, non smettere di crederci mai!

... E soprattutto grazie a te Chiara! Grazie per essermi stata accantoin questo periodo di lontananza con una tenerezza indescrivibile(lo so benissimo quanto questo ti sia costato), grazie per essere lapersona che condivide la sua vita con me da quasi 8 anni, e grazieper essere tutto cio che in una donna si puo desiderare. Sei unica nondimenticarlo mai!

1 “Mettiti nei miei panni che lo sopporto per 600 km alla settimana!!!”2 ...c’e ancora gente che si chiede che cosa ci facesse quel ragazzo con degli hamburger in tasca...3 Beh, da un certo punto di vista e un bene essersi conosciuti dopo, altrimenti staremmo ancora

tutti studiando impianti..

v

IND ICE

1 introduzione . . . . . . . . . . . . . . . . . . . . . . . . 1

1.1 Scopo e struttura del lavoro di tesi . . . . . . . . . 4

2 teoria . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

2.1 Modello fisico-matematico . . . . . . . . . . . . . . 7

2.2 Equazioni aggiunte laminari . . . . . . . . . . . . . 9

2.3 Equazioni aggiunte turbolente . . . . . . . . . . . 11

2.4 Cost function . . . . . . . . . . . . . . . . . . . . . 13

2.5 Calcolo della sensitivity . . . . . . . . . . . . . . . 13

2.6 Approssimazione alla Lohner . . . . . . . . . . . . 15

3 implementazione . . . . . . . . . . . . . . . . . . . . . 17

3.1 Presentazione del software OpenFOAM . . . . . . 17

3.2 Solver per il sistema aggiunto . . . . . . . . . . . . 20

3.3 Script Python . . . . . . . . . . . . . . . . . . . . . 20

3.4 Loop di ottimizzazione in Bash . . . . . . . . . . . 21

3.5 Simulazioni e criteri di arresto . . . . . . . . . . . 22

4 casi 2d . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

4.1 Cost function . . . . . . . . . . . . . . . . . . . . . 25

4.2 Caso 2D laminare . . . . . . . . . . . . . . . . . . . 26

4.2.1 Mesh . . . . . . . . . . . . . . . . . . . . . . 26

4.2.2 Condizioni al contorno . . . . . . . . . . . 27

4.2.3 Simulazione diretta . . . . . . . . . . . . . . 27

4.2.4 Simulazione aggiunta . . . . . . . . . . . . 30

4.2.5 Deformazione del profilo . . . . . . . . . . 31

4.2.6 Loop . . . . . . . . . . . . . . . . . . . . . . 33

4.3 Caso 2D turbolento . . . . . . . . . . . . . . . . . . 38

4.3.1 Mesh . . . . . . . . . . . . . . . . . . . . . . 38

4.3.2 Modello di turbolenza . . . . . . . . . . . . 38




4.3.6 Deformazione del profilo . . . . . . . . . . 42

4.3.7 Loop . . . . . . . . . . . . . . . . . . . . . . 44

5 casi 3d . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49

5.1 snappyHexMesh . . . . . . . . . . . . . . . . . . . 49

5.2 ANSA . . . . . . . . . . . . . . . . . . . . . . . . . . 50

5.3 checkMesh . . . . . . . . . . . . . . . . . . . . . . . 50

5.4 Schemi numerici . . . . . . . . . . . . . . . . . . . 52

5.5 Caso 3D con wing e mesh poliedrica . . . . . . . . 53




5.6 Caso 3D con wing e trim-mesh . . . . . . . . . . . 61


5.7 Rear Model . . . . . . . . . . . . . . . . . . . . . . 64

vii

Indice



6 conclusioni e sviluppi futuri . . . . . . . . . . . . . 67

a coefficienti di forza aggiunti . . . . . . . . . . . . 71

b algoritmo di openfoam . . . . . . . . . . . . . . . . . 75

Bibliografia . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81

viii

ELENCO DELLE F IGURE

Figura 1 Andamento tipico in direzione upstream dellavelocita aggiunta in un caso 2D . . . . . . 11

Figura 2 Struttura di un caso in OpenFOAM . . . 18

Figura 3 Dominio nel caso bidimensionale . . . . . 26

Figura 4 Zoom della mesh attorno al bordo d’attacco delprofilo NACA 0012 . . . . . . . . . . . . . 26

Figura 5 Coefficiente di portanza adimensionale nel caso2D laminare . . . . . . . . . . . . . . . . . 28

Figura 6 Coefficiente di resistenza adimensionale nelcaso 2D laminare . . . . . . . . . . . . . . 28

Figura 7 Campo di velocita attorno al profilo NACA0012 nel caso 2D laminare . . . . . . . . . 29

Figura 8 Cp attorno al profilo NACA 0012 nel caso 2Dlaminare . . . . . . . . . . . . . . . . . . . 29

Figura 9 Residui della simulazione aggiunta nel caso 2Dlaminare . . . . . . . . . . . . . . . . . . . 30

Figura 10 Andamento della sensitivity sul profilo nel caso2D laminare . . . . . . . . . . . . . . . . . 31

Figura 11 Deformazione del profilo nel caso 2D lamina-re . . . . . . . . . . . . . . . . . . . . . . . 32

Figura 12 Deformazione del profilo, zoom al bordo d’at-tacco nel caso 2D laminare . . . . . . . . . 32

Figura 13 Deformazione del profilo, zoom al bordo d’u-scita nel caso 2D laminare . . . . . . . . . 33

Figura 14 Aumento del CL durante l’ottimizzazione nelcaso 2D laminare . . . . . . . . . . . . . . 33

Figura 15 Riduzione del CD durante l’ottimizzazione nelcaso 2D laminare . . . . . . . . . . . . . . 34

Figura 16 Evoluzione dell’aumento del CL nel caso 2Dlaminare . . . . . . . . . . . . . . . . . . . 35

Figura 17 Evoluzione della riduzione del CD nel caso 2Dlaminare . . . . . . . . . . . . . . . . . . . 35

Figura 18 Evoluzione dell’aumento dell’efficienza nel ca-so 2D laminare . . . . . . . . . . . . . . . 36

Figura 19 Confronto del Cp sul dorso del profilo nel caso2D laminare . . . . . . . . . . . . . . . . . 37

Figura 20 Confronto del Cp sul ventre del profilo nel caso2D laminare . . . . . . . . . . . . . . . . . 37

Figura 21 Riduzione della funzione obiettivo nel caso 2Dlaminare . . . . . . . . . . . . . . . . . . . 37

Figura 22 Coefficiente di portanza adimensionale nel caso2D turbolento . . . . . . . . . . . . . . . . 39

Figura 23 Coefficiente di resistenza adimensionale nelcaso 2D turbolento . . . . . . . . . . . . . 40

Figura 24 Campo di velocita attorno al profilo NACA0012 nel caso 2D turbolento . . . . . . . . 40

ix

Elenco delle figure

Figura 25 Cp attorno al profilo NACA 0012 nel caso 2Dturbolento . . . . . . . . . . . . . . . . . . 40

Figura 26 Residui della simulazione aggiunta nel caso 2Dturbolento . . . . . . . . . . . . . . . . . . 41

Figura 27 Andamento della sensitivity sul profilo nel caso2D turbolento . . . . . . . . . . . . . . . . 42

Figura 28 Deformazione del profilo nel caso 2D turbolen-to . . . . . . . . . . . . . . . . . . . . . . . 42

Figura 29 Deformazione del profilo, zoom al bordo d’at-tacco nel caso 2D turbolento . . . . . . . . 43

Figura 30 Deformazione del profilo, zoom al bordo d’u-scita nel caso 2D turbolento . . . . . . . . 43

Figura 31 Aumento del CL durante l’ottimizzazione nelcaso 2D turbolento . . . . . . . . . . . . . 44

Figura 32 Riduzione del CD durante l’ottimizzazione nelcaso 2D turbolento . . . . . . . . . . . . . 44

Figura 33 Evoluzione dell’aumento del CL durante l’otti-mizzazione nel caso 2D turbolento . . . . 45

Figura 34 Evoluzione della riduzione del CD durante l’ot-timizzazione nel caso 2D turbolento . . . 45

Figura 35 Evoluzione dell’aumento dell’efficienza duran-te l’ottimizzazione nel caso 2D turbolento 46

Figura 36 Confronto del Cp sul dorso del profilo nel caso2D turbolento . . . . . . . . . . . . . . . . 46

Figura 37 Confronto del Cp sul ventre del profilo nel caso2D turbolento . . . . . . . . . . . . . . . . 47

Figura 38 Riduzione della funzione obiettivo nel caso 2Dturbolento . . . . . . . . . . . . . . . . . . 47

Figura 39 Definizione di skewness . . . . . . . . . . . 51

Figura 40 Definizione di nonOrthogonality . . . . . . 51

Figura 41 Mesh poliedrica, zoom attorno all’ala . . . 53

Figura 42 Mesh poliedrica, zoom attorno all’estremita ala-re . . . . . . . . . . . . . . . . . . . . . . . 54

Figura 43 Coefficiente di portanza adimensionale nel caso3D poliedrico . . . . . . . . . . . . . . . . 55

Figura 44 Coefficiente di resistenza adimensionale nelcaso 3D poliedrico . . . . . . . . . . . . . 55

Figura 45 Residui della simulazione aggiunta nel caso 3Dpoliedrico . . . . . . . . . . . . . . . . . . . 57

Figura 46 ∆Ux iterazioni 2500− 3000 . . . . . . . . . 58

Figura 47 ∆Ux iterazioni 3000− 3500 . . . . . . . . . 58

Figura 48 Andamento dell’integrale della sensitivity nelleprime 150 iterazioni . . . . . . . . . . . . . 59

Figura 49 Andamento dell’integrale della sensitivity tral’iterazione 1000 e 4000 . . . . . . . . . . . 59

Figura 50 Andamento della sensitivity sul dorso del pro-filo - mesh poliedrica . . . . . . . . . . . . 60

Figura 51 Andamento della sensitivity sul ventre del pro-filo - mesh poliedrica . . . . . . . . . . . . 60

Figura 52 Residui della simulazione aggiunta nel caso 3Dcon trim-mesh . . . . . . . . . . . . . . . . 62

x

Elenco delle figure

Figura 53 Andamento dell’integrale della sensitivity tral’iterazione 1 e 3500 . . . . . . . . . . . . . 63

Figura 54 Andamento della sensitivity sul dorso del pro-filo - trim-mesh . . . . . . . . . . . . . . . . 63

Figura 55 Andamento della sensitivity sul ventre del pro-filo - trim-mesh . . . . . . . . . . . . . . . . 63

Figura 56 Andamento dei residui nel caso 3D con meshpoliedrica . . . . . . . . . . . . . . . . . . . 72

Figura 57 Coefficiente di “resistenza” aggiunto nel caso3D con mesh poliedrica . . . . . . . . . . . 72

Figura 58 Coefficiente di “portanza” aggiunto nel caso3D con mesh poliedrica . . . . . . . . . . . 73

xi

1 INTRODUZ IONE

L’operatore aggiunto (o duale) relativo ad un certo operatore diffe-renziale e uno strumento molto comune in analisi funzionale, ed ilprocedimento piu semplice per ricavarlo e analogo ad un’integrazioneper parti (per ulteriori dettagli si veda ad esempio [1]). Esso assumeinoltre particolare rilievo nel momento in cui viene utilizzato nell’ambi-to dell’ottimizzazione, nonche in relazione al sistema di equazioni cheregge il moto di un fluido, ossia le ben note equazioni di Navier-Stokes.

La ricerca dell’ottimo e un argomento che puo essere affrontatoin piu campi anche molto diversi fra loro tuttavia, focalizzandosinell’ambito aerodinamico, essa puo essere intesa come la ricerca di unageometria capace di minimizzare (o massimizzare) una certa funzioneobiettivo relativamente all’imposizione di alcuni vincoli; ci sono quindimolti modi per perseguire questo scopo e, tra i vari procedimenti ditipo “gradient-based”, il metodo che sfrutta l’operatore aggiunto e unotra quelli piu accurati e computazionalmente economici.

Una rassegna storica su alcuni degli approcci che sono stati utiliz-zati nel corso degli anni per ottenere l’andamento della sensitivitadi un corpo rispetto alla sua deformazione e relativamente ad unafunzione obiettivo (anche chiamata sensitivity) puo essere trovata nel-l’articolo [2]: storicamente, il primo metodo che e stato utilizzato equello relativo alle differenze finite; in questo caso viene valutata lavariazione del campo di moto e delle quantita ad esso correlate otte-nute tramite una modifica arbitraria della geometria del corpo. Taleprocedimento viene ripetuto per ogni variabile di design considerata,e quindi la sensitivity viene assemblata. Il prezzo della semplicita delmetodo viene pagato con il costo computazionale correlato al calcolodel gradiente, stimato nell’ordine di 2 • nα soluzioni delle equazioni diNavier-Stokes (non lineari), dove α e proprio il numero di variabili didesign che viene utilizzato; si capisce quindi che, per casi complessi,questo metodo risulti in pratica inutilizzabile. Un altro punto criticomenzionato e correlato alla scelta dell’ordine di grandezza delle va-riazioni delle variabili di design δα: se tali valori sono troppo piccoliallora i cambiamenti nella geometria saranno poco significativi, maqualora essi siano troppo grandi, invece, la linearizzazione tramite l’e-spansione di Taylor che viene introdotta per il calcolo della sensitivitaviene invalidata. Questo metodo e stato utilizzato negli anni ′70 peril progetto di profili (i primi studi sono stati effettuati da Hicks [3] eVanderPlaats [4]) ma, con l’avvento dell’ottimizzazione applicata aicasi 3D, ha lasciato gradualmente posto ad altri, computazionalmentepiu efficienti.

Un secondo approccio (piu recente, siccome la sua introduzionerisale agli anni ′80, quando e stato applicato per la prima volta daBristow e Hawk [5]) consiste nel cosiddetto discrete direct method, doveil residuo delle equazioni di Navier-Stokes viene differenziato rispettoalle stesse variabili α del caso precedente: in questo caso il calcolo del

1

introduzione

gradiente riesce ad essere effettuato tramite un’unica soluzione nonlineare del flusso e ad nα soluzioni di sistemi lineari per ricavare lavariazione delle equazioni di Navier-Stokes rispetto alle variabili didesign. Qui, il costo di calcolo del gradiente risulta essere piu bassorispetto al caso precedente, ma ancora dipendente da α.

Viene dunque introdotto negli anni ′90 il metodo dell’operatore ag-giunto discreto, calcolato a partire dalle equazioni di Navier-Stokes giadiscretizzate (discretize-then-derive), e viene dimostrato come, in questocaso, il calcolo del gradiente sia indipendente dal vettore contenente levariabili di design α, la cui presenza affliggeva invece i due approcciprecedenti (le prime applicazioni di questo concetto sono dovute aShubin e Franck, nel 1991 [6]).

Parallelamente cominciano gli studi anche sull’operatore aggiuntocontinuo, che e quello utilizzato in questa tesi, ricavato a partire dalleequazioni dirette continue e successivamente discretizzato (derive-then-discretize) (i primi lavori con questo metodo si devono a Pironneau [7]e Jameson [8]). Anche in questo caso, come per l’operatore aggiuntodiscreto, l’ottenimento del suddetto gradiente passa attraverso il solocalcolo di una simulazione diretta (ossia con le equazioni di Navier-Stokes canoniche) e di una simulazione aggiunta (con le equazionidi Navier-Stokes aggiunte, ricavate a partire dalle corrispondentiequazioni dirette). Come gia precedentemente annunciato, la virtuprincipale di questo approccio risiede nell’indipendenza del suo costocomputazionale rispetto al numero di variabili di design e quindi,soprattutto per problemi complessi con un elevato numero di gradidi liberta, questo metodo risulta un’ottima scelta. Invece, il costo diuna simulazione aggiunta e dell’ordine di quello della corrispondentesimulazione diretta.

E percio possibile, al prezzo di due sole simulazioni, ottenere unacompleta mappa di sensitivita alla deformazione del corpo consideratorispetto ad una specifica funzione obiettivo. Quest’ultima formulazio-ne risulta, come ricavato in [9], meno costosa da implementare rispettoa quella discreta (in termini computazionali), sebbene l’andamentoesatto del gradiente lo si ottenga solamente nel primo caso. Tutta-via, sempre in questo articolo, vengono confrontati gli andamenti ele differenze tra le sensitivita ottenute tramite i due approcci, dimo-strando come sia possibile utilizzarle indifferentemente, in quanto leinformazioni estratte sono pressoche identiche nei due casi.

Una volta ricavato il gradiente sulla superficie del corpo (il cui an-damento per i casi considerati nella tesi e legato, come sara dettagliatomeglio in seguito, alle condizioni al contorno per il sistema aggiunto),e possibile quindi modificare la geometria sulla base dell’informazioneottenuta e, successivamente, far partire un secondo ciclo di ottimiz-zazione, composto da una nuova risoluzione delle equazioni direttee aggiunte, nonche dall’ottenimento di un nuovo andamento dellasensitivita sul corpo deformato.

Nell’ambito dell’ottimizzazione dei condotti interni, particolare rilie-vo ha la tecnica della topological optimization, ove l’operatore aggiuntoe utilizzato al fine di determinare quali celle nel dominio fluido sianofavorevoli ad essere mantenute oppure no rispetto alla funzione costoscelta. Tale argomento e stato ampiamente trattato da vari autori, tra

2

introduzione

cui in particolare occorre menzionare C. Othmer [10]. Nell’articolocitato viene messo in luce come, in ambito industriale, il processo didesign e ottimizzazione sia spesso composto da una collaborazionefra reparti diversi con un approccio di tipo “trial and error”; questometodo e ovviamente poco efficiente e percio, soprattutto per quantoriguarda l’ottimizzazione di flussi interni (un esempio nell’ambitoautomobilistico sono i condotti per l’aria condizionata piuttosto chegli scarichi), le tecniche piu consolidate comprendono l’impiego, adesempio, di algoritmi genetici. Questo modo di procedere, sebbenemolto robusto, ha fra i suoi punti deboli di nuovo un elevato costocomputazionale, che ne limita percio l’utilizzo. Con il metodo del-l’operatore aggiunto, al contrario, e possibile ottenere la sensitivitaalla funzione obiettivo indipendentemente dalla parametrizzazionedella mesh utilizzata. Uno degli obiettivi che viene ricercato, e chesi traduce nella formulazione della relativa funzione obiettivo, e adesempio la minimizzazione della perdita di pressione tra monte evalle del condotto stesso. Viene dunque introdotto il concetto dellatopological optimization, dove il criterio di ottimizzazione viene legatoall’introduzione, sia nel sistema di equazioni dirette che aggiunte,di un termine di “porosita artificiale di Darcy” (chiamato α), il cuivalore in ogni cella alla fine del ciclo iterativo permettera di capiredove e necessario agire. Infatti, ad equazioni aggiunte risolte, vienevisualizzato l’andamento di α nel dominio: quelle zone dove esso havalori piccoli si comportano “bene” relativamente alla cost functionscelta, e possono percio essere lasciate invariate, mentre le parti incui α assume valori elevati hanno un comportamento “sbagliato” evanno percio rimosse deformando la geometria a ridosso di questecelle. Tuttavia, i limiti di questa applicazione risiedono soprattutto indue punti, ossia (1) nel fatto che l’ottimizzazione e fatta a livello divolume, e quindi questo approccio ha poco senso in casi di aerodina-mica esterna dove cio che importa e la deformazione della superficie,e che (2) questo metodo porta ad una modifica del corpo abbastanzagrossolana e non sicuramente di dettaglio.

L’approccio che risulta tuttavia piu adatto per i casi di aerodinami-ca esterna, come quelli presi in considerazione, e quello della shapeoptimization, leggermente diverso dal caso precedente principalmenteper l’assenza del termine di porosita nelle equazioni, per l’espressionedella sensitivity stessa, e per le funzioni obiettivo che vengono utiliz-zate. E quindi possibile, risolvendo i due set di equazioni, ricavarel’andamento del gradiente direttamente sul corpo da ottimizzare, edeseguirne la deformazione con un accurato livello di dettaglio. Altresıin questo caso e reperibile una buona letteratura, come ad esempio [11].In questo articolo di Jameson viene ulteriormente sottolineata la ne-cessita di utilizzare metodi automatici per la determinazione di ungradiente che indichi la direzione in cui effettuare l’ottimizzazione, equesta volta in ambito specificatamente aeronautico, dove altrimentisarebbero necessarie continue prove pratiche su molte configurazionidiverse in galleria del vento. Da qui l’utilizzo dell’operatore aggiuntonell’ambito della shape optimization: esso viene ricavato correlandoloalla teoria del controllo, e in particolare con l’introduzione dei molti-plicatori di Lagrange (che e il procedimento utilizzato in questa tesi);

3

introduzione

la variabile di controllo diventa dunque la forma del corpo, mentreviene introdotto un funzionale contenente sia la cost function da mi-nimizzare che le equazioni di Navier-Stokes espresse sotto forma diresidui. Dalla minimizzazione di tale funzionale vengono dunquericavate sia le equazioni aggiunte che la sensitivita, per le equazionidi Eulero e per quelle di Navier-Stokes viscose. Particolare attenzioneviene data alla decomposizione di tale funzione obiettivo, espressatramite la somma di contributi volumetrici e superficiali, e su comeciascuno di questi due termini intervenga nelle equazioni: sebbene ilprimo vada a modificare il sistema aggiunto stesso, diventando unasorta di termine sorgente, il secondo agisce soltanto sulle condizionial contorno del corpo: in quest’ultimo caso, le equazioni risultanoindipendenti dalla cost function, ed e percio possibile utilizzare il mede-simo solutore anche per ottimizzazioni di tipo diverso semplicementemodificando, da simulazione a simulazione, il tipo di condizione alcontorno che viene imposta. Da qui si puo dedurre come particolareattenzione sia riposta in funzioni obiettivo che comportano integrali diuna data quantita solo sulla superficie al fine di mantenere il sistemadi equazioni invariato: cio, nell’ambito aerodinamico, trova la sua piunaturale applicazione nell’utilizzo di grandezze legate all’integraledella pressione e degli sforzi viscosi sul corpo, potendo percio consi-derare direttamente l’ottimizzazione rispetto alle forze di portanza eresistenza.

1.1 scopo e struttura del lavoro di tesi

Lo scopo del lavoro svolto e quello di estrarre l’andamento dellasensitivity avvalendosi dell’ausilio del programma OpenFOAM, e dianalizzarne l’andamento su corpi di varia geometria. L’obiettivo prin-cipale e quello di riuscire ad evidenziare, tramite il gradiente stesso,le parti di un oggetto piu sensibili ad una modifica per ottimizzarela funzione di costo scelta, e cio anche con l’intenzione di mettere inluce zone che difficilmente sarebbero prese in considerazione (cosache accade soprattutto quando si considerano corpi con geometriacomplessa). La scelta di sfruttare OpenFOAM come software di rife-rimento risiede nel fatto che, ad oggi, pochi programmi commercialiper le simulazioni CFD contengono al loro interno un solutore per leequazioni aggiunte nell’ambito dell’ottimizzazione di forma. Solamen-te l’ultima versione di ANSYS FLUENT ne contiene una, come pureil programma SU2 (Stanford University Unstructured), sviluppatodall’omonimo ateneo.

La tesi viene percio strutturata nei seguenti capitoli:

2. Teoria in questo capitolo vengono derivate in modo completole equazioni di Navier-Stokes aggiunte a partire dalle corrispon-denti equazioni dirette, insieme all’espressione della sensitivity.Viene inoltre ricavata l’approssimazione dello stesso sistemaaggiunto introdotta da Lohner.

3. Implementazione il software OpenFOAM viene presentato, in-

4

1.1 scopo e struttura del lavoro di tesi

sieme alla struttura generale di una simulazione eseguita conquesto programma. Viene poi introdotto il solutore per le equa-zioni aggiunte, nonche l’algoritmo di deformazione utilizzato nelcaso bidimensionale e il loop di ottimizzazione. Infine, vengonoanalizzati i criteri di arresto utilizzati nelle simulazioni.

4. Casi 2D il metodo e stato quindi applicato in questo capitolo adun profilo NACA 0012, sia con flusso laminare che turbolento, eper cui e stato sviluppato anche un algoritmo di deformazionedel corpo stesso, e cio al fine di validare il procedimento in casidapprima bidimensionali.

5. Casi 3D successivamente lo stesso metodo e stato applicato inun caso tridimensionale, senza pero la scrittura di un algoritmodi deformazione. Al fine di mettere in luce i due problemi princi-pali che sono presenti nelle equazioni aggiunte, ossia un’elevatadipendenza dalla griglia nonche dagli schemi numerici utilizzatiper la risoluzione del sistema stesso, viene proposto un confrontotra due mesh volumetriche differenti, una composta da poliedrie una puramente costituita da esaedri, e prendendo come riferi-mento la stessa geometria alare. Il metodo e infine applicato adun caso piu complesso di cui e riportato il procedimento, ossiaalla parte posteriore di una macchina di F1.

6. Conclusioni e sviluppi futuri

7. Appendici

5

2 TEOR IA

In questo capitolo vengono ricavate le equazioni aggiunte a partiredalle analoghe equazioni dirette laminari, e sfruttando il procedimentoche consiste nell’utilizzo dei moltiplicatori di Lagrange. Il trattamentodella turbolenza nel sistema aggiunto tramite il concetto della turbo-lenza congelata (frozen-turbulence) e quindi introdotto e discusso, esuccessivamente viene ricavata l’espressione della sensitivity. Infine,viene derivato il sistema aggiunto ottenuto con l’approssimazioneintrodotta da Lohner.

2.1 modello fisico-matematico

Il modello che viene utilizzato e quello delle equazioni di Navier-Stokes aggiunte continue, che sono ricavate a partire dalle analogheequazioni dirette stazionarie, viscose e incomprimibili nella loro formaanalitica, e solo successivamente alla loro derivazione vengono discre-tizzate; per contro viene anche ricordato l’approccio duale, ossia quellodelle equazioni di Navier-Stokes aggiunte discrete, ove il processo diderivazione delle stesse parte dalle analoghe equazioni dirette gia informa discreta; quest’ultimo approccio non verra seguito in questa tesi.Esistono diverse modalita per ricavare il set di equazioni in esame e ilprocedimento qui esposto sfrutta i moltiplicatori di Lagrange.

Innanzitutto occorre scegliere quale set di equazioni dirette utilizza-re come punto di partenza; per i casi analizzati, nonche per il regimedi moto che viene simulato, si e scelto di ricavare l’aggiunto delleequazioni di Navier-Stokes stazionarie, viscose e incomprimibili, chepossono essere espresse nel seguente modo [12]:

(v •∇)v = −∇p +∇ • (2ν D(v))

∇ • v = 0(1)

dove:

• D(v) =12(∇v + (∇v)T)

Le equazioni appena scritte sono valide a rigore solo in caso di regimedi moto laminare. Nel momento in cui si deve rappresentare inveceuna corrente turbolenta, tra i vari modelli a disposizione si e utilizza-to quello delle equazioni RANS (Reynolds-Averaged Navier-Stokes)unitamente all’ipotesi di Boussinesq, in cui tutta la modellazione dellaturbolenza e racchiusa all’interno di una viscosita spaziale non unifor-

7

teoria

me νt, che permette di riformulare le precedenti equazioni nel modoseguente: (v •∇)v = −∇p +∇ • (2(ν + νt(x))D(v))

∇ • v = 0(2)

dove e stata volutamente messa in evidenza la dipendenza dallecoordinate spaziali della viscosita turbolenta, la quale viene ricavataa partire dall’introduzione di uno a scelta tra i modelli di turbolenzache sono stati sviluppati.

La successiva derivazione delle equazioni aggiunte seguira lo svilup-po proposto da Othmer in [12]; viene introdotta quindi una funzionecosto J(β, v, p), dove:

• β e la deformazione in direzione normale al corpo;

• v e la velocita diretta;

• p e la pressione diretta;

che esprime, nella sua accezione piu generale, l’andamento di unacerta quantita sulla superficie del corpo e all’interno del dominio;l’obiettivo e la sua minimizzazione (o massimizzazione), unitamenteal rispetto del vincolo dettato dalle equazioni di Navier-Stokes stesse,che vengono espresse nella seguente forma residuale:(R1, R2, R3)T = (v •∇)v +∇p−∇ • (2νD(v))

R4 = −∇ • v(3)

dove viene definito R = (R1, R2, R3, R4)T, che e uguale a 0 per defini-

zione.

Viene percio introdotto il seguente funzionale:

L = J +∫

Ω(u, q) • R dΩ (4)

in cui i moltiplicatori di Lagrange u e q prendono il nome rispet-tivamente di velocita e pressione aggiunta. Calcolando dunque lavariazione totale di tale funzionale rispetto alle variabili β, v e p siotterra sia la sensitivity alla funzione costo J che le equazioni di Navier-Stokes aggiunte.

Nella scrittura della variazione di L, la viscosita ν verra trattatacome una costante, non considerando il termine che nascerebbe dallarelativa variazione. Questo passaggio e a rigori corretto solamente nelcaso di flusso laminare, in cui essa e effettivamente costante.

Qualora invece le equazioni aggiunte ricavate dovessero rappresen-tare l’evoluzione di un flusso turbolento, l’aver trascurato tale contri-buto portera all’approssimazione della cosiddetta “frozen-turbulence”,(in cui il campo di viscosita turbolenta e ritenuto lo stesso costante)ampiamente utilizzata in letteratura. Un approccio piu corretto al pro-blema imporrebbe di tener conto altresı di tale termine, e cio dovrebbe

8

2.2 equazioni aggiunte laminari

essere fatto tramite l’introduzione di un modello di turbolenza ancheper l’operatore aggiunto; questo approccio, tuttavia, esula dagli scopidella tesi e non verra trattato.

La variazione totale di L viene espressa nel modo seguente:

δL = δβL + δvL + δpL (5)

ed e quindi composta da un termine dipendente dalla deformazionedel corpo (δβL) e da due termini dipendenti dalla variazione dellevariabili dirette (δvL e δpL).

I moltiplicatori di Lagrange u e q sono ricavati in modo da far valerela condizione:

δvL + δpL = 0 (6)

Questa scelta viene fatta affinche, per il calcolo della sensitivity, non sianecessario conoscere la variazione del funzionale rispetto alle variabilidirette v e p, di fatto rendendo indipendente l’ottenimento del gra-diente dal calcolo di un numero di soluzioni del flusso pari al numerodi variabili di design considerate. Dalla risoluzione del sistema che siricava dalla condizione appena imposta, si otterra il sistema delle equa-zioni di Navier-Stokes aggiunte con relative condizioni al contorno,mentre la sensitivity rimarra semplicemente data dal termine:

δL ≡ δβL = δβ J + δβ

∫Ω(u, q) • R dΩ (7)

2.2 equazioni aggiunte laminari

Le equazioni aggiunte si ricavano, come gia detto, dall’imposizionedella condizione

δvL + δpL = 0 (8)

Le variazioni dei termini di (3) rispetto alle variabili di stato v e p siesprimono come:

δv(R1, R2, R3)T = (δv •∇)v + (v •∇)δv−∇ • (2νD(δv))

δvR4 = −∇ • δv

δp(R1, R2, R3)T = ∇δp

δpR4 = 0

(9)

e percio la condizione (8) viene tradotta nella seguente equazione:

δv J + δp J +∫

Ωu • ((δv •∇)v + (v •∇)δv−∇ • (2νD(δv))dΩ+

−∫

Ωq(∇ • δv)dΩ +

∫Ω

u • (∇δp)dΩ = 0(10)

A questo punto, decomponendo J in contributi di volume e disuperficie, ossia

J =∫

ΓJΓ dΓ +

∫Ω

JΩ dΩ (11)

9

teoria

e integrando per parti ciascun termine contenuto in (10), al fine dispostare le variazioni delle variabili di stato al di fuori degli operatoridifferenziali, le equazioni di Navier-Stokes aggiunte con gli opportunitermini di bordo possono essere ricavate.

Il sistema cosı ottenuto e:−∇u • v− (v •∇)u = −∇q +∇ • (2νD(u))− ∂JΩ

∂v

∇ • u =∂JΩ

∂p

(12)

con:

• D(u) =12(∇u + (∇u)T),

mentre le condizioni al contorno per le variabili aggiunte possonoessere dedotte annullando singolarmente i termini seguenti, partendodall’analisi delle condizioni al contorno che vengono imposte nellevarie zone del dominio per le variabili dirette:∫

Γ(n(u • v)+u(v • n) + 2νn • D(u)− qn +

∂JΓ

∂v) • δv dΓ+

−∫

Γ2νn • D(δv) • u dΓ = 0

∫Γ(u • n +

∂JΓ

∂p)δp dΓ = 0

(13)

Come si puo notare, il sistema appena ricavato e molto simile aquello delle equazioni di Navier-Stokes dirette, se non per tre aspettiprincipali:

• e lineare, infatti il trasporto della velocita aggiunta u avvienetramite la velocita diretta v e non tramite se stessa;

• la convezione dell’informazione avviene in direzione opposta,cioe upstream, rispetto al sistema diretto, come si puo vederedal segno − che precede i due termini del primo membro dellaprima equazione, e come esemplificato dalla figura 1. E evidenteda quest’immagine come l’orientamento del vettore velocitaaggiunta sia localmente in prevalenza opposto alla direzione delflusso;

• e presente un termine aggiuntivo in entrambe le equazioni, lega-to a JΩ, mentre il contributo superficiale della funzione obiettivo,JΓ, compare unicamente nelle condizioni al contorno.

10

2.3 equazioni aggiunte turbolente

Figura 1: Andamento tipico in direzione upstream della velocita aggiunta inun caso 2D

Un punto particolarmente importante e che, qualora J non contienetermini dipendenti dal volume, le equazioni aggiunte ne risultanoeffettivamente del tutto indipendenti, e quindi il sistema da risolvererimane invariato anche considerando funzioni di costo (purche solosuperficiali) diverse fra loro. In questo caso, che e quello preso inesame, l’unica dipendenza del sistema aggiunto rispetto a J rimane,come gia dimostrato in precedenza, nelle condizioni al contorno perle variabili aggiunte.

Per quanto riguarda i dati di bordo del sistema aggiunto per i casiconsiderati, che consistono in un corpo immerso in un dominio fluido,si puo trovare in letteratura che la condizione al contorno principaleche determina l’andamento delle equazioni e quella che si impone perla velocita aggiunta sul corpo da deformare, e che puo essere intesacome una condizione di “traspirazione” sul corpo stesso; per le condi-zioni sul farfield, invece, si ricava che sono sufficienti delle condizionidi Neumann sia per la velocita che per la pressione aggiunta (si vedaad esempio [13] per ulteriori dettagli).

2.3 equazioni aggiunte turbolente

Nel momento in cui si volesse calcolare l’aggiunto a partire da equazio-ni di Navier-Stokes che descrivono l’evoluzione di un flusso turbolento,sarebbe necessario tenere in considerazione, come gia detto, anche lavariazione spaziale della viscosita ν che in questo regime di moto ecomposta dalla somma di due contributi, ossia la viscosita molecolareνm e la viscosita turbolenta νt.

Cio viene fatto tramite l’approssimazione della frozen-turbulence, checonsta nell’ipotesi di considerare “congelato” il campo ν non uniformederivante dal solutore diretto, e di inserirlo cosı com’e nelle equazioni

11

teoria

aggiunte [13], che percio vengono modificate nel termine viscoso nelseguente modo:−2(D(u))v = −∇q + 2ν(∇ • D(u)) + 2∇ν • D(u)− ∂JΩ

∂v

∇ • u =∂JΩ

∂p

(14)

dove questa volta ν = νm + νt.

La validita dell’utilizzo di tale approssimazione e giustificata in [14],dove vengono proprio messi a confronto l’implementazione e l’erroreche si commette qualora si utilizzi una viscosita costante piuttosto cheun modello di turbolenza nelle equazioni aggiunte: cio che si ricavanei casi considerati e che l’informazione sulla sensitivity ottenuta con lafrozen-turbulence differisce di poco da quella ricavata tramite il modellocompleto. Piu nello specifico viene dimostrato come, con il modelloapprossimato, sia rispettato l’andamento qualitativo del gradienterelativamente a quello ottenuto con il sistema completo, e cio sia intermini di zone del corpo messe in evidenza, che di direzione in cuieffettuare la deformazione dell’oggetto; quello che non e conservatoesattamente e il valore puntuale della sensitivita, e il rapporto tragli ordini di grandezza della sensitivity stessa nei vari punti. Lavalidita dell’approssimazione e infine ulteriormente verificata con ilparallelo tra i valori del gradiente ottenuti con i due approcci e quelloricavato con le differenze finite, confronto che sottolinea la fondatezzadell’approssimazione tramite viscosita costante (in quanto i valoriottenuti nei 3 casi differiscono di poco fra loro). L’eventuale differenzatra i due andamenti della sensitivita risulta percio tradursi solamentein un’ottimizzazione “piu lenta” nel caso dell’approssimazione conviscosita congelata, seppur sempre alla stessa soluzione asintotica altermine dell’ottimizzazione, ma in questo caso e solamente necessarioeseguire piu cicli per ottenere lo stesso risultato del modello completo.

Nel caso con viscosita costante, tuttavia, i cicli di ottimizzazionesono molto piu veloci, essendo il sistema aggiunto composto da menoequazioni da risolvere.

Un possibile problema nasce nel momento in cui le equazioni direttevengono risolte tramite l’introduzione di un modello di turbolenzache sfrutta l’utilizzo delle wall functions, ossia di quelle funzioni chegestiscono automaticamente l’andamento delle quantita turbolentea parete (k, ε, νt ad esempio), permettendo di evitare un’eccessivarisoluzione della griglia in questa zona. La controparte “aggiunta” ditali wall functions e a tutt’oggi un argomento di ricerca, e percio lesoluzioni che possono essere attuate per ovviare a questa mancanzasono ad oggi sostanzialmente due:

• Low-Re turbulence models: e sicuramente l’approccio piu correttoe consiste nell’impiego di una mesh e relativo modello di turbo-lenza che risolvano totalmente lo strato limite, con un valore cioedi y+ attorno al corpo dell’ordine dell’unita; con questo metodoviene eliminato totalmente il suddetto problema, in quanto none necessaria l’introduzione delle wall functions neanche per leequazioni dirette. Tuttavia questo approccio porta, soprattutto

12

2.4 cost function

per casi tridimensionali, ad un elevato aumento del numero dicelle nel dominio, con conseguente innalzamento dei tempi dicalcolo della simulazione;

• L’altro criterio prevede di non utilizzare le wall functions perl’operatore aggiunto mentre per il diretto si. In questo caso sipossono ancora ricavare delle “buone” sensitivities utilizzandouna mesh che non sia troppo rifinita alla parete.

2.4 cost function

L’introduzione nel problema aggiunto di funzioni costo per l’ambitospecificatamente aeronautico e la derivazione del sistema in questocaso e presente ad esempio nell’articolo [15]: Soto e Lohner intro-ducono dapprima la forma generale della suddetta cost function chepermettera una riduzione della resistenza oppure un aumento dellaportanza sviluppato dal corpo. Come funzione obiettivo viene dunquepresa proprio la forza aerodinamica sviluppata che, nel caso in cui siabbia a che fare con fluidi viscosi, e composta dalla somma di duecontributi, uno legato alla pressione e uno agli sforzi tangenziali, epercio viene espressa come un integrale solamente sulla superficie delcorpo che subisce la deformazione, nella seguente forma:

J ≡ JΓ =∫

Γ(pn− 2νD(v)n)e f dΓ (15)

dove n e la normale rispetto al corpo assunta positiva uscente,e f la direzione della forza che si vuole ridurre e Γ la superficie sucui si effettua l’ottimizzazione. Analizzando dunque la variazionedi tale funzione obiettivo rispetto alle variabili v e p, e possibileesprimere esplicitamente i termini δv J e δp J contenuti nell’equazione(10) traducendoli, nello sviluppo del sistema, in condizioni al contornoda imporre per la velocita aggiunta sul corpo, che assumono percio laseguente forma:

u|Γ = −e f (16)

Come gia annunciato in precedenza, per cost function di questotipo, l’unica dipendenza del sistema aggiunto da J e nella condizioneappena ricavata.

2.5 calcolo della sensitivity

Il calcolo della sensitivity, sempre seguendo l’articolo [12], passaattraverso il calcolo della seguente quantita:

δL ≡ δβL = δβ J + δβ

∫Ω(u, q) • R dΩ (17)

che, dopo l’ottenimento delle equazioni aggiunte, e l’unico termine adessere rimasto diverso da zero.

Qualora la funzione di costo J abbia la forma analizzata in prece-denza, essa non ha dipendenza esplicita dalla deformazione, e quindi

13

teoria

il termine δβ J e identicamente nullo. Per quanto riguarda il secondocontributo esso risulta esprimibile, sempre nel caso in cui J ≡ JΓ,come un integrale esteso solo alla superficie del corpo di cui si attual’ottimizzazione (per dettagli sulla derivazione si veda nuovamenteanche [15]).

Viene dunque sfruttata la seguente proprieta: innanzitutto sonoripresi i residui delle equazioni di Navier-Stokes che, per definizione,sono identicamente nulli, ossia vale che:

R = 0 (18)

Da questa proprieta deriva direttamente la seguente:

δR = δβR + δvR + δpR = 0 (19)

ed e percio possibile esprimere δβR in funzione della variazionedelle variabili dirette, ossia:

δβR = −δvR− δpR (20)

L’equazione (17) puo quindi essere riformulata nel seguente modo:

δL = −∫

Ω(u, q) • δvR dΩ−

∫Ω(u, q) • δpR dΩ (21)

ed e facilmente verificabile come lo sviluppo di tali termini porti adun risultato analogo a quello contenuto nell’equazione (10), a menodei termini δv J e δp J. Entra a questo punto in gioco l’ipotesi fonda-mentale della dipendenza di J unicamente dalla superficie: solo inquesto caso, infatti, dall’integrazione per parti dei termini contenutiin (21), si ottiene esattamente il sistema di equazioni di Navier-Stokesaggiunte (12) (a meno dei contributi legati a JΩ), permettendo quindidi cancellare i termini relativi agli integrali di volume e lasciando lacost function (17) composta da contributi dipendenti unicamente dallasuperficie del corpo.

Quindi, facendo l’approssimazione di deformare il corpo solo nelladirezione localmente normale (come fatto ad esempio da Othmer), eomettendo tutti i passaggi, che sono comunque riportati negli articolicitati, e possibile riformulare la sensitivity nel seguente modo:

δL = −∫

Γν((n •∇)u) • ((n •∇)v)dΓ (22)

Come si puo vedere, essa dipende dal prodotto dei gradienti norma-li tra la velocita diretta e aggiunta a parete; questa quantita e perciodirettamente collegata al comportamento attorno al corpo della velo-cita aggiunta: dove appariranno, in relazione alla cost function scelta, isuoi gradienti piu elevati, allora in tali zone ci sara da aspettarsi unasensibilita maggiore alla modifica del corpo stesso.

Tale campo scalare viene dunque ricavato sul corpo considerato e,ove esso assume valori positivi, allora sara necessario deformare local-mente in direzione normale uscente per ottenere un miglioramentodella prestazione richiesta, e viceversa laddove avesse valori negativi.

Si sottolinea infine come la sensitivity sia una quantita che puo esserecalcolata puramente a livello di post-processing, una volta cioe che idue campi di moto siano gia stati ricavati tramite la risoluzione deirelativi sistemi di equazioni.

14

2.6 approssimazione alla lohner

2.6 approssimazione alla lohner

Sempre nell’articolo [15], e introdotta una variante nel sistema aggiun-to: viene infatti fatto notare come le equazioni ottenute siano statericavate a partire da un’approssimazione del second’ordine della va-riazione di R rispetto alle variabili di stato e, seguendo la derivazionepresentata nell’articolo, viene suggerito di trascurare il terzo terminecontenuto in (10) (che e appunto del second’ordine), allo scopo diimplementare un’approssimazione delle equazioni al prim’ordine.

Cosı facendo il sistema che si ricava e:−(v •∇)u = −∇q +∇ • (2νD(u))− ∂JΩ

∂v

∇ • u =∂JΩ

∂p

(23)

In questo caso il problema aggiunto, sebbene lineare, risulta averela medesima caratterizzazione convettiva delle equazioni di Navier-Stokes dirette, (ossia la convezione avviene tramite un solo termine enon tramite due come nel caso completo) e quindi viene sottolineatocome tale sistema possa venire risolto con pressoche le medesimetecniche numeriche utilizzate per l’analogo diretto.

Anticipando poi l’implementazione numerica del problema all’in-terno del software OpenFOAM, viene qui sottolineato come l’avertrascurato nelle equazioni il contributo:

−∇u • v (24)

porti anche ad una stabilizzazione nel comportamento stesso delsistema aggiunto.

15

3 IMPLEMENTAZ IONE

Una breve introduzione al software open source OpenFOAM e di segui-to proposta, con particolare enfasi sia sulla facilita di implementazionedelle equazioni alle derivate parziali che sulla struttura di una simu-lazione eseguita con questo programma. Viene quindi introdotto ilsolver per le equazioni aggiunte che e stato utilizzato, nonche il pro-gramma di deformazione del profilo nel caso bidimensionale e quelloper gestire il loop di ottimizzazione. Infine vengono presentati i criteriutilizzati per valutare la convergenza delle simulazioni.

3.1 presentazione del software openfoam

Il software utilizzato e la versione 2.1.1 di OpenFOAM. Questo pro-gramma sta acquisendo sempre piu popolarita, anche negli ambientiindustriali, soprattutto per la sua caratteristica di essere open source,ossia con codice sorgente disponibile, e modificabile da parte di chiun-que. E un codice a volumi finiti, molto versatile, il cui utilizzo spaziainfatti attraverso un’ampia gamma di problemi, passando dalle classi-che simulazioni RANS, LES e DNS, arrivando tutt’oggi a solver chegestiscono simulazioni chimiche, elettrostatiche e per l’analisi deglisforzi nei solidi.

Uno dei suoi punti di forza risiede sicuramente nella facilita discrittura delle equazioni alle derivate parziali che si vogliono risolvere,la cui sintassi ricorda molto quella che viene utilizzata quando si trattail problema dal punto di vista matematico.

In OpenFOAM, l’equazione

∂ρU∂t

+∇ • (φU)−∇ • (µ∇U) = −∇p (25)

viene infatti espressa nella seguente forma:

f vm :: ddt(rho, U)

+ f vm :: div(phi, U)

− f vm :: laplacian(mu, U)

==

− f vc :: grad(p)

(26)

Come si puo ben vedere, a parte i prefissi fvm ed fvc, ogni terminediscreto e ben riconducibile all’analogo continuo. Un altro punto afavore del programma risiede nel linguaggio di programmazione stes-so, ossia il C++, che con la sua innata caratteristica di object orientationpermette di definire dei contenitori, le cosiddette classi, che rendonomolto semplici eventuali modifiche e aggiunte al programma stesso.

17

implementazione

In un linguaggio orientato ad oggetti, particolare rilevanza assumela caratteristica del polimorfismo, estesamente utilizzata in OpenFOAM,con la quale e possibile definire oggetti appartenenti a classi diversecon lo stesso nome, potendo poi chiamare quello desiderato tramite ilriferimento alla sua classe di appartenenza.

In OpenFOAM l’esecuzione di una simulazione puo essere sostan-zialmente suddivisa in tre fasi:

• Una fase di pre-processing, in cui viene generata la mesh: per quan-to riguarda questo punto, il programma mette a disposizionenativamente due tool: il primo e blockMesh, che permette di gene-rare delle griglie strutturate, mentre il secondo e snappyHexMesh,con cui e possibile creare mesh non strutturate, semplicementedando in ingresso un file in formato stl contenente il corpo euna griglia di background;

• Una fase di risoluzione vera e propria in cui, chiamando unospecifico solutore, vengono risolte le equazioni ad esso associate,specificando le condizioni iniziali e al contorno, la mesh, glischemi numerici per i vari termini nonche alcuni parametri dicontrollo della simulazione;

• Una fase finale di post-processing, dove il programma ParaViewpermette di visualizzare i risultati ottenuti, consentendo di ana-lizzare i vari campi (vettoriali, scalari e tensoriali) all’interno deldominio.

La struttura di un caso e quella riportata nella figura 2:

Figura 2: Struttura di un caso in OpenFOAM

Come si puo vedere sono tre le cartelle principali, ossia 0, constant esystem. Esse contengono, rispettivamente:

18

3.1 presentazione del software openfoam

• 0: contiene le condizioni iniziali e al contorno per ogni variabilepresente nel problema; e disponibile all’interno del programmauna vasta scelta di condizioni gia pre-impostate (come quel-le di Dirichlet oppure di Neumann); nell’eventualita in cui lacondizione al contorno ricercata non fosse presente, e possibilescriverne una ex-novo e implementarla nel solver desiderato;

• constant: contiene al suo interno alcuni file per la definizione delcampo di moto, ossia transportProperties, turbulenceProperties eRASProperties, che definiscono il modello di turbolenza usato e ilvalore di alcuni coefficienti (come la viscosita), nonche la cartellapolyMesh che contiene la griglia vera e propria, rappresentatadai file points, cells, faces e boundary.

Qualora la mesh venisse generata tramite la utility blockMesh,allora occorrerebbe compilare il file blockMeshDict (come e statofatto per il caso 2D) nelle seguenti parti:

– Vertices: in questa prima parte vengono elencate le coordi-nate di tutti i nodi che compongono i vertici dei blocchicostituenti la griglia, nel formato (x y z);

– Blocks: in questa sezione vengono effettivamente composti ivari blocchi, specificandone i vertici in uno specifico ordine,il numero di celle in ogni direzione, nonche il relativograding; un tipico elemento di questa parte e espresso nellaforma:

hex (0 1 2 3 4 5 6 7) (10 20 1) simpleGrading (5 3 1) (27)

– Edges: questa sezione e una tra le piu importanti, poiche equi che ad esempio, nel caso 2D, sono elencate le coordinateche compongono il profilo, tramite la funzione polyLine;

– Patches: qui, infine, vengono costruite e assemblate tutte lepatch che compongono il dominio.

• system: e composta almeno da tre file, che sono quelli fondamen-tali, ossia:

– controlDict: contiene i parametri che governano la simula-zione, come ad esempio il numero di iterazioni oppure il∆t per le simulazioni instazionarie;

– fvSchemes: questo file e molto importante, poiche contienegli schemi numerici che sono utilizzati per ogni terminecontenuto nelle equazioni da risolvere: per le simulazioniaggiunte e stata riscontrata una grande sensibilita soprat-tutto dai termini convettivi e dai gradienti: l’utilizzo, perqueste due classi di termini, di alcuni schemi piuttostoche di altri, cambia radicalmente il comportamento delleequazioni;

– fvSolution: qui, infine, sono elencati i solutori e i criteri ditolleranza utilizzati per ogni singola variabile. Vengonoanche specificati altri parametri quali i relaxation factors,nonche eventuali criteri di arresto basati sui residui.

19

implementazione

3.2 solver per il sistema aggiunto

Il solutore utilizzato per risolvere le equazioni aggiunte prende spuntodal gia presente adjointShapeOptimizationFoam, solver implementato daC. Othmer per l’ottimizzazione topologica utilizzata, come gia detto,soprattutto nello studio dei flussi interni. Nel tool citato, ad ognipasso vengono risolte in cascata dapprima le equazioni dirette e poiquelle aggiunte, e viene fatto l’update del parametro α, che svolge ilruolo di “porosita” di ogni singola cella: esso funziona da parametrodi controllo della geometria e, dove piu tale valore e elevato, allorapiu in tale zona il fluido si comporta in realta come se fosse un solido,e quindi quella e la parte in cui l’ottimizzazione rispetto alla funzionecosto impostata sara piu efficace.

La risoluzione e attuata tramite un algoritmo di tipo SIMPLE (Semi-Implicit Method for Pressure-Linked Equations), che e quello utilizzatoanche nel solver simpleFoam per risolvere le canoniche equazioni diNavier-Stokes stazionarie e incomprimibili. L’approccio sfruttato daquesto solutore e quello segregato, ossia viene risolta dapprima l’equa-zione della quantita di moto e solo successivamente l’equazione dicontinuita. Piu nello specifico, ad ogni iterazione si risolve dapprimauna pseudo-equazione per il momentum, da cui si ricava una velocitaapprossimata; dopodiche viene risolta un’altra equazione, ricavatadalla divergenza della quantita di moto, e inserita nell’equazione dicontinuita al fine di determinare la pressione “corretta”, che e poiutilizzata per aggiornare altresı il valore di velocita precedentemen-te trovato (che solo ora rispetta il vincolo di continuita). Siccome ilsistema aggiunto sfrutta una modifica di SIMPLE, viene anch’essorisolto in modo segregato: questo e uno dei motivi principali in cuirisiede la cattiva convergenza delle equazioni aggiunte che, seppurlineari, venendo risolte per via iterativa risultano mal condiziona-te. Tuttavia, un solutore accoppiato all’interno di OpenFOAM none ad oggi presente, e quindi il problema non riesce ad essere eliminato.

Per adattare il codice alla shape optimization e stato percio necessarioeliminare il termine di porosita artificiale dal sistema aggiunto noncheprovvedere al calcolo della sensitivity, discretizzando l’equazione (22)nel modo seguente:

∆Li = −νAi((n •∇)u)i • ((n •∇)v)i (28)

Essa viene quindi calcolata sul corpo ed e un output del programma,quindi utilizzabile a posteriori, piu nello specifico nell’algoritmo dideformazione. Si e anche deciso di separare la risoluzione dei dueset di equazioni, e quindi verranno portate a convergenza dapprimale equazioni dirette, e successivamente i relativi campi di velocita epressione saranno sfruttati nel sistema aggiunto.

3.3 script python

Questo script, scritto tramite il linguaggio Python, provvede a defor-mare il profilo nel caso bidimensionale; esso riceve in input ad ognistep le coordinate del corpo nonche l’andamento della sensitivity, di

20

3.4 loop di ottimizzazione in bash

cui viene effettuato uno smoothing (per uniformarne l’andamento) eun rescaling rispetto al valore massimo. Quest’ultimo passaggio rendepiu agevole la scelta del fattore di “guadagno” per cui va moltiplicatoil gradiente ottenuto: siccome i campi di moto aggiunti sono (sia intermini di andamento e soprattutto di valori) molto dipendenti siadalla cost function scelta che dal campo di moto diretto, a parita di meshun’operazione di riscalatura del gradiente permette di limitare confacilita la massima deformazione ottenibile per ciascuna iterazione,che rimane collegata solamente alla scelta del guadagno. In questomodo e possibile evitare che, da uno step al successivo, il profilo sideformi eccessivamente e venga inficiata cosı la qualita dei risultatiottenuti. Dopodiche, viene calcolato l’andamento delle normali sulcorpo, ed ogni punto del profilo viene spostato in quella direzione,di una distanza pari al prodotto tra la sensitivity locale e il guadagnointrodotto in precedenza, costante per tutto il profilo. Accorgimentiparticolari sono presi al bordo d’uscita, siccome dopo la deformazionetale zona potrebbe risultare “aperta” oppure “intrecciata”. Come out-put sono infine create le nuove coordinate del profilo. Esse vengonoquindi importate da un secondo programma, che ha il compito digenerare in modo automatico il nuovo file di blockMesh.

3.4 loop di ottimizzazione in bash

Il procedimento sopra elencato e quindi automatizzato, nel caso 2D,tramite uno script in bash. Esso e, nei sistemi Unix e Unix-like, unashell (o terminale) che permette di eseguire sia comandi (copiare file,creare cartelle) che programmi, interfacciandosi direttamente con ilsistema operativo.

Quindi, e possibile scrivere un file che contenga la sequenza delleoperazioni da eseguire in un’iterazione del loop, e poi farle ancheripetere ciclicamente.

Le varie fasi che vengono effettuate ad ogni step sono:

1. Creazione cartelle per contenere il caso di OpenFOAM;

2. Risoluzione delle equazioni dirette tramite la chiamata al solversimpleFoam;

3. Passaggio dei file di velocita, pressione e viscosita turbolentaalla cartella contenente la struttura del caso per le equazioniaggiunte;

4. Risoluzione delle equazioni aggiunte con chiamata al relativosolver;

5. Passaggio dei file contenenti la sensitivity e le coordinate delprofilo allo script di deformazione;

6. Deformazione del profilo con lo script in Python introdotto inprecedenza;

7. Generazione del nuovo file di blockMesh e ripartenza dal punto1;

21

implementazione

3.5 simulazioni e criteri di arresto

Tutte le simulazioni dirette sono state svolte in modo canonico, ri-solvendo iterativamente le equazioni di Navier-Stokes stazionarie,viscose e incomprimibili fino al raggiungimento di una condizione diconvergenza per i coefficienti di forza e anche fino alla stabilizzazionedei residui delle varie componenti di velocita. Particolare attenzione,invece, e stata posta sullo studio del criterio di arresto da utilizzarenelle simulazioni aggiunte: tale concetto e, per questo set di equazioni,un argomento abbastanza complesso.

In primo luogo, occorre far presente che i residui del sistema ag-giunto sono risultati essere sempre piu elevati rispetto a quelli dellecorrispondenti equazioni dirette. Inoltre, va tenuto in considerazioneche essi potrebbero anche non essere l’indicatore piu credibile pergiudicare una simulazione aggiunta come “arrivata a convergenza”.

Per finire, e necessario mettere in evidenza che non e possibiledefinire a priori un set di schemi numerici che garantiscano la miglio-re convergenza delle equazioni aggiunte, in quanto esse dipendonomolto dal caso in esame, in termini di regime e tipologia del campodi moto, e in particolar modo dalla mesh utilizzata: un set di schemiottimi per una griglia composta da poliedri, ad esempio, potrebbeinvece portare alla divergenza se applicato ad una mesh composta daesaedri e viceversa.

Percio, per ogni simulazione, vengono messi a confronto tre indica-tori, ossia:

1. Residui delle componenti della velocita aggiunta;

2. Variazione della norma L2 della velocita aggiunta;

3. Variazione del valore massimo della velocita aggiunta.

Piu nel dettaglio ciascun criterio e definito come:

1. Sono i valori calcolati automaticamente da ogni solutore all’inter-no di OpenFOAM; essi sono adimensionali, e ottenuti dividendoad ogni iterazione la differenza tra il Right Hand Side e il LeftHand Side del sistema lineare per un fattore di normalizzazione;

2. Questo criterio analizza la variazione della velocita aggiuntanel dominio tra due generiche iterazioni, nel momento in cui iresidui si sono stabilizzati; la norma viene calcolata tramite unasommatoria su tutto il dominio dei moduli della differenza tra levelocita aggiunte nelle due iterazioni considerate, e adimensio-nalizzata rispetto alla sommatoria dei moduli su tutto il dominiocalcolata in una delle due iterazioni:

NL2 =

√∑ [(U1x −U2x)2 + (U1y −U2y)2 + (U1z −U2z)2]√

∑ (U21x + U2

1y + U21z)

(29)

dove con U1 si intende la velocita aggiunta ad una iterazione, econ U2 la velocita aggiunta N iterazioni dopo; e percio un valoreadimensionale, espresso in termini percentuali rispetto ad U1;

22

3.5 simulazioni e criteri di arresto

3. Per finire, viene analizzata la variazione del valore massimo divelocita aggiunta, sempre tra due generiche iterazioni, e unavolta che i residui si siano stabilizzati; tale valore e definitocome la differenza tra i valori massimi nel dominio della velocitaaggiunta, diviso per il massimo in una delle due iterazioni:

∆Umax =abs(max(|U1|)−max(|U2|))

abs(max(|U1|)(30)

Questo criterio, a differenza del precedente, e “locale”, in quantonon tiene in conto di tutti i valori della velocita nel dominio (chepotrebbero anche cambiare), ma solamente di quelli di picco;sebbene non globale, il valore massimo della velocita aggiunta sie verificato essere un buon indicatore per prevedere l’evoluzionedella simulazione: in runs divergenti, infatti, esso cresce senzalimiti, mentre in simulazioni stabili satura ad un certo valore eal limite oscilla.

23

4 CAS I 2D

In questo capitolo l’operatore aggiunto e stato applicato a casi bidi-mensionali. Sono state quindi effettuate due differenti simulazioni,una in regime di moto laminare e una in regime turbolento: questascelta e stata fatta dal momento che, per il secondo caso, e necessariointrodurre l’ulteriore approssimazione della gia citata frozen-turbulence,con conseguente modifica delle equazioni da risolvere. I casi sonostati svolti prendendo come riferimento in entrambe le simulazioni unprofilo NACA 0012 a 2 di incidenza e il setup, in termini di dimensio-ni della mesh, e il medesimo. Con questi studi si e voluto dimostrarela validita del metodo in entrambi i regimi di moto, insieme all’imple-mentazione di un vero e proprio ciclo di ottimizzazione comprendenteanche la deformazione del profilo.

4.1 cost function

Con lo scopo di ricercare un setup che fosse ingegneristicamenterilevante, si e percio deciso di implementare una cost function cheprovvede ad effettuare un aumento di portanza nonche una riduzionedi resistenza contemporaneamente. La quantita di cui si ricerca laminimizzazione e percio la seguente:

J = δ • CD − CL (31)

Sfruttando la linearita del problema aggiunto, e possibile dimostrareche la condizione al contorno da imporre per la velocita aggiunta sulprofilo con questa funzione costo “combinata” e pari alla sommavettoriale delle condizioni al contorno che si otterrebbero con le costfunction di “aumento di portanza” e “riduzione di resistenza” presesingolarmente (per un rimando al calcolo del tipo di condizione alcontorno da usare con funzioni obiettivo di questo tipo si guardi(16)). Facendo un esempio, supponendo che la corrente di fluidosia nella direzione [ 1 0 0 ], e che la portanza e la resistenza sianorispettivamente dirette in direzione [ 0 1 0 ] e [ 1 0 0 ], si dovrebbeimporre, come condizione al contorno per la velocita aggiunta sulcorpo:

• Condizione al contorno per la riduzione di resistenza: [−1 0 0 ]

• Condizione al contorno per l’aumento di portanza: [ 0 1 0 ]

• Condizione al contorno combinata:

δ [−1 0 0 ] + [ 0 1 0 ] = [−δ 1 0 ] (32)

dove δ e un fattore il cui scopo e quello di scalare opportuna-mente i due termini di portanza e resistenza, qualora avesseroordini di grandezza diversi. Un riferimento a riguardo puoessere trovato in [16]

25

casi 2d

4.2 caso 2d laminare

4.2.1 Mesh

Per questi casi, la griglia e stata costruita manualmente, generandolatramite la scrittura di un file blockMeshDict, e tramite l’esecuzione dellautility blockMesh contenuta in OpenFOAM. Essa e composta da seiblocchi, tre al di sotto del profilo e tre al di sopra, e ciascuno di essi e asua volta suddiviso in esaedri, che sono gli unici elementi componentila mesh, e a cui viene dato nelle varie direzioni il grading opportunoper rappresentare lo strato limite e la scia dietro al profilo. SiccomeOpenFOAM risolve equazioni tridimensionali, per rappresentare labidimensionalita del problema viene messa solamente una cella nellaterza dimensione, ossia quella perpendicolare al corpo. La mesh “a C”strutturata cosı ottenuta e composta totalmente da 292000 esaedri, e siestende per 10 corde davanti al profilo e per 40 dietro di esso.

Nelle figure 3 e 4 vengono riportati rispettivamente il dominio nelcaso bidimensionale e uno zoom rappresentante la griglia attorno albordo d’attacco del profilo:

Figura 3: Dominio nel caso bidimensionale

Figura 4: Zoom della mesh attorno al bordo d’attacco del profilo NACA 0012

26


La griglia cosı ottenuta viene suddivisa in 4 patch, che sono asse-gnate nel modo seguente:

• Inlet: e composta dalla C, nonche dalla parte inferiore della mesh,al fine di simulare un angolo di attacco maggiore di 0;

• Outlet: e composta dalla zona superiore della mesh e dalla patchverticale a destra;

• FrontAndBack: sono le due facce parallele tra loro e perpendico-lari all’asse z;

• Airfoil: e il profilo vero e proprio.

4.2.2 Condizioni al contorno

Le condizioni al contorno che vengono imposte per le variabili direttesono le seguenti:

• Velocita:

– Inlet: valore fissato, pari al modulo della velocita, in dire-zione x;

– Outlet: condizione di Neumann;

– Airfoil: condizione di adesione, ossia velocita nulla in tuttee tre le componenti;

• Pressione:

– Inlet: condizione di Neumann;

– Outlet: valore fissato a 0;

– Airfoil: condizione di Neumann;

Mentre per le variabili aggiunte si impone:

• Velocita aggiunta:



– Airfoil: la suddetta condizione, che puo essere intesa comeuna condizione di “traspirazione”;

• Pressione aggiunta:



– Airfoil: condizione di Neumann;

4.2.3 Simulazione diretta

La simulazione diretta e stata svolta con il solutore simpleFoam, cherisolve le equazioni di Navier- Stokes stazionarie e incomprimibili. Ilprofilo considerato e posto, come gia detto, a 2 di incidenza, ha cordadi 1 metro ed e investito da una corrente di aria con modulo di velocitapari a 0.01 m/s, in modo da avere un numero di Reynolds basato sulla

27

casi 2d

corda pari a 680, con un flusso totalmente in regime laminare.

Sfruttando la utility forceCoeffs di OpenFOAM e possibile inoltreottenere in output dal solutore i valori di CL e CD del profilo, e percioritenere che la simulazione sia a convergenza nel momento in cui talicoefficienti hanno raggiunto un valore stazionario.

Questo e l’andamento dei coefficienti di forza al primo step dellasimulazione, ossia con il profilo indeformato, a titolo esemplificativodella convergenza raggiunta ad ogni passo. Nei grafici seguenti,(figure 5 e 6) i valori di CL e CD sono stati adimensionalizzati rispettoal valore a convergenza:

1000 2000 3000 4000 5000 6000Iterazioni

0.95

1.00

1.05

1.10

1.15

1.20

1.25

C L

Figura 5: Coefficiente di portanza adimensionale nel caso 2D laminare

1000 2000 3000 4000 5000 6000Iterazioni

0.98

0.99

1.00

1.01

1.02

1.03

1.04

1.05

C D

Figura 6: Coefficiente di resistenza adimensionale nel caso 2D laminare

28


Le immagini seguenti riportano invece l’andamento dei campi dellavelocita (figura 7) e del coefficiente di pressione (figura 8), sempre allostep iniziale:

Figura 7: Campo di velocita attorno al profilo NACA 0012 nel caso 2Dlaminare

Figura 8: Cp attorno al profilo NACA 0012 nel caso 2D laminare

Il profilo viene investito da una corrente inclinata di 2 dal bassoverso l’alto e percio si puo notare come, sia nel valore della velocitache nel Cp, ci siano dei picchi rispettivamente positivi e negativi suldorso del profilo.

29

casi 2d

4.2.4 Simulazione aggiunta

Il campo di velocita ottenuto viene dunque inserito, come gia det-to, all’interno delle equazioni aggiunte e comincia percio la lororisoluzione.

Inizia ora un’analisi sui criteri di arresto citati all’inizio del capitolo,i quali possono essere utilizzati come indicatori della convergenzadella simulazione.

Il primo che viene esposto e quello relativo ai residui, come esem-plificato dalla figura 9:

0.0001

0.001

0.01

0.1

1

0 1000 2000 3000 4000 5000 6000 7000 8000

Res

idu

i

Iterazioni

"Uax""Uay"

Figura 9: Residui della simulazione aggiunta nel caso 2D laminare

Come si puo vedere, i valori si stabilizzano attorno ad un ordinedi grandezza di 10−3 per le componenti della velocita aggiunta gia apartire dall’iterazione 2000 in avanti: questi non sono valori “piccoli”se paragonati alle simulazioni canoniche, ma il punto sicuramentepiu importante da mettere in luce e il fatto che la simulazione cosıottenuta non diverge. Questo aspetto e quello di maggior rilievo perle simulazioni aggiunte, unitamente al fatto che anche la sensitivityraggiunge un andamento costante al variare delle iterazioni nel mo-mento in cui i residui si stabilizzano.

Inoltre, vengono riportati i valori di N2 e ∆Umax: il calcolo di questidue indicatori e stato effettuato una volta che i residui si fossero stabi-lizzati, e piu precisamente tra le iterazioni 7000 e 8000.

I valori cosı ottenuti sono:

• N2 = 0.027 (2.7%)

• ∆Umax = 0.00017 (0.017%)

Percio, da questi due criteri si puo avere un’idea del raggiungi-mento della convergenza. Gli indicatori cosı ottenuti sono ritenutisoddisfacenti: quindi la simulazione viene fermata all’iterazione 8000,

30


e il relativo gradiente viene calcolato, sfruttando la conoscenza dientrambi i campi di velocita. Si passa infine al punto centrale delprocedimento, ossia la parte di deformazione del corpo.

I tre indicatori, in termini di valori e andamenti, seppur riportatiper un solo step della simulazione aggiunta, sono stati ottenuti adogni iterazione del loop.

4.2.5 Deformazione del profilo

La sensitivity e ottenuta come output del solutore aggiunto e vienepassata al programma di deformazione scritto in Python. Qui neviene effettuato uno smoothing (per eliminare eventuali forti gradientipresenti nel suo andamento) e un rescaling rispetto al valore massimo,in modo che sia sempre compresa tra i valori −1 e 1 (per il motivo giaspiegato nel capitolo dell’implementazione).

Per il caso in esame l’andamento della sensitivity ottenuto e quelloriportato in figura 10:

1.0 0.5 0.0 0.5 1.0Corda del profilo

1.0

0.5

0.0

0.5

1.0

Sens

itivi

ty

Figura 10: Andamento della sensitivity sul profilo nel caso 2D laminare

dove, scorrendo l’asse delle ascisse da −1 a +1 si puo analizzarel’andamento della sensitivity a partire dal bordo d’uscita e ruotando insenso orario attorno al profilo. Le zone piu sensibili ad una modificasono proprio quelle vicine al bordo d’uscita (dove sono concentratii valori in assoluto piu elevati) e al bordo d’attacco, mentre, come sipuo vedere, il resto del profilo risulta abbastanza insensibile. Unaspiegazione al fatto che il gradiente abbia i suoi valori piu elevati albordo d’uscita del profilo e da ricercarsi nel comportamento del flussodiretto attorno al corpo: come si puo vedere dalla figura 7, in questocaso laminare e presente una scia molto consistente, e la sensitivityevidenzia come zona piu propensa ad una modifica proprio questa,dove il flusso e meno attaccato al corpo.

31

casi 2d

Un punto particolarmente importante risiede nella scelta dell’impo-sizione di uno o piu vincoli geometrici sul corpo una volta deformato,al fine di rendere possibile un confronto significativo tra step diversi.Il criterio scelto prevede di vincolare la corda del profilo affinche, adogni iterazione, sia sempre unitaria (e pari cioe a quella del profiloindeformato). Percio, ad ogni step, sia le ascisse che le ordinate del-l’oggetto deformato vengono divise per il valore della “nuova” corda,riscalandone le dimensioni ogni volta.

Al termine del ciclo di ottimizzazione il profilo risulta essere defor-mato come riportato in figura 11:

0.0 0.2 0.4 0.6 0.8 1.0

0.10

0.05

0.00

0.05

0.10

Profilo indeformatoProfilo deformato

Figura 11: Deformazione del profilo nel caso 2D laminare

Con uno zoom sul bordo d’attacco e sul bordo d’uscita (rispettiva-mente figure 12 e 13):

0.05 0.00 0.05 0.10 0.15 0.200.10

0.05

0.00

0.05

0.10Profilo indeformatoProfilo deformato

Figura 12: Deformazione del profilo, zoom al bordo d’attacco nel caso 2Dlaminare

32


0.80 0.85 0.90 0.95 1.00 1.050.10

0.05

0.00

0.05


Figura 13: Deformazione del profilo, zoom al bordo d’uscita nel caso 2Dlaminare

Come si puo vedere, e proprio la zona del bordo d’uscita a subire lamaggiore deformazione, diventando addirittura a doppia curvatura.

4.2.6 Loop

Come gia accennato, i passaggi analizzati in precedenza vengono ripe-tuti ciclicamente, e cosı il profilo si deforma sempre di piu seguendol’informazione ricavata dalla sensitivity. Di seguito, vengono presentatii grafici che mostrano la convergenza dei coefficienti nelle equazionidirette per il CL e per il CD (figure 14 e 15), a titolo esemplificativosolo al primo e all’ultimo step:

2500 3000 3500 4000 4500 5000 5500 6000 6500 7000Iterazioni

0.99

1.00

1.01

1.02

1.03

1.04

C L

Step_0Step_22

Figura 14: Aumento del CL durante l’ottimizzazione nel caso 2D laminare

33

casi 2d

2500 3000 3500 4000 4500 5000 5500 6000 6500 7000Iterazioni

0.990

0.995

1.000

1.005

C D

Step_0Step_22

Figura 15: Riduzione del CD durante l’ottimizzazione nel caso 2D laminare

Come si puo facilmente constatare, relativamente allo step 0, ossiaquello iniziale, si ottiene un aumento nel coefficiente di portanza dipiu del 3% e una riduzione del coefficiente di resistenza di piu dell’1%.

La simulazione viene bloccata nel momento in cui un’ulteriore stepdi modifica del profilo non genera piu l’effetto ricercato (aumentodi portanza e contemporaneamente riduzione di resistenza), ancheriducendo di alcuni ordini di grandezza il fattore di guadagno.

34


Di seguito viene infine riportato l’intero andamento dei coefficientidi forza attraverso le iterazioni (figure 16 e 17), ciascuno adimensiona-lizzato rispetto al valore del profilo indeformato:

0 5 10 15 20 25Iterazioni

1.000

1.005

1.010

1.015

1.020

1.025

1.030

1.035C L

Figura 16: Evoluzione dell’aumento del CL nel caso 2D laminare


0.986

0.988

0.990

0.992

0.994

0.996

0.998

1.000

C D

Figura 17: Evoluzione della riduzione del CD nel caso 2D laminare

35

casi 2d

In aggiunta a cio, se si combinano i due coefficienti, si puo anchenotare come venga aumentata ulteriormente l’efficienza del profilostesso, infatti si ha:


1.00

1.01

1.02

1.03

1.04

1.05

Effic

ienz

a

Figura 18: Evoluzione dell’aumento dell’efficienza nel caso 2D laminare

Il valore dell’efficienza e normalizzato rispetto a quello dello stepiniziale. Si puo percio vedere dalla figura 18 che, alla fine dell’ottimiz-zazione, essa risulta aumentata di piu del 4%.

Viene riportato inoltre l’andamento del coefficiente di pressione,ossia del Cp, tra la prima e l’ultima iterazione (figure 19 e 20), perdimostrare come venga ad essere effettivamente modificato il com-portamento del flusso, nel profilo deformato rispetto a quello iniziale.Come si puo notare, l’andamento del Cp nel grafico relativo al dorso, acirca x = 0.1, e inusuale: cio e dovuto al tipo di algoritmo di deforma-zione che e stato implementato, e alla presenza di piccole irregolaritache vengono a generarsi sul corpo stesso.

36


Figura 19: Confronto del Cp sul dorso del profilo nel caso 2D laminare

Figura 20: Confronto del Cp sul ventre del profilo nel caso 2D laminare

Infine, come ulteriore validazione del caso, viene riportato con lafigura 21 altresı l’andamento della funzione obiettivo, il cui valore,(adimensionalizzato rispetto a quello di partenza) giustamente, siriduce da un’iterazione alla successiva, siccome ne viene richiesta laminimizzazione:


0.55

0.60

0.65

0.70

0.75

0.80

0.85

0.90

0.95

1.00

Funz

ione

obi

ettiv

o

Figura 21: Riduzione della funzione obiettivo nel caso 2D laminare

37

casi 2d

4.3 caso 2d turbolento

In questa seconda simulazione viene ricercata un’ottimizzazione delprofilo nel caso in cui il regime di moto sia completamente turbolento.Cio porta all’introduzione dell’approssimazione della frozen-turbulencenelle equazioni aggiunte, e quindi anche il campo di viscosita turbo-lenta, oltre a quelli di velocita e pressione, viene passato e sfruttatodal sistema aggiunto.

Nel caso turbolento il parametro δ, introdotto in precedenza perpesare i contributi di portanza e resistenza nella cost function, haassunto un valore pari a 10, al fine di tenere in considerazione i diversiordini di grandezza delle due forze.

4.3.1 Mesh

La griglia utilizzata per il caso turbolento e, nelle dimensioni, identicaa quella del caso laminare. Essa viene infatti generata con il mede-simo file blockMeshDict in termini di vertices, blocks e patches, mentrele uniche differenze tra le due griglie risiedono nel numero di celle,che passa da 292000 a 528000, e nel grading (ossia nel rapporto tra ledimensioni di celle adiacenti) in direzione perpendicolare al profilo, inmodo da ottenere un valore di y+ ragionevole in relazione al modellodi turbolenza utilizzato, attorno a tutto il corpo.

4.3.2 Modello di turbolenza

Il modello di turbolenza che si e utilizzato e quello di Spalart-Allmaras:esso e un modello “ad un’equazione” completo, che permette di tro-vare direttamente l’andamento della viscosita turbolenta νt tramite larisoluzione di un’unica equazione differenziale in cui essa e la solaincognita presente. L’equazione, che non deriva da alcuna legge fisica,viene modellata tramite l’introduzione di 3 funzioni e 8 coefficien-ti empirici, e il modello trova la sua piu naturale applicazione nelcontesto aeronautico, per il quale e stato sviluppato.


Per quanto riguarda le variabili v e p, come pure le variabili aggiunte,le condizioni al contorno sono le stesse utilizzate nel caso laminare senon per il valore del modulo della velocita diretta all’inlet del dominio.Per quanto riguarda invece la viscosita turbolenta, il suo valore efissato anch’esso all’inlet, con l’imposizione di un valore ricavato conla seguente legge (per ulteriori dettagli si puo consultare [17]):

νt =

√32

v I l (33)

dove:

• v e la velocita all’infinito a monte;

38


• I e l’intensita della turbolenza, assunta pari all’1%;

• l e la lunghezza della scala turbolenta, presa pari al 5% dellacorda del profilo.

Lo stesso valore viene utilizzato per inizializzare il campo di vi-scosita turbolenta all’interno del dominio, mentre sul profilo vieneimposto un valore molto piccolo (10−20), siccome imponendo il valorenullo vero e proprio si genererebbero dei problemi numerici.


La simulazione diretta e molto simile a quella svolta per il casolaminare, se non per due punti, ossia:

• il valore della velocita all’infinito a monte viene imposto paria 30 m/s, al fine di avere un numero di Reynolds basato sullacorda del profilo pari a 2 • 106;

• viene simulato un flusso turbolento, utilizzando il modello diSpalart-Allmaras appena introdotto.

L’angolo di attacco viene invece sempre mantenuto pari a 2 comenel caso precedente e il valore della y+ attorno al profilo ha un valoremedio di 35.

Anche in questo caso le simulazioni sono state svolte fino a conver-genza guardando, oltre che l’andamento dei residui, la stabilizzazionedei coefficienti di forza.

Viene dunque riportato l’andamento del CL (figura 22) e del CD

(figura 23) rispetto al numero di iterazioni necessarie per la lorostabilizzazione, e adimensionalizzati rispetto al valore di convergenza:

1000 2000 3000 4000 5000 6000Iterazioni

0.86

0.88

0.90

0.92

0.94

0.96

0.98

1.00

1.02

C L

Figura 22: Coefficiente di portanza adimensionale nel caso 2D turbolento

39

casi 2d

1000 2000 3000 4000 5000 6000Iterazioni

0.995

1.000

1.005

1.010

1.015

1.020

1.025

C D

Figura 23: Coefficiente di resistenza adimensionale nel caso 2D turbolento

Ugualmente per questa seconda simulazione vengono riportati degliscreenshots sugli andamenti di velocita (figura 24) e pressione (figura25):

Figura 24: Campo di velocita attorno al profilo NACA 0012 nel caso 2Dturbolento

Figura 25: Cp attorno al profilo NACA 0012 nel caso 2D turbolento

40


Il profilo viene sempre investito da una corrente deflessa di 2 dalbasso ma, in questo caso ad alto numero di Reynolds, si puo vederecome la scia sia molto piu ridotta rispetto alla simulazione laminare.


La simulazione e stata quindi condotta fino a che i residui si fosserostabilizzati. Come si puo vedere in figura 26, il numero di iterazioninecessarie alla stabilizzazione della simulazione e molto maggiorerispetto al caso precedente, e cio puo essere imputato alla differen-za della mesh nonche del solutore, che contiene ora il termine divariazione spaziale della viscosita.

1e-06

1e-05

0.0001

0.001

0.01

0.1

1

0 2000 4000 6000 8000 10000 12000 14000 16000 18000 20000

Res

idu

i

Iterazioni

"Uax""Uay"

Figura 26: Residui della simulazione aggiunta nel caso 2D turbolento

Questa volta i residui si assestano a dei valori piu bassi rispetto alcaso laminare, attorno all’ordine di grandezza di 10−5.

I valori di NL2 e ∆Umax raggiunti a questo numero di iterazioni sonopiu piccoli rispetto al caso precedente:

• NL2 = 0.0082 (0.82%)

• ∆Umax = 0.0056 (0.56%)

E tali indicatori sono stati calcolati addirittura tra un intervallo mag-giore rispetto al caso laminare, ossia tra l’iterazione 17500 e l’iterazione20000, e dunque , con questi valori, la simulazione aggiunta vieneassunta come “arrivata a convergenza”.

Ugualmente in questo caso viene fatto presente che e stato riportatosolamente l’andamento dei residui di uno step siccome tutti gli altrisono molto simili, in quanto le equazioni e la geometria del profilodifferiscono di poco da un’iterazione alla successiva.

41

casi 2d

4.3.6 Deformazione del profilo

Come nella simulazione precedente, viene a questo punto acquisito enormalizzato l’andamento della sensitivity, che viene percio riportatonella figura 27:

1.0 0.5 0.0 0.5 1.0Corda del profilo

1.0

0.8

0.6

0.4

0.2

0.0

0.2

0.4

Sens

itivi

ty

Figura 27: Andamento della sensitivity sul profilo nel caso 2D turbolento

Come si puo vedere, a differenza del caso laminare, la zona piusensibile e sicuramente quella del bordo d’attacco, e cio puo dipenderedal diverso andamento sia del campo diretto che di quello aggiunto:nella fattispecie, questa volta la scia al bordo d’uscita e meno elevatarispetto al caso laminare, infatti il flusso e quasi del tutto attaccato,e quindi tale zona risulta meno propensa ad una modifica rispettoinvece al bordo d’attacco.

Il profilo viene dunque deformato seguendo l’indicazione del sud-detto gradiente e, anche in questo caso, viene imposto il vincologeometrico relativo alla corda, mantenuta sempre pari a quella delprofilo iniziale. Ecco quindi cio che si ottiene dopo 33 iterazioni (figura28):

0.0 0.2 0.4 0.6 0.8 1.0

0.10

0.05

0.00

0.05

0.10

Profilo indeformatoProfilo deformato

Figura 28: Deformazione del profilo nel caso 2D turbolento

42


Con, ancora una volta, uno zoom su bordo d’attacco (figura 29) ebordo d’uscita (figura 30), ossia sulle zone piu sensibili:

0.05 0.00 0.05 0.10 0.15 0.200.10

0.05

0.00

0.05


Figura 29: Deformazione del profilo, zoom al bordo d’attacco nel caso 2Dturbolento

0.80 0.85 0.90 0.95 1.00 1.050.10

0.05

0.00

0.05


Figura 30: Deformazione del profilo, zoom al bordo d’uscita nel caso 2Dturbolento

43

casi 2d

4.3.7 Loop

Infine, anche per questo caso viene riportata la storia della convergen-za dei coefficienti di forza nella prima e ultima iterazione, al fine dievidenziare la convergenza delle equazioni dirette che e stata raggiun-ta ad ogni passo:

1000 2000 3000 4000 5000 6000Iterazioni

0.8

0.9

1.0

1.1

1.2

1.3

1.4

1.5

1.6

C L

Step_0Step_33

Figura 31: Aumento del CL durante l’ottimizzazione nel caso 2D turbolento

1000 2000 3000 4000 5000 6000Iterazioni

0.97

0.98

0.99

1.00

1.01

1.02

1.03

1.04

1.05

C D

Step_0Step_33

Figura 32: Riduzione del CD durante l’ottimizzazione nel caso 2D turbolento

44


Come si puo constatare dalle figure 31 e 32, in questo caso il coef-ficiente di portanza aumenta di piu del 50% rispetto al suo valoreiniziale, mentre la riduzione del coefficiente di resistenza si attestaattorno al 2%. In piu, combinando i due valori, si ha che l’efficienzadel profilo e aumentata del 58% circa al termine dell’ottimizzazione.Anche qui il loop e stato fermato nel momento in cui, da uno stepal successivo, solo una delle due richieste in termini di aumento diportanza o riduzione di resistenza non veniva piu rispettata.

Di seguito viene quindi riportato l’andamento sia dei singoli CL

(figura 33) e CD (figura 34) che dell’efficienza (figura 35), durantel’ottimizzazione (tutti i valori sono adimensionalizzati rispetto a quelliiniziali):

0 5 10 15 20 25 30 35Iterazioni

1.0

1.1

1.2

1.3

1.4

1.5

1.6

C L

Figura 33: Evoluzione dell’aumento del CL durante l’ottimizzazione nel caso2D turbolento

0 5 10 15 20 25 30 35Iterazioni

0.975

0.980

0.985

0.990

0.995

1.000

C D

Figura 34: Evoluzione della riduzione del CD durante l’ottimizzazione nelcaso 2D turbolento

45

casi 2d

0 5 10 15 20 25 30 35Iterazioni

1.0

1.1

1.2

1.3

1.4

1.5

1.6

Effic

ienz

a

Figura 35: Evoluzione dell’aumento dell’efficienza durante l’ottimizzazionenel caso 2D turbolento

E anche qui sono riportati gli andamenti del Cp tra la prima e l’ul-tima iterazione (nelle figure 36 e 37), per evidenziare ulteriormentecome cambia l’andamento della pressione nel profilo deformato rispet-to a quello iniziale. In questo caso e molto evidente come l’aumento diportanza ottenuto sia legato sostanzialmente ad un valore di Cp mas-simo piu basso sul dorso del profilo, al termine dell’ottimizzazione:

Figura 36: Confronto del Cp sul dorso del profilo nel caso 2D turbolento

46


Figura 37: Confronto del Cp sul ventre del profilo nel caso 2D turbolento

In questi grafici l’irregolarita nell’andamento del Cp, presente so-prattutto sul dorso del profilo a x = 0.26, e legata ancora una voltaall’algoritmo di deformazione usato.

Per finire, anche in queste simulazioni e stata calcolata la funzionecosto per validare il metodo; in questo caso, come in quello laminare,il suo valore si riduce nel susseguirsi delle iterazioni:

0 5 10 15 20 25 30 35Iterazioni

25

20

15

10

5

0

Funz

ione

obi

ettiv

o

Figura 38: Riduzione della funzione obiettivo nel caso 2D turbolento

L’andamento della funzione obiettivo (come quello riportato infigura 38) e un buon criterio per validare i due casi trattati, in quantoil suo comportamento evidenzia se, effettivamente, l’ottimizzazionesta procedendo nel modo corretto oppure no. Non e stato invecepossibile effettuare ulteriori confronti con casi gia esistenti, a causadelle diversita in termini di geometria, settings, condizioni al contorno,e soprattutto parametrizzazione della forma.

47

casi 2d

Con questi due casi e stato possibile verificare l’efficacia dell’opera-tore aggiunto nell’ambito bidimensionale sotto le ipotesi sia di flussolaminare che turbolento, al fine di verificare separatamente l’ipotesidella frozen-turbulence, e unitamente alla scrittura di un algoritmo dideformazione per il profilo considerato. Sono state inoltre evidenziatele differenze tra il comportamento delle equazioni aggiunte nelle duesimulazioni, soprattutto in termini di andamento dei residui dellecomponenti di velocita. Tutto il procedimento di ottimizzazione estato inoltre automatizzato al fine di avere un tool versatile. A questopunto, dopo aver applicato il metodo in quest’ambito, si e percio pron-ti ad applicare il procedimento di risoluzione del sistema aggiunto nelpiu complesso caso 3D.

48

5 CAS I 3D

In questo capitolo vengono presentati i risultati relativi ai casi tridi-mensionali: lo studio qui effettuato ha permesso di ottenere un buonlivello di convergenza delle equazioni aggiunte, e con l’analisi svolta irisultati ricavati possono essere presi come punto di partenza per esse-re applicati a casi piu complicati. La stessa geometria e stata utilizzatacome base per due mesh volumetriche differenti, ossia una trim-meshottenuta con la utility snappyHexMesh e una griglia completamentecomposta da poliedri, generata con il programma ANSA. In questomodo, si e potuto studiare il comportamento dell’operatore aggiuntorelativamente alla sua dipendenza dalla mesh, e anche rispetto aglischemi numerici utilizzati.

Verranno ora trattate brevemente le tipologie di generazione di unagriglia per il caso tridimensionale; i metodi utilizzati sono due:

5.1 snappyhexmesh

Con questo tool incluso in OpenFOAM e possibile generare una trim-mesh 3D a partire da un file in formato stl (ossia STereoLithography)contenente la geometria del corpo, e da una griglia di background:la mesh che viene costruita in questo modo e composta puramenteda esaedri, se non a ridosso dell’oggetto, dove sono presenti deglielementi di forma tetraedrica per rappresentare il contorno del corpostesso. Il procedimento per la sua creazione e composto da 3 passaggi:

1. Castellated mesh: in questa prima fase viene generata una grigliaabbastanza grossolana, con delle zone di rifinitura della stessaattorno al corpo sulla base dei parametri specificati nel sottodizionario castellatedMeshControls; l’effetto cosı ottenuto e quellodi una griglia “a gradini” che verra in seguito rifinita, e succes-sivamente vengono anche rimosse le celle all’interno del corpostesso in quanto inutili;

2. Snapping process: a questo punto viene presa la mesh della faseprecedente e ogni vertice delle celle costituenti la superficie delcorpo “a gradini” e mosso al fine di ottenere una griglia il piupossibile simile al file fornito in input all’inizio, nel rispetto deiparametri di qualita imposti;

3. Layers: in quest’ultimo passaggio viene estruso, a partire dallasuperficie, un certo numero di prism layers, al fine di simulare lostrato limite attorno al corpo;

Le griglie generate in questo modo sono quasi sempre di “qualita”sufficiente per le simulazioni; tuttavia, bisogna far presente che edifficile far crescere uno strato di prismi uniforme attorno a geometrie

49

casi 3d

complesse in quanto, soprattutto attorno ad angoli affilati, i layerscreati nell’ultimo step tendono a collassare e cio puo effettivamen-te creare problemi nella risoluzione delle equazioni. Per ovviare aquesto problema e necessario agire sui parametri contenuti nel sottodizionario addLayersControls, nonche rilassare i parametri di qualitain meshQualityControls, a scapito della “bonta” della mesh ottenuta: equesta, infatti, una delle cause principali del comportamento sopracitato, ossia il fatto che snappyHexMesh rimuove gli strati di prismidalle zone che non rispettano le caratteristiche specificate.

5.2 ansa

Il programma industriale che viene utilizzato e prodotto dall’aziendaBETA CAE Systems, e permette di partire da un file CAD di uncorpo e di generare sia la mesh superficiale che quella volumetrica;al fine di creare una griglia di sufficiente qualita per OpenFOAM, estato necessario utilizzare come elementi costituenti la mesh stessa itetraedri oppure i poliedri. C’e inoltre da citare la totale integrazionetra i due programmi siccome e possibile, una volta completato ilprocesso di generazione della griglia volumetrica, creare ed esportareanche le altre cartelle che compongono un caso di OpenFOAM (legia citate 0, constant e system); all’interno di ANSA e possibile altresıgenerare direttamente alcuni file (controlDict, fvSchemes ed fvSolutionad esempio) che vengono esportati anch’essi, in modo da avere gia adisposizione il caso completo e pronto per la simulazione.

Con le griglie generate tramite questo programma si sono ugualmen-te ottenuti buoni risultati, rispetto a quelle create con snappyHexMesh.

5.3 checkmesh

Tra le varie utilities fornite con l’installazione di OpenFOAM, ne epresente una di particolare importanza, ossia checkMesh. Questo tool ein grado di controllare i parametri della mesh generata, in termini di:

• statistiche della griglia, ossia viene prodotto un riepilogo dellecaratteristiche della mesh in termini di numero di punti, di faccee di celle; viene anche dettagliata la natura delle celle stesse,siano esse hexaedra, prisms, tetrahedra e cosı via;

• dopodiche viene controllata la topologia, ossia e fatto un checksulla validita dei contorni, sull’utilizzo dei punti, e viene stam-pato un riepilogo delle patch che compongono il dominio;

• infine viene analizzata la geometria, e vengono calcolati, perogni cella, dei parametri che ne caratterizzano la qualita. I valoripiu critici sono due, ossia la skewness e la nonOrthogonality:

– skewness: con riferimento alla figura 39, questo valore edefinito nel modo seguente: prese due celle adiacenti sitraccia un linea, chiamata d, che unisce i due centri cellaP ed N. Viene poi trovata l’intersezione di d con la facciadella cella, ossia fi. E infine calcolata la distanza tra fi e il

50

5.3 checkmesh

centro della faccia f , che viene chiamata m. La skewness equindi definita come il rapporto tra m e d.

Questo valore viene segnalato come critico, nell’output dicheckMesh, qualora superi il valore di 4.

Figura 39: Definizione di skewness

– nonOrthogonality: riferendosi invece alla figura 40, questoindicatore viene definito come l’angolo in gradi formatotra la normale ad una generica faccia, nell’immagine sottoriportata s, e il vettore che unisce i centri cella P ed N, ossia∆.

Figura 40: Definizione di nonOrthogonality

Questo valore e ottimale qualora sia inferiore a 70; nelcaso stia tra 70 e 90 e necessario prendere alcuni accorgi-menti nella soluzione (come un aumento del parametro dinNonOrthogonalCorrectors situato nel file fvSolution), che por-tano inevitabilmente ad un aumento del tempo di calcolo,pur di garantire la stabilita della soluzione; infine qualorala nonOrthogonality fosse superiore a 90, sarebbe necessa-rio assolutamente rifare da capo la mesh perche altrimenti

51

casi 3d

la stabilita della simulazione verrebbe compromessa. Leimmagini sopra riportate sono tratte da [18]

Qualora dunque il controllo della griglia segnalasse la presenza dicelle in cui i due valori citati superano quelli massimi, allora bisogne-rebbe tenere in considerazione il fatto che la stabilita della soluzionepotrebbe essere messa in crisi. Gli accorgimenti da prendere in questocaso potrebbero consistere semplicemente nell’adozione di schemi piuconservativi, oppure nell’aumento del gia citato valore di nNonOrtho-gonalCorrectors, a scapito della velocita di convergenza della soluzione,oppure la decisione di generare una nuova mesh.

5.4 schemi numerici

Nei casi 3D viene esaltata la dipendenza dell’aggiunto anche daglischemi numerici: per non incorrere nella divergenza delle simulazionisono serviti, come sara meglio dettagliato in seguito, dei limitatoriper i termini di gradiente nonche degli schemi upwind per i terminiconvettivi.

Oltre a cio, e stato necessario agire sui valori dei relaxation factorsper le variabili aggiunte: la tecnica dell’under-relaxation della soluzio-ne viene infatti utilizzata per migliorarne la stabilita, e consta nellalimitazione del cambiamento delle variabili da uno step al successivo,controllando il valore di un parametro α, con 0 < α ≤ 1, definito nelmodo seguente:

φn+1 = φn + αφ(φ∗ − φn) (34)

dove φn+1 e il valore che la variabile assumera nell’iterazione successi-va, αφ e il relaxation factor associato a φ, e φ∗ e il valore assunto dallavariabile alla fine del ciclo iterativo al passo n + 1.

Piu tali valori sono bassi, (si puo infatti scegliere un valore del para-metro per ciascun campo) e piu sara garantita la stabilita del calcoloa scapito della velocita di convergenza della simulazione; tuttavia,almeno nelle prime fasi, e necessario tenere quei coefficienti a valorinon troppo elevati per evitare la divergenza del run. Successivamente,una volta che i residui si sono stabilizzati, e possibile aumentare α,per far arrivare a convergenza la soluzione piu velocemente.

52

5.5 caso 3d con wing e mesh poliedrica


In questo caso la mesh e stata generata con l’ausilio del programmaANSA. I passaggi necessari alla sua creazione sono sostanzialmentetre:

1. generazione della mesh superficiale, sia per il corpo che per ildominio esterno;

2. estrusione di un numero prefissato di strati di prismi su tuttal’ala, per simulare la presenza di uno strato limite;

3. creazione della mesh volumetrica.

La griglia cosı ottenuta e composta completamente da poliedri, perun totale di circa 1 milione di celle. L’ala viene ruotata a prioridell’angolo di attacco scelto (pari a 3) affinche il dominio esterno,un parallelepipedo rettangolo, possa essere suddiviso nelle seguentizone:

• un’unica faccia del dominio, quella davanti all’ala, costituisce lapatch di inlet;

• la faccia parallela alla precedente costituisce l’outlet;

• le restanti quattro facce sono piani di simmetria, per simulare ilpiu possibile il campo di moto indisturbato.

L’ala ha una corda costante di 0.16 m, e il dominio esterno si estendeper circa 10 corde davanti e dietro al corpo, come pure sopra, sotto elateralmente ad esso.

Di seguito, nelle figure 41 e 42, viene riportata sia un’immaginedella mesh utilizzata, che un ulteriore zoom attorno all’estremita alare:

Figura 41: Mesh poliedrica, zoom attorno all’ala

53

casi 3d

Figura 42: Mesh poliedrica, zoom attorno all’estremita alare

Essa e stata creata affinche rispettasse i parametri di qualita pre-cedentemente menzionati, che in questo caso assumono i seguentivalori:

• skewness ≈ 1.2 (e percio minore del valore critico di 4);

• nonOrthogonality ≈ 48 (e percio minore del valore critico di 70);


Il setup di questa simulazione in termini di condizioni al contornonon differisce molto da quello del caso bidimensionale.

Per la viscosita turbolenta viene utilizzata la condizione al contornodi tipo freestream sia all’inlet che all’outlet, e il valore che viene im-posto e ricavato dalla stessa legge che e stata introdotta nel capitoloprecedente. La boundary condition appena introdotta si comporta comeuna condizione di Neumann nel momento in cui il flusso esce daldominio, mentre si trasforma in una condizione di Dirichlet in caso diflusso inverso.

• velocita:

– Inlet: valore fissato, pari a 50 m/s in direzione x;


– Symmetry: condizione di simmetria;

– Profilo: condizione di adesione;

• pressione:


– Outlet: valore fissato pari a 0;


– Profilo: condizione di Neumann;

• viscosita turbolenta:

54


– Inlet: freestream con valore fissato pari a 0.012;

– Outlet: freestream con valore fissato pari a 0.012;


– Profilo: valore fissato, pari a 10−20;


Anche in questo caso ci si e avvalsi del solutore stazionario simpleFoam,senza particolari modifiche rispetto al caso bidimensionale. Il valoremedio di y+ ottenuto attorno al corpo e di circa 19 e, nuovamente,la simulazione e stata svolta fino alla stabilizzazione di residui ecoefficienti di forza (figure 43 e 44):

100 200 300 400 500Iterazioni

0.998

1.000

1.002

1.004

1.006

1.008

1.010

1.012

C L

Figura 43: Coefficiente di portanza adimensionale nel caso 3D poliedrico

100 200 300 400 500Iterazioni

0.997

0.998

0.999

1.000

1.001

1.002

1.003

1.004

1.005

1.006

1.007

C D

Figura 44: Coefficiente di resistenza adimensionale nel caso 3D poliedrico

55

casi 3d

Come si puo vedere dai grafici precedenti, la simulazione arriva aconvergenza molto rapidamente; dopo circa 200 iterazioni, infatti, icoefficienti di forza hanno gia raggiunto il valore di regime.


In questa simulazione, la cost function scelta mira all’aumento dellaportanza e questa richiesta viene tradotta, come gia illustrato, in unacondizione al contorno di “traspirazione” per la velocita aggiunta sulcorpo.

Inoltre, la ricerca di un set di schemi che garantisse la convergenzadel problema aggiunto si e rivelata essere abbastanza critica, ed e statonecessario eseguire diverse simulazioni prima di ottenerne una chenon presentasse problemi di divergenza.

Piu nello specifico, partendo dal seguente setup:

• gradSchemes: Gauss linear;

• divSchemes: Gauss linear;

• laplacianSchemes: Gauss linear limited 1;

si e potuto verificare la necessita sia di utilizzare un limitatore per iltermine del gradiente, nonche degli schemi upwind sia per i gradSche-mes che per i divSchemes. Cio che si e riscontrato e stato soprattuttoun comportamento delle equazioni aggiunte molto dipendente dal-lo schema del gradiente, in particolare dalla presenza o meno dellimitatore faceLimited, il cui scopo e quello di effettuare un clip, ossiaun taglio, delle componenti stesse del termine, al fine di limitarne ilcambiamento da un’iterazione alla successiva. Il cattivo comporta-mento della simulazione si e rivelato essere correlato particolarmenteall’assenza di tale limiter. Il valore numerico 1 che viene specificato staad indicare che tale limitatore e del tutto attivo mentre, se tale numeroscendesse progressivamente verso 0, la sua efficacia andrebbe via viadiminuendo. Per evitare la divergenza delle simulazioni aggiunte estato necessario imporre il suo valore ad 1.

Percio la configurazione migliore che e stata ottenuta sfrutta iseguenti schemi per i vari termini:

• gradSchemes: faceLimited Gauss upwind phia 1;

• divSchemes: Gauss upwind;


56


In questo caso, l’andamento ottenuto per i residui e quello riportatonella figura 45:

0.0001

0.001

0.01

0.1

1

0 200 400 600 800 1000

Res

idu

i

Iterazioni

"Uax""Uay""Uaz"

Figura 45: Residui della simulazione aggiunta nel caso 3D poliedrico

Analizzando l’immagine, si puo vedere come tutti i residui scen-dano pressoche monotonicamente, assestandosi stabilmente a valoricompresi tra 10−3 e 10−4, gia a partire dall’iterazione 300 in avanti.

I valori degli altri due indicatori sono, considerando la coppia diiterazioni 1000− 2000:

• NL2 = 0.2396 (≈ 24%);

• ∆Umax = 0.049 (≈ 5%);

Come si puo notare, il valore di NL2 non e piccolo, e cio puo volerdire che il sistema di equazioni aggiunte non converge, mentre ilvalore di ∆Umax ha un valore inferiore. Anche provando a cambiarela coppia di iterazioni considerate, nonche riducendo l’intervallo tracui tale indicatore e calcolato, esso non diminuisce. E percio statasvolta un’analisi di questo campo di moto, riuscendo a verificare comeesso in realta “oscilli”, soprattutto nella zona a ridosso del bordo diuscita dell’ala; la possibilita che ci siano delle oscillazioni nel campodi moto, e piu precisamente vicino al corpo, viene menzionata in [19],seppur non e fatto nessun ulteriore studio a riguardo. Di seguitoviene cercato di visualizzare la presenza di tale comportamento per lavelocita aggiunta ed e stato fatto quanto segue: una volta che i residuisi sono stabilizzati, cosa che accade dall’iterazione 300 in avanti, sonostati calcolati dei nuovi campi scalari, che contengono la differenza, adistanza di 500 iterazioni, tra componenti omologhe di velocita aggiun-ta (e anche tra gli stessi moduli). Le immagini seguenti, puramentequalitative, vogliono dimostrare l’effettiva presenza di oscillazioni nelcampo di moto, e a titolo esemplificativo viene considerata la coppiadi iterazioni 2500− 3000 (figura 46) e la 3000− 3500 (figura 47): esserappresentano uno zoom attorno al bordo d’uscita dell’ala ad unagenerica sezione, e si riferiscono alla componente x di velocita; le zone

57

casi 3d

rosse visualizzano una differenza di velocita positiva tra la prima ela seconda iterazione (cioe velocita alla prima iterazione maggiore diquella alla seconda), mentre, ove il campo di moto e blu, allora lastessa differenza e invece negativa. Se rapportate al valore massimo divelocita aggiunta nel dominio in una delle due iterazioni, tali differen-ze sono attorno al 30%. Questo comportamento si ripete inoltre dalmomento in cui i residui si stabilizzano in avanti e avviene per tutte etre le componenti di velocita.

Figura 46: ∆Ux iterazioni 2500 −3000

Figura 47: ∆Ux iterazioni 3000 −3500

Come gia precedentemente riportato, l’oscillazione e concentratasoprattutto al bordo d’uscita dell’ala e, per cercare di quantificarla intermini percentuali, ogni valore di questo campo di moto e stato anchediviso per il modulo della velocita aggiunta puntuale in una delle dueiterazioni considerate. Cio che si ricava in questo caso e che, allonta-nandosi dal profilo, il campo scende rapidamente al di sotto dell’1%ma, sempre nella zona vicina al bordo d’uscita, i valori arrivano adoscillare di alcune volte il proprio modulo locale; questo comporta-mento e legato al fatto che, in alcuni punti, esso e di un ordine digrandezza inferiore rispetto alla differenza delle componenti al nume-ratore. Cambiando inoltre l’iterazione presa come riferimento per ilcalcolo del modulo della velocita i risultati non si modificano. Vienetuttavia messo ulteriormente in evidenza come, in letteratura, non siapossibile trovare nessuna analisi di questo genere a riguardo, e anche ivalori stessi, in termini di ampiezza delle oscillazioni, sono difficili daquantificare non essendo il campo di moto aggiunto un campo fisico.Cio che si cerca di ottenere in quest’ambito e, invece di una soluzionepriva di oscillazioni, un andamento costante della sensitivita sul corpo.

Sulla simulazione diretta e stata comunque svolta la medesima anali-si, al fine di verificare se ci fossero eventuali oscillazioni ugualmente inquesto campo, che potessero venire trasferite alla simulazione aggiun-ta. Tuttavia questo comportamento non si verifica, siccome il flussodiretto e completamente attaccato all’ala e le oscillazioni al bordo diuscita sono di entita trascurabile, e quindi il problema sembra nonattribuibile a questa simulazione.

Di conseguenza, oltre ai campi di moto aggiunti, anche la sensitivityha un comportamento oscillatorio; per dimostrare cio viene introdot-ta un’ulteriore quantita, ossia l’integrale dello stesso gradiente sulcorpo considerato: viene analizzato dunque l’andamento di questoparametro per verificare quando esso si stabilizza: l’intento e quello

58


di poter a questo punto decidere di ricavare un andamento “medio”della sensitivity, facendo un’average su tutte le iterazioni considerate.Viene dunque riportata, in figura 48, l’evoluzione di questo indicatorenelle prime 150 iterazioni, dove si puo vedere come esso raggiunga inbreve tempo una condizione stazionaria; i valori riportati sono adimen-sionalizzati rispetto alla media del parametro stesso nelle iterazioniconsiderate:

0 20 40 60 80 100 120 140Iterazioni

0.0

0.2

0.4

0.6

0.8

1.0

1.2

1.4Integrale della sensitivity

Figura 48: Andamento dell’integrale della sensitivity nelle prime 150iterazioni

Il grafico seguente (figura 49) riporta invece l’andamento dellostesso parametro una volta che i residui si sono stabilizzati, e piuprecisamente tra l’iterazione 1000 e 4000, dimostrando che vi e al piuuno scostamento di ± 2% rispetto al valor medio stesso (calcolatoanche questa volta nelle iterazioni considerate):

1000 1500 2000 2500 3000 3500 4000Iterazioni

0.96

0.98

1.00

1.02


Figura 49: Andamento dell’integrale della sensitivity tra l’iterazione 1000 e4000

59

casi 3d

L’obiettivo principale che qui e stato raggiunto e il conseguimento diuna simulazione stabile che non presentasse problemi di divergenza;tuttavia, non e stato possibile ottenere risultati migliori in termini diconvergenza delle equazioni, ne cambiando gli schemi numerici, nei valori dei relaxation factors, e cio puo venire imputato alla gia citatadifficolta dell’operatore aggiunto a convergere.

E possibile dunque visualizzare l’andamento della sensitivity media-ta (dall’iterazione 1000 in avanti) cosı ottenuta sul corpo, dapprimasul dorso del’ala (figura 50) e poi sul ventre della stessa (figura 51):

Figura 50: Andamento della sensitivity sul dorso del profilo - mesh poliedrica

Figura 51: Andamento della sensitivity sul ventre del profilo - mesh poliedrica

L’andamento del gradiente riportato e stato adimensionalizzato ri-spetto al suo valore massimo e percio ha valori massimi compresi tra−1 e 1. Analizzando le immagini proposte, si nota come gran partedel corpo sia poco sensibile ad una deformazione per l’aumento dellaportanza, mentre le zone che risultano piu sensibili sono quelle delbordo d’uscita e dell’estremita alare. Come si puo vedere, nella zonavicina al trailing edge sul dorso vi e un’elevata sensitivity negativa: ciosignifica che viene indicato di deformare l’ala in direzione normaleverso l’interno del corpo, mentre nell’analoga zona sul ventre, oveinvece la sensitivity e positiva, vi e l’indicazione di deformare normal-mente verso l’esterno. Queste due informazioni sono “concordi” e illoro significato e quello di far sı che il bordo d’uscita venga modificatoal fine di aumentare la curvatura dell’ala stessa.

60

5.6 caso 3d con wing e trim-mesh


Per lo stesso corpo e stata poi costruita una trim-mesh con il toolsnappyHexMesh. In questo caso, per generare una griglia di buonaqualita, occorre focalizzare l’attenzione sulla creazione di uno stratodi prismi uniforme attorno al corpo. Questo punto e quanto maicritico, se si pensa al fatto che la presenza o meno di tale zona attornoall’ala influenza il valore della y+, che anche per questa simulazioneha un valore medio sul corpo di circa 19, e allo stesso modo questo el’aspetto piu complesso della generazione di una griglia con snappy,come gia citato nell’introduzione.

La mesh cosı costruita e composta da 7 milioni di celle circa, per-lopiu esaedri, mentre la suddivisione delle patch, come pure la parterimanente del setup della simulazione, e mantenuta invariata rispettoal caso precedente, e quindi verra trattata solamente la simulazioneaggiunta, mettendo in luce analogie e differenze rispetto alla meshpoliedrica.

I parametri di qualita della mesh sono, in questo caso:

• skewness ≈ 1.12 (anche in questo caso minore del valore criticodi 4);

• nonOrthogonality ≈ 69 (minore del valore di soglia di 70);

E possibile quindi verificare come la qualita di questa griglia (nellafattispecie in termini di nonOrthogonality) sia piu bassa rispetto al casopoliedrico, seppur comunque ancora sufficiente per non far ritornarealcun messaggio di errore dal controllo di checkMesh.

Cio e legato al fatto che, come gia precedentemente introdotto,snappyHexMesh tende a far collassare gli strati di prismi, e percioa rimuoverli in quelle parti in cui non sono rispettati i parametridi qualita pre-impostati: con l’intento di ridurre al minimo questaevenienza, e stata data al programma la massima “liberta” possibile,e dunque i due indicatori di skewness e nonOrthogonality sono statiimpostati a valori pari a quelli limite, oltre ai quali potrebbero esserciproblemi nella simulazione. Proprio per questo motivo la griglia,durante la sua creazione, e libera di assumere come parametri diqualita quelli massimi ammessi.


La simulazione aggiunta e dapprima stata eseguita sfruttando lostesso set di schemi numerici del caso con la mesh poliedrica. Inquesta situazione, si e potuto verificare quanto preannunciato piuvolte, ossia la simulazione cosı impostata diverge rapidamente (nelgiro di circa 150− 200 iterazioni). Guardando la qualita della grigliautilizzata, si puo supporre che il problema non provenga da qui,come pure dalla simulazione diretta, che e stata anch’essa svolta senzaproblemi. Quindi, la grande differenza che rimane tra questo caso e ilprecedente, e che puo essere imputata come una delle cause principaliper la divergenza, risiede proprio nel tipo di mesh utilizzata. Percioe stato necessario modificare gli schemi numerici e successivamente,

61

casi 3d

con il setup riportato, e stato possibile evitare la divergenza dellasimulazione:

• gradSchemes: faceLimited Gauss upwind phia 1;

• divSchemes: Gauss upwind;


Come si puo vedere, si e cambiato il termine all’interno del dizio-nario laplacianSchemes rispetto alla precedente simulazione. Comeriportato nel manuale di OpenFOAM [20], alla scelta del tipo didiscretizzazione relativa al laplaciano viene fortemente correlato ilcomportamento numerico della simulazione: la scelta che qui e statonecessario attuare utilizza, tramite la sintassi limited 0, schemi delprimo ordine, mentre nel caso poliedrico, con la sintassi limited 1 perl’analogo termine, veniva utilizzato uno schema del second’ordine.

I residui cosı ottenuti sono i seguenti (figura 52):

0.0001

0.001

0.01

0.1

1

200 400 600 800 1000 1200 1400

Res

idu

i

Iterazioni

"Uax""Uay""Uaz"

Figura 52: Residui della simulazione aggiunta nel caso 3D con trim-mesh

I valori di NL2 e ∆Umax, calcolati qui tra l’iterazione 1000 e 2000,sono:

• NL2 = 0.169 (≈ 17%);

• ∆Umax = 0.021 (≈ 2%);

Anche qui si ripresenta l’oscillazione al bordo d’uscita nel campo dimoto aggiunto gia riscontrata nel caso precedente, e valgono parimen-ti le stesse considerazioni gia tratte. Per quanto riguarda l’integraledella sensitivity, esso ci mette piu iterazioni a stabilizzarsi, e cio puoessere correlato alla differente tipologia di griglia utilizzata. Cio nontoglie che venga lo stesso raggiunta una condizione in cui il valoreriportato ha un’andamento oscillante al variare delle iterazioni, comeesemplificato dalla figura 53:

62


0 500 1000 1500 2000 2500 3000 3500Iterazioni

0.0

0.2

0.4

0.6

0.8

1.0

1.2

1.4


Figura 53: Andamento dell’integrale della sensitivity tra l’iterazione 1 e 3500

E percio la media della sensitivity viene ricavata in questo caso apartire dall’iterazione 1000 in avanti, ottenendo i seguenti andamenti(figure 54 e 55):

Figura 54: Andamento della sensitivity sul dorso del profilo - trim-mesh

Figura 55: Andamento della sensitivity sul ventre del profilo - trim-mesh

63

casi 3d

Come si puo vedere, anche in questo caso le zone messe in lucedal gradiente sono le medesime rispetto al caso precedente, percio sisuppone che l’informazione ottenuta sia valida poiche confermata dadue simulazioni diverse fra loro.

Cio che si e voluto dimostrare con questo confronto e questa geome-tria, e la difficolta che e correlata all’operatore aggiunto in termini diconvergenza della simulazione; e infatti difficilmente possibile definirea priori un setup di schemi numerici che sia in grado di evitare ladivergenza della simulazione nel momento in cui si ha a che farecon griglie diverse ma anche, talvolta, con la stessa griglia e con costfunction differenti. E stato inoltre difficile far scendere i residui al disotto di una certa soglia, pur utilizzando differenti setup di scheminumerici, e riducendo sia le tolleranze che permettono di passare daun’iterazione alla successiva che i relaxation factors (i cui valori sonostati abbassati di molto rispetto a quelli delle simulazioni canoniche);si presentano infine delle oscillazioni nel campo di moto aggiuntoche non e stato possibile eliminare, ma la cui presenza e gia statamenzionata in altri lavori.

Queste considerazioni, unitamente ai risultati paragonabili in termi-ni di sensitivity ottenuta nei due casi, hanno permesso sia di ritenerevalida l’informazione ricavata tramite il gradiente, che di procedereall’applicazione del metodo stesso ad un caso piu complicato, i cuipassi principali saranno spiegati nel paragrafo seguente.

5.7 rear model

Vengono quindi presentati i passaggi che sono stati svolti nella simu-lazione piu complessa che e stata affrontata. Non e stato possibileriportare risultati relativamente a questo caso in cui l’operatore ag-giunto e stato applicato ad una configurazione recente della macchinadi F1 della Scuderia Toro Rosso, le cui geometrie non sono pubbliche.

Piu nello specifico, cio che e stato fatto consiste nell’ottenimentodella sensitivity alla parte posteriore della suddetta vettura. Per que-sto caso, le problematiche relative all’operatore aggiunto piu voltemenzionate sono state esaltate, e questo fatto e legato soprattutto allacomplessita della geometria da simulare che presenta molti angoli affi-lati, punti dove le equazioni aggiunte sono risultate essere piu soggettead andare in crisi. La mesh tetraedrica utilizzata e stata ottenuta con ilprogramma ANSA, mentre la cost function implementata nelle equa-zioni aggiunte mira ad aumentare il downforce sviluppato dal corpo: alcontrario di una normale simulazione aeronautica, dove l’ala sviluppaportanza, quindi una forza perpendicolare al vento relativo e rivolta ingenere verso l’alto, in ambito automobilistico la richiesta principale stası nell’aumento della forza perpendicolare al vento relativo, ma questavolta rivolta verso il suolo: da qui la necessita e la richiesta di cer-care di evidenziare le zone piu sensibili all’aumento di questa quantita.

64

5.7 rear model


Questo caso non e stato “banale” da trattare neanche nella simulazionediretta; essa si e infatti rivelata gia complessa di per se, poiche il corpoe letteralmente tagliato dalla patch di inlet, e quindi si e reso necessarioricercare una condizione iniziale su questa parte del dominio che fossesensata.

E stato percio deciso di utilizzare un campo di velocita “estratto” daun altro programma, STAR-CCM+, che e il solutore delle equazioni diNavier-Stokes comunemente utilizzato in azienda. Il campo di motosfruttato e dunque quello relativo ad una simulazione gia arrivata aconvergenza, e risulta essere sicuramente un’ottima scelta. Le coordi-nate dei punti, nonche i valori delle componenti di velocita sulla patchconsiderata, sono importati da tale software sotto forma di tabella,vengono poi manipolati e riarrangiati per essere espressi nel correttoformato richiesto da OpenFOAM e, infine, vengono utilizzati comecondizione al contorno di Dirichlet non uniforme sull’inlet.

Un altro punto particolarmente rilevante e quello relativo alla simu-lazione della rotazione della ruota posteriore della macchina: per farecio essa e stata definita come una patch a se stante, indipendente daquella della macchina, e su di essa e stata applicata una condizione alcontorno per la velocita che permette, dato un asse di rotazione, unpunto che funge da origine e una velocita angolare, di simulare uncorpo che ruota.

Per quanto riguarda invece tutti gli altri aspetti, essi non differisconorispetto a quelli dei casi gia trattati. La simulazione diretta e perciofatta partire e anche qui e stata condotta fino alla stabilizzazione deicoefficienti di forza e dei residui.


In questo caso l’obiettivo menzionato precedentemente, ossia quello diottenere una sensitivity il cui andamento fosse stabile sul corpo, e statoveramente difficile da conseguire. Pur ricercando e ottenendo unabuona qualita della mesh iniziale, la sensibilita dell’operatore aggiuntosoprattutto rispetto agli angoli affilati presenti nella geometria hacompromesso la stabilita della simulazione utilizzando le equazioninella loro forma completa. Neppure il setup in termini di scheminumerici precedentemente utilizzato e stato sufficiente per evitarequesto inconveniente, tanto quanto la riduzione ulteriore dei valoridei relaxation factors per le variabili aggiunte. Una possibile soluzioneal problema, qualora non fosse possibile evitare la divergenza del-l’operatore aggiunto, viene proposta da Lohner, tramite la modificadelle equazioni da risolvere nel modo gia trattato nel capitolo di teoria.Con tale accorgimento e stato possibile stabilizzare il comportamentodel sistema aggiunto e ottenere l’andamento della sensitivita anche inquesto caso.

Dopodiche, al fine di stabilizzare il sistema aggiunto nella sua formacompleta, e stato introdotto un limitatore manuale per la velocitaaggiunta: ogni componente, sia all’interno del dominio che sullepatch che lo delimitano, e stata confrontata con un threshold; qualorala variabile considerata lo ecceda, allora il suo valore viene posto

65

casi 3d

uguale a quello di soglia. Il valore numerico di tale limitatore per ognicomponente di velocita aggiunta e stato scelto sulla base dell’ordine digrandezza delle componenti di velocita nella simulazione “alla Lohner”ricavata in precedenza, che non divergeva. Un diverso controllo dellalimitazione poteva consistere nel confronto tra i moduli di velocitaaggiunta nei due casi tuttavia, cosı facendo, occorreva comunquerisalire in seguito alle singole componenti di velocita; percio si edeciso, in prima istanza, di non utilizzare questo criterio.

Molto schematicamente, e stato percio implementato il seguentealgoritmo, ad ogni iterazione:

• Risoluzione del sistema di equazioni aggiunte e ottenimento delcampo di velocita aggiunta ad una certa iterazione;

• Valutazione di ciascuna componente di velocita aggiunta in ognipunto del dominio, patch incluse;

• if componente > soglia:

→ componente = soglia.

• if componente < −soglia:

→ componente = −soglia.

A questo punto una nuova simulazione e stata fatta partire, questavolta con il solver completo, e si e potuto verificare come il limitatoreentrasse in gioco, agendo in quelle celle che creavano problemi inprecedenza, questa volta pero evitando la divergenza del sistema. Lamappa di sensitivita cosı ricavata e qualitativamente molto simile aquella ottenuta nel caso con le equazioni “incomplete”, in terminisoprattutto di segno del gradiente che viene rispettato in entrambe lesimulazioni quasi ovunque.

La sensitivity, come si e potuto constatare, mette in evidenza siaalcune zone che potrebbero apparire “ovvie” da deformare (come adesempio l’aumento o la riduzione della curvatura di un profilo peraumentarne la deportanza sviluppata) ma anche altre che, ad una pri-ma impressione, non sembrerebbero cosı particolarmente soggette allamodifica per un miglioramento della prestazione richiesta. L’obiettivoche ci si era prefissati all’inizio della tesi e stato percio ottenuto, e ilpasso successivo sara quindi quello di provare a deformare le zoneindicate nella direzione suggerita, per verificare la validita dell’infor-mazione ottenuta. Una complicazione ulteriore, e che non puo esseretenuta in considerazione dall’operatore aggiunto (per lo meno in fasedi ottenimento del gradiente), risiede nella presenza di vincoli geo-metrici, nella fattibilita delle deformazioni suggerite, piuttosto che inlimiti derivanti da misure genericamente “fissate” da un regolamento,che percio riducono l’efficacia, oppure addirittura impediscono, lemodifiche di alcune delle zone evidenziate dalla sensitivity.

66

6 CONCLUS ION I E SV I LUPP IFUTUR I

Il concetto dell’operatore aggiunto nell’ambito dell’ottimizzazione diforma, e l’ottenimento della relativa sensitivity, sono stati applicati incasi 2D e 3D, prendendo come riferimento il sistema delle equazionidi Navier-Stokes stazionarie, viscose e incomprimibili. Il metodo estato dunque applicato su un profilo NACA 0012, imponendo una costfunction che permettesse simultaneamente di aumentarne la portanzae ridurne la resistenza, di fatto aumentando l’efficienza del corpo stes-so. Sono state percio svolte due simulazioni, una in regime di motolaminare e una in regime di moto turbolento, al fine di analizzare sepa-ratamente l’introduzione dell’approssimazione della frozen-turbulence,e la modifica nelle equazioni di moto aggiunte che essa comporta. Perquesto caso e stato sviluppato un algoritmo di deformazione dellageometria per validare il metodo proposto avvalendosi dell’interpretePython, nonche tutto il ciclo di ottimizzazione e stato automatizzatotramite l’implementazione di un loop in Bash. Per il caso tridimensio-nale e stata invece concentrata l’attenzione nell’ottenimento di unasimulazione che non presentasse problemi di divergenza, affrontandole problematiche legate alle equazioni aggiunte, dovute principalmen-te ad una difficile convergenza e ad un’elevata dipendenza da scheminumerici e qualita della griglia. Proprio in quest’ottica sono statecondotte due differenti simulazioni con lo stesso corpo, ma con meshvolumetriche composte da elementi diversi. I risultati ottenuti hannomesso in luce come la stessa simulazione aggiunta, in termini di setupe schemi numerici, se applicata a griglie diverse, puo portare alladivergenza in uno dei due casi. Una volta che e stato ottenuto unsetup soddisfacente, e stata quindi sviluppata la stessa analisi sullaparte posteriore di una macchina di F1, dovendo affrontare i problemisopra citati in misura ancora maggiore, e riuscendo ad ottenere anchequi l’andamento del gradiente, mettendo in risalto le zone piu sensibiliad una modifica della forma.

Gli studi effettuati sono stati finalizzati ad un’applicazione preli-minare dell’operatore aggiunto in correlazione al programma Open-FOAM al fine di ottenere un codice industrialmente utilizzabile, eall’argomento della shape optimization: a tal proposito e stato modifica-to un solutore gia presente nel programma affinche potesse risponderealle esigenze richieste: si e quindi implementato il calcolo della sensiti-vity per questo caso particolare, e sono state utilizzate delle funzionidi costo per ottimizzare le forze aerodinamiche agenti sul corpo, conle relative condizioni al contorno.

Ulteriori sviluppi possono comportare sia un ulteriore approfondi-mento delle debolezze del sistema aggiunto, con l’implementazionedi tecniche piu raffinate per riuscire ad arginarle, cosı come la scrit-tura di un programma per la deformazione del corpo anche nel caso3D. La strada e promettente, e il punto chiave che rende degno di

67

conclusioni e sviluppi futuri

studio questo metodo e la relativa facilita con cui e possibile otteneresul corpo in esame una mappa completa delle zone piu sensibili aduna deformazione: infatti, una volta fatto il setup del caso, ossia so-stanzialmente nel momento in cui si sia ottenuta una configurazionestabile per quanto riguarda gli schemi numerici nonche una grigliadi buona qualita, il costo legato all’ottenimento della sensitivity e sicu-ramente trascurabile, e paragonabile all’incirca all’esecuzione di duesimulazioni canoniche con le equazioni dirette di Navier-Stokes.

68

A COEFF IC IENT I D I FORZAAGG IUNT I

In questa appendice viene presentata l’idea relativa ad un ulteriorecriterio di arresto per le equazioni aggiunte oltre a quelli gia sfruttatiin precedenza: si tratta dell’utilizzo di “pseudo” coefficienti di forzadefiniti analogamente a quelli delle equazioni dirette, con il propositodi verificare se, in corrispondenza della “stabilizzazione” di questivalori, si abbia contemporaneamente anche l’assestamento dei residuidelle equazioni aggiunte, e percio se si abbia la possibilita di usarequesti coefficienti come indicatori di convergenza o stabilita della si-mulazione.

Immaginando che le equazioni di moto aggiunte descrivano l’e-voluzione di un campo di moto fisico come le classiche equazionidi Navier-Stokes, possiamo definire il coefficiente di forza aggiuntocome:

CFa =Fa

12

ρv2S(35)

dove:

• Fa e la forza esercitata sul corpo;

• v e una velocita di riferimento;

• S e una superficie di riferimento;

Dopodiche, con l’intento di “legare” il significato di questo valore aquello delle equazioni canoniche, esso viene proiettato in direzioneparallela e perpendicolare al flusso diretto, dando luogo al coefficientedi “resistenza” aggiunto e al coefficiente di “portanza” aggiunto;si ribadisce ulteriormente che i termini resistenza e portanza sonoutilizzati solo perche le loro direzioni sono le medesime di quelle deicoefficienti classici.

Questi valori non hanno alcun senso fisico, essendo correlati alleequazioni aggiunte e non a quelle dirette; possono pero indicare ilraggiungimento, almeno in termini integrali, di una condizione distazionarieta. Essi vengono adimensionalizzati rispetto alle medesimequantita utilizzate per la simulazione diretta, ossia una superficie S ingenere unitaria e una velocita v uguale a quella del flusso indisturbatodelle equazioni dirette. Anche in questo caso questa e una scelta deltutto arbitraria, fatta solamente al fine di mantenere un legame con leanaloghe quantita fisiche. Questi coefficienti vengono quindi ricavatiattraverso la stessa utility forceCoeffs che viene utilizzata per il calcolodi CL e CD.

71

coefficienti di forza aggiunti

Per rappresentare il loro significato, viene ripreso il caso 3D con lamesh poliedrica. L’andamento dei residui in questa simulazione era ilseguente (figura 56):

0.0001

0.001

0.01

0.1

1

0 200 400 600 800 1000

Res

idu

i

Iterazioni

"Uax""Uay""Uaz"

Figura 56: Andamento dei residui nel caso 3D con mesh poliedrica

Percio, come si puo vedere, per questa simulazione viene raggiuntauna stabilizzazione dei residui attorno a 300 iterazioni.

Parallelamente, l’andamento di dei coefficienti aggiunti e riportatonelle figure 57 e 58:

0 200 400 600 800 1000Iterations

0.0004

0.0002

0.0000

0.0002

0.0004

C Da

Step_0

Figura 57: Coefficiente di “resistenza” aggiunto nel caso 3D con meshpoliedrica

72

coefficienti di forza aggiunti

0 200 400 600 800 1000Iterations

0.004

0.002

0.000

0.002

0.004

C La

Step_0

Figura 58: Coefficiente di “portanza” aggiunto nel caso 3D con meshpoliedrica

Per questi indicatori si raggiunge un andamento stazionario, seppuroscillante, in circa 300 iterazioni. Il fenomeno delle oscillazioni, cheera presente nel campo di moto, risulta esserci anche qui. Per questocaso, sembra dunque che ci sia una corrispondenza tra l’andamentodei residui e quello dei coefficienti di forza aggiunti, seppur questoargomento necessiti ancora di ulteriori approfondimenti prima dipoterne trarre delle conclusioni di validita generale.

73

B ALGOR ITMO DI OPENFOAM

Viene brevemente analizzato e riportato il solver che e stato utiliz-zato per la risoluzione delle equazioni aggiunte, modificando il giapresente adjointShapeOptimizationFoam: il solutore sfrutta l’algoritmodi tipo SIMPLE come anche l’analogo solver diretto simpleFoam, e ditipo segregato, e puo essere suddiviso sostanzialmente in due macroblocchi:

• Nel primo e risolta la prima equazione di (12), da cui e ricavatoun’andamento di velocita aggiunta approssimato;

• Nel secondo e sfruttato il campo appena ricavato per risolvereuna versione modificata della seconda equazione del medesimosistema citato, che permette di modificare il campo di velocita alfine di soddisfare effettivamente il sistema;

Come output di questa prima parte si ottengono i campi di velocitae pressione aggiunti; essi vengono dunque sfruttati nella parte suc-cessiva del codice, ove e stato implementato il calcolo della sensitivitysulla superficie, discretizzando l’equazione (22) come:

∆Li = −νAi((n •∇)u)i • ((n •∇)v)i (36)

Infine, il codice riportato contiene l’implementazione del limitatoreche e stato necessario utilizzare per il caso del modello posterioredella macchina.

Qualora si volesse effettuare l’approssimazione di Lohner, infine, sa-rebbe necessario togliere dal codice, commentandoli, sia la dichiarazio-ne del vettore adjointTransposeConvection sia la sua presenza all’internodel solutore.

// * * * * * * * * * * * * * * * * * * * * * * * * * * //

========= |

\\ / F ield | OpenFOAM: The Open Source CFD Toolbox

\\ / O peration |

\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation

\\/ M anipulation |

References of the original solver,

adjointShapeOptimizationFoam:

"Implementation of a continuous adjoint for topology optimization

of ducted flows"

C.Othmer,

E. de Villiers,

H.G. Weller

AIAA-2007-3947

// * * * * * * * * * * * * * * * * * * * * * * * * * * //

#include "fvCFD.H"

75

algoritmo di openfoam

#include "singlePhaseTransportModel.H"

#include "RASModel.H"

#include "simpleControl.H"

#include <fstream>

// * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])

#include "setRootCase.H"

#include "createTime.H"

#include "createMesh.H"

#include "createFields.H"

#include "initAdjointContinuityErrs.H"

simpleControl simple(mesh);

// * * * * * * * * * * * * * * * * * * * * * * * * * * //

Info<< "\nStarting time loop\n" << endl;

while (simple.loop())

Info<< "Time = " << runTime.timeName() << nl << endl;

// Adjoint Pressure-velocity SIMPLE corrector

// Adjoint Momentum predictor

volVectorField

adjointTransposeConvection((fvc::grad(Ua) & U));

dimensionedScalar NU(laminarTransport.lookup("nu"));

volScalarField nutot(NU + nut);

tmp<fvVectorMatrix> UaEqn

(

fvm::div(-phi, Ua)

- adjointTransposeConvection

- fvm::laplacian(nutot, Ua)

- (fvc::grad(Ua) & fvc::grad(nutot))

);

UaEqn().relax();

solve(UaEqn() == -fvc::grad(pa));

pa.boundaryField().updateCoeffs();

volScalarField rAUa(1.0/UaEqn().A());

Ua = rAUa*UaEqn().H();

76


UaEqn.clear();

phia = fvc::interpolate(Ua) & mesh.Sf();

adjustPhi(phia, Ua, pa);

// Non-orthogonal pressure corrector loop

while (simple.correctNonOrthogonal())

fvScalarMatrix paEqn

(

fvm::laplacian(rAUa, pa) == fvc::div(phia)

);

paEqn.setReference(paRefCell, paRefValue);

paEqn.solve();

if (simple.finalNonOrthogonalIter())

phia -= paEqn.flux();

#include "adjointContinuityErrs.H"

// Explicitly relax pressure for adjoint momentum corrector

pa.relax();

// Adjoint momentum corrector

Ua -= rAUa*fvc::grad(pa);

Ua.correctBoundaryConditions();

// Limiter of velocity over the domain and on the patches

scalar threshold = value;

forAll(Ua,cellI)

Ua[cellI].component(0)=min(Ua[cellI].component(0), threshold);



forAll(Ua,cellI)

Ua[cellI].component(0)=max(Ua[cellI].component(0), -threshold);



forAll(Ua.boundaryField(),patchI)

forAll(Ua.boundaryField()[patchI],faceI)

77


Ua.boundaryField()[patchI][faceI].component(0)=

min(Ua.boundaryField()[patchI][faceI].component(0)

, threshold);



, threshold);



, threshold);

forAll(Ua.boundaryField(),patchI)

forAll(Ua.boundaryField()[patchI],faceI)


max(Ua.boundaryField()[patchI][faceI].component(0),

-threshold);



-threshold);



-threshold);

// Sensitivity calculation

word patchName = "wing";

label patchID = mesh.boundary().findPatchID(patchName);

scalarField temp =

-NU * (U.boundaryField()[patchID].snGrad()

& Ua.boundaryField()[patchID].snGrad());

sensitivity.boundaryField()[patchID] =

mesh.boundary()[patchID].magSf() * (temp);

turbulence->correct();

runTime.write();

Info<< "ExecutionTime = "

78


<< runTime.elapsedCpuTime()

<< " s\n\n" << endl;

Info<< "End\n" << endl;

return(0);

// * * * * * * * * * * * * * * * * * * * * * * * * * * //

79

B IBL IOGRAF IA

[1] A. Quarteroni. Modellistica numerica per problemi differenziali.Springer, 2008.

[2] J.E.V: Peter e R.P. Dwight. “Numerical sensitivity analysis for ae-rodynamic optimization: a survey of approaches”. In: Computers& Fluids 39 373-391 (2010).

[3] R. Hicks, E. Murman e G. VanderPlaats. An assessment of airfoildesign by numerical optimization. Rapp. tecn. TMX 3092, NASA,1974.

[4] G. VanderPlaats e R. G. Hicks. Numerical airfoil optimization usinga reduced number of design coordinates. Rapp. tecn. TMX 73151,NASA, 1976.

[5] D. Bristow e J. Hawk. Subsonic panel method for the efficient analysisof multiple geometry perturbations. Rapp. tecn. CR 3528, NASA,1982.

[6] G. Shubin e P. Frank. A comparison of implicit gradient approachand the variational approach to aerodynamic design optimization.Rapp. tecn. AMS-TR-163, Boeing Computer Service, AppliedMathematics e Statistics, 1991.

[7] O. Pironneau. “On optimum profiles in Stokes flow”. In: Journalof Fluid Mechanics 59 (1) (1973).

[8] A. Jameson. “Aerodynamic design via control theory”. In: Jour-nal of Scientific Computing (3) (1988).

[9] S. K. Nadarajah e A. Jameson. “A comparison of the conti-nuous and discrete adjoint approach to automatic aerodynamicoptimization”. In: AIAA paper 2000-0667 (2000).

[10] C. Othmer. “CFD topology and shape optimization with adjointmethods”. In: VDI-Berichte (2006).

[11] A. Jameson. “Aerodynamic shape optimization using the adjointmethod”. In: Lectures at the Von Karman Institute (2003).

[12] C. Othmer. “A continuous adjoint formulation for the computa-tion of topological and surface sensitivities of ducted flows”. In:Int. J. Num. Meth. Fluids (2007).

[13] E. de Villiers e Petropoulou S. “An efficient adjoint-based me-thod applied on vehicle aerodynamic optimization”. In: Open-FOAM international conference, 2007.

[14] C.A. Marta e S Shankaran. “On the handling of turbulenceequations in RANS adjoint solvers”. In: Computer & Fluids 74102-113 (2013).

[15] O. Soto e Lohner R. “On the computation of flow sensitivitiesfrom boundary integrals”. In: AIAA Paper 2004-0112 (2004).

[16] D. J. Mavriplis. “A discrete adjoint-based approach for optimi-zation problems on three dimensional unstructured meshes”. In:AIAA Paper 2006-0050 (2006).

81

Bibliografia

[17] CFD-online. Turbulence free-stream boundary conditions. 2013. url:http : / / www . cfd - online . com / Wiki / Turbulence _ free -

stream_boundary_conditions.

[18] H. Jasak. “Error analysis and estimation for the finite volumemethod with applications to fluid flow”. Tesi di dott. ImperialCollege, 1996.

[19] A. Stuck. “Adjoint Navier-Stokes methods for hydrodynamicshape optimization”. Tesi di dott. Technische Universitat Hamburg-Harburg, 2012.

[20] Open∇FOAM, The Open Source CFD Toolbox.

[21] T.D. Economon, F. Palacios e J.J. Alonso. “A coupled-adjointmethod for aerodynamic and aeroacoustic optimization”. In:AIAA Paper 2012-5598 (2012).

[22] C. Othmer, E. Villiers e H.G. Weller. “Implementation of a con-tinuous adjoint for topology optimization of ducted flows”. In:AIAA Paper 2007-3947 (2007).

[23] K. Fridolin. CFD for air induction systems with OpenFOAM. 2012.

[24] C. Othmer e T. Grahs. “Approaches to fluid dynamic optimiza-tion in the car development process”. In: EUROGEN. 2005.

[25] O. Soto e R. Lohner. “On the boundary computation of flowsensitivities”. In: AIAA Paper 2004-112 (2004).

[26] R. Lohner, O. Soto e C. Yang. “An adjoint-based design methodo-logy for CFD optimization problems”. In: AIAA Paper 2003-0299(2003).

82

http://www.cfd-online.com/Wiki/Turbulence_free-stream_boundary_conditions

http://www.cfd-online.com/Wiki/Turbulence_free-stream_boundary_conditions

Date post:	12-Jan-2020
Category:	Documents
Upload:	others
View:	3 times
Download:	0 times

Sensitivita aerodinamica nelle autovetture` da competizione: uno … · 2017-07-13 · corpi con...

Documents