Aero_2009_APPENDICI.pdf

C.GOLIA: Aerodinamica 1 Pagina |APP. 1

Aero_APP_a.docx 23/01/2009

Appendici

Scopo

Di solito l’appendice in libri e articoli scientifici, è una parte aggiunta che contiene dati, tabelle ecc.., certamente secondaria rispetto a quella principale. Lo scopo delle appendici di queste note è alquanto differente, più che fornirvi dati (tranne poche) vogliono svolgere un’azione di tutoraggio per aiutarvi a capire chi, cosa, come, perchè e quindi ad aiutarvi a fare, a calcolare, ed a presentare quanto riuscite a produrre, autonomamente, a questo stadio di evoluzione del vostro apprendimento. In conclusione queste appendici non sono secondarie al testo, e certamente sono pariteticamente principali, come o forse più del testo perché dovrebbero darvi quel qualcosa di professionalità che vi rimarrà più a lungo.

Indice Appendice 1 Atmosfera Standard Appendice 2 Progetto AEROFORCES Appendice 3 Profili Alari Appendice 4 Le equazioni del moto di un corpo rigido in un fluido Appendice 5 Codice NACAPROF Appendice 6 Codice Pannelli Sorgenti Appendice 7 Metodo dei Vortici Concentrati Appendice 8 Metodo Pannelli Portanti Appendice 9 Ala Finita Lineare Dritta Appendice 10 Strato Limite Completo Appendice 11 Schema di rapporto tecnico

C.GOLIA: A

Aero_APP_

Aerodinamica

_a.docx

Unità S.I.

Unità Ing

1

.

glesi

Pag

23

Ap

Atmosfer

Ap

Progetto

gina |APP. 2

3/01/2009

ppendice 1

ra Standard

ppendice 2

Aeroforces

2

9

d

2

s



PROGETTO n° 1 Scrivere, ed utilizzare un programma di calcolo scritto in FORTRAN e chiamato AEROFORCES capace di leggere dati (analitici/numerici/sperimentali) delle velocità e dei coefficienti di attrito superficiali per ricavarne i coefficienti aerodinamici risultanti. Il calcolo del sistema di forze che il fluido genera sul corpo (o viceversa) può essere fatto mediante l'integrazione degli sforzi superficiali (forze su superfici unitarie), cioè mediante integrazione superficiale degli sforzi normali (pressione) e degli sforzi tangenziali (viscosi, sforzi di attrito). Ne derivano, per la forza aerodinamica, le seguenti componenti: • nella direzione dell'asse X:

( ) dStnpF xxS

X τ+−= ∫

• nella direzione dell'asse Y:

( ) dStnpF zzS

Z τ+−= ∫

