In partea a patra sunt prezentate concluziile acestei lucrări, contribuțiile personale si direcțiile de cercetare ulterioare in acest domeniu. [301414]
Introducere
Dezvoltarea rapida a zonelor urbane a [anonimizat], a condus la suprasolicitarea rețelelor existente de canalizare. [anonimizat]. [anonimizat] a apei se realizează ținând cont de depunerile de sedimente in decursul timpului. [anonimizat], estimat pentru diverse materiale de construcție a canalului. O estimare cat mai precisa a coeficientului de rugozitate din formula lui Manning pentru un colector orășenesc este necesara pentru a evita incidentele in exploatarea rețelelor de canalizare urbana. Acest lucru se poate determina pe cale experimentala (măsurări PIV) sau cu ajutorul modelarii numerice (CFD).
Prezenta lucrare este compusa din patru parți. In prima parte, s-a realizat o cercetare documentara privind situația actuala si cercetările realizate pana in prezent pe modele asemănătoare cu cel propus de noi.
In partea a doua a lucrării s-a [anonimizat] a Instalațiilor din București. Pe acest model fizic s-a realizat o [anonimizat] s-au determinat o serie de parametri atât pentru curgerea peste o frontiera lisa cat si pentru curgerea peste o frontiera rugoasa. Rezultatele s-au comparat cu cele prezentate in literatura.
In urma prelucrării datelor experimentale s-au calculat profilele de viteză și mărimile turbulente ale curgerii, s-au comparat diverse metode pentru calculul vitezei de frecare, s-a realizat unui studiu de sensibilitate privind estimarea vitezei de frecare pornind de la distribuția logaritmică a profilului de viteze si s-a estimat coeficientului de rugozitate într-o conductă în care au apărut deja depuneri de sedimente și frontiera curgerii nu mai este netedă.
In partea a treia a lucrării, s-a realizat un model numeric pornind de la caracteristicile modelului fizic de laborator. Astfel se prezintă elaborarea unor modele numerice al curgerii studiate în partea a doua, utilizând o largă paletă de posibilități oferite de tehnica CFD (Computational Fluid Dynamics) și de codul comercial Ansys Fluent. Rugozitățile artificiale studiate pe cale experimentală au fost reproduse în modelul geometric utilizat pentru simulările numerice.
[anonimizat].
In partea a [anonimizat].
Studiu bibliografic
Viteza de frecare
Viteza de frecare, u*, [anonimizat]. [anonimizat], un factor esential in aplicatii fizice sau legate de mediul inconjurator. Acesta este o variabila fundamentala si un parametru de scalare a [anonimizat]. Estimarea pragurilor critice ale eroziunii sau ale depunerii sedimentelor necesita determinarea fortelor hidrodinamice aplicate stratului de sedimente. Estimarea pragului critic al transportului de sedimente este strans legat de determinarea efortului de forfecare.
Efortul de forfecare τ este strans legat de viteza de frecare si se determina cu formula:
(1)
ρ – densitatea apei
Viteza de frecare este strans legata de curgerea turbulenta in apropierea frontierei si este esentiala in intelegerea dezvoltarii turbulentei in apropierea peretelui. In cazul curgerilor peste un pat de pietris, efortul de forfecare nu se poate masura direct si este determinat pe cale indirecta prin estimarea lui u*. Datorita faptului ca efortul de forfecare variaza direct proportional cu patratul vitezei de frecare, este necesar ca determinarea valorilor lui u* sa se face cu o precizie cat mai ridicata.
Viteza de frecare se determina utilizand profile de viteze, in special pentru cele masurate in stratul interior. In cazul curgerilor peste un strat de pietris liber, se dezvolta un strat rugos in imediata vecinatate a patului de sedimente (Rauph, 1981; Nikora si Smart, 1997), afectand profilul de viteze pentru stratul inferior si cel interior. Dinamica curgerii in stratul rugos este direct influentata de scara ce se utilizeaza pentru a defini rugozitatea. Aceasta curgere este, de regula, tridimensionala. Nu exista o marime unic definita pentru a caracteriza dimensiunile rugozitatii. In literatura, se regasesc, in schimb, mai multe marimi pentru a carcateriza rugozitatea:D50 – Townsend (1976); doua pana la cinci diametre D, Raupach et al. (1991); trei diametre D, Wilcock (1996); σD – abataerea medie patratica a inaltimii rugazitatilor, Nikora si Goring (2000); 0,05 h (h – adancimea apei), Smart (1999). Deasupra stratului rugos, curgerea se poate considera peste frontiera neteda. Stabilirea tipului de curgere depinde, de asemenea, de rugozitatea relativa D/h. Daca rugozitatea relativa creste, inaltimea de curgere poate sa nu fie suficienta pentru ca profilul de viteze sa respecte legile curgerii peste o frontiera neteda. Katul et al. (2002) sugereaza ca aceste legi sunt valabile pana la limita h<10D.
Viteza de frecare in curgerile peste un strat rugos se determina utilizand profile de viteze detaliate pe intreaga inaltime de curgere, pentru ca aceasta sa nu fie influentata de rugozitatea a carei inaltime nu este cunoscuta. In cazul raurilor cu pat sedimentar, Wilcock (1996) si Smart (1999) au utilizat metoda profilului logaritmic pentru a determina viteza de frecare, Nikora si Goring (2000) au aplicat metoda efortului Reynolds, iar MacVicar si Roy (2007) au aplicat metoda energiei cinetice turbulente (TKE). Rowinski et al. (2005) au comparat aceste metode pe un model experimental cu pat de pietris ancorat si au ajuns la concluzia ca metoda profilului logaritmic nu se preteaza in acest caz. S-a subliniat faptul ca, datorita neregularitatii stratului de pietris brut, este dificil sa se foloseasca masurari intr-un singur punctpentru a determina viteza de frecare cu ajutorul metodelor utilizate.
In anumite conditii, poate fi dificil de obtinut un profil real de viteze, in specia in campul curgerii. In acest caz, viteza de frecare trebuie estimata pornind de la musrari intr-un singur punct. Aceasta abordare este des utilizata in cazul studiilor oceanografice (Stapelon si Huntley, 1995). Kim et al. (2005) au utilizat aceleasi metode ca si Rowinski et al. (2005), cat si metoda spectrala in cazul curgerii intr-un estuar peste un strat de mal. Biron et al. (2004) au realizat o serie de experimente de laborator privind curgerea peste un strat de nisip fin utilizand aceleasi metode si au ajuns la concluzia ca metoda profilului logaritmic pentru determinarea vitezei de frecare a condus la obtinerea unor valori mai mari decat in cazul celorlalte metode. Ei sugereaza ca masurarile trebuie realizate intr-un punct situat la distanta 0,1h, dar deasupra stratului rugos. Pentru validarea acestui punct, studiul trebuie continuat pentru mai multe tipuri de curgere.
Majoritatea studiilor au fost efectuate pentru straturi de sedimente fine (D50 <0,5cm). Bagherimiyab si Lemmin (2013) propun un studiu de laborator pentru curgerea peste un strat rugos a carei dimensiune caracteristica este D50 = 1,5cm.
Tehnici de estimare a vitezei de frecare
Viteza de frecare a fost initial definita cu ajutorul teoriei curgerii peste stratul limita. Metodele utilizate de obicei se bazeaza pe ipoteza existentei un efort de forfecare constant pe toata inaltimea curgerii, ceea ce nu se respecta in curgerea intr-un canal deschis. Cu toate acestea, teoria se aplica pentru sudiile de laborator privind curgerea cu suprafata libera. Pentru curgerea cu suprafata libera peste o frontiera rugoasa, Nezu si Nakagawa (1993) indica patru metode pentru calculul vitezei de frecare si a efortului de forfecare: (1) pornind de la inclinarea patului de sedimente in conditiile unei curgeri uniforme, (2) metoda profilului logaritmic, (3) metoda efortului Reynolds, (4) masurari directe.
In cazul curgerilor uniforme bi-dimensionale, viteza de frecare se poate estima pe baza echilibrului fortelor si poate fi data ca valoare de referinta:
()
g – acceleratia gravitationala
R – raza hidraulica
I – panta stratului sedimentar
Totusi, aceasta metoda ofera o valoare general estimata si nu este adecvata pentru a caracteriza curgerea. In cazul curgerilor rugoase, valorile masurate ale vitezei de frecare difera semnificativ fata de cele calculate cu ajutorul formulei de mai sus, in special datorita variatiei stratului rugos.
Metoda profilului logaritmic de viteze
Katul et al. (2002) sugereaza faptul ca poate exista un profil logaritmic de viteze in cazul curgerii peste suprafete rugoase, acesta manifestandu-se in partea inferioara a curgerii (mai exact in cincimea inferioara a curgerii), pentru cazul in care h > 10D. Distributia logaritmica a vitezei este descrisa cu ajutorul ecuatiei von Karman-Prandtl equation (Schlichting, 1987):
(3)
k – contanta von Karman, cu valoarea utilizata k = 0,41
Dependenta lui k fata de tipul curgerii este prezentata de Ferreira et al. (2002). z0 reprezinta lungimea hidraulica catacteristica a rugozitatii, iar u reprezinta viteza medie pe directia curgerii, la inaltimea z deasupra startului rugos. Monin si Yaglom (1971) definesc z0 ca fiind inaltimea la care viteza medie de curgere devine zero, daca regula logaritmica s-ar putea aplica de la acesta inaltime in jos. In cazul curegrii peste o frontiera rugoasa, atat marimea relativa a lui z0 cat si lungimea scalata a rugozitatii sunt importante pentru determinarea limitelor de aplicabilitatea a legii logaritmice. Pentru un strat rugos format din nisip omogen, Monin si Yaglom (1971) au stabilit relatia z0/D50 =1/30. In cazul rugozitatii cu dimensiuni neregulate, se poate aplica regula z0/D50 =1/10 sau chiar z0/D50 =1/5 (Monin si Yaglom, 1971; Townsend, 1976). Acest coeficient poate, insa, sa varieze, tinand cont de faptul ca el depinde de forma rugozitatii.
Metoda profilului logaritmic de viteze este adesea utilizata in studiul curgerii cu suprafata libera (in canale deschise si in rauri) (Nezu si Nakagawa, 1993). Aceasta prezinta avantajul ca valoarea estimata a lui z0 nu intervine in calculul lui u* , intrucat aceasta depinde doar de panta profilului. In cazul curgerii peste frontiere rugoase, profilul logaritmic de viteze se dezvolta deasupra acestui strat.
Metoda efortului Reynolds
Cand se pot efectua masurari ale turbulentei, viteza medie locala de frecare se poate determina cu ajutorul valorilor distributiei efortului Reynolds:
(4)
in care u’ si v’ reprezinta fluctuatiile vitezei pe directia curgerii, respectiv pe directia perpendiculara. Ele sunt mediate in functie de timp. Datorita efortului interior volumului in care se realizeaza masurarile, instrumentele de masura cu efect Doppler au o acuratete mai scazuta in straturile in care gradientul de viteza este mai accentuat, cum este cazul zonei din apropierea patului sedimentar (Lhermitte si Lemmin, 1994; Dombroski si Crimaldi, 2007). Aceasta metoda este, de asemenea, sensibila la diviatiile 2D in cazul curgerilor neuniforme (Nezu si Nakagawa, 1993; Kim et al., 2000; Nikora si Goring, 2000; Albayrak si Lemmin, 2011) si este necesara o aliniere stricta a sondelor. In cazul curgerilor cu suprafata libera peste frontiere rugoase, efortul Reynolds variaza liniar in stratul exterior (Nezu si Nakagawa, 1993). Acesti autori propun o extrapolare a efortului Reynolds peste frontiera rugoasa:
(5)
In plus, aceasta metoda permite verificarea cinditiilor de curgere 2D printr-o distributie liniara a efortului Reynolds.
Metoda energiei cinetice turbulente (TKE)
Efortul de forfecare se poate obtine din fluctuatiile vitezei turbulente cu ajutorul TKE. TKE se defineste ca:
(6)
u’, v’ si w’ reprezinta componentele vitezei pe cele trei axe. Dependenta liniara intre TKE si efortul de forfecare a fost prezentata de Townsend, 1976. Soulsby (1980) a dat relatia efortului de forfecare in functie de TKE:
(7)
de unde rezulta:
(8)
unde C1 reprezinta o constanta de proportionalitate. Valoarea C1 0,19 a fost utilizata de MacVicar si Roy (2007) in cazul curgerii peste un strat de pietris intr-un rau, de Rowinski et al. (2005) in cazul curgerii cu suprafat libera penste o frontiera rugoasa si de Pope et al. (2006) in cazul curgerii in rauri si in conditiile studiilor de laborator. Wolf (1999) a propus valoarea C1 0,19 bazandu-se pe studii ocenografice. Kim et al. (2000) a propus valoarea C1 0,21 in cazul curgerii in estuare pentru a se adapta la valorile obtinute, sugerand necesitatea studiilor ulterioare pentru a valida aceasta valoare ca valoare hidraulica universala.
Kim et al. (2000) au presupus o relatie liniara a TKE si variabile (denumite in continuare w’) si a sugerat faptul ca efortul de forfecare poate fi corelat cu componenta verticala a variabilei:
(9)
(10)
Kim et al. (2000) propune C2 0,9 comparand rezultatele obtinute pentru TKE w’ cu rezultatele obtinute cu ajutorul celorlalte metode.
Nezu si Nakagawa (1993) au demonstrat faptul ca efortul Reynolds si TKE sunt corelate in cazul curgerii cu suprafta libera. In stratul interior, coeficientul de corelare are o valoare apropiata de 0,1. Acesta lucru conduce la . In urma utilizarii metodei efortului Reynolds, extrapoland profilul TKE in zona stratului rugos, se obtine:
(11)
Metoda similitudinii la perete
Acest concept de similitudine se regaseste in cazul curgerilor turbulente in apropierea frontierei, pentru numere Reynolds ridicate si implica existenta unei zone din curgere in care energia turbulenta produsa si cea disipata sunt in echilibru relativ, iar difuzia este neglijabila, netinandu-se cont de conditiile impuse de curgerea in apropierea unui strat rugos. In stratul de echilibru (0,15 ≤ z/h ≤ 0,6) exista un echilibru intre energia produsa ρ si energia disipata ε (Townsend, 1976):
(12)
iar termenul difuziei energiei turbulente este:
(13)
Astfel, fluxul vertical al TKE in stratul de echilibru. Hurther si Lemmin (2000) au aratat faptul ca fluxul vertical al TKE normalizat cu viteza de frecare la puterea trei este dat de relatia:
(14)
Aceasta relatie tine cont de efectul stratului rugos.
(15)
Fk este o constanta cu valoarea de 0,3 pentru intervalul 0,15 ≤ z/h ≤ 0,6 (Lopez si Garcia, 1999; Hurther si Lemmin, 2009). In consecinta, o estimare a valorii vitezei de frecare se poate obtine cu ajutorul fluxului energiei turbulente in orice punct in intervalul 0,15 ≤ z/h ≤ 0,6. Cum aceasta plaja de valori este destul de mare, nu este necesar sa se cunoasca exact pozitia punctului de masurarepentru a determina valoarea lui u*. Aceasta metoda nu este sensibila la erori in ceea ce priveste determinarea fluxului de energie. De aceea, se preteaza la masurari realizate in natura.
Modelare experimentala
Utilizarea modelului fizic in laborator
In cadrul laboratorului de Instalatii al Facultatii de Inginerie a Instalatii din cadrul Universitatii Tehnice de Constructii Bucuresti, s-a realizat un model exprimental cu ajutorul caruia s-a studiat curgerea cu suprafata libera peste rugozitati artificiale (Figura 1b). Determinarile au presupus masurarea campurilor de viteza si calcularea diversilor parametri statistici caracteristici curgeii turbulente (abaterea medie pătratică pentru fluctuațiile vitezei în lungul curgerii respectiv pe direcția perpendiculară pe sensul curgerii și covarianța componentelor vitezei). Pentru aceasta s-a utilizat un sitem de masura de tip PIV (Particle Image Velocimetry), prelucrarea datelor obtinute fiind realizata in cadrul centrului de cercetare CAMBI.
In cele ce urmeaza, este descris standul experimental si principiul de functionare al acestuia.
a) b)
Figura 1 a) sedimente într-o conducta de canalizare reala, b) rugozități artificiale pentru modelul fizic
Forma depunerilor in straturi sedimentare in interiorul conductelor este complexa si variaza de la un caz la altul. Pentru modelul de laborator s-au utilizat rugozitati sub forma de sfera, din material plastic, prezentate in figura 1b . Aceasta configuratie a fost aleasa pentru a putea valida rezultatele obtinute in urma masurarilor prin compararea lor cu rezultate obtinute pe modele similare si prezentate in literatura [1], [2]. Tipul rugozitatii, in functie de criteriile prezentate de Perry [3], este prezentat in aceasta lucrare. Alte tipuri de rugozitati artificiale prezentate in literatura de specialitate sunt glaspapir [4], bare cu diverse sectiuni [5], impletitura de sarma, cuburi [6], etc.
Criterii de similitudine colector real/model experimental
Datorita fatului ca efectuarea de masurari PIV nu se poate realiza direct intr-un colector orasenesc, s-a realizat un model fizic in Laboratorul de Instalatii din cadrul Facultatii de Inginerie a Instalatiilor. Pentru a putea valida rezultatele obtinute, trebuie indeplinite criterii de similitudine intre cele doua colectoare, atat din punct de vedere al caracteristicilor geometrice, cat si din punct de vedere al parametrilor cinematici si dinamici.
In privinta similitudinii geometrice, rapoartele lungimilor caracteristice ale celor doua colectoare, sunt egale.
Prin similitudine cinematica se intelege faptul ca rapoartele vitezelor si ale acceleratiilor sunt constante. Prin similitudine dinamica se intelege faptul ca rapoartele fortelor corespunzatoare sunt constante. In cazul curgerii apei in conducte, prezinta interes raportul fortelor de inertie si vascozitate, eprimat prin criteriul Reynolds (16), si raportul fortelor de inertie si gravitatie, exprimat prin criteriul Froude (17).
(16)
(17)
este viteza medie pe secțiunea curgerii (m/s);
este coeficientul de vâscozitate cinematică (mp/s);
este o lungime ce caracterizează curgerea (m);
este accelerația gravitațională (m/s2);
Conform E. Iatan [1], indeplinirea ambelor criterii de similitudine implica si scalarea proprietăților fizice ale fluidelor ceea ce ar conduce la costuri foarte ridicate in realizarea modelului fizic, ori la utilizarea unor fluide care pur si simplu nu exista. Chanson (2009) propune ca in cazul curgerilor cu suprafața libera, criteriul Froude sa fie cel semnificativ datorita importantei forțelor gravitaționale. Pentru a contracara efectul forțelor de vâscozitate, Gill și Pugh (2009) propun utilizarea criteriului Froude coroborat cu o valoare minimă a numărului Reynolds de 2000 pentru curgeri cu suprafața libera.
Se notează:
(18)
raportul lungimilor considerate caracteristice, indicele n este pentru natură, indicele m este pentru model iar cu r se va nota raportul mărimilor caracteristice considerate.
Pornind de la criteriul Froude pentru similitudinea natură – model, adică se impune și se pot deduce următoarele relații pentru viteze și debite:
(19)
în care este raportul vitezelor medii pe secțiunile curgerilor din natură respectiv din modelul considerat,
(20)
în care este raportul debitelor natură – model.
Figura 2 – Caracteristici geometrice ale colectorului real
Pentru colectorul orasenesc ave urmatoarele caracteristici (Figura 2):
– diametrul interior al colectorului = 400 mm
– grosimea depunerilor de sedimente consolidate pe fundul canalului este de 20 mm,
– debitul de apă 10 l/s,
– raza hidraulică 0,068 m,
– panta de montaj a conductei 0,5%
Pentru colectorul modelului fizic, se alege diametrul maxim ce se poate prelucra pentru conducte de plexiglas, si anume diametrul interior de 144 mm. Debitul de apa evacuat prin colectorul modelului fizic se poate scala in functie de raportul ales și de criteriul Froude impus. S-a considerat ca lungime caracteristică raza hidraulică și se poate calcula raportul lungimilor caracteristice (coeficient de scară):
(21)
Utilizând relația de mai sus, stabilita pe baza criteriului Froude, se poate calcula debitul de apa evacuat prin colectorul modelului fizic, pornind de la coeficientul de scară Qr pentru debit.
Qr = Lr5/2 = 2,775/2 = 12,8 (22)
Qm = Qn / Qr = 10/12,8 = 0,78 l/s (23)
Figura 3– Caracteristici geometrice ale conductei utilizată ca model.
Se obțin următoarele caracteristici ale modelului:
– diametrul interior al conductei de plexiglas = 144 mm,
– grosimea depunerilor de sedimente consolidate pe fundul conductei = 7,5 mm,
– debitul de apă 0,78 l/s,
– raza hidraulică 0,0245 m,
– panta de montaj a conductei 0,5%.
În de mai jos sunt prezentate centralizat datele propuse pentru colectorul real și pentru model, valorile calculate și necunoscutele pe care dorim să le estimăm.
Tabel 1 – Datele propuse pentru colectorul real și pentru model
Standul experimental propus pentru dezvoltarea în cadrul tezei de doctorat
Descrierea standului experimental
Se considera un sistem in care se simulează curgerea gravitaționala cu suprafața libera intr-o conducta de canalizare cu depuneri, de secțiune circulara. Conducta are diametrul Dn150 (Di=144mm) si o lungime de 3,9m, fiind compusa din doua tronsoane egale. Îmbinarea celor două tronsoane s-a realizat cu o flanșă confecționată special tot din plexiglas și garnitură de cauciuc. Suprafața garniturii de cauciuc a fost pe cât posibil netezită pentru ca prezența acestei îmbinări să nu introducă perturbații suplimentare în curgere. Circuitul hidraulic, cuprinzând rezervor, pompă, debitmetru și conducte este reprezentat pe scurt în Figura 4.
Figura 4– Circuitul hidraulic al modelului experimental
Depunerile de sedimente din interiorul colectorului orasenesc au fost simulate prin realizarea unor placute din plastic montate la partea inferioara a conductei de evacuare a modelului fizic. Rugozitatile sunt simlate cu ajutorul unor sfere din plastic montate pe aceste placute. Sferele au diametrul de 4,5 mm si simuleaza depunerile de sedimente consolidate de la partea inferoara a colectorului. Sferele sunt calibrate si sunt montat in perforatii ale placutelor (Figura 5). Pentru fixarea sferelor s-a utilizat un adeziv compus din solutie de acetona si di-cloretan. Intrucat s-a observat ca solutia adeziva topeste suprafata din plastic a placutelor, gaurile au fost realizate cu un diametru de 4,2 mm. In acest fel, sferele reusesc sa fie ingolbate in placutele din plastic pana in dreptul diametrului acestora. Configuratia aleasa pentru simularea rugozitatii cu sfere din plastic se datoreaza faptului ca, pentru validare, rezultatele obtinute au fost comparate cu rezultate din literatura de specialitate, obtinute pe configuratii asemanatoare.
Figura 5 – Placute perforate cu sfere de diametru 4,5mm
Plăcuțele au fost aliniate la partea inferioară a conductei modelului din plexiglas astfel încât generatoarea inferioară a conductei să se suprapuna cu linia mediană trasată la nivelul plăcuțelor. În acest scop s-a realizat dispozitivul din Figura 6. De asemenea, s-a încercat realizarea îmbinării la acelasi nivel între plăcuțele succesive, pentru a nu influenta curgerea apei.
Figura 6 – Dispozitiv utilizat pentru alinierea plăcuțelor de plexiglas la partea inferioară a conductelor din model
Principiul de functionare al dispozitivuluide aliniere ce cuprinde o nivelă fixată pe un disc rigid de plastic cu diametrul foarte apropiat de diametrul interior al conductei de plexiglas este urmatorul:
-pe acest disc este marcat punctul cel mai de jos; in momentul în care cele două ″bule″ ale nivelei sunt centrate se marchează pe conductă un punct care corespunde punctului cel mai de jos marcat pe discul rigid; punctele succesive astfel desenate pe conductele de plexiglas ne-au ajutat să materializăm cu o sfoară generatoarea inferioară a conductei și să aliniem cu precizie plăcuțele de plexiglas, a se vedea în Figura 7, pașii descriși în acest paragraf.
Cele opt plăcuțe au fost prevăzute cu ″piciorușe″, poziționate la extremități, care au fost decupate astfel încât să urmărească curbura conductei. În zona acestor ″piciorușe″ s-a utilizat silicon pentru a realiza fixarea de conducta de plexiglas. Pentru anularea efectului introdus în curgere de muchia ce apare în lateralul plăcuțelor, la îmbinarea cu conducta de plexiglas s-a rotunjit muchia cu ajutorul unui strat de silicon, relativ rezistentă la apă. [1]
Figura 7 – Alinierea plăcuțelor de plexiglas relativ la generatoarea inferioară a conductelor
Din cauza curburii conductei de plexiglas și datorită diferenței de mediu, aer la exteriorul conductei respectiv apă la interiorul acesteia, în zona în care se realizează măsurătorile conducta a fost introdusă într-un “acvariu” umplut cu apă. Astfel, în zona de măsură, conducta de plexiglas a avut apă și în exteriorul ei.
Pentru asigurarea unui debit constant in circuitul standului experimental, s-a montat o pompa de circulație marca DAB, tipul A 110/180XM. Pompa cu turație fixa are curba caracteristica de funcționare conform Figura 8.
Figura 8– Pompa circulatie DAB model A110/180XM
Rezervorul tampon montat pe aspirația pompei de circulație este de tip Valrom, cu un volum total de 300 l. Dimensiunile acestuia sunt L x l x H = 880mm x 680mm x 950mm.
Figura 9– Rezervor tampon deschis 300 l
Pentru reglarea debitului s-a utilizat o vana de reglaj cu debitmetru, de diametru Dn50 (2”), marca CALEFFI 132, prezentata in fig.16. Aceasta vana s-a montat pe circuitul conductelor din Cu. Se urmărește realizarea măsurărilor pentru diferite valori ale debitului.
Vana de reglaj cu debitmetru are o plaja de fixare a valorilor debitului cuprins intre 50 si 200 l/min.
Figura 10 – Vana de regalj cu fluxometru CALEFFI
Figura 11– Pozitia vanei de regalj cu fluxometru CALEFFI in cadrul modelului experimental
Pentru măsurarea exacta a debitului, s-a folosit un debitmetru cu ultrasunete de tip SIEMENS SITRANS F (Figura 12). Acesta s-a montat pe conducta din cupru Dn40, in amonte de vana de reglaj cu fluxometru, la o distanta de peste șase diametre fata de orice fitting sau organ de reglaj pentru a evita erorile de măsurare (Figura 13).
Figura 12– Debitmetru cu ultrasunete SIEMENS SITRANS F
Figura 13– Montaj debitmetru cu ultrasunete
Pentru a simula curgerea gravitaționala intr-un sistem de canalizare s-a folosit un cămin de liniștire. Acesta este de tip Valrom 200/160, cu o inaltime H=1130 mm. Curgerea se face pompat din rezervorul tampon pana in căminul de liniștire, iar de aici, mai departe, gravitațional prin conductele de plexiglas.
Figura 14– Camin din polietilena
Principiul de funcționare al modelului experimental
Modelul experimental simulează curgerea cu suprafața libera intr-un sistem de canalizare cu depuneri. Depunerile sunt simulate cu ajutorul unor sfere cu diametrul de 4,5mm fixate pe placi lise din material plastic.
Rezervorul tampon de 300 l si căminul de liniștire se umplu cu apa pana când aceasta ajunge la nivelul cotului Dn150 din partea superioara a căminului de liniștire.
Figura 15– Camin de linistire
Standul experimental prezintă un circuit închis al apei. In conductele din Cu Dn40 apa curge sub presiune, debitul fiind asigurat de pompa de circulație. Pentru a trece de la curgerea sub presiune la o curgere gravitaționala cu suprafața libera, se montează un cămin de liniștire. In momentul in care nivelul apei depaseste nivelul inferior al cotului de la intrarea in sistemul de conducte din plexiglas, curgerea devine gravitaționala. Debitul evacuat prin conducta de plexiglas este egal cu debitul vehiculat in interiorul conductei de Cu Dn40.
Pentru a regla si măsura debitul s-a folosit o vana de regalare cu fluxostat Dn50. Acestea funcționează după următorul principiu: se stabilește debitul dorit si se închide/deschide treptat vana de reglaj pana in momentul in care valoare citita a debitului corespunde valorii setate.
Prin utilizarea standului experimental descris mai sus se obține o curgere gravitaționala stabila in interiorul conductelor. Pentru a nu influenta in vreun fel curgerea gravitaționala (eventuale vibrații datorate funcționarii pompei de circulație) sistemul de conducte gravitaționale este izolat din punct de vedere al suporților si susținerilor fata de restul conductelor din Cu si echipamentelor montate.
Figura 16– Stand experimental – elemente de susținere a componentelor
Pentru determinarea câmpurilor de viteze se injectează la partea inferioara a conductei un amestec de apa cu particule trasoare. Acestea sunt sfere cu pereți subțiri din sticla, imbracate in argint, cu diametrul D50 egal cu 10 μm. Pentru o evaluare a înregistrărilor PIV este necesara prezenta unui număr suficient de mare de particule trasoare in zona de interes.
Principii de măsură și echipamente disponibile
Majoritatea curgerilor întâlnite în domeniul ingineresc sunt turbulente. Mărimile utilizate pentru caracterizarea lor, atât scalare cât și vectoriale prezintă de cele mai multe ori repartiții spațiale puternic tridimensionale influențate de prezența structurilor de tip vârtej caracteristice – începând cu vârtejurile primare de scară mare, coerente și până la scările turbulente cele mai mici. În aceste situații posibilitatea rezolvării ecuațiilor Navier-Stokes rămâne destul de restrânsă, chiar și în prezent, în condițiile dezvoltării importante a tehnicii de calcul și a modelelor de tip CFD. Importanța investigării experimentale a curgerilor rămâne deci primordială, cu atât mai mult cu cât orice model nou sau extindere a aplicării lui va necesita întotdeauna o validare experimentală.
La fel ca și în cazul tehnicilor CFD, mijloacele metrologice de astăzi au cunoscut o dezvoltare importantă, o adevărată explozie de posibilități, facilitată de miniaturizarea circuitelor integrate și de scăderea prețului lor de fabricație. Dacă până acum aproximativ zece ani principalele mijloace de măsurare a câmpului de viteză într-un fluid erau reprezentate de tubul Pitot și anemometrul cu fir cald [1, 2], în prezent mijloacele optice de măsura s-au „democratizat”. Într-adevăr, măsurările de Viteza cu ajutorul Imaginilor de Particule – Particle Image Velocimetry (PIV) – sau cu ajutorul Efectului Doppler Laser – Laser Doppler Velocimetry (LDV) nu mai reprezintă astăzi mijloace rezervate celor mai mari și mai bogate în fonduri dintre laboratoarele de cercetare. Dezvoltarea captorilor CCD și CMOS și miniaturizarea lor au avut drept rezultat faptul că mai ales tehnicile de vizualizare (Particle Image Velocimetry, tomografie, Particle Traking Velocimetry, holografie, etc.) au devenit foarte accesibile pentru experimentatorii din domeniul mecanicii fluidelor. Acest lucru a fost posibil și datorită faptului că pe lângă dezvoltarea senzorilor optici, a scăzut prețul de producție a laserilor de tip pulsat, de putere medie și mare. Consecința este legată de oferirea unor noi perspective de investigare a câmpurilor bidimensionale și tridimensionale ale vitezei de dimensiuni cu mult mărite față de ceea ce era posibil în urma cu câțiva ani numai [3], frecvențele de achiziție maxime apropiindu-se de rezolvarea unor fenomene întâlnite în aplicații practice pe o scară din ce mai largă.
De exemplu, studii experimentale cu ajutorul tehnicii de măsurare PIV, cu mijloace standard pentru acea vreme, erau aplicate în perioada 2006-2007 [4-7] în laboratorul francez LEPTIAB unor curgeri de tip jet de dimensiuni reduse, câmpurile de viteză instantanee având dimensiuni medii de câteva zeci de cm², o rezoluție de aproximativ 1-2 pixel²/mm² și o frecvență de înregistrare de maxim 15Hz. Acești parametri nu permiteau „rezolvarea” curgerii din punct de vedere spațial și cu atât mai puțin din punct de vedere temporal, frecvențele naturale ale stratului limită de tip jet având valori de câteva zeci de Hz pentru numere Reynolds relativ mici [8, 9]. Astăzi la UTCB avem un sistem PIV ce permite achiziția unor câmpuri de viteză instantanee de până la aproximativ 60cm x 60cm [10] iar la LEPTIAB un sistem stereoscopic PIV permite înregistrarea a până la 500 de câmpuri instantanee de viteză cu toate cele trei componente de dimensiuni de peste 100cm² [11]. Putem spune astfel că anumite metode experimentale de investigare a curgerilor « rezervate » investigațiilor la scară mică « de laborator », devin astăzi posibile și în cazul studiilor experimentale realizate în condiții cât mai reale.
Cercetările anterioare legate de dinamica fluidelor, în care măsurarea vitezelor se realiza cu ajutorul sondelor de presiune sau a anemometrelor cu fir cald, ridicau problema intruzivității senzorilor în curgerile studiate. Astfel, problemele ce apar în folosirea acestor metode intruzive sunt cele legate de introducerea unor perturbații ce pot genera la rândul lor instabilități, vârtejuri și deteriorează condițiile de reproducere a fenomenelor naturale [2, 12-15]. Mai mult, sondele intruzive necesită calibrare și sunt sensibile la factorii externi (temperatură, umiditate etc.) [16].
Odată cu apariția laserelor în a doua jumătate a secolului XX, metodele non-intruzive de măsurare au devenit mult mai practice. După introducerea laserelor cu gaz, tehnica de măsură cu ajutorul efectului Laser Doppler – Laser Doppler Velocimetry (LDV )- a fost dezvoltată în primul rând de către cercetătorii americani Yeh și Cummins [17] . Acesta a fost unul dintre cele mai importante progrese pentru investigarea curgerii fluidului. Rezultatul este liniar, nu sunt necesare etalonări, nivelul de zgomot este cu mult mai redus, răspunsul este de frecvență înaltă și viteza este măsurată independent de alte variabile. În ultimele trei decenii, tehnica LDV a cunoscut progrese semnificative în ceea ce privește metodele optice, cum ar fi fibra optică, precum și tehnicile avansate de procesare a semnalului și dezvoltarea de software. În plus, metoda LDV a fost extinsă, devenind Doppler Phase Technique pentru măsurarea dimensiunilor particulelor și bulelor, odată cu viteza.
Dezvoltarea uimitoare a tehnologiilor laser și a senzorilor CCD și CMOS a deschis posibilitatea de investiga calitativ mișcarea fluidului (vizualizarea curgerii) și mai târziu de a analiza din punct de vedere cantitativ câmpul de curgere. Prin dezvoltarea tehnicii PIV a devenit una dintre cele mai populare instrumente pentru investigarea curgerii de fluid, având numeroase aplicații.
Ulterior, injectarea de particule fluorescente în lichide sau de vapori de acetonă în aer a făcut posibilă măsurarea concentrațiilor și temperaturilor în paralel cu măsurarea câmpului de viteză. Astfel a apărut tehnica Laser Induced Fluorescence (LIF) .
Fluidul este însămânțat cu particule fine și cu un marker fluorescent pentru a transmite semnalele PIV și respectiv LIF. Camera PIV, echipată cu un filtru ce corespunde lungimii de undă de excitație a moleculelor fluorescente de către lumina laser, detectează deplasarea particulelor, iar camera LIF, echipată cu un filtru ce corespunde lungimii de undă a semnalului fluorescent, detectează concentrația de molecule fluorescente.
Astfel, cu ajutorul acestor tehnici, se pot determina câmpuri de viteză, de concentrație, de temperatură, de turbulență, dimensiuni ale particulelor, iar cu ajutorul programelor de post-procesare a datelor se pot determina și alte mărimi ce derivă din aceste cantități măsurabile.
Tehnicile de investigare a curgerilor își găsesc aplicabilitatea într-o gamă largă de domenii din mecanica fluidelor. Pornind de la curgeri în spații de câțiva microni (micro-canale) și ajungând la ingineria vântului sau la curenții din oceane, întâlnim echipamente de măsură dedicate fiecărui caz.
Datorită răspunsului de mare frecvență (chiar sute de kHz) a senzorilor cu fir cald (Hot Wire Anemometry – HWA) acest tip de investigare devine unealta ideală pentru studiul mișcărilor turbulente sau măsurări de analiză spectrală. Astfel, aplicațiile din domenii legate de fenomenele naturale, protecția mediului, de turbulență în general folosesc cu predilecție măsurările cu fir cald.
Investigațiile cu Imagini de Particule, PIV, sunt preferate în domeniile în care trebuie măsurați sute sau chiar mii de vectori simultan și unde sunt redate câmpuri mari de viteze. Ca și metoda cu laser Doppler, LDV, fiind metode de investigare neintruzive, sunt folosite atunci când curgerea fluidului nu trebuie să fie influențată de factori exteriori: în aeronautică, domeniul biomedical, combustie, științe naturale, studiul spray-urilor etc. Aceste metode sunt folosite și pentru validarea modelelor CFD (Computational Fluid Dynamics).
Probleme legate de difuzia luminii și particule
Injectarea controlată de particule în curgerea studiată permite o analiză cantitativă și calitativă în ceea ce privește talia, distribuția și concentrația acestora. În general, aceste particule are trebui să fie îndeajuns de mici pentru a urmări fidel curgerea și îndeajuns de mari pentru a reda un semnal luminos suficient de bun pentru captor. De asemenea trebuie să fie ne-toxice, ne-corozive și inerte din punct de vedere chimic dacă este posibil. Melling [18] trece în revistă majoritatea particulelor trasoare care au fost utilizate în experimente PIV, precum și metodele de generare și de introducere ale acestora în curgerile de fluid.
Alegerea modalității și a particulelor injectate depind de o serie de parametri. În primul rând natura particulelor ar trebui să fie aleasă în funcție de curgerea ce urmează să fie analizată și de echipamentul laser disponibil. În general, particulele trebuie să fie alese cât mai mari posibil, cu scopul de a difuza cât mai multă lumină, însă dimensiunea lor este limitată deoarece o talie prea mare a particulei nu va urmări curgerea în mod corespunzător. Astfel, mărimea maxim admisibilă a particulelor scade cu creșterea vitezei de curgere, turbulenței și gradientului de viteză.
Ideal, materialul de injectare cu particule ar trebui în primul rând să fie ales astfel încât
acestea să fie neutre din punctul de vedere al inerției față de fluidul transportor, acest aspect fiind de multe ori lăsat în plan secund. În al doilea rând ar trebui să fie luat în considerare modul de injectare a particulelor.
În cazul curgerilor de lichide, comparativ cu cele gazoase, spectrul materialelor folosite pentru însămânțare este mult mai mare, fapt datorat incintei închise în care se realizează curgerea. În cazul curgerilor de aer sunt folosite: fum pentru spectacole, diferite tipuri de ulei pulverizat, apă, dioxid de titan (TiO2), oxid de aluminiu(Al2O3) etc. Generatoarele de fum pentru spectacole sunt ieftine și generează suficiente particule. Uleiul poate fi pulverizat folosind dispozitive cum ar fi duza Laskin, generatoare de particule de talie de un micron până la sub un micron, fiind deosebit de utile pentru aplicații de mare viteză. Dioxidul de titan (TiO2) și oxidul de aluminiu (Al2O3) sunt utile pentru aplicații la temperaturi ridicate, cum ar fi combustia și măsurări ale flăcării.
Concentrația naturală de particule foarte mici, este adesea mult
mai mare decât cea a particulelor utile pentru măsurare. În unele cazuri, în general pentru lichide, aceasta cauzează o creștere nedorită a nivelului de zgomot ca urmare a semnalelor incoerente ale particulelor de talie mică. În general, se recomandă, ori de câte ori este posibil, să se controleze mărimea și concentrația particulelor prin filtrarea fluidului și, ulterior, prin adăugarea de particule determinate ca dimensiune.
Talia particulelor influențează de asemenea și lumina difuzată de acestea. Când o undă luminoasă străbate un mediu, câmpul electromagnetic al undei interacționează cu particulele aflate în mediul respectiv, energia undelor fiind absorbită de acestea și apoi re-emisă, lumina fiind astfel împrăștiată (difuzată) în toate direcțiile (Figura 23). Teoria clasică a difuziei luminii a fost fondată de Lordul Rayleigh, fiind numită difuzie de tip Rayleigh. Aceasta teorie este însă aplicabilă particulelor mici, adică având o dimensiune cu mult mai mică față de lungimea de undă a luminii difuzate (diametrul este mai mic decât raportul λ /10 ). Intensitatea luminii difuzate în acest caz este invers proporțională cu puterea a patra a lungimii de undă (I ~ 1/λ4). Pentru particule mari (din punct de vedere optic), se aplică teoria Rayleigh-Gans, intensitatea luminoasă variind în același mod, adică este invers proporțională cu puterea a patra a lungimii de undă (I ~ 1/λ4). Dacă o particulă este mai mare decât lungimea de undă, lumina poate fi difuzată diferit în funcție de unghiul de observație. Difuzia se numește de tip Mie, iar intensitatea luminii difuzate este invers proporțională cu puterea a doua a lungimii de undă : I ~ 1/λ2.
Figura 17 – Distribuția intensității luminii difuzate (lungime de unda λ) pentru o sferă de rază: a) r<< λ ; b) r= λ; c) r> λ [19]
Particulele folosite în vizualizările laser au diametre între 1 μm și 10 μm, cum ar fi fumul sau particule de praf, produc difuzie de tip Mie când sunt iluminate în spectrul vizibil.
Figura 18 – Vizualizarea unei curgeri cu ajutorul unei pânze de lumină laser
Diverse particule, cum ar fi fumul, pot fi injectate în curgere pentru a urmări liniile de curent. Aceste particule pot fi iluminate cu o pânză laser (Figura 24) pentru a vizualiza o secțiune din curgerea complexă a fluidului: presupunând că particulele urmăresc îndeaproape liniile de curent, se pot realiza investigări în general calitative. Pentru a obține și o investigare cantitativă se folosesc tehnicile ce vor fi descrise ulterior. Gratie unei cadențe ridicate de înregistrare a imaginilor, aceasta tehnică reprezintă o abordare calitativă foarte utilă în observarea curgerilor turbulente și permite direcționarea către o analiză cantitativă a acestora.
Mai mult, detectarea contururilor prin proceduri de procesare a imaginilor permite accesul la informații cantitative ale dinamicii turbionare, prelevate de altfel cu firul cald.
Figura 19 – Schema principiului de funcționare pentru vizualizări cu tomografia laser după [20]
Principiul este simplu: se iluminează o secțiune a curgerii cu ajutorul unei pânze de lumină laser foarte puternică (Figura 25). Această secțiune este mai apoi filmată cu ajutorul unei camere echipate cu o memorie internă suficient de mare și captori de tip CMOS cu frecvență înaltă de înregistrare.
Figura 20 – Dispozitiv experimental pentru vizualizarea curgerilor prin tomografie rapidă
Măsurări de viteză cu Imagini de Particule plană cu două componente PIV 2D
Măsurările cu ajutorul Imaginilor de Particule (PIV) reprezintă o metodă optică de investigare a curgerilor, mai precis a câmpurilor de viteză (2D în cazul de față). Principalul avantaj față de alte tehnici de investigare a curgerilor (ex. LDV) este capacitatea sistemului PIV de a măsura un câmp întreg de viteze într-un singur pas.
Această tehnică se bazează pe intercorelarea de imagini ale unei curgeri, înregistrate de senzori de tip CCD sau CMOS. Curgerea este însămânțată în prealabil cu particule fine solide sau lichide. Principiul esențial al acestei metode de măsură este determinarea vitezelor locale ale curgerii pornind de la deplasările locale ale particulelor.
Dacă în timpul unui interval de timp foarte scurt ∆t, o particulă se deplasează din poziția la poziția , viteza locală de deplasare poate fi exprimată:
(25)
În acest scop se înregistrează semnalul Mie difuzat de particule, pe două imagini succesive, separate în timp cu ∆t. Se aplică apoi un tratament statistic spațial de intercorelare asupra imaginilor digitalizate în funcție de nivelul lor de gri.
Pentru a putea obține vectorii locali de viteză, se împart imaginile în ferestre mici numite rețele de discretizare. Cu cât discretizarea este mai fină, cu atât câmpul de viteză este mai bine rezolvat în spațiu.
Figura 21 – Auto-corelarea după [21]
Deplasarea este evaluată în subdomeniile definite de grila de calcul, prin determinarea poziției maximului de corelare (vârf maxim de corelare) (Figura 27).
Pentru ca deplasarea cea mai probabilă calculată să fie corectă, trebuie ca fiecare celulă din rețeaua de discretizare să conțină un număr suficient de particule. Numărul critic de particule este de 8 într-o celulă. Măsurările sunt corecte atunci când deplasarea medie, în pixeli, nu depășește sfertul numărului de pixeli ai dimensiunii celulei grilei de calcul.
Vectorul viteză pentru celula j este vectorul viteză cel mai probabil al acestei celule. Factorul de probabilitate depinde de raportul Semnal / Zgomot (SNR – Signal to Noise Ratio). Cu cât valoarea lui SNR va fi mai mare cu atât vectorul calculat va fi mai probabil. Digitalizarea imaginii induce o discretizare a definiției vârfului de corelare.
Figura 22 – Schema de principiu PIV, după DANTEC
Valoarea deplasării este egală cu un număr întreg multiplicat cu pasul de pixeli ai imaginii digitale (peak locking sau « blocaj de vârf »). Această deplasare influențează dinamica măsurărilor și limitează accesul micilor deplasări. Originea « blocajului de vârf » provine din variația nivelului de gri pe un singur pixel. Una din soluții constă în mărirea taliei particulelor. Se arată că un număr de cel puțin 3 pixeli este necesar pentru reprezentarea corectă a unei particule [22]. O altă soluție pentru evitarea « blocajului de vârf » este optimizarea determinării poziției maximului vârfului de corelare prin interpolare. Aceasta se realizează cu ajutorul metodelor așa-numite de precizie « sub-pixel » dintre care cele mai utilizate folosesc ajustarea vârfului cu o funcție gausiană (cea mai apropiată de vârful de corelare).
Tehnica de măsură PIV este perfect adaptată pentru vizualizarea câmpurilor instantanee de viteză a curgerii. Totuși, ținând cont de frecvențele limitate ale laserelor Yag, sistemele PIV clasice nu permit accesul la rezoluția temporală a curgerii turbulente, motiv pentru care câmpurile de viteză 2D măsurate sunt necorelate temporal între ele. Alegerea dimensiunilor particulelor utilizate pentru măsurări PIV reprezintă adesea un compromis între necesitatea unui timp scăzut de relaxare dat de o talie redusă a particulei și o difuzie suficientă a luminii date de iluminarea în retro-difuzie sau de o sursă luminoasă de mare intensitate. Pentru a respecta cel de-al doilea criteriu este mai puțin costisitoare mărirea calității particulelor față de mărirea intensității sursei luminoase.
Echipamente optice de măsură disponibile la Laboratorul de Instalații de la UTCB
Pentru realizarea măsurărilor PIV se propune următoarea configurație a echipamentelor de măsura:
Figura 23 – Schema de principiu pentru masurarile PIV
La Laboratorul de Instalații din cadrul Facultății de Inginerie a Instalațiilor, a Universității Tehnice de Construcții București sunt disponibile la ora actuală două sisteme optice de măsură:
Un sistem PIV clasic, compus dintr-o cameră de înalta sensibilitate FlowSense MKII 4M, cu captor CCD, având o rezoluție de 4 * 106 pixeli și dintr-un laser Litron de 200mJ, ce produce un plan luminos cu lungimea de undă de 532nm. Frecvența de achiziție a sistemului este de 7.5Hz.
Un sistem PIV rapid, compus dintr-o cameră Nanosense, cu captor CMOS, de rezoluție 512 * 512 pixeli, si dintr-un laser Nanopower ce produce un plan luminos cu o lungime de 765nm. Frecvența maximă de achiziție a acestui sistem este de 3.5 KHz.
Sistemul PIV utilizat
Sistemul PIV utilizat a fost cel cu camera de tip Nanosense MKII, cu un obiectiv Nikon cu o distanta focala de 50 mm și o deschidere de 1.2. Această cameră are o frecvență de achiziție de 5000 Hz pentru o imagine de 512 x 512 pixels. Imaginile sunt realizate pentru frecvențe de 300, 500 și 1000 Hz. Numărul imaginilor înregistrate variază între 1000 și 3000. Imaginile stocate pe memoria internă a camerei sunt transferate la finalul fiecărei achiziții pe PC. Pentru iluminarea particulelor trasoare din cadrul curgerii (particulele sunt sfere de diametru 10 μm) se folosește un laser cu infraroșu Nanopower, având puterea de 4W producând o lungime de undă de 795 nm.
Alegerea noastră s-a îndreptat către acest sistem dată fiind posibilitatea de surprindere a dinamicii curgerii studiate datorită frecvențelor de achiziție importante. Sistemul PIV clasic disponibil la centrul de cercetare CAMBI, a fost utilizat în cadrul tezei de doctorat al dnei. Dr. Ing. Elena Iatan [5].
Cea mai comună metodă a tehnicilor PIV este de a înregistra două imagini consecutive asupra curgerii fluidului. În lichidul de lucru se introduc particule trasoare care urmăresc mișcarea locală a fluidului, fără o influență majoră a acestuia. Între cele două imagini există un timp de decalare, Δt, foarte mic, ce poate varia intre câteva zeci si câteva sute de μs. Imaginile ce conțin câmpul de particule trasoare sunt împărțite în zone elementare (de obicei pătrate cu latura de 32 pixeli) pentru o determinare probabilistică a deplasării locale a particulelor. Prima estimare a vitezei locale a fluidului este obținută astfel:
In cadrul măsurărilor s-au folosit ca particule trasoare sfere cu pereți subțiri din sticla, imbracate in argint, cu diametrul D50 egal cu 10 μm. Pentru o evaluare a înregistrărilor PIV este necesara prezenta unui număr suficient de mare de particule trasoare in zona de interes.
În general, pentru măsurările PIV, controlul cantității de particule introduse în curgere reprezintă un aspect delicat pentru că pe lângă omogenitatea însămânțării trebuie asigurată o densitate de particule de minim 8 particule pe celula de interogare. Mai mult, întotdeauna trebuie examinat numărul mediu de pixeli asociați unei particule. Astfel, odată ce obiectivul camerei și câmpul de măsură sunt fixate, trebuie verificat dacă particula măsoară în jur de 2 x 2 pixeli. Într-adevăr, cu cât numărul de pixeli este mai mare, vârful de corelare este amplificat și precizia de deplasare medie a particulelor poate fi alterată. Cu cât numărul este mai mic, fenomenul de « blocaj de vârf » poate apărea, ceea ce introduce erori de calcul în metodele de evaluare sub-pixel [6, 7].
Figura 24– Sistem utilizat pentru masurarile PIV
Sistemul PIV este folosit în cadrul lucrării pentru determinarea distribuțiilor de viteze, identificarea punctului de desprindere a stratului limită și determinarea liniei suprafeței libere.
Pentru măsurarea distribuțiilor de viteze s-a folosit un număr de 350 perechi de imagini înregistrate cu frecvența de 1000 Hz.
Etapele în post-procesarea înregistrărilor PIV sunt următoarele:
Achiziționarea de imagini succesive la un interval de timp Δt;
După achiziționare este necesară îndepărtarea imaginilor foarte luminoase și a imaginilor cu puține particule;
Zonele de neinteres se vor ”masca”, astfel se va oferi procesarea cu mai puține erori;
Analiza se bazează pe mai multe imagini (350 de perechi de imagini);
Imaginile sunt împărțite în zone egale. Aceste zone au o rezoluție de 16 x 16 pixeli, unde sunt urmărite deplasările particulelor;
Se obțin 360 de câmpuri instantanee de viteza;
Se calculează media statistica a câmpurilor de viteza, obținându-se vectorii medii precum si abaterile pătratice (rms) ale acestora.
Etalonarea imaginilor a rezultat într-o rezoluție spațială de 0.07mm/pixel ceea ce corespunde unui câmp de vizualizare de 36mm 29mm. Cele 350 perechi de imagini au fost procesate cu ajutorul unui algoritm de corelare adaptativă pe mai multe niveluri de grilă ținând cont de deformarea celulelor și deplasarea sub-pixel. Grila finală este compusă din celule de 1616 pixeli cu 75% acoperire ceea ce poate fi tradus ca fiind un vector la fiecare 4.12 pixeli sau, o rezoluție spațială de 0.36 mm.
În ambele cazuri metoda de corelare a fost validată la fiecare nivel dacă valoarea SNR a corelării este mai mare față de un prag fixat la 1.1. În medie, mai puțin de 2% din vectori au fost detectați ca fiind invalizi. Aceștia au fost corectați printr-o interpolare biliniară. O inspecție sistematică a histogramelor de deplasare a particulelor arată o distribuție bimodală bine discretizată (fără prezența « blocajului de vârf ») [8].
Precizia algoritmilor de corelare din în cadrul software-ului Dynamic Studio utilizat, a fost evaluată de Calluaud [9] plecând de la o imagine reală a unei curgeri în care au fost injectate particule solide cu traiectorie impusă. Pentru algoritmii “convenționali” de corelare ce au fost utilizați în cadrul acestui studiu, s-a găsit o rezoluție de aproximare sub-pixel de ordinul a 1/64 pixeli. Pentru estimarea erorii sistematice pentru măsurările PIV cu acest sistem a fost utilizată ca referință o tehnică LDV (Laser Doppler Velocimetry) a cărei precizie a fost prealabil validată. Acest studiu a fost realizat în cadrul unor lucrări anterioare [10]. Astfel eroarea relativă nu depășește 3% pentru component de viteza în lungul direcției principale a curgerii. Pentru cealaltă component de viteza, date fiind valorile foarte mici, eroarea poate atinge 35%.
Protocol de măsură
S-au realizat măsurări în condiții de laborator utilizând standul experimental și sistemul de măsură descrise în paragraful precedent. Configurațiile studiate și planurile în care s-au realizat măsurătorile PIV sunt prezentate în Figura 25. Sensul curgerii este în lungul axei X iar axa Y este perpendiculară pe sensul curgerii.
configurația 1, p/k = 4 configurația 2
Figura 25. Configurații studiate pentru poziția rugozităților artificiale
Tabelul 2
În tabelul 2 sunt prezentate valorile celor 3 debite la care s-au realizat măsurările PIV.
Validarea rezultatlor experimentale
Campuri de viteza
Așa cum precizam mai devreme, pentru determinarea câmpurilor de viteza prin metoda PIV este necesară introducerea particulelor de argint în volumul de apa. Metoda de măsura PIV se bazează în mod indirect pe vizualizările obținute prin iluminarea acestor particule cu ajutorul sursei laser. Aceasta sursă este prevăzută cu un modul optic realizat din lentile speciale ce generează o pânză fină de lumină.
Figura 26. – Principul de injecție a particulelor în domeniul de măsură
Pentru realizarea vizualizărilor au fost injectate particule la partea inferioara a colectorului așa cum este arătat în Figura 5. Pentru configurația cu rugozități artificiale, au fost studiate cele două planuri longitudinale din Figura 4 pentru cele trei debite din Tabelul 1, respectiv 0,7 l/s, 0,9 l/s si 1,1 l/s.
Planul luminos generat de sursa laser a fost amplasat, inițial în planul 1 de măsură, adică în planul median al canalului, ceea ce revine la o configurație de măsură între două rânduri de semisfere. Intr-o a doua etapă planul a fost repoziționat în configurația 2 de măsură pe primul rând de semisfere adiacent planului median.
Figura 27. – Subdivizarea zonei de măsură în două regiuni cu parametri de corelare diferiți
Așa cum era de așteptat, din literatură și pe baza studiile anterioare ale dr. ing. Elena Iatan [5], dinamica curgerii peste semisfere este una foarte complexă. Sistemul PIV utilizat, având o rezoluție mai mică decât cel din studiile anterioare, necesită pentru operațiunile de corelare a imaginilor de particule, subdivizarea câmpului vizualizat în două regiuni în care parametrii de corelare să fie ușor diferiți. Această împărțire este prezentată în Figura 33.
Prezentăm în cele ce urmează (Figurile 34 – 39) distribuții tipice ale celor două componente medii ale vitezei în planurile de măsură pentru diferitele configurații studiate. Vectorii sunt colorați în funcție de magnitudinea lor.
Rezultatele obtinute sunt comparate cu date din literatura. Agelinchaab si Tachie (2006) au realizat un studiu privind curgerea apei intr-un canal deschis peste o frontiera rugoasa. In studiul de laborator s-au utilizat trei configuratii cu sfere montate pe fundul colectorului, dupa cum sunt prezentate in fugura [ ]:
Comparatia datelor obtinute in cadrul tezei se realizeaza pentru valorile rezultate in configuratia p/k=4, asemanatoare cu configuratia modelului de laborator. Mai jos este prezentata configuratia rugozitatilor pentru modelul fizic utilizat pentru masurarile din cadrul prezentei teze.
In imaginile urmatoare, se pot observa vectorii de viteza ce caracterizeaza curgerea, atat pentru zona mare, cat si pentru zona intre sfere.
In figurile anterioare se pot observa vectorii viteza pentru diferite viteze de curgere intr-un canal cu rugozitati. Liniile de câmp, in zona de deasupra sferelor, sunt paralele in sensul de curgere al apei. In zona dintre sfere se observa formarea vârtejurilor si o circulație in sens opus curgerii in canal.
Pentru o validare vizuala a modelului de laborator, s-au comparat imaginile obtinute in urma prelucrarii datelor cu imagini din literatura. Astfel, se observa ca profilurile de viteze ce se formeaza in sensul curgerii sunt asemanatoare cu cel obitnute in cadrul studiului Agelinchaab si Tachie (2006). De asemenea, se poate urmari si formarea vartejurilor intre doua sfere succesive. Zona de formare a vartejurilor corespunde cu rezultatele studiilor similare.
Figura 28. – Configurația 1 – Plan 2 – Vectori de viteza pentru debitul de 0,7l/s in zona de deasupra semisferelor
Figura 29 – Configurația 1 – Plan 2 – Vectori de viteza pentru debitul de 0,7l/s in zona dintre semisfere
Profiluri de viteza
Pentru a tesata daca punctul din curgere in care sunt extrase datele utilizate are o influenta asupra acestora, s-au realizat profilurile de viteza pe doua sfere succesive si s-a constatat faptul ca cele doua profiluri de viteza se suprapun aproape in totalitate.
Profilurile de viteza s-au realizat pentru toate cele trei debite, respectiv 0,7 l/s, 0,9 l/s si 1,1 l/s.
Având in vedere faptul ca profilurile de viteza sunt aproape identice, nu trebuie tinut cont de sfera in vârful căreia se extrag datele experimentale, putând extrapolând rezultatele către toate sferele din lungul curgerii.
Abaterea medie patratica
Sistemul de măsură PIV înregistrează câmpuri de viteză instantanee într-un plan paralel la planul median al conductei. Pornind de la aceste câmpuri de viteză instantanee se pot deduce câmpuri medii ale curgerii și apoi se calculează parametrii statistici ai curgerii turbulente utilizând descompunerea Reynolds descrisă în relația:
În care: este viteza instantanee, este viteza medie iar reprezintă fluctuațiile vitezei.
Se notează abaterile medii pătratice pentru fluctuațiile vitezei în sensul curgerii :
și pentru fluctuațiile de viteză perpendiculare pe sensul curgerii.
În, și …… este prezentată variația abaterii medii pătratice și a covarianței în lungimea curgerii, la diferite adâncimi ale curgerii și pentru diverse configurații ale frontierei studiate. Datele obținute sunt comparate cu date similare din literatură, Agelinchaab și Tachie (2006) si E. Iatan (2013).
Se notează cu „” abaterile medii pătratice normalizate cu viteza de frecare:
,
Configuratia este asemanatoare cu ce utilizata de Agelinchaab și Tachie (2006), dupa cum se poate vedea si din figura urmatoare.
In figura este prezentat planul de masura pentru care s-au prelucrat datele, in timp ce in figurile… sunt prezentate variatiile abaterii medii patratice si ale covariantei conform studiului Agelinchaab și Tachie (2006), respectiv valorile obtinute pentru aceiasi parametri in studiul E. Iatan (2013).
Dupa cum se poate vedea si in graficele de mai jos, valorile obtinute atat pentru abaterea medie patratica, cat si pentru covarianta sunt apropiate de cele din literatura.
In cazul studiului Agelinchaab și Tachie (2006), valoarea maxima pentru abaterea medie patratica (u) se obtine pentru raportul y/k = 1, si anume 0,047 m/s. In cazul nostru, pentru masurarile la debitul de 0,7 l/s, valoarea maxima se obtine pentru raportul y/k = 1, 5, cu o valoare maxima de 0,045 m/s. Spre comparatie, datele obtinute in cadrul studiului E. Iatan (2013), indica o valoare maxima pentru u de 0,087 m/s, pentru raportul p/k = 0,75.
Se observa ca, in cazul masurarilor la debitele de 0,9 l/s si 1,1 l/s tendinta se pastreaza, valorile maxime ale abaterii medii patratice obtinandu-se pentru raportul p/k = 1,25.
Comparand datele obtinute cu cele din literatura, se observa valori apropiate studiul Agelinchaab și Tachie (2006), la valori similare ale vitezei de curgere. In cazul studiului E. Iatan (2013), valorile sunt mai mari. Desi modelele de laborator sunt similare, diferenta dintre valorile obtinute se poate datora simularii curgerii gravitationale. In cazul prezentului studiu, curgerea gravitationala este controloata prin montarea unui camin de linistire. Acesta este instalat petru a ne asigura de o curgere gravitationala, cu suprafata libera, apa evacuandu-se prin conducta de plexiglas dupa ce nivelul apei din acest rezervor depaseste nivelul inferior al conductei. Debitul se regleaza cu ajutorul fluxostatului, pe conducta de alimentare cu apa a rezervorului.
Si in privinta covariantei, se observa valoarea maxima obtinuta pentru configuratia p/k = 1.
Pentru validarea modelului experimental, s-au comparat valorile abaterii medii patratice si ale covariantei cu date din literatura (E. Iatan, 2013) pentru planul 1 (planul median).
In figura….. etse prezentat planul de masura pentru care s-au prelucrat datele.
Dupa cum se poate observa, din datle obtinute pentru debitul de 0,7 l/s, pentru valoarea raportului p/k =0,75, distributia valorilor lui u este una parabolica, cu o valoare maxima de 0,13 m/s, similar cu datele obtinute in cadrul studiului E. Iatan (2013). Pentru valori ale lui p/k > 0,75, distributia valorilor lui u este aproape liniara, tendinta ce se observa si in datele din literatura.
In privinta covarintei, distributia valorilor obtinute la raportul p/k = 0,75 este tot una parabolica, cu o valoare maxima de 0,0012 m2/s2. Pentru celelalte valori ale lui p/k, distributia lui uv este una aproape liniara.
Aceleasi tendite se pastreaza si pentru celelalte masurari la debitele de 0,9 l/s si 1,1 l/s.
Calculul vitezei de frecare utilizand metoda profilului logaritmic de viteze
Zona din curgere în care distribuția vitezei este considerată logaritmică diferă de la o situație la alta, de exemplu: Bridge și Jarvis, (1982), pentru curgeri în râuri au considerat valori între 15 și 20 % din adâncimea curgerii, Ferguson și Ashworth, (1992), au propus valori de peste 50 % din adâncimea curgerii pentru distribuția logaritmică. Ne interesează pentru cazul de față, în care frontiera curgerii prezintă rugozități, realizarea unui studiu de sensibilitate în funcție de poziția primului punct introdus în profilul logaritmic față de frontiera inferioară a curgerii. Un studiu mai complex privind sensibilitatea metodei ce utilizează distribuția logaritmică a vitezei este prezentat de Biron (1997).
Pentru calculul vitezei de frecare se va utiliza distribuția logaritmică a profilului de viteze, metoda de calcul fiind prezentată în Capitolul ……
In cele ce urmeaza, este prezentata determinarea vitezei de frecare prin metoda profilului logaritmic de viteze, pentru cele doua planuri de masura, la toate cele trei valori ale debitelor debite 0,7 l/s, 0,9 l/s si 1,1 l/s.
In figurile……, viteza de frecare este determinata prin metoda profilului logaritmic de viteze pentru planul situat pe varful unei sfere.
Aceeasi metoda s-a folosit pentru determinarea vitezei de frecare in planul median, conform figurii….., prezentata mai jos.
Pentru setul al treilea de masurari, s-a inlocuit placa cu sfere, cu o placa lisa (figura…).
In mod similar, s-a determinat viteza de frecare prin metoda profilului logaritmic de viteze.
Valorile obtinute pentru viteza de frecare in cazul masurarilor pentru cele trei debite, pe planul de masura prezenta in figura…., se regasesc in tabelul…., de mai jos:
Calculul vitezei de frecare utilizand metoda extrapolarii efortului Reynolds
Pentru utilizarea acestei metode sunt necesare măsurători pe toată adâncimea curgerii astfel încât să identificăm corect zona liniară din distribuția efortului turbulent. În această lucrare am dispus de măsurători pe toată adâncimea curgerii doar pentru planul 1 la debitul de 0,7 l/s.
Calculul vitezei de frecare s-a realizat pentru planuri diferite si anume pe sfera, pe marginea a doua sfere succesive si in spatiul median dintre acestea, acest lucru fiind detaliat in figurile….
Valorile vitezei de frecare obtinute prin metoda extrapolarii efortului Reynolds pentru cele patru planuri de prelucrare a datelor, sunt prezentate in tabelul….
Calculul vitezei de frecare utilizand metoda abaterilor medii patratice ale fluctuatiilor vitezei
In figurile….. este prezentat determinarea vitezei de frecare prin metoda abaterilor patratice medii ale ale flucatuaiilor vitezei. Viteza de frecare s-a obtinut prin incercari succesive, modificand valoarea vitezei de frecare pana cand profilurile experimentale s-au potrivit cu cele teoretice.
Rezultatele obtinute pentru viteza de frecare, utilizand metoda abaterilor medii patratice ale fluctuatiilor vitezei, sunt prezentate in tabelul….
In tabelul…. sunt prezentate valorile vitezei de frecare obtiunte prin cele trei metode, la debitul de 0,7 l/s. Se observa ca toate cele trei valori sunt apropiate ca ordin de marime si sunt similare cu valori obtinute in literatura. Vitezele obtinute prin metoda extrapolarii efortului Reynolds si prin metoda abaterilor medii patratice ale fluctuatiei vitezelor sunt apropiate ca valori. Acest lucru se poate datora si faptului ca metoda logaritmica depinde foarte mult de profilul obtinut, in special de punctul de pornire al acestui profil. Aceasta metoda se aplica pentru valorile vitezei masurate in zona rugazitatilor, si anume in primele 20% din inaltimea curgerii, in timp ce metodele ce se bazeaza pe extrapolarea efortului Reynolds, respectiv pe reprezentarea abaterilor medii patratice ale fluctruatiei vitezei, se aplica pe toata inaltimea curgerii.
Abordarea numerică a studiului
Modelarea de tip CFD implică utilizarea unui sistem de ecuații cu derivate parțiale format din:
Ecuația de continuitate ce exprimă conservarea masei de fluid.
Ecuațiile de mișcare Navier-Stokes, ce exprimă conservarea cantității de mișcare.
Ecuația energiei Fourier-Kirchhoff ce exprimă conservarea energiei.
În cazul curgerilor și fenomenelor de transfer termic din domeniul nostru introducem și următoarele ipoteze simplificatoare: fluidul este considerat newtonian, de obicei monofazic, incompresibil, supus câmpului gravitațional și cu o viscozitate constantă. În aceste condiții, ecuațiile enumerate mai sus pot fi exprimate în felul următor:
Fie componenta vitezei definite într-un reper cartezian pe direcția (i=1,2,3), cu axele x, y, și z. Ecuația de continuitate se poate scrie pentru un volum elementar de fluid în modul următor:
(8)
Iar ecuația de conservare a cantității de mișcare va avea forma :
(9)
Unde : t este timpul, sunt coordonatele carteziene x, y, z (i=1,2,3), este componenta vitezei pe directia , p este presiunea, ρ densitatea, μ viscozitatea dinamică, iar termenul ține cont de eventuale surse.
Ecuația de conservare a energiei pentru un volum elementar de fluid se scrie în modul următor:
(10)
Unde este coeficientul de difuzie, este numărul Schmidt pentru fluid, μ viscozitatea dinamică, Cp căldura specifică, λ conductivitatea termică, T temperatura și un termen sursă.
Modelarea de tip CFD se poate realiza cu diferite grade de finețe și aproximare a variațiilor temporale și spațiale ale parametrilor fluidului, astfel putem distinge între modele DNS (Direct Numerical Simulation), LES (Large Eddy Simmulation) și modele statistice tip RANS (Reynolds Averaged Navier- Stokes).
Modelarea turbulenței
Curgerile de fluide sunt prezente în jurul nostru fie că este vorba de natură sau de aplicații tehnice. În cadrul acestora turbulența este o caracteristică dominantă a curgerilor. Ea este o proprietate a curgerii, nu a fluidului în cauză. Turbulența nu are o definiție specifică, ci se caracterizează mai degrabă prin proprietățile sale [25, 26].
Atunci când o curgere este turbulentă, mărimile fizice precum viteza și presiunea, variază rapid și aleatoriu, iar temperatura și concentrația sunt caracterizate printr-o difuzivitate crescută.
Multă vreme știința nu a dat un răspuns concret dacă turbulența este sau nu aleatorie. Și dacă nu, care este setul de ecuații care caracterizează această curgere, cum se întâmplă în restul fenomenelor din natură. S-a ajuns apoi la concluzia că numai sistemele neliniare pot caracteriza o mișcare haotică și întâmplătoare [27]. Cu toate acestea, turbulența nu este pe deplin înțeleasă, acest lucru rămânând marea provocare a oamenilor de știință.
Turbulența este o stare a mișcării unui fluid caracterizată de structuri spațiale ce se dezvoltă în timp, denumite vârtejuri. Acestea au diferite mărimi caracteristice, cele mai mari dintre ele fiind de același ordin de mărime cu lungimea caracteristică a curgerii (de exemplu, diametrul unei conducte, înălțimea unei încăperi, diametrul unui difuzor, etc.). Aceste vârtejuri au și o viteză caracteristică bine definită în funcție de scara lor spațială și de viteza curgerii. Când întâlnim un regim turbulent, acesta domină de obicei orice alte fenomene de curgere, rezultând o creștere a disipării energiei, a amestecului și transferului de căldură. Caracteristica de a genera noi vârtejuri din vechile vârtejuri este esențiala pentru turbulență. Acest fenomen se explică datorită instabilității structurilor mari care sub acțiunea forțelor de forfecare se întind, generând vârtejuri mai mici ce preiau la rândul lor energia vârtejurilor din care provin. Aceste vârtejuri mici vor fi supuse la rândul lor aceluiași proces și vor da naștere unor vârtejuri și mai mici, continuând până când este atinsă cea mai mica scară, aceasta corespunzând unui echilibru între forțele de inerție datorate turbulenței și forțele de viscozitate moleculară. Aceasta scară se numește ”micro-scara Kolmogorov”. Astfel, energia cinetică este transportată de la scări mari de curgere către scările cele mai mici posibile, fenomen cunoscut sub termenul de cascadă de energie.
Scările de lungime relevante pentru interacțiunile fizice ce au loc într-o curgere turbulentă sunt: L – macro-scara (lungimea caracteristică a fenomenului studiat); lT – macro-scara Taylor (dimensiunea caracteristică a celor mai mari și mai energetice structuri), lλ – micro-scara Taylor (dimensiunea medie a vârtejurilor), – micro-scara Kolmogorov (dimensiunea celor mai mici vârtejuri). Aceasta din urmă este definită în felul următor :
(11)
În care ε este rata de disipare turbulentă, iar ν viscozitatea cinematică a fluidului.
Rata de disipare este la rândul ei definită în felul următor:
(12)
Unde ui' reprezintă fluctuațiile vitezei, iar bara superioară un operator de mediere.
Micro-scara lui Kolmogorov este legată de numărul Reynolds prin relația :
(13)
În care Re este numărul lui Reynolds :
(14)
U fiind viteza caracteristică a curgerii și L lungimea caracteristică a curgerii.
Revenind la tipurile de modele enunțate anterior, fiecare dintre acestea prezintă particularități în ceea ce privește modelarea turbulenței. Astfel, modelele DNS permit descrierea fină a fazei fluide, cu reprezentarea celor mai fine scări temporale și spațiale, prin rezolvarea directă a sistemului de ecuații prezentat în acest paragraf. Pentru ca ecuațiile Navier-Stokes să poată permite obținerea unei soluții cu acuratețe este necesară o discretizare cu celule ale căror dimensiuni să fie de ordinul de mărime al celor mai mici scări temporale și spațiale. Principalul inconvenient este deci legat de faptul că această metodă necesită o discretizare foarte fină ceea ce conduce la un timp de calcul ridicat. Ca și exemplu, numărul necesar de noduri N, din grila de discretizare pentru o direcție spațială poate fi determinat cu relația:
(15)
Unde lK reprezintă dimensiunea celor mai mici structuri turbulente (scara de lungime Kolmogorov)
Pentru curgeri turbulente tridimensionale putem astfel ajunge cu rapiditate la o grilă de discretizare cu N3 noduri. Pe de altă parte, rezoluția temporală a scării de timp trebuie să fie de același ordin de mărime dacă se dorește surprinderea fenomenelor dinamice ale curgerii.
Deocamdată este dificil a folosi astfel de modele pentru curgeri complexe din cauza limitărilor impuse de tehnica de calcul disponibilă. În cadrul domeniului nostru de interes, al curgerilor la scara clădirilor, putem afirma cu certitudine că modelele DNS nu sunt adaptate pentru abordarea numerică a fenomenelor.
În ceea ce privește modelele LES, acestea sunt în plină dezvoltare. Principiul ce stă la baza lor constă în rezolvarea ecuațiilor de conservare ce guvernează curgerea utilizând o discretizare spațială mai puțin fină decât în cazul modelării DNS. Vârtejurile de dimensiuni superioare celor ale celulelor discretizării sunt determinate prin calcul, iar transferul de energie dinspre acestea spre structurile mai fine este reprezentat printr-un model dezvoltat de Smagorinsky [28].
Figura 30 –Extinderea modelării pentru diferite tipuri de modele aplicate la studiul curgerilor turbulente, după [29]
Abordarea statistică (RANS), ce antrenează mai puțin timp de calcul decât cele precedente, este des folosită în codurile de calcul industriale și pare adaptată pentru curgerile specifice la interiorul clădirilor. Obiectivul acestei metode este de a neglija mișcarea turbulentă instantanee a fluidului, mult prea complexă, și de a căuta ecuații care să prezică simplu evoluția câmpurilor medii. Astfel, ecuațiile de mișcare sunt mediate pentru a reduce termenii fluctuanți, iar noile necunoscute sunt luate în calcul în modelele de închidere sau de turbulență [30]. Cu cât aceste modele sunt mai elaborate, cu atât reprezentarea fizică a curgerii este mai reală, dar timpul de calcul se mărește considerabil. Astfel, pentru fiecare tip de curgere, trebuie găsit cel mai bun compromis între precizie și complexitatea modelului utilizat.
Astfel, principiul fundamental al modelarii clasice a curgerilor turbulente se bazează pe descompunerea tuturor variabilelor de curgere într-un termen fluctuant și unul mediu. Pentru determinarea mărimilor medii sunt întâlnite mai multe metode: media spațială, medie temporală, medie statistică sau medie stohastică [31, 32].
Pentru simularea caracterului aleatoriu al turbulenței, se introduce descompunerea Reynolds în ecuațiile de bază. Aceasta abordare presupune descompunerea fiecărui termen scalar u într-un termen mediu ū și un termen fluctuant u’:
U= ū +u’ (16)
Dacă termenul ū este mediat în raport cu timpul, avem:
(17)
Δt este o scară de timp mult mai mare decât cea mai mare scară de timp a fluctuațiilor turbulente. Acestea sunt presupuse a fi aleatoare, deci putem deduce că media unei fluctuații este nulă:
(18)
Astfel pentru o curgere medie staționară rezultă condiția:
(19)
Rescriind ecuațiile de bază, ținând cont de descompunerea Reynolds, obținem:
Conservarea masei:
(20)
Conservarea cantității de mișcare, în care singura forță exterioară care influențează curgerea este cea gravitațională:
(21)
Conservarea energiei:
(22)
În ecuațiile de mai sus întâlnim necunoscute suplimentare (eforturile Reynolds sau turbulente) și (flux de căldură turbulent). Acești termeni reprezintă influența fluctuațiilor diverselor mărimi pentru curgerea medie.
În literatura de specialitate, o abordare statistică poate conduce la rezultate satisfăcătoare în cadrul modelarii curgerii în interiorul clădirilor. În paragrafele ce urmează, vom prezenta principalele modele existente.
Viscozitatea turbulentă
Prin analogie cu legea de comportament care leagă tensorul eforturilor datorate viscozității, Boussinesq propune în 1877 relaționarea eforturilor turbulente cu gradientul de viteză din curgerea medie prin conceptul de viscozitate turbulentă:
(23)
Unde reprezintă viscozitatea turbulentă. Aceasta este proporțională pe de-o parte cu o scară de viteză v și pe de altă parte cu scara lungimilor L, caracterizând curgerea turbulentă:
, Cµ constantă (24)
În raport cu viscozitatea cinematică υ a cărei valoare este dependentă de natura fluidului considerat, viscozitatea turbulentă depinde doar de caracteristicile turbulente locale ale curgerii.
Prin analogie cu transportul turbulent al fluctuațiilor vitezei, termenul de difuzivitate termică turbulentă at (proprietate locală a curgerii) asigură dependența între fluxul de căldură turbulent și gradientul de temperatură în curgerea medie:
(25)
Ținând cont de ipoteza Reynolds care spune că fluctuația temperaturii T’ se comportă la fel ca și în cazul vitezei,, difuzivitatea termică turbulentă poate fi scrisă în funcție de viscozitatea turbulentă și numărul Prandtl turbulent:
(26)
Astfel, ecuația (27) se scrie:
(27)
Pentru a "închide" modelul de turbulență, mai trebuie determinată o singură necunoscută, în acest caz viscozitatea turbulentă. Există mai multe modele de închidere a turbulenței:
Modele de viscozitate turbulentă bazate pe ipoteza Boussinesq;
Modele de transport pentru eforturile Reynolds;
În funcție de ordinul corelațiilor introduse pentru rezolvarea sistemului nedeterminat de ecuații putem avea:
Modele de ordinul I care calculează eforturile Reynolds în funcție de curgerea medie cu sau fără ecuații suplimentare;
Modele de ordinul II care tratează anizotrop tensorul Reynolds (ecuații suplimentare care iau în considerare toate eforturile turbulente precum și fluxurile de căldură turbulente);
Astfel, putem întâlni modele de ordin I fără ecuație de transport, în care deducem câmpuri medii în cazuri simple, modele de ordin I cu ecuație/ecuații de transport, în care deducem câmpuri medii mai complexe și mărimi turbulente caracteristice, modele de ordinul II cu modele de închidere, în care deducem câmpuri medii și câmpuri fluctuante mediate.
În paragrafele următoare vom trece în revistă și vom detalia caracteristicile diferitelor tipuri de modele de turbulență propuse în cadrul codului comercial Fluent și utilizate pe parcursul studiului numeric. Nu vom detalia celelalte modele existente și nici cele propuse de Fluent și care nu au fost utilizate în acest studiu. Dintre acestea amintim modelul de ordinul II RSM ce nu a putut fi implementat în cadrul lucrărilor de teză din cauza limitelor impuse de resursele de calcul.
Modele de ordinul I
Modele fără ecuație de transport
Modelele fără ecuații suplimentare sunt integral bazate pe ipotezele Boussinesq în ceea ce privește conceptul de viscozitate turbulentă. Astfel, viscozitatea turbulentă este calculată cu ajutorul unei relații algebrice.
(28)
Prandtl a propus în 1925 modelul fără ecuație de transport, urmărind un raționament inspirat de teoria cinetică a gazului, model aplicabil curgerilor pentru curgerile cu gradient semnificativ.
(29)
Unde lm reprezintă lungimea de amestec caracterizată de scara de turbulentă în punctul ales.
Aceasta lungime este determinată prin intermediul relațiilor empirice. Pentru o curgere liberă, de exemplu un jet, se presupune că lm este proporțională cu dimensiunea transversală a jetului. Prandtl propune pentru lungimea de amestec în apropierea unui perete, ca și pentru curgerea pe placă plană sau o curgerea în conductă : lm=ky, unde k este o constantă de proporționalitate (constanta von Karaman, de valoare 0.41), iar y este distanța normală la perete.
Utilizarea acestei metode, în ciuda ușurinței de a o aplica, depinde foarte mult de alegerea parametrului lm. Un astfel de model nu poate fi corect utilizat pentru curgeri complexe în care determinarea valorii exacte a lui lm este dificilă.
Astfel, această abordare algebrică nu este adaptată în domeniul clădirilor, pentru că la interior, fenomenele de transport turbulent sunt complexe, cu precădere în camerele în care întâlnim zone de recirculare.
Modele cu o ecuație de transport
Modelele cu o ecuație de transport, întotdeauna bazate pe ipoteza lui Boussinesq, nu mai iau în calcul viscozitatea turbulentă în funcție de parametrii exteriori, temporali și spațiali, pentru că acest termen este determinat de data aceasta cu ajutorul a două necunoscute suplimentare (energie cinetică turbulentă și rata de disipare a acesteia).
Astfel, prima necunoscută, energia cinetică turbulentă k este determinată cu ajutorul unei ecuații de transport în timp ce rata de disipare a energiei cinetice este calculată algebric:
(30)
Ținând cont de definiția lui k, care dă scara de viteza v și de definiția viscozității turbulente din relația (30) avem:
(31)
Unde Cµ este o constantă și L scara lungimilor.
Ecuația de transport pentru energia cinetică turbulentă poate fi scrisă astfel:
(32)
Unde :
transportul energiei cinetice turbulente datorat curgerii medii;
ține cont de redistribuirea în spațiu a energiei cinetice turbulente datorată activității moleculare (difuzie);
= producerea de energie cinetică turbulentă legată de gradientul de viteză medie;
=> transfer prin corelația presiune-viteză; reprezintă difuzia turbulentă pentru k și este modelată cu ajutorul numărului Prandtl relativ la k;
= producere/disipare prin forțele gravitaționale;
rata de disipare vâscoasă a energiei cinetice turbulente, ε, dinspre vârtejurile mari spre cele mici;
Rata de disipare vâscoasă a energiei cinetice turbulente poate fi obținută și cu ajutorul analizei dimensionale:
(33),
Unde este o constantă determinată empiric și L scara lungimilor ce depinde de curgerea studiată.
Întâlnim rezultate bune pentru configurații simple, ca acelea pentru straturi limită sau zone de recirculare, curgeri în conducte, dar nu sunt adaptate pentru calculul curgerilor complexe, cum ar fi un jet tridimensional sau o curgere în jurul unui obstacol.
Modelul Spalart – Allmaras
Majoritatea modelelor cu o ecuație de transport iau în calcul ecuația de transport a energiei cinetice turbulente k. Spalart și Allmaras [33] au propus să se rezolve direct o ecuație de transport pentru viscozitatea turbulentă (modelul Spalart-Allmaras). Față de restul modelelor cu o ecuație, modelul S-A găsește soluții independente de soluțiile calculate pentru celulele alăturate și este compatibil cu orice structură a discretizării. Literatura de specialitate arată că acest model, printre puținele modele cu o ecuație aplicabile în domeniul curgerilor de aer din clădiri, este fiabil și destul de larg utilizat. În 2006, Torano et al. [34] au simulat curgerea aerului în tunel aerodinamic cu o viscozitate turbulentă constantă, comparând modelul k-epsilon cu cel S-A, obținând rezultate performante pentru ambele cazuri. Mai mult, modelul S-A a fost implementat într-una dintre cele mai noi metode de modelare a turbulenței – detached eddy simulation (DES).
Modelul S-A este un model cu o ecuație relativ simplu, ce rezolvă ecuația de transport pentru viscozitatea cinematică turbulentă. Astfel nu mai este necesar calculul scării de lungime pentru grosimea stratului de forfecare. În solverul Fluent, modelul S-A, față de forma sa originală, a fost modificat pentru implementarea funcțiilor de perete, atunci când discretizarea în acea zona nu este suficient de fină. Mai mult, în zonele limită, gradienții variabilelor din ecuațiile de transport sunt mult mai mici față de cei din modele k-epsilon sau k-omega. Acest lucru poate duce la o perturbare mult mai mică a calculului pentru zonele limită unde celulele nu sunt stratificate. Totuși, pentru curgerile complexe, modelul nu este atât de potrivit pentru calcule.
Modelul propus de Spalart și Allmaras rezolvă ecuația de transport pentru o variabilă care este o formă modificată a viscozității cinetice turbulente.
Astfel, , este identică viscozității cinematice turbulente, cu excepție în zonele limită (afectate de viscozitate). Ecuația de transport este:
(34)
Unde:
reprezintă producerea de viscozitate turbulentă;
disiparea viscozității turbulente care apare în zonele limită din cauza obstacolului și amortizării vâscoase;
, sunt constante, iar este viscozitatea cinematică;
Viscozitatea turbulentă se determină astfel:
, , ( 1)
, , (35)
, , (36)
Constantele modelului:
,
Modele cu două ecuații de transport
Modelele cu două ecuații de transport țin cont de o ecuație de transport a energiei turbulente, k, și o ecuație de transport pentru scara de lungime turbulentă sau un alt parametru legat de această mărime. Odată ce ecuația lui k este rezolvată, orice variabilă de forma z=ka/Lb poate fi utilizată pentru că este cunoscut k.
Variabila cel mai des folosită este rata de disipare a energiei turbulente, ε, care apare explicit în termenul sursă al ecuației de conservare a energiei cinetice turbulente k.
Majoritatea calculelor CFD pentru curgerile de aer din interiorul clădirilor sunt bazate pe modelul de turbulență de tipul k-epsilon.
O gamă foarte largă a variantelor acestui model este prezentată cu caracteristicile și sugestiile de îmbunătățire de către David C. Wilcox [35].
Multe alte variante ale modelului k-epsilon au fost dezvoltate ulterior, cum ar fi modelul RNG k-epsilon [36] bazat pe teoria grupului de normalizare, sau modelul k-epsilon „realizable” bazat pe o noua ecuație a lui ε, [37].
În continuare vom prezenta modelul standard k-epsilon, acesta fiind o bună aproximare matematică a curgerilor de aer întâlnire în interiorul clădirilor și este constituit din două ecuații de transport.
Avantajul unui astfel de model este că rezultatul este un bun compromis între calitatea datelor obținute și puterea de calcul implicată.
Modelul k-epsilon standard
Acest model de turbulență cu două ecuații nu ține cont de ipoteza lungimii de amestec. Pentru curgerile cu număr Reynolds mare, viscozitatea turbulentă se scrie de forma:
(37)
Unde rata de disipare a energiei cinetice turbulente este de forma:
(38)
Pentru calculul viscozității turbulente în orice punct al curgerii, la ecuațiile de conservare a masei și cantității de mișcare, trebuie adăugate două ecuații de transport suplimentare:
O ecuație pentru energia cinetică turbulentă:
(39)
O ecuație pentru rata de disipare a energiei cinetice turbulente, obținută prin analogie cu ecuația anterioară:
(40)
Unde Gk este rata de producere a energiei turbulente rezultate din interacțiunea eforturilor turbulente cu mișcarea medie și Gb corespunde producerii de energie datorată forțelor arhimedice. Astfel, avem:
(41)
Constantele utilizate în model sunt determinate empiric, valorile uzuale fiind:
=0.09; =1.44; =1.92; =1; =1.217;=1
Termenul reprezintă impactul forțelor de presiune asupra ratei de disipare a energiei cinetice turbulente [38].
Modelul k-epsilon este superior modelelor cu lungime de amestec, fiind unul dintre cele mai folosite modele în codurile de calcul comerciale. Totuși, acest model prezintă și câteva inconveniente:
Ecuația ratei de disipare este de formă aproximată, obținută prin metode intuitive;
Coeficienții de închidere sunt ajustați într-o manieră empirică pentru reprezentarea fizică a curgerii;
Modelul nu aduce decât informații globale asupra mecanismelor de transfer între diferite scări de turbulență.
Acest model simplu, dar robust, reprezintă un bun compromis între fiabilitate, complexitate și performanțe. Acest fapt explică de ce acest model a devenit o alegere uzuală în simulările numerice ale mișcării turbulente [32, 39].
Modelul k-epsilon reprezintă nivelul minim de modelare fizică acceptat. Acest model este adaptat pentru curgerile cu eforturi tangențiale mari, care se dezvoltă prin intermediul unei guri de refulare, în interiorul clădirilor. Prezintă dezavantaje, în cazul modelarii în apropierea stratului limită și în calculul curgerii în jurul obstacolelor, atunci când apar zone de recirculare sau vârtejuri de scară mare. În aceste cazuri, modele de ordin superior pot fi utilizate cu rezultate mai bune, introducând totuși un timp de calcul cu 50-60% mai mare față de modelul k-epsilon.
Modelul k-epsilon realizable este o variantă a modelului de bază, ce pornește de la ipoteza universalității scărilor mici de turbulență, ducând în final la rezultate mai bune într-un număr mare de aplicații și fiind o bună alternativă și pentru curgerile complexe.
Modelul k-epsilon realizable
Denumirea acestui model provine din faptul că modelul satisface anumite restricții matematice legate de fizica curgerii turbulente. Prezintă o performanță sporită pentru studiul curgerilor ce implică zone de recirculare, strat limită. De exemplu, eforturile normale ale tensorului Reynolds, care sunt pozitive prin definiție, pot deveni negative (ne-realiste) atunci când deformarea devine semnificativă, având pentru curgerea medie următoarea expresie:
(42)
Această inegalitate poate fi obținută după prelucrări cu ajutorul ipotezei lui Boussinesq și expresiei viscozității turbulente cinematice.
În aceste condiții, pentru a se asigura pozitivitatea eforturilor normale, constanta Cµ ia valori variabile. Valoarea potrivită pentru Cµ pentru un sub-strat limită inerțial este de 0.09. Invers, valoarea convenabilă a lui Cµ într-o curgere puternic forfecată este de 0.05. Astfel, modelul k-epsilon realizable propune o expresie de calcul pentru valoarea lui Cµ necesară pentru calcul viscozității turbulente:
(43)
Unde termenii ce apar în expresia de mai sunt definiți astfel:
(44)
(45)
(46)
Unde reprezintă tensorul ratei de deformație medii si tensorul vitezelor medii de deformație exprimat prin relația:
(47)
Celelalte mărimi care apar în ecuații, ε și ω, reprezintă rata de disipare a energiei cinetice turbulente, respectiv rata de disipare specifică sau frecvența specifică. Aceasta poate fi definită prin:
(48)
Constantele AT și AS sunt definite astfel:
și (49)
, , (50)
Plecând de la ecuația (46) ajungem la concluzia că depinde de mărimile turbulente k și ε și de ratele de deformare.
Analizând modelul k-epsilon standard s-a ajuns la concluzia că modelând ecuația ratei de disiparea a energiei cinetice turbulente k este cauza principală a anomaliei comportamentului jeturilor, pe când modelul k-epsilon realizable propune o nouă expresie pentru aceasta ecuație.
Ecuațiile de transport ale lui k și ε în model se scriu astfel:
Ecuația conservării energiei cinetice turbulente:
(51)
Ecuația conservării ratei de disipare a energiei cinetice:
(52)
Dacă în modelul standard, ecuația lui ε se bazează pe termeni de producție și disipare proporționali celor din ecuația lui k, în modelul realizable, modelarea lui ε este modificată prin introducerea fluctuațiilor de frecvență medie, [37], în ipoteza că la numere Reynolds mari se verifică relația următoare:
( 2)
Termenul care desemnează generarea ratei de disipare a energiei cinetice turbulente nu mai este legat de generarea energiei cinetice turbulente k. De altfel, termenul Gk conține o altă expresie în raport cu alte modele de tip k-epsilon.
Parametrul C1 este determinat după expresia:
în care (53)
Unde S reprezintă modulul tensorului ratei de eforturi medii:
(54)
Constanta C3ε care cuantifică influența forțelor arhimedice asupra lui ε, este calculată cu expresia:
(55)
Unde W reprezintă componenta vitezei paralelele la vectorul acceleratei gravitaționale, iar U este componenta perpendiculară pe acest vector.
Celelalte constante au următoarele valori:
=1.9; =1.44; =1.2; =1
Modelul k-epsilon RNG
Modelul k-epsilon RNG a fost creat utilizând metode de tipul Re-Normalisation Group (RNG) [36] pentru a normaliza ecuațiile Navier-Stokes, în scopul de a contoriza și efectul celor mai mici scări de mișcare. În modelul standard k-epsilon viscozitatea turbulentă este determinată cu o singură scară de lungime turbulentă, astfel că difuzia turbulentă este cea care apare doar la scara de lungime specificată, pe când în realitate toate scările de mișcare contribuie la această difuzie turbulentă. Abordarea RNG, o tehnică matematică menită să determine un model similar lui k-epsilon, redă în final o ecuație a lui epsilon modificată pentru termenul de generare de energie.
Similar ca forma cu modelul standard k-epsilon, dar include:
Un termen adițional în ecuația lui care îmbunătățește modelarea jeturilor puternic deformate;
Acuratețe pentru curgerile turbionare;
O formulă analitică pentru numărul Prandtl turbulent, față de valorile constante din modelul standard;
În timp ce modelul standard k-epsilon este un model cu numere mari Reynolds, teoria RNG propune o formulă dedusă pe cale analitică pentru viscozitatea ce apare în urma efectelor numerelor Reynolds mici. Totuși această caracteristică depinde de o abordare adecvată în straturile limită;
Nu redă corect curgerea unui jet circular;
Ecuația conservării energiei cinetice turbulente:
(56)
Ecuația conservării ratei de disipare a energiei cinetice:
(57)
Termenul R din ecuație este un termen adițional legat de deformarea medie și mărimile turbulente.
Constantele sunt determinate utilizând teoria RNG.
Modelul k-omega standard
Modelul k-omega din Fluent se bazează pe modelul k-omega formulat de Wilcox, 1998 (Turbulence Modeling for CFD), în care sunt incorporate modificările pentru efectele produse de numerele Reynolds mici și alte dificultăți de calcul.
Modelul standard k-omega este un model empiric bazat pe ecuațiile de transport ale energiei cinetice turbulente, k, și ratei specifice de disipare, ω, care determină scara de turbulență. De-a lungul anilor modelul a fost modificat, astfel încât termenul de producere a fost adăugat în ambele ecuații, k și ω, fiind îmbunătățită acuratețea în ceea ce privește modelarea curgerilor libere.
Energia cinetică turbulentă, k și rata specifică de disipare, ω, se obțin din următoarele ecuații de transport:
(58)
(59)
În aceste ecuații, termenii sunt explicitați astfel:
Gk reprezintă generarea energiei cinetice turbulente, k, datorată gradienților medii de viteză;
Gω reprezintă generarea ratei specifice de disipare, ω;
Γk și Γω reprezintă difuzivitatea efectiva a lui k și ω, respectiv;
Yk și Yω reprezintă disiparea lui k și ω din cauza turbulenței;
Sk și Sω sunt termeni sursă definiți de utilizator;
Difuzivitatea efectivă pentru modelul k – omega se reprezintă prin relațiile:
(60)
(61)
Unde sunt numerele Prandtl turbulente pentru k și ω. Viscozitatea turbulentă este calculată astfel:
(62)
Corecția pentru numerele Reynolds mici
Coeficientul diminuează viscozitatea turbulentă, fiind nevoie de o corecție pentru numerele Reynolds mici:
(63)
Unde:
(64) ; (65) ; (66) ; (67)
Pentru numerele Reynolds mari,
. (68)
Modelarea generării turbulenței
Generarea lui k
Termenul reprezină generarea de energie cinetică turbulentă. Din ecuația de transport a lui k, acesta poate fi definit:
(69)
Ținând cont de ipoteza lui Boussinesq, avem:
(70)
Unde S este modulul tensorului mediu al ratei de deformare, definit ca și pentru modelul k-epsilon:
(71)
Generarea lui ω
(72)
Coeficientul α este redat de relația:
(73)
Unde Rω= 2.95; α* și Ret sunt definite mai sus în ecuațiile (64) și (65), respectiv.
Pentru numere Reynolds mari, în modelul k-omega, α=α∞=1.
Modelarea disipării turbulenței
Disiparea lui k
(74)
Unde:
(75)
(76)
(77)
(78)
Disiparea lui ω:
(79)
Unde
(80)
(81)
(82)
(83)
Constantele modelului:
; ; ; ;
; ; ;
Modelul k-omega SST
Modelele k-omega sunt fără îndoială o bună abordare pentru curgerile de aer aflate la interiorul clădirilor, prezentând o bună acuratețe și stabilitate numerică. Multe studii de profil indică modelul k-omega SST ca având o performanță superioară modelului k-epsilon standard sau RNG.
Modelul k-ω SST (Shear-Stress Transport) a fost dezvoltat de către Menter [40], pentru a obține acuratețea modelului standard k-ω în zonele parietale și independența curgerii în zonele îndepărtate redată de modelul k-ε. Pentru a obține acest lucru, modelul k-ε a fost transformat într-o formulare de tip k-ω. Modelul k-ω SST este similar celui standard, dar include următoarele modificări:
În cadrul acestui model există o funcție care în zona parietală are o formă (modelul k-omega) și în zona liberă are o altă formă (modelul k-epsilon);
Modelul SST incorporează un termen derivativ de difuzie în ecuația lui ω;
Definiția viscozității turbulente este modificată astfel încât să ia în considerare transportul eforturilor tangențiale;
Constantele de modelare sunt diferite;
Aceste caracteristici dau modelului k-omega SST o mai mare acuratețe și aplicabilitate pentru o categorie mai largă de curgeri de fluid.
Ecuațiile de transport din modelul k-ω SST
Modelul k-ω SST este similar celui standard în ceea ce privește ecuațiile de transport:
(84)
În aceste ecuații întâlnim următorii termeni:
reprezintă generarea de energie cinetică turbulentă datorată gradienților medii de viteză;
reprezintă generarea lui ω;
și reprezintă difuzivitatea pentru k și ω;
și reprezintă disiparea lui k și ω datorată turbulenței ;
reprezintă difuzia încrucișată ;
și sunt termeni definiți de către utilizator.
Modelarea difuzivității efective
Difuzivitatea efectiva pentru modelul k-omega SST este calculată astfel:
Unde sunt numerele Prandtl turbulente pentru k și ω. Viscozitatea turbulentă este calculată astfel:
(85)
Unde S este mărimea ratei eforturilor și :
(86)
Funcțiile de amestec, F1 și F2 sunt :
(87)
(88)
(89)
(90)
(91)
Unde y este distanța până la suprafața vecină, iar este termenul pozitiv al difuziei transversale.
Modelarea generării turbulenței
Generarea lui k
Termenul Gk reprezintă generarea energiei cinetice turbulente definită astfel:
(92)
Unde Gk este definit la fel ca în modelul k-omega standard.
Generarea lui ω
Termenul Gω reprezintă generarea lui ω. Față de modelul standard, această formulare diferă prin faptul că în modelul de față α este redat printr-o formulă și nu are valoare constantă:
(93)
Unde :
(94)
(95)
În care k are valoarea de 0.41.
Modelarea disipării turbulenței
Disiparea lui k
Termenul Yk reprezintă disiparea energiei cinetice turbulente și este definit de o manieră similară ca în modelul standard k-omega. Diferența constă în faptul că termenul fβ* este în modelul SST o constantă egală cu 1. Deci:
(96)
Disiparea lui ω
Termenul Yω reprezintă disiparea lui ω și este definit similar ca în modelul standard, diferența fiind în modalitatea de definire a termenilor βi și fβ : fβ este constant egal cu 1, iar βi este variabil. Disiparea lui ω se definește astfel:
(97)
În loc să fie o valoare constantă, βi este dat de relația:
(98)
Termenul de difuzie transversală
Modelul SST k-omega are la bază deopotrivă componente din modelul k-epsilon și k-omega standard. Pentru a realiza acest model, modelul standard k-epsilon a fost modificat, utilizându-se ecuații pentru k și ω, ceea ce a dus la apariția termenului de difuzie transversală Dω :
(99)
Constantele modelului
; ; ; ;;
;
Toate celelalte constante care apar au aceleași valori ca pentru modelul k-omega standard.
Modele de ordinul II
Față de modelele de ordinul I, bazate pe conceptul de viscozitate turbulentă, modelele de ordinul II se bazează pe ecuațiile de transport pentru eforturile Reynolds. Pentru curgerile caracterizate de tensiuni de forfecare mari, avantajul modelelor de ordin II nu este întotdeauna vizibil, însă situația poate fi diferită pentru curgeri mai complexe. Aceste modele sunt mereu mai puțin disipative față de modelele cu viscozitate turbulentă (care furnizează rezultate mai confuze față de realitate), caracterizează mai bine aspectele neliniare, instabilitățile, permițând un schimb de energie între mișcările fluctuante și mișcările medii.
Ca modele de ordinul II se remarcă modelul RSM (Reynolds Stress Model) și modelul ASM (Algebraic Stress Model). Programul Fluent propune modelul RSM ce presupune modelarea difuziei turbulente, a corelațiilor presiune – tensiuni și a ratei de disipare a energiei cinetice turbulente, prin rezolvarea eforturilor Reynolds. În modelul ASM se pot face două ipoteze: transportul eforturilor turbulente se presupune a fi proporțional cu energia cinetică turbulentă și efectele convective și difuzive sunt considerate neglijabile, ținându-se cont de aceste aspecte în alegerea modelului.
Aceste modele permit o mai bună aproximare a fizicii curgerii și redau o descriere mai corectă a turbulenței. Totuși, nu sunt foarte utilizate, pentru că necesită mult timp și spațiu de calcul și introduc foarte multe dificultăți numerice în utilizare:
Fiecare nouă ecuație introduce un număr de necunoscute din ce în ce mai mare, pentru care trebuie formulate ipoteze;
În ecuațiile de mișcare medie, rolul eforturilor turbulente apare sub formă de termeni sursă.
Comparativ cu modelele tip k-epsilon sau k-omega, modelul RSM necesită cu15-20% mai multă memorie alocată și cu 50-60% mai mult timp de calcul din cauza numărului mare de ecuații de transport adiționale, motiv pentru care acest model nu a făcut parte din opțiunile luate în calcul în cadrul acestui studiu.
Am realizat în Tabelul 2 o trecere în revistă sintetică a modelelor de turbulență prezentate până acum.
Tabel 2 Comparație între modelele de turbulență RANS
Modelarea stratului limită
Majoritatea modelelor de turbulență au fost elaborate cu ipoteza unei curgeri libere și unor numere Reynolds turbulente mari. Astfel, curgerile nu sunt bine caracterizate în apropierea obstacolelor. Mai mult, acestea au o importanță deosebită asupra curgerii fluidului, determinând de exemplu diminuarea scării de lungime turbulentă, anizotropia semnificativă a turbulenței în această zonă și mai ales apariția unei regiuni în care viscozitatea moleculară este preponderentă. Astfel sunt necesare unele corecții în modelele de turbulență sau introducerea unor noi modele chiar în apropierea frontierelor solide.
Regiunea stratului limită poate fi împărțită în 3 zone, pornind de la variația vitezei în raport cu distanta față de perete, y:
Zona interioară foarte aproape de perete unde eforturile vâscoase predomină (substratul vâscos);
Zona de tranziție sau zona tampon unde efectele moleculare și cele turbulente sunt de același ordin de mărime;
Zona exterioară îndepărtată de perete unde întâlnim stratul de frecare turbulentă constantă;
Modelarea curgerii în zona de perete se poate realiza cu ajutorul a mai multe metode:
Metoda funcțiilor de perete: regiunea în care viscozitatea moleculară este luată în calcul nu este rezolvată, modelele de turbulență pentru numere Reynolds mari sunt cuplate la o formulare globală a stratului limită;
Metoda cu două straturi: regiunea apropiată de perete este rezolvată cu modele specifice, modelele de turbulență pentru numere Reynolds mari sunt cuplate cu modele mai simple;
Metoda cu numere Reynolds mici (LRN- Low Reynolds Number): rezolvarea se face până la perete prin introducerea unor funcții de amortizare în ecuațiile de transport ale lui k și ε;
Metoda funcțiilor de perete
Pentru realizarea unui calcul care să redea curgerea cât mai fidel, este necesară discretizarea foarte fină a domeniului de calcul din regiunea substratului vâscos, loc în care apar gradienți importanți ai mărimilor caracteristice. Astfel sunt antrenate resurse și timp de calcul foarte mari.
Astfel, metoda de față propune evitarea discretizării zonei vâscoase și racordarea prin așa numitele legi de perete a condițiilor la limită la condiții impuse în interiorul curgerii, situate la o graniță imaginară, acolo unde efectele moleculare domină. Mărimea modelată cu ajutorul unei funcții de perete este simulată pe un domeniu de calcul care exclude zona vâscoasă, primul nod nemaifiind cel de pe frontiera propriu-zisă, iar condițiile la limită fiind impuse în acest prim nod.
Frontiera reală de calcul
Frontiera imaginară unde se impun condițiile la limită
Figura 32 –Reprezentare schematică a metodei funcțiilor de perete (după[23])
Funcțiile de perete pentru viteză și temperatură sunt obținute plecând de la ecuațiile simplificate de conservare a cantității de mișcare și a energiei, scrise pentru stratul vâscos de-a lungul unei plăci plane.
Metoda dublu-strat
O altă modalitate de a trata curgerile la nivelul straturilor limită constă în utilizarea unor modele de turbulență diferite pentru zonele cu număr Reynolds mic, respectiv pentru cele cu număr Reynolds ridicat. Un model de turbulență uzual este folosit în zonele neinfluențate de frontierele solide (zone turbulente cu număr Reynolds ridicat), iar în stratul vâscos din apropierea pereților se aplică un model de ordin inferior.
Delimitarea între cele două zone este determinată de criteriul Reynolds turbulent, construit cu distanța y de la nodul de calcul la peretele cel mai apropiat:
(100)
Modelele folosite de regulă în zona vâscoasă se clasifică în:
Modele bazate pe introducerea a două scări de lungime caracteristică: difuzia turbulentă și disiparea vâscoasă, această abordare fiind cel mai des întâlnită în simulările CFD pentru curgeri interioare;
Modele bazate pe corecții aduse scării de lungime în zona stratului limită;
Modele bazate pe corectarea expresiei de calcul a viscozității turbulente conform fenomenelor de turbulență redusă specifice stratului vâscos;
Figura 33 –Modelarea stratului limită în Fluent [43]
Metoda cu numere Reynolds mici
Metoda a fost dezvoltată pornind de la modelele de turbulență cu două ecuații de transport (în general modelul k-epsilon).
Metoda constă în introducerea de funcții de amortizare și de termeni adiționali în ecuațiile de transport pentru energia cinetică turbulentă și rata de disipare a acesteia în zonele de strat limită unde numărul Reynolds turbulent local are valori scăzute. Scopul acestor modificări este de a atenua comportamentul turbulent în aceste regiuni în care predomină efectele vâscoase.
Discretizarea domeniului de calcul
Unul dintre cele mai importante aspecte ale modelarii numerice de tip CFD este reprezentat de discretizarea domeniului de calcul. Aceasta reprezintă de fapt prima etapă în simularea numerică a modelelor fizice bazate pe sisteme de ecuații cu derivate parțiale. Discretizarea domeniului de calcul este importantă din cauza faptului că nerealizarea sa în mod adecvat poate duce fie la rezultate eronate fie la blocarea completă a procesului numeric.
Soluțiile analitice ale modelului cu ecuații cu derivate parțiale au o formă continuă în tot domeniul de calcul. Spre deosebire de acestea, soluțiile numerice sunt date într-o formă discretă. Astfel, domeniul continuu de calcul este înlocuit printr-o mulțime finită de puncte denumită rețea de discretizare.
Există mai multe metode discretizare dar cele mai utilizate sunt: diferențe finite, volume finite, elemente finite și elemente de frontieră. Fiecare tip de metodă conduce la aceeași soluție dacă rețeaua de discretizare este suficient fină, dar fiecare dintre ele este de preferat pentru o anumită categorie de probleme. În aplicațiile inginerești obișnuite, codurile comerciale CFD folosesc pe scară largă metoda volumelor finite.
Implementarea Metodei Volumelor Finite în Fluent
În metoda volumelor finite, punctele rețelei se numesc noduri sau vârfuri. Elementele de bază formate din mai multe noduri unite sau conexe se numesc ochiurile sau celulele rețelei.
Solverul Fluent folosește ca metodă de discretizare metoda volumelor finite.
Pașii de implementare a metodei volumelor finite într-un cod CFD sunt următorii:
Discretizarea geometriei de calcul în volume finite pentru calculul principalelor fenomene de transport: difuzie, convecție și termenii sursă;
Proceduri de discretizare pentru fenomene nestaționare;
Procese iterative pentru cuplarea corectă între toate variabilele curgerii;
Algoritmi de calcul pentru sistemele de ecuații discretizate;
Implementarea condițiilor la limită;
Metoda presupune o tehnică de calcul bazată pe volume de control (CV) pentru a converti o ecuație scalară de transport într-o ecuație algebrică ce poate fi rezolvată numeric. Conservarea a unei variabile Φ într-o curgere, de exemplu o componentă a vitezei, într-un volum de control poate fi exprimat ca un echilibru între diferite procese. Cu alte cuvinte, variația variabilei Φ raportată la timp în volumul de control este egală cu fluxul net al variabilei Φ datorat convecției la care se adaugă fluxul net al variabilei Φ datorat difuziei și termenul sursă. Se pornește de la ecuația conservării scrisă în formă integrală, pentru regim staționar:
(101)
Unde este viteza normală la suprafața S a volumului , iar si fiind termenul de difuzie și termenul sursă pentru variabila.
a) b)
Figura 34 –Stabilirea volumelor de control: a) Metoda Cell-Centred b) Metoda Vertex-Centred
Domeniul de calcul este împărțit într-un număr finit de mici volume de control cu ajutorul unei grile definește limitele volumelor de control (Figura 28a) în comparație cu nodurile computaționale, care sunt stabilite în acest caz în centrul de greutate al celulei de calcul (Metoda Cell-Centred ).
Cu această metodă fluxurile sunt calculate pe o suprafață a celulei prin interpolarea valorilor din nodurile adiacente. Pentru anumite grile de calcul, este posibil să se stabilească întâi nodurile de calcul, pentru ca apoi să se construiască în jurul lor celulele care definesc volumele de control prin trasarea fețelor la jumătatea distanței între două noduri (Figura 28b), astfel încât fluxurile pe o suprafață sunt calculate cu ajutorul nodurilor care se învecinează cu suprafața respectivă (Metoda Vertex-Centred ).
Figura 35 –Volum de control tipic și notația utilizată pentru o grilă carteziană 3D [44]
Prima metodă de calcul are avantajul unui consum de resurse computaționale mai redus și aproximează mai bine valoarea pe întreg volumul de control, aceasta fiind dată de nodul central, dar cea de-a doua metodă calculează mai precis valoarea integrată pe o suprafață atunci când aceasta se află la mijlocul distanței între două noduri.
Primul pas în discretizare este de a integra ecuația conservării pe fiecare volum de control și se vor însuma toate ecuațiile pentru obținerea ecuației de conservare globală.
Totul se traduce într-un set de ecuații liniare care trebuie rezolvate pentru obținerea unei valori a fiecărei variabile pentru fiecare celulă, fiind necesară implementarea unei scheme iterative eficiente de calcul.
Rezolvarea ecuațiilor constă în aproximarea integralelor de suprafață și a integralelor volumice. Determinarea integralelor de suprafață conduce la valori pe suprafața CV-ului, acestea fiind calculate prin interpolare cu ajutorul valorilor aflate în centrul celulei. Metoda volumelor finite implică două niveluri de aproximare a valorilor:
Calculul valorilor variabilelor pe suprafața volumului de control – interpolare;
Calculul integralelor volumice și de suprafață – integrare;
Interpolare ,, upwind”
În mod implicit, solverul Fluent stochează valorile discrete ale scalarului Φ în mijlocul celulei. Valorile pe fiecare suprafață a fiecărei celule sunt necesare pentru calculul termenului convectiv din ecuația (102) și trebuie interpolat din valorile aflate în interiorul celulelor. Această interpolare se realizează cu ajutorul unei scheme tip “upwind”. Acest termen implică faptul că valoarea de pe o față Φe este derivată dintr-o celulă din amonte relativ la direcția vitezei normale, . Solverul Fluent dă posibilitatea de a alege dintre mai multe scheme “upwind” :” first order upwind”, ” second order upwind”, ”power law” și” QUICK”. Utilizarea unor scheme de interpolare de ordin mai mare de II sunt dificil de dezvoltat în curgerile tridimensionale.
(103)
Unde H reprezintă termenii de ordin superior.
Când se dorește o acuratețe de ordin I, valorile de pe suprafețele celulei sunt determinate considerând că valorile din centrul acesteia () reprezintă o valoare medie și sunt valabile pentru toată celula – valorile fiecărei fețe sunt identice cu valorile aflate în celula din amonte, ținând seama de direcția de curgere.
Acesta este un mod de aproximare care satisface condițiile la limită necondiționat, nu dă soluții oscilatorii, însă este difuzivă numeric. Dezvoltarea în serie Taylor este în acest caz utilizată doar până la primul termen din partea dreaptă a ecuației de mai sus.
Termenul eroare de trunchiere este de tip difuziv și se aseamănă cu un flux difuziv:
(104)
Unde reprezintă difuzia falsă, numerică.
Această difuzie falsă pune probleme mai ales în cazurile multidimensionale, unde curgerea este oblică pe grilă, eroarea de trunchiere producând o difuzie pe direcția normală pe curgere și pe direcția acesteia. Apariția unor valori false în variabile vor duce la o soluție eronată, rezolvarea acestei probleme fiind generarea unor grile foarte fine.
În cazul interpolării de ordin II, valorile situate pe fețele celulei sunt calculate utilizând o abordare liniară prin interpolarea dintre două noduri apropiate:
Unde este valoarea variabilei pe fața e, iar sunt valorile în nodurile P, W, E și EE; sunt factori de interpolare.
Din analiza erorilor acestei scheme de interpolare se obține o precizie de ordinul II, fiind atins un nivel de precizie mai bun. Dezavantajul constă în faptul că soluția are nevoie de mai mult timp de convergență. O metodă des întâlnită este aceea ce a începe simulările cu o schemă de ordin I, iar după convergența soluției să se înceapă iterațiile cu o schemă de ordin II, acești pași ducând mai rapid la o soluție convergentă.
Aproximarea integralelor pe suprafețe
Fluxul net prin suprafața limită a unui volum de control se determină prin suma integralelor pe fețele acestuia:
(105)
Unde f reprezintă componenta vectorului convectiv () sau difuziv () în direcția normală la fața volumului de control (Figura 29).
Pentru a fi îndeplinită condiția de conservare a unei cantități, este necesar ca volumele de control să nu se suprapună, fiecare față a unui CV fiind unică și aparținând doar a două CV, care se situează de-o parte și de alta.
Pentru a se calcula corect integrala suprafeței, trebuie să se cunoască integrantul f în orice punct al suprafeței S. Cum această informație nu există pentru că numai valorile nodale din centrul CV sunt calculate, trebuie făcută o aproximare. Această aproximare este de două feluri:
Integrala este aproximată în funcție de valorile variabilei într-una sau mai multe locații pe fața în cauză;
Valorile fețelor celulei sunt aproximate în funcție de valorile nodale.
Aproximarea integralei de volum
Unii termeni din ecuația de transport cer integrarea pe volumul CV. Cea mai simplă aproximare de ordin II este să se înlocuiască integrala de volum cu un produs între o valoare medie și volumul CV. Această valoare medie este stabilită ca fiind valoarea nodală:
(106)
Unde este valoarea în centrul CV, această valoare fiind ușor de determinat având toate variabilele disponibile în nodul P, nefiind necesară interpolarea.
Aproximarea este exactă fie dacă este constant fie dacă are o variație liniară în CV, neîndeplinirea acestor două condiții ducând la o eroare de gradul II. O aproximare de grad mai mare implică valori ale lui în mai multe locații, nu doar central. Aceste valori trebuie obținute prin interpolarea valorilor nodale.
Alegerea domeniului de analiză. Generarea frontierelor
Conectivitatea unei rețele definește forma geometrică a elementelor sale. De exemplu, un triunghi este compus din trei noduri, iar un patrulater sau un tetraedru din patru noduri.
În Eroare! Fără sursă de referință. este redată o porțiune dintr-o rețea de discretizare cu elemente de tip patrulater, în planul xy. Distanțele dintre punctele rețelei în direcția lui x, notate Δx, sau în direcția lui y, notate cu Δy, pot fi sau nu constante. Vom reveni ulterior asupra acestui aspect.
Figura 36 –Rețea discretă de puncte [45]
a)b)
c) d)
Figura 37 –Exemple de diverse tipuri de rețele: a) structurată, plană, cu elemente de tip dreptunghi [45], b) nestructurată, tridimensională cu elemente de tip paralelipiped [46], c) nestructurată, plană cu elemente de tip triunghi [47], d) nestructurată, plană cu elemente de tip patrulater [47]
Elementele rețelei pot avea diferite forme geometrice. De exemplu în cazul problemelor bidimensionale, putem avea elemente de tip triunghi, patrulater sau hexagon, iar în cazul problemelor tridimensionale putem avea element de tip tetraedru, paralelipiped sau prisma hexagonală.
O rețea de discretizare se numește structurată dacă conectivitatea sa este de tip diferență finită, adică distanțele dintre nodurile sale sunt constante după cele două sau trei axe ale unui reper cartezian. În mod contrar, pentru o rețea nestructurată distanțele dintre noduri sunt diferite ceea ce înseamnă că conectivitatea este oarecare. Rețelele structurate sunt formate din elemente de tip dreptunghi în plan și din prisme hexagonale în probleme tridimensionale. Rețelele nestructurate folosesc în mod frecvent triunghiuri în plan și tetraedre în spațiu. Alte combinații de elemente geometrice sunt de asemenea posibile (Figura 32). Rețelele hibride sunt compuse din cele doua tipuri de rețele structurate și nestructurate. Există și rețele multibloc și adaptative.
a)b)
Figura 38 –a) Elemente de tip tetraedru (la stânga), diverse prisme și poliedre, b) rețea cu poliedre [45]
Construcția unei rețele de discretizare trebuie să țină cont de geometria domeniului de calcul. Construirea geometriei domeniului și generarea rețelei este pe departe cea mai consumatoare de timp în raport cu întregul proces CFD. Timpul consumat constă în definirea geometriei și introducerea acestor informații în modulul software care generează rețeaua în mod automat. Întâlnirea unei rețele inadecvate problemei tratate, din cauza prea puținelor puncte sau a distribuției necorespunzătoare a acestora, poate conduce în mod frecvent la reconstituiri multiple ale rețelei pentru problema dată, pentru ca simularea curgerii să fie optimizată. Discretizarea spațială a domeniului trebuie să se obțină fără discontinuități ale spațiilor rețelei și fără introducerea de elemente sau celule cu deformări mari. Scopul este acela de genera o rețea cât mai netedă cu putință corespunzând cât mai bine frontierelor fizice ale problemei.
Dacă în trecut, tehnicile de discretizare foloseau rețele de tip structurat ceea ce limita domeniul de aplicare a codurilor CFD la geometrii relativ simple, la ora actuală, datorită dezvoltării metodelor de generare a rețelelor de discretizare, este posibilă reprezentarea unor domenii cu geometrii din ce în ce mai complexe cu ajutorul elementelor nestructurate. Rețelele nestructurate prezintă avantajul unei capacități de adaptare și de automatizare ridicate.
În domeniul construcțiilor, ca și în alte domenii, geometriile problemelor ce trebuie tratate pot fi destul de complexe, de aceea vom prezenta în continuare diferite metode de generare a rețelelor de discretizare nestructurate.
În general rețelele de discretizare nestructurate sunt compuse din triunghiuri în probleme plane și din tetraedre în probleme tridimensionale. Metodele automate de generare a rețelelor de discretizare nestructurate se bazează pe aceste elemente deoarece acestea permit adaptarea facilă a rețelei de discretizare la geometrii complexe ale domeniului de calcul [23].
Există și rețele nestructurate cu elemente de tip patrulater sau hexaedru, iar în ultimul timp se utilizează tot mai mult rețele tridimensionale cu elemente de tip poliedru [23]. Acestea din urmă prezintă o serie de avantaje în comparație cu rețelele ce folosesc tetraedre. Astfel o celulă de tip poliedru are mai multe celule învecinate decât o celulă de tip tetraedru. Acest lucru se traduce printr-o aproximare mai corectă a gradienților mărimilor vectoriale și scalare calculate. Acest lucru asigură, de asemenea, evitarea formării de direcții preferențiale artificiale a curgerii simulate. În același timp o rețea cu poliedre asigură un număr mai mic de noduri ceea ce implică un timp de calcul de câteva ori mai redus. Unul dintre cele mai importante avantaje ale acestui tip de elemente este legat de evitarea apariției de elemente de tip alungit care poate ridica o serie de probleme numerice.
Chiar dacă avantajele enumerate mai sus sunt extrem de importante, utilizarea acestui tip de elemente în codurile comerciale CFD are o aplicare destul de restrânsă. Acest tip de geometrii este mai degrabă folosit pentru probleme fundamentale de cercetare din mecanica fluidelor. Pentru aplicațiile inginerești se folosesc pe scară largă rețele nestructurate cu tetraedre, de aceea vom prezenta mai departe noțiunile de bază legate de generarea de rețele de discretizare ce conțin astfel de elemente.
Orice metodă de generare automatizată a unei rețele de discretizare cuprinde următorii pași de bază:
definirea frontierelor domeniului de calcul;
specificarea funcției de distribuție a dimensiunilor celulelor rețelei;
generarea rețelei de discretizare interioare respectând discretizarea frontierei;
optimizarea rețelei de discretizare.
Primul pas spre generarea rețelei de discretizare este constituit de definirea precisă din punct de vedere matematic a geometriei domeniului de calcul sau a frontierelor sale. Cea mai fericită situație este aceea în care curbele sau suprafețele frontierelor domeniului sunt disponibile ca o funcție analitică de (x,y) în probleme plane, sau de (x,y,z) în probleme tridimensionale. Acest lucru se întâmplă rareori în practică și de aceea se generează curbe sau suprafețe din coordonate discrete sau din alte surse ce conțin asemenea informații [48].
Următoarea etapă poate fi realizată fie în mod implicit – mărimea unei celule interne depinde de maniera în care este discretizată frontiera, sau într-un mod explicit – dimensiunile elementelor de discretizare sunt controlate local.
Ultima etapă este opțională, dar ea asigură obținerea unor rețele de bună calitate.
Metode de tip frontal
Aproximativ 40 % dintre metodele de discretizare cu rețele nestructurate folosesc metode de tip frontal [23, 49, 50]. Unul dintre cei care au contribuit la dezvoltarea acestor metode este Lohner [51].
Acest tip de discretizare se realizează plecând de la un front inițial, de exemplu un plan, un punct sau o latură a frontierei domeniului de calcul.
Frontul de discretizare se deplasează de la interior prin inserarea de puncte în funcție de punctele deja existente. Noile puncte sunt legate de elementele frontului de discretizare și formează astfel noi celule ale rețelei. Frontul activ este menținut acolo unde sunt formate noile celule.
a)
b)
Figura 39 –a) Etapele generării frontale a unei rețele de discretizare plane [49] b) Rețea de discretizare tridimensională de tip frontal cu diferite surse de plecare : vârful conului, cercul delimitând baza conului, planul de bază al conului [45]
Algoritmul de discretizare este iterativ, un element este selecționat și se determina amplasarea unui nou punct astfel încât dacă rezultă un nou element ce poate fi acceptat din punct de vedere al calității formei sale, punctul este păstrat și inserat în rețea.
In Figura 33a este prezentat un exemplu simplu [49] în două dimensiuni, de discretizare de tip frontal. Pe măsură ce algoritmul de generare a rețelei este aplicat, frontul avansează astfel încât aria domeniului rămasă nediscretizată să poată fi acoperită cu triunghiuri. Pentru fiecare latură poziționată pe frontul activ, algoritmul caută amplasarea optimă a unui vârf și de asemenea verifică nodurile existente care ar putea forma un triunghi de formă optimă cu latura respectivă. Vom reveni ulterior asupra noțiunii de formă optimă. Algoritmul selectează fie un punct nou, fie un punct existent pentru a forma cel mai bun triunghi posibil. Se verifică de asemenea posibilele intersecții între laturi astfel încât să nu apară situația unei suprapuneri de triunghiuri. Frontul este o structură de date dinamică ce se schimbă în mod continuu pe măsură ce generarea de elemente progresează, la fiecare pas frontul activ este reînnoit. În această structură de date, fiecare segment disponibil de a forma un nou element este memorat în timp ce fiecare segment ce a fost deja integrat într-un element este șters din memorie. Generarea rețelei este terminată atunci când frontul este gol. Acest tip de discretizare este ușor de aplicat dar în cazul problemelor tridimensionale pot apărea uneori probleme de convergență.
Metode de descompunere spațială
Sunt metode de discretizare ce se bazează pe o structură ierarhică de tip arborescent. Discretizarea constă în reconstituirea domeniului de calcul prin asamblarea celulelor prealabil definite în mod recursiv pentru a satisface condiții locale de dimensiune. Celulele finale fiind supuse apoi unui proces de triangularizare.
Aceste metode sunt robuste datorită simplității de punere în practică, dar pot ridica probleme de conectivitate în regiunile sensibile cum ar fi colțuri ale domeniului de calcul. Ele nu sunt foarte răspândite în cadrul codurilor comerciale CFD, fiind folosite sub 10% dintre acestea[23].
Metode de tip Delaunay
Aceste metode sunt cele mai populare pentru generarea de rețele de discretizare nestructurate triangulare sau tetraedrale. Ele sunt bazate pe aplicarea criteriului lui Delaunay.
În rețelele de discretizare structurate, conexiunile dintre puncte sunt definite în mod automat având în vedere reperul cartezian considerat, după ordinea (i, j, k). În rețelele de discretizare nestructurate nu există o astfel de ordonare a punctelor. Astfel, conexiunile dintre puncte trebuie și ele definite și memorate pe lângă amplasarea în spațiu a respectivelor puncte. Metodele de triangulație Delaunay utilizează un criteriu deosebit de simplu pentru conectarea punctelor pentru a forma elemente conforme nesuprapuse. Acest tip de construcție geometrică a fost cunoscut de multă vreme dar a fost doar de curând pentru generarea de rețele de discretizare în codurile CFD. Criteriul geometric folosit furnizează doar un mecanism de conectare a punctelor, generarea lor trebuie realizată independent. Generarea unei rețele de discretizare prin metoda Delaunay implică, prin urmare, două probleme distincte: crearea punctelor și conectarea lor. În 1850 Dirichlet a propus o metodă de descompunere a unui domeniu dat într-un spațiu arbitrar, într-un ansamblu de regiuni convexe [52]. Pentru o mulțime de puncte P, domeniul este împărțit în regiuni astfel încât fiecare regiune este mai apropiată de un punct P decât de oricare alt punct. Acest tip de descompunere geometrică este cunoscut ca tesselarea lui Dirichlet. Această operațiune aplicată unui domeniu închis are ca rezultat o mulțime de regiuni convexe denumită diagrama lui Voronoi, sau regiunile lui Voronoi. Acestea reprezintă mulțimea punctelor cele mai apropiate de un punct Pi decât față de orice alt punct [53].
Figura 40 –Diagrama Voronoi și triangulare Delaunay (linii punctate)[49]
Figura 41 –Criteriul lui Delaunay: a) respectarea criteriului, b) nerespectarea criteriului.
Figura 34 ilustrează faptul că în plan, laturile unui poligon Voronoi situate în jurul unui punct P este construit din medianele segmentelor ce unesc punctul P cu toate punctele vecine lui. Dacă toate perechile de puncte ce aparțin a două poligoane adiacente sunt reunite se obține o triangularizare Delaunay. În trei dimensiuni, granițele subdomeniilor sunt reprezentate de fețele poliedrelor Voronoi echidistant amplasate în raport cu perechile de puncte. Triangularizarea are ca rezultat în acest caz obținerea unei mulțimi de tetraedre. O proprietate interesanta a triangulație Delaunay este criteriul Delaunay care constă în verificarea proprietăților următoare: cercul (sau sfera) înscris elementului de triangularizare nu conține nici un vârf al ansamblului de puncte, pe de altă parte cercul circumscris nu conține decât vârfurile elementului respectiv (Figura 35).
Operații de optimizare și adaptare
Chiar dacă metodele de generare a rețelelor de discretizare sunt automatizate, în majoritatea cazurilor, simpla lor aplicare nu poate garanta obținerea unor rețele de o calitate suficientă pentru folosirea lor directă în calculele de simulare. De cele mai multe ori este necesară o intervenție din partea utilizatorului pentru corectarea problemelor apărute, pentru adaptarea și optimizarea rețelei. Calitatea unei rețele nu poate fi exprimată folosind o unică definiție. Fiecare problemă în parte are nevoie de o rețea de discretizare dedicată. În general, apreciem calitatea unei rețele de discretizare în funcție de calitatea rezultatelor obținute. Acest lucru nu este întotdeauna evident dacă nu există o metodă de verificare a acestor rezultate. Este de dorit realizarea unei validări experimentale pentru a putea aprecia calitatea simulării și implicit a discretizării. La un nivel superior, o rețea de bună calitate, asigură nu numai obținerea de rezultate corecte, ci și costuri numerice minime – timp redus de calcul ceea ce este echivalent cu un număr redus de elemente de discretizare fără afectarea preciziei rezultatului.
La nivel local, pentru fiecare element în parte, se apreciază factorul de formă al unei celule, ce constă în devierea de la un volum ideal. Pentru probleme tridimensionale CES sau cell equivolume skewness în engleză se exprimă prin raportul:
(107)
Unde reprezintă volumul unei celule de formă geometrică regulată (de exemplu tetraedru) înscris într-o sferă de aceeași rază ca și elementul pentru care se apreciază calitatea, de volum . Notăm că pentru problemele plane, volumele sunt înlocuite de arii.
Cu cât valoarea indicelui CES de calitate a celulei este mai apropiată de zero, cu atât elementul este mai bun. Cu cât valorile sunt mai apropiate de unitate, cu atât mai defectuos este elementul. În această situație, elementele sunt aproape coplanare ceea ce induce dificultăți de ordin numeric.
Optimizarea unei rețele de discretizare constă în ameliorarea globală a calității elementelor sale, și reprezintă ultima etapă din cadrul procesului de generare. Există două metode principale de realizare a optimizării unei rețele de discretizare: metode cu poziția fixă a vârfurilor pentru care se poate schimba conectivitatea elementelor, și metode cu conectivitate fixă pentru care se pot schimba pozițiile vârfurilor.
Metodele de optimizare cu conectivitate fixă implică repoziționarea iterativă a vârfurilor pentru a ameliora local calitatea elementelor. În general varietatea de tehnici existente reprezintă varianta unei tehnici de bază propuse de Field [54, 55], în care poziția unui punct al rețelei este înlocuită printr-o amplasare medie obținută în raport cu pozițiile vârfurilor care împart același element ca și punctul ce urmează a fi deplasat. Diferitele metode se deosebesc prin procedeul prin care se calculează această nouă poziție medie a punctului. Se pot adăuga și condiții legate de criterii de calitate a elementelor nou create.
Metodele de optimizare cu poziție fixă păstrează amplasarea vârfurilor prin ameliorarea formei sau a topologiei elementelor de discretizare. Astfel pentru probleme plane se poate schimba latura comună a două elemente sau pentru probleme tridimensionale, fața comună a două elemente.
Figura 42 –Schimbarea laturilor în două dimensiuni (după [23])
Metodele cu schimbarea topologiei se bazează pe un concept de grad optim la vârf, înțelegând prin gradul unui vârf numărul de laturi incidente în acel punct. În două dimensiuni, valoarea gradului optim este de 6, în timp ce în trei dimensiuni această valoare este de 12.
O rețea de discretizare chiar și după ce este supusă procedurilor descrise anterior poate necesita alte modificări. Aceste operații finale privesc mai ales densitatea nodurilor (adaptare sau rafinare a rețelei) și urmăresc adaptarea rețelei la particularitățile problemei tratate. Anumite fenomene fizice pot necesita rafinarea locală a rețelei, de aceea trebuie să se determine unde și cum se adaptează rețeaua. Regiunile unde se adaptează rețeaua sunt strâns legate de natura problemei. Este necesară uneori cunoașterea a priori a unor regiuni particulare din curgerea studiată, de exemplu: locurile unde gradienții sunt importanți (în stratul limită), sau unde pot apărea desprinderi (puncte de schimbare a curburii suprafețelor solide) (Figura 37). Aceste regiuni particulare se pot comporta diferit în funcție de regimul curgerii (Figura 38)
a) b)
Figura 43 –a) Rafinarea rețelei în stratul limită, b) rafinarea rețelei în funcție de curbura suprafețelor (după [45])
Figura 44 –Stânga – rețea de discretizare pentru curgerea pe placă plană la număr mare Reynolds, Dreapta – rețea de discretizare pentru curgerea pe placă plană la număr mic Reynolds, (după [45])
Există mai multe metode principale de adaptare și anume: subdivizarea laturilor elementelor de discretizare, subdiviziunea directă a elementelor de discretizare, inserarea de noi puncte în rețeaua de discretizare.
Subdivizarea unei laturi a elementelor duce la decuparea fiecărui element ce conține latura respectivă în alte două elemente de discretizare (în plan aceasta se poate traduce de exemplu prin decuparea fiecăruia dintre cele două triunghiuri ce conțin latura respectivă în alte două triunghiuri, în spațiu prin decuparea fiecărui tetraedru în alte două tetraedre.
Inserarea de puncte noi reprezintă o metodă simplă de rafinare a rețelei de discretizare. Această operație poate consta în introducerea unui punct în centrul de greutate al fiecărui element, ceea ce conduce de exemplu la decuparea unui triunghi în trei și a unui tetraedru în patru. O problemă legată de acest procedeu este afectarea calității a noilor elemente rezultate. În acest caz este necesar de exemplu să se recurgă la o nouă optimizare printr-un algoritm de tip Delaunay pentru ștergerea locală a anumitor elemente și pentru conectarea noului punct respectând criteriul Delaunay.
Deseori în modulul de generare a rețelelor de discretizare din codurile comerciale CFD se folosește importarea directă a unor geometrii create cu ajutorul software-urilor de proiectare CAD ce reprezintă domeniul sau o parte a domeniului de calcul. În acest caz, uneori se poate întâmpla ca o data importata acestei geometrii să nu fie potrivita pentru generarea directă a unor rețele de discretizare. Acest lucru se poate datora de exemplu diferenței de precizie spațială dintre programul de discretizare și programul de proiectare [49]. Mai multe probleme pot apărea astfel încât discretizarea domeniul format din geometria importată să nu poată fi făcută în mod direct:
o față a domeniului poate fi definită de muchii care nu sunt coincidente (Figura 39);
un volum poate fi alcătuit din fețe vecine cu muchii ce nu coincid (Figura 40);
o geometrie poate fi alcătuită din volume cu fețe comune ce nu sunt coincidente;
Aceste imperfecțiuni creează probleme de discretizare și pot duce la blocarea modulului de discretizare. Uneori pot fi importate elemente de discretizare având proprietăți nepotrivite (CES de calitate mică), având de exemplu una dintre laturi foarte mici (Figura 41). Acest lucru poate duce la soluții mai puțin precise și chiar și la divergențe în procesul de soluție [45].
Figura 45 –Suprafețe cu muchii care nu coincid [45]
Figura 46 –Suprafețe învecinate cu muchii care nu coincid [45]
Figura 47 –Element de tip triunghi având un CES neadecvat [45]
Din cauza imperfecțiunilor menționate anterior, o geometrie importată dintr-un soft CAD trebuie verificată și dacă este nevoie, aceste imperfecțiuni trebuie remediate. Astfel elementele cu laturi foarte mici pot fi înlăturate prin unirea vârfurilor sau prin contopirea lor cu un element vecin. Muchiile adiacente pot fi contopite.
Alegerea modelului numeric
Măsurările experimentale sunt scumpe și consumatoare de timp. Acesta este motivul pentru care s-a decis încă de la începutul acestor lucrări de cercetare efectuarea de simulări numerice (mai puțin scumpe și mai ușor de executat) și validarea experimentală a acestora doar pentru diferite modele inițial propus. Așa cum a fost arătat în Referatul al doilea din cadrul programului de pregătire doctorală, alegerea noastră a fost orientată către o abordare de tip CFD.
Studiile efectuate în aceasta teză nu includ o evaluare a programelor comerciale de tip CFD, motivația alegerii programului Ansys FLUENT ca unealtă de modelare fiind dată de larga sa răspândire în momentul de față dar și disponibilitatea sa la UTCB.
În afară de modelul propus în referatul anterior, am continuat în paralel dezvoltarea unui alt model mai simplu, conținând o singură semisferă elementară, pentru validarea unor ipoteze de dinamică a structurilor tipice curgerilor studiate pe baza unor rezultate experimentale din literatură. Acest nou model elementar ne va ajuta să calibrăm anumite modele teoretice și să propunem o analiză fină a fenomenelor dinamice la scara stratului limită pentru o singură semisferă.
Figura 48 – Complexitatea dinamicii instabilităților în strat limită la nivelul unei semisfere [25]
Figura 49 – Configurație experimentală studiată de Dixen et al. [25]
Figura 50 – Dinamica vârtejurilor în jurul unei semisfere obținută din modelare numerică [25]
Fenomenele globale vor fi studiate cu ajutorul modelului cu mai multe semisfere descris anterior, ce va fi validat pe baza masurarilor PIV prezentate în capitolul precedent.
În Figura 31 este prezentat domeniul de calcul în cazul unei singure semisefere elementare și discretizarea inițială, urmând să continuăm cu un studiu de dependență a soluției numerice fată de nivelul de discretizare al grilei și să vedem și pe această cale dacă modelul de turbulență spre care ne-am orientat anterior (raportul 2 de cercetare) este potrivit și la micro-scara unei singure semisfere, pentru analiza dinamicii feneomenelor de instabilitate ale stratului limită. O abordare de tip LES este de asemenea prevăzută.
Figura 51 – Domeniul de calcul pentru cazul cu o singură semisferă elementară
Figura 52 – Dicretizarea inițială cu 1.25mil elemente
Reluăm în cele ce urmează prezentarea modelului cu mai multe semisfere prezentat în raportul de cercetare anterior.
Acuratețea simulărilor CFD depinde de redarea cu finețe a geometriei care definește domeniul. În cazul simulărilor ce presupun suprafețe complexe, crearea și generarea grilei de calcul reprezintă o provocare.
Datorită simetriei geometrice a cazului studiat pentru simularea numerică a fost utilizat ½ din domeniul investigat experimental.
În urma alegerii formei geometrice corespunzătoare simulărilor s-a trecut la pasul următor, alegerea discretizării. Determinarea numărului necesar de celule pentru calculul soluției se realizează în urma unui studiu de dependență a soluției față de numărul de elemente de discretizare. Acest studiu este necesar pentru obținerea unei soluții cât mai corecte, aceasta fiind proporțională cu numărul de celule din grila de calcul. Totuși, odată cu creșterea acestui număr, timpul de lucru respectiv resursele de calcul trebuie să crească corespunzător. De la un anumit număr de celule, pentru numere superioare de celule, diferența între soluțiile găsite, nu mai variază în așa mare măsură astfel încât se poate alege o grilă de calcul care să rezolve compromisul soluție corectă – resurse computaționale.
Figura 53 – Domeniu investigat prin simulare numerică și condițiile la limită
1.6 million elements a)
3 million elements b)
7 million elements c)
Figura 54 – Discretizările inițiale cu elemente tetraedrale
Numărul minim de elemente ale grilei de calcul de la care soluția nu mai variază foarte mult, se determină realizând simulări numerice succesive cu diferite discretizări spațiale. Pentru studiul independenței rezultatelor obținute prin simulare numerică față de discretizarea geometrică au fost realizate trei grile diferite de: 1.6, 3 și 7 milioane elemente tetraedrale.Cele trei discretizări geometrice au fost importate în programul Fluent unde celulele tetraedrale din domeniul investigat au fost convertite în celule poliedrale. Prin acest mod numărul de elemente s-a micșorat numărul de elemente în cazul discretizării de la 1.6 milioane elemente tetraedrale la 0.34 milioane elemente poliedrale, de la 3 milioane elemente tetraedrale la 0.64 milioane elemente poliedrale, iar în cazul discretizării geometrice de 7 milioane elemente tetraedrale, numărul de elemente s-a micșorat la 1.5 milioane elemente poliedrale.
0.34 milioane elemente poliedrale a)
0.64 milioane elemente poliedrale b)
1.5 milioane elemente poliedrale c)
Figura 55 – Nivelul de discretizare pentru cazurile studiate, utilizând elemente poliedrale
Pentru realizarea studiului de dependența a soluției în funcție de discretizarea aleasă, ne-am hotărât să alegem un model de turbulență – cel care ni se părea cel mai adaptat pentru utilizarea sa ulterioară. Timpul limitat pentru realizarea acestui studiu nu ne-a permis realizarea analizei mai multor cupluri model de turbulență – discretizare așa cum este recomandat de către specialiștii CFD [26].
Bazându-ne pe lucrări din literatura recentă, modelul de turbulență k-omega SST se dovedește cel mai fiabil dintre modelele cu două ecuații atunci când este dorită reproducerea unor curgeri relativ complexe, caracterizate de valori ale numărului Reynolds relativ mici [27, 28]. Am ales deci acest model pentru abordarea numerică, și testarea discretizării spațiale a fost realizată pentru acest model de turbulență.Pentru fiecare grilă în parte s-a calculat valoarea factorului de simetrie (skewness) a celulelor. O celulă cu valoarea factorului de simetrie egală cu 0, descrie o celulă structurată regulată, în timp ce o valoare a factorului de simetrie egală cu 1 descrie o celulă de cea mai proastă calitate.
În continuare sunt prezentate câmpuri de viteză și vorticitate în planul median al primului rând de semisfere paralel cu planul median al conductei.
Numărul minim de elemente de discretizare a fost determinat după efectuarea de calcule pentru fiecare dintre cazurile enumerate precedent. Am urmărit evoluția câmpurilor de viteză și vorticitate din figurile anterioare în încercarea de a remarca cazul pentru care nu mai sesizăm o schimbare vizibilă a acestor distribuții. Am ales geometria de 0.63 milioane de celule pentru a fi utilizată mai departe în restul cazurilor studiate, deoarece îndeplinește condițiile necesare unor soluții independente de discretizare utilizată.
1a)1b)
1a) 2b)
3a)3b)4a)4b)
Figura 56 – Distribuțiile valorilor medii ale componentei longitudinale a vitezei (a) și ale vorticității transversale (b): 1) PIV, 2) k-omega SST 0.34 milioane de elemente, 3) k-omega SST 0.64 milioane de elemente, 4) k-omega SST 1.5 milioane de elemente
În literatură se sugerează posibilitatea utilizării modelului de turbulență k-epsilon RNG pentru modelarea curgerii cu suprafață liberă în canal circular cu rugozități artificiale (Gao 2007). Am testat și noi posibilitatea utilizării acestui model prin comparație cu distribuțiile valorilor medii ale componentei longitudinale a vitezei. În figura 104 sunt prezentate aceste distribuții împreună cu datele experimentale. După cum se poate remarca cu ușurință în această figură rezultatele obținute cu ajutorul modelului k-epsilon RNG nu sunt satisfăcătoare în cazul nostru.
1a)1b)
2a)2b)
Figura 57 – Distribuțiile valorilor medii ale componentei longitudinale a vitezei (a) și ale vorticității transversale (b): 1) PIV, 2) k-epsilon RNG -0.64 milioane elemente
Pentru aceeași grilă de discretizare am testat și posibilitatea utilizariisimularii numerice cu vârtejuri de scară mare (LES). În Figura 105 și Figura 106 sunt prezentate distribuțiile valorilor medii, respectiv instantanee ale componentei longitudinale a vitezei și a vorticității transversale. După cum se poate remarca în Figura 10, zonele de vorticitate negativă asociate vârtejurilor normale generate de suprafețele semisferelor sunt surprinse destul de bine de modelul LES.
1a)1b)
Figura 58 – Distribuțiile valorilor medii ale componentei longitudinale a vitezei (a) și ale vorticității transversale (b): 1) PIV, 2) LES
1a)1b)
2a)2b)
Figura 59 – Distribuțiile valorilor instantanee ale componentei longitudinale a vitezei (a) și ale vorticității transversale (b): 1) PIV, 2) LES
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: In partea a patra sunt prezentate concluziile acestei lucrări, contribuțiile personale si direcțiile de cercetare ulterioare in acest domeniu. [301414] (ID: 301414)
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.