Modelarea matematică si numerică a curgerilor [617044]
UNIVERSITATEA POLITEHNICA BUCURE STI
FACULTATEA DE INGINE RIE AEROSPATIALA
1
Modelarea matematică si numerică a curgerilor
compresibile in regim laminar si turbulent
Absolvent: [anonimizat] : Prof. d r. ing. Dănăilă Sterian
BUCUREȘTI
2017
UNIVERSITATEA POLITEHNICA BUCURE STI
FACULTATEA DE INGINE RIE AEROSPATIALA
2
CUPRINS
1 Introducere ………………………….. ………………………….. ……………….. 4
2 Aspecte fizice ale mișcărilor compresibile ……………………… 5
3 Modele matematice ………………………….. ………………………….. ….. 8
3.1 Ecuațiile Navier -Stokes ………………………….. ………………………….. …….. 8
3.1.1 Legi de conservare ………………………….. ………………………….. ………………………….. …………………….. 8
3.1.2 Teorema de transfe r ………………………….. ………………………….. ………………………….. …………………… 9
3.1.3 Conservarea masei ………………………….. ………………………….. ………………………….. ………………………. 10
3.1.4 Conservarea impulsului ………………………….. ………………………….. ………………………….. ……………. 11
3.1.5 Conservarea energiei ………………………….. ………………………….. ………………………….. ……………….. 12
3.1.6 Formularea conservativă locală ………………………….. ………………………….. ………………………….. .. 14
3.1.7 Formularea integrală conservativă ………………………….. ………………………….. ………………………… 15
3.1.8 Formularea neconservativă locală ………………………….. ………………………….. ………………………… 15
3.1.9 Formularea tensorială ………………………….. ………………………….. ………………………….. ………………. 16
3.1.10 Formularea adimensională ………………………….. ………………………….. ………………………….. ……….. 16
3.1.11 Comentarii privind ecuațiile Navier -Stokes ………………………….. ………………………….. ……………. 18
3.2 Ecuațiile mediate Reynolds ………………………….. …………………………. 18
3.2.1 Medierea ecuațiilor Navier -Stokes ………………………….. ………………………….. ………………………… 18
3.2.2 Comentarii privind ecuațiile mediate Reynolds ………………………….. ………………………….. ……… 20
3.3 Ecuațiile Euler ………………………….. ………………………….. …………………. 21
3.3.1 Proprietăți matematice ale ecuațiilor Euler ………………………….. ………………………….. ……………. 21
3.3.1.1 Formularea integrală ………………………….. ………………………….. ………………………….. 22
3.3.1.2 Formularea diferențială co nservativă ………………………….. ………………………….. …. 22
3.3.1.3 Matricele iacobiene în formularea conservativă ………………………….. ……………… 23
3.3.1.4 Formularea în varibile primitive ………………………….. ………………………….. ………….. 24
3.3.1.5 Legătura dintre formularea conservativă și cea în variabile primitive …………… 24
3.3.1.6 Variabilele caracteristice ………………………….. ………………………….. …………………….. 25
3.3.1.7 Condiții la limită ………………………….. ………………………….. ………………………….. ……… 28
3.3.1.8 Precondiționarea sistemului Euler ………………………….. ………………………….. ……… 28
3.3.1.9 Problema Riemann ………………………….. ………………………….. ………………………….. … 32
3.3.1.10 Matricele iacobiene în variabile primitive ………………………….. ………………………… 34
3.3.1.11 Valorile proprii ale sistemului Euler ………………………….. ………………………….. …….. 34
4 Metode numerice pentru rezolvarea ecuațiilor
Navier -Stokes ………………………….. ………………………….. ……………… 35
4.1 Metode cu volume finite ………………………….. ………………………….. ….. 36
4.1.1 Discretizări conservative ………………………….. ………………………….. ………………………….. …………… 36
4.1.2 Scheme cu volume finite ………………………….. ………………………….. ………………………….. ………….. 38
4.2 Discretizarea în regim compresibil ………………………….. …………………… 40
4.2.1 Discretizarea fluxurilor convective ………………………….. ………………………….. ………………………… 40
4.2.1.1 Scheme upwind de tip Godunov ………………………….. ………………………….. ………… 40
4.2.2 Discretizarea fluxurilor d ifuzive ………………………….. ………………………….. ………………………….. …. 42
UNIVERSITATEA POLITEHNICA BUCURE STI
FACULTATEA DE INGINE RIE AEROSPATIALA
3 5 Modele de turbulență ………………………….. ………………………….. … 44
5.1 Descrierea statistică a turbulenței ………………………….. ……………….. 44
5.2 Modelul
k standard ………………………….. ………………………….. ……… 46
6 Programul NSC2KE ………………………….. ………………………….. … 49
6.1 Ecuațiile guvernatoare ………………………….. ………………………….. …….. 50
6.2 Tehnici numerice ………………………….. ………………………….. …………….. 51
6.2.1 Rezolvitorul Navier -Stokes ………………………….. ………………………….. ………………………….. ……….. 52
6.2.2 Proced ura de integrare în timp ………………………….. ………………………….. ………………………….. …. 53
6.3 Alegerea schemei numerice ………………………….. ………………………… 54
6.4 Funcțiile subrutinelor ………………………….. ………………………….. ………. 55
7 Rezultate numerice ………………………….. ………………………….. ……. 56
8 Concluzii ………………………….. ………………………….. ………………….. 56
9 Bibliografie ………………………….. ………………………….. ……………….. 56
UNIVERSITATEA POLITEHNICA BUCURE STI
FACULTATEA DE INGINE RIE AEROSPATIALA
4
1 Introducere
Progresele majore ale tehnologiilor computerizate și numerice au făcut posibilă propunerea unei
alternative sau cel puțin o abordare complementară tehnicii clasice a experimentelor de laborator. Din ce
în ce mai mult, Dynamic FLuyd Computational, CFD devine parte a procesului de proiectare. De fapt,
chiar dacă nu elimină necesitatea de a face experimente prin realizarea de prototipuri, CFD -ul ar contribui
la reducerea numărului de cazu ri de testat .
Scopul acestei lucrări este acela de a studia interacțiunea undă de șoc -strat limită cu ajutorul
programului de calcul NSC2KE, având diferite date de intrare și diferite geometrii rezultând în final , ca și
concluzie cu o comparație între pol ara calculată a profilului aerodinamic NACA0012 și o polară calculată
cu alt program .
Programul de calcul NSC2KE a fost implemenentat de către Bijan Mohammadi -Stephane Lanteri în
anul 1994 . Acesta este un rezolver 2D și axialsimetric a ecuațiilor Euler și Navier -Stokes.
Când curgerea peste un corp are viteză supersonică, se formează undele de șoc, cauzate de : o
schimbare în panta suprafețe i, contrapresiunea care constrânge curgerea să devină subsonică, etc. În
aerodinamica modernă, se pot cita un număr ma re de circumstanțe în care undele de șoc sunt prezente [7].
Întâlnirea unei unde de șoc cu stratul limită rezultă într -un fen omen complex. Interacțiunea din tre
unda de șoc și stratul limtă, este adesea critică pentru performan ța vehiculelor. Șocul supune stratul limită
la un gradient de presiune negativ care poate altera puternic profilul de viteze. În același timp, în curgeri
turbulente, producția turbulentă este sporită lucru care amplifică disipația vâscoasă și conduce la sc ăderi
grave ale peroformanță.
O parcurgere a cuprinsului lucrării ilustrează cele mai sus mentționate. Așadar, după un capitol
introductiv privind ceea ce se va face până la finalul prezentei lucrări Capitolul doi reprezintă un capitol
introductiv privind aspecte fizice ale mișcărilor compresibile în cadrul căruia se introduc o serie de
noțiuni de bază și unele exemple clasice. În capitolul trei se prezintă modelele matematic e necesare
rezolvării problemei: deducerea ecuațiilor Navier -Stokes în diferite for me utile, mai ales pentru abordarea
numerică a aplicațiilor, ecuațiile mediate Reynolds pentru regimul compresibil și ecuațiile Euler unde s -au
prezentat proprietățile matematice ale acestora. Metodele numerice pentru rezolvarea ecuațiilor Navier –
Stokes s unt prezentate în capitolul patru. În cadrul acestui capitol sunt deasemenea prezentate metodele
cu volume finite și discretizarea în regim compresibil a fluxurilor convective și a celor difuzive. Necesar
calculului numeric este și un model de turbulență. Prin urmare, modelul de turbulență
k este
prezentat în capitolul cinci. Capitolul șase cuprinde caracteristicile programului de calcul NSC2KE (
ecuații gurvernatoare, tehnici numerice, scheme numerice, schema logică, etc). Ultimel e două capitole
sunt constituite din rezultatele numerice obținute și concluziile la care s -au ajuns.
La sfârșitul lucrării s -au introdus datele bibliografice. Am încercat, pe cât posibil să fac apel la
indicațiile bibliografice de origine.
Cuvinte chei e: ecuațiile Navier -Stokes, Modelul de turbulență
k , NSC2KE, ecuațiile mediate
Reynolds, ecuațiile Euler.
UNIVERSITATEA POLITEHNICA BUCURE STI
FACULTATEA DE INGINE RIE AEROSPATIALA
5 2 Aspecte fizice ale mișcărilor compresibile
Pentru mișcările în jurul unui corp in fluid nelimitat, domeniul tr ansonic poate fi limitat la un anumit
interval al numărului Mach,
M . Numărul Mach critic , este numărul Mach
M pentru care apare viteza
suntetului (M=1) în câmpul mișscării. Așadar, regimul de curgere este în regim subsonic daca
c MM ,
iar transonic dacă
c MM . Dacă
c MM dar
1 M atunci viteza este subsonic ă. În regiunile cu
1 M
va fi reprezentat ă o incluzicune supersonică adiacentă suprafeței corpului. Un astfel de regim
uneori poartă denumirea de regim transnic -subsonic .
Când
1 M apare o undă de șoc frontală întotdeauna detașată (fig.2.1 ., fig2.2 .). În avalul zonei
centr ale a undei frontale apare o zona subsonică, urmată de o zonă supersonică. În regiunea cu
1 M se
prezintă o incluziune subsonică , curentul fiind din nou supersonic departe de corp. Dacă profilul este
subțire și cu bord de atac ascuț it, la creșterea numărului Mach, unde de șoc se atașează la un anumit
număr Mach,
cM . Sub condițiile prezentate, pentru regimul de curgere devine în totalitate supersonic .
Figure 2.1.Unda de șoc frontală produsă de un obstacol rotunjit [2]
Figure 2.2.Curgerea în vecinătatea bordului de atac al unui profil supersonic[ 2]:
a. bord d e atac rotunjit
b. bord de atac ascuțit, diedru mare
c. bord de atac ascuțit, diedru mic
d. bord de atac puțin rotunjit[]
Rezultă deci că, pentru mișcările unui curent uniform la infinit, în prezența unui corp, regimul
transonic se întinde în vecinătatea va lorii
1M .
ccM M M
(2.1)
UNIVERSITATEA POLITEHNICA BUCURE STI
FACULTATEA DE INGINE RIE AEROSPATIALA
6 Regimul transonic conține atât particularități ale regimului subsonic cât și ale regimului
supersonic. Studiul regimului transonic este legat de o serie de dificultăți precum:
– Caracterul hiperbolic, eliptic al ecuațiilor mișcării in regiunile subsonice si supersonice;
– Decelerarea de la vite za supersonică la viteză subsonică se face brusc sub forma unor unde de
șoc;
– Influența efectelor de frecar e atât în mod direct prin producerea desprinderii stratului limită
datorită saltului de presiune produs de unda de șoc cât și în mod indirect prin modificarea
localizărilor liniilor sonice si îndeosebi a undelor de șoc;
Aceste dificultăți fac dificilă des crierea fenomenelor, facând aproape imposibilă elaborarea de
metode analitice. Așadar suntem conduși aproape exclusiv la necesitatea dezvoltării de metode numerice.
Considerații fizice
Pentru un profil simetric plasat în curent fără incidență, în regim incompresibil distribuția de viteze
si presiuni este simetrică față de mijlocul corzii iar la bordul de atac si bordul de fugă, viteza se anulează.
În regim subsonic, sunt menținute aceleași particularități dar cu deosebirea că efectele sunt mai
pronunțat e. Odată cu depășirea numărului Mach critic apare o regiune supersonică iar revenirea la
subsonic se face brusc printr -o undă de șoc. Pe măsură ce numărul Mach crește, incluziunea supersonică
se mărește ca întindere iar unda de șoc de deplasează către bord ul de fugă si crește în intensitate. Atunci
când numărul Mach depășește unitatea apare o undă de șoc frontală cu o incluziune subsonică în avalul
acesteia. În acest timp viteza crește devenind din nou supersonică iar la bordul de fugă se formează o
undă de șoc oblică. La numere Mach
c MM
, unda de șoc se atașează bordului de atac ascuțit iar
regimul de curgere poate fi considerat pur supersonic [2].
În regim transonic se observă modificarea ditribuției pres iunilor astfel încât, de la una simetrică în
regim subsonic la una antisimetrică în regim supersonic. Mișcarea în regim subsonic are loc pentru a
exista o interdependență între curgerea de pe intrados și cea de pe extrados. În schimb, în regim
supersonic c urgerea de pe intrados este complet dependentă de curgerea de pe extrados.
În concluzie particularitățile specifice regimului transonic afectează foarte mult caracteristicile
aerodinamice (rezistența la înaintare, portanța, etc).
Efecte de frecare
Efectele de frecare sunt întotdeauna prezente. În regimurile subsonice si supersonice au o pondere
realtiv redusă afecând distribuția presiunilor dacă stratul limită care se formează nu se desprinde și
producând o rzistență la înaintare deasemenea, redusă. To todată, rezistența la înaintare are ca și
componente rezistența de frecare și rezistența de formă . Conform paradoxului lui d ’Alambert rezistența
de formă este zero pentru un fluid perfect. În cazul fluidului real, aceasta poate fi importantă dacă stratul
limită se desprinde de suprafața corpului. În condițiile în care stratul limită nu se desprinde, curgerea din
exteriorul corpului este foarte puțin afectată de existența stratului limită. Așadar, rezistența de formă este
circa 20-40% din cea de frecare. Rez istența totală la înaintare poate include și rezistența indusă și
rezistența de undă [2].
Concluzionând, efectele de frecare ră man aceleași și în regim transonic deosbirile făcându -se
simțite datorită formării undelor de șoc.
Interacțiunea dintre o undă de șoc si stratul limită din apropierea peretelui este extrem de
complexă datorită faptului că viteza scade la zero pe corp devenind subsonică pe o anumită porțiune, deci,
unda de șoc trebuie sa se stingă undeva în interiorul s tratului limită.
În fig.2.3. sunt prezentate experimental două situații extreme astfel încât dacă stratul limită este
turbulent, unda de șoc se îngroașă în interiorul acestuia și se stinge în apropierea peretului. Dacă stratul
limită este laminar atunci, apare un efect al undei principale și puțin amonte în stratul limită, având
aspectul literei grecești
.
UNIVERSITATEA POLITEHNICA BUCURE STI
FACULTATEA DE INGINE RIE AEROSPATIALA
7 Datorită creșterii entropiei care duce la creșterea rezistenței, apariția undelor de șoc conduce la
pierderi. În regimuril e cu
1 M , creșterea rezistenței nu este explicabilă decât parțial. Aceaste devine
explicabilă doar dacă ținem cont de efectul undei de șoc asupra desprinderii stratului limită. Într-adevăr
unde de șoc este întotdeauna însoșită de un salt de presiune care dacă este îndeajuns de puternic, provoacă
desprinderea locală a stratului limită. Dacă condițiile din avalul undei de șoc (apar dacă gradientul de
presiune este negativ, constant sau pozitiv dar mic în valoarea absolută) sunt favorab ile acesta se poate
atașa din nou iar efectele asupra rezistenței la ănaintare rămân reduse.
În stratul limită turbulent stratul limită se desprinde mai greu așadar, in regim transonic, aceasta
apare în zona decelerată fiind urmată de o decelerare a cure ntului. Odată cu desprinderea stratului limită
sub aceste condiții nefavorabile, asupra rezistenței la înaintare apar efecte negative notabile. Deci, în
apariția undelor de șoc în regim transonic cu
1 M cel mai important este fenome nul de desprindere a
stratului limită.
O altă consecință a efectelor de frecare este alterarea mișcării exterioare a corpului. Condițiile la
limită sunt alterate și pot fi considerate de o manieră aproximativă fie prin alterarea conturului fie prin
alterarea efectivă a condițiilor la limită pe contur. Dacă stratul limită nu se desprinde atunci moficiările
menționate sunt mici și pot fi neglijabile. În regim transonic, aceste alterări pot modifica localizarea undei
de șoc. Odată cu aceasta, ditribuția pres iunilor este modificată cu efecte asupra portanței și asupra
rezistenței la înaintare.
Incluziuni supersonice
Regimul transonic -subsonic, din punct de vedere fizic se particularizează în principal prin apariția
unei regiuni în care viteza este supersoni că. O asemenea incluziune supersonică este mărginită în partea
de amonte de linia sonică iar în partea de aval de o undă de șoc (fig.2.4.a.).
În cele ce urmează se va demonstra posibilitatea ca la un anumit număr Mach să existe o frontieră
pentru care in cluziunea supersonică este continuu mărginită de o linie sonică, fără unde de șoc (fig.2.4.b.).
Această problemă a fost tratată matematic de către C. Morawetz care a arătat că o soluție continuă poate
exista pentru un anumit număr Mach dar că o asemenea soluție se comportă ca o soluție izoltă. Acest
rezultat explică de ce în practică, în mod obișnuit nu se întâlnesc soluții fără unde de șoc.
Evantaie de expansiune
Un evantai de expansiune supersonic, Prandtl -Meyer, este un proces de extindere central izat care
are loc atunci când o curgere supersonică se întoarce în jurul unui colț convex. Evantaiul constă dintr -un
număr infinit de unde Mach. Atunci când o curgere se întoarce în jurul unui colț neted, circular, aceste
unde pot fi extinse înapoi pentru a se întâlni într -un punct. Fiecare undă din evantaiul de expansiune
transformă treptat curgerea. Este imposibil din punct de vedere fizic ca curgerea să treacă printr -o singură
undă de "șoc", deoarece aceasta ar încălca a doua lege a termodinamicii. Prin intermediul evantaiului de
expansiune, curgerea se accelerează iar numărul Mach crește, în timp ce presiunea statică, temperatura și
densitatea scad. Deoarece procesul este izentropic, proprietățile de stagnare (presiunea totală și
temperatura totală) răm ân constante în lungul evantaiului [4].
Prima linie Mach este la un unghi
1
11arcsinM
iar ultima linie Mach este la un unghi
2
21arcsinM
Numărul Mach de după întoarcere
2M are legătură cu numărul Mach inițial
1M iar unghiul
de întoarcere
este:
21MM
(2.2)
UNIVERSITATEA POLITEHNICA BUCURE STI
FACULTATEA DE INGINE RIE AEROSPATIALA
8 unde,
M este funcția Prandtl -Meyer . Această funcție determină unghiul prin care o curgere sonică
trebuie să se întoarcă pentru a ajunge la un anumit număr Mach și are următoarea expresie:
2
22
21 1 1arctan 1 arctan 11 1112M dMM M MMM
(2.3)
max1121
(2.4)
iar unghiul maxim de întoarcere are următoarea expresie:
max max 1 M
(2.5)
3 Modele matematice
3.1 Ecuațiile Navi er-Stokes
Modelul Navier -Stokes este format din ecuația de continuitate, ecuația de impuls și ecuația energiei
completate cu alte ecuații.
În continuare vor fi analizate câteva reprezentări ale sistemului Navier -Stokes fiind introdu se
representăr i utile pentru analiza teoretică a sistemului de ecuații diferențiale.
3.1.1 Legi de conservare
Mișcarea unui fluid se cunoaște dacă la orice moment de timp pentru un fluid perfect
incompresibil se cunoaște o singură proprietate statică și anume presiunea, pe ntru un fluid perfect
compresibil in echilibru termodinamic se cunosc două proprietăți statice, presiunea si densitatea
și dacă se cunoaște câmpul de viteze [3].
Expresia matematică a principiilor generale de conservare (conservarea masei, impulsului
și energiei) cărora li se aduc completări cu relații constitutive constituie ecuațiile dinamicii
fluidelor .
În continuare este necesară introducerea unor noțiuni cu rolul de instrumente în mecanica
mediului continuu :
– Mărime a extensivă reprezintă orice mărime care își dublează valoarea dacă se
adiționează două volume identice(masa).
– Mărimea intensive păstrează aceeași valoare la adiționarea a două volume
identice(densitatea).
– Volumul de control constit uie o zonă mărginită de o suprafață perfect permeabilă(
suprafață de control ) la schimb de masa, impuls și energie și asupra căruia pot
acționa forte exterioare. În funcție de sistemul de referință la care îl raportăm,
volumul poate fi în repaus sau în mișcare.
– Volumul de material (sistem termodinamic simplu) se obține în cazul în care frontiera
domeniului spațial este perfect impermeabilă la orice transfer de masa, impuls sau
energie.
UNIVERSITATEA POLITEHNICA BUCURE STI
FACULTATEA DE INGINE RIE AEROSPATIALA
9 3.1.2 Teorema de transfer
Ecuația generală de transport sau transfer [3]:
PFt
(3.1)
P
– producția (sau sursa) mărimii extensive
pe unitatea de timp;
F – fluxul în unitatea de timp prin suprafața
;
Producția poate avea loc atât în interiorul volumului
cât și pe suprafața acestuia
.
P
poate reprezenta o su rsă care generează (pozitivă) sau una care cedează(negative).
Producția totală va fi suma contribuțiilor
q (sursă de volum) si
q (sursă de suprafață):
,, P q t x d q t x nd
(3.2)
Fluxul
F reprezintă schimbul mărimii
dintre volumul
și mediu. Așadar, acesta trebuie
sa conțină informații asupra geometriei
:
( , ) F J t x nd
(3.3)
Fluxul vector,
( , )J t x al mărimii
prin frontiera domeniului apare justificat. Datorită
procesului convectiv datorat mișcării macroscopice a fluidului și unui proces difuziv datorat mișcării
submacroscopice fluxul vector poate fi descompus in două componente.
( , ) ( , ) ( , )cd J t x J t x J t x
(3.4)
unde
( , )cJ t x reprezint ă componenta convectivă a fluxului iar
( , )dJ t x reprezintă componenta difuzivă
a acestuia.
Odată cu introducerea densității volumice
și considerării vitezei de transport macroscopică
V
fluxul convectiv se poate scrie:
,,cJ t x t x Vd
(3.5)
Fluxurile difuzive, de obicei se pot exprima prin intermediul unor legi empirice care exprima
fluxul conductiv în funcție de parametri macro scopici direct măsurabili.
Prin descompunerea (3 .4) se pune în evidență o superpoziție liniară a două fenomene fizice:
convecția macroscopică si agitația moleculară. Cele două fenomene sunt cuplate iar fluxul difuziv trebuie
considerat ca fiind o funcție i mplicită de câmpul de viteze.
Înlocuind relațiile (3.2) -(3.5) în ecuația generală de transfer (3 .1) obținem [3]:
, , , , ,d t x d q t x d q t x nd t x Vnd J t x ndt
(3.6)
Cu normala
n orientată spre interiorul volumului
, aplicăm teorema divergenței si obținem:
,, , , , 0dtxq t x q t x t x V J t x dt
(3.7)
Relația (2.7) poartă denumirea de forma integrală a teoremei de t ransfer .
Volumul
a fost ales arbitrar așadar, din relația (3 .7) rezultă:
UNIVERSITATEA POLITEHNICA BUCURE STI
FACULTATEA DE INGINE RIE AEROSPATIALA
10
,, , , , 0dtxq t x q t x t x V J t xt (3.8)
Relația (3 .8) reprezintă forma locală a teoeremei de transfer. În continu are vom particulariza ecuația (3 .8)
pentru următoarele mărimi extensive: masă, impuls și energie totală.
3.1.3 Conservarea masei
Notăm:
m și
–
m este msa;
–
este densitatea fluidului;
Înlocuind în relația
3.8 obținem:
,, , , , 0m m m
dtxq t x q t x t x V J t xt
(3.9)
Conform principiului conservării masei sursele de masă sunt nule :
,0mq t x
,0mq t x
,, , 0m
dtxt x V J t xt
(3.10)
Vor trebui scrise ecuații de conservare pentru fiecare componentă a amestecului și respectiv, pentru
amestec în totalitate. Pentru amest ec fluxul difuziv va fi nul deși, pentru o anumită componentă poate fi
nenul [3].
,0m
dJ t x
(3.11)
Odată cu înlocui rea acestui termen in relația (3 .10) obținem form a locală conservativă a ecuației de
continuitate :
,, , 0txt x V t xt
(3.12)
Forma neconservativă a ecuației de continuitate:
0DVDt
(3.13)
unde:
0t
Pentru mișcări incompresibile:
0V
(3.14)
UNIVERSITATEA POLITEHNICA BUCURE STI
FACULTATEA DE INGINE RIE AEROSPATIALA
11
3.1.4 Conservarea impulsului
Vom considera in ecuația (3 .8) ca mărime extensivă impulsul I[3]:
I Vd
(3.15)
0III
d V q q V V Jt
(3.16)
unde , mărimile
,I
q t x și
,I
dJ t x sunt mărimi tensoriale deoarece depinde de două mărimi vectoriale
(normala la suprafața volumului de control ,
n și viteza fluidului,
V ).
Conform legii fundamentale a mecanicii forța de impuls este cea ca re acționează asupra
sistemului.
Vom nota cu f densitatea forțelor exterioare sau de volum pe unitatea de masă, așadar:
,Iq t x f
(3.17)
Sursele de impuls distribuite pe suprafață vor fi repr ezentate de tensiunile distribuite pe suprafață:
,I
q t x
(3.18)
Conform ipotezei anterior admise,
,0I
dJ t x așadar, ecuația (2.17) devine:
0 V f V Vt
(3.19)
Ecuațiile de echilibru Navier vor fi completate cu relații constitutive numite postulatele lui Stokes.
Pentru un fluid newtonian expresia matematică a celor trei postulate este:
2s p V I
(3.20)
1; ; , 1,2,32j i
ij ij
jiu us s s i jxx
(3.21)
s
– tensorul vitezelor de deformație;
p
– presiunea statică;
I
-tensor unitar de ordinul doi;
– vâscozitate dinamică;
Pentru determinarea ecuației de transfer a impulsului, vom introduce următoarele relații [3]:
Se face apel la ipoteza Stokes pent ru determinarea parametrului
din ecuația (2.21):
1
3p Tr
(3.22)
unde
Tr este urma tensorului
. Așadar se obțin e o legatură între vâscozitatea dinamică și
:
2 3 0
(3.23)
UNIVERSITATEA POLITEHNICA BUCURE STI
FACULTATEA DE INGINE RIE AEROSPATIALA
12 Pentru un caz mai general, ipoteza Stokes se poate înlocui cu:
2
3
(3.24)
unde
reprezintă al doilea coeficient de vâscozitate.
Așadar, ecuația de t ransfer a impulsului este:
2203V f V V s p V It
(3.25)
Forma conservativă [3] se obține substituind tensorul tensiunilor totale
și tensorul tensiunilor
vâscoase,
:
pl
(3.26)
223s V I
(3.27)
V V V f plt
(3.28)
Vom înlocui in relația (2.29)
V pentru obținerea formei neconservative [3]:
DV VV V f pDt t
(3.29)
sau în coordonate carteziene:
223j i i i k
j i ij
j i j j i ku u u u u puft x x x x x x
(3.30)
Se mai p oate scrie și in formă echivalentă :
2
2VVV V f pt
(3.31)
Ecuațiile (3.28) -(3.31 ) poartă denumirea de ecuațiile Navier -Stokes .
3.1.5 Conservarea energiei
Suma dintre energia internă pe unitatea de masă a sistemului(e) și energia cinetică pe unitatea de
masă
2
2V
este energia totală pe unitatea de masă a unui sistem termodinamic, E [3].
2
2VEe
(3.32)
UNIVERSITATEA POLITEHNICA BUCURE STI
FACULTATEA DE INGINE RIE AEROSPATIALA
13 Ecuația (2.8) scrisă pentru energia totală a sistemului devine:
0EE
d E q q EV Jt
(3.33)
Conform principiului intâi al termodinamicii sursele volumice vor fi reprezentate ca fiind suma
dintre lucrul mecanic al forțelor exterioare ,
fV și sursele de caldură pe unitatea de masă ,
vq:
Eq fV q
(3.34)
unde lucrul mecanic al tensiun ilor reprezintă sursele distribuite pe suprafață:
EqV Vpl
(3.35)
Fluxul difuziv este reprezentat de Legea Fourier în absența difuziei masice:
,dJ t x k T
(3.36)
unde:
– k este conductivitatea termică;
– T este temperatura absolută;
Așadar eccuația de conservare a energiei se va scrie [3]:
E EV fV q k T Vt
(3.37)
Odată cu introducerea entalpiei totale , H și a entalpiei statice , h obținem următoarea formă a energiei:
E EV fV q k T Vt
(3.38)
2
2VpH h E
(3.39)
phe
(3.40)
Având în vedere ecuația (3.33) vom pune în evidență variația energiei interne a fluidului:
2
2DE De D V
Dt Dt Dt
(3.41)
Pentru a obține o ecuație pentru energia cinetică înmulțim scalar cu V ecuația de impuls (3.29):
DVV fV V p VDt
(3.42)
Dacă se introduc ecuațiile (3.41) și (3.42) în ecuaț ia (3.38) se va obține:
DeV q k T V VDt
(3.43)
unde suma ultimilor doi termeni reprezintă disipația
.
UNIVERSITATEA POLITEHNICA BUCURE STI
FACULTATEA DE INGINE RIE AEROSPATIALA
14
V V (3.44)
Ca atare, ecuația (3.43) se poate scrie:
DeV q k TDt
(3.45)
Disipația este întotdeauna pozitivă și pume în evidență ireversibilitatea transformării energiei
cinetice în căldură sub e fectul vâscozității.
În cazul în care fluidul este incompresibil iar conductivitatea termică este constantă, ecuația
anterioară se poate scrie:
Deq k TDt
(3.46)
3.1.6 Formularea conservativă locală
Ecuația de continuitate:
0 Vt
(3.47)
Ecuațiile de impuls:
VV V f pt
(3.48)
223d V I
(3.49)
Ecuația de energie:
E HV fV q k Tt
(3.50)
Pentru obținerea formei locale conservative în sistem cartezian, vom particulariza ecuația [3]:
y x zF F F UQt x y z
(3.51)
– U- vectorul variabilelor conservative;
TU u v w E
(3.52)
–
,,Fx Fy Fz – fluxuri;
2T
x xx xy xz xx xy xzTF u u p uv uw uH u v w kx
(3.53)
2T
y yx yy yz xy yy yzTF v uv v p uw vH u v w ky
(3.54)
2T
z zx zy zz zx zy zzTF w uw vw w p wH u v w kz
(3.55)
UNIVERSITATEA POLITEHNICA BUCURE STI
FACULTATEA DE INGINE RIE AEROSPATIALA
15
2 2 22 , 2 , 23 3 3xx yy zzu v w v u w w v u
x y z y x z z y x (3.56)
,,xy yx xz zx yz zyu v u w w v
y x z x y z
(3.57)
O scriere utilă pentru aplicațiile numerice a ecuatiei (3.51) se obține punând în evidență
componentele convectivă și d ifuzivă a fluxului:
x x y y z zUF G F G F G Qt x y z
(3.58)
unde
,,x y zG G G reprezintă componentele fluxului difuziv.
De obicei forțele masice și termenul sursă Q, sunt nule.
3.1.7 Formularea integrală conservativă
Aceasta se obține prin integrarea ecuației
UFQt pe un domeniu
urmată de aplicarea
teoremei de divergență pentru integrarea fluxului [3]:
Ud Fnd Qdt
(3.59)
unde,
reprezintă frontiera domeniului, iar n reprezintă normala la suprafața acestei frontiere.
3.1.8 Formularea neconservativă locală
Ecuația de continuitate:
0Vt
(3.60)
Ecuațiile de impuls:
DVfpDt
(3.61)
Ecuația de energie:
vDep V q k TDt
(3.62)
unde
reprezintă funcția de disipație.
UNIVERSITATEA POLITEHNICA BUCURE STI
FACULTATEA DE INGINE RIE AEROSPATIALA
16
3.1.9 Formularea tensorială
Ecuația de continuitate:
0i
iutx
(3.63)
Ecuațiile de impuls:
, 1,2,3i
i j i ij ij
jjuu u f it x x
(3.64)
2
3j ik
ij ij
j i ku uu
x x x
(3.65)
unde,
ij reprezintă tensorul tensiunilor vâscoase.
1
0ijij
ij
(3.66)
Ecuația de energie:
j j j v ij i
j j j jTe u H f u q k ut x x x x
(3.67)
3.1.10 Formularea adimensională
În continuare se va prezenta forma adimensionalizată a ecuațiilor Navier -Stokes coresp unzătoare
relației (3.58) [3].
, , , , , , ,tV x y z u v wx y z t u v wL L L L V V V
(3.68)
2 2, , , ,p T ep T eTV V
(3.69)
În ipoteza neglijării forțelor ex terioare ecuația (3.58) se va scrie:
0cx vx cy vy cz vzUE E E E E E
t x y z
(3.70)
unde :
–
U reprezintă vectorul variabilelor conservative adimensionalizate:
TU u v w E
(3.71)
UNIVERSITATEA POLITEHNICA BUCURE STI
FACULTATEA DE INGINE RIE AEROSPATIALA
17
–
,,cx cy czE E E sunt fluxurile convective adimensionalizate:
2
2
2T
cx
T
cy
T
czE u u p uv uw uH
E v uv v p uw uH
E w uv vw w p wH
(3.72)
– iar
,,vx vy vzE E E sunt fluxurile difuzive adi mensionalizate:
0
0
0T
xx xy xz xx xy xzvx x
T
xy yy yz xy yy yzvy y
T
xz yz zz xz yz zzvz zE u v w q
E u v w q
E u v w q
(3.73)
Entalpia totală adimensionalizată are următoarea formă:
pHE
(3.74)
Valorile adimensiona le ale energiei totale rezultă în forma:
2 2 2
2u v wEe
(3.75)
Tensiunile vâscoase adimensionale au următoarea formă:
2 2 22 , 2 , 23Re 3Re 3Rexx yy zz
L L Lu v w v u w w v u
x y z y x z z y x
(3.76)
,,Re Re Rexy yx xz zx yz zy
L L Lu v u w w v
y x z x y z
(3.77)
Parametrii adimensionali care apar în expresiile de mai sus:
ReLV L V L
(3.78)
VVMa RT
(3.79)
Prpc
k
(3.80)
UNIVERSITATEA POLITEHNICA BUCURE STI
FACULTATEA DE INGINE RIE AEROSPATIALA
18
3.1.11 Comentarii privind ecuațiile Navier -Stokes
Aspecte fizice
Un rol esențial în mișcarea unui fluid îl are raportul a două tipuri de forțe: forțele de frecare și
forțele de inerție . Numărul Re pune în evidență acest raport. În funcție de acesta se pot distinge trei
regimuri de curgere: laminar, turbulent, tranzitoriu [1].
Aspecte matematice
Ecua țiile Navier -Stokes constituie un sistem de ecuații cu de rivate parțiale de ordinul doi.
Termenul de inerție
VV introduce neliniaritatea . Acesta introduce interacțiuni complexe între
structurile de diferite scări, așadar este o sursă esențială de turbulență. Neliniaritatea poate fi considerată
slabă atâta timp cât forțele de inerție sunt mici în raport cu cele de frecare.
În cazul mișcării unui fluid guvernată de ecuații cu derivate parțiale neliniare, pornind de la
ecuațiile gene rale Navier -Stokes, integrarea este imposibilă deoarece numărul de variabile este infinit .
Din punct de vedere practic, în aerodinamică și în mecanica fluidelor, suntem interesați de
caracteristici ale curgerii cum ar fi: fluxul de căldură, coeficientul d e frecare, câmpul de presiuni medii,
vitezele medii, etc. Putem media statistic ecuațiile Nivier -Stokes, rezultând astfel, ecuațiile Reynolds care,
folosind modele de turbulență pot fi integrate numeric [1].
Aspecte numerice
Problema, în calculul numeric al curgerilor turbulenta constă în coexistența unor structuri
caracterizate prin scări de timp, de lungimi și viteze de diferite ordine de mărime. Pentru soluționare
numerică rețeaua de calcul va trebui sa aibă o rezoluție s pațială si temporală mai fină decât scările fizice
[1]. Din cauza numărului de pași ai grilei de calcul și ai pașilor de timp foarte mici, se ajunge la
concluzia că se depășește capacitatea celor mai puternice calculatoare din punct de vedere al timpului de
calcul și al memoriei necesare.
3.2 Ecuațiile mediate Reynolds
Principalul model de calcul utilizat pentru calculul mișcărilor în regim turbulent, îl
constituie, ecuaț iile mediate Reynolds . Termenii suplimentari ( tensiuni apa rente, fluxuri termice aparente )
rezultați din medierea ecu ațiilor Navier -Stokes sunt exprimați în funcție de parametrii medii prin
intermediul modelelor de turbulență. Pentru închiderea sistemului, modelele de turbulenșă introduc
ipoteze suplimentare.
Obținerea ecuațiilor mediate Reynolds rezultă din descompunerea parametrilor fluidului într -o
valoare medie și o fluctuație urmată de medierea sistemului de ecuații .
3.2.1 Medierea ecuațiilor Navier -Stokes
Fluide compresibile
Pentru obținerea ecuațiilor mediate R eynolds, vom utiliza operatorul de mediere ponderat masic.
Se vor defini următoarele valori medii:
UNIVERSITATEA POLITEHNICA BUCURE STI
FACULTATEA DE INGINE RIE AEROSPATIALA
19
uiui
(3.81)
HH
(3.82)
EE
(3.83)
și fluctuațiile:
,,i i iu u u H H H E E E
(3.84)
,,p p p T T T
(3.85)
Se evidențiază faptul că
0 pT dar mediile fluctuațiilor nu sunt nule așadar,
obținem:
0, 0, 0,iu H E
(3.86)
,,i
iuHEu H E
(3.87)
După introducerea descompunerilor în ecuațiile (3.63), (3.64), (3.67) și după mediere (statistică
sau temporală), obținem :
– Ecuația mediată de continuitate:
0i
iut t x
(3.88)
– Ecuațiile mediate de impuls:
, 1,2,3i i j ij ij i j
j i jpu u u u u it x x x
(3.89)
22,,33i j k j ikijij ij ij
k j i j i ku u u u u u
x x x x x x
(3.90)
În ecuația (3.74) termenul
ijuu are semnificația unei tensiuni aparente numită tensiune
Reynolds sau tensiune turbulentă
– Ecuația mediată a energiei:
j ij i iij ij i j
j j jTE u H u u u u H kt x x x
(3.91)
unde:
UNIVERSITATEA POLITEHNICA BUCURE STI
FACULTATEA DE INGINE RIE AEROSPATIALA
20
1 1 1 1,,2 2 2i i i iii E e u u k H h u u k k u u
(3.92)
Ecuația (3.94) se poate scrie sub forma:
1
2j i ii j j i i ij i ij
j j j jTE u H k u u u u u u u ut x x x x
(3.93)
Dacă se alege ca și necunoscută entalpia mediată masic , ecuația anterioară devine:
,i
jjj ij j
j j j j j jp p T u h h pu u u k u ht x t x x x x x
(3.94)
se poate alege ca necunoscută și temperatura mediată masic .
'
' i i
ij ij
jju u
xx
(3.95)
În concluzie, ecuațiile mediate Reynolds se vor scrie:
Ecuația de continuitate:
( ) 0i
iutx
(3.96)
Ecuațiile de impuls:
''( ) ( u ) ( )i i j ij i j
i i jpu u u ut x x x
(3.97)
Ecuația energiei:
''( ) ( u H) ( ' k )j ij i i ij j
i j jTE u u u Ht x x x
(3.98)
În aceste ecuații
, ,p,E,H,T,i iju reprezintă valorile medii ale mărimilor respective, iar cu „prim”
s-au notat fluctuațiile acestora. În cazul compresibil, sunt mediile ponderate masic, iar în c el
incompresibil, mediile statistice. De asemenea, fluctuațiile vor fi calculate întotdeauna în raport cu media
adoptată. Această scriere va pune în evidență asemănarea cu ecuațiile Navier -Stokes pentru valori
instantanee, precum și toți termenii supliment ari care apar.
3.2.2 Comentarii privind ecuațiile mediate Reynolds
Introducerea valorilor medii conduce la o netezire a parametrilor mișcă rii care permite calcul
numeric;
Datorită faptului că există o infinitate de valori instantanee care realizează aceeași med ie, valorile
medii nu mai prezintă o semnificație fizică bine determinată deci, este imposibilă regasirea configurației
reale a curgerii (valorile instantanee) daca se dispune doar de valorile medii;
Datorită apariției mai multor necunoscute ale sistemului de ecuații, sistemul de ecuații Reynolds nu
mai este un sistem închis iar pentru soluționare vom folosi modele de turbulență . În ecuațiile de impuls
apare divergența unor tensiuni aparente (turbulente) datorate transportului de impuls prin fluctuații
UNIVERSITATEA POLITEHNICA BUCURE STI
FACULTATEA DE INGINE RIE AEROSPATIALA
21 turbulente:
ij
juux , unde,
ijuu este definit ca fiin d tensorul tensiunilor turbulente . În
ecuația ecergiei apare o difuzie termică aparentă (turbulentă). Această difuzie turbulentă se adișionează la
difuzia moleculară , calculată cu valorile medii, facând schimbul de căldură într -o mișcare turbulentă față
de una laminară. Tot în ecuația energiei este prezentă o disipație turbulentă care pune în evidență
caracterul disipativ al mișcărilor turbulente.
Totodată se pot evi denția termenii care reprezintă interacțiunile complexe care se manifestă într -o
mișcare turbulentă: interacțiunea dintre câmpul fluctuațiilor vitezei și câmpul fluctuațiilor presiunii prin
termenul
j
jpux
și interacțiunea dintre câmpu l fluctuațiilor și vitezele de deformație prin termenul
ij iu
. Pentru integrarea ecuației energiei, este necesară modelarea acestor termeni.
3.3 Ecuațiile Euler
În aplicațiile practice inginerești prezintă interes curgerile la care term enii de convecție sunt cei
dominanți ( numere Reynolds mari). Acești termeni sunt incluși în ecuațiile Euler . Pentru partea
convectivă a ecuațiilor Navier -Stokes vor fi valabile aceleași metode numerice de integrare pentru
sistemul Euler, urmând ca termenii de difuzie sa fie discretizați separat.
Una dintre probleme actuale în calculul științific este aproximarea numerică a mișcărilor
compresibile multidimensionale. Mișcarea este modelată prin sistemul Euler, admițâmd ipoteza fluidului
fără vâscozitate. Ob iectivul de bază al calculului numeric la ora actuală este simularea numerică a
prezenței undelor de șoc, prin variații continue dar foarte rapide, fără a introduce oscilații fără
semnificație fizică în soluție.
3.3.1 Proprietăți matematice ale ecuațiilor Euler
Ecuațiile Euler formează un sistem neliniar de ecuații diferențial cu derivate parțiale de ordinul
întâi. Din formularea matematică a principiilor de conservare a fluidului fără vâscozitate pot fi deduse
ecuațiile Euler deoarece interpretarea fizică a aces tora o reprezintă chiar aceste principii.
În funcție de alegerea variabilelor independente sunt posibile mai multe reprezentări matematice
ale ecuațiilor Euler:
1. dacă se aleg ca necunoscute densitatea, componentele impulsului și energia totală atunci este o
reprezentare conservativă . Forma conservativă este esențială pentru a calcula corect viteza de
propagare si intensitatea undelor de șoc și a discontinuităților de contact, în fluidul fără
vâscozitate;
2. reprezentarea neconservativă constă în alegerea ca și necunoscute a densității, compo nentelor
vitezei și a presiunii. Forma neconservativă este de obicei aleasă în analiza teoretică a
proprietăților matematice ale modelului.
3. reprezentările în variabile caracteristice sunt posibile dacă necunoscutele sunt rep rezentate de
mărimi ce se conservă în lungul unor suprafețe caracteristice;
UNIVERSITATEA POLITEHNICA BUCURE STI
FACULTATEA DE INGINE RIE AEROSPATIALA
22 3.3.1.1 Formularea integrală
Pentru un volum de control
, mărginit de suprafața
ecuațiile de conservare se vor scrie:
1. Conserva rea masei:
0 d Vndt
(3.99)
2. Conservarea impulsului:
e Vd V V p nd f dt
(3.100)
unde
ef reprezintă forțele exter ioare.
3. Conservarea energiei:
e Ed HVnd f dt
(3.101)
Sistemul poate fi scris și intr -un sistem de referință mobil aflat intr -o mișcare de rotație cu viteza
unghiulară constantă
. În acest caz se va înlocui peste tot viteza V cu viteza relativă w unde:
w V x
(3.102)
Energia totală se va înlocui cu
E iar entalpia totală se va înlocu i cu rotalpia I.
Totodaă sistemul poate fi scris intr -o formă compactă introducând vectorul varibilelor
conservative U, vectorul flux conservativ, F și vectorul sursă Q. În final sistemul poate fi reprezentat în
forma:
Ud Fnd Qdt
(3.103)
Această ecuație trebuie completată cu ecuașia de stare care definește proprietățile termodinamice ale
fluidului.
3.3.1.2 Formularea diferențială conservativă
Reprezentare diferențială se obține aplicând teore ma diverg enței în sistemul (3.103 ):
y x zF F F UQt t t t
UFQt
(3.104)
Sistemul Euler în formă liniarizată :
UA U Qt
(3.105)
unde A, matricea iacobiană:
UNIVERSITATEA POLITEHNICA BUCURE STI
FACULTATEA DE INGINE RIE AEROSPATIALA
23
FAU (3.106)
Având în vedere componentele carteziene ale vectorului flux, sistemul Euler poate fi scris în
forma:
x y zUA U A U A U Qt x y t
(3.107)
dar,
0y x zA A AU U Ux y z
(3.108)
3.3.1.3 Matricele iacobiene în formularea conservativă
Pentru a putea calcula matricele iacobiene
,,x y zA A A este necesară exprimarea componentelor
vectorului flux
,,x y zF F F în raport cu vectorul conservativ U. Așadar, trebuie admisă o relație
constitutivă pentru fluid deci, vom considera că fluidul este perfect și satisface ecuația de stare
,vpRT e c T
.
1 2 3 4 5TTU u u u u u u v w E
(3.109)
2
3 2 4 2
1 2 3 4 5 2 2 2 5
1 1 1 1T
T T
x x x x x xu u u uF F F F F F u p u u u pu u u u
(3.110)
Similar se pot rescrie și celelalte componente ale vectorului flux
,yzFF . Pe baza acestora se pot calcula
elementele matricelor iacobiene:
1 1 1
1 2 5
2 2 2
1 2 5
5 5 5
1 2 5…
…
… … … …
…k k k
k k k
k
k
k k kF F F
u u u
F F F
Fu u u AU
F F F
u u u
(3.111)
– Mișcarea unidimensională. Vom considera că mișcarea are loc în direcția x:
2
320 1 0
3 3 12
1132kuA A u
u uE E u u
(3.112)
– Mișcarea bidimensională. Considerăm o mișcare în planul (x,y):
UNIVERSITATEA POLITEHNICA BUCURE STI
FACULTATEA DE INGINE RIE AEROSPATIALA
24
22
2 2 20 1 0 0
13 1 12
0
11 2 12xu V u v
Auv v u
u E V E V u uv u
(3.113)
22
2 2 20 0 1 0
00
11 3 12
11 1 22xuv v
A v V u v
v E V u E V v v
(3.114)
3.3.1.4 Formularea în varibile primitive
În continuare vom prezenta formularea în variabile primitive a ecuațiilor Euler deși pentru
dezvoltarea schemelor numerice se preferă exprimarea conservativă. Această prezentare a ecuațiilor Euler
este ultilă pentru obținerea unor scheme numerice de ordin superior de precizie. Așadar, expresi
sistemului Euler este:
– Ecuația de continuitate:
0 VVt
(3.115)
– Ecuația de impuls:
eVpV V ft
(3.116)
– Ecuația energiei:
0DepVt
(3.117)
Dacă se alege ca variabilă primitivă presiunea, ecuația energiei trebuie transformată intr -o ecu ație pentru
presiune:
20pV p c Vt
(3.118)
unde, c este viteza locală a sunetului.
Sistemul, în formă compactă are următoarea expresie:
pqB q Qt
(3.119)
în care, q reprezintă vectorul variabilelor primitive iar B este matricea iacobiană.
3.3.1.5 Legătura dintre formularea conservativă și cea în variabile primitive
Anterior s -a specificat faptul că pentru integrarea numerică a sistemului Eu ler se preferă
reprentarea conservativă, singura capabilă să capteze discontinuitățile care apar intr -o curgere de fluid
perfect. Dificultatea analizei proprietăților matematice a matricilor iacobiene este dată chiar de structura
acestora corespunzătoare r eprezentării conservative. În schimb, în reprezentarea în variabile primitive
acestea au o structură mult mai avantajoasă unor calcule analitice. Așadar, este de preferat să putem
UNIVERSITATEA POLITEHNICA BUCURE STI
FACULTATEA DE INGINE RIE AEROSPATIALA
25 transfera reultatele calculului analitic de la formularea în variabile prim itive la formularea în varibile
conservative. Acest lucru este posibil dacă vom considera o transformare neliniară de forma
() U U q .
Matricea iacobiană a transformării are următoarea formă:
UMq
(3.120)
care, în urma calculului analitic are forma:
15
15
552
151 0 0 0 0
… 0 0 0
0 0 0… … …0 0 0
… 1
21uuuvvvMwuu
Vvv u v w
(3.121)
M este o matrice nesingulară, iar inver sa ei pune în evidență că transformarea neliniară
() U U q este
ireversibilă :
21 0 0 0 0
10 0 0
10 0 0
10 0 0
1 1 1 1 12u
v
M
w
Vu v w
(3.122)
În variabile primitive , ecuațiile Euler au următoarea formă:
pqB q Qt
(3.123)
Ținând cont de matricea iacobiană a transformării putem scrie:
1
i i iq q U UMx U x x
(3.124)
Înmulțind cu matricea M la stânga, obținem:
1
pUMBM U MQt
(3.125)
Reprezentarea cvasiliniară a sistemului Euler în coordonate conservative este:
UA U Qt
(3.126)
În mod evident comparând relațiile (3.121) și (3.122) putem spun e ca matricea B are o structură
mult mai simplă decât matricea sistemului conservativ.
3.3.1.6 Variabilele caracteristice
Variabilele caracteristice ale sistemului Euler sunt date de relația:
1W L q
(3.127)
UNIVERSITATEA POLITEHNICA BUCURE STI
FACULTATEA DE INGINE RIE AEROSPATIALA
26 unde:
1 2 3 2
451, , ,
11,x z y x
ddw p w d w d u w d w d uc
w i V p w i V pcc
(3.128)
Variabila caracteristică
1w descrie perturbațiile de entropie , variabilele
23,ww descriu
undele de vârte juri. Ultimele două varibile caracteristici
45,ww corespund cu undele acustice.
Variabilele caracteristice W nu pot fi integrate în cazul general dar variația
W poate fi definită
întotdeauna.
Un mod mai simp lu de a exprima varibilele carcteristice constă în exprimarea directă a varibilelor
primitive prin cee conservative în relația (3.124).
În cazul unidimensional, varibilele caracteristice pot fi efectiv determinat și vor fi definite de:
1 2 3 21 1 1,, w p w u p w u pc c c
(3.129)
Sistemul caracteristic:
1 1 1 1
p L M A U L M Qt
(3.130)
devine:
11
22
330
0
0wwutx
wwuctx
wwuctx
(3.131)
Fig. 3.1 caracteristicile și conul Mach pentru mișcarea unidimenisonală [1]
În Fig.3.1. sunt prezentate caracteristicile
C și
C poartă numele de linii Mach și reprezintă înfășurarea
caracteristicilor iar
0C reprezintă o linie de curent care se propagă în lungul caracteristicii
C .
Putem interpreta starea fluidului într -un punct din domeniul mișcării unidimensionale ca fiind
rezultată din suprapunerea mărimilor transportate în lungul caracteristicilor. Așadar, intr -un punct P(x,t)
UNIVERSITATEA POLITEHNICA BUCURE STI
FACULTATEA DE INGINE RIE AEROSPATIALA
27 starea fizică este determinată de entropia transportaă în lungul caracteristicii
0C precum și de mărimile
transportate în lungul caracteristicii
C (fig.3.2.).
Fig. 3.2 caracteristicile într -o mișcare unimensională pentru o mișcare subsonică și o mișcare supe rsonică [1]
Dacă
0u
x și
0 ucx , în punctul M se intersectează două caracteristici de tipul
C
,situație prezentată în fig.3.3 .
Fig. 3.3 reprezintă formarea undelor de șoc în mișcarea unidimenională [1]
În această situație, în punctul M trebuiesc îndeplinite următoarele condiții:
112 2 2 2,1 1 1 1M P M Mu c u c u c u c
(3.132)
Acest lucru este imposibil, deoarece:
1122
11PMu c u c
(3.133)
Imposibilitatea de mai sus conduce la o discontinuitate de tip undă de șoc . Viteza C a undei de
șoc satisface relația:
11PMu c C u c
(3.134)
UNIVERSITATEA POLITEHNICA BUCURE STI
FACULTATEA DE INGINE RIE AEROSPATIALA
28 Această relație semnifică, din punct de vedere fizic, faptul că înformația care este transportată în lungul
caracteristicilor se propagă în lungul discontinuității.
În situația inversă :
11PMu c C u c
(3.135)
Informația se propagă de la discontinuitate în câmpul curgerii generând „unde de șoc de expansiune ”.
Pentru a se elimina acestea dintr soluțiile modelu lui matematic se va impune o restricție și anume condiția
de entropie.
3.3.1.7 Condiții la limită
Considerăm o suprafață constantă S și un punct al suprafeței acesteia P(x,y,z,t). Numărul
condițiilor la limită care necesită a fi impuse pe frontiera care are normal a n, vor coincide cu numărul
caracteristicilor asociate vectorului de undă d=n.
Pentru o curgere unidimensională, numărul de condiții la limită care trebuiesc impuse depinde de
modalitatea în care curbele caracteristice corespunzătoare normalei la supraf ață, intersectează frontierele.
Caracteristicile
0,CC , în secțiunea de intrare au pantele întotdeauna pozitive deci, vor transporta
informația de la frontieră către interiorul domeniului. Așadar, aceste mărimi transportate vor fi i mpuse
drept condiții la limită pe frontieră. Caracteristica,
C care depinde de numărul Mach la intrare, pentru o
intrare subsonică va avea panta negativă iar informația transportată iese din domeniu și mărimea
caracteristică asoci ată, nu trebuie dată ca și condiție la limită pe frontieră. Pentru o intrare supersonică va
avea panta pozitivă situația va sta exact invers.
Schemele numerice cer specificarea valorilor tuturor necunoscutelor pe frontiere așadar, pentru
realizarea sistem ului Euler se vor introduce condiții la limită numerice care se vor adăuga la condițiile
fizice fără a le influența și fără a intra în contradicție cu acestea.
3.3.1.8 Precondiționarea sistemului Euler
Precondiționarea are drept obiectiv definirea unor operatori matriceali sau diferențiali care să
îmbunătațească proprietățile formulării inițiale. Sistemul precondiționat reprezintă o aproximație a
formulării exacte.
Ecuațiile (3.115), (3.116), (3.117 ) reprezintă f orma neconservativă a sistemului Euler dar ecuația
energiei (3.117 ) mai poate fi scrisă sub forma:
phV p V htt
(3.136)
Introducând vectorii
M și
L, sistemul Euler poate fi scris ca:
0 LM
(3.137)
unde,
M poate fi exprimat în funcție de fluxul conservativ F. Astfel, ecuațiile de impuls scrise matriceal
sunt de forma:
M Z F
(3.138)
11 0 0 0 0
1 0 0 0
0 1 0 0
0 0 1 0
1u
Z v
w
H u v w
(3.139)
Sistemul neconservativ mai poate fi reprezentat sub forma:
UNIVERSITATEA POLITEHNICA BUCURE STI
FACULTATEA DE INGINE RIE AEROSPATIALA
29
0 L Z F (3.140)
Deoarece ,
1 UZLt
(3.141)
unde U reprezintă vectorul variabilelor constitutive, rezultă ca sistemul neconstitutiv poate fi scris sub
forma:
0UZ Z Ft
(3.142)
Variabilele primitive dependente vor fi alese ca fiind: presiunea, temperatura și componentele
vitezei. În urma definirii vectorului
Q și a matricii
, pent ru evaluarea matricii
se va considera
ecuația de stare intr -o formă în care sa includă atât regimul incompresibil cât și pe cel compresibil.
,pT
(3.143)
Așadar matricea iacobiană
calculată explicit are următoarea formă:
0 0 0
00
00
00
1pT
pT
pT
pT
p T puu
vv
qw
H u v w H c
(3.144)
Precondiționarea sistemului neconservativ
Ecuația de stare pentru un gaz p erfect este:
pRT
(3.145)
unde, derivata
p este:
21
pRT c
(3.146)
Pentru fluidul incompresibil
0p ceea ce înseamnă că viteza perturbațiilor acustice se
propagă la infinit iar din punct de vedere matematic anularea derivatei, duce la pierderea caracterului
hiperbolic al ecuației de continuitate. În vederea extinderii metodelo r numerice de integrare a sistemelor
hiperbolice și în cazul compresibil vom altera valoarea exactă a derivatei
p și cu aceasta, vom putea
controla viteza de propagare a perturbațiilor acustice.
Sistem ul Euler precondiționat se va scrie:
0PQS Z Ft
(3.147)
unde, matricea
S a fost înlocuită cu matricea
PS .
UNIVERSITATEA POLITEHNICA BUCURE STI
FACULTATEA DE INGINE RIE AEROSPATIALA
30 Reprezentarea conservativă precondiționată în variabile primitive
Această reprezentare se obține prin înmulțire a sistemului (3.147 ) la stânga cu matricea
1Z
(3.135). Ca urmare, sistemul va avea forma:
0PQZFt
(3.148)
unde,
1 PPZ Z S
(3.149)
Așadar sistemul (3.148 ) se poate scrie sub forma:
10P QZFt
(3.150)
Reprez entarea cvasiliniară a sistemului precondiționat
Reprezentarea cvasiliniară a sistem ului precondiționat (3.150 ) se poate obține dacă relația:
1F Z M
(3.151)
se înlocuiește în (3.146):
10P QSMt
(3.152)
unde,
11
1
1
1
1P
pV V p
pVux
pVv SMy
pVwz
V V Tc
(3.153)
După înlocuirea produsului (3.153) în sistemul (3.151 ) se obține:
0 x y zQ Q Q QA A At x y z
(3.154)
unde,
,,x y zA A A se pot identifica din relația (3.153 ).
UNIVERSITATEA POLITEHNICA BUCURE STI
FACULTATEA DE INGINE RIE AEROSPATIALA
31 Valorile proprii ale sistemului precondiționat
Matric ea iacobiană a sistemului (3.154 ) are următoarea formă :
0
0 0 0
0 0 0
0 0 0
0y xz
x
y
z
y xz
p p pd dddV
ddV
ddV K
ddV
d dddVc c c
(3.155)
Soluțiile ecuației:
det 0KI
(3.156)
reprezintă valorile proprii ale matricei
K . Rezolvând această ecuație vom avea o rădacină multiplă de
ordinul trei:
1 2 3 dV
(3.157)
iar radăcinile ecuației:
2
20ddV
(3.158)
sunt
4 și
5 care reprezintă vitezele de propagare a undelor acustice în fluid.
Putem calcula parametrul
pentru cazul compresibil:
21 1 1P T
p
ppc TR Tc TR TR c
(3.159)
Înlocuind în relația (3.158 ) rezultă că:
4,5 ; dV dc
(3.160)
Există tendința de a se aproxima regimul incompresibil cu o mișcare compresibilă la vi teze mici.
În această situație
4,5 rămân finite dar cu valori absolute foarte mari față de valorile proprii. În mod
obișnuit se obține
1,2,3 2
4,510
care ne spune că undele acustice devin dominante în comparație cu cel e
de entropie și vorticitate sau condiția de stabilitate a algoritmilor conduc la pași de timp foarte mici
(număr foarte mare de iterații). Această situație poate apărea chiar în cazul mișscărilor compresibile( de
exemplu, în vecinătatea undelor de șoc).
Pentru a limita viteza de propagare a perturbațiilor acustice vom considera o referință maximă de
propagare a acestora, astfel încât:
()
10 ( )
mod ( ). c la viteze mari regim transonic
V la viteze mici regim incompresibil
Vd la viteze erate regim compresibil subcritic
(3.161)
Astfel, ecuația (3.154), devine:
UNIVERSITATEA POLITEHNICA BUCURE STI
FACULTATEA DE INGINE RIE AEROSPATIALA
32
2 220 Vd d V (3.162)
Vom introduce un factor de compresibilitate
și o funcție comutatoare
pentru a scrie unitar
cele trei cazuri de mai sus:
21,
0,T
p
pîn regim compresibilccîn regim incompresibil
(3.163)
2 0,1
12 ,2în regim compresibilV
în regim incompresibil
(3.164)
Putem exprima valorile proprii
4,5 sub forma:
22
4,521 Vd V V d
(3.165)
În regim compresibil, ăentru cazul vitezelor moderate în care se consideră
V Vd putem lua
0
.
3.3.1.9 Problema Riemann
Problema Riemann mai este numită și problema tubului de șoc și are o mare importanță pentru
intregrarea sistemu li Euler. Putem construi o soluție ex actă a modelului Euler unidimensional care este
intens folosită. Cazul prezentat de aceasta ne permite să analizăm comportarea schemelor nume rice în
prezența undelor de șoc, a evantaielor de expansiune și a discontinuităților de contact, simultan iar baza
unor scheme numerice o constituie soluția exactă sau aproximativă pentru mișcări multidimensionale.
Fig. 3.4 Tubul de șoc la momentul inițial [1]
Odată cu spargerea diafragmei la momentul t=0, se formează o undă de șoc , care se propagă în
zona cu presiune scăzută (1) simultan cu producerea evantaiului de expansiune în zona cu presiune
ridicară (4).
UNIVERSITATEA POLITEHNICA BUCURE STI
FACULTATEA DE INGINE RIE AEROSPATIALA
33
Fig. 3.5 Schema mișcării în tubul de șoc [1]
Dim moment ce unda de șoc se deplasează către dreapta, aceasta va accelera fluidul în sens
pozitiv până la viteza
2u . Evantaiul de expansiune, se deplasează către stânga și accelerează fluidul până
la viteza
3u . Așadar frontul evantaiulii de expansiune se va deplasa cu viteza suntetului
4c , iar ultima
undă din evantai se va deplasa cu viteza
3c . Deoarece entropia fluidului este diferită în fața undei si în
fața evantaiului de expansiune, deosebim două zone separate printr -o discontinuitate de contact .
În urma prezentării analiti ce a problemei Riemann (Anderson, Saad), obtinută considerându -se o
curgere unidimensională de fluid fără vâscozitate care satisface ecuația de stare a gazului perfect, se obțin
următoarele relații pentru:
– Densitatea, temperatura, viteza undei de șoc, vite za în spatele undei de șoc:
12
1
111112pCcp
(3.166)
1 2 1 1 2
2
1 1 1 1 121111c p pupp
(3.167)
11
2 2 1 2 1 2 2 1 2 1 2
1 1 1 1 1 1 1 1 1 1 1 11 1 1 11 , 1 ,1 1 1 1sT p p p p p CMT p p p p p c
(3.168)
Pentru calculul curgerii în interiorul evantaiului de expansiune vom reaminti ecuația
caracteristicii
C :
dxucdt
x u c t
(3.169)
de unde rezultă:
42
1xuct
(3.170)
Se consideră gazul perfect ca fiind
RTc rezultă:
2
44112Tu
Ta
(3.171)
UNIVERSITATEA POLITEHNICA BUCURE STI
FACULTATEA DE INGINE RIE AEROSPATIALA
34 iar cu ajutorul relațiilo r evoluției izentropice obținem:
22
11
4 4 4 4111 , 122u p u
a p a
(3.172)
3.3.1.10 Matricele iacobiene în variabile primitive
Vom scrie explicit sistemul Euler în variabile primitive, într -un sistem cartezian de referință :
x y z pq q q qB B B Qt x y z
(3.173)
unde,
, , ,x y zB B B se deduc prin identificare.
Pentru cazul unidimensional vectorul conservativ și matricea iacobiană poat avea următoarea
formă:
, qu
p
(3.174)
20
0 1/ .
0xu
Bu
cu
(3.175)
Pentru cazul bidimensional, acestea vor fi:
,uqv
p
(3.176)
200
0 0 1/,0 0 0
00xu
uBu
cu
(3.177)
200
0 0 0.0 0 1/
00yv
vBv
cv
(3.178)
3.3.1.11 Valorile proprii ale sistemului Euler
Se caută soluția de forma u nei unde:
ˆ ,i dx tq qe
(3.179)
unde,
ˆq este amplitudinea complexă a undei, d vectorul de undă,
pulsația undei iar x este vectorul de
poziție.
Faza undei care s epropagă în direcția d cu pulsația
este reprezentată de relația:
UNIVERSITATEA POLITEHNICA BUCURE STI
FACULTATEA DE INGINE RIE AEROSPATIALA
35
,S x t dx t (3.180)
Matricea iacobiană a sistemului este de forma:
,p x x y y z zK Bd B d B d B d
(3.181)
iar
, , ,x y zd d d sunt parametrii nespecificați ai vectorului de undă d.
Valorile și vectorii proprii ai matricei iacobiene sunt:
ˆˆ ,
.pK q q
(3.182)
Ecuația caracteristică, rezultă din condiția ca sistemul să fie omogen și să admită o soluție diferită
de cea banală:
det 0.pKI
(3.183)
4 Metode nume rice pentru rezolvarea ecuațiilor Navier -Stokes
Discretizarea termenilor din ecuațiile Navier -Stokes se face diferențiat în raport cu semnificația
fizică a acestora pentru curgerile compresibile turbulente, unde efectele de inerție sunt dominante.
Ecuați ile Navier -Stokes pentru problemele care conțin interacțiuni complexe, zone de separare sau
domenii de calcul ce conțin zone subsonice și supersonice, în ceea ce privește soluția inițială, aceasta se
alege arbitrar iar ecuațiile sunt integrate în timp pent ru obșinerea unei solușii staționare. Așadar, soluțiile
intermediare, nu au semnificație fizică acestea constituind doar condițiile inițiale pentru iterația care
urmează. Pentru mărirea eficienței calculului, se impune ca pasul de timp impus de condiția de stabilitate
sa furnizeze cea mai mare valoare posibilă. În cazul obținerii unei soluții nestaționare, pasul te timp, este
determinat nu numai din impunerea condiției de stabilitate ci și de condiția ca acesta să fie compatibil cu
fenomenul fizic dar și fu rnizarea unor condițiiinițiale cu semnificație fizică.
Pentru rezolvarea sistemului Navier -Stokes se vor enumera câteva observații privind cerințele ce
trebuiesc îndeplin ite:
– pentru integrarea numerică a sistemului Navier -Stokes, grila de cal cul treb uie sa aibă o rezoluție
spațială mult mai mare față de cea utilizată pentru integraea sistemului Euler. Gradientul câmpului de
viteze și temperaturi, în vecinătatea peretelui, dacă numărul Reynolds este mare, este mai mare pe direcția
normalei la perete de cât în lungul acestuia. Pentru a se putea reproduce fenomenul fizic de către soluția
numerică, pasul grilei în direcția normale la perete trebuie să fie extrem de mic. În regim turbulent, grila
va fi mult mai fină în vecinătatea pereților decât în regim la minar. Această deformare a grilei, are un efect
negativ asupra schemelor numerice precum condițiile de stabilitate sunt foarte dure iar ordinul de precizie
este alterat drastic. Prin urmare, se impune utilizarea unor grile de tip structurat, pretabile tran sformării
geometrice pentru integrarea sistemului.
– fluxurile nevâscoase (convective) datorită faptului că sunt cuprinse în întregime în sistemul
Euler, discretizarea acestora se poate face cu scheme centrate, scheme upwind, etc. Fluxurile difuzive,
pentru a fi în concordanță cu simetria spațială a fenomenului fizic pe care -l reprezintă se vor discretiza
centrat.
– termenii de transport convectiv, pentru integrarea sistemului Navier -Stokes se vor discretiza cu
scheme de cel puțin oridnul doi de precizi e. În discretizarea fluxurilor convective, dacă se utilizează
scheme de ordin superior, este necesară introducerea unei vâscozități artificiale.
– schemele implicite sunt preferate în cazul în care dorim obținerea unor soluții staționare dar ele
sunt mul t mai stabile decât cele explicite. În probleme multidimensionale, pentru mărirea eficienței
calculului numeric, se va aplica o procedură de decuplare a direcțiilor spațiale sau de alternare a
direcțiilor.
UNIVERSITATEA POLITEHNICA BUCURE STI
FACULTATEA DE INGINE RIE AEROSPATIALA
36 4.1 Metode cu volume finite
La bază, aceste metode au integrarea directă, în spațiul fizic, a formulărilor conservative ale
ecuațiilor și sistemelor de ecuații diferențiale cu derivate parțiale [1]. Acestea pot fi aplicate atât pe grile
structurate cât și pe grile nestructurate . Necunoscutele, reprezintă o valoare medie pe celula asociată unui
anumit punct.
Metodele cu volume finite, pe de o parte, pentru cazul grilelor nestructurate pot fi privite ca o
extindere a metodelor cu diferențe iar pe de altă parte, ca o variantă a form ulării integrale din metoda
reziduurilor ponderate [1].
4.1.1 Discretizări conservative
Prin integrarea formală pe un volum de control în jurul unui nod, al rețelei, se caută ca soluția să
respecte condițiile la limită asociate și să aibă un flux diferențiabil . Acest lucru smnifică faptul că fluxul
este continuu la trecerea dintr -un volum de control în altul asigurând păstrarea caracterului conservativ al
ecuației discrete. Flexibilitatea metodei este dată de posibilitatea de a se modifica forma și localizarea
spațială a volumelor de control asociate unui punct din grilă dar și posibilitatea de a mări precizia
evaluării fluxurilor prin interfețele volumelor de control . Proprietatea necesară pentru orice schemă
numerică, conservarea la nivel discret a masei, impulsului și energiei este asigurată prin discretizarea
directă a legilor de conservare) [ 1].
Surse numerice
Ecuația conservativă de transport:
()0.u f u
tx
(4.1)
Forma conservativă este echivalentă cu forma cvasilinară:
( ) 0,uuautu
(4.2)
unde:
( ) .fuauu
(4.3)
Discretizarea cu diferențe finite centrate a ecuației (4.1) în trei puncte succesive este:
1/ 2 1/ 2 3/ 2 1/ 2 1/ 2 3/ 2
110, 0, 0.i i i i i i
i i iu f f u f f u f f
t x t x t x
(4.4)
Însumând, obținem:
3/ 2 3/ 2
110.ii
i i iu u u f f
t t t x
(4.5)
Discretizând centrat pe celula
3/2 3/2,iixx , obținem:
1 1 3/2 3/20.3i i i i i
iu u u f f
tx
(4.6)
Relațiile (4.5) și (4.6) sunt identice deoarece fluxurile în punctele interioare s -au anulat.
Aceeași discretizare dar pentru forma neconservativă rezultă:
UNIVERSITATEA POLITEHNICA BUCURE STI
FACULTATEA DE INGINE RIE AEROSPATIALA
37
1/ 2 1/ 2 3/ 2 1/ 2 3/ 2 1/ 2
11
110, 0, 0.i i i i i i
i i i
i i iu u u u u u u u ua a at x t x t x
(4.7)
Dacă însumăm relațiile obținem:
1 3 3/ 2 3/ 2 3/ 2 3/ 2
3/ 2 3/ 2 1/ 2 3/ 2
1/ 2 3/ 2
3/ 2 1/ 2( ) ( )3 6 6
.6i i i i i i i
i i i i
ii
iiu u u u u u ua a a at x x
uuaax
(4.8)
Prin discretizarea directă pe celula
3/2 3/2,iixx se obține:
1 1 3/ 2 3/ 2
3/ 2 3/ 2( ) 0.36i i i i i
iiu u u u uaatx
(4.9)
Aceste două relații nu mai sunt identice. Semnificația unui termn datorat doar discreti zării numerice se
regăsește în membrul drept al ecuației (4.8).
Testele numerice pun în evidență că rezultatele obținute pe formulări conservative au o acuratețe
superioară celor obținute pe formulări neconservative [1].
Fluxul numeric și convergen ța
Se consideră
1/ 2if
o formulă numerică de aproximare a funcției
fu iar dependența o
considerăm de forma:
**
1/ 2 1 ( , , , ),i i i i kf f u u u
(4.10)
Așadar, schema conservativă (4.4) va fi implementată numeric sub forma:
**
1/2 1/20.ii
iu f f
tx
(4.11)
Fluxurile numerice îndeplinesc condiția:
1 *( , , , ) ( ).i i i ku u u u f u u u f u
(4.12)
dacă schema este consistentă, iar:
11
001
0( , ) ( , ) .nnnn
nnxx tt
nn
n
xx ttu dx u dx f x t dt f x t dt
(4.13)
Această teoremă garantează că atunci când soluția converge, aceasta o va face către o ecuație de
conservare pe continuu. Ea nu garantează că schema numerică este convergentă sau consistentă [1].
Aproximații integrale
Schema numerică (4.11) poate fi obținută și cu volume finite:
**
1/2 1/210.i
iiduffdt x
(4.14)
Totuși, între cele două scheme există o deosebire. În schema (4.11),
iu reprezintă valoarea
funcției u în nodul i iar în schema (4.14),
iu reprezintă o valoare medie a funcției u pe celul asociată
punctului, i fără a mai fi atribuită unui anumit punct al celulei [1].
UNIVERSITATEA POLITEHNICA BUCURE STI
FACULTATEA DE INGINE RIE AEROSPATIALA
38
Schemele numerice cu volume finite sunt mult mai ușor de obținut decât schemele cu diferențe
finite doarece este suficient s ă integrăm formal ecuația de conservare pe un subdomeniu al domeniului de
integrare. În dinamica fluidelor, această simplitate de obținere este mult mai avantajoasă tinând cont de
faptul că forma naturală a ecuațiilor de mișcare este cea integrală. Cu toa te că principala problemă care
apare într -o discretizare cu volume finite este estimarea valorii fluxului la interfața dintre două celule
vecine, mai pot apărea și alte cerințe suplimentare (pentru cazul multidimensional, asigurarea conservării
fluxurilor la nivelul interfețelor) [1].
4.1.2 Scheme cu volume finite
Dacă vom diviza domeniul, într -o serie de subdomenii:
) ( ( ) ,
jj j j jdU F U nd Qdt
(4.15)
unde:
11,,
jjjj
jjU Ud Q Qd
(4.16)
Schema numerică cu volume finite are forma generică:
*
) ( ( , ) ,
jj j jL jR j jdU F U U nd Qdt
(4.17)
iar,
,jL jRUU reprezintă valorile
jU de o parte și de alta a frontierei.
Definirea grilei și a volumelor de control
Metodele cu diferențe finite pot fi aplicate în grile structurat e sau în grile nestructurate. Pentru
cazul grilelor nestructurate, se impun insă o serie d e limită ri suplimentare . După ce am selectat rețeaua,
avem două variante de a defini volumele finite : centrat pe celula sau centrat pe nod .
Fig. 4.1. Tipuri de volume finite [6]
Datorită faptului că schema numerică în varianta centrată pe nod furnizează valorile medii ale
necunoscutelor pe un volum finit, considerate reprezentative pentru toate punctele acestuia, obținem
valorile necunoscutelor în nodurile grilei de bază . Valori le necunoscutelor în varianta centrat pe nod se
vor interpola funcție de valorile calculate în celulele care îl înconjoară [1].
UNIVERSITATEA POLITEHNICA BUCURE STI
FACULTATEA DE INGINE RIE AEROSPATIALA
39 Fluxul numeric prin fețele celulei
În metodele cu volume finite, principala problemă o consti tuie evaluarea fluxului numeric la
interfața dintre două celule de integrare. Fluxul numeric la nivelul interfeței se poate calcula prin două
variante: centrat si upwind.
Pentru introducerea unor expresii ale fluxului numeric, considerăm relația de conser vare locală pe
un domeniu bidimensional:
( ) ,x x y yij ABCD ijdUdxdy F n F n ds Qdxdydt
(4.18)
,.xyF Fi F Fj
(4.19)
( ) ,xyij ABCD ijdUdxdy F dy F dx Qdxdydt
(4.20)
Dacă aplicăm teorema pentru calculul integrarelor flux vom avea:
,, ( ) ,x AB AB y AB AB ij ij
ABCDdU F y F x Qdt
(4.21)
Sheme centrate
Evaluarea centrată a fluxului presupune adaptarea unei discretiz ări numerice care să
îndeplinească condițiile de simetrie.
Pentru o discretizare centrată putem lua în considerare variantele următoare:
– Medierea fluxurilor calculate cu valorile necunoscutelor din celulele vecine:
, 1,( ) / 2,AB i j i jF F F
(4.22)
– Medierea necunoscutei între celulele vecine:
, 1,.2i j i j
ABUUFF
(4.23)
Sheme upwind
Principiul de bază constă în considerarea fluxului numeric prin una dintr e fețele celulei egal cu
fluxul generat de „celula din amonte” [1]. Principiul upwind intuiește că fluxul printr -o interfață a celulei
este complet determinat de valoarea transportată în direcția vitezei de convecție.
Aplicân d schema (4.21) pe o grilă carteziană vom avea:
, , , , , , 1, 1,
, , , , , , , 1 , 1,,
,,x AB AB x i j x i j x CD CD x i j x i j
y BC BC y i j y i j y DA DA y i j x i jF y F y a u y F y F y a u y
F y F x a u x F y F x a u x
(4.24)
iar ecuația discretă este:
, , , , 1, , , , , 10.i j x i j x i j y i j y i jdu F F F F
dt x x
(4.25)
Gradientul necuno scutelor
Chiar dacă în schemele cu volume finite necunoscuta numerică reprezintă o valoare medie pe o
celulă de integrare a necunoscutei fizice, pot exista și situații în care este necesară estimarea „gradientulul
mediei”. Totodată, și în schemele de or din superior pentru ecuația de transport
convectiv, se impune determinarea gradienților medii pe o celulă de integrare.
În discretizarea cu volume finite, medierea numerică a derivatelor variabilelor
poate fi făcută prin intermediul teoremei lui Gauss:
Ud Und
(4.26)
UNIVERSITATEA POLITEHNICA BUCURE STI
FACULTATEA DE INGINE RIE AEROSPATIALA
40
4.2 Discretizarea în regim compresibil
Ecuațile Navier -Stokes pot fi soluționate pentru a găsi soluții staționare sau nestaționare. Soluția
inițială se alege de obicei arbitrar p entru probleme care presupun zone de separare, interacțiuni complexe
sau domenii de calcul ce conțin zone subsonice și hipersonice. Pentru obținerea unei soluții staționare,
ecuațiile sunt integrate în timp. Ațadar, în acest caz pentru mărirea eficienței d e calcul, se va urmări ca
pasul de timp impus de condiția de stabilitate să capete valoarea cea mai mare posibilă.
Pentru a obține soluții nestaționare, condițiile inițiale, tebuiesc formulate cu o semnificație fizică
iar determinarea pasului de timp va r eieși nu numai din impunerea condiției de stabilitae ci și din condiția
impusă ca acest sa fie compatibil cu fenomenul fizic [1].
4.2.1 Discretizarea fluxurilor convective
Discretizarea fluxurilor convective se poate face prin:
– decuplarea vectorului flux;
– scheme centrate cu adționarea unei vâscozități artificiale;
– scheme TVD, ENO sau WENO;
În cazul utilizării unor scheme de ordin superior pentru discretizarea fluxurilor convectiveeste
necesară introducerea unei vâscozități artific iale.
4.2.1.1 Scheme upwind de tip Godunov
Formularea generală a schemelor de tip Godunov
Forma generală (pentru un sistem hiperbolic unidimensional) a unei scheme de bază Godunov
este:
1
111 1 1 1
2 2 2 2, 1, , , 1, ,n n n n n nR R R R R R
i i i i i i
i i i itU U F U F U U U U U U U U Ux
(4.27)
Datorită faptului că soluționarea exactă a problemelor Riemann numerice este costisitoare iar
detalii incluse în soluția exactă dispar prin mediere, se poate înlocui soluția exactă cu o soluție
aproximativă. Solvitorii Riemann aproximativi sunt metode de construcție a solțiilor aproximative ale
problemelor Riemann. Cele mai utilizate sunt metode Osher și metoda Roe.
Schema Roe
Rezolvitorul Roe are la bază decuplarea după direcțiile caracteristice a variațiilor fluxurilor.
Metoda Roe poate fi o altern ativă (baz ată pe o liniarizare locală dup ă direcțiile caracteristice) de calcul a
fluxului numeric. Schema Roe este una dintre cele mai utlizate scheme pentru integrarea sistemului Euler.
Reprezentarea conservativă asociată schemei upwinnd pentru sisteme hiperbolice liniare
este:
1
1/2 1/2nn
i i i iU U F F
(4.28)
În continuare, vom prezenta algoritmul general de implementare a rezolvitorului Roe intr –
o discretizare cu volume finite.
Exprimarea conservativă integrală a sistemului Euler, scrisă pe o celulă de calcul
k ,
conduce la sistemul diferențial:
kkdU tF U nddt
(4.29)
UNIVERSITATEA POLITEHNICA BUCURE STI
FACULTATEA DE INGINE RIE AEROSPATIALA
41 unde:
1
kk
eU t U t d
(4.30)
Următorul sistem difer ențial ordinar a rezultat în urma aplicării ecuației (4.28) fiecărei
celule:
1, 1,2…k
k
kdUR U kdt
(4.31)
iar,
kRU este inte grala fluxului F pe forntiera celulei k.
Frontiera unei celule în cazul bidimensional este reprezentată de segmente de dreaptă, așadar,
kRU
se aproximează:
11kk
ll
lNN
kn
llR U F U nd F U d
(4.32)
Vom nota:
lllF F U d
(4.33)
Pentru a estima fluxul
lnF putem aplica variante de decuplare sau schema Roe. Fluxul F este o
funcție omogenă de gradul întâi de vecto rul conservativ U:
F AU
(4.34)
unde, A este matricea iacobiană iar componentele carteziene:
nF Fn nAU KU
(4.35)
Observația pr ecum că local, la nivelul unei interfețe apare o similitudine cu cazul unidimensional
dacă se identifică o direcție x din cazul unidimensional cu direcția normalei n. Această observație ne
permite înlocuirea fluxului numeric pe interfață cu fluxul numeric Roe:
,
l lRLnF U Fn U U
(4.36)
unde
,
lRL Fn U U
poate fi calculat cu :
lk
knk R
kF Fn r
(4.37)
iar,
k reprezintă valorile proprii,
kr reprezintă vectorii proprii și
k reprezintă amplitudinile
vectorilor proprii ale matrice Roe. Matricea Roe, va avea aceeași formă cu matricea iacobiană K, în care
entalpia, densitatea, și componentele vitezei sunt înlocuite cu valorile acestora mediate masic.
Relația (4.31) devine:
11,kk
lNN
l RLk n l
llR U F F U U
(4.38)
Pentr u integrare putem folosi o metodă de tip Euler, care conduce la schema:
1
1,.k
lNnn
k k R Lnl
lktU U F U U
(4.39)
UNIVERSITATEA POLITEHNICA BUCURE STI
FACULTATEA DE INGINE RIE AEROSPATIALA
42 Schema Osher
Osher, propune, pentru fluxul numeric Godunov, o formulare exactă [1]. Schema numerică
Godunov în funcție de fluxul numeric, poate fi:
1
11
1/ 2
1( ),
( ),min
maxii
iiii
u u u
i
ii
u u uF u pentru u u
F
F u pentru u u
(4.40)
Dacă funcția
Fu este convexă atunci relația de mai sus arată că fluxu l numeric este monoton
pe intervalul
1,iiuu .
Ambiguitatea schemei de bază la traversarea punctului critic este rezolvată atunci când relația
anterior menționată este valabilă chiar dacă derivata
dF
dx se anulează . Luând o referință
u pentru care
derivata se anulează, fluxul numeric Godunov va lua forma [1]:
22
1/ 2 1 111max ,22i i iF u u
(4.41)
4.2.2 Discretizarea fluxuril or difuzive
Discretizarea fluxurilor difuzive se va face prin utilizarea unor scheme centrate în spațiu.
În continuare se va exemplifica pe cazul unei mișcări bidimensionale, pentru care ecuațiile
Navier -Stokes au forma:
U F F G G
t
(4.42)
Forma explicită a fluxurilor difuzive este:
1 3 1 3
3 2 4 2
40
Rea u a v c u c v
Ga u a v c u c v J
G
(4.43)
1 4 1 3
3 2 3 2
40
Rec u c v b u b v
Gc u c v b u b v J
G
(4.44)
unde:
2 2 2 2 2 2 2 2
1 2 3 1 2 34 4 1 4 4 1, , , , ,3 3 3 3 3 3x y x y x y x y x y x y a a a b b b
(4.45)
2 2 2 2
4 4 1 244, , , ,33x y x y x x y y x x y y a b c c
(4.46)
3 4 422,,33y x x y x y y x x x y y c c c
(4.47)
UNIVERSITATEA POLITEHNICA BUCURE STI
FACULTATEA DE INGINE RIE AEROSPATIALA
43 Scheme cu decuplarea vectorului flux
Pentru discretizarea sistemului Navier -Stokes ne vom referi la elementele necesare discretizării
sistemului.
GG
(4.48)
Vom considera termenul care pune în evidență existența unor termeni de forma
1au și
3av care
conțin derivatele după direcția
și a unor termeni de forma
1cu ,
3cv care conțin derivatele în lungul
direcției
:
1 3 1 32 ReG a u a v c u c vJ
(4.49)
Un astfel de termen îl vom nota cu
LM unde
1ReL a uJ
iar
Mu
Pentru derivata
G
, discretizând centrat, obținem:
1, 1,1, 1, 2
,1
2i j i ji j i j
ijGML L M L M
(4.50)
Similar , pentru derivata
G
vom obține o discretizare a fluxurilor vâscoase de forma:
1, , 1,21, , 1,
,,
1, , 1, 21, , 1,1
2
1
2i j i j i j
i j i j i j
i j i j
i j i j i ji j i j i jGGe G e G e G
f G f G f G
(4.51)
Figure 4.1. Compo nentele vectorului
G și
G din ecuația (4.15) [1]
UNIVERSITATEA POLITEHNICA BUCURE STI
FACULTATEA DE INGINE RIE AEROSPATIALA
44
5 Modele de turbulență
5.1 Descrierea statistică a turbulenței
Modelul matematic cel mai frecvent folosit în aplicațiile ing inerești pentru calculul mișcărilor
turbulent îl constituie ecuațiile mediate Reynolds.
Medierea ecuațiilor Navier -Stokes conduce la apariția unor termeni suplimentari care sunt
interpretați ca tensiuni aparente și fluxuri termice aparente asociate cu miș carea turbulentă. Acești
termeni sunt exprimați în funcție de parametrii medii prin intermediul modelelor de turbulență .
Scări de turbulență
Noțiunea de scară reprezintă un instrument foarte util în dinamica fluidelor deoarece , permite
evaluarea importn ței relative a diferitelor fenomene participante la mișcare determinând factorii secundari
ce pot fi neglijați ca și regiunile din spațiu unde anumite fenom ene pot fi luate în considerare [3]. Această
notiune se aplică atât lun gimilor caracteristice cât și timpilor caracteristici.
Timp de difuzie și timp de convecție
Timpul necesar pentru ca efectul peretelui să fie transmis la distanța
, se poate considera ca
timp caracteristic de difuzie :
2
dt
(5.1)
Timpul în care o particulă de fluid parcurge o distanță L, este timpul caracteristic fenomenului de
convecție :
cLtU
(5.2)
Conform ipotezei stratului limită scările de timp de difuzie și de convecție sunt de același ordin
de mărime:
1/2 21,.
ReL UL
UL
(5.3)
Într-o curgere turbulentă proceselor prezente în cazu l mainar (convecția și difuzia vâscoasă) li se
adaugă ca efect al fluctuțiilor vitezei fenomenul de difuzie turbulentă . Direcția dominantă a difuziei
turbulente o constituie direcția normalei la perete. Structurile turbulente de talie inferioară sunt gener ate
în mod continuu de către structurile turbulente de talie superioară, prin fenomenul de cascadă energetică.
Așadar este normală definirea unor scări distincte asociate fluctuațiilor turbulente de diferite amplitudini
[3].
UNIVERSITATEA POLITEHNICA BUCURE STI
FACULTATEA DE INGINE RIE AEROSPATIALA
45
Figure 5.1-Convecția și difuzia într -un strat limită turbulent [3]
Scara de timp care reprezintă timpul necesar unei particule de a se deplasa sub efectul fluctuțiilor
turnulent e pe direcția normalei pe o distanță
ˆl este:
ˆˆ
ˆkltu
(5.4)
Timpul caracteristic al mișcării de convecție și cel caracteristic al turbulenței sunt de același
ordin:
ˆ
ˆuU
Ll
(5.5)
Totodată această relație mai poate fi reprezentată astfel: scara de timp a turbulenței este impusă
de mișcarea de convecție cu viteza medie.
Scara de timp a turbulenței este impusă de viteza de forfecare:
1
ˆ ˆkUU
ty l
(5.6)
Difuzie moleculară și difuzie turbulentă
Dacă aerul este în mișcare datorită scăderii densității prin încălzire, timpul caracteristic. Dacă
există o mișcare turbulentă cu cea mai mare scară, L și viteza caracteristică,
ˆu , timpul caracteristic este:
ˆ
ˆtLtu
(5.7)
Scări laminare și scări turbulente
În comparație cu cazul laminar, fluctuațiile turbulente modifică în mod semnificativ toate scările.
Reprezentarea grosimii în cazul turbulent pe placa plană rezultă experimental (legea lui Blasius):
56
1/50.38, 5 10 Re 4 10Ret
x
xpentrux
(5.8)
0,30.076 Ret
x
l
(5.9)
Raportul dintre scara de timp a difuziei moleculare și dscara de timp a fluctuațiilor turbulente este
invers proporțional cu numărul Reynolds asociat:
2ˆˆˆ ˆˆRe ,ˆ ˆl
t
tt l u ul
t l
(5.10)
Scara de timp a difuziei turbulente este cu cel puțin două ordine de mărime mai mică decât scara
de timp a difuziei vâscoase.
UNIVERSITATEA POLITEHNICA BUCURE STI
FACULTATEA DE INGINE RIE AEROSPATIALA
46 Scări mici de turbulență
Fiind rezultatul un ei superpoziții de mișcări cu scări foarte diferite, turbulența eare o structură
foarte complexă. Structurile de amplitudini mari, se formează în interiorul structurilor turbulente. Acestea
se formează prin întindere care, la rândul lor formează mișcări la scări și mai mici. Datorită faptului că în
această zonă se găsește deficitul maxim de impuls și energie, acest proces se petrece simultan cu difuzia
structurilor turbulente succesive către perete. Limita inferioară a mărimilor caracteristice pentru
struct urile turbulente este legată de vâscozitatea fluidului. Pe măsură ce numărul Reynolds scade,
vâscozitatea crește și are rolul de a transforma energia cinetică în căldură. Pentru structurile mici,
disipația este foarte importantă și devine preponderentă [3].
Deși vâscozitatea are efecte neglijabile în zona turbulentă dezvoltată, ea, în domeniul scărilor
mici este în continuare prezentă și provoacă disiparea energiei cinetice a mișcării medii. Este de preferat
să considerăm că rata disipației trebuie să fie egală cu rata aportului de energie pentru susținerea
mișcării fluctuante staționare [3].
Este normală admiterea faptului că viteza de forfecare este proporțională cu rata disipației
energiei cinetice instantanee. Deoarece transformarea energiei cinetice a fluidului în căldură este un efect
al vâscozității moleculare, este evident ca rata dispației tr ebuie să depindă de vâscozitate [3].
d ij ij ij ijR C s s C s s
(5.11)
Această relație pune în evidenșă faptul că media ratei de disipație a unei mișcări turbulente apare
ca suma a termenilor:
ij ij C s s care reprezintă rata de disipație asociată mișcă rii medii și
ij ij C s s
reprezintă rata de disipație asociată câmpului de viteze fluctuante.
Parametrii cei mai importanți pentru a defeni cele mai mici scări sunt: rata disipației pe unitatea
de masă și vâscozitatea cinematică. Cu aceșt i parametrii se pot construi scările de timp, viteză și lungime
denumite microscări Kolmogorov.
1/ 4 1/ 2 3
1/ 4ˆ ˆ ,,tu
(5.12)
Este utilă corelarea ca ordin de mărime a scărilor cu scările legate de scări le mari de turbulență:
3/4 1/2 1/4ˆˆˆ 1 1 1,,ˆˆ ˆ Re Re Ret t tt u u
u ll
(5.13)
ˆˆRetul
(5.14)
Aceste relații arată că discrepanța dintre scările pentru cele mai mici nuclee și scările pentru
nucleele mari, crește pe măsură ce numărul Reynolds al turbulenței crește.
5.2 Modelul
k standard
Modelul
k este cel mai popular model de turbulență. Referința de bază pentru acest model
rămane lucrarea lui Jones și Launder. Acum, este cunoscut sub numere de modelul
k standard. Acesta
a fost dedus pentru mișcări la numere Reynolds mari. Impunerea condițiilor la limită pentru disipația
turbulentă, pentru zona din vecinătatea peretelui se face utilizând legea la perete coroborată cu ipoteza
echilibrului dintre disipație și energia cinetică turbulentă. Acest lucru nu permite extinderea calculului
până la suprafață așadar, se pierd o serie de detalii legate de: va riația energiei turbulente, variația vitezei și
a disipației turbulente.
Ecuația energiei cinetice turbulente, poate fi obținută din ecuațiile Navier -Stokes prin
descompunerea Reynolds:
211
2i i i k k i i
k i k
k k k k k k k kK K U U U U U P K U UU U Ut x x x x x x x x
(5.15)
UNIVERSITATEA POLITEHNICA BUCURE STI
FACULTATEA DE INGINE RIE AEROSPATIALA
47
Exceptând ultimii doi termeni, toți ceilalți termeni trebuiesc modelați.
Modelarea termenului de producție
Dacă admitem ipoteza lui Boussinesq atunci termenul de producție poate fi modelat prin
înlocuirea tensiunii tur bulente cu valoarea dată de ipoteza Boussinesq.
2
3j ii
i k t ij i k
k k kU UUU U K U Ux x x
(5.16)
2
3j i i i
i j t ij
j i j jU U U UU U Kx x x x
(5.17)
unde:
t
tq
(5.18)
reprezintă vâscozitatea de vârtej.
Termenul de transport din punct de vedere fizic este un termen în care transportul se realizează ca
o difuzie sub efectul mișcării fluctuante, așa cum difuzia moleculară se reali zează sub efectul agitației
moleculare.
1
2j i i k
kjKU UUU
xx
(5.19)
Pentru a modela acest termen se introduce o ipoteză similară difuziei moleculare care este cunoscută sub
numele de ipoteza gradientului .
t
j
kjKKUx
(5.20)
1
2i i k t
k j k jU U U K
x x x
(5.21)
Termenul de interacțiune viteză -presiune
Din punct de vedere fizic, este un termen de difuzie turbulentă produsă de interacțiunea dintre
fluctuațiile vitezei și cele ale presiunii. Acest termen este modelat tot cu ipoteza gradientului. Variația
presiunii, va amplifica sau va diminua difuzia turbulentă. Dacă presiunea crește în lungul curger ii difuzia
scade iar dacă presiunea scade în lungul curgerii atunci difuzia este amplificată.
Ecuația modelată a disipației este:
2
12j i i t
jt
j j i j j jU UUU C Ct x x x x K K x x
(5.22)
În general, difuzia vâscoasă pentru turbu lența dezvoltată este neglijabilă pentru că începând cu
zona logaritmică difuzia turbulentă devine dominantă.
Termenul de producție pentru
este scris prin raportare la energia cinetică. Producția pozitivă și
negativă fiin alter ate cu două constante empirice. Rata de variație a
energiei cinetice
turbulente Termen
de
producție Termen
de
transport Interacțiu
ne viteză –
presiune Termen
de difuzie
vâscoasă Disipația
turbulentă
UNIVERSITATEA POLITEHNICA BUCURE STI
FACULTATEA DE INGINE RIE AEROSPATIALA
48 Raportul
k
care apare în expresia termenului de producție are semnificație unei scări de timp.
Idee exploatată într -o serie de variante ale modelului numite modele de tip număr Reynolds mic .
121,44, 1,92, 1,3 CC
(5.23)
Vâscozitatea turbulentă se calculează cu ipoteza Bradshaw :
2
0,09tKC
C
(5.24)
Constantele de închidere
În modelu l
k standard constatele de empirice care apar, trebuie să asigure proprietatea
modelului de a fi un model complet (model care se aplică fără a fi nevoie de a schimba mărimi). În
deducerea valorilor acestor constante, Lauder și Jone s au apelat la rezultate experimentale valabile pentru
o serie de probleme test de curgeri turbulente.
Pentru deducerea
C s-au făcut următoarele ipoteze:
– curgerea turbulentă este complet dezvoltată;
– pentru turbulența dezvoltată bilanțul energetic pune în evidență faptul că producția de
turbulență este echilibrată de disipația turbulentă, transportul este neglijabil;
Celelalte constante,
12,CC din ecuația disipației sunt obținute din condiția ca modelul să
reproducă cazul turbulenței de grilă.
Dacă avem un grilaj plasat într -un curent atunci la curgerea în jurul fiecărei bare se formează un
siaj turbulent.
Dacă avem o grilă amplasată într -o suflerie odată cu curgerea spre aval, turbulența scade, chiar
dispare și se numește rata de atenuare a turbulenței.
Experimental, în avalul grilei, energia cinetică turbulentă
mKx ,
1,25 0,06m . În
lucrările care abordează turbulența de grilă se specifică două intervale de variație a lui m. Valori mai mici
se găsesc în zona de început a grilei numită zona inițială și valori ceva mai mari la depărtare față de grilă
numită faza finală a turbulenței de grilă.
Pentru aplicarea modelului la turbulența de grilă se mențin ipotezele anter ioare privind turbulența
dezvoltată dar nu se neglijează difuzia turbulentă din ecuația lui
. Din acest caz rezultă o relație de
legătură între constantele
1C și
2C .
2
12 1/2xCCC
(5.25)
Adoptând o valoare medie pentru
1,25m se obține
21,92 C .
Dacă
0,41x iar
1,3 atunci
11,44 C .
Singurele constante care au fost alese, având ca referință lucrări publicate anteriro de către
Kolmogorov și Nokoyama sunt
1k și
1,3 .
Condiții la limită
Una dintre deficiențele majore ale modelului standard, este legată de faptul că ecuația disipației
turbulente
nu este valabilă la perete și în zona imediat adiacentă a peretelui (substrat vâscos).
Din acest motiv trebuie să impunem condiții la limită pentru mod elul
k standard nu pe
suprafața peretelui ci într -un punct din curgere care este suficient de îndepărtat de perete pentru a putea fi
considerat în zona inerțială.
UNIVERSITATEA POLITEHNICA BUCURE STI
FACULTATEA DE INGINE RIE AEROSPATIALA
49
În punctul P, făcând ipoteza turbulenței în echilibru, putem ded uce valoarea energiei cinetice
turbulente.
1/2 2
pKCU
(5.26)
3
, 30 40P K P
PUPyxy
(5.27)
Viteza de forfecare din condițiile la limită ale modelului
k standard se obține din punerea
legii la perete:
30 40
11ln lnP
PP
Py
U U yy C CU
(5.28)
Această ecuație devine o ecuație transcendentă din care se calculează viteza de frecare.
Dacă peretele nu este plan se procedează absolut similar doar că în toate relațiile vitezei
PU este
înlocuită cu viteza în direcție tangentei locale la perete.
Avantaje ale folosirii modelului:
modelul este relative simplu de implementat ;
conduce la calcule stabile care converg relativ ușor ;
Dezavantaje ale folosirii modelului:
predicții deficitare pentru:
– curgerile turbionare și rotative;
– curger ile cu o separare puternică;
– jeturile axial simetrice;
– curgerile pe depli n dezvoltate în canalele ne -circulare
valabil numai pentru curgerile complet turbulente ;
ecuația ε este simplistă .
6 Program ul NSC2KE
Software -ul dinamicii fluidelor NSC2KE folosește o combinație de metode cu volum e finite și
metode cu elemente finite pent ru a simula o gamă largă de câmpuri de curgere care merg de la regim
subsonic la regim hipersonic pentru config urații bidimensionale și axialsi metrice. Deoarece configurația
UNIVERSITATEA POLITEHNICA BUCURE STI
FACULTATEA DE INGINE RIE AEROSPATIALA
50 reală implică adesea geometrii complexe, gri lele nestructurate par a fi cel mai bu n instrument pentru a
aborda corect realitatea.
Pentru a rezolva parte a Euler a ecuațiilor, sunt disponibili trei rezolvitori: un Roe, un Osher și un
Cinetic. Pentru a calcula curgerile turbulente un model
k este disponibil . Turbu lența din apropierea
peretelui este calculată utilizând legea pereților sau abordearea în două straturi. Problemele dependente de
timp sunt rezolvate cu ajutorul solverului Runge -Kutta de ordinul 4.
Codul este aplicat într -o forma adimensională. Pentru a d efini configurația calculului, utilizatorul
va trebui să definească două variabile adimensionale, mai exact, numarul Mach și numărul Reynolds, care
vor descrie câmpul curgerii.
NSC2KE simulează o gamă largă de curgeri, inclusiv:
– curgeri interne și extern e;
– curgeri nevâscoase de la subsonic la hipersonic ;
– curgeri nevâscoase de la subsonic la hipersonic la numere Reynolds mici ;
– curgeri staționare și nestaționare;
6.1 Ecuațiile guvernatoare
densitatea,
12, u u u viteza, T temperatura,
2
2uET energia totală,
1 pT
presiunea,
,ij uu gradientul lui u,
,ij Du este divergența,
2
3tS u u DI tensorul
deformației.
Ecuațiile Na vier-Stokes
0 ut
(6.1)
tuu u p St
(6.2)
ttEE p u Su k k Tt
(6.3)
Mod elul
k la numere Reynolds mari
,tkkuk k St
(6.4)
.t u c St
(6.5)
unde termenul de producție este:
2
3ktS P kD
(6.6)
iar termenul de disipație este:
2
1
122
3cS c kP D cck
(6.7)
UNIVERSITATEA POLITEHNICA BUCURE STI
FACULTATEA DE INGINE RIE AEROSPATIALA
51 Formularea globală
Vom considera următoarea formă conserva tivă a ecuațiilor Navier -Stokes și a modelului de
turbulență fără termeni sursă:
0WF W N Wt
(6.8)
unde W este vectorul variabilelor conservative iar N reprezintă operatorii vâscoși și cei advec tivi.
Condiții la limită
Condițiile la limită disponibile sunt destul de simple. Datorită faptului că s -a încercat să nu se
facă un cod prea complicat, este disponibil un set complet, dar nu exhaustiv, de condiții la limită care
permit utilizatorului s ă calculeze o gamă largă de domenii de curgere , inclusiv curgeri interne și externe
[5].
Corpuri solide
Pentru calcule laminare condiția la limită clasică pentru viteză es te condiția fără alunecare
0u
în
timp ce pentru temeperatură poate fi f olosită o condiție la limită Neu mann sau Dirichlet, depinde de
problemă.
Pentru calcule turbulente utilizatorul poate folosi abordarea cu legile la perete sau abo rdarea în două
straturi (two -layer). Pentru primul caz este folosită o condiție la limită slabă pentru viteză, o condiție
Newmann pentru temperatură iar pentru k și
se folosește o condiție la limită neomogenă Dirichlet .
Pentru cel de -al doilea caz sunt aplicate aceleași condiții pe ntru u și T, pentru k este folosită o condiție
omogenă Dirichlet și Neumann pentru
.
Limitele slabe ale curgerii interioare și exterioare
NSC2KE tratează limitele infinite printr -o tehnică caracteristică. Acest lucru înseamnă că fluxurile
sunt împărțite în parți pozitive și parți negative , urmând semnul valorilor proprii ale matricei iacobiene;
Profilul de admisie
În calculul curgerilor trebuiesc formultate condiții in secțiunea de intrare. Aceste date pot fi obținute
fie din experimente, fie dintr -o altă simulare numerică.
Simetria și condițiile de alunecare
Condițiile de nepenetrare (impermeabilitate) (
0 un
) sunt necesare pentru calculele Euler și uneori
în cazurile vâscoase. Această condiție este aplic ată pe limite solide pentru calculele Euler și pe liniile
simetrice, în funcție de fizica problemei.
6.2 Tehnici numerice
– Algoritmul este explicit în timp;
– Pentru a ajunge la soluția stabilă este utilizată o schemă iterativă;
– Ecuațiile Navier -Stokes și k -epsilon sunt rezolvate printr -o tehnică Upwind cu Volume FInite –
Galerkin, folosind un solver Osher, Roe sau Cinetic pentru partea convectivă a ecuațiilor;
– Termenii vâscoși sunt tratați folosind o tehnică standard Galerkin.
UNIVERSITATEA POLITEHNICA BUCURE STI
FACULTATEA DE INGINE RIE AEROSPATIALA
52 6.2.1 Rezolvitorul Navier -Stokes
Fie
h j jU discretizarea în triunghiuri a domeniului Ω și
h i iUC împarțirea în celule.
Figure 6.1 Discretizarea domeniului
Putem asocia fiecărui
hhwV , unde
hV este un set de funcții continue
1P în triangulația noastră, și
/
hw
o funcție constantă pe o porțiune pe celulă:
/ 1
ih i hC
iw C wC
(6.1)
În schimb, știind că
/
hw este constanta unei porțiuni,
hw este obținut prin:
/()h i h iw S w C
(6.2)
Dacă presupunem că F variază liniar pe fiecare triunghi, formularea slabă a ecuației globale (3.19),
luând în calcul
hV este:
Găsim
6()hhWV astfel încât,
hhV
( )( ) ( ) ( ) 0h
h h h h h h h hWF N W F N nt
(6.3)
este echivalentă cu urmatoarea formulare slabă, obținută prin scoaterea părții convective a ecuației
anterioare pentru
h fiind funcția caracteristică a lui
iC și folosind un timp explicit de integrare:
1
()
inn
n
idCWWC F W n RHSt
(6.4)
Datorită faptului că folosim o schemă centrată , putem calcula partea din dreapta a problemei:
. . . ( ) ( ) ( )
hnn
nh R H S N W N W n
(6.5)
Știind că,
( ) ( )n
dhF W F W pentru
iC , și în rest
dF este o porțiune constantă aproximată
a lui F(W) verificând:
UNIVERSITATEA POLITEHNICA BUCURE STI
FACULTATEA DE INGINE RIE AEROSPATIALA
53
1( ,v) ( (u) F(v)) d(u,v)2uF
(6.6)
unde d este difuzia numerică. În NSC2KE utilizatorul poate alege diferite proporții pentru d în funcție d e
schema aleasă (Roe, Osher și C inetic).
Formularea precedentă este corectă doar pentru ecuații de gradul I. Pen tru cele de gradul II
precizia se obține folosind o extensie de tip MUSCL (Schemă Monotonă Amonte -Centrată pentru Legile
de conservare ), implicând combinații de scheme upwind și gradienți centrați. Mai precis, definim
iW ca
fiind o aproximare a gradientului W la nodul i. Definim pe segmentul
,ij cantitățile următoare:
0.5 ( ( ) ,(1 )( ))
0.5 ( ( ) ,(1 )( ))ij i i i j
ji j j j iW W Lim W ij W W
W W Lim W ij W W
(6.7)
cu limitator de tip Van Albada:
22
22(a ) ( )( ,b) 0.5(1 sgn(ab)) 1 02b b aLim a cuab
(6.8)
conține cantitatea de urcare de la un punct la altul. C el de -al doilea ordin de acuratețe în spațiu este
obținut prin inlocuirea lui
/
iW și
/
jW în (6.7 ) prin
ijW și
jiW .
Totuși , această abordare nu garantează pozitivitatea termenilor
k și
, prin urmare fluxu rile
convective pentru ec uațiile turbulente sunt calculate folosind schema propusă pentru conservarea
pozitivității speciilor chimice.
Mai precis, odată ce densitatea este calculată, curgerile turbulente sunt deduse cu:
0( . 0) . ( . ) . .
i i i i i iijC C C C C Cuk n k resp k u n dacã u resp n
(6.9)
Termenii sursă pentru modelul k – ε au fost luați în vedere în mod explicit.
6.2.2 Procedura de integrare în timp
WRHS Wt
(6.9)
În această ecuație integrarea în timp s-a făcut cu o schemă explicită. NSC2KE folosește mai degrabă
următoare schemă Runge -Kutta de ordinul patru:
0 nWW
(6.10)
01, 1,…,4kk
k W W tRHS W k
(6.11)
14nWW
(6.12)
unde, alegerea optimă pentru
k este:
1 2 3 40,11, 0,2766, 0,5, 1,0.
(6.13)
Această schemă ne permite deasemenea să calculăm curgerile dependente de timp.
UNIVERSITATEA POLITEHNICA BUCURE STI
FACULTATEA DE INGINE RIE AEROSPATIALA
54
Pasul de timp local . NSC2KE folosește următoarea formulă pentru a calcula pasul de timp local
pentru un nod dat,
i .
2Prmin ,2i
txxtuc
(6.14)
unde,
x este înălțimea minimă a triunghiurilor având nodul
i în comun.
Pentru calculul Navier -Stokes merge deasemenea și pasul de timp Euler, înseamnând că putem
lua doa r primul termen din ecuația mai sus menționată. Acest lucru nu este întotdeauna adevărat așa că,
trebuie să fim atenți când alegem stategie pentru pasul de timp. Așadar, pentru cazul stașionar, trebuie să
folosim mereu o strategie pentru pasul de timp loca l în schimb, pentru cazul nestaționar este necesar. un
pas de timp global.
Corecții axialsimetrice
Pentru a adapta rezolvitorul la calcul axialsimetric, multiplicăm zonele celulare
iC , zonele
triunghiulare
iT și lung imile de frotieră
iE cu raza r .
,i i iC C r
(6.15)
unde
ir este a doua coordonată a nodului
i .
1 2 3,3iir r rTT
(6.16)
1 2 3,,r r r
sunt razele corespunzătoare nodurilor triunghiului.
12,2iirrEE
(6.17)
iar,
12,rr sunt razele extremităților
iE .
6.3 Alegerea schemei numerice
Utilizatorul are posibilitatea de a defini tipul de solver pe care îl folosește. Solverii Roe, Osher și
Cinetic sunt disponibili pentru partea Euler a ecuațiilor , în timp ce pentru partea vâscoasă este disponibilă
doar o schemă centrată . Cei trei solvers pot fi utilizați pentru regimul subsonic sau transonic. Utilizatorul
poate alege solverul dorit prin parametrul iflux din fișierul DATA.
1;
2;
3;iflux rezolvitor Roe
iflux rezolvitor Osher
iflux rezolvitor Cinetic
(6.18)
Precizia schemei . NSC2KE propune scheme cu ordinul întâi și doi de precizie. Putem alege
ordinul de precizie al schemei. Acest lucru se face prin parametru nordre .
1;
;
2 lim ;nordre ordinul întâi
nordre ordinul doi
nordre ordinul doi itat
(6.19)
Regula este aceea de a începe cu scheme cu ordinul întâi și de a schimba la ordinul doi după un
timp.
UNIVERSITATEA POLITEHNICA BUCURE STI
FACULTATEA DE INGINE RIE AEROSPATIALA
55 Modelul de turbulență. Utilizatorul poate alege abordarea clasică a legii pereților sau abordarea
în două straturi. Cea mai bună alegere este abordarea în două straturi în acest caz deoarece descrie
curgerea în apropierea peretelui aproape fara empirism. Totuși rezolvarea cu acest model va dura mai
mult deoarece utilizatorul trebuie sa folosească o grilă m ult mai detaliată.
Alegerea tipului de abordare se va face alegând numărul aferent abordării dorite.
1 –> =0 two -layer technique, =1 wall laws
Rețeaua de calcul. Într-un calcul, primul punct constă în definirea domeniului de calcul. Acest
lucru este făcut printr -un fișier ce conține punctele grilei și ordinea acestora. Este nevoie și de un tabel de
conectivitate, ce descrie cum sunt aceste puncte conectate între ele. Programul „NSC2KE” citește fișierul
rețelei folosind următorul format al acest eia:
Figure 6.2. Citirea fișierului rețelei
Condiții inițiale . Pentru a inițializa calculele, programul NSC2KE propune două strategii.
Utilizatorul poate începe fie de la o stare inițială uniformă defini tă de condițiile curgerii ( numărul Mach,
respectiv temperatura din interiorul curgerii din fișierul DATA) sau să reînceapă de la o soluție obținută
anterior definită în fișierul INIT_NS.
Pentru a folosi modelul de turbulență k – ε, utilizatorul poate de asemenea să înceapă de la o stare
inițială uniformă sau să folosească o soluție anterioară definită în fișierul INIT_KE.
6.4 Funcțiile subrutinelor
Tabel 6.1 Schema logică a programului NSC2KE
Programul conț ine 27 de surse dintre care un fișier de parametrii și 26 de subrutine. Fișierul
param2d conține parametrii și definiții comune pentru rezolverul 2D Navier -Stokes. Subrutina nsc2ke
constituie subrutina principală în care se apelează următoarele subrutine c are au următoarele funcții:
subrutina config în care sunt citiți parametrii dependenți ai aplicației . Se citesc parametri ca: parametrii
UNIVERSITATEA POLITEHNICA BUCURE STI
FACULTATEA DE INGINE RIE AEROSPATIALA
56 gazului perfect, numărul Pandtl pentru regim laminar și turbulent, numărul Mach, unghiul de incidență,
parametrii proce sului de integrare Runge -Kutta, parametrii aferenți definirii tipului de rezolver (Roe,
Osher, Cinetic). Încărcarea structurii grilei globale se face prin apelarea subrutinelor mailla și geowall . În
subrutina mailla se definesc variabilele locale și se fac e strategia de renumerotare a nodurilor. Calculul
volumelor de control și a zonelor triunghiulare se face prin apelarea subrutinei aires . Așadar se
cazlculează zona triunghiulară curentă iar aceasta se adună în volumul de control. În cadrul subrutinei
clhaut se calculează altitudinile volumului de control. Acest calcul rezultă din definirea varibilelor locale
și inițializarea altitudinilor de vârf. Construcția frontierelor triunghiulării se face prin apelarea
subrutinelor axigeo și seg2d . În cadrul acesteia se definesc variabilele locale, se calculează limitele
normale ale volumului de control, se normalizează limitele normale și se construiesc marginile limită.
Prin apelarea subrutinei init se inițializează soluția fizică după care se schimbă. Calculul pres iunilor
inițiale se realizează prin apelarea subrutinei calprc . Prin apelarea subrutinei caldtl se calculează
vâscozitaea sutherland. În continuare, se apelează subrutina result prin care se salvează în fișiere
rezultatul. La rândul ei, această subrutină a pelează subrutina isovat.
Ultima subrutină apelată în subrutina nsc2ke este resexp. În cadrul acesteia se calculează procesul
de integrare în timp Runge -Kutta . Pentru rezoluția explicită, subrutina resexp apelează subrutinele fluroe,
fluosh, flucin, keps2 d, gravity, vitfrot și sa iar pentru tratamentul condițiilor la limită se apelează subrutina
cdl. Astfel, în subrutina fluroe se definesc varibilele locale după care : se indexează local nodurile marginii
curente , se face adresarea indirectă pe noduri a stă rilor fizice, se calculează entalpia, se calculează valorile
medii Roe, valorile proprii, partea centrată a fluxului, partea difuzivă a fluxului iar în final se face
colectarea fluxurilor elementare în cele globale. În subrutinele fluosh și flucin se face aceleași calcule ca
în subrutina fluroe cu diferențe la schimbarea rezolverului. Subrutina keps2d face calculele aferente
modelului de turbulență. Aceasta la rândul ei, apelează subrutinele ke_law și ke_two în care se fac calulele
pentru cazul legii la per ete și cel al stratului dublu. Subrutina gravity definește forțele gravitaționale. În
ultima subrutină apelată vitfrot , prin apelarea subrutinei loglaw
7 Rezultate numerice
8 Concluzii
9 Bibliografie
1. Dănăilă, S., Berbente, C., “Metode numerice în dinamica flu idelor”, București, Editura
Academiei Romane, 2003 .
2. Carafoli E., Constantinescu V.N., “Dinamica fluidelor compresibile ”, București, Editura
Academiei Republicii socialiste Romania, 1984.
3. Constantinescu V.N., Dănăilă, S., Găletușe S,, “D inamica fluidelor în regim turbulent ”,
București, Edit ura Academiei Romane, 2008.
4. www.wikiepdia.org
5. INRIA, Bijan Mohammadi, “Fluid Dynamics Computation with NSC2KE , An User -Guide ,
Release 1.0 ”, Franța, 1994.
6. www.creeaza.com/tehnologie/aeronautica
7. J Délery, J -P Dussauge , “Some physical aspects of shock wave/boundary layer
interactions ”, 26 July 2009 .
Copyright Notice
© Licențiada.org respectă drepturile de proprietate intelectuală și așteaptă ca toți utilizatorii să facă același lucru. Dacă consideri că un conținut de pe site încalcă drepturile tale de autor, te rugăm să trimiți o notificare DMCA.
Acest articol: Modelarea matematică si numerică a curgerilor [617044] (ID: 617044)
Dacă considerați că acest conținut vă încalcă drepturile de autor, vă rugăm să depuneți o cerere pe pagina noastră Copyright Takedown.