• per il Momento (componente sull'asse Y), Mo, rispetto al bordo di attacco ( x=0 , z=0)

considerato positivo se cabrante (naso in su):

( ) ( ) dStnpxtnpzM zzxxS

o τ+−−τ+−= ∫

dove: n = versore normale, nx ed nz componenti del versore rispetto agli assi x e z, t = versore tangente, tx ed tz componenti del versore rispetto agli assi x e z,

La forza aerodinamica è normalmente decomposta: • nella direzione normale alla velocità asintotica (Portanza): α+α−= cosFsinFL zx • nella direzione della velocità asintotica (Resistenza): α+α= sinFcosFD zx Ovviamente queste forze dipendono dalla forma del corpo, dalla sua posizione rispetto alla corrente e dal regime del campo di moto [che dipende dalla velocità, dimensioni del corpo e dall'ambiente di volo (quota dell’atmosfera che determina p, ρ, ν)]. Si preferisce usare una rappresentazione adimensionale di queste forze e del momento; ne discendono i coefficienti aerodinamici che dipenderanno da altri numeri dimensionali [angolo d’attacco, numeri di Reynolds, Mach…] :

n

τp

X

L

α

Z

V∞

Dt

dSz

x

dS

c = corda

xcp

F

Mo



Coefficiente di Portanza: Sq

LSV

LC 221L

∞=

ρ≡

Coefficiente di Resistenza: Sq

DSV

DC 221D

∞=

ρ≡

Coefficiente di Momento (b.d.a): cSq

McSV

MC o2

21

oM o

∞=

ρ≡

Posizione del Centro di Pressione: CL

Ccx oM

cp −=

In estrema sintesi, per la determinazione dei coefficienti aerodinamici si devono effettuare le seguenti integrazioni che si derivano facilmente dalle formule sopra date:

( ) ( ) ⎥⎦

⎤⎢⎣⎡ +−α+⎥⎦

⎤⎢⎣⎡ +−α−=α+α−=≡ ∫∫

∞∞∞dStcnccosdStcncsincos

SqFsin

SqF

SqLC zfzp

Sxfxp

Szx

L

( ) ( ) ⎥⎦⎤

⎢⎣⎡ +−α+⎥⎦

⎤⎢⎣⎡ +−α=α+α+=≡ ∫∫

∞∞∞dStcncsindStcnccossin

SqFcos

SqF

SqDC zfzp

Sxfxp

Szx

D

( ) ( ) ⎥⎦⎤

⎢⎣⎡ +−−⎥⎦

⎤⎢⎣⎡ +−=≡ ∫∫

∞dStcncxdStcncz

ScqMC zfzp

Sxfxp

So

Mo

Suggerimenti per la scrittura del codice Supponiamo un codice di calcolo in cui i dati vengono letti da un file esterno che fornisce:

• i numeri di punti sul dorso (Nd) e sul ventre (Nv) • e successivamente le ascisse del punto, la velocità e il coefficiente di attrito locale:

do i = 1, Nd xd(i) zd(i) vd(i) cfd(i) ……… enddo do j = 1, Nv xv(i) zv(i) vv(i) cfv(i) ……… enddo

dove: • xd(i), zd(i), vd(i), cfd(i) sono rispettivamente le coordinate x e z del punto e la velocità

ed il coefficiente di attrito sul dorso, • xv(i), zv(i), vv(i), cfv(i) sono rispettivamente le coordinate x e z del punto e la velocità

ed il coefficiente di attrito sul ventre. Ovviamente i coefficienti di pressione sono ricavabili immediatamente dalle velocità superficiali dalle relazioni:

cpd(i) =1- vd(i)2, cpv(i) =1- vv(i)2



Riferendoci alla figura:

è facile riscontrare che valgono le seguenti relazioni:

)i(x)1i(x)i(z)1i(z

xztan;sin;cos;cosn;sinn

zxzx −+−+

≅∂∂

=θθ=τθ=τθ=θ−=

e che inoltre, volendo effettuare gli integrali nella sola variabile "x": dx

xzdxtandx

cossinsindsdz;

cosdxds

∂∂

−=θ−=θθ

−=θ−=θ

=

da cui:

dxtandssinds;ddscosdsn;dxdscosds;dxtandssindsn zzxx θ=θ=τΔ=θ==θ=τθ−=θ−= in questo modo le formule di integrazione diventano:

( )⎥⎥⎦

⎤

⎢⎢⎣

⎡⎟⎟⎠

⎞⎜⎜⎝

⎛∂∂

+∂∂

++−= ∫ ∫c

0

c

0 ventreventre,f

dorsodorso,fventre,pdorso,pF dx

xzC

xzCdxCC

c1C

X

( )⎥⎥⎦

⎤

⎢⎢⎣

⎡++⎟⎟

⎠

⎞⎜⎜⎝

⎛∂∂

−∂∂

= ∫ ∫c

0

c

0ventre,fdorso,f

ventreventre,p

dorsodorso,pF dxCCdx

xzC

xzC

c1C

Z

( ) ( )

( )

⎥⎥⎥⎥⎥⎥

⎦

⎤

⎢⎢⎢⎢⎢⎢

⎣

⎡

⎟⎟⎠

⎞⎜⎜⎝

⎛+

∂∂

−+⎟⎟⎠

⎞⎜⎜⎝

⎛+

∂∂

+⎟⎟⎠

⎞⎜⎜⎝

⎛∂∂

+∂∂

−−

=

=τ+−−τ+−=

∫∫

∫∫

∫

c

0ventreventre,f

ventreventre,pdorsodorso,f

dorsodorso,p

c

0

c

0 ventreventre,f

dorsodorso,fventre,pdorso,p

c

02

zzxxS

M

dxzCxzCdxzC

xzC

dxxxzC

xzCdxxCC

c1

dStnpxtnpzCo

ed quindi:

α+α−= cosCsinCC

zX FFL

x

z

i+1

i

xi xi+1

zi

zi+1

n

nx

nz

θ

ds

x

z

i+1

i

xi xi+1

zi

zi+1τ

τx

τz θ

ds



α+α= sinCcosCCzX FFD

L'integrazione può effettuarsi, in prima battuta, con la semplicissima regola del trapezio. Per ogni stazione:

( ) ( )11ii1i

x

xi xx

2)x(f)x(fdxxfInt

1i

i

−+

≈≡ ++∫

+

l'integrale totale (se N è il numero di punti) sarà quindi:

i

1N

1i

c

0

Intdx)x(f ∑∫−

=

=

! ! ATTENZIONE ! ! Quanto segue è una serie di considerazioni primordiali su come si

scrive un codice di calcolo in FORTRAN. Quanto detto è in forma arcaica e per niente ottimizzata.

Se ne sconsiglia la lettura a quanti hanno esperienza nella programmazione FORTRAN

Facciamo riferimento al FORTRAN in formato fisso (contrapposto al formato libero disponibile negli ultimi compilatori). In questo formato ogni linea di programma è divisa in quattro sezioni.

• Le colonne da 1 a 5 sono usate per labels di riferimento (per GOTO, o per FORMATS). Una C nella prima colonna, farà ignorare tutto quanto segue nella riga; questo fatto sarà quindi usato per inserire commenti.

• Una qualsiasi lettera nella 6.sta colonna collegherà quanto scritto alla riga precedente (i.e. continuazione).

• Gli statements di istruzione devono essere compresi tra la 7ma e la 72.ma colonna. • Tutto quanto compreso dalla 73.ma posizione in poi sarà ignorato.

L'uso di un compilatore moderno faciliterà il riconoscimento delle quattro zone, in quanto colorerà differentemente quanto ivi inserito. Inoltre tutte le funzioni interne del linguaggio (nomi riservati) saranno colorati in un modo diverso dagli statements. Il linguaggio FORTRAN ignora la differenza tra lettere maiuscole e quelle minuscole: per cui le variabili ICAT, icat, Icat, iCat, ecc saranno considerati come simili. E' da dire che il codice che stiamo considerando è, in effetti, molto semplice, ma approfittiamo di questa occasione per illustrare l'uso di particolari tecniche (COMMON, SUBROUTINES, FUNCTIONS, ecc), che saranno usate quindi ad libitum al solo scopo di rendere lo studente

dx)x(f1i

i

x

x∫

+

)xx(2

)x(f)x(fi1i

1ii −+

++

f(x)

x

f(xi+1)

f(xi)

xi+1xi



familiare al loro utilizzo [questo causerà gonfiamento e complicazione, non necessario, del programma]. Un programma FORTRAN inizia con una riga del tipo "PROGRAM XYZ" e finisce con una riga "END" In genere la sua struttura è la seguente: • subito dopo la prima riga occorre mettere le dichiarazioni (singola/doppia precisione, variabili

reali, variabili interi, variabili caratteri, variabili logiche, aree comuni, vettori, matrici, ecc) (capitolo 0)

• dopo mettere gli statement di esecuzione che, in via generale, comprendono una fase di lettura (input) di informazioni (capitolo 1), una fase di elaborazione delle informazioni (capitolo 2), ed una fase di scrittura (output) dei risultati (capitolo3).

Anche se non necessario, nello specifico caso usiamo delle subroutine per i vari capitoli. 1234567890123456789012345678901234567890123456789012345678900123456789012

program AEROFORCES c-----capitolo 0: dichiarazioni inserire qui le vostre dichiarazioni …………………………… ………………………… c------capitolo 1 call reading c------capitolo 2 call calcola c------capitolo 3 call printout c------fine stop end Capitolo.0...DICHIARAZIONI Una variabile se non dichiarata è IMPLICITAMENTE considerata dal FORTRAN come

• intera se inizia con una lettera che va da I ad N, • reale, a singola precisione, se la variabile inizia con una lettera che va da A ad H ovvero da

O a Z. Per cui la variabile Icat sarà considerata con intera, la Ciat come reale a precisione singola. Se si volesse la Icat come reale occorrerà inserire lo statement: REAL Icat analogamente se si volesse la variabile Ciat come intera occorrerà inserire lo statement: INTEGER Ciat Vettori e matrici devono essere specificati nelle DIMENSION. Nello specifico preferiamo specificarli nelle aree COMMON etichettate. Nota abbiamo separato le COMMON di variabili intere da quelle contenenti variabili reali.



Questo non è necessario se si usa precisione singola, ma è consigliabile/necessario se si usa precisione doppia.

NOTA: Nella fattispecie, anche se non ve ne è bisogno, si usa precisione doppia in modo implicito (cosa categoricamente sconsigliata da tutti i manuali di programmazione); l’uso della doppia precisione è invero raccomandato qualora si devono fare calcoli di tipo massivo (i.e. se, in genere, si devono risolvere sistemi con dimensioni maggiori di 50 ovvero si devono usare algoritmi iterativi pesanti).

In conclusione il capitolo 0 potrebbe essere come segue: 12345678901234567890123456789012345678901234567890123456789012 c-----capitolo 0: dichiarazioni implicit double precision (a-h,o-z) common /main1/ nd,nv,alfa

common /main2/ xd(401),zd(401),vd(401),cpd(401),cfd(401) & xv(401),zv(404),vv(401),cpv(401),cfv(401),cf(401) Capitolo.1...LETTURA DATI Nello specifico i dati saranno letti da un file esterno. A tale scopo si usa una SUBROUTINE chiamata READING, che comunica alle altre subroutine, tutti i dati tramite le aree COMMON. Da notare che occorre mantenere la doppia precisione implicita, ed occorre specificare che le variabili FILEIN e FILEHEAD sono delle stringhe di caratteri. Nel seguito il nome del file da leggere viene richiesto su video (unità *), e si dovrà immettere da tastiera (unità *). Dopo di che si apre una unità di lettura (nella fattispecie 1), si legge quanto dovuto, dopo di che la si chiude. Nota si usa una funzione PAUSE che blocca la esecuzione, questa continuerà dopo che si preme il tasto della tastiera indicato sullo schermo. L'uso della pause può essere importante nella fase di debugging del programma, nella quale si fanno delle stampe per verificare la correttezza dell'esecuzione di casi extra-semplici. In conclusione il capitolo 1 dovrebbe contenere una SUBROUTINE come quella che segue SUBROUTINE READING implicit double precision (a-h,o-z) character*64 filein, filehead common /main1/ nd,nv,alfa

common /main2/ xd(401),zd(401),vd(401),cpd(401),cfd(401), & xv(401),zv(404),vv(401),cpv(401),cfv(401),cp(401),cf(401) common/filenames/ filein,filehead c.....leggi dati da file esterno print*,' Inputta percorso e nome file di uscita' read*, filein open(1,file=filein,status='unknown') c.....leggi/stampa titolo del file read (1,'(a)') filehead print*, ' ' write(*,'(a)') filehead



c.....leggi numeri di punti (sul dorso e sul ventre) et angolo alfa READ(1,*) np,nv, alfa c print*,np,nv,alfar,alfa c.....leggi dati sul dorso do i=1,np read(1,*) xd(i),zd(i),vd(i),cfd(i) print*,i,xd(i),zd(i),vd(i),cfd(i) end do print*,' sopra i dati ventre dorso' pause c.....leggi dati sul ventre do i=1,nv read(1,*) xv(i),zv(i),vv(i),cfv(i) print*,i,xv(i),zv(i),vv(i),cfv(i) end do print*,' sopra i dati sul ventre' close(1) RETURN END c================fine capitolo di lettura dati======================== Capitolo 2: questa fase contiene le fasi di calcolo. Nella fattispecie si è fatta la scelta (pessima!) di mettere tutti i calcoli in una SUBROUTINE SUBROUTINE CALCOLA implicit double precision (a-h,o-z) common /main1/ nd,nv,alfar,alfa

common /main2/ xd(401),zd(401),vd(401),cpd(401),cfd(401), & xv(401),zv(404),vv(401),cpv(401),cfv(401),cp(401),cf(401) common/prnt/CL,CD,CM c pig = DATAN(1.d0)*4. Alfar = alfa*pig/180 print*, 'angolo di attacco (in radiante)=', alfar print*, 'angolo di attacco (in gradi )=', alfa pause c c....converti velocitat in cp: do i=1,nd cpd(i)= 1.-vd(i)*vd(i) end do c do i=1,nv cpv(i)= 1.-vv(i)*vv(i) end do c xo= xd(1)



zo= zd(1) corda= xd(np) - xd(1) c.....inizializza Fzd= 0. Fxd= 0. CMod= 0. Fzv= 0. Fxv= 0. CMov= 0. c:::::::::::::::::::::DORSO::::::::::::::::::: print*,' valori medi i dati dorso: dzdx,x,z,cp,cf' do i=1,nd-1 dxd= xd(i+1)-xd(i) dzd= zd(i+1)-zd(i) xmd= DMEDIA(xd(i+1),x(i)) zmd= DMEDIA(zd(i+1),z(i)) cpdm= DMEDIA (cpd(i+1),cpd(i)) cfdm= DMEDIA (cfd(i+1),cfd(i)) dzdxd= dzd/dxd c Fzd= Fzd +(-cpdm + cfdm*dzdxz)*dxd/corda

Fxd= Fxd +(+cpdm*dzdxd + cfdm)*dxd/corda CModp= CModp + (cpdm*(xmd + zmd*dzdxd)+ & cfdm*(zmd - xmd*dzdxd))*dxd/corda/corda end do c pause C::::::::::::::::::::VENTRE::::::::::::::: do i=1,nv-1 dxv= xv(i+1)-xv(i) dzv= zv(i+1)-zv(i) xmv= DMEDIA(xv(i),xv(i+1)) zmv= DMEDIA(zv(i),zv(i+1)) cpvm= DMEDIA(cpv(i+1),cpv(i)) cfvm= DMEDIA(cfv(i+1),cfv(i)) dzdxv= dzv/dxv c Fzv= Fzv + (cpvm+ cfvm*dzdxv)*dxv/corda c Fxv= Fxv +(- cpvm*dzdxv + cfvm)*dxv/corda c CMovp= CMovp – (cpvm*(xmv+zmv*dzdxv)+ & cfvm*(zmv-xmv*dzdxv))dxv/corda/corda

end do c------componi le forze Fx= Fxd+Fxv Fz= Fzd+Fzv CL= -Fx*sin(alfar)+Fz*cos(alfar) CD= Fx*cos(alfar)+Fz*sin(alfar)



c RETURN END c============fine calcoli================== Capitolo 3: Outputs L'output di stampa [su video (*) e su file (5)] viene realizzato (pessimamente) mediante questa subroutine: SUBROUTINE PRINTOUT

implicit double precision (a-h,o-z) character*64 filein, filehead, fileout common/filenames/ filein,filehead common/prnt/CL,CD,CM c print*,' vuoi output su file ? (0=NO) (1=SI)' read*, intprint

if(intprint .eq. 1) then c print*,' Inputta percorso e nome file di uscita dati es. dati.txt' read*, fileout open(5,file=fileout,status='unknown') c write(5,'(a)') 'file di lettura: ','filein' write(5,'(a)') 'intestazione file: ','filehead' write(5,1000) np,nv,alfa 1000 FORMAT(1x,'np = ',I3,1x,' nv = ',I3,2x,' alfa = ',f7.3) write(5,*)' ' write(5,*)' coefficiente di portanza =', CL write(5,*)' ' write(5,*)'coeff.resist. =', CD,CDd,CDv write(5,*)' ' write(5,*)'coeff.mom.(bda)=', CMo,CMod,CMov write(5,*)' ' write(5,*)'xcp/c= ', CMo/CL

close(5) c end if c write(*,'(a)') 'file di lettura: ','filein' write(*,'(a)') 'intestazione file: ','filehead' write(*,1000) np,nv,alfa write(*,*)' ' write(*,*)' coefficiente di portanza =', CL write(*,*)' ' write(*,*)'coeff.resist. =', CD,CDd,CDv write(*,*)' ' write(*,*)'coeff.mom.(bda)=', CMo,CMod,CMov



write(*,*)' ' write(*,*)'xcp/c= ', CMo/CL RETURN END c=======fine capitolo 3================== Per ragioni squallidamente dimostrative si fa (squallido) uso di una funzione esterna che fa la media di due costanti: c========== funzione MEDIA========= FUNCTION DMEDIA(f1,f2)

implicit double precision (a-h,o-z) DMEDIA= (f1+f2)/2.d0 RETURN END c==== fine funzione MEDIA========

Ne approfittiamo per chiederci: perché la funzione è dichiarata DMEDIA invece che MEDIA?

ATTENZIONE

Quanto riportato in grassetto può essere copiato in un listato di programma FORTRAN con copia/incolla. Ma bada che: • il FORTRAN accetta file text e non doc! • vi sono nascosti degli errori (alcuni diabolicamente voluti altri sciaguratamente presenti) che il

vostro compilatore dovrebbe rilevarvi, • così avrete la possibilità di esercitarVi al debugging del codice, i.e. a trovare gli errori ed a

correggerli (niente di più istruttivo il capire gli errori che si fanno!)

Buon lavoro

Dopo che il programma funziona ricordate di:

• testarlo • riscriverlo in modo semplice e pulito, eliminando tutte le astrusità messe soltanto per

dimostrazione

Il testing (validazione) di un codice è fondamentale per il suo uso professionale. Per il testing si usano diverse tecniche (alternative e/o consecutive):

• run di casi semplici per cui sono evidenti soluzioni elementari algebriche, che vengono concepiti in modo da isolare logicamente alcune particolarità (geometriche ovvero qualitative)

• re-run di casi complicati di cui si ha soluzione certa, per confrontarne i dati. Nel primo caso potrebbe essere “di utilità” considerare un file dati molto semplice a 5+5 punti, facilmente modificabile, come quello che segue:



file: aeroforces_valid_01.dat file di validazione del codice aeroforces : lastra piana CL=1, CD=2 ? 5 5 0. 0. 0. 1. 1. 0.25 0. 1. 1. 0.5 0. 1. 1. 0.75 0. 1. 1. 1.0 0. 1. 1. 0. 0. 0. 1. 0.25 0. 0. 1. 0.5 0. 0. 1. 0.75 0. 0. 1. 1.0 0. 0. 1. Secondo Voi, • questo file cosa rappresenta? • usando il cervello è corretto prevedere CL=1, CD=2? • Provate ad usare il vostro codice per vedere cosa esce per CL e CD! • Così avrete fatto un primo passo di validazione! Sbizzarritevi a creare, in modo intelligente, altri casi che significativamente possono testare tutte le particolarità del vostro codice. Nel secondo caso potrete modificare il file di cui sopra per simulare il caso di un cilindro con circolazione la cui velocità periferica sia [cap.2.3.5.1]:

Vθ= 2 sin(θ) + K/2 , cf=0. Per cui si dovrebbe verificare:

CL=2 π k , CD=0 Provate a fare questa validazione!

D=1

x

y

θ

V∞=1

K



Appendice 3

Profili Alari NACA-0006 NACA-0009 NACA-0012 NACA-2412 NACA-4412 NACA-23012 NACA-651-012 GA(W)-1

Aero

o_APP_a.docx

C.GOLIA: Ae

Naca

erodinamica 1

0006

23/01/2009

Pagina

9

a |APP. 16

Aero

o_APP_a.docx

C.GOLIA: Ae

Naca

erodinamica 1

a 0009

23/01/2009

Pagina

9

a |APP. 18

Aero

o_APP_a.docx

C.GOLIA: Ae

Naca 0

erodinamica 1

0012

23/01/2009

Pagina

9

a |APP. 19

Aero

o_APP_a.docx

C.GOLIA: Ae

Nac

erodinamica 1

ca 2412

23/01/2009

Pagina

9

a |APP. 20

Aero

o_APP_a.docx

C.GOLIA: Ae

Naca 4

erodinamica 1

4412

23/01/2009

Pagina

9

a |APP. 21

Aero

o_APP_a.docx

C.GOLIA: Ae

Naca

erodinamica 1

a 23012

23/01/2009

Pagina

9

a |APP. 22

Aero

o_APP_a.docx

C.GOLIA: Ae

Nac

erodinamica 1

ca 651-012

23/01/2009

Pagina

9

a |APP. 24

Aero

o_APP_a.docx

C.GOLIA: Aeerodinamica 1

GA(W)-1

23/01/2009

Pagina

9

a |APP. 25



Appendice 4

Le equazioni del moto di corpi rigidi in un fluido

Consideriamo un corpo in moto in un fluido viscoso, e che la estensione della sua traiettoria risulti trascurabile rispetto al raggio della Terra. In tale caso si potrà assumere la terra piatta e che quindi il vettore gravità sia costante (non cambia intensità e direzione) agendo sempre secondo l’asse verticale (in un riferimento cartesiano). Per semplificazioni didattiche considereremo campi di moto piani, e usiamo un riferimento cartesiano con gli assi x e z nel piano del moto. Non considereremo quindi le componenti (delle forze, delle accelerazioni, delle velocità e degli spostamenti) nella direzione (y) normale a tale piano. Ovviamente quanto detto sarà facilmente estensibile a moti non planari, con opportune implementazioni che non altereranno però i concetti fondamenti per le equazioni alla traslazione, ma che richiederanno una dinamica più complicata per le equazioni alla rotazione. Considereremo inoltre corpi rigidi, in tal caso le equazioni del moto del baricentro del corpo e quelle del moto attorno al baricentro (rotazione) saranno indipendenti. Nel caso di corpi a massa variabile (ma configurazione rigida) faremo ulteriormente l’ipotesi semplificativa che il momento d’inerzia ICG attorno al baricentro (CG) sia rappresentabile mediante un raggio di girazione costante . E’ abbastanza ovvio che conviene porre l’origine degli assi nel baricentro del corpo. Per l’orientazione degli assi (fermo restando che l’asse delle y è costante e normale (entrante) al piano delle figure che elaboreremo, abbiamo delle scelte (tutti riferimenti orto-normali):

Assi vento (s,n) Assi terra (x,z) Assi corpo (ξ,ζ) In realtà se non si considera il moto attorno al baricentro (dinamica alla rotazione, i.e. solo il moto del punto massa) l’una scelta equivale alle altre. Nell’altro caso l’uso degli assi corpo porta a semplificazioni, anche notevoli, in quanto è abbastanza ovvio che è preferibile scegliere gli assi principali d’inerzia, rispetto ai quali il tensore d’inerzia è diagonale. Gli angoli sono:

α= angolo d’attacco; θ= angolo di assetto; σ=angolo di spinta

MCG è il momento delle forze aerodinamiche rispetto al baricentro LE EQUAZIONI DEL MOTO RISPETTO AGLI ASSI VENTO [VARIABILI INDIPENDENTI V, θ]:

( ) variabilem costante,mI 2CG =

x

zn

s

ζ ξ

θ

θ

θ

σ

α

α

α

l

m g

V

L

D

S

MCG

B



lungo la traiettoria (asse “s”): ( ) ( ) θ−−−α+σ=+ sinBgmDcosS

dtdV)mm( app

lungo la normale alla traiettoria (asse “n”): ( ) ( ) θ−−+α+σ=

θ+ cosBgmLsinS

dtdV)mm( app

equilibrio alla rotazione (nel piano del moto): CGrotapp

2 MsinSldtd)mm(

dtd

+σ=⎥⎦⎤

⎢⎣⎡ θ

+

equazione dello spazio percorso: V

dtds

=

equazione della traiettoria: θ= cosVdtdx ; θ= sinV

dtdz

NOTA: sviluppando le derivate temporali dell'equazione alla rotazione si ricava:

CG

rotapp2

2

2rotapp

2 MsinSldtd

dt)mm(d

dtd)mm( +σ=⎥⎦

⎤⎢⎣⎡ θ+

+θ

+

compare così un termine proporzionale alla velocità di rotazione ed alla velocità di variazione della massa (termine negativo, pari alla portata di propellente) che genera uno smorzamento alla rotazione (fenomeno del Jet damping)

LE EQUAZIONI DEL MOTO RISPETTO AGLI ASSI TERRA [VARIABILI INDIPENDENTI U,W: U=VCOSθ , W=VSINθ ] SONO: lungo la orizzontale (asse “x”): ( ) θ−θ−α+θ+σ=+ sinLcosDcosS

dtdu)mm( app

lungo la verticale (asse “z”): ( ) ( )BmgcosLsinDsinSdtdw)mm( app −−θ+θ−α+θ+σ=+

angolo di assetto: ⎥⎦⎤

⎢⎣⎡=θ −

uwtan 1

equilibrio alla rotazione: ( ) CGrotapp

2 MsinSldtdmm

dtd

+σ=⎥⎦⎤

⎢⎣⎡ θ

+

equazione dello spazio percorso: 22 wuVdtds

+==

equazione della traiettoria: θ== cosVudtdx ; θ== sinVw

dtdz

LE EQUAZIONI DEL MOTO RISPETTO AGLI ASSI CORPO (ξ,ζ) [VARIABILI INDIPENDENTI V,θ]: lungo l’asse “ξ”): ( ) [ ] ( ) ( )α+θ−−α+α−σ=

α+ sinBgmsinLcosDcosS

dtcosVdmm app

lungo l’asse “ζ”): ( ) [ ] ( ) ( )α+θ−−α+α+σ=α

+ cosBgmcosLsinDsinSdt

Vsindmm app

equilibrio alla rotazione: ( ) CGrotapp

2 MsinSldtdmm

dtd

+σ=⎥⎦⎤

⎢⎣⎡ θ

+



equazione dello spazio percorso: Vdtds

=

equazione della traiettoria: θ= cosVdtdx ; θ= Vsin

dtdz

EQUAZIONE DELL'ENERGIA (TERMICA DEL CORPO)

( ) VDTSTThSdtdmcT

dtdTcm 4

oirradiantecorpoo

convettivacorpo

.equvscaricoscarico

.equvcorpo +εσ+−+−=

dove: T = temperatura media del corpo h = coefficiente globale di scambio termico (funzione della forma, del regime di moto e quindi

dei numeri di Re, Gr, Mach) cequiv = calori specifici equivalenti, rispettivamente per il corpo e per i gas di scarico σo = 5.67 108 [W/(m2K4)] costante di Stefan-Boltzman ε = coefficiente di emissione del corpo (=1 per corpo nero) D = resistenza aerodinamica V = modulo della velocità

Per l’integrazione del sistema di equazioni differenziali occorre, ovviamente, fornire le equazioni ausiliarie:

Portanza: L=0.5 ρ V2 A CL Resistenza: L=0.5 ρ V2 A CD Momento aerodinamico: MCG=0.5 ρ V2 A L CMCG

• Forma di galleggiamento: B= ρfluido Volumecorpo ggravità Massa apparente alla traslazione: mapp = K ρfluido Volumecorpo

(K costante che dipende dalla forma e dall'asse di traslazione) Massa apparente alla rotazione: mrot

app = Krot ρfluido Volumecorpo (Krot costante che dipende dalla forma e dall'asse di rotazione)

Nota: le masse apparenti e la forza di galleggiamento sono in genere trascurabili se la densità del corpo è molto maggiore di quella del fluido in cui si muove I valori dei coefficienti aerodinamici:

CL, CD e CM = funzioni di (α, Re, Mach, eventualmente Fr) Le leggi di controllo (variazioni con il tempo):

α=α(t) ; S=S(t) ; m=m(t)

Nota che nel caso di propulsione a razzo si verifica:



I parametri termodinamici del fluido (ambiente):

ρ=ρ(z) p=p(z) T=T(z)

Risoluzione delle equazioni Una volta assegnati i valori iniziali le equazioni possono essere agevolmente integrate numericamente con una qualsiasi routine (ad es. RungeKutta del 4° ordine).

Esercizi che consentono una soluzione analitica Esercizio A.4. 1 Un miscelatore è composto da due coppe (semi-sfere) di diametro 10 cm, poste ad una distanza di 0.6 m dall'asse di rotazione. Stimare la potenza necessaria per miscelare una sostanza acquosa, con una velocità di rotazione di 60 rpm Esercizio A. 4. 2 Un miscelatore è composto da due semi-tubi di diametro 10 cm e lunghi 30 cm, posti ad una distanza di 0.6 m dall'asse di rotazione. Stimare la potenza necessaria per miscelare una sostanza acquosa, con una velocità di rotazione di 60 rpm Esercizio A. 4. 3 Si vuole dimensionare un paracadute capace di rallentare un corpo del peso di 120 kg ad una velocità di discesa terminale di 5 m/s. Determinare il diametro in condizioni di atmosfera standard. Esercizio A. 4. 4 Un viscosimetro è basato sul fatto che il coefficiente di resistenza di una sfera, per Re<1 , è CD=24/ReD.

( )

dtm-mmsarà caso talein

scarico divelocità c

solida) epropulsion in costante solito di tempo,del funzione (possibile scarico di gasportata m

ppAcmS

t

0piniziale

jet

p

jetexitjetp

∫•

•

•

=

=

=

−+=



Una sfera di acciaio (ρsfera = 7800 kg/m3) con diametro D = 6 mm, lasciata cadere in un olio lubrificante (olio SAE 30: ρolio = 933 kg/m3), ha una velocità terminale di 70 cm/s ad una certa temperatura. • Determinare la viscosità dinamica (μ) dell'olio a quella temperatura. • Verificare se l'ipotesi ReD<1 è verificata. • SE l'ipotesi non fosse verificata cosa fareste per misuare la μ? Esercizio A. 4. 5 Un anemometro è composto da due coppe (semi-sfere) aventi diametro di 5 cm, poste a 15 cm di distanza da un asse di rotazione collegato ad una generatore elettrico a tensione costante (12 Volts) che opera con un rendimento di 0.85. Calcolare la curva caratteristica V-Ampere Esercizio A.4. 6 Un fumaiolo di forma cilindrica , diametro 1 m , altezza 25 m, è investito da un vento uniforme di 50 km/h (assenza di turbolenza e raffiche). Assumendo condizioni di atmosfera standard e trascurando gli effetti di bordo, determinare: • La forza resistente, • Il momento flettente alla base, • La frequenza (rad/s) delle vibrazioni indotte dal vento sul fumaiolo. Esercizio A. 4. 7 Su di un'isola deserta un discendente di R.Crusoe, con i potenti mezzi a sua disposizione, costruisce un rudimentale generatore eolico per alimentare l'impianto stereo che è riuscito a salvare dal naufragio della sua nave. Ha trovato un cassa rettangolare (lato 0.80 m, H=1.20 m), ed ha pensato di tagliarla diagonalmente, e di accoppiare, con un moltiplicatore, i due pezzi ad un asse ruotante verticale posto sull'alternatore che ha recuperato dal motore della scialuppa di salvataggio (ovviamente naufragata). Assumendo un rendimento globale meccanico pari a 0.8, quanto deve essere la velocità del vento per poter alimentare lo stereo al massimo (150 w) della potenza per ascoltare il suo CD preferito, fortunatamente trovato sulla spiaggia? In tali condizioni a che velocità gira la macchina? Quale deve essere il rapporto di moltiplicazione (stabilizzato) se l'alternatore (a 12Volts) ha una caratteristica Ampere=.010 x Nrpm Esercizio A. 4. 8 Un aereo ultra leggero acquistato in USA, ha le seguenti caratteristiche: massa a vuoto =3000 lb , superficie alare 300 ft2 , velocità di decollo 100 ft/s , potenza installata 45 Hp Polare sperimentale: CL= 0.35 (0.2 α +1) CD=0.008 (α2+1) (α in gradi) Con due persone a bordo (160 kg) e pieno di carburante (30 kg) l'aereo riesce a decollare? Quanto vale la potenza minima al decollo?



Esercizio A. 4. 9 La vostra automobile ha un peso a vuoto di 800 kg, un'area frontale di 1.8 m2 , un coefficiente di resistenza (finestrini e prese d'aria chiusi) Cx=0.5 ed un coefficiente di attrito al rotolamento di 0.015 [(forza resistente)/(forza peso)] La potenza massima del motore è 110 HP, il consumo (medio) 300 gr/(HP ora), rendimento di trasmissione 0.85. In condizioni standard ( 2 passeggeri + 40 litri di combustibile), e senza aria condizionata: • determinare la velocità massima (supponendo che il rapporto al ponte è ottimale) • determinare il consumo di carburante (litri/100km) per la crociera (in piano) , a 70 km/h, a 130

km/h ed alla velocità massima • quale è la velocità per cui la resistenza aerodinamica è pari a quella di attrito al rotolamento Quindi: • ripetere l'analisi in condizione di massimo carico (5 x (passeggeri +15 kg di bagaglio), pieno di

carburante=80 litri) Quindi bis: • ripetere le due analisi su di un falsopiano in salita

10%. Esercizio A. 4. 10 Progetto preliminare delle superfici alari di un aliscafo. Dati di progetto: Velocità di splash off 15 nodi velocità di crociera = 40 nodi. peso totale 600 tonnellate larghezza 7 metri Determinare le potenze necessarie in condizioni di splash-off e quelle di crociera, determinare le superfici aerodinamiche, i profili, i calettamenti e discutere i criteri per la scelta del/dei flap di controllo Esercizio A. 4. 11 Una schiera di palette di un compressore assiale usa profili NACA 66-212 con una corda di 10 cm, ed un passo pari a 5 cm. Le palette sono poste a σ=20° e la schiera è investita da una corrente con velocità V=400 m/s, α1=30°. Assumendo Rrif=0.8 m e che tali condizioni si realizzano con una velocità angolare ω=250 rad/s, determinare il momento torcente, la potenza e la spinta assiale sulla schiera usando la polare del profilo.

12 8

600 ton.



Esercizi che richiedono una soluzione numerica Esercizio A. 4. 12 Un corpo di forma sferica, densità 8000 kg/m3, volume 0.1 m3, area frontale 0.25 m2, CD=1.1 è accelerato in acqua ad una velocità iniziale di 13 m/s. Determina la traiettoria. Esercizio A. 4. 13 Un rais di un paese sottosviluppato vuole acquistare un cannone gasdinamico capace di sparare proiettili di 30 cm di diametro e peso fino a 200 kg ad una velocità iniziale di Mach = 4, [ ]RT4V.e.i o γ= Supponendo che il peso del proiettile sia 150 kg , che non vi sia portanza, e che il coefficiente di resistenza del proiettile vari con Mach secondo la formula:

⎥⎥⎦

⎤

⎢⎢⎣

⎡ +−+=

86MexpM25004.0C 2

D

determinare le gittate e la massima quota (S ed zmax) al variare dell'angolo di alzata θo.

Assumere per semplicità atmosfera isoterma [T=288 k], per cui discende che ⎥⎦⎤

⎢⎣⎡=

ρρ

=Hhexp

pp

o0,

dove h è l'altezza dal suolo (metri) e per l'atmosfera terrestre H=7050 (metri). Esercizio A. 4. 14 Un aereo d'attacco al suolo vola orizzontalmente ad una quota 10000 m. con una velocità di 800 km/h, e sgancia una bomba sferica da 300 Kg avente un diametro di 20 cm, ed un CD=0.08. Assumendo atmosfera isoterma ed assenza di vento, determinare la traiettoria della bomba (i.e. dove va a sbattere relativamente al punto di sgancio?) Esercizio A. 4. 15 Il solito rais del problema precedente si è esaltato e vuole costruire un razzo ICBM per festeggiare il Capodanno lanciando su un'isola lontana una bomba di carta (sic!). Il razzo pesa alla partenza, 2 tonnellate il cui 90 % è costituito da propellente. Il propellente solido usa una camera di combustione ed un ugello che realizza una velocità di scarico di cjet=3000 m/s,. Memo: trascurando la contro-pressione allo scarico, la spinta del motore a razzo è S = - cjet (dm/dt). [di solito in un razzo a propellente solido dm/dt è costante dall'accensione allo spegnimento]

θo

Vo

x

z

s

zmax

σ

mg

S

V

θ

α

D

L

θo

y

x

traiettoria



Data la rudimentalità della tecnologia disponibile, il razzo segue una traiettoria a gravità, cioè una traiettoria per cui sia portanza nulla e la spinta allineata con la direzione della velocità. Per questo caso le equazioni generali del moto di possono semplificare: nella direzione della velocità: ( ) θ−−= singMDS

dtdVM

nella direzione normale alla traiettoria: θ−=

θ cosMgdtdMV

Equazioni della traiettoria: θ= Vsin

dtdy θ= cosV

dtdx

Trascurando l'inerzia alla rotazione, determinare le traiettorie y(t), x(t) , y(x) per diversi valori dell'angolo di lancio θo=80°,60°,40°. Assumere che • il coefficiente di resistenza del missile vari con Mach secondo la formula:

⎥⎥⎦

⎤

⎢⎢⎣

⎡ +−+=

86MexpM25004.0C 2

D

• atmosfera isoterma [T=288 k], per cui discende che ⎥⎦⎤

⎢⎣⎡=

ρρ

=Hhexp

pp

o0, dove h è l'altezza dal

suolo (metri) e per l'atmosfera terrestre H=7050 (metri).

P.S. non seguite l'esempio e non sparate i botti né a Capodanno, né quando vince il Napoli.

Valori del CD per vari corpi BIDIMENSIONALI semplici (profondità unitaria)

forma del Corpo Rapporti di forma CD Numero di

Reynolds

circolare: 1.2 104 → Cilindro 105 Cilindro ellittico: a:b=2:1 0.6 4 104 0.46 105 4:1 0.32 2.5 104 → 105

a

b



8:1 0.29 2.5 104 0.20 2.5 104 1:2 1.8 < 5 105 0.5 >1 106

Forme aerodinamiche: c/t=1 0.35 107 1.5 0.12 107 2 0.06 107 3 0.04 107 4 0.04 107

5 0.04 107

6 0.045 107 7 0.048 107

8 0.052 107 Cilindri quadrati: 2.0 3.5 104

1.6 3.5 104 Cilindri triangolari: 2.0 104

1.72 104 2.15 104 1.60 104

120

120

60°

90°

90°

c

t



2.20 104 1.39 104 1.80 104 1.0 104 Cilindri semi- tubolari 2.3 104 1.12 104

Valori del CD per corpi TRI-DIMENSIONALI semplici

forma del Corpo Rapporti di forma CD Numero di Reynolds

Cubo: 1.1 104 Coni: 60° 0.5 104 Emisfere: 1.7 104

60°

30°

30°



0.4 104 Piastre rettangolari: b/h=1 1.15 104 5 1.2 104 10 1.3 104 20 1.5 104 Disco: 1.28 104 Cilindri allineati: h/D=0.5 1.15 104 1 0.9 104 2 0.85 104 4 0.87 104 8 0.99 104

Fusi: b/h=1 0.070 107 1.5 0.032 107 2 0.030 107 3 0.033 107 4 0.041 107

5 0.050 107

6 0.060 107 7 0.068 107

8 0.078 107 Listato di un codice in FORTRAN che integra un sistema di 4 equazioni differenziali del primo ordine con un metodo Runge Kutta 4.to ordine PROGRAM PROVARK01 C demo di un programma di prova tipico per un'integrazione C alla Runge Kutta del 4.to ordine DIMENSION Y(15),F(15) c print*,'******************************************************' print*,'* Fludodinamica *' print*,'* programma di prova per integrazione ODE *' print*,'* RUNGE KUTTA 4.o ordine *' print*,'* alunno Mormile Luigi a.a. 2000/2001 *'

b

h

h

D

L

D



print*,'******************************************************' C condizioni iniziali X=0.0 XLIM=30.0 H=0.1 ISTEP=20 ICOUNT=0 M=0 N=4 Y(1)=0.0 Y(2)=0.0 Y(3)=100.0 Y(4)=100.0 open(1,file='Dati.dat') print*,' Tempo x z dx/dt dz/dt ' WRITE(1,1000) X,(Y(I),I=1,N) WRITE(*,1000) X,(Y(I),I=1,N) C inizia l'integrazione C verifica se est superato il limite 8 IF (X-XLIM) 6,6,7 C si incrementa il contatore 6 ICOUNT=ICOUNT+1 9 CALL RUNGE(N,Y,F,X,H,M,K) GO TO (10,20),K C K=1 calcola il valore delle derivate 10 CALL FNC(N,Y,F,X) C ritorna a RUNGE GO TO 9 C K=2 completato il passo di integrazione 20 CONTINUE c verifica se si vuole stampare IF(ICOUNT/ISTEP*ISTEP.EQ.ICOUNT) THEN write(1,1000) x,(y(i),i=1,n) write(*,1000) x,(y(i),i=1,n) END IF C continua integrazione GOTO 8 CLOSE(1) 7 STOP 1000 FORMAT(1X,F8.2,1x,5(E10.4,1X)) END C....sssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss SUBROUTINE RUNGE(N,Y,F,X,H,M,K) C usa Runge-Kutta del 4.to ordine con i coefficienti'dí Gill C C integra il sistema di N equazioni ordinarie del prima ordine:



C dY(i)/dx=F(i) ; i=1, ... N ; con passo dx=H C M deve essere impostato alla prima chiamata pari a zero C K è un interruttore di switch,per K=2 l'integrazione è completa DIMENSION Y(15),F(15),Q(15) M= M+1 GOTO (1,4,5,3,7), M 1 DO 2 I= 1,N 2 Q(I)= 0. A= 0.5 GOTO 9 3 A=1.707107 4 X= X + 0.5*H 5 DO 6 I=1,N Y(I)=Y(I)+A*(F(I)*H-Q(I)) 6 Q(I)=2.*A*H *F(I)+(1.0-3.0*A)*Q(I) A= 0.2928932 GOTO 9 7 DO 8 I=1,N 8 Y(I)=Y(I)+H*F(I)/6.-Q(I)/3. M=0 K=2 GOTO 10 9 K=1 10 RETURN END C.....FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF SUBROUTINE FNC(N,Y,F,X) DIMENSION Y(15),F(15) c Definisce il sistema di 2 equazioni del secondo ordine c come un sistema di 4 equazioni del primo ordine c Y(1)=x c Y(2)=z dove x è l'ascissa, z è la profondità c Y(3)=u mentre le u,w sono rispettivamente le c Y(4)=w componenti della velocità lungo x e z g=9.81 Massa=800.0 Rho=1000.0 A=0.25 CL=0.0 CD=0.01 F(1)= Y(3) F(2)= Y(4) F(3)=(-0.5*Rho*A*CD*Y(3)*sqrt(Y(3)**2+Y(4)**2) & +0.5*Rho*A*CL*Y(4)*sqrt(Y(3)**2+Y(4)**2))/Massa F(4)=(-0.5*Rho*A*CD*Y(4)*sqrt(Y(3)**2+Y(4)**2) & +0.5*Rho*A*CL*Y(3)*sqrt(Y(3)**2+Y(4)**2))/Massa - g RETURN END



C=================================== Esercizio A. 4. 16 capire questo codice cosa sta integrando (back-engineering: a volta serve!) …………………………………………….. Esercizio A. 4. 17 vedere se il codice sopra riportato vi funziona Esercizio A. 4. 18 adattare il codice ad un altro problema vedi esercizi da A.4.12 a A.4.15

C.GOLIA: A

Aero_APP_

Considera

……………Esercizio ACapire il li…………… Esercizio AProva ad in……………

Aerodinamica

_a.docx

a il seguente

……………A. 4. 19 istato di cui ……………

A. 4. 20 ntegrare lo s……………

1

e file in MA

………………

sopra che p………………

stesso probl

………………

APLE:

….

problema ris……………

lema dell’es………..

solve ……………

sercizio 4.16 con MAPPLE

Pagi

23

ina |APP. 40

3/01/2009

0

9



Appendice 5

Codice NACAPROF Il codice NACAPROF genera i dati delle superfici per i profili NACA a 4 e 5 cifre. Alla partenza il codice chiede l’immissione

• del numero di nodi: sul ventre “Nventre” e sul dorso “Ndorso”, ovviamente il numero di nodi totali sarà pari a Nventre+Ndorso

• del codice del profilo NACA (4 o 5 cifre) che viene letto come intero: nacan= abc d e • del percorso e del nome del file di uscita (ad es. n4412.dat)

Sfruttando il fatto che nacan è un numero intero (la divisione tra interi genera un intero, si scartano gli eventuali decimali), con divisioni successive si estraggono i parametri del profilo:

• Si tratta di un profilo della famiglia a 5 cifre se [nacan/1000] > 10 • Spessore: {nacan-1000*[nacan/1000]-100*[nacan/100-10* nacan/1000]}*0.1

Se della famiglia a 4 cifre:

• Posizione max camber: [nacan/100-10* nacan/1000]*0.1 • Camber: [nacan/1000]*0.01

Se della famiglia a 5 cifre: • Posizione max camber: 0.2025 • Camber: 2.6595*(0.2025)3

I punti sulla corda sono individuati usando la legge del coseno:

2cos1

cx θ+

= per π=θ 2,...0 con

l’anomalia θ equamente ripartita. Sicché la numerazione dei punti avviene in senso orario:

• inizia da 1 nel bordo di uscita,

1

Nventre

Nventre+Ndorso

x/c

y/c

zz

camber(zz)/c

thick(zz)/c

Inclinazione Linea media

β=dcamber/dx

Linea media



• prosegue (sul ventre) fino ad Nventre (nel bordo di uscita) , • successivamente prosegue sul dorso fino ad arrivare a Ntot=Nventre+Ndorso nuovamente sul

bordo di uscita. In tale modo si realizza un addensamento dei nodi sui bordi. Su video appaiono i seguenti dati:

x(i)/c z(i)/c csl(i)/c ……………....ecc………….... dove csl(i)/c è lo sviluppo del profilo a partire dal b.d.u adimensionalizzato rispetto alla corda. Nel file di uscita vengono stampati i seguenti dati: Ntot x(i)/c z(i)/c ……………....ecc………….... Il codice genera, nella presente versione, un profilo chiuso: i.e. z(1) = z(Ntot) = 0.0 Nota in genere per molti profili i punti del bordo di uscita sono distinti per il ventre e per il dorso, per tenere conto dello spessore del lamierino usato per costruire, in pratica, l’ala. I dati del profilo vengono generati costruendo prima la linea media, e riportando successivamente su di essa la distribuzione degli spessori inclinata dell’inclinazione locale della linea media. Per ogni ascissa zz, i dati di spessore, camber ed angolo sono generati dalla

subroutine NACA45 (zz,thick,camber,beta,nacan.tau,epsmax,ptmax) che accetta in ingresso: zz, nacan,tau,epsmax,ptmax e fornisce in uscita: thick, camber, beta Ovviamente il codice permette di trattare famiglie diverse da quelle NACA cambiando semplicemente questa subroutine.



Appendice 6

Codice Pannelli Sorgenti Nel paragrafo 3.12.3, avevamo sottolineato che una delle metodologie più semplici per determinare il campo di moto non portante attorno a corpi simmetrici (e posti ad angolo d’attacco nullo) è di ricorrere ad un metodo a pannelli lineari con distribuzioni di sorgenti costanti su ogni pannello. A tale scopo si suddivide il contorno del corpo in N pannelli, che di solito, per successive implementazioni portanti, si fanno iniziare dal bordo d’uscita e si numerano in senso orario. Le N incognite σi=dQi/ds (che sono di per sé, funzioni Laplaciane) vengono determinate imponendo che il corpo sia linea di corrente. In una formulazione alla Neumann, questo viene fatto imponendo l’annullamento delle velocità normali al corpo in N punti di controllo, assunti essere medi degli N pannelli.

Ponendo ij

n

1jjinn,n NcosUVVV ∑

=∞∞ σ+β=+=

ne segue un sistema di equazioni lineari del tipo: ijji bN =σ Dove:

( )⎪⎩

⎪⎨

⎧

≠π

=

= ∫ jisedsdn

Rlnd21

jise2/1

Nj

i

ji

jpannello

ji

ii cosUb β−= ∞

• Nij è detto coefficiente di influenza del pannello j sul pannello i e rappresenta la velocità indotta

dal pannello j nel punto di controllo del pannello i e normalmente al pannello stesso, • Rji è il raggio vettore dal generico punto del pannello j al punto di controllo del pannello i, • ni la normale al pannello i, • βi l’angolo orientato della normale al pannello i rispetto alla corrente asintotica (asse X). Una volta risolto il sistema di equazioni lineari per le σj , le velocità tangenziali negli N punti di controllo i sono calcolate tramite la:

ij

n

)ji(1j

jist,t TsinVVVV ∑≠

=∞∞ σ+β−=+=

dove Tij che rappresenta il coefficiente di influenza del pannello “j” sulla velocità tangenziale del punto di controllo del pannello “i” è dato da:

α o

o

o

o

o

o

o

o

o s

estremi pannello punti di controllo

P(x,y)

θ

ni

βi

pannello i-esimo

pannello j-esimo

o

o S j n j

s j

ds j

(X ,Y )

j+1(X ,Y )

j j

j+1

β φ jj

pannello j-esimo

pannello i-esimo

φ β

(X ,Y )

(X ,Y )R

x

y

ij

(x ,y )

i i

i i

i i

i+1 i+1

n i



( )∫ ∂

∂

π=

jj

i

ijij ds

sRln

21T

Da notare che è scomparso il contributo alla velocità tangenziale di un pannello su se stesso. A parte la prova matematica (3.12.3) , si può facilmente intuire, da un punto di vista fisico, che un pannello sorgente emette flusso di volume dalla sua superficie verso l'esterno e verso l'interno del corpo, in direzione perpendicolare alla tangente al pannello stesso e che quindi il contributo nella direzione tangenziale è nullo. Una volta calcolate le velocità superficiali, il coefficiente di pressione superficiale, per ogni

pannello, si determina facilmente come: 2

ii,p V

V1c ⎟⎟

⎠

⎞⎜⎜⎝

⎛−=

∞

Un codice di calcolo, rudimentalmente semplice, capace di realizzare tale formulazione è il PANSOR Il codice è composto da un main e da 5 subroutines chiamate di seguito dal main. La subroutine INDATA(N) provvede a leggere o generare i nodi estremi dei pannelli del corpo che sono conservati nei vettori XP(i) e YP(i), posti nell’area COMMON. A tal fine chiede un input per tre modalità:

1. lettura da file esterno in questo caso chiede di immettere il percorso ed il nome del file di input che deve essere strutturato in modo tale da contenere in prima linea il numero N dei nodi totali (<400) e successivamente N linee che riportano la X(i) e la Z(i) del i-esimo nodo che rappresenta una delle estremità dei pannelli

2. generazione automatica di ellissi/cilindro in questo casi chiede di immettere il rapporto degli assi b/a dell’ellisse (1 per cilindro) e successivamente il numero totale dei nodi dei pannelli (<400) i cui nodi vengono generati con la legge del coseno con anomalie Δθ costanti

3. generazione automatica di profili NACA simmetrici in questo caso chiede di immettere il valore dello spessore percentuale (i.e. 12 per NACA 0012) e quindi in ) e successivamente il numero totale dei nodi dei pannelli (<400) i cui nodi vengono generati con la legge del coseno con anomalie Δθ costanti

La subroutine COEFF rappresenta il cuore del programma in quanto provvede a generare i coefficienti di influenza Nij e Tij che per ragioni di linguaggio sono contenuti nelle matrici AN(i,j) e AT(i,j). All’uopo:

• genera i punti di controllo come mezzerie dei pannelli e li conserva nei vettori XC(i), YC(i) • calcola la lunghezza di ogni pannello che viene conservata nel vettore S(i) ed il perimetro

del profilo SUP • calcola i coseni direttori dei versori normali [ DTXY(i,1-2) ] e tangenziali [ DNxy(i,1-2) ] ai

pannelli (1 componente asse x, 2 componente asse y)

Poiché questo codice risolve l’equazione di campo moltiplicando entrambe i termini per 2π:

ijji b2N2 π=σπ ne risulta, per i termini diagonali delle matrici di influenza: A(i,i)= ½ * 2 π = π , T(i,i)=0



Le altre componenti delle matrici di influenza sono calcolate, facendo riferimento alla figura, come segue: Si vogliono usare le formule di integrazione delle distribuzioni di sorgenti (2.2.6.1) per le componenti della velocità indotta in un punto,

( )( ) ⎥

⎥⎦

⎤

⎢⎢⎣

⎡

+−++

πσ

= 22

22

sorg z2/xz2/xln

21

2)z,x(u ⎥

⎦

⎤⎢⎣

⎡

−+πσ

= −222

1sorg )2/(zx

ztan2

)z,x(w (2.1)

poiché tali componenti sono relative ad un sistema di riferimento solidale al pannello inducente, le x,z delle formule devono essere le xq e yq della figura e le u e w delle formule diventano la Vx e Vy della figura. Le DX e DY, calcolabili nel sistema assoluto sono le componenti del vettore Rji che può essere proiettato sugli assi locali del pannello “j” per ottenere:

xq = R ji • t(j) yq = R ji • n(j)

e quindi ( σ(j) intensità unitaria)

( )( ) ⎥

⎥⎦

⎤

⎢⎢⎣

⎡

+−

++= 2

q2

q

2q

2q

y2/)j(sx

y2/)j(sxln

21Vx

⎥⎥⎦

⎤

⎢⎢⎣

⎡

−+= −

22q

2q

q1

)2/)j(s(zx

)j(sytanVy

che sono le velocità (unitarie) indotte dal pannello “j” nel punto di controllo del pannello “i” ma orientate secondo gli assi del pannello “j”. Per ottenere le componenti di velocità nelle direzioni tangenziale e normali al pannello “i”, occorre ripassare per il sistema assoluto: prima calcolare la velocità nel sistema assoluto: V ji,ass=Vx t(j)+Vy n(j) e proiettare questa sul pannello “i” : V ji,tangente al pannelli i = V ji,ass• t(i) V ji,normale al pannelli i

= V ji,ass• n(i) ovvero:

[ ] ⎥⎦

⎤⎢⎣

⎡ •+•==TNIJTTIJ

i,t )i(t)j(nVY)i(t)j(tVxV)j,i(AT [ ] ⎥⎦

⎤⎢⎣

⎡ •+•=≡DNIJDTIJ

i,n )i(n)j(nVY)i(n)j(tVxV)j,i(AN

Poiché si assume sempre velocità asintotica unitaria il termine noto è semplicemente b(i)= - n(i) • i ≡ A( I , N+1) La subroutine assembla il sistema nella matrice “aumentata” a destra [ ])1N,I(A|)J,I(ANA += Subroutine GAUSS(1) risolve il sistema di equazioni algebriche e pone la soluzione il A( I , N+1)

X

Y

x

yq

Pannello “j”

Pannello “i”

xq

t(j)

t(i)

n(j)

n(j)

y

DX

DY

Vx

Vy

R j i



Subroutine AEROCOF(N) calcola le velocità tangenziale nei punti di controllo secondo la

iij

n

)ji(1j

jt sinVTV β−σ= ∞

≠=∑

Do i=1,N VT(i)=0. Do j=1,N VT(i)=VT(i)+A(j,Nplus)*AT(I,J) Enddo VT(i)=VT(i)+DTXY(I,1) Enddo

La Subroutine OUTDATA(N) provvede all’output. Chiede se si vuole un output su file: se si (input1) chiede il percorso ed il nome del file di output in cui scrive (su file e su video) in prima linea il perimetro del profilo e quindi di seguito I,XC(i),YC(i),VT(i),CP(i) se NO (input 0) scrive solo su video



Appendice 7

Metodo dei Vortici Concentrati La risoluzione numerica di linee medie con una metodologia vorticosa porterebbe naturalmente a discretizzare il corpo con N pannelli lineari su cui disporre una distribuzione di vorticità γ(i). Tale metodologia, possibile, è però ostacolata da due problematiche: l’imposizione della condizione di Kutta al bordo di uscita e il non buon condizionamento nel sistema algebrico che ne deriva per la determinazione delle incognite γ(i), derivante dal fatto che la velocità normale indotta da un pannello vorticoso su se stesso è nulla (2.x.x.), sicché la diagonale principale della matrice è identicamente nulla. Ovviamente l’uso di vortici concentrati risolve entrambi i problemi: pannelli lineari sono da considerarsi come tante lastre piane, per le quali il loro valore della vorticità Γ , se concentrata nel punto neutro anteriore (1/4 della corda del pannello), è stimabile esattamente imponendo l’annullarsi della velocità totale (asintotica + indotta) normalmente al pannello nel punto neutro posteriore (3/4 della corda del pannello). E’ dimostrato [vedi paragrafo 4.8] che tale metodologia soddisfa la condizione di Kutta per il pannello. Inoltre è facilmente intuibile che l’influenza del pannello su se stesso è senz’altro maggiore delle influenze degli altri pannelli, sicché la matrice risultante è certamente diagonale dominante. L’implementazione di tale metodo è particolarmente semplice in quanto se il corpo è discretizzato in N pannelli, per ognuno di questi si conoscono: • le coordinate del fuoco XP(i),ZP(i) e quelle del punto di controllo XC(j),ZC(j), • i coseni direttori della normale al pannello nel punto di controllo DNX(j),DNZ(j) misurati

rispetto ad una corda di riferimento di lunghezza CRIF, • l’angolo di attacco “α” della velocità asintotica misurato rispetto alla stessa corda di

riferimento, • e quindi le componenti RXij ed RZij del vettore Rij che va dal fuoco del generico pannello i al

punto di controllo del generico pannello j, Le vorticità incognite concentrate nei punti focali GAMMA(j) sono soluzione del sistema:

ijij bA =Γ dove: • i coefficienti di influenza Aij (velocità indotta da un vortice unitario nel fuoco di j, nel punto di

controllo del pannello "i" nella direzione normalmente al pannello stesso) sono dati da

y

2ij

iijij

R

nR

21A ⎟

⎟

⎠

⎞

⎜⎜

⎝

⎛ ∧

π=

• il termine noto (componente della velocità asintotica nella direzione normale al pannello i) è dato da: ii nVb •−= ∞

ovvero più semplicemente: 2π 2

IJiijiijij R/)DNX*RzDNZ*RX(A −= ; -2π )DNZ*VZDNX*VX(*2b iii +π−= laddove:

Rxij = XP(j)-XC(i) Rzij = ZP(j)-ZC(i) Rij2 = RXij

2+RZij2



VX = cosTETA VZ = sinTETA

Una volta note le intensità dei vortici e la loro collocazione è immediato ricavare i coefficienti del sistema di forze:

Crif2C

i

N

1il

Γ

=∑

= 2

xp2C

Crif

ii

N

1i0M

Γ

=∑

= l

0Mcp

CC

Crifx

−= 25.0*CCC l0M)4/1(M +=

Tale metodologia è usata dal programma VORCON, che consiste da un MAIN e da una SUBROUTINE e da una FUNCTION. Il MAIN chiede in input: 1. leggere dati da un file esterno

in questo caso chiede di digitare il percorso ed il nome del file dati che deve contenere nella prima linea:

il numero di pannelli N ed il valore della corda di riferimento Crif e su ognuna delle N righe successive

XP(i), ZP(i), XC(i), ZC(i), DNX(i), DNZ(i) 2. generare una linea media NACA

in questo caso si assume corda unitaria (Crif=1) e si chiede di immettere nell’ordine: • numero dei pannelli (N) • camber massimo (RM) (ex. per NACA 2412 immettere 0.02) • posizione del camber massimo RP (ex. per NACA 2412 immettere 0.40)

il programma genererà automaticamente N pannelli di pari lunghezza di cui calcolerà i parametri necessari all’algoritmo richiamando la FUNCTION ZNACA(X,RM,RP). E’ ovvio che se si volessero esaminare linee medie di famiglie diverse da quelle della NACA basta sostituire questa function.

Il programma richiede quindi di immettere l’angolo di attacco α (in gradi) e ad abbundantiam fa una stampa su video per la verifica dei dati immessi, dopo di che chiede di immettere un numero qualsiasi se i dati vanno bene o di immettere ESC se si vuole uscire. Il programma forma il sistema di equazioni (matrice aumentata con il termine noto) che viene risolto con la solita subroutine GAUSS, con soluzione nella colonna N+1. La soluzione XP(i), ZP(i), G(i) viene stampata su video e su di un file di output (per eventuali plottaggi) di cui il programma chiede l’immissione del percorso e del nome. Il programma quindi calcola i parametri aerodinamici sopra descritti che vengono stampati solamente su video.



Appendice 8

Metodo Pannelli Portanti Il metodo più semplice per la risoluzione di problemi portanti attorno a profili arbitrari (senza la limitazione dei piccoli disturbi) è senza dubbio quello ideato in casa DOUGLAS e dovuto a Smith & Hess. Basato sulla idea di discretizzare il profilo (di qualsiasi forma, anche non simmetrici) con N pannelli, e volendo mantenere una metodologia alla Neumann, l’idea si basa sulla estensione dei metodi a pannelli sorgenti (di cui in App. 6) a problemi portanti mediante la introduzione di vorticità. Poiché il problema non portante è determinato dalle condizioni di tangenza nei punti di controllo degli N pannelli e quello portante è determinato dalla condizione di Kutta, vi sono a disposizione N+1 equazioni. La soluzione più semplice, per chiudere il problema, è di introdurre una sola incognita per la vorticità, che Smith&Hess pensarono di distribuire uniformemente su tutto il contorno. Quindi su ogni pannello del profilo sussiste sia una distribuzione di sorgenti, σi , che varia da pannello a pannello, sia una distribuzione di vorticità γ che è la stessa per tutti i pannelli. Ne deriva che le velocità, nei punti di controllo, dell’i-esimo pannello possono porsi come:

• componente normale: .vortj,i

N

1jj

.sorgj,i

N

1ji

.vorti,n

.sorgi,ni,ni,n ANANnVVVVV ∑∑

==∞

∞ γ+σ+•=++=

• componente tangenziale: .vortj,i

N

1jj

.sorgj,i

N

1ji

.vorti,t

.sorgi,ti,ti,t ATATtVVVVV ∑∑

==∞

∞ γ+σ+•=++=

dove: • ANsorg, ATsorg sono rispettivamente i coefficienti di influenza normale e tangenziale delle

distribuzioni di sorgenti • ANvort, ATvort sono rispettivamente i coefficienti di influenza normale e tangenziale delle

distribuzioni di vortici La condizione di parete impermeabile porta ad N equazioni: Vn,i=0 , ∀i=1, .. ,N ovvero:

N,..,1inVANAN i.vort

j,i

N

1jj

.sorgj,i

N

1j

=∀•−=γ+σ ∞==

∑∑ (A8.1)

La condizione di Kutta, che teoricamente dovrebbe essere imposta nel bordo di uscita, non può essere implementata in modo esatto, e viene approssimata imponendo eguale pressione nella parte del ventre e del dorso immediatamente adiacente al bordo di uscita, cioè imponendo [teorema di Bernoulli] che le velocità tangenziali nei punti di controllo dei due pannelli di coda (nella nostra notazione il primo e l’N-esimo) siano uguali. Questa approssimazione è tanto più valida quanto più prossimi sono i due punti di controllo al bordo di uscita, e soprattutto quanto più simile sia la loro distanza dal b.d.u (questo richiede che il pannello 1 ed il pannello N abbiano la stessa lunghezza e che questa sia quanto più piccola possibile).



bordo d’uscita

1

2

NN-1

Vt1>0

VtN>0

Ovviamente a causa della numerazione (oraria a partire dal b.d.u.) nel primo pannello la tangente sarà diretta verso il bordo di attacco, cioè avrà direzione contraria alla corrente asintotica, nel mentre quella dell’ennesimo pannello sarà concorde alla corrente asintotica. Ne segue che l’imposizione di iso-baricità nei pressi del b.d.u risulta approssimata dall’imposizione di

N,t1,t VV =− ovvero

0VV N,t1,t =+ ovvero:

0ATATtVATATtV .vortj,N

N

1jj

.sorgj,N

N

1jN

.vortj,1

N

1jj

.sorgj,1

N

1j1 =γ+σ+•+γ+σ+• ∑∑∑∑

==∞

==∞

ovvero:

[ ]N1.vort

j,N

N

1j

.vortj,1

N

1jj

.sorgj,N

N

1j

.sorgj,1

N

1j

ttVATATATAT +•−=⎥⎥⎦

⎤

⎢⎢⎣

⎡+γ+σ

⎥⎥⎦

⎤

⎢⎢⎣

⎡+ ∞

====∑∑∑∑ (A8.2 )

Le (A8.1 e A8.2) rappresentano un sistema di N+1 equazioni delle N+1 incognite σj + γ, che può rappresentarsi in forma matriciale-partizionata come segue:

⎥⎦

⎤⎢⎣

⎡=⎥

⎦

⎤⎢⎣

⎡γσ

⎥⎥⎦

⎤

⎢⎢⎣

⎡

2

1

2221

1211bb

AAAA

dove: A11 è una matrice quadrata di dimensioni N x N simile a quella del programma PANSOR

⎥⎥⎥⎥⎥

⎦

⎤

⎢⎢⎢⎢⎢

⎣

⎡

=

sorgN,N

sorgN,1

sorgN,1

sorg1,1

11

ANAN

ANAN

A

A12 è una matrice-vettore colonna di dimensioni N x 1

⎥⎥⎥⎥⎥⎥⎥⎥

⎦

⎤

⎢⎢⎢⎢⎢⎢⎢⎢

⎣

⎡

=

∑

∑

=

=

vortj,N

N

1j

vortj,1

N

1j

12

AN

AN

A



A21 è una matrice-vettore riga di dimensioni 1 x N avente componenti non nulle solo nelle posizioni 1 ed N

[ ] [ ]⎥⎥⎦

⎤

⎢⎢⎣

⎡++= ∑∑

==

.sorgj,N

sorgj,1

N

1j

.sorgj,N

sorgj,1

N

1j21 ATAT00ATATA

A22 è una matrice scalare di dimensioni 1x1 (un sol componente)

[ ]⎥⎥⎦

⎤

⎢⎢⎣

⎡+= ∑

=

vortj,N

vortj,1

N

1j22 ANANA

b1 un vettore colonna di dimensioni N x 1

⎥⎥⎥⎥⎥

⎦

⎤

⎢⎢⎢⎢⎢

⎣

⎡

•−

•−

=

∞

∞

N

1

1

nV

nV

b

b2 un vettore scalare di dimensioni 1x1 (un solo componente)

[ ][ ]111 ttVb +•−= ∞ Nota che a causa della simmetria tra sorgenti (campi radiali) e vortici (campi circonferenziali)

ANvort =- ATsorg e ATvort=ANsorg e che quindi non occorre determinare niente di nuovo rispetto al caso del PANSOR. In pratica, quindi, per implementare un codice portante occorre soltanto portare alcune modificazioni al codice a pannelli sorgente di cui all’Appendice 6. Il codice SMITH_HESS è composta da un main e da 8 subroutines. Il Main provvede a chiamare nell’ordine le varie subroutines e a richiedere al momento giusto l’immissione dell’angolo d’attacco ALFA (in gradi); Nota che se si immette un valore maggiore di 90 il programma si chiude. La Subroutine INDATA chiede di inserire separatamente il numero di punti sul ventre (Nventre) e poi quello sul dorso (Ndorso) con l’avvertenza che la loro somma deve essere minore di 400. Chiede successivamente di immettere la sigla del profilo naca, da cui deriva i dati del profilo secondo la classificazione della Naca (idem Appendice 5) La subroutine SETUP fissa le coordinate dei punti sulla corda (unitaria) con la legge del coseno e i punti del dorso/ventre che vengono generati chiamando la subroutine BODY(Z,SIGN,X(I),Y(I)). Calcola le pendenze dei pannelli e conserva i seni [SINTHE(I)]ed i coseni [COSTHE(I)] degli angoli orientati dei pannelli rispetto alla corda, stampa su video il numero del pannello e le sue coordinate [I,X(I),Y(I)] La subroutine BODY(Z,SIGN,X(I),Y(I)) genera automaticamente le coordinate dei punti dei profili NACA chiamando la subroutine NACA45(Z,THICK,CAMBER,BETA)



La Subroutine COFISH(SINALFA,COSALFA) definisce le matrici del metodo e forma la matrice di risoluzione del problema A(N+1,N+2) aumentata del termine noto. La Subroutine GAUSS(1) risolve il sistema e pone la soluzione nella N+2-esima colonna di A La Subroutine VELDIS(SINALFA,COSALFA) rileva la soluzione dalla N+2-esima colonna di A e la stampa su video [I,Q(I)], calcola, nel punto mediano di ogni pannello, la velocità tangenziale (che coincide con la velocità totale essendo sul corpo la componente normale nulla) e il coefficiente di pressione e li stampa su video [XMID(I),VEL(I),CP(I)] La Subroutine FANDM(SINALFA,COSALFA,alfa) calcola e stampa su video i dati aerodinamici: i coefficienti di portanza, resistenza indotta e momento (rispetto al bordo d’attacco): [CL,CD.CM] e poi, se CL≠0, stampa su video la posizione del centro di pressione [Xcp/c], e poi il coefficiente di momento a c/4 (molto prossimo al coefficiente di momento focale) [CM(c/4)] e chiede se si vuole una stampa su file: • in caso positivo chiede percorso e nome del file, in cui scrive [NODTOT,ALFA], e poi

[XMID(I),YMID(I),CP(I),VEL(I)] • chiede poi se si vuole un output per file per analisi di Strato Limite

• in caso positivo richiede percorso e nome del file, cerca il punto di ristagno, separa i dati del dorso aerodinamico e quelli del ventre aerodinamico e stampa nel file indicato prima [Ndorso, Nventre,ALFA] e poi successivamente Ndorso righe che, iniziando dal punto di ristagno, riportano i dati del dorso aerodinamico, [X(i),Y(i),VEL(i)] seguiti da Nventre righe che iniziando dal punto di ristagno, riportano i dati del ventre aerodinamico [X(i),Y(i),VEL(i)].

• infine chiede se si vuole anche un file di output per la geometria del profilo • in caso positivo chiede percorso e nome del file, in cui stampa [NODTOT] e poi [X(i),Y(i)]

dopo di che si torna al main che chiede, per un nuovo caso, di immettere il valore dell’angolo d’attacco ALFA (in gradi): se ALFA > 90 il programma viene terminato. Nota il programma, anche se non lo sembra, è molto semplice ed orientato alla determinazione della polare di profili NACA. Nel caso si debbano usare altri profili occorre variare la subroutine NACA45, ovvero implementare il programma con una opzione per permettere di leggere i dati del profilo da un file esterno

Appendice 9

Ala Finita Lineare Dritta Le ali lineari a grande allungamento e senza freccia in pianta (rispetto alla linea dei fuochi) sono modellate da teorie vorticose che usano un segmento vorticoso aderente alla linea retta dei fuochi dell’ala e da una superficie vorticosa di vortici liberi che si estende nella direzione dell’asse delle x (ali poco caricate!). Ne deriva che la distribuzione di vorticità aderente e quindi libera, deve soddisfare l’equazione integro-differenziale di Lanchester-Prandtl. La risoluzione di tale equazione viene fatta seguendo l’analisi di Glauert :



a) scambiando l’ascissa di apertura “y” con una anomalia θ, [ ] θ−=π∈θ∀ cos2by:,0

b) rappresentando la Vorticità aderente con uno sviluppo in serie di Fourier di solo seni dell’anomalia,

