POLITECNICO DI TORINO

Interazione fluido-struttura di hydrofoil in materiale composito

Tesi di laurea magistrale di: Pietro Casalone

Relatori: Prof.ssa Giuliana Mattiazzo
Prof. Ugo Icardi
Prof. Andrea Urraci

ANNO ACCADEMICO 2018-2019

1
RINGRAZIAMENTI

Desidero qui ringraziare tutti coloro che mi hanno aiutato nella stesura di questo lavoro o che per
me rappresentano un punto di riferimento.

La professoressa Giuliana Mattiazzo, Relatore della tesi e responsabile del Polito Sailing Team, ha
reso possibile questa esperienza in team e mi ha guidato nella definizione della mia figura
professionale.

Il professor Ugo Icardi, Co-relatore, che mi ha seguito nell' approfondimento dei vari aspetti teorici
analizzati. Il suo supporto e quello del Prof. Andrea Urraci sono stati per me preziosi.

Il professor Luca Ridolfi che per primo mi ha introdotto alla fluidodinamica e mi ha guidato alla sua
scoperta.

Un ringraziamento particolare va a tutti i colleghi ed agli amici del Sailing Team con cui sono
cresciuto come ingegnere e soprattutto come persona e con cui ho condiviso un'esperienza
indimenticabile.
Ringrazio i miei amici di una vita, che per me rappresentano una seconda famiglia, senza di loro
sarebbe tutto completamente diverso.
Ringrazio infine le persone a me più care, i miei genitori, che mi hanno sempre permesso di
scegliere, nello studio e nella vita, quello che più mi appassiona.

2
INTRODUZIONE

Il lavoro della presente tesi è nato da una fortunata coincidenza: l’aver scoperto per caso, quattro
anni fa, che esisteva all’interno del PoliTo, un sailing team che si proponeva di progettare,
realizzare e condurre prototipi di skiff (sail keep it fast and flat).
Appassionato da sempre alla pratica della vela ho deciso con entusiasmo di far parte del progetto.
Seguire i propri sogni ed assecondare le proprie passioni è la cosa più importante nella vita, non
importa se uno è costretto a lasciare porti sicuri o a mettere in discussione quello che ha fatto fino
a quel momento.

È così che è nata questa tesi, dalla volontà di formarsi e reinventarsi in un altro campo, cosa in cui
gli ingegneri, come ci insegna la storia, sono maestri.
Ho così ridefinito la mia figura professionale e personale, ho aggiunto al mio piano di studi
approfondimenti specifici legati alla fluidodinamica e oggi esco dal Politenico con una laurea in
Ingegneria Civile e tanta voglia di mettermi in gioco in un campo che non è quello specificamente
previsto dal percorso di studi che ho seguito.
Per quattro anni, infatti, ho conciliato il mio percorso scolastico con l’attività del Sailing Team,
cercando di crearmi conoscenze e competenze necessarie per lavorare nel mondo della vela.

“L’obiettivo principale del gruppo è la formazione pratica dell’ingegnere di domani: progettare è

mediare tra un’idea e ciò che è realmente realizzabile, dunque la preparazione teorica non è
sufficiente”.

Questa è l’ottica in cui il nostro gruppo si è migliorato di anno in anno, ottimizzando i processi di
lavorazione, sperimentando sandwich sempre diversi, modificando le geometrie di scafo e vele in
base alle condizioni di mare e di vento attese dal campo di regata.

Coerentemente con l’esperienza descritta, la mia tesi tratta dell’ interazione fluido-struttura degli
hydrofoil, le appendici dell’imbarcazione che le permettono di “volare” sopra l’acqua esattamente
come un aereo.
Il problema è molto complesso e di interesse sia per il mondo universitario che industriale: negli
ultimi anni, infatti, la diffusione degli hydrofoil ha avuto un forte boom e sono nate tantissime

3
nuove imbarcazioni che sfruttano questo principio per una navigazione più veloce ed efficiente.
La modellazione di tali appendici risulta un compito estremamente complesso: il fenomeno è
multi-fisico e l’interazione tra la struttura e il fluido introduce una serie di problematiche che
devono essere accuratamente modellate per ottenere una soluzione corretta.
Per effettuare l’interazione fluido struttura è possibile utilizzare alcuni software commerciali che
permettono di fare analisi multi-fisiche, oppure si possono accoppiare due solutori differenti, uno
che si occupa solamente della parte fluidodinamica e uno che si occupa della sola parte
strutturale, andando poi a definire un algoritmo per il trasferimento delle informazioni tra i due
software.
Nell’ambito di questa tesi ho deciso di sfruttare la seconda opzione e di utilizzare software open
source: questo tipo di licenza permette infatti di accedere al cosiddetto codice sorgente del
programma e di andare a modificarlo nel caso ce ne sia bisogno. La possibilità di utilizzare questo
tipo di software senza dover pagare una licenza risulta allettante per la maggior parte degli
utilizzatori medio-piccoli e quindi degna di considerazione e di studio. Da parte mia l’interesse per
questo argomento non si esaurisce con lo studio presentato e spero che il futuro mi offra la
possibilità di ulteriori approfondimenti sul tema.

4
INDICE:
1. Sommario pag. 7

2. Concetti base sui foil pag. 9

3. Concetti base dell’interazione fluido struttura pag. 14
3.1 Modellazione pag. 14

4. Modellazione fluidodinamica pag. 18

4.1 XFLR5, un codice 3D ai pannelli pag. 18
4.2 Mesh pag. 19
4.3 Solutori pag. 20

4.4 Lifting Line Theory (LLT) - Non Lineare pag. 20

4.5 Vortex Lattice Method (VLM) pag. 21
4.6 Metodo 3D ai pannelli lineare pag. 22
4.7 Considerazioni ulteriori pag. 23

5. Confronto tra Star-ccm+ e XFLR5 pag. 25

5.1 Modello su XFLR5 pag. 25
5.2 Modello in Star-ccm+ pag. 29
5.3 Conclusioni pag. 33

6. Modellazione strutturale pag. 34

6.1 Sandwich a facce laminate multistrato pag. 34

6.2 Modellazione dei pannelli sandwich pag. 35

6.3 Code Aster pag. 37
6.4 Elemento DST pag. 38
6.5 Cinematica pag. 38
6.6 Effetti del taglio e delle tensioni trasversali pag. 39

6.7 Considerazioni ulteriori pag. 41

7. Benchmark effettuati su CA pag. 43

7.1 test 1: piastre isotropa semplicemente appaggiata pag. 43
7.2 test2: piastre isotropa multilayer semplicemente appaggiata pag. 46
5
 7.3 test3:piastre isotropa asimmetrica semplicemente appoggiata pag. 48
7.4 test4: piastra in sandwich semplicemente appoggiata pag. 50

8. Correzione modello strutturale pag. 55

8.1 Introduzione pag. 55
8.2 Metodo dei moltiplicatori di Lagrange pag. 55
8.3 Metodo di Riz-Rayleigh pag. 56

8.4 Modello di spostamento ZZA pag. 57

8.5 ZZA pag. 59

8.6 Costruzione di soluzioni analitiche pag. 60

9. Correzione del test 4 con applicazione della ZZA pag. 62

9.1 Conclusioni pag. 66

10. Ala incastrata pag. 68

10.1 Conclusioni pag. 72

11. Trasferimento dei dati tra le due diverse mesh pag. 74

11.1 Principio dei lavori virtuali pag. 74
11.2 XFLR5 -> Code Aster pag. 74
11.3 Code Aster -> XFLR5 pag. 78

11.4 Updated Lagrangian approach pag. 81

12. Conclusioni finali pag. 83
13. Sviluppi futuri pag. 85
14. Elenco foto e figure pag. 86
15. Lista Sigle pag. 88

15. Bibliografia pag. 89

6
 1. SOMMARIO

Lo scopo della seguente tesi è quello di sviluppare, utilizzando software open source ed un

approccio innovativo, uno strumento in grado di eseguire un ciclo di interazione fluido-struttura.

Tale interazione è necessaria ogni qual volta si abbia a che fare con elementi sottili e deformabili

sottoposti a carichi fluidodinamici esterni, variabili quindi in funzione della geometria stessa della

struttura. I campi di applicazione di questo tipo di analisi sono molteplici, così come gli approcci

testati negli ultimi decenni, sia dal mondo universitario che da quello industriale.

L' analisi condotta consiste nel trovare la soluzione ad un problema fluidodinamico per poi

applicarla ad un modello strutturale come condizione di carico. Ottenuta la geometria deformata

la si importa nel solutore fluidodinamico e si inizia un ciclo iterativo finché due diverse deformate

non differiscono tra loro meno di una certa soglia prestabilita. In quel momento anche i carichi

fluidodinamici risultano variare in maniera trascurabile tra un ciclo e l’altro e si considera

raggiunta la convergenza.

Quando si va a progettare una struttura deformabile, soggetta a forze di tipo fluidodinamico, è

necessario effettuare un’analisi di interazione fluido struttura, in quanto la forma disegnata a CAD

e non soggetta ai carichi sarà sensibilmente diversa dalla geometria finale deformata, come

saranno diverse le forze applicate a parità di condizioni al contorno esterno.

Nei casi di interazione fluido struttura la potenza di calcolo necessaria è sempre un problema di

primaria importanza. La maggior parte degli approcci esistenti risultano infatti

computazionalmente molto pesanti: bisogna iterare più volte entrambi i solutori, sia quello

fluidodinamico che quello strutturale, prima di raggiungere la convergenza sotto certe condizioni

al contorno, per poi eventualmente passare alle condizioni di carico successive.

7
È quindi necessario, volendo tempi e potenza di calcolo contenuti, che il singolo solutore, sia dal

lato strutturale che da quello fluidodinamico, utilizzi una teoria semplificata e veloce. Con questa

idea in mente, ho deciso di utilizzare un elemento piastra per la parte strutturale e un modello a

pannelli dal lato fluidodinamico. I dati di output vengono poi corretti per una migliore accuratezza

dei risultati.

Un approccio di questo tipo permette di avere un ridotto costo computazionale e consente quindi

di esplorare rapidamente lo spazio di design; il tool è pensato infatti per una fase preliminare di

progetto dove è importante riuscire ad indagare su una casistica ampia, al prezzo di una minore

accuratezza dei risultati.

Per la correzione degli output del solutore ad elementi finiti, si è utilizzata la teoria di spostamento

ZZA applicandola e adattandola al caso in questione. Questa permette di ottenere risultati simili ai

modelli 3D più raffinati, con un tempo di calcolo ridotto.

La complessità dell’obiettivo prefissato e del problema fisico in esame, hanno portato a trascurare

alcuni effetti che saranno oggetto di un futuro approfondimento; nel capitolo finale si descrivono

tali approssimazioni e si evidenzia cosa sarebbe da implementare per migliorare alcuni aspetti

dell’analisi.

8
 2. CONCETTI BASE SUI FOIL

Negli ultimi anni, nel mondo della vela, l’utilizzo dei foil si è molto diffuso su barche di ogni tipo e

concezione: dalle piccole derive per singolo o doppio atleta fino alle imbarcazioni oceaniche di

oltre 20 metri.

Nel momento in cui i materiali compositi ad alta prestazione e le tecnologie realizzative del

processo sono diventati accessibili ad un’utenza più vasta, sono nate molte nuove classi che

sfruttano il principio degli hydrofoil per ottenere una navigazione più stabile e veloce.

Questo principio è ormai molto diffuso anche sulle tavole wind e kite surf ed è ormai chiaro come

il mondo della nautica e della vela in particolare stia spingendo in questa direzione.

Nonostante il grande sviluppo sia avvenuto negli ultimi 20 anni, già pochi anni dopo la creazione

del primo aeroplano, nel 1906 Enrico Forlanini sul lago di Como, era riuscito per la prima volta a

far “volare” l'idroplano, progenitore dei moderni aliscafi grazie appunto ad un sistema di foil.

Foto 1: Forlanini sul lago di Como nel 1906

9
Il funzionamento di un hydrofoil è infatti lo stesso di un’ala di un aereo che genera portanza in

direzione verticale, con la differenza di essere immersa in un liquido mille volte più pesante. Le

superfici e soprattutto le velocità saranno di conseguenza ridotte rispetto a quelle di un aereo per

garantire la portanza di progetto.

Si può in effetti notare come una barca durante la navigazione con foil sia effettivamente molto

più simile ad un aereo che ad una classica imbarcazione dislocante o planante; si perdono i classici

fenomeni legati all’interazione tra lo scafo e l’acqua e quindi legati al numero di Froude e alla

lunghezza, larghezza e forma dello scafo, mentre diventa più importante l’equilibrio al beccheggio

che, se non soddisfatto, può causare conseguenze non volute e a volte catastrofiche.

Foto 2: Team New Zealand ingavona durante l’ultima edizione dell’America’s Cup alle Bermuda

La velocità massima raggiungibile dalla barca non è più limitata dall’attrito generato dallo scafo

che, raggiunta una certa velocità, esce completamente dall’acqua, bensì dall’attrito generato dai

10
foil. Un’ala ha un’efficienza molto superiore a quella di uno scafo che, oltre a generare la forza

verticale richiesta, deve anche avere una forma dove sia possibile muoversi, ed avrà quindi un

attrito di forma e di parete molto elevato.

La funzione dei foil è quindi quella di sollevare l’imbarcazione dall’ acqua in maniera da avere

meno superficie bagnata possibile.

In questo modo l’attrito è ridotto al minimo, la forza di Archimede e la portanza generate dal

volume immerso dello scafo, che servono a tenere sollevata la barca sull’acqua, sono sostituite

dalla portanza dei foil.

Nelle due successive foto si osserva come due imbarcazioni che hanno velocità di avanzamento e

dislocamento simili vadano a creare delle scie completamente diverse.

Nel caso del Volvo Ocean Race sotto, la scia è pronunciata e si nota come per raggiungere tali

velocità sia necessario un vento molto sostenuto

Foto 3: Volvo Ocean Race in vento forte

11
 Foto 4: AC 50 in vento leggero

