Dinamica e modellistica della turbolenza - Esercitazioni del corso
-
Upload
alfredo-patrizi -
Category
Documents
-
view
196 -
download
2
description
Transcript of Dinamica e modellistica della turbolenza - Esercitazioni del corso
Universita degli Studi di Roma Tor Vergata
Facolta di Ingegneria
Corso di Laurea Magistrale in Ingegneria Meccanica
DMT a.a. 2012-2013
Esercitazioni dal corso di Dinamica e
Modellistica della Turbolenza
Giulio Topaziomatr. 0197585
Alfredo Patrizimatr. 0194329
27 settembre 2013
INDICE 2
Indice
1 Esercitazione 1 - Attrattore di Lorenz 51.1 Metodo esplicito: Eulero forward . . . . . . . . . . . . . . . . 5
1.1.1 Simulazioni numeriche e risultati . . . . . . . . . . . . 61.2 Metodo implicito: mid-point . . . . . . . . . . . . . . . . . . . 9
1.2.1 Simulazioni numeriche e risultati . . . . . . . . . . . . 101.3 Dominio di calcolo della soluzione . . . . . . . . . . . . . . . . 131.4 Distribuzione di probabilita . . . . . . . . . . . . . . . . . . . 13
2 Esercitazione 2 - Equazione di Burgers 182.1 Integrazione nel dominio dei numeri d’onda . . . . . . . . . . . 18
2.1.1 Simulazioni numeriche e risultati . . . . . . . . . . . . 182.2 Integrazione nello spazio fisico (x,t) . . . . . . . . . . . . . . . 21
2.2.1 Simulazioni numeriche e risultati . . . . . . . . . . . . 22
Elenco delle figure
1 Soluzione con metodo esplicito, dt = 0.02. . . . . . . . . . . . 62 Coordinata x. . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 Coordinata y. . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 Coordinata z. . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 Attrattore di Lorenz. . . . . . . . . . . . . . . . . . . . . . . . 86 Coordinata x. . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 Coordinata y. . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 Coordinata z. . . . . . . . . . . . . . . . . . . . . . . . . . . . 129 Attrattore di Lorenz. . . . . . . . . . . . . . . . . . . . . . . . 1210 Funzione densita di probabilita. . . . . . . . . . . . . . . . . . 1511 Ordine 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1612 Ordine 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1713 Massimo ordine ammissibile. . . . . . . . . . . . . . . . . . . . 1714 Andamento dei vari modi. . . . . . . . . . . . . . . . . . . . . 1915 Trasferimento di energia ai modi successivi. . . . . . . . . . . . 2016 Dissipazione dell’energia. . . . . . . . . . . . . . . . . . . . . . 2017 Dissipazione dell’energia. . . . . . . . . . . . . . . . . . . . . . 2118 Dissipazione dell’energia: ν = 10−2m2
s. . . . . . . . . . . . . . 23
19 Dissipazione dell’energia: ν = 10−3m2
s. . . . . . . . . . . . . . 23
ELENCO DELLE FIGURE 4
I risultati esposti di seguito sono stati ottenuti attraverso l’implementa-zione degli algoritmi di calcolo necessari in codice C++, programmato uti-lizzando il compilatore Visual C++ 2010 su piattaforma Microsoft Windows(x64).
L’analisi statistica, ove richiesta, e stata condotta sempre attraverso ap-posite funzioni del listato C++. La sola rappresentazione grafica e statarealizzata attraverso software commerciali (Excel, Matlab) o freeware (gnu-plot).
1. Esercitazione 1 - Attrattore di Lorenz 5
1 Esercitazione 1 - Integrazione numerica del-
le equazioni dell’attrattore di Lorenz
L’esercitazione si propone di risolvere, previa integrazione numerica, ilseguente sistema di equazioni:
x = σ(y − x)
y = ρx− y − xzz = −βz + xy
(1)
Successivamente, si richiede di computare i momenti della distribuzionedi probabilita per i valori di x, y, z, fino al generico ordine k.
Si riportano di seguito i risultati ottenuti utilizzando per l’integrazio-ne numerica, rispettivamente, metodo esplicito (Eulero forward) e implicito(midpoint).
1.1 Metodo esplicito: Eulero forward
Il metodo si basa sull’approssimazione della derivata prima della generi-ca funzione u(x) a variabile indipendente x mediante lo sviluppo di Taylor(arrestato al primo ordine) nel punto x0; in particolare, introducendo unadiscretizzazione per x del tipo x = i∆x, vale:
δu
δx|i =
ui − ui−1
∆x+ o(∆x) (2)
Le equazioni dell’attrattore di Lorenz divengono dunque:xi = xi−1 + σ(yi−1 − xi−1)∆t
yi = yi−1 + (ρxi−1 − yi−1 − xi−1zi−1)∆t
zi = zi−1 + (−βzi−1 + xi−1yi−1)∆t
(3)
Segue l’estratto dal listato di codice che implementa la funzione perl’integrazione numerica delle equazioni attraverso il metodo in esame.
void explicitIntegration(double *x, double *y, double *z,double x0, double y0, double z0)
{x[0] = x0;y[0] = y0;
1.1 Metodo esplicito: Eulero forward 6
z[0] = z0;for(int i = 1; i<gI_steps; i++){
x[i] = x[i-1]+ gD_deltaT*gD_sigma*(y[i-1]-x[i-1]);y[i] = y[i-1]+
gD_deltaT*(gD_rho*x[i-1]-y[i-1]-x[i-1]*z[i-1]);z[i] = z[i-1]+ gD_deltaT*(-gD_beta*z[i-1]+x[i-1]*y[i-1]);
}}
Listato 1: La funzione explicitIntegration().
1.1.1 Simulazioni numeriche e risultati
La soluzione ottenuta con metodo esplicito non e incondizionatamen-te stabile; si puo dimostrare infatti che il limite di stabilita si ha per undt = 0.0192625. Si riporta un esempio computato utilizzando il suddettometodo di Eulero “in avanti” con uno step d’integrazione pari a 0.02.
Figura 1: Soluzione con metodo esplicito, dt = 0.02.
A scopo dimostrativo sono state calcolate le soluzioni a partire da con-dizioni iniziali uguali a meno di una variazione di ordine 10−5 su y0: i corri-spondenti grafici mostrano nei valori di x, y e z variazioni di ordine 10.
Si riporta inoltre un plot 3d dell’attrattore di Lorenz integrato con e senzaperturbazione iniziale.
1.2 Metodo implicito: mid-point 9
1.2 Metodo implicito: mid-point
Utilizzando la medesima simbologia adottata nel paragrafo precedente, ilmetodo si basa sull’approssimazione della derivata prima ui nel punto i− 1
2
(punto medio tra i − 1 e i). Il metodo e detto implicito in quanto presumela conoscenza del valore della derivata nel punto di arrivo (i); cio comportal’esigenza di risoluzione per via iterativa, a vantaggio, tuttavia, di un ordinedi accuratezza pari a 2.
δu
δx|i− 1
2=ui − ui−1
∆x+ o(∆x2) (4)
Nel listato, per il generico istante temporale i, ad ogni iterazione k vieneassunto, come stima per ciascuna delle variabili x, y, z, il valore medio trail valore calcolato nel ciclo immediatamente precedente ed il valore iniziale,essendosi preliminarmente posto quest’ultimo pari a quello computato nelpasso temporale i− 1.
xpredict = 0.5(xtemp + xi−1) (5)
In particolare nela precedente e possibile calcolare xtemp, mediante l’in-troduzione di un opportuno fattore di sottorilassamento α, come segue:
xtemp = (1− α)xki + αxk−1i (6)
L’iterazione prosegue fintanto che non si giunge a convergenza secondouna prestabilita soglia di tolleranza.
Il sistema da risolvere numericamente diventa quindi:xi−xi−1
∆t= σ(yp − xp)
yi−yi−1
∆t= ρxp − yp − xpzp
zi−zi−1
∆t= βzp + xpyp
(7)
void midPoint(double *x, double *y, double *z, double x0,double y0, double z0)
{char tempString[1000];double xP, yP, zP, underRelaxationFactor = 0;x[0] = x0;y[0] = y0;z[0] = z0;
1.2 Metodo implicito: mid-point 10
double xT, yT, zT, conv = 0.00000001;for(int i = 1; i<gI_steps; i++){
x[i] = x[i-1];y[i] = y[i-1];z[i] = z[i-1];do {
xT = (1-underRelaxationFactor)*x[i] +underRelaxationFactor*xT;
yT = (1-underRelaxationFactor)*y[i] +underRelaxationFactor*yT;
zT = (1-underRelaxationFactor)*z[i] +underRelaxationFactor*zT;
xP = (x[i-1]+xT)*0.5;yP = (y[i-1]+yT)*0.5;zP = (z[i-1]+zT)*0.5;x[i] = x[i-1] + gD_deltaT*gD_sigma*(yP - xP);y[i] = y[i-1] + gD_deltaT*(gD_rho*xP - yP - xP*zP);z[i] = z[i-1] + gD_deltaT*(-gD_beta*zP + xP*yP);
} while (abs((x[i]-xT)) > conv || abs((y[i]-yT)) > conv ||abs((z[i]-zT)) > conv);
}}
Listato 2: La funzione midPoint().
1.2.1 Simulazioni numeriche e risultati
Analogamente al caso della soluzione calcolata con metodo esplicito, siriportano gli andamenti di x, y, z a partire da condizioni iniziali imperturbatee perturbate. Ancora, a una variazione di ordine 10−5 su y0 corrispondonovariazioni di ordine 10 in x, y, z.
1.3 Dominio di calcolo della soluzione 13
1.3 Dominio di calcolo della soluzione
Sono state eseguite anche diverse runs del listato andando a modificarela sola estensione del dominio di integrazione (fino 500, 5000 e 50000 unitatemporali, con dt = 0.005), a parita di ciascuno degli altri parametri.
Seguono le tabelle con le rispettive medie e deviazioni standard.
Tabella 1: Valori medi e deviazioni standard per t = 500, dt = 0.005
time = 500.000000, dt = 0.005
< x > = -0.430321, SDx = 9.027859< y > = -0.430707, SDy = 10.454852< z > = 30.609365, SDz = 9.685412
Tabella 2: Valori medi e deviazioni standard per t = 5000, dt = 0.005
time = 5000.000000, deltaT = 0.005
< x > = 0.102233, SDx = 9.035895< y > = 0.102108, SDy = 10.451678< z > = 30.613749, SDz = 9.665703
Tabella 3: Valori medi e deviazioni standard per t = 50000, dt = 0.005
time = 50000.000000, dt = 0.005
< x > = -0.052813, SDx = 9.033357< y > = -0.052792, SDy = 10.460466< z > = 30.594805, SDz = 9.686376
1.4 Distribuzione di probabilita
Sia per le soluzioni ottenute con il metodo di Eulero in avanti che perquelle computate implicitamente, e stata condotta un’analisi statistica sulladistribuzione di valori per le coordinate x, y e z, calcolandone le rispettivefunzioni di distribuzione di probabilita ed i momenti teorici centrati nel valore
1.4 Distribuzione di probabilita 14
atteso, fino al generico ordine k, attraverso gli appositi metodi probDensity()e probMomentum().
double *probDensity(double *x, int startValue, int steps, intnRanges)
{double *probFunction = new double[nRanges]();
double max = maximumValue(x, startValue, steps), min =minimumValue(x, startValue, steps);
for(int i = startValue; i < steps; i++){
bool tempCheck = false;for(int k = 0; k < nRanges && !tempCheck; k++)
if(x[i]>=min+k*(max-min)/nRanges && x[i] <min+(k+1)*(max-min)/nRanges)
{probFunction[k] += 1;tempCheck = true;
}}return probFunction;
}double **probMomentum(double *x, double *probX, int startValue,
int steps, int nRanges, int order){
double **probMomentum = new double*[order];double max = maximumValue(x, startValue, steps), min =
minimumValue(x, startValue, steps), avg = average(x,startValue, steps) ;
for(int j = 0; j < order; j++){
probMomentum[j] = new double[nRanges];for(int i = 0; i < steps; i++)
for(int k = 0; k < nRanges; k++)if(x[i]>=min+k*(max-min)/nRanges && x[i] <
min+(k+1)*(max-min)/nRanges)probMomentum[j][k] += pow(x[i]-avg, j);
}return probMomentum;
}
Listato 3: Le funzioni probDensity() e probMomentum().
Segue il grafico dell’andamento della funzione distribuzione di probabilita; ivalori delle ascisse risultano centrati nei rispettivi valori medi.
1.4 Distribuzione di probabilita 15
Figura 10: Funzione densita di probabilita.
Si riportano invece di seguito, ancora a titolo dimostrativo, alcuni gra-fici delle funzioni integrande per il calcolo dei momenti di distribuzione diprobabilita:
� ordine 1: il grafico e esemplificativo dell’andamento della funzione in-tegranda per momenti di ordine dispari; la distribuzione risulta sim-metrica rispetto allo zero ed il momento centrato e dunque nullo perciascuna coordinata;
� ordine 2: il grafico e esemplificativo dell’andamento della funzione inte-granda per momenti di ordine pari; i momenti centrati in x, y, z risul-tano pari alle rispettive varianze, come si deduce confrontando i valoridei primi con i valori delle deviazioni standard calcolate con funzioniappositamente dedicate;
1.4 Distribuzione di probabilita 16
Tabella 4: Valori medi, deviazioni standard e momenti del primo e secondo ordine per t =50000, dt = 0.005
time = 50000.000000, dt = 0.005
< x > = -0.059689, SDx = 9.201369< y > = -0.059824, SDy = 10.552230< z > = 31.485401, SDz = 9.193244M I
x = -0.000020 M IIx = 84.665358
M Iy = -0.000031 M II
y = 111.349276M I
z = -0.000031 M IIz = 84.515324
Figura 11: Ordine 1.
1.4 Distribuzione di probabilita 17
Figura 12: Ordine 2.
Si e infine valutato il massimo ordine di validita della statistica graficandoi valori delle funzioni integrande dei momenti statici di ordini superiori. Ilgrafico di seguito riportato evidenzia che la statistica e da non ritenersi piuvalida a partire dall’ordine 27.
Figura 13: Massimo ordine ammissibile.
2. Esercitazione 2 - Equazione di Burgers 18
2 Esercitazione 2 - Integrazione numerica del-
l’equazione di Burgers
L’esercitazione ha come scopo l’integrazione numerica della seguente equa-zione nota come equazione di Burgers:
δu
δt+ u
δu
δx= ν
δ2u
δx2(8)
La relazione monodimensionale appena esposta e di notevole importanza inquanto presenta le principali caratteristiche delle equazioni di Navier-Stokesfatta eccezione per il termine di pressione.
L’integrazione verra condotta in due modi differenti: dapprima si cercherala soluzione nel dominio dei numeri d’onda, limitando cioe l’∞ ad un numerodiscreto di modi, procedendo poi all’integrazione nello spazio fisico (x,t).
2.1 Integrazione nel dominio dei numeri d’onda
Ipotizzando una opportuna soluzione serie di seni
u(x, t) =∞∑k=1
Ak(t)sin(kx) (9)
si giunge, mediante l’utilizzo delle proprieta di ortogonalita delle funzioniseno, ad una soluzione del tipo:
Ak +∞∑m=1
m(AmAk−m
2+AmAk+m
2− AmAm−k
2) = −νk2Ak (10)
2.1.1 Simulazioni numeriche e risultati
L’equazione (10) e stata integrata limitando il numero di infiniti a 3 modie a 150 modi. Il grafico seguente mostra l’andamento dei suddetti modi neltempo nella circostanza in cui l’energia venisse immessa al solo modo A(1)al tempo iniziale ( tutti gli altri valori di A(i, t0) sono nulli).
2.1 Integrazione nel dominio dei numeri d’onda 19
(a) Infinito a 3 modi.
(b) Infinito a 150 modi.
Figura 14: Andamento dei vari modi.
2.1 Integrazione nel dominio dei numeri d’onda 20
Si riporta anche un grafico che mostra il trasferimento di energia peristanti successivi a i modi caratterizzati da numeri d’onda maggiori (cascatadiretta dell’energia). La soluzione proposta riguarda il caso di “infinito ”a150 modi ma ne vengono riportati, per ovvie ragioni, soltanto 14.
Figura 15: Trasferimento di energia ai modi successivi.
La validazione del metodo e stata condotta andando a computare l’ef-fettiva dissipazione di energia durante il fenomeno. L’energia in questione ecalcolabile dalla somma dei quadrati delle ampiezze dei singoli modi.
Figura 16: Dissipazione dell’energia.
2.2 Integrazione nello spazio fisico (x,t) 21
Viene infine valutata la soluzione u(x, t) serie di seni. Il grafico che segueriporta l’andamento della suddetta a 3 diversi intervalli temporali. Si osservaun’appiattimento della curva attorno al valore x = π, corrispondente allameta del dominio spaziale(D = [0, 2π]). Cio e giustificabile ancora una voltadal meccaniscmo di cascata dell’energia: andando avanti nel tempo l’energiaviene trasferita a strutture sempre piu piccole cui corrispondono gradientisempre piu ripidi.
Figura 17: Dissipazione dell’energia.
2.2 Integrazione nello spazio fisico (x,t)
In questa sezione l’equazione (8) verra risolta direttamente nello spaziofisico utilizzando il metodo di integrazione esplicito di Adams-Bashforth. Ilmetodo qui proposto e un metodo detto multistep poiche, nella scrittura delleequazioni alle differenze che lo caratterizza, sono coinvolti piu valori, calcolatiagli istanti precedenti, della serie che rappresenta la soluzione: per il genericosistema di equazioni differenziali
{y} = {f ({y} , t)} (11)
al passo n + k la soluzione e approssimata attraverso l’uso di combinazionilineari delle funzioni y ed f valutate ai passi n+ j con j = 0, 1, . . . , k dove kindica il numero di passi del metodo.
2.2 Integrazione nello spazio fisico (x,t) 22
Nel caso in esame e stato utilizzato un metodo al secondo ordine la cuiforma generale e la seguente:
yn = yn−1 + h3
2f(xn−1, yn−1)− h1
2f(xn−2, yn−2) (12)
La soluzione dell’equazione di Burgers viene definita nell’intervallo [0, 2π]e si ipotizza che la soluzione sia periodica in x con media nulla. Come nelcaso di integrazione nel dominio dei numeri d’onda la soluzione u(x, t) si puoespandere con una serie di seni del tipo:
u(x, t) =∞∑k=1
Ak(t)sin(kx) (13)
la cui dinamica temporale della soluzione e tenuta in conto dai coefficien-ti Ak(t). Si pongano nulle le soluzioni per x = 0 e x = 2π (condizioni acontorno).
2.2.1 Simulazioni numeriche e risultati
Per la prima integrazione si e costruita una griglia di calcolo spazio-temporale uniforme. Il passo temporale di integrazione fissato e di ∆t =0.001s, mentre il passo di integrazione spaziale e ∆x = 2π
stepscon steps =
100000.Il grafico che segue riporta, come nel caso dell’integrazione dell’equazione
di Burgers nel dominio dei numeri d’onda, l’andamento della suddetta a 5diversi intervalli temporali. Si osserva un’appiattimento della curva attornoal valore x = π, corrispondente alla meta del dominio spaziale(D = [0, 2π]).Cio e giustificabile ancora una volta dal meccaniscmo di cascata dell’energia:andando avanti nel tempo l’energia viene trasferita a strutture sempre piupiccole cui corrispondono gradienti sempre piu ripidi.
Il primo grafico riporta l’andamento della soluzione per valori di viscositacinematica ν = 10−2m2
sevidentemente troppo elevati in quanto si osserva
un forte smorzamento (riduzione dell’ampiezza dell’oscillazione) per tempibrevi.
2.2 Integrazione nello spazio fisico (x,t) 23
Figura 18: Dissipazione dell’energia: ν = 10−2m2
s
Si e poi diminuito il valore di viscosita per osservare una dinamica “piulenta ”. Se ne riporta l’andamento nel grafico che segue.
Figura 19: Dissipazione dell’energia: ν = 10−3m2
s
Infine si e variato, a parita di tutti gli altri parametri, l’intervallo tempo-