( ) ( )∑∞

=∞ θ=θΓ

1nn nsinAbV2

c) usando le formule di Glauert per eliminare gli integrali impropri presenti, d) per arrivare ad una equazione per gli infiniti coefficienti An della serie che deve valere per

tutti i punti dell’apertura [-b/2≤y≤b/2, i.e. per tutti i valori dell’anomalia θ, [ ]π∈θ∀ ,0 ]

[ ] { }[ ]∑∞

=

θ+μθ=θα−αμ1n

n0,L sinn)n(sinAsin (A9.1)

dove:

• si è posto: )(4

θμμ ==b

ac ass

• aass è il coefficiente della retta di portanza del tronco di ala in esame, considerato operare nel tratto lineare:

( ) ( ) ( )iLass

Lieff

ass

LieffassL d

dCddCaC ααα

ααα

ααα −−⎟

⎠⎞

⎜⎝⎛=−⎟

⎠⎞

⎜⎝⎛=−= 0

La risoluzione dell’equazione (A9,1) viene fatta con la tecnica della collocazione: assumendo un numero finito N* di armoniche ovvero di An , [i.e. ∀ n=1,..,N* ] e imponendo che l’equazione (approssimata) sia verificata in N* punti dell’apertura dell’ala.

In parole povere, dovreste aver capito che lo scopo di tutta l’ingegneria è di

riuscire ad avere N* equazioni per le N* incognite con una matrice di sistema ben condizionata.

Questo processo di troncamento/collocamento è chiamato nella letteratura inglese, risoluzione dell’ equazione del monoplano.

Ovviamente se si considera un’ala simmetrica (in forma e carico, i.e. senza imbardata e rollio) basterà considerare soltanto i termini dispari [i.e. n=1,3,5,7,....] e una sola semiala [ ]π∈θ∀ ,0 e raddoppiare i valori dei coefficienti trovati:

( ) ( ) ( ) ( ) ( )[ ]Γ θ θ θ θ θ= + + + +∞2 3 5 71 3 5 7bU A sin A sin A sin A sin ........ Una volta noti i valori dei coefficienti della serie ( An ) si possono calcolare tutti i parametri di interesse:

( ) ( )L b U A C U SL= ≡∞ ∞ ∞ ∞2 1

2 112ρ π ρ

C AbS

AL = = ℜ1

2

1π π

( )α θθθ

θθ

θθi A A

sinsin

Asin

sinA

sinsin

= + + + +1 3 5 733

55

77( ) ( ) ( )

.....



C A nAA

C AA

AA

AAD i

n

n

NL

, .......= ℜ +⎛⎝⎜

⎞⎠⎟

⎡

⎣⎢⎢

⎤

⎦⎥⎥

⎡

⎣⎢⎢

⎤

⎦⎥⎥

=ℜ

+⎛

⎝⎜

⎞

⎠⎟ +

⎛

⎝⎜

⎞

⎠⎟ +

⎛

⎝⎜

⎞

⎠⎟ +

⎡

⎣⎢⎢

⎤

⎦⎥⎥=

∑ππ1

2

1

2

2

232

12

52

12

72

121 1 3 5 7

Il codice MONOPLANcf risolve ali lineari, a grande allungamento, senza freccia in pianta, con rastremazione, svergolamento (geometrico ed aerodinamico) e con flaps. Consiste da un main che chiama, di volta in volta, quattro subroutines. La subroutine GEOMONOP legge o genera la geometria e l’angolo d’attacco, all’uopo propone l’opzione: 1. di immettere da tastiera i dati geometrici. In questo caso chiede di immettere da tastiera:

• l’apertura alare B (in metri), • la corda alla radice CR (in metri), • la corda all’estremità CE (in metri), • l’angolo di svergolamento alla radice ASVR (in gradi), • l’angolo di svergolamento all’estremità ASVE (in gradi), • per la sezione di radice: il coefficiente delle retta di portanza (1/rad) e l’angolo di

portanza nulla (in gradi positivo verso il basso) CLAR,ALOGR • per la sezione di estremità: il coefficiente delle retta di portanza (1/rad) e l’angolo di

portanza nulla (in gradi positivo verso il basso) CLAE,ALOGE 2. di leggerli da un file. In questo caso chiede di immettere da tastiera il percorso ed il nome di un

file che deve contenere 9 righe: 1. l’apertura alare B (in metri), 2. la corda alla radice CR (in metri), 3. la corda all’estremità CE (in metri), 4. l’angolo di svergolamento alla radice ASVR (in gradi positivo verso il basso), 5. l’angolo di svergolamento all’estremità ASVE (in gradi positivo verso il basso), 6. l’angolo di portanza nulla (in gradi positivo verso il basso) per la sezione di radice ALOGR, 7. l’angolo di portanza nulla (in gradi positivo verso il basso) per la sezione di radice ALOGE, 8. il coefficiente delle retta di portanza (1/rad) per la sezione di radice CLAR, 9. il coefficiente delle retta di portanza (1/rad) per la sezione di estremità CLAE,

Successivamente: • chiede di immettere l’estensione del flap in % della semiapertura alare PFALP (se 0 no flap)

• se PFLAP ≠0 chiede di immettere l’incremento dell’angolo di portanza nulla causato dall’abbassamento del flap nelle sezioni (in gradi, negativo se flap down) DFALO

• chiede di immettere: • il numero di punti lungo la semi-apertura N, che posizione con la regola del coseno, • l’angolo di attacco ALFAG (in gradi misurata rispetto alla corda di riferimento (di solito la

corda media dell’ala). Essendo un’ala lineare, i dati delle N sezioni alle varie stazioni intermedie vengono interpolati linearmente tra quelli di radice ed estremità e stampati su video: I, TETA(I),Y(I),C(I),SVERG(gradi), Clalfa,Alfa0(gradi)



La subroutine COFINMONOP genera i coefficienti del sistema di equazioni del monoplan, ne deriva una matrice:

A(I,J)= SIN(Jdispari*TETA(I))+(RMU*Jdispari+SIN(TET(I)) ed un termine noto che viene collocato nella matrice A aumentata a destra:

A(I,N+1)=RMU*(ALFAR-ASV(i)-ALO(I)*SIN(TET(I)) La subroutine GAUSS(NRHS) risolve il sistema e restituisce la soluzione nella (n+1)-esima colonna di A:

An(i)=A(I,N+1) La subroutine COEFFMONOP calcola i coefficienti aerodinamici dell’ala (raddoppia quelli calcolati per la semiala) e le distribuzioni lungo l’apertura e provvede alle stampe su video e su file. All’uopo:

• recupera la soluzione dall’output di GAUSS per ogni stazione Jdispari, calcola: • il coefficiente di portanza locale CL(I) come somma rispetto a Jdispari di

AN(j)*SIN(J*TET(I)*2*4*BH/C(I), • l’angolo di incidenza indotta locale AI(I) come somma rispetto a Jdispari di

J*AN(j)*SIN(J*TET(I) /SIN(TET(I)), • il coefficiente di resistenza indotta locale CD(I) come prodotto di AI(I)*CL(I), • il parametro di resistenza indotta locale SCDT come somma rispetto a Jdispari di

J*(AN(j)/AN(1))**2 e ne stampa su video la distribuzione lungo l’apertura: Y(I), CL(I), Y(I)/BH, CL(I)/CLradice, CD(i) calcola i coefficienti totali:

• di portanza CLT= PIG*AR*AN(1) • di resistenza indotta CDT=PIG*AR*AN(1)**2*SCDT

quindi stampa su video • Apertura alare • Angolo d’attacco, alfa(gradi) • Estensione del flap, Delta(alfal0) derivante dall’abbassamento del flap • rastremazione di corde • angoli di svergolamenti (radice et estremità) • Angoli di portanza nulla (radice et estremità) • coefficienti totali di Portanza e di resistenza indotta chiede se si vuole output su file? • se si chiede di immettere percorso e nome del file su cui stampa:

• eventualmente il nome del file di input usato per l’ala • angolo d’attacco (in grado) e coefficienti totali di Portanza e di resistenza indotta • le distribuzioni lungo l’apertura Y(I), CL(I), Y(I)/BH, CL(I)/CLradice,

CD(i)



USO DEL PROGRAMMA MONOPLAN PER L’ANALISI (di primissimo livello) DELL’ALA.

Tutti questi casi prevedono l’utilizzo di molti run del programma, per cui è preferibile usare input mediante un file esterno (dal nome preferibilmente corto e posto nella stessa directory dell’eseguibile). Sia ala0.txt il file che contiene i dati dell’ala, proveniente dagli aerodinamici. 1. CALATTAMENTO DELL’ALA L’ufficio di progetto deve calettare (posizionare con un dato angolo) l’ala su di un velivolo e vuole farlo con qualche criterio, ad esempio; se si tratta di un velivolo civile, è pensabile che per la stragrande maggioranza del tempo il velivolo (di massa nota) volerà ad una certa quota, con una certa velocità. Da questi dati è facile ricavare il CLalare di regime. Infatti dall’eguaglianza PESO=PORTANZA:

( ) ala.rif2

21

alare SVquota alla*CLPortanzag*carico)fluidi e carburante vuotoa massa(PESO ρ==++= si deriva:

( ) ala.rif2

21alare SVquota alla

g*carico)fluidi e carburante vuotoa massa(CLρ

++=

Ne esce un valore che generalmente è compreso tra 0.4 e 0.5, assumiamo ad esempio CLalare=0.45.

Ora è abbastanza ovvio che in questo caso si vuole che l’ala sia posizionata rispetto al velivolo in modo tale che quando l’ala genera CLalare=0.45, la fusoliera sia posizionata rispetto alla direzione della corrente ad angolo d’attacco nullo (per fare in modo che la fusoliera generi una resistenza minima!).

Il problema è quindi di calettare opportunamente la corda di riferimento dell’ala rispetto all’asse di fusoliera.

Molto semplice:

• Fare vari runs del codice MONOPLAN usando il file ala0.txt a vari angoli d’attacco (che saranno rispetto alla corda di riferimento dell’ala) e trovare quello per cui si raggiunge CLalare=0.45; sia ad esempio:

αcorda rif ala=+1.5(°) che rappresenta l’angolo di calettamento

• Andare nel file alfa0.txt e sottrarre ad entrambi gli angoli di svergolamento , di radice

ASVR e di estremità ASVE, il valore trovato di αcorda rif ala=+1.5(°). • Salvare il file con un altro nome, ad es. ala0_calettato_1.5.txt • Controllare con un nuovo run che usando questo file e ponendo α=0 risulti CLalare=0.45.

Se ciò avviene abbiamo trovato l’angolo di calettamento dell’ala sul nostro velivolo. Evidentemente questo è solo la stima di primo tentativo, dovrete in seguito considerare altri fattori quali: le interferenze, i down-wash … le condizioni di off-design….eccetera….



2. POLARE DELL’ALA Gli aerodinamici avranno certamente generato, si suppone seriamente, le curve di polare dell’ala. Da buon seguaci di S.Tommaso, potrebbe venirvi voglia di ri-calcolarle, diciamo per sfizio o per una verifica semplicistica, a volo d’uccello, sfruttando il codice MONOPLAN. Usate il file dell’ala calatteta, ed assumete che l’ala continui ad avere comportamento scalabile linearmente, sezione per sezione, ad ogni angolo d’attacco, fino allo stallo. Dovete avere a disposizione:

• le polari dei profili di radice ( R ) e di estremità ( E ), i.e. dovrete, per ogni valore dei rispettivi coefficienti di portanza CLR e CLE, conoscere i valori dei coefficienti di resistenza parassita CDpR, CDpE ed i valori dei coefficienti di momento focale CMFR,CMFE,

• le corde dei profili di radice CR e di estremità CE e la corda media Cm,

Fate diversi run con il codice MONOPLAN, facendo incrementare l’angolo d’attacco e per ogni run abbiate cura di creare una tabella come quella che segue: alfa 0 2 4 6 8 9 10 CL_ala CDind_ala CLR CLE CDpR = f(CLR) CDpR = f(CLR) CMFR = f(CLR) CMFE = f(CLR) CDp_ala=(CDpR*CR + CDpE*CE)/CM

CD_ala = CDp_ala + CDind_ala CMF_ala = (CMFR*CR + CMFE*CE)/Cm

NOTA :

• le caselle bianche sono riempite direttamente dagli outputs del codice MONOPLAN • le caselle grigio chiare devono essere calcolate, sfruttando le polari dei profili di radice e di

estremità per determinare i coefficienti di resistenza di profilo e di momento focale in funzione dei coefficienti di portanza locali,

• le caselle grigio più scuro forniscono la prima grossolana valutazione del coefficiente di resistenza parassita CDp_ala e quindi del coefficiente di resistenza totale CD_ala e del Momento focale CMF_ala (MEMO per l’ala diritta tutte le sezioni hanno lo stesso fuoco!) in funzione di CL_ala e di alfa

Ergo: da questa tabella si può tracciare una prima grossolana stima delle curve polari

dell’ala!



y = 0 y = b/2

CLmax,radice

CLmax,estremità

Curva CLmaxlungo l’apertura

3. INIZIO E CAMMINO DI STALLO Dovrete plottare l’andamento del CLmax dei profili significativi della vostra ala, lungo l’ apertura (y). Questo è semplice per profili non flappati. Infatti nello spirito della linearità basta individuare dai data-base aerodinamici disponibili (in funzione dei numeri di Reynolds calcolati in base alle corde locali) il CLmax del profilo di radice, il CLmax del profilo di estremità e unirli con una linea retta.



y = 0 y = b/2

CLmax,radiceCLmax,estremità

CLmax,radice

con flap CLmax,estremità

con flap

zona di estensionedei flap

Curva CLmaxlungo l’aperturacon flap azionati

y = 0 y = b/2

CLmax,radice

CLmax,estremità

stazione diinizio stallo

α2

α1

α3

α4= angoloinizio stallo

Curva CLMAXlungo l’apertura

per ogni α:curve CL(y)e un valore

di CLala

α5 >α4

Più complicato il caso di un’ala flappata, in questo caso dovrete individuare anche il CLmax del profilo di radice flappato ed il CLmax del profilo di estremità flappato, unire i dati dei CL con flap e dei CL senza flap, ed, in funzione della reale estensione del flap, assumere l’una o l’altra linea. Spiegazione troppo complicata,

meglio vedere l’esempio in figura.

Aumentando ad ogni run l’angolo

d’attacco dell’ala, avrete distribuzioni di carico CL(y) lungo l’apertura che variano con α, che plotterete idealmente su questi diagrammi, per verificare se la curva disegnabile dai dati di output del codice CL(y), tocca la curva CLmax che avete ricavato dai dati aerodinamici dei profili di radice e di estremità. Quando ciò accade avrete fatto una prima stima:

• dell’angolo d’attacco di inizio stallo, ed il corrispondente valore di CLstallo dell’ala • la stazione di inizio stallo.

Bada Bene: in queste condizioni il comportamento aerodinamico dell’ala è fortemente non-lineare, quindi quella fatta è solo una stima preliminare. Ad assetti ancora maggiori lo stallo interesserà altre sezioni, si potrà così stimare la forma del cammino dello stallo. Ma non illudetevi, sarà solo una stima qualitativa, non sostenuta da alcuna razionalità (dopo l’inizio stallo dovreste escludere le zone stallate altrimenti potrebbe accadere che il CLala continui ad aumentare!).

Aero_2009_APPENDICI.pdf

Documents

Transcript of Aero_2009_APPENDICI.pdf