Sopra invece due catamarani volanti durante l’America’s Cup del 2013 che veleggiano in acque

calme e vento molto leggero. Nonostante nella foto il vento sia minimo le barche avanzano ad una

velocità sorprendente essendo in grado di procedere a circa tre volte la velocità del vento, finché

questo non aumenta troppo. La scia, indice della resistenza della barca, è molto ridotta e

scompare dopo pochi metri.

In tutte le barche a vela, al contrario degli idrovolanti e aliscafi, non esiste una velocità di regime

che viene raggiunta e mantenuta, infatti il vento è variabile e la velocità della barca funzione di

esso. La portanza da generare è però sempre la stessa ed equivalente al peso della barca e

dell’equipaggio.

Osservando la classica formula della portanza prodotta da un’ala, risulta quindi ovvio che non

potendo andare a modificare l’area ed essendo la velocità condizione imposta ambientale, deve

12
essere modificato il 𝐶𝐿 cambiando l’angolo di attacco per rendere possibile la navigazione con i foil

in un range di condizioni atmosferiche abbastanza ampio.

È per questo motivo che tutte le imbarcazioni fornite di foil sono anche munite di un sistema di

controllo (che può essere attivo, passivo o sfruttare entrambi i principi) per modificare in funzione

della velocità della barca, l’angolo di attacco del foil.

13
 3. CONCETTI BASE DELL’INTERAZIONE FLUIDO STRUTTURA

3.1 Modellazione

Il problema dell'interazione tra fluidi e solidi si pone, sotto diverse vesti, in molti campi applicati

dell’ingegneria, dalla biomedica alla nautica, dall’aereospaziale alla civile; i concetti alla base dei

modelli, rimangono comunque gli stessi.

Sebbene sia ampiamente studiata in letteratura, a causa della sua complessità, l’interazione fluido-

struttura è ancora oggetto di un'indagine intensiva. Una revisione completa della letteratura

esistente va ben oltre lo scopo di questo lavoro; qui vorrei semplicemente menzionare i diversi

approcci che possono essere scelti e le ragioni alla base di queste scelte.

La prima classificazione è legata al fatto che si sia interessati solo alla soluzione di equilibrio o

all'intero comportamento transitorio. Nel caso della FSI di un foil non siamo interessati allo studio

del transitorio, perché le frequenze proprie principali del foil sono molto diverse dalle frequenze

principali della forzante. Il problema può quindi essere studiato con un approccio quasi-statico

dove la convergenza è raggiunta sotto determinate condizioni al contorno prima di passare a

quelle successive, leggermente diverse dalle precedenti.

Inoltre, la soluzione finale a cui converge il foil sotto certe condizioni al contorno, non dipende in

maniera forte dalle condizioni iniziali imposte, purché queste siano verosimili, come potrebbe

accadere invece nel caso di una bandiera al vento (che rappresenta un caso di interazione fluido

struttura forte).

Un altro aspetto importante da prendere in considerazione è l'ampiezza della deformazione e

l'effetto della massa aggiunta: nel caso di soluzione non stazionaria, infatti, il flusso esercita sul

corpo una forza proporzionale alla sua accelerazione. È come se una massa virtuale di fluido

14
dovesse essere spostata per consentire al corpo di muoversi: per tale motivo si parla di effetto di

massa aggiunta (che va a modificare i modi propri della struttura). Questo determina la strategia di

accoppiamento da utilizzare per ottenere simulazioni stabili.

In caso di piccole deformazioni e / o alto rapporto di massa tra solido e fluido, un semplice

accoppiamento può essere utilizzato, il cosiddetto approccio “partitioned weakly coupled”, in cui i

due risolutori possono essere eseguiti separatamente (cioè senza mai assemblare la matrice del

sistema completo): questo consente la possibilità di accoppiare due codici esistenti.

Nel nostro caso, con la trattazione quasi statica, andiamo a trascurare completamente l’effetto

della massa aggiunta che nel caso del foil rimane comunque di importanza secondaria.

Sfortunatamente, quando l'effetto massa aggiunto diventa più importante e il comportamento

transitorio deve essere preso in considerazione, è stato dimostrato che l’interazione “debole”

diventa incondizionatamente instabile ed è necessario introdurre delle sub iterazioni utilizzando

un fattore di rilassamento (under relaxation factor).

L'approccio “strongly coupled” consiste nel diminuire il passo temporale ed introdurre fattori di

rilassamento ad ogni ciclo tra i due risolutori. Lo svantaggio di questo approccio è che il numero di

sub-iterazioni può diventare grande, dell'ordine di qualche decina a seconda del caso, e quindi il

costo computazionale può aumentare molto rapidamente. Questa rimane comunque la strada da

percorrere se siamo di fronte ad un importante effetto di massa aggiunta come per esempio nel

caso dell’interazione fluido-struttura di una vela per andature portanti (gennaker, spinnaker...).

Un'alternativa è lo schema monolitico in cui risolviamo il problema accoppiato non lineare tutto in

una volta.

In questo caso, un sistema lineare globale viene assemblato e risolto per tutte le incognite del

problema: variabili fluide, strutturali e tutte le informazioni geometriche (compresa la posizione

dell'interfaccia).

15
L'intero sistema fluido-solido deve essere risolto nel suo complesso e simultaneamente.

Potrebbero tuttavia sorgere problemi legati alla memoria a causa delle dimensioni della matrice

globale. Inoltre, è necessario un solutore unico ad hoc e non è possibile accoppiare due software

che trattano ognuno una parte separata del problema.

Quando si fa un’analisi con interazione fluido struttura, se non si utilizza un approccio monolitico,

le mesh relative al modello strutturale e fluidodinamico sono solitamente diverse in quanto,

descrivendo i modelli fenomeni diversi, è consigliabile utilizzare anche una diversa discretizzazione

per mantenere un buon compromesso tra accuratezza e costo computazionale.

Questo è il motivo per cui deve essere condotta un'analisi sull'interpolazione dei dati da mesh non

conformi. Quando si passano le forze da due diverse mesh infatti, è importante conservare almeno

l'energia dell'intero sistema in modo che sia possibile rispettare il Principio dei lavori virtuali e non

venga generato o eliminato lavoro nel processo di trasferimento.

Per diminuire il tempo di calcolo, è necessario considerare quale livello di accuratezza sia

necessario per la simulazione delle parti fluide e solide e se qualsiasi modello semplificato possa

essere preso in considerazione.

Ad esempio, nel caso di un foil che lavora a bassi angoli di attacco non è solitamente

indispensabile risolvere le equazioni complete di Navier-Stokes; anzi, modelli più semplici, come

quello basato sul flusso potenziale, sono abbastanza precisi da ottenere campi di flusso accurati.

Considerazioni simili possono essere applicate alla parte strutturale, in cui assunzioni come il

modello di elasticità lineare o di piccole deformazioni possono portare a metodi numerici

semplificati e molto più veloci.

16
 Figura 1: flow chart che descrive il processo di interazione fluido-struttura.

Date le ipotesi sulla tipologia di interazione fluido-struttura in analisi ho optato per accoppiare due

solutori open source esistenti, in maniera che sia possibile, volendo, accedere al codice sorgente e

modificarne delle parti nel caso sia necessario.

17
 4. MODELLAZIONE FLUIDODINAMICA

4.1 XFLR5, un codice 3D a potenziale

Come solutore fluidodinamico ho scelto XFLR5, software open source rilasciato con licenza GNU.

Basato su una teoria bidimensionale e su un metodo a potenziale, XFLR5 permette di svolgere

analisi viscide e inviscide anche 3D; è previsto un modello addizionale a due equazioni integrali per

la trattazione dello strato limite. Inoltre è possibile indagare sul punto di transizione turbolenta del

profilo e sulla bolla di separazione.

XFLR5 è di semplice utilizzo ed ha un costo computazionale irrisorio se paragonato a quello di un

normale RANS solver. Questo lo rende un candidato ideale per essere accoppiato con un solver

strutturale: si parla infatti di mesh di decine di migliaia di pannelli contro mesh di milioni di

elementi necessari per la risoluzione delle equazioni di Navier-Stokes ai volumi finiti.

Lavorando con angoli di attacco bassi che garantiscono un flusso attaccato al profilo per la maggior

parte di esso, un metodo a potenziale risulta essere accurato nella predizione del lift e anche nella

distribuzione delle pressioni lungo la corda dell’ala, nonostante le forti ipotesi semplificative

adottate.

Risulta invece completamente errato il calcolo sul 𝐶𝑑 , che viene regolarmente sottostimato

proprio a causa delle ipotesi che andrò ad esporre nelle righe seguenti.

Nell’ottica di voler sviluppare un tool per un progetto preliminare, dove non interessa eseguire

un’ottimizzazione sull’efficienza dell’ala, la componente di attrito ricopre un ruolo secondario in

quanto la sua influenza sulla deformata finale è completamente trascurabile.

Queste considerazioni mi hanno portato ad abbandonare l’approccio ai volumi finiti e ad utilizzarlo

solamente come strumento di validazione dei modelli a potenziale.

18
Segue una rapida descrizione delle ipotesi semplificative adottate nel metodo a potenziale, delle

teorie utilizzate e dei limiti di applicabilità delle stesse.

4.2 Mesh

L’ala è meshata in un numero finito di pannelli e un vortice o una doppietta e sorgente sono

associati ad ogni pannello. Sul quarto anteriore di ogni pannello è presente un punto di controllo

sul quale verranno applicate le forze fluidodinamiche.

È raccomandabile aumentare la densità della mesh alle estremità di ogni sezione geometrica

dell’ala e soprattutto lungo la corda avvicinandosi al bordo di attacco e di fuga del profilo per

meglio descrivere i gradienti di pressione. Per evitare errori numerici e migliorare la stabilità del

solutore è consigliabile usare pannelli con forma in pianta quadrata o quasi; se le dimensioni del

pannello sono molto differenti nelle due direzioni si va incontro ad errori numerici.

Figura 2: Discretizzazione dell’ala e disposizione dei vortici e dei control point nel VLM method.

19
 4.3 Solutori

In XFLR5 sono presenti tre tipi differenti di solutori:

Il primo è un metodo derivato dalla teoria della linea portante di Prandlt.

Il secondo è un VLM, Vortex Lattice Method.

Il terzo è un metodo 3D ai pannelli.

La novità in XFLR5 sta nel fatto di aggiornare questi algoritmi classici (rendendoli non lineari) con i

risultati di analisi 2D viscose effettuate precedentemente su XFOIL (noto software open source per

l’analisi di profili 2D, a partire dal quale è stato sviluppato XFLR5 che ne rappresenta l’estensione

3D) col fine di stimare anche la parte viscosa dell’attrito e non solo la componente di pressione.

4.4 Lifting Line Theory (LLT) - Non Lineare

Nella classica teoria della linea portante la relazione 𝐶𝐿 − 𝛼 è lineare e gli effetti viscosi non sono

tenuti in conto. In XFLR5 è stata implementata una LLT non lineare: la teoria si basa sul fatto che

un’ala possa essere sostituita da una linea portante e da un distacco incrementale di vortici lungo

il bordo di fuga dell’ala stessa lungo linee dritte e nella direzione della velocità in ingresso al

contorno.

La forza di questi vortici di estremità, che si staccano dal bordo di fuga, è proporzionale al

gradiente della portanza lungo l’ala. Questi vortici inducono una velocità aggiuntiva normale alla

direzione del flusso. Di conseguenza l’angolo di attacco di ogni sezione dell’ala è differente dal

geometrico angolo di attacco di una quantità che è appunto l’angolo indotto dalla velocità.

20
L’angolo di attacco effettivo è quindi legato alla distribuzione della portanza tramite l’angolo

indotto. In più l’angolo di attacco effettivo è legato al Cl della sezione tramite interpolazione su

dati 2D di analisi effettuate precedentemente su XFOIL.

I vincoli imposti dalla LLT lineare e l’interpolazione sulle polari precedentemente generate su

XFOIL devono essere soddisfatte contemporaneamente; si procede quindi con un loop iterativo al

fine di tenere conto del comportamento non lineare.

Limitazioni della LLT:

Si nota come la LLT abbia due importanti limitazioni fondamentali:

- le ipotesi alla base della teoria sono le stesse della teoria di Prandlt, quindi calcoli eseguiti su ali

con basso allungamento o freccia importante potrebbero non essere accurati.

- la forma in pianta dell’ala dovrebbe giacere nel piano XY, quindi con un dietro assente o molto

ridotto.

4.5 Vortex Lattice Method (VLM):

Il metodo VLM è stato implementato come un’alternativa per l’analisi delle geometrie che

ricadono al di fuori delle limitazioni imposte dalla LLT.

Il principio alla base del VLM è quello di modellare la perturbazione generata dall’ala con una

somma di vortici distribuiti sulla geometria della stessa. La forza di ogni vortice è calcolata in

maniera che le seguenti condizioni al contorno vengano soddisfatte:

• condizione di non penetrabilità del pannello dove la velocità su di esso deve essere

tangente e non può avere componente normale.

21
 • condizione di Kutta che impone velocità e pressioni uguali al bordo di fuga per due

particelle che arrivano rispettivamente dal dorso e dal ventre dell’ala.

Il metodo è applicabile a qualsiasi geometria di ala, comprese quelle con forte freccia, basso

allungamento e dietro pronunciato (con possibilità di modellare quindi winglet o hydrofoil con

forme a L o J)

La risoluzione del VLM richiede l’inversione di una matrice quadrata (effettuata con il metodo di

Gauss) che è della dimensione del numero dei pannelli; come sempre in questi casi il costo

computazionale del problema (che va circa col quadrato del numero dei pannelli) può essere

ridotto enormemente facendo considerazioni di simmetria.

4.6 Metodo 3D ai pannelli lineare:

L’obiettivo principale del metodo ai pannelli è il seguente:

rifinire i risultati della LLT e del VLM, tenendo di conto il reale spessore dell’ala e non soltanto la

linea media del profilo nonché dare una migliore descrizione della distribuzione del 𝐶𝑝 sulle

superfici dell’ala.

Il principio alla base del metodo ai pannelli sta nel modellare la perturbazione generata dall’ala

con una somma di doppiette e sorgenti distribuite lungo l’ala.

La forza di queste doppiette e sorgenti è calcolata in maniera da andare a rispettare le giuste

condizioni al contorno che possono essere di Dirichlet o di Neumann.

Una descrizione completa delle teorie alla base di XFLR5 è oltre lo scopo di questa trattazione.

22
 Figura 3: Curve di lift di XFLR5 e dati sperimentali

4.7 Considerazioni ulteriori

L’assunzione più discutibile delle teorie sovraesposte sta nell’ usare i risultati di XFOIL, trovati per

ali di lunghezza infinita anche per ali con allungamento finito.

Il metodo consiste nell’interpolare, utilizzando le polari generate da XFOIL, il 𝐶𝐿 di output del 3D,

per trovare la componente viscosa dell’attrito associata a tal coefficiente; ovviamente questa

risulta un’approssimazione senza nessun tipo di sostegno teorico, in quanto gli effetti 3D sono

importanti, che può però portare a risultati non lontani dalla realtà quanto più l’ala si avvicina ad

un’ala infinita.

Come accade solitamente quando si trasferiscono dei risultati dal 2D al 3D l’attrito, e in particolare

la sua parte viscosa, viene ampiamente sottovalutato, e questo tipo di analisi può portare ad

assunzioni troppo ottimistiche con un’ala che ha un’efficienza nella realtà molto minore.

23
Tra le altre limitazioni è da notare come XFLR5 non possa predire il comportamento di un’ala

quando si utilizzano angoli di attacco vicini allo stallo (anche se questo non è il caso di un hydrofoil

dove ci si tiene ampiamente nella parte lineare del grafico 𝐶𝐿 − 𝛼), quando è presente

interferenza tra piano principale e piano di coda o volendo indagare su fenomeni come la

cavitazione o la ventilazione.

Si evidenzia come sia possibile effettuare su XFLR5 un’analisi di stabilità del volo: il programma è

infatti pensato per analizzare interi aeromodelli, e non solamente ali isolate.

Questo tipo di analisi risulta particolarmente utile nell’ottica di voler descrivere il reale

comportamento di un’imbarcazione che naviga con foil; è possibile fissare una portanza che deve

essere generata, nel nostro caso uguale al peso della barca, e il programma cambia

automaticamente velocità o angolo di attacco al fine di soddisfare tali condizioni.

Inoltre, assegnando alle ali dei valori di massa concentrata che permettano di rappresentare il

reale baricentro del modello, si valuta l’angolo di attacco da conferire al piano principale e al piano

di coda al fine di avere una navigazione stabile.

24
 5. CONFRONTO TRA STAR-CCM+ E XFLR5

5.1 MODELLO XFLR5

Nell’ottica di indagare l’accuratezza dei vari metodi disponibili su XFLR5, tenendo conto delle varie

limitazioni imposte dalle teorie sopracitate, si sono eseguite differenti analisi su un modello

geometrico liberamente ispirato alle forme ad L utilizzate nelle ultime edizioni della America’s Cup

a bordo dei catamarani volanti.

È stata effettuata un’analisi di confronto utilizzando il solutore RANS Star-ccm+, ritenuto un

riferimento affidabile e accurato, anche se computazionalmente molto più pesante, per verificare

quale fosse il reale grado di accuratezza di XFLR5 e in quali situazioni fosse lecito utilizzare il

modello a potenziale (nel caso per esempio si voglia indagare su alti angoli di attacco vicini allo

stallo, su interazione tra deriva e timone o su fenomeni di cavitazione o ventilazione è necessario

l’utilizzo di un software ai volumi finiti viscoso).

Foto 5: Team USA nell’edizione dell’ AC del 2013. In particolare, si vede il foil sopravvento fuori dall’acqua mentre
quello sottovento è completamente immerso.

25
La geometria scelta è un’ala con forte diedro e quindi una superficie portante che genera sia una

forza verticale che laterale.

La forza verticale deve essere in ogni momento uguale al peso della barca in maniera che questa

possa navigare sempre sollevata dall’acqua con un attrito che viene ridotto al minimo; la forza

laterale deve essere paragonabile alla forza laterale prodotta dalle vele al fine di mantenere una

rotta dritta.

Nonostante nel caso di ali simmetriche sia consigliabile studiare solamente metà dominio per

ridurre il costo computazionale dell’analisi, in questo caso non è stato possibile.

Oltre ad avere un angolo di attacco rispetto al piano XY imposto con dei controlli in maniera attiva,

i foil hanno un angolo di attacco, dello stesso ordine di grandezza, anche rispetto al piano YZ,

dovuto al moto proprio dell’imbarcazione. Infatti, il flusso in ingresso che sentono i foil non è

dovuto alla sola velocità di avanzamento della barca ma anche alla componente di velocità

laterale, chiamata scarroccio, che fa sì che anche la parte verticale del foil abbia un angolo di

attacco e vada a generare portanza.

La distribuzione di carico sui due foil sarà quindi diversa in quanto nell’ala a sinistra vediamo il

campo di bassa pressione sempre dal lato “interno”, mentre in quella a destra il campo di bassa

pressione è interno per la parte orizzontale ma esterno per la parte verticale; la seconda

configurazione risulta meno efficiente rispetto alla prima in quanto nella parte centrale e inclinata

dell’ala c’è l’inversione del segno del diagramma di pressione.

L’angolo di attacco imposto nell’analisi è stato di 3 gradi sia per la parte orizzontale che verticale

del foil.

26
 Figura 4: Modello completo su XFLR5

Figura 5: Dettaglio della mesh: distribuzione coseno della densità dei pannelli lungo la corda dell’ala

Per il paragone con Star-ccm+ ho preso in considerazione l’ala con la bassa pressione tutta sul lato

interno, che è la configurazione di navigazione standard.

27
Ho plottato il 𝐶𝐿 per i tre differenti pezzi in cui è divisa l’ala ed ho eseguito il paragone con il 𝐶𝐿 di

Star; inoltre in alcune sezioni significative è stato eseguito il confronto del 𝐶𝑝 in maniera da essere

sicuri non solo di avere la stessa risultante di carico ma anche una distribuzione delle pressioni

vicina a quella reale.

Mi sono concentrato su questi due parametri che vanno a definire il carico da applicare al solutore

strutturale. Purtroppo, date le ipotesi semplificative su cui si basano le varie teorie utilizzate dal

metodo ai pannelli, la componente di attrito risulta in tutti i casi ampiamente sottostimata e si

sconsiglia di affidarsi a XFLR5 per uno studio sull’efficienza dell’ala (almeno nel range di Reynolds

indagato).

Si nota comunque come le stesse approssimazioni vengano applicate nello stesso modo a tutti i

modelli studiati; risulta quindi possibile effettuare studi di tipo comparativo (soprattutto se le

condizioni al contorno e le geometrie sono simili) anche se la precisione di questi ultimi rimane

adeguata ad un progetto di tipo preliminare.

Seguono alcune foto dal solutore XLFR5

Figura 6: Distribuzione del 𝐶𝑝 adimensionale per Star-ccm+ e per XFLR5 (in rosso Star)

28
 Figura 7: Distribuzione della portanza calcolata su XFLR5

Tabella 1: Confronto dei risultati dei due modelli

5.2 Modello in Star

Va oltre lo scopo di questa tesi descrivere accuratamente il modello realizzato sul solutore Star-

ccm+ in quanto è stato utilizzato solamente come fonte di validazione e non è parte attiva

dell’intero tool di simulazione.

Si riportano comunque di seguito alcune foto che mostrano come sia stata ottenuta una buona

convergenza della soluzione e che quindi sia lecito proseguire utilizzando i valori restituiti da Star-

ccm+ per validare gli output di XFLR5.

29
Figura 8: Coefficiente di portanza all’aumentare delle iterazioni: da questa e dalla successiva immagine si può notare
come la simulazione effettuata sia andata a convergenza.

Figura 9: Andamento dei residui delle equazioni del modello su Star-ccm+.

30
Figura 10: Scena di pressione con tracciamento delle linee isobariche

Figura 11: Scena di pressione con tracciamento delle linee isobariche

31
 Figura 12: Dettaglio della mesh fluidodinamica

Figura 13: Valori di y+ che vanno a determinare il tipo di legge a parete utilizzata nel solutore.

32
 5.3 CONCLUSIONI:

Osservando i risultati appena mostrati in tabella, si conclude che il solutore XFLR5 possa essere

usato, nel range di condizioni al contorno di interesse, per trovare la distribuzione di pressioni

lungo un’ala anche in caso di bassi allungamenti e marcati angoli di dietro, la differenza sul 𝐶𝐿 è

risultata essere < 3% e la distribuzione del 𝐶𝑝 praticamente identica, rispetto agli output di Star-

ccm+.

Quindi, nonostante il calcolo sull’ attrito risulti sottostimato (spesso il 𝐶𝑑 calcolato è la metà di

quello reale), si è deciso di utilizzare comunque XFLR5 in quanto per trovare una deformata

dell’ala, è sufficiente la componente di portanza e non serve quella di attrito.

A livello di solutore è stato scelto il VLM2 in quanto garantisce una buona accuratezza con un peso

computazionale molto ridotto; la distribuzione di pressioni lungo l’ala è risultata essere

praticamente la stessa del 3D panel method ma con un tempo di calcolo inferiore.

Le limitazioni imposte sulla LLT mi hanno portato ad escludere tale modello che non permette di

tenere conto dell’angolo di diedro.

Si evidenzia come i due differenti approcci utilizzati, i volumi finiti ed il metodo a pannelli,

necessitino di tempi di calcolo molto differenti tra di loro: XFLR5 conclude un’analisi nel giro di

decine di secondi o pochi minuti mentre Star-ccm+ impiega due o tre ore per la stessa analisi.

Il risparmio computazionale risulta quindi enorme a fronte di un calcolo sottostimato della

resistenza fluidodinamica.

33
 6. MODELLAZIONE STRUTTURALE

6.1 Sandwich a facce laminate multistrato:

Per pannello a sandwich (o struttura a sandwich) si intende un elemento costituito da due strati

resistenti, detti pelli o facce, distanziati tra loro e collegati rigidamente ad un elemento connettivo

che prende il nome di core: la struttura così composta ha un comportamento statico

notevolmente migliore delle singole parti da cui è costituita.

Il core è in genere una schiuma leggera e poco resistente oppure una struttura a nido d’ape, che

permette di distanziare le pelli, composte di materiale nobile e di spessore ridotto.

Le pelli sono adibite alla distribuzione dei carichi nel piano, la presenza del core è invece utile ad

aumentare il valore della rigidezza flessionale del pannello, che dipende dalla distanza delle lamine

dal piano medio, ovvero dall’inerzia della sezione. L'impiego di tale struttura è quindi paragonabile

al concetto della trave con sezione a L, dove l'anima serve ad aumentare la rigidezza flessionale

nella direzione della stessa.

Distanziando le pelli si ottiene un incremento notevolissimo della rigidezza rispetto a un pannello

costituito soltanto da uno spessore di materiale pari a quello delle due facce, con un grande

miglioramento del momento flessionale massimo rispetto al peso della sezione.

Per queste ragioni negli ultimi quaranta anni nell'industria aerospaziale e navale si è sempre più

consolidato l'impiego di pannelli sandwich.

34
 Figura 14: Pannello sandwich realizzato in fibra di carbonio e core in pvc espanso

6.2 Modellazione dei pannelli sandwich:

La modellazione di strutture in materiale composito è spesso affrontata con l’utilizzo di elementi

finiti tridimensionali e bidimensionali, soprattutto quando a comporre la sezione partecipano

materiali con caratteristiche meccaniche profondamente diverse, come nel caso di un sandwich.

Fintanto che si vanno a modellare pannelli laminati monolitici composti quindi da strati con

orientazione variabile delle fibre ma assenza di core, la teoria della sezione omogenizzata

equivalente (ESL) riproduce fedelmente il comportamento strutturale. Tale teoria consiste nel

modellare l’intero pannello con un solo elemento piastra e nel trovare la rigidezza equivalente che

dovrebbe avere un materiale omogeneo per ottenere la stessa energia di deformazione e, se

possibile, la stessa freccia massima del composito sottoposto al carico.

Nel momento che si introduce un core con scarse caratteristiche meccaniche tra le due pelli, la

teoria della sezione omogenizzata diventa imprecisa, da correggere e comunque “case sensitive”.

35
Attualmente, l’approccio più corretto, anche se computazionalmente pesante, consiste

nell’andare a modellare il core con elementi solidi 3D e le pelli esterne con elementi piastra 2D. I

gradi di libertà coinvolti da questi elementi permettono infatti di tenere conto del complesso

campo deformativo e tensionale che caratterizza questo tipo di strutture. Si nota comunque che

usando elementi standard dispacement-based, non sono soddisfatte le condizioni al contorno sulla

continuità della tensione alle interfacce, che comunque diventano sempre meno importanti al

crescere del numero di elementi.

Si evidenzia quindi come un approccio che vada a utilizzare solamente elementi 3D risulti inadatto

alla descrizione meccanica del sandwich in quanto gli elementi 3D non hanno tutti i gradi di libertà

necessari a descrivere il comportamento delle pelli (infatti, solitamente, gli elementi 3D hanno più

nodi ma ammettono solamente le tre traslazioni mentre i nodi degli elementi piastra ammettono

anche le rotazioni intorno a x e y).

Durante i miei studi ho deciso comunque di provare ad utilizzare la teoria della sezione

equivalente omogenizzata, in quanto la soluzione del sistema strutturale risulta estremamente

veloce.

Gli elementi piani, infatti, presentano solamente 4 nodi per elemento con 5 gradi di libertà per

ogni nodo ed è possibile utilizzare un solo elemento lungo lo spessore della struttura rimanendo

entro ottimi margini di errore. Il risparmio computazionale di un approccio 2D è molto forte in

quanto si ottengono sistemi con un numero di gradi di libertà che è di circa un ordine di grandezza

inferiore.

Andando ad utilizzare però la teoria della sezione omogenizzata si ottengono risultati sbagliati e

spesso troppo ottimistici, inoltre si ha una descrizione del campo deformativo e tensionale

completamente distorta ed è necessario passare a modelli di ordine superiore volendo indagare in

36
questa direzione, soprattutto nel caso di pannelli sandwich che presentano campi tensionali e di

spostamento complessi come si vede sotto in figura.

Figura 15: Spostamento in direzione x lungo lo spessore per sandwich simmetrico composto da pelli in carbonio e core
in pvc espanso. Il campo di spostamento in figura è stato trovato dopo l’applicazione della ZZA al pannello sandwich
del test 4 (pg. 50 )

6.3 Code Aster

Come solutore FEM ho scelto Code_Aster, software libero (licenza GNU GPL) di simulazione

numerica dei materiali e delle strutture meccaniche, sviluppato principalmente dal dipartimento

«Analyses Mécaniques et Acoustiques» del servizio Recherche et Développement (R&D) di EDF -

Électricité de France.

Il suo codice sorgente è liberamente scaricabile dal sito del progetto.

Segue la descrizione, riportata dalla documentazione di Code Aster, dell’elemento finito utilizzato

nei diversi test.

37
 6.4 Elemento DST

Questo tipo di modellazione con elementi piastra è pensato per essere usato in campo di piccole

deformazioni e spostamenti. I DST sono elementi piani che non tengono conto della curvatura

geometrica della struttura essendo elementi con solo 4 nodi, al contrario degli elementi shell che

possono descrivere una curvatura utilizzando 9 nodi.

La conseguenza è che potrebbero nascere delle tensioni parassite che possono essere ridotte

usando un numero maggiore di elementi in maniera da poter approssimare la curvatura

geometrica correttamente. La formulazione risulta quindi semplificata e il numero di gradi di

libertà ridotto.

6.5 CINEMATICA

L’ipotesi cinematica fondamentale è la seguente: la sezione perpendicolare alla superficie media

rimane perpendicolare alla stessa a deformazione avvenuta. Da questo approccio deriva un campo

di spostamenti che varia linearmente lungo lo spessore della piastra. Questi elementi utilizzano un

modello cinematico che è quello di Hencky-Mindlin e si possono scrivere le seguenti equazioni:

Campo di spostamento:

𝑢 (𝑥,𝑦,𝑧) 𝑢 (𝑥,𝑦) 𝛽𝑦 (𝑥,𝑦,)

( 𝑣 (𝑥,𝑦,𝑧) ) = ( 𝑣 (𝑥,𝑦) )+ 𝑧 (𝛽𝑥 (𝑥,𝑦,))
𝑤 (𝑥,𝑦,𝑧) 𝑤(𝑥,𝑦) 0

Campo deformativo 3D e della superficie media:

ℰ𝑥𝑥 = 𝑒𝑥𝑥 + 𝑧𝑘𝑥𝑥 2ℰ𝑥𝑧 = 𝛾𝑥

ℰ𝑦𝑦 = 𝑒𝑦𝑦 + 𝑧𝑘𝑦𝑦 2ℰ𝑦𝑍 = 𝛾𝑦

2ℰ𝑥𝑦 = 𝛾𝑥𝑦 = 2𝑒𝑥𝑦 + 2𝑧𝑘𝑥𝑦

38
 𝜕𝑢 𝜕𝛽𝑦
𝑒𝑥𝑥 = 𝑘𝑦𝑦 =
𝜕𝑥 𝜕𝑦

𝜕𝑣 𝜕𝛽𝑥 𝜕𝛽𝑦
𝑒𝑦𝑦 = 2𝑘𝑥𝑦 = +
𝜕𝑦 𝜕𝑦 𝜕𝑥

𝜕𝑣 𝜕𝑢 𝜕𝑤
2𝑒𝑥𝑦 = + 𝛾𝑥 = 𝛽𝑥 +
𝜕𝑥 𝜕𝑦 𝜕𝑥

𝜕𝛽𝑥 𝜕𝑤
𝑘𝑥𝑥 = 𝛾𝑦 = 𝛽𝑦 +
𝜕𝑥 𝜕𝑦

dove le 𝑒𝑥𝑥 sono deformazioni membranali della superficie media, le 𝛾 sono le deformazioni

trasversali associate al taglio verticale, le K sono le curvature della superficie media e le 𝛽 sono le

rotazioni della superficie media. In questo modello cinematico le gamma risultano costanti lungo

lo spessore, e di conseguenza i tagli; ciò viola le condizioni al contorno sulle tensioni a taglio

trasversali in quanto sulle facce libere non caricate avremo comunque degli stress non nulli,

mentre lo dovrebbero essere. Questo implica che compiendo un lavoro sulle deformazioni si ha

un’energia associata errata che genera a sua volta una deflessione scorretta; devono pertanto

essere definiti dei coefficienti correttivi.

6.6 EFFETTI DEL TAGLIO E DELLE TENSIONI TRASVERSALI

Date le ipotesi semplificative appena annunciate, per descrivere in maniera più corretta gli sforzi di

taglio cioè le 𝜎𝑥𝑧 e le 𝜎𝑦𝑧 , è stato a priori determinato un fattore di correzione K tramite

39
equivalenze energetiche con modelli 3D più raffinati, in maniera che la rigidezza trasversale del

modello composto da elementi piastra sia il più vicino possibile alla realtà.

Teoria di Reissner – FSDT (first order shear deformation theory)

La teoria di Reissner è sviluppata partendo dalle tensioni:

la variazione degli sforzi membranali è supposta lineare lungo lo spessore della piastra nel caso

della teoria di Hencky, dove questo deriva dalla variazione lineare delle 𝜀𝑥𝑥 e 𝜀𝑦𝑦 lungo lo

spessore.

Nell’ambito del modello cinematico di Hencky le 𝛾𝑥𝑧 e 𝛾𝑦𝑧 sono costanti lungo lo spessore ma

questo, come già detto, va a violare le condizioni al contorno; è però possibile usare le equazioni di

equilibrio per determinare attraverso queste la variazione degli stress trasversali lungo lo spessore

della piastra.

Si utilizza un modello che, indagando sulla conservazione dell’energia interna, dopo un confronto

con una teoria 3D, restituisce, per un materiale elastico lineare, una relazione tra gli sforzi

risultanti, le rotazioni e la freccia media.

Nel caso di comportamento membranale e flessionale non accoppiati valgono quindi le seguenti

equazioni costitutive:

40
Si introduce il coefficiente moltiplicativo K col valore di 5/6 nella relazione che lega gli sforzi di

taglio trasversali con le distorsioni gamma in maniera da assicurare almeno l’eguaglianza

energetica tra il modello 3D e quello semplificato 2D con distorsioni costanti.

6.7 CONSIDERAZIONI ULTERIORI

Date le ipotesi semplificative applicate nelle teorie degli elementi appena descritti sembrerebbe

possibile utilizzarli solamente in alcune situazioni dove gli effetti 3D non sono predominanti.

Per testare le reali possibilità di applicazione di questi elementi al caso in esame, sono stati

eseguiti sul software Code Aster (attraverso la piattaforma Salome) alcuni benchmark di cui è nota

la soluzione analitica con estrema precisone.

I risultati ottenuti sono stati confrontati con lo stesso modello 2D creato sul software commerciale

MSC Nastran e con la soluzione corretta nota da letteratura o fornita da un FEM mixed 2D/3D

universitario ritenuto affidabile.

Nei benchmark analizzati gli effetti tensionali e deformativi 3D, tipici di un sandwich, vengono

studiati utilizzando geometrie e condizioni vincolari e di carico tali che detti effetti vengano

evidenziati.

Le piastre studiate infatti, soprattutto le più spesse, presentano una deformazione tridimensionale

importante e permettono di capire quando i classici modelli 2D vanno in errore con previsioni

spesso troppo ottimistiche.

Nel caso di una reale ala questi effetti tridimensionali saranno minori rispetto alle piastre in

esame, anche se, nel caso di rapporti di allungamento alare bassi, sicuramente non trascurabili. Ho

comunque preso in considerazione casi particolarmente complessi in maniera da poter testare le

reali capacità del modello nelle condizioni peggiori, all’interno delle quali sicuramente ricadono i

casi classici più semplici.

41
 7. BENCHMARK EFFETTUATI SU CODE ASTER

Test 1: piastra isotropa semplicemente appoggiata

Test 2: piastra ortotropa multilayer semplicemente appoggiata

Test 3: piastra ortotropa asimmetrica semplicemente appoggiata

Test 4: piastra in sandwich semplicemente appoggiata

Test 1: piastra isotropa semplicemente appoggiata

Geometria: Piastra quadrata

L=100mm;

Materiale: dural

E=73000 MPa,

ν=0.3,

G=E/(2(1+ν))=28076 MPa.

Mesh usata: ogni lato è stato suddiviso in 120 elementi, quindi i nodi sono in totale

121*121=14641.

42
Figura 16: mesh usata nei test 1, 2, 3

Condizioni di vincolo: i quattro lati della piastra hanno bloccato soltanto lo spostamento lungo z,

di contro i vertici del quadrato sono vincolati anche in x e y per far sì che la matrice di rigidezza

non sia labile per traslazioni orizzontali.

Carico: ho simulato il carico distribuito discretizzandolo con piccoli carichi concentrati, con valore

variabile lungo x e y, con la seguente legge:

f(x,y)=0.625*sin(π x /L)*sin(π y /L)

* il valore 0.625 è stato trovato utilizzando il PLV per trasformare delle pressioni in forze nodali

equivalenti; lo stesso principio si utilizzerà poi per passare dalle forze di pressione del solutore

fluidodinamico alle forze nodali da applicare sul FEM (vedi appendici)

43
RISULTATI:

Tabella risultati test 1: freccia massima

Figura 17: Spostamento in direzione perpendicolare al piano della piastra per L/H=10

Nel test 1, l’elemento finito piastra è stato in grado di descrivere perfettamente il problema anche

in caso di piastre particolarmente spesse, L/h=4, dove il taglio gioca un ruolo importante.

44
L’elemento DST, infatti, tenendo in conti degli sforzi trasversali tramite il fattore K=5/6 permette di

descrivere accuratamente il comportamento di piastra isotrope omogenee anche spesse.

Si nota inoltre con piacere come il software Code Aster, a differenza del Nastran, sia in grado di

cogliere gli spostamenti nelle direzioni x e y con grande precisione; una piastra appoggiata caricata

perpendicolarmente al proprio piano, si inflette ma si contrae anche trasversalmente. (come si

vede nelle figure a pagina 48)

Questo effetto di contrazione diventa importante nel caso di piastre sottili dove gli sforzi

membranali sono dello stesso ordine di grandezza di quelli flessionali; andando a trascurare

questo tipo di effetto si ottengono risultati molto errati anche per quel che riguarda la freccia.

Test 2: piastra ortotropa multilayer semplicemente appoggiata

Geometria: Piastra quadrata

L=100mm;

Materiale singola lamina: Ortotropo

E11=172400 MPa

E22=E33=6890 MPa

ν12= ν13= ν23=0.25

G12=G13=3450 MPa

G23= 1378MPa

Laminazione: [0/90/0] -> gli angoli sono calcolati rispetto l’asse x; i 3 strati hanno lo stesso

identico spessore, pari ad h/3, dove h è l’altezza totale del laminato.

Mesh usata: ogni lato è stato suddiviso in 120 elementi, quindi i nodi sono in totale

45
121*121=14641.

Carico: ho simulato il carico distribuito discretizzandolo con piccoli carichi concentrati, con valore

variabile lungo x e y, con la seguente legge: f(x,y)=0.625*sin(π x /L)*sin(π y /L)

RISULTATI:

Tabella risultati test 2: freccia massima

Figura 18: Spostamento in direzione perpendicolare al piano della piastra per L/H=10 (mesh visibile)

46
Anche nel caso di laminato simmetrico la modellizzazione utilizzata si è dimostrata in grado di

cogliere una freccia e degli spostamenti in x e y corretti per tutti i rapporti L/H testati.

Test 3: piastra ortotropa asimmetrica semplicemente appoggiata

Geometria: Piastra quadrata

L=100mm;

Materiale singola lamina: Ortrotropo

E11=172400 MPa

E22=E33=6890 MPa

ν12= ν13= ν23=0.25

G12=G13=3450 MPa

G23= 1378MPa

Laminazione: [0/90] -> gli angoli sono calcolati rispetto l’asse x; i 2 strati hanno lo stesso identico

spessore, pari ad h/2, dove h è l’altezza totale del laminato.

Mesh usata: ogni lato è stato suddiviso in 120 elementi, quindi i nodi sono in totale

121*121=14641.

Carico: ho simulato il carico distribuito discretizzandolo con piccoli carichi concentrati, con valore

variabile lungo x e y, con la seguente legge:

f(x,y)=0.625*sin(π x /L)*sin(π y /L)

47
RISULTATI:

Tabella risultati test 3: freccia massima

Anche in questo caso la modellizzazione utilizzata si è dimostrata in grado di cogliere una freccia

corretta per tutti i rapporti L/H testati. Si può concludere quindi che finché si analizzano laminati

multistrato, anche spessi e asimmetrici, con orientazione qualsiasi delle fibre, i risultati rimangono

accurati, a patto di utilizzare sempre lo stesso materiale per tutti gli strati.

Ciò è dovuto alla semplificazione che applica il codice di considerare la sezione composita come

“omogenizzata”.

48
Test 4: piastra in sandwich semplicemente appoggiata

Geometria: Piastra quadrata L=100mm;

Figura 20: sandwich del test 4 composto da core di 8mm e pelli di 1mm.

Materiale pelli: Ortotropo

E11=132380 MPa

E22=E33=10760 MPa

ν12= ν13= 0.25

ν23=0.49

G12=G13=5650 MPa

G23= 3610 Mpa

Materiale core: Isotropo

E=35 Mpa

G=12.3 Mpa

ν=0.4

49
Laminazione: [0/ISO/90] -> gli angoli sono calcolati rispetto l’asse x; le due pelli hanno lo stesso

identico spessore, pari ad 1mm, l’altezza del core è di 8mm, per un’altezza totale della sezione

composita di 10 mm.

Freccia massima corretta (letteratura e software mixed universitario) = 3.93 mm

Nel momento che si va a modellizzare un composito con caratteristiche dei materiali che

compongono i vari strati molto differenti tra di loro, il modello cade drasticamente in errore.

Si sono infatti realizzati tre differenti modelli utilizzando elementi finiti diversi: un modello 2D

classico, un modello 2D con offset e un modello ibrido con facce 2D e core 3D.

Figura 21: Mesh ibrida del test 4

Il primo modello, come nei test case precedenti, utilizza un solo elemento piastra per tutto lo

spessore e per tutti i materiali coinvolti: va a creare una sezione omogenea. Ovviamente tale

metodo risulta un’approssimazione troppo forte che non permette di indagare il comportamento

50
di un sandwich spesso, soprattutto al variare dello spessore del core si ottengono risultati molto

lontani dal reale e poco affidabili anche per un design preliminare. La freccia massima è 2.7 mm

contro i 3.93 della soluzione corretta.

Il secondo modello utilizza elementi piastra sovrapposti uno sopra l’altro. Ogni strato è

rappresentato da un solo elemento piastra: si realizza quindi la mesh definita nei casi precedenti e

si procede poi a copiarla per due volte e ad imporre un offset delle pelli rispetto al core. Si ottiene

un totale di circa 45000 nodi, cioè il triplo del caso 1.

Si impone che gli spostamenti verticali e le rotazioni associate ai nodi che condividono le

coordinate x, y siano coincidenti; questo corrisponde ad assumere una condizione di incollaggio

perfetto e di continuità tra i vari strati, nonché ad imporre una deformazione trasversale nulla nel

core. Risulta quindi che le due pelli rimangono sempre distanti di quanto è stato imposto con

l’offset come se tra di loro vi fosse un materiale estremamente rigido cioè indeformabile.

Ciò produce notevoli errori per quanto riguarda la freccia sotto carico.

L’utilizzo dell’offset implica quindi un aumento dell’inerzia della sezione che va con il cubo

dell’altezza della stessa, in quanto il materiale resistente si trova sulle facce; ovviamente questo

comportamento è molto lontano dalla realtà in quanto una sezione composta da materiali con

caratteristiche meccaniche altamente differenti non si deforma come una piastra sottile

omogenea. Inoltre, nella teoria degli elementi piastra si fa l’ipotesi di 𝜎𝑧 trascurabili, nonché di

stato di deformazione piano: non è possibile quindi cogliere la compressibilità trasversale del

sandwich e in particolare del core che diventa un fenomeno importante. Nella realtà il distacco

delle due pelli tramite l’inserimento di un core migliora particolarmente le performance del

sandwich (soprattutto se si guarda al rapporto peso-momento resistente), ma non permette di

aumentare la resistenza della sezione con il cubo dell’altezza. Intervengono infatti fenomeni di
51
tipo 3D, come appunto lo schiacciamento del core, che vengono completamente trascurati nella

modellazione fatta sin ora.

Per tenere conto degli effetti dovuti alla presenza di un core con caratteristiche meccaniche

estremamente inferiori a quelle delle pelli, si è descritto nel terzo approccio il core utilizzando

elementi 3D. Questa modellazione risulta estremamente più precisa anche se

computazionalmente molto più pesante. Se infatti per gli elementi piastra si possono avere

elementi con forme anche molto lontane dalla piastra (es. nei test precedenti per L/H=4

l’elemento piastra era un parallelepipedo a base quadrata di lato 0.8mm e alto 25 mm), nel caso di

elementi solidi 3D è necessario avere forme che non si distaccano troppo dal cubo per mantenere

una buona accuratezza nei risultati.

Ciò porta a discretizzare il core con 10 elementi lungo lo spessore, per mantenere una forma

cubica dell’elemento, e ad aumentare quindi drasticamente il numero di nodi che passa da 45.000

a circa 150.000 con un conseguente aumento del tempo computazionale e della memoria richiesta

dal software.

In questo modo è però possibile indagare a fondo il comportamento del materiale lungo lo

spessore della sezione in quanto si dispone di molti più punti per imporre le condizioni fisiche del

problema.

Anche in questo caso ibrido 2D/3D è stata usata la tecnica di offset imponendo che spostamenti

dei nodi con stessa x e y fossero coincidenti. Queste condizioni sono però state imposte solamente

tra la pelle superiore e lo strato superiore di elementi 3D, e allo stesso modo tra la pelle inferiore e

lo strato inferiore di elementi solidi. A differenza di prima ora esistono altri 8 elementi solidi a

descrivere il core che non sono vincolati a muoversi come le pelli e possono quindi rappresentare

52
un comportamento tridimensionale. La massima freccia ottenuta con questo approccio è stata di

4,3 mm.

Si conclude quindi che non risulta possibile all’interno di Code Aster utilizzare solamente elementi

piastra se bisogna descrivere il comportamento di un sandwich dove esiste un core con

caratteristiche meccaniche profondamente diverse da quelle delle pelli esterne; in questi casi è

necessario utilizzare elementi solidi per il core oppure un sistema di post processing che permetta

di rielaborare i dati per fornire una soluzione più accurata.

Si evidenzia come sia invece possibile descrivere laminati ortotropi multilayer anche spessi purché

il materiale degli strati sia lo stesso; il problema è dato dall’insorgere di effetti 3D che mal vengono

descritti dagli elementi piastra che hanno alla base una teoria 2D.

Nell’ ottica quindi di utilizzare comunque elementi 2D, data la convenienza dal punto di vista

computazionale, si è deciso di testare una metodologia innovativa che è stata sviluppata

precedentemente da Icardi e Urraci, per correggere gli output del FEM.

Questa teoria partendo dagli spostamenti del modello semplificato del FEM, restituisce,

imponendo ulteriori condizioni, una deformata e degli stress più corretti.

53
 8. CORREZIONE MODELLO STRUTTURALE

8.1 Introduzione

Al fine di indagare accuratamente il campo deformativo e tensionale a cui è soggetto un sandwich,

mantenendo tempi computazionali ridotti, si introducono allora nuove condizioni (rispetto a

quelle delle teorie classiche come la FSDT prima esposta, che corrisponde all’elemento DST).

Queste una volta soddisfatte, permettono di determinare in forma chiusa e ben approssimata le

grandezze coinvolte nell'analisi.

Si riportano brevemente due metodi tramite i quali è possibile implementare e calcolare le

ulteriori incognite del problema strutturale coinvolto.

8.2 Metodo dei moltiplicatori Lagrange

I moltiplicatori di Lagrange rappresentano nuovi gradi di libertà del sistema. Con un aumento del

costo computazionale permettono di soddisfare delle condizioni al contorno tramite un approccio

di tipo variazionale; si studiano i massimi e i minimi vincolati di una funzione a più variabili in

riferimento ad un vincolo espresso tramite una o più equazioni. La potenza di tale metodo consiste

nell'ottenere come risultato il valore numerico dei moltiplicatori, anche se questi non

rappresentano delle incognite fisiche del problema.

Si osserva infine che tale metodo fornisce condizioni necessarie ma non sufficienti per la

determinazione dei punti di stazionarietà del funzionale: permette infatti di identificarli ma non di

classificarli, ed è quindi necessario verificare di aver trovato realmente un minimo o un massimo

locale.

54
 8.3 Metodo di Rayleigh-Ritz

Il metodo di Rayleigh-Ritz è un metodo di tipo variazionale che permette di approssimare la

soluzione esatta quando questa non può essere calcolata per via analitica. È un metodo molto

usato che si può applicare ogni volta che le equazioni da risolvere si possono porre in forma

variazionale. Il principio del metodo è il seguente: sia ψ un elemento di uno spazio di Hilbert H a

cui appartengono le soluzioni dell’equazione che vogliamo studiare. Si suppone che tali soluzioni

siano elementi di H per cui un certo funzionale E[ψ] sia stazionario. Si utilizzano delle funzioni di

prova per approssimare ψ, le cosiddette trial function.

Nel caso di un problema a carattere strutturale il funzionale va a rappresentare la definizione del

potenziale totale del sistema (pertanto tale tipo di formulazione prende nome anche di metodo

energetico, o del minimo dell’energia potenziale), somma del potenziale totale dei carichi applicati

e dell'energia di deformazione elastica.

55
 8.4 Modello di spostamento ZZA

Si ipotizza la seguente descrizione del campo di spostamenti, attraverso lo spessore dove

𝛼 = 𝑥, 𝛽 = 𝑦, 𝜁 = 𝑧:

u ( x, y, z ) =  u 0 ( x, y ) + z (0 ( x, y ) − w0 ( x, y ), ) 0 +  Fu ( z ) i

ni n

   ( x, y)( z − zk ) H k ( z ) +   Cuj ( x, y) H j ( z ) c
k

k =1 j =1

ni
u ( x, y, z ) =  w ( x, y )
0
0 +  
F (z) i +    k ( x, y)( z − zk )H k ( z ) +
k =1
ni n

  ( x, y)( z − z ) H ( z) +  C ( x, y)H ( z) 
k =1
k
k
2
k
j =1
j
j c

Si osserva che la soluzione è espressa come somma di tre termini.

Il primo termine [… ]0 , ha contributo lineare per 𝑢𝛼 e costante per 𝑢𝜁 ,ed è lo stesso della teoria di

Reissner-Mindlin e introduce gli unici cinque gradi di libertà funzionali della ZZA: 𝑢0 , 𝑣 0 , 𝑤 0 , 𝛤𝑥0 , 𝛤𝑦0

ovvero le tre traslazioni e le rotazioni rispetto agli assi x e y. Il secondo contributo [… ]𝑖 contiene

termini di ordine superiore, assunti essere un'espansione di potenze attraverso lo spessore:

 Fu ( z) i =  Ci ( x, y) z 2 + Di ( x, y) z 3 + (Oz 4 ...) i =  3 (.~) i +  (Oz 4 ...) i

 F  ( z) i =  bi ( x, y) z + c i ( x, y) z 2 + d i ( x, y) z 3 + ei ( x, y) z 4 + (Oz 5 ...) i =  4
~
(. ) i +

 (Oz 5 ...) i

Questo tipo di formulazione abbandona l'interpolazione unica lungo lo spessore, a favore di

un'interpolazione che valga solo all'interno di ogni singolo strato della struttura. Si mantengono i

cinque gradi di libertà funzionali classici e si calcolano i rimanenti coefficienti del modello di

56
spostamento in funzione dei primi, tramite calcolo simbolico. L’ appellativo "adattivo" deriva dalla

possibilità di espandere la serie di potenze attraverso lo spessore, sia per la componente di

spostamento u che w.

Le espressioni di 𝐶𝛼𝑖 , 𝐷𝛼𝑖 e dei coefficienti da 𝑏𝑖 a 𝑒 𝑖 sono determinate tramite l'imposizione delle

condizioni al contorno alle tensioni, sia per la faccia superiore che inferiore. Alle già descritte

condizioni, considerando i termini di ordine ancora maggiore, si aggiungono il soddisfacimento

delle equazioni di equilibrio locale sempre per punti scelti attraverso lo spessore.

  , +   , +  z , z = b ;

  , +   , +   , = b

Il terzo termine […]c, rappresenta dei contributi ulteriori, detti layerwise, che determinano la

continuità delle tensioni fuori dal piano e del gradiente della tensione trasversale normale, lungo

lo spessore, attraverso le interfacce.

Questi ultimi sono espressi come prodotto di funzioni zig-zag (𝑧 − 𝑧𝑘 )𝐻𝑘 (𝑧) e (𝑧 − 𝑧𝑘 )2 𝐻𝑘 (𝑧) e

di ampiezze zig-zag 𝛷𝛼𝑘 , 𝜓 𝑘 e 𝛺𝑘 .

Tali ampiezze vengono determinate imponendo le seguenti condizioni al contorno: continuità della

tensione a taglio trasversale, continuità della tensione trasversale normale e del suo gradiente:

  (( k ) z + ) =   (( k ) z − ) ;

  ((k ) z + ) =   ((k ) z − ) ;

  , (( k ) z + ) =   , (( k ) z − )

57
 𝑗 𝑗
Le condizioni di continuità degli spostamenti alle interfacce sono legate ai termini 𝐶𝛼 𝐶𝑢 e 𝐶𝜁 i

quali sono indipendenti dai gradi di libertà; tali condizioni possono essere espresse come

u ( ( k ) z + ) = u ( ( k ) z − ) ;

u (( k ) z + ) = u (( k ) z − )

Il termine (𝑧 − 𝑧𝑘 ) è la funzione zig-zag di Di Sciuva mentre il termine (𝑧 − 𝑧𝑘 )2 è la funzione zig-

zag parabolica di Icardi; il termine 𝐻𝑘 (𝑧) rappresenta la Heaviside unit step function. Le ampiezze

zig-zag provvedono a garantire il giusto cambio di pendenza degli spostamenti alle interfacce degli

strati con materiali e/o orientazione differente.

La potenza di tale modello consiste nel fornire risultati accurati senza ricorrere a tecniche di post-

processing facendo utilizzo delle sole equazioni costitutive; è quindi possibile affermare che il

modello ben descrive l'energia di deformazione elastica.

8.5 ZZA*

È inoltre possibile dimostrare come, ridefinendo i coefficienti delle funzioni che descrivono il

campo di spostamento lungo lo spessore, non ci sia più bisogno che i suddetti contributi ZZ

appaiano in maniera esplicita: è solamente necessario un polinomio di ordine minimo n che sia

sufficiente a descrivere il modello di spostamento lungo lo spessore.

È quindi possibile sostituire tali funzioni ZZ con dei polinomi e trovare i coefficienti tramite

l’imposizione di condizioni al contorno lungo lo spessore.

58
 8.6 Costruzione di soluzioni analitiche

Si ottengono soluzioni in forma chiusa per l'equazione di governo, indipendentemente dalla teoria

utilizzata, esprimendo i gradi funzionali di libertà, attraverso lo spessore, come un'espansione di

serie troncata di ampiezze sconosciute e funzioni di prova 𝔑𝑖 (𝑥, 𝑦) che soddisfano

individualmente le condizioni al contorno prescritte:

m
 =  Ai i ( x, y)
i =1

Le ampiezze restanti sono determinate derivando il funzionale rispetto alle ampiezze ancora

incognite e imponendolo uguale a zero; all’interno del funzionale è considerato il lavoro delle forze

di inerzia:

 [ − b u ]dv =  [ −  u u ]dv

V
i i
V
i i

In questo modo, si ottiene un sistema algebrico, la cui soluzione fornisce il valore numerico di

ciascuna ampiezza, quindi i campi di spostamento, deformazione e sollecitazione possono essere

calcolati. Per definire le funzioni di forma, a seconda del reale caso in analisi, si vanno poi ad

imporre diverse condizioni al contorno.

Per esempio, in corrispondenza del bordo incastrato delle travi a sbalzo, vengono applicate le

seguenti condizioni:

u 0 (0,0) = 0; w0 (0,0) = 0; w0 (0,0), x = 0; 0x (0,0) = 0;

u (0, z ), z = 0; u (0, z ), z = 0; u (0, z ), xz = 0

Per garantire che la forza risultante dalla forza di taglio trasversale sia uguale alla forza di vincolo,

si impongono le seguenti condizioni meccaniche sul taglio agli estremi della trave:

h/2
 −h / 2
 xz (0, z)dz = T

59
 h/2
−h / 2
 xz ( L, z )dz = 0

Si noti che le precedenti condizioni al contorno vengono applicate utilizzando i metodi dei

moltiplicatori di Lagrange; non si applicano ulteriori condizioni ai momenti flettenti ma, se

necessario, è sufficiente scegliere un ordine di espansione sufficiente ed applicare le seguenti:

h/2
−h / 2
 xz (0, z) zdz = M

h/2
−h / 2
 xz ( L, z ) zdz = 0

Qualsiasi altra condizione al contorno potrebbe essere applicata allo stesso modo, semplicemente

scegliendo funzioni di prova che soddisfino individualmente le condizioni al contorno prescritte.

Nel caso in cui questo non sia possibile, o per soddisfare le condizioni meccaniche, si applica il

metodo dei moltiplicatori di Lagrange. In questo modo è possibile ottenere una soluzione in forma

chiusa con uno sforzo computazionale ridotto, anziché ricorrere alla FEA come spesso accade per

tali casi.

Gli spostamenti allora devono essere descritti da funzioni di classe C0, continue attraverso le

interfacce del materiale, atte a simulare e riprodurre in modo efficace l'effetto zig-zag, invece che

utilizzare funzioni di classe C1 continue.

La scelta delle prime rende possibile la descrizione dell'inversione di pendenza degli spostamenti

alle interfacce, soddisfacendo così le condizioni di equilibrio locale per il campo di tensioni, la

continuità del taglio trasversale, degli sforzi normali trasversali e del gradiente di questi ultimi.

60
 9. CORREZIONE DEL TEST 4 CON APPLICAZIONE DELLA ZZA

Nell’intento di migliorare la soluzione della modellazione con elementi piastra, e sfruttando i

precedenti studi di Icardi e Urraci, si utilizzano gli spostamenti nodali del FEM degli elementi

appartenenti ad un asse della piastra per ricreare, con un polinomio di ordine 6, le funzioni di

forma del problema.

Si è notato infatti come queste funzioni differiscano dalla reale deformata solamente in ampiezza

e non in forma, questo può essere corretto imponendo ulteriori condizioni al contorno fisiche.

Si utilizzano quindi gli spostamenti nodali del FEM per ricostruire le funzioni di forma nei casi in cui

queste non siano note a priori.

Andando però a plottare tali funzioni ci si accorge come queste siano leggermente disturbate al

contorno da effetti di bordo, soprattutto risultano “sporcate” le derivate prima e seconda che

rappresentano rispettivamente rotazione e curvatura.

Fig 22: w(x) del Test 4 (normalizzata con la freccia massima della modellazione 2D): estratta dal

FEM in grigio, soluzione analitica in giallo

61
È possibile alterare le funzioni di forma, estratte con l’interpolazione polinomiale, con dei termini

aggiuntivi che impongono il rispetto delle condizioni al contorno meccaniche; così si può

correggere questo errore nel comportamento dovuto alla modellazione ad elementi finiti e,

ottenendo un campo di spostamento corretto, raggiungere l’accuratezza delle teorie più avanzate.

Nel caso in questione (est 4) di piastra semplicemente appoggiata soggetta a carico bisinusoidale si

sono imposte le seguenti condizioni che hanno permesso di modificare il polinomio che

rappresenta le funzioni di forma.

Condizioni al contorno legate all’appoggio semplice e al carico simmetrico:

𝑤(0, 𝑦) = 𝑤 (𝐿𝑥 , 𝑦) = 0 𝑤(𝑥, 𝐿𝑦 ),𝑦𝑦 = 0

𝑤(𝑥, 0) = 𝑤(𝑥, 𝐿𝑦 ) = 0 𝐿𝑥
𝑤,𝑥 ( , 𝑦) = 0
2
𝑤(0, 𝑦),𝑥𝑥 = 0
𝐿𝑦
𝑤,𝑦 (𝑥, )=0
𝑤(𝑥, 0),𝑦𝑦 = 0 2

𝑤(𝐿𝑥 , 𝑦),𝑥𝑥 = 0

Si nota come, in questo caso semplice, a causa delle condizioni geometriche, di vincolo e carico, la

deformata e quindi le condizioni aggiuntive da imporre fossero note a priori: infatti, ancor prima di

svolgere l’analisi, era possibile immaginare una deformata sinusoidale per detta piastra.

Nei casi reali con geometrie o condizioni di vincolo e carico variabili, la funzione di forma è

interpolata dagli spostamenti nodali del FEM e corretta con termini aggiuntivi per risolvere il

problema al bordo imponendo il rispetto delle condizioni al contorno meccaniche.

Procedendo in tal senso risulta che, modificando la formulazione della soluzione, è possibile che si

aggiungano ulteriori gradi di libertà fisici ma non strutturali per i quali è necessario definire

ulteriori condizioni da imporre.

62
 Figura 23: stratigrafia di riferimento per la formulazione delle condizioni al contorno

Per determinare i coefficienti di ZZA si impongono le seguenti condizioni al contorno sulle tensioni

trasversali, sulle facce superiori e inferiori:

ℎ
𝜎𝑧 (𝑥, 2) = 𝑞 (𝑥, 𝑦)
ℎ
𝜎𝑧 (𝑥, − 2) = 0
ℎ
𝜎𝑧,𝑧 (𝑥, ± 2) = 0
ℎ
𝜏𝑥𝑧 (𝑥, ± 2) = 0

alle interfacce per le componenti di tensione trasversali e per le componenti di spostamento si

impone la continuità:

ℎ ℎ
𝜎𝑧𝑖 (𝑥, 2𝑖 ) = 𝜎𝑧𝑖+1 (𝑥, − 𝑖+1
2 )
ℎ ℎ
𝑖
𝜏𝑥𝑧 (𝑥, 2𝑖 ) = 𝜏𝑥𝑧
𝑖+1
(𝑥, − 𝑖+1
2 )
ℎ𝑖 ℎ
𝑢𝑖 (𝑥, ) = 𝑢𝑖+1 (𝑥, − 𝑖+1 )
2 2
ℎ ℎ
𝑤 𝑖 (𝑥, 2𝑖 ) = 𝑤 𝑖+1 (𝑥, − 𝑖+1
2 )

𝑝𝑒𝑟 𝑖=1,…..,𝑛 𝑖𝑛𝑡𝑒𝑟𝑓𝑎𝑐𝑐𝑒

63
E gli equilibri:

𝜎𝑥,𝑥 + 𝜎𝑥𝑧,𝑧 = 0
𝜎𝑧,𝑧 + 𝜎𝑥𝑧,𝑥 = 0

Si riportano di seguito nelle immagini tensioni e spostamenti lungo lo spessore di una sezione di

mezzeria della piastra.

Si nota come la freccia massima, parametro di riferimento per tutti i benchmark, sia pressoché la

stessa del fem mixed 2D/3D universitario, nonché dei casi riportati in letteratura.

𝐹𝑖𝑔. 24: σx lungo lo spessore 𝐹𝑖𝑔. 25: σxy lungo lo spessore

Fig. 26: Spostamento in direzione y lungo lo spessore 𝐹𝑖𝑔 27: σyz lungo lo spessore
64
 Fig. 28: Spostamento verticale lungo lo spessore

9.1 Conclusioni

Si è dimostrato come la ZZA sia in grado di cogliere il campo tensionale e deformativo anche nel

caso di sandwich spessi; l’accuratezza di questa modellazione è simile a quella del FEM mixed 2D/3D.

La difficoltà nel trovare una soluzione in forma chiusa e ben accurata risiede nel soddisfare le

condizioni al contorno meccaniche le quali comportano problemi soprattutto nel caso in cui siano

presenti vincoli di semplice appoggio. Si nota, infatti, come in corrispondenza di tale vincolo gli

spostamenti estratti dal FEM e soprattutto le derivate prima e seconda risultino sporcate da effetti

di bordo e quindi da correggere con termini aggiuntivi.

Lo studio svolto si è mosso in una direzione che rispetti intrinsecamente le condizioni al contorno

meccaniche comprese quelle con condizioni di vincolo relativamente arbitrarie; i coefficienti delle

polinomiali in z, variabili strato per strato, soddisfano le condizioni fisiche quali continuità,

condizioni al bordo ed equilibri locali lungo lo spessore.

Si osserva come l'utilizzo del calcolo simbolico in ambiente MATLAB abbia permesso un approccio

più rapido al problema posto: non c’è bisogno di imporre restrizioni al campo di spostamento,

65
tensioni e deformazioni al fine di ovviare alle difficoltà algebriche, poiché queste sono a carico del

simbolico. In tal modo si perviene ad una formulazione generalizzata dove l’utente può inserire

arbitrariamente le condizioni di interesse; il calcolo viene poi svolto automaticamente a fronte di

errori numerici, limitazioni di memoria e velocità.

La programmazione in MATLAB risulta semplice ed immediata, ma lo stesso codice ben scritto in

un linguaggio base è di più rapida esecuzione.

Una volta trovati, tramite l’applicazione della ZZA, i campi tensionali e deformativi corretti è

possibile andare a creare un elemento ESL (equivalent single layer) che restituisca la stessa freccia

di ZZA. Quest’ultima è dovuta all’energia di deformazione flessionale, membranale e legata alle

distorsioni trasversali.

Andando quindi a tarare i coefficienti della matrice di rigidezza di un elemento ESL, in maniera che

si realizzino due modelli energeticamente equivalenti, è possibile ottenere risultati corretti almeno

per quanto riguarda il campo di spostamento. Non è comunque possibile ottenere maggiori

informazioni sul campo tensionale utilizzando un elemento ESL che non ha i gradi di libertà per

descrivere il campo tensionale relativo ad un sandwich.

Il costo computazionale derivante da tale formulazione è ridotto, paragonabile ad esempio a

quella di una FSDT, fornendo comunque dei risultati accurati.

66
 10. ALA INCASTRATA

Si è deciso di eseguire un test per una geometria assimilabile ad un hydrofoil, cioè una trave

incastrata ad un estremo e libera sull’altro; il carico, uniformemente distribuito, agisce sulla faccia

superiore. Non volendo in questa sede analizzare un reale hydrofoil ma semplicemente testare

l’affidabilità della modellazione usata, si è scelta una stratigrafia asimmetrica che permettesse di

evidenziare gli effetti 3D e layerwise e di mettere in crisi la ZZA: una teoria in grado di descrivere il

complesso campo deformativo e tensionale tridimensionale associato ad un sandwich

asimmetrico, sarà sicuramente anche in grado di cogliere dei valori corretti nel sotto caso

semplificato di stratificazione simmetrica.

Si riporta nella seguente tabella la stratificazione utilizzata e le condizioni di vincolo e carico:

Si è proceduto ad imporre le diverse condizioni, le quali coinvolgono il giusto valore del taglio

all'incastro poiché, essendo in quella sezione nulli gli spostamenti e le rotazioni, risulterebbe nullo

anche il taglio stesso; stessa cosa per il momento. Inoltre, si impone il valore del taglio = 0

sull’estremo libero.

Si osserva che le condizioni meccaniche sono state implementate nel calcolatore tramite il metodo

dei Moltiplicatori di Lagrange.

• condizioni al contorno meccaniche

ℎ/2
∫ 𝜎𝑥𝑧 (0, 𝑧)𝑑𝑧 = 𝑇
−ℎ/2

67
 ℎ/2
∫ 𝜎𝑥𝑧 (𝐿, 𝑧)𝑑𝑧 = 0
−ℎ/2

ℎ/2
∫ 𝜎𝑥 (0 , 𝑧)𝑧 𝑑𝑧 = M
−ℎ/2

ℎ/2
∫ 𝜎𝑥 (𝐿𝑥 , 𝑧)𝑧 𝑑𝑧 = 0
−ℎ/2

Le condizioni meccaniche vanno a sostituire 3 condizioni imposte sul soddisfacimento

dell'equazione di equilibrio locale. 𝜐1

Vanno aggiunte anche le condizioni al contorno legate alla geometria e ai vincoli:

𝑢(0, 𝑧) = 𝑤(0, 𝑧) = 𝑤(0, 𝑧),𝑥 = 0

Inoltre, per assicurare che le precedenti condizioni valgano per tutto lo spessore si impone anche:

𝑢(0, 𝑧),𝑧 = 𝑢(0, 𝑧),𝑧 = 𝑢(0, 𝑧),𝑥𝑧 = 0

Si sono inoltre aggiunte le seguenti condizioni alle interfacce per determinare i coefficienti di ZZA:

-condizioni di continuità alle interfacce:

ℎ ℎ
𝜎𝑧𝑖 (𝑥, − 2𝑖 ) = 𝜎𝑧𝑖+1 (𝑥, 𝑖+1
2 )
ℎ ℎ
𝑖
𝜏𝑥𝑧 (𝑥, − 2𝑖 ) = 𝜏𝑥𝑧
𝑖+1
(𝑥, 𝑖+1
2 ) 𝑝𝑒𝑟 𝑖 = 1, … . . , 𝑛𝑖𝑛𝑡𝑒𝑟𝑓𝑎𝑐𝑐𝑒
𝑖 (𝑥, − ℎ𝑖 ) = 𝜎 𝑖+1 (𝑥, ℎ𝑖+1 )
𝜎𝑧,𝑧 2 𝑧,𝑧 2
ℎ ℎ
𝑢𝑖 (𝑥, − 2𝑖 ) = 𝑢𝑖+1 (𝑥, 𝑖+1
2 )
ℎ ℎ
𝑤 𝑖 (𝑥, − 2𝑖 ) = 𝑤 𝑖+1 (𝑥, 𝑖+1
2 )

68
-Condizioni al bordo sulle tensioni:

ℎ
𝜎𝑧 (𝑥, + 2) = 𝑞(𝑥 )
ℎ
𝜎𝑧 (𝑥, − 2) = 0
ℎ
𝜎𝑧,𝑧 (𝑥, ± 2) = 0
ℎ
𝜏𝑥𝑧 (𝑥, ± 2) = 0

- condizioni di equilibrio locale in determinati punti lungo lo spessore:

𝜎𝑥,𝑥 + 𝜎𝑥𝑧,𝑧 = 0
𝜎𝑧,𝑧 + 𝜎𝑥𝑧,𝑥 = 0

Si riportano nelle immagini seguenti alcune delle grandezze che caratterizzano lo stato tensionale

e deformativo della struttura: in blu il risultato del FEM 2D/3D mixed, in rosso la ZZA.

Si nota come i risultati siano stati normalizzati per renderne più immediata la lettura ed ottenere

una freccia unitaria; l’accuratezza di ZZA è praticamente la stessa del FEM mixed restituendo

campi deformativi e tensionali corretti.

𝐹𝑖𝑔. 29: 𝜎𝑥𝑧 lungo lo spessore per x =L/2

69
 Fig. 30: U(z) per x =L/2

Fig. 31: W(x) normalizzato con la freccia massima del FEM mixed 2D/3D

70
10.1 Conclusioni

I risultati ottenuti tramite il modello attuale sono praticamente coincidenti con quelli ottenuti

tramite la FEA-2D/3D mixed la quale è stata già dimostrata essere una soluzione molto accurata.

La FSDT risulta essere inappropriata per la descrizione di tale problema poiché la propria

definizione del campo di spostamenti comporta risultati che non rappresentano in modo accurato

gli effetti tridimensionali, soprattutto la compressibilità del core. Risulta quindi conveniente

applicare la ZZA, quando possibile, poiché il risparmio computazionale è forte.

Nello studio effettuato si è partiti da un solutore ad elementi finiti utilizzando varie teorie

associate a diversi tipi di elementi. Si è dimostrato come in alcuni casi (definiti da geometria

semplice come per esempio una piastra, una trave, un’ala o un foil) sia possibile determinare ed

utilizzare una soluzione analitica in forma chiusa invece di un modello ad elementi finiti,

ottenendo oltre ad un’ottima accuratezza dei risultati, anche un tempo computazionale ridotto.

Questo implica però di dover definire tale soluzione per parti se la geometria è data dalla

successione di componenti elementari (tante travi incastrate una dopo l’altra) come nel caso del

nostro hydrofoil.

Infatti, la soluzione analitica è funzione della geometria e delle condizioni al contorno di partenza

che vanno a determinare le funzioni di forma.

Volendo velocizzare e rendere più flessibile tale tipo di modellazione sarebbe opportuno andare a

definire un nuovo elemento finito basato sul modello analitico implementato. In questo modo la

discretizzazione di geometrie complesse risulterebbe immediata e non bisognerebbe andare a

comporre soluzioni analitiche ognuna definita su una parte differente del dominio.

71
La definizione di un elemento finito che vada a riprodurre i risultati del modello analitico

comporterebbe un numero di gradi di libertà associati all’elemento tale che sia possibile rispettare

le condizioni al contorno imposte alle interfacce.

Tale elemento dovrebbe avere un numero di nodi lungo la verticale atto a garantire una corretta

descrizione del campo di spostamento in direzione trasversale alla piastra.

Un’alternativa è utilizzare la modellazione ZZA solamente per tarare i termini della matrice di

rigidezza di un elemento ESL equivalente al composito in analisi; questo garantisce una freccia ed

un’energia associata alle deformazioni corretta ma non permette di indagare sul campo

tensionale.

Per soddisfare i complessi criteri di rottura inerenti ad un sandwich sarà comunque necessario

passare a modelli di ordine superiore che non considerano la sezione omogenizzata.

72
 11. TRASFERIMENTO DEI DATI TRA LE DUE DIVERSE MESH

(XFLR5 -> Code Aster / Code Aster -> XFLR5)

11.1 Principio dei lavori virtuali

Il principio dei lavori virtuali (spesso indicato come PLV) afferma che il lavoro virtuale svolto dalle

forze esterne su un corpo solido è uguale a quello svolto dalle forze interne. Il termine "virtuale"

indica che il teorema è valido per lavori calcolati per un dato sistema qualsiasi di forze esterne

(forze di superficie e di volume) equilibrato con tensioni unitarie, e per un qualsiasi campo

di spostamenti congruente con le deformazioni unitarie, ma non necessariamente conseguente al

sistema di forze esterne applicate.

Il teorema dei lavori virtuali può essere esteso a sistemi discreti di corpi (internamente continui)

tra loro vincolati.

11.2 XFLR5 -> Code Aster

Solitamente, nei casi di interazione fluido struttura si lavora con due mesh differenti, per i due

differenti solutori. Infatti, i due diversi problemi fisici, strutturale e fluidodinamico, richiedono, al

fine di mantenere l’errore di discretizzazione sotto una certa soglia prestabilita, una diversa mesh.

Nel caso fluidodinamico è necessario infittire la mesh nei pressi del bordo di attacco e di fuga di un

profilo aereodinamico in quanto in quelle zone si hanno i gradienti più forti. Nel caso strutturale

invece, non avendo forze concentrate ma solamente pressioni, è necessario infittire

principalmente nei pressi del vincolo e non nella lunghezza dell’ala.

73
La mesh strutturale è quindi una mesh uniforme (che si infittisce nei pressi del vincolo), mentre

quella fluidodinamica segue un andamento uniforme nella direzione dell’ala e un andamento

coseno nella direzione della corda andandosi ad infittire sul bordo di attacco e di fuga.

Inoltre, anche supponendo una mesh uniforme per entrambi i casi, è difficile ottenere lo stesso

numero di celle al fine di mantenere approssimativamente uguali gli errori di discretizzazione del

modello, in quanto si sta descrivendo un diverso fenomeno fisico.

Quando si vanno a trasferire le forze di pressione dalla mesh fluidodinamica alla mesh strutturale

bisogna far sì che i due sistemi di carico siano equivalenti dal punto di vista energetico cioè che

non si sia creato o eliminato lavoro trasferendo le forze da una mesh all’altra.

Figura 32: Output di XFLR5

Si va quindi ad applicare il principio dei lavori virtuali e deve valere la seguente equivalenza

W1=W2

W1→ è il lavoro delle forze della mesh fluidodinamica per gli spostamenti nodali dei nodi della

stessa. Il problema è che dalla soluzione del FEM non si ottengono gli spostamenti nodali della

mesh fluidodinamica ma solamente di quella strutturale.

Esistono comunque varie soluzioni per trovare, dati gli spostamenti nodali della mesh strutturale,

gli spostamenti nodali della mesh fluidodinamica. Nello specifico ne sono stati indagati due:

74
Il primo consiste nell’ utilizzare le funzioni di forma degli elementi finiti per sapere, una volta

fissato il sistema di riferimento, lo spostamento rispetto al nodo adiacente appartenente all’altra

mesh.

Nel secondo, in alternativa, è possibile far passare un polinomio interpolante f(x,y) (definito nel

sistema di riferimento locale appartenente al piano del pannello) per la deformata della mesh

strutturale. Trovato il polinomio, lo spostamento lungo z di qualsiasi punto appartenente al

dominio viene di conseguenza fornendo x e y.

Si è notato come questo secondo approccio risulti migliore e fornisca risultati più precisi a patto di

utilizzare ai fini dell’interpolazione un dominio più grande di quello a cui si è realmente interessati:

infatti lungo i bordi del dominio la soluzione del polinomio risulta leggermente sporcata da effetti

di bordo.

W2 → è il lavoro delle forze equivalenti a quelle del solutore fluidodinamico, trasferite sulla mesh

strutturale, per gli spostamenti dei nodi della mesh strutturale.

Quando si vanno a trasferire le forze da una mesh all’altra è necessario introdurre davanti ad esse

un coefficiente moltiplicativo che tenga conto di due cose: la differente area dell’elemento (sia nel

caso di forze di pressione che di forze nodali) e il diverso lavoro che viene compiuto a

deformazione avvenuta. Infatti, pur utilizzando due sistemi di forze la cui risultante è equivalente,

il lavoro compiuto da esse è differente se la discretizzazione è diversa.

L’applicazione del PLV permette proprio di trovare il coefficiente moltiplicativo che devono avere

le forze nodali affinché il sistema di carico sia energeticamente equivalente e non si vada ad

aggiungere o togliere del lavoro quando si trasferiscono le forze rispetto alla reale condizione di

carico rappresentata dal solutore fluidodinamico.

75
Lo stesso processo è alla base della definizione del coefficiente 0.625 applicato ai benchmark

effettuati su Code Aster per passare da forze di pressione a forze nodali equivalenti. Si è utilizzato

tale coefficiente dopo aver fatto uno studio sulla conservazione del lavoro svolto per un caso

semplice di piastra isotropa di cui è nota la soluzione con precisione.

La funzionalità PROJ_CHAMP di Code Aster può trasferire delle forze da una mesh ad un’altra

affidandosi a considerazioni energetiche basate sul principio dei lavori virtuali e potrebbe essere

un’alternativa automatica agli approcci appena proposti.

Ultima via da testare è usare una mesh strutturale uguale alla mesh fluidodinamica definita su

XFLR5. Andando infatti a correggere gli output del solutore ad elementi finiti si riesce a mantenere

una buona accuratezza della soluzione anche utilizzando una mesh che non risulta ideale o

ottimizzata per un problema di tipo strutturale.

È quindi da valutare la bontà di questa alternativa con ulteriori test.

76
 11.3 Code Aster -> XFLR5

La creazione di un’ala con relativa mesh in XFLR5 è un processo abbastanza semplice che si può

intuire guardando la sottostante tabella che definisce, per pezzi, la geometria dell’ala in analisi.

Figura 33: definizione per pezzi dell’ala su XFLR5

77
Il procedimento consiste nell’assemblare una serie di pannelli rettangolari piani uno dopo l’altro.

Tra un pannello e il successivo può esserci un cambiamento di corda, freccia, diedro, angolo di

attacco e profilo aereodinamico. In questo modo si va a ricreare una geometria molto simile a

quella reale con l’utilizzo di pochi punti fondamentali (basta infatti utilizzare i 4 vertici di ogni

rettangolo per definire una superficie tridimensionale in questo modo).

Dopodichè è sufficiente assegnare il numero di pannelli e la distribuzione di questi ultimi per ogni

rettangolo e il processo di mesh è ultimato in meno di un secondo.

Si nota come sia anche possibile importare in XFLR5 una geometria in formato .stl che viene poi

trasformata dal programma in una successione di rettangoli come quella in tabella.

Ho deciso comunque di non sfruttare questa funzionalità in quanto la modellazione sul FEM è

basata su un modello 2D che non tiene conto del reale profilo aereodinamico; sarebbe quindi

inutile esportare un .stl deformato da Code Aster in quanto sarebbe la deformata di un modello

2D e non del reale modello geometrico.

Data la modellazione 2D utilizzata su Code Aster, risulta conveniente realizzare un algoritmo su

Matlab che dia gli input ad XFLR5 per costruire la tabella che definisce l’ala a deformazione

avvenuta.

Gli unici parametri che risultano fondamentali per definire l’ala deformata sono: l’angolo di diedro

che corrisponde alla deformazione flessionale e l’angolo di attacco che corrisponde alla

deformazione torsionale dell’ala.

Per definire l’ala in questo modo ci si può concentrare solamente sugli spostamenti nodali dei

vertici che compongono i macrorettangoli che definiscono per parti la geometria.

78
L’algoritmo sviluppato è il seguente:

- Per ogni rettangolo si tracciano le coordinate dei 4 vertici che lo definiscono

- A deformazione avvenuta si guarda lo spostamento nelle tre componenti di ogni nodo.

- Si trasformano le tre componenti dello spostamento rispetto alla terna giacente e

perpendicolare al piano del rettangolo.

- L’angolo di dietro è l’arcotangente della differenza degli spostamenti perpendicolari al

piano di due nodi appartenenti ad un lato lungo del rettangolo, fratto la lunghezza del

rettangolo stesso.

- L’angolo di attacco è l’arcotangente della differenza di spostamento perpendicolare al

piano, per due punti appartenenti al lato corto del rettangolo, fratto la corda dell’ala.

79
 11.4 UPDATED LAGRANGIAN APPROACH

I foil sotto carico sono soggetti a grandi deformazioni e sarebbe quindi necessario utilizzare un

modello non lineare per tenere conto di questo comportamento.

Un’alternativa è utilizzare il cosidetto approccio “Updated Lagrangian”.

Questo approccio consiste nel meshare ad ogni iterazione la geometria per ripartire da un modello

scarico da tensioni ma deformato. Applicando nuovamente i carichi, questa volta solamente la

differenza con la condizione precedente e non il contributo intero, non si esce dal campo lineare

delle piccole deformazioni e ciò permette di usare una teoria semplificata.

Nell’ approccio “Total Lagrangian” infatti, le equazioni discrete sono formulate nella

configurazione di riferimento indeformata. Per l’approccio “Updated Lagrangian” invece, ad ogni

iterazione si aggiorna il sistema di riferimento e le equazioni discrete vengono scritte rispetto alla

geometria deformata corrente, che diventa la configurazione di riferimento.

Utilizzando questo tipo di approccio è necessario controllare, rispetto alla configurazione iniziale,

se qualche criterio di rottura sia stato soddisfatto: guardando rispetto all’ultima configurazione

deformata ed applicando solamente una differenza di carichi, non si arriverebbe infatti mai a

rottura.

Si usa quindi una teoria lineare e poi si controlla se le deformazioni e tensioni rispetto alla

configurazione a riposo sono tali da provocare la rottura del materiale o si può proseguire con

l’analisi.

80
Figura 34: Flow chart che descrive il ciclo iterativo completo

81
 12. CONCLUSIONI FINALI

• Si può affermare che tale metodologia risulta una scelta efficace ed accurata per la

descrizione del comportamento di un hydrofoil. Il campo di pressione che agisce

sull’ala è descritto fedelmente e ciò è considerato sufficiente nell’ottica di voler

trovare una deformata della struttura e testare le condizioni che generano le

condizioni critiche nelle sezioni più sollecitate.

Si sottolinea come uno studio che utilizzi come base un metodo a potenziale sia

inadatto ad uno studio di ottimizzazione basato sull’efficienza: l’attrito infatti risulta

sempre sottostimato anche nel caso di angoli di attacco bassi ed ali allungate e

rettilinee.

• La teoria ZZA utilizzata ha prodotto risultati praticamente coincidenti con quelli del

FEM mixed 2D/3D utilizzato come riferimento preciso.

• Si nota come la soluzione strutturale restituisca, con una metodologia innovativa, la

giusta freccia e le giuste tensioni legate ad una condizione di carico di tipo

flessionale. Durante questo studio tutte le condizioni al contorno imposte sono

legate ad un comportamento di tipo flessionale; la correzione adottata tiene conto

solamente di questo tipo di comportamento, Infatti non si sono imposte ulteriori

condizioni sul torcente.

82
• Il modello pone comunque le basi per la costruzione di una soluzione che consideri

anche il comportamento torsionale di piastre o travi; è sufficiente infatti aggiungere

delle ulteriori condizioni sul momento torcente per uscire dalla ZZA con una

deformata che tenga conto anche di questo comportamento.

• Si evidenzia come nel nostro caso questo tipo di implementazione risulti essere di

fondamentale importanza in quanto è necessario tenere correttamente conto delle

variazioni di angolo di attacco tra una deformata e l’altra, essendo queste ultime le

maggiori responsabili del cambiamento delle forze fluidodinamiche rispetto alla

geometria indeformata.

• È possibile utilizzare la ZZA per tarare i termini della matrice di rigidezza K di un

elemento ESL per ottenere un modello equivalente almeno dal punto energetico

che restituisca quindi dei valori di freccia corretti.

83
 13. SVILUPPI FUTURI

Purtroppo il tempo necessario per sviluppare per intero ed efficacemente un tool del genere

eccede quello che è possibile dedicare ad una tesi magistrale. Ritengo tuttavia auspicabili ulteriori

approfondimenti dello studio condotto e mi auguro di riuscire in futuro a lavorare ai vari aspetti

del processo elaborato per incrementarne l’accuratezza e la facilità di utilizzo. Magari per far sì che

possa essere utilizzato anche dal mondo industriale e non solo da quello della ricerca universitaria.

Le implementazioni più importanti da realizzare prima che il processo venga utilizzato per

progettare o verificare un hydrofoil sono, secondo me, le seguenti:

• modificare la modellazione analitica di ZZA per far si che possa tenere conto del

comportamento torsionale

• andare a creare un elemento finito da aggiungere ad un FEM che sia basato sulla teoria

analitica ZZA, in maniera da poter modellare qualsiasi geometria che non sia

descrivibile da una funzione e in maniera da poter più facilmente assegnare condizioni

di vincolo e carico particolari come quelle uscenti da XFLR5.

• Testare in maniera più approfondita i diversi modi che sono stati precedentemente

descritti per definire il coefficiente correttivo per trasferire le forze da una mesh

all’altra.

• Tradurre i codici MATLAB utilizzati in un linguaggio base più adeguato all’ HPC (high

performance computing)

• Scrivere un codice (C++ o Python) che permetta di lanciare e concludere il ciclo in

automatico senza il bisogno di aprire differenti software e mandare manualmente le

analisi una dopo l’altra.

84
14. ELENCO FOTO E FIGURE

Foto 1: Forlanini sul lago di Como nel 1906 pag. 9

Foto 2: Team New Zealand ingavona durante l’ultima edizione dell’America’s Cup alle Bermuda pag. 10

Foto 3: Volvo Ocean Race in vento forte pag. 11

Foto 4: AC 50 in vento leggero pag. 12

Figura 1: flow chart che descrive il processo di interazione fluido-struttura. pag. 17

Figura 2: Discretizzazione dell’ala e disposizione dei vortici e dei control point nel VLM method pag. 19

Figura 3: Curve di lift di XFLR5 e dati sperimentali pag. 23

Foto 5: Foto 5: Team USA nell’edizione dell’AC del 2013. In particolare, si vede il foil sopravvento

fuori dall’acqua mentre quello sottovento è completamente immerso. pag. 25

Figura 4: Modello completo su XFLR5 pag. 27

Figura 5: Dettaglio della mesh: distribuzione coseno della densità dei pannelli lungo la

corda dell’ala pag. 27

Figura 6: Distribuzione del 𝐶𝑝 adimensionale per Star-ccm+ e per XFLR5 (in rosso Star) pag. 28

Figura 7: Distribuzione della portanza calcolata su XFLR5 pag. 29

Figura 8: Coefficiente di portanza all’aumentare delle iterazioni: da questa e dalla successiva

immagine si può notare come la simulazione effettuata sia andata a convergenza. pag. 30

Figura 9: Andamento dei residui delle equazioni del modello su Star-ccm+ pag. 30

Figura 10: Scena di pressione con tracciamento delle linee isobariche pag. 31

Figura 11: Scena di pressione con tracciamento delle linee isobariche pag. 31

Figura 12: Dettaglio della mesh fluidodinamica pag. 32

Figura 13: Valori di y+ che vanno a determinare il tipo di legge a parete utilizzata nel solutore pag. 32

Figura 14: Pannello sandwich realizzato in fibra di carbonio e core in pvc espanso pag. 35

85
Figura 15: Spostamento in direzione x lungo lo spessore per sandwich simmetrico composto

da pelli in carbonio e core in pvc espanso pag. 37

Figura 16: mesh usata nei test 1, 2, 3 pag. 44

Figura 17: Spostamento in direzione perpendicolare al piano della piastra per L/H=10 pag. 45

Figura 18: Spostamento in direzione perpendicolare al piano della piastra

per L/H=10 (mesh visibile) pag. 47

Figura 19: Spostamento in direzione x per L/H=10 pag. 49

Figura 20: sandwich test 4 pag. 50

Figura 21: Mesh ibrida del test 4 pag. 51

Figura 22: w(x) del test 4 pag. 62

Figura 23: stratigrafia di riferimento per la formulazione delle condizioni al contorno pag. 64

Figura. 24: σx lungo lo spessore pag. 65

Figura 25: σxy lungo lo spessore pag. 65

Figura 26: Spostamento in direzione y lungo lo spessore pag. 65

Figura 27: σyz lungo lo spessore pag. 65

Figura 28: Spostamento verticale lungo lo spessore pag. 66

Figura 29: σxz lungo lo spessore per x =L/2 pag. 70

Figura 30: U(z) per x =L/2 pag. 71

Figura 31: W(x) normalizzato con la freccia massima del FEM mixed 2D/3D pag. 71

Figura 32: Output di XFLR5 pag. 75

Figura 33: definizione per pezzi dell’ala su XFLR5 pag. 78

Figura 34: Flow chart che descrive il ciclo iterativo completo pag. 82

86
15. LISTA SIGLE

AC = America’s Cup

𝑪𝑫 = Coefficient di drag

𝑪𝑳 = Coefficient di lift

𝑪𝑷 = Coefficiente di pressione

DST = discrete shear theory

ESL = Equivalent single layer

FEA = Finite element analysis

FEM = Finite element method

FSDT = first order shear deformation theory

HPC = High performance computing

LLT = Lifting line theory

PLV = Principio dei lavori virtuali

RANS = Reynolds Averaged Navier Stokes

URF = Under relaxation factor

Y+ = distanza adimensionalizzata a parete

VLM = Vortex lattice method

ZZA = Zig zag adaptive theory

87
 16. BIBLIOGRAFIA

1) FLEXIBLE HYDROFOIL OPTIMIZATION FOR THE 35th AMERICA’S CUP WITH CONSTRAINED EGO METHOD
M. Sacher, Naval Academy Research Institute, France, [email protected]
M. Durand, K-Epsilon, Groupama Team France, Sirli Innovations, France
E. Berrini, MyCFD, Universit´e Cˆote d’Azur, INRIA, CNRS, France
F. Hauville, Naval Academy Research Institute, France
R. Duvigneau, Universit´e Cˆote d’Azur, INRIA, CNRS, France
O. Le Maitre, LIMSI - CNRS, France
J. A. Astolfi, Naval Academy Research Institute, France

https://hal.archives-ouvertes.fr/hal-01583591/file/20-Sacher_flexible_foil_optimisation.pdf

2) An experimental investigation of high aspect-ratio rectangular sails

By A. Crook, M. Gerritsen y AND N. N. Mansour – Center forTurbulence Research Annual Research
Briefs 2002

https://web.stanford.edu/group/ctr/ResBriefs02/crook2.pdf

3) A VELOCITY PREDICTION PROCEDURE FOR SAILING YACHTS WITH A HYDRODYNAMIC MODEL BASED ON
INTEGRATED FULLY COUPLED RANSE-FREE-SURFACE SIMULATIONS
Christoph BÖHM

https://pdfs.semanticscholar.org/32b7/0815195513a61ba7df58fcc4e457d1c64bdd.pdf

4) Computational Methods for Hydrofoil Design - A composite analysis using panel code and RANS
Gilberto Besana
University of Southampton - Thesis · January 2015

5) Dynamic Behaviour of a Flexible Yacht Sail Plan.

Benoit Augier, Patrick Bot, Frédéric Hauville, Mathieu Durand.
Ocean Engineering, Elsevier, 2013, 66, pp.32-43.

6) Geometric model for automated multi-objective optimization of foils

Elisa Berrini: National Institute for Research in Computer Science and Control
Bernard Mourrain: National Institute for Research in Computer Science and Control
Régis Duvigneau: National Institute for Research in Computer Science and Control
https://hal.archives-ouvertes.fr/hal-01524287/document

88
7) Numerical simulation of sailing boats: dynamics, FSI, and shape optimization
Matteo Lombardi, Nicola Parolini, Alfo Quarteroni and Gianluigi Rozza

https://infoscience.epfl.ch/record/175879/files/PaerErice-Lombardi-parolini-quarteroni-Rozza.pdf

8) SAIL TRIMMING FSI SIMULATION - COMPARISON OF VISCOUS AND

INVISCID FLOW MODELS TO OPTIMISE UPWIND SAILS TRIM
M. Sacher1, [email protected]
F. Hauville1, P. Bot1, M. Durand2
5th High Performance Yacht Design Conference Auckland, 10-12 March, 2015

https://www.k-epsilon.com/uploads/pdf/papiers/HPYD5_Presentation.pdf

9) Mechanics of laminated composite plates and shells: theory and analysis.

J.N. Reddy CRC Press, Boca Raton, 2003.

http://fn.iust.ac.ir/files/fnst/ssadeghzadeh_52bb7/files/EB__Mechanics_of_Laminated_Composite
_Plates_and_Shells_-JN_Reddy.pdf

10) Development of an efficient zig-zag model with variable representation of

displacements across the thickness.

U. Icardi, F. Sola. J. of Eng. Mech. 140 (2014), 531-541,

https://doi.org/10.1061/(ASCE)EM.1943-7889.0000673.

11) Free and forced vibration of laminated and sandwich plates by zig-zag theories differently accounting
for transverse shear and normal deformability.

U. Icardi, A. Urraci. Aerosp. MDPI. 5 (2018),108,

https://doi.org/10.3390/aerospace5040108

12) Novel HW mixed zig-zag theory accounting for transverse normal

deformability and lower-order counterparts assessed by old and new elastostatic benchmarks.

U. Icardi., A. Urraci. Aer. Sci. &amp;

Tech. 80 (2018), 541-571, https://doi.org/10.1016/j.ast.2018.07.040.

13) Documentazione di Code Aster

https://www.code-aster.org/spip.php?rubrique19

89
14) Documentazione di XFLR5

https://sourceforge.net/projects/xflr5/files/

15) Documentazione di Star-ccm+

https://www.plm.automation.siemens.com/global/it/products/simcenter/STAR-CCM.html

16) Wikipedia: pannello a sandwich

https://it.wikipedia.org/wiki/Pannello_a_sandwich

17) I Kreja, "A literature review on computational models for laminated composite and sandwich panels,"
Central European Journal of Engineering, vol. 1, pp. 59 - 80, 2011.

18) NASA, Vortex-lattice utilization. NASA SP-405, NASA-Langley, Washington, 1976.

https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19760021075.pdf

19) Wikipedia: Teorema dei lavori virtuali

https://it.wikipedia.org/wiki/Teorema_dei_lavori_virtuali

90