Universitatea Politehnica Din Bucuresti

Alexandru Sorin Patulea Prof.dr.ing. Diana Robescu

-2010-

Cuprins:

1.Generalitaţi.

1.1. Curgerea irotaţională a lichidelor în (x,y)1.2. Curgerea printr-un canal îngustat

2.Curgerea viscoasă.

2.1. Determinarea regimurilor de curgere a lichidelor în conducte 2.2. Curgerea viscoasă în canal 2.3 Curgere viscoasă printr-un canal îngustat2.4. Injecţia la viteză uniformă2.5. Similaritate dinamică

1.Generalităţi

Curgerea fluidelor are ca obiect studiul fluidelor sub aspectul comportării lor mecanice. Mai exact, sunt studiate echilibrul (statica) şi mişcarea (dinamica) fluidelor, precum şi interacţiunile dintre acestea şi suprafeţele solide cu care sunt în contact. Este o ramură a mecanicii mediilor continue, domeniu care modelează materia la nivel macroscopic, făcând abstracţie de comportarea la nivel atomic si nuclear. Dinamica fluidelor, constituie un domeniu de cercetare activ cu multe probleme nerezolvate sau rezolvate parţial. Curgerea fluidelor poate fi formulată printr-un formalism matematic avansat bazat pe teoria ecuaţiilor diferenţiale şi algebra complexă. Modelul matematic este obţinut şi prin întrebuinţarea calculului numeric implementabil pe diverse programe CAE sau PDE de simulare. Deasemenea, folosind proprietatea vizibilităţii deosebite a curgerii, fluidele pot fi analizate comportamental prin metoda vizualizării traiectoriilor particulelor.

1.1 Curgerea irotaţională a lichidelor în (x,y)

Cum densitatea normală a unui lichid în domeniul de presiune întâlnit în aplicaţiile practice este putin variabila putem considera ca aceasta este constantă. Mai putem presupune că lichidul alunecă liber pe suprafaţa solidă şi că forţele vâscoase sunt semnificativ mai mici decât forţele ineţiale. Aceste presupuneri sunt valabile în majoritatea cazurilor practice. Deci legea conservării de masă se poate exprima astfel:

∇ ∙ ( ρ0V )=−∂ ρ0

∂t

(1.1)

Unde ρ0 este densitatea masei şi v este vectorul viteză. Presupunând că densitatea este constantă putem deduce că are loc o conservare a volumului după formula:

∇ ∙ V =0 (1.2)

sau în forma explicită:∂ v x

∂ x+

∂ v y

∂ y=0 (1.3)

Cum în FlexPDE programul cu care lucrăm nu putem scrie formula sub această formă apelăm la urmatoarea convenţie. Din fericire putem nota componentele vitezei ca derivate ale funcţiei φ. Funcţiile fiind definite în continuare:

vx=∂∅∂ x

, v y=∂∅∂ y

(1.4)

Conveţie ce conduce mai departe la relaţia:

∂2∅∂2 x

+ ∂2∅∂2 y

=0 (1.5)

Care este o bine cunoscută ecuaţie Laplace.

Până acum sa folosit doar principiul conservării masei dar este important de stiut că orice solutie PDE trebue să fie şi irotaţională (∇ × v = 0), deoarece:

(∇×V )z=∂ vv

∂ x−

∂ vx

∂ y= ∂2∅

∂ x ∂ y− ∂2∅

∂ y ∂ x=0 (1.6)

Dar legea conservării energiei ne conduce mai departe la ecuaţia de miscare a lui Benoulli :

12

ρ0 v2+ p+ρ0 gy=constant (1.7)

v=√v x2+v y

2

Unde:v- este viteza;p- este presiunea g- acceleraţia gravitaţională (considerând axa y veticală)

1.2. Curgerea printr-un canal îngustat

Prima aplicaţie folosind ecuaţiile de la punctul 1 este curgerea printr-un canal orizontal, limitat de o suprafaţă plană perpendiculară pe un domeniu dat.

Următoarea descriere defineşte problema şi introduce în PDE potenţialul vitezei φ. Dupa aflarea lui φ trebuie să derivam pentru a putea afla componentele vitezei. Cunoscând aceste componente putem afla viteza (v). Presupunând o curgere orizontală, unde forţa gravitaţională se anulează, din relatia lui Bernoulli rezultă:

12

ρ0 v2+ p=12

ρ0 v02+p0 (1.8)

Deci în final obţinem expresia presiuni p pe segmentul definit.Pentru optimizarea rezultatului alegem ngrid=1(specifică numărul rândurilor meshări în

fiecare dimensiune, standard ngrid în versiunea student 2d este 10).La definirea Boundary Conditions se specifică viteza de intrare vx0. Pentru iesire se

impune o constantă pentru potentialul φ. Valoarea sa absolută este arbitrară deorece numai derivatele parţiale vor fi folosite, dar specificând valoarea constantei asupra domeniului se intueste vy=dy(φ) şi se forţează ieşirea lichidului pe direcţia x.

În limbajul FlexPDE cele da mai sus se scriu sub forma programului de mai jos:

TITLE 'Curgere printr-un canal îngustat' { fex1.pde }COORDINATES cartesian2 { Sistemul de coordonate, 1D,2D,3D, etc }SELECT errlim=1e-5 ngrid=1 spectral_colors { Culori Curcubeu }VARIABLES { Definitre variabile in PDE } phi { Potentialul vitezei }

DEFINITIONS { Definire parametri } Lx=1 Ly=1 k=0.5 { Coeficientul de comprimare } vx0=3.0 { Viteza de intrare } p0=1e5 {Presiunea atmosferica } dens=1e3 { Densitate} vx=dx(phi) vy=dy(phi) { Componentele vitezei } v=vector( vx,vy) vm=sqrt( vx^2+ vy^2) { Viteza } p=p0+ 1/2*dens*(vx0^2-vm^2) { Presiunea } div_v=dx( vx)+ dy( vy) { Divergenta, or div( v) } curl_z=dx( vy)- dy( vx) { Vorticitate, or curl( v) }EQUATIONS dxx( phi)+ dyy( phi)=0 { Sau div( grad( v)) }BOUNDARIES { Definirea domeniului }region 'domeniu' start 'exterior' (0,Ly)natural( phi)=-vx0 line to (0,-Ly) { intrare }natural( phi)=0 line to (Lx,-Ly) to (2*Lx,-Ly*k) to (3*Lx,-Ly*k)value( phi)=0 line to (3*Lx,Ly*k) { iesire }natural( phi)=0 line to (2*Lx,Ly*k) to (Lx,Ly) to closeMonitorsPLOTScontour( phi) vector( v) norm contour( vm) paintedcontour( p) painted contour( p) zoom(1.5*Lx,0, Lx,Ly)surface( p) zoom(1.5*Lx,0, Lx,Ly)elevation( vm) on 'exterior' { verificarea conditiilor }contour( div_v) contour( curl_z)END

Vectorul este modificat astfel încât să urmarească normala datorită comenzi comform conturului definit. Asta înseamnă că săgeţile vor fi definite ca lungimi egale, culoarea lor însă va indica amploarea si anume viteza. Graficul indică o viteză constantă la iesire şi că urmează o crestere a acesteia fată de valoarea de intrare pe perioada când are loc îngustarea canalului şi în imediata apropiere din faţa acesteia. Comanda “elevation plot” pe figura definită arată mult mai clar acest lucru.

Vectori schiţaţi mai jos arată, reprezintă câmpul definit al vitezei. De asemenea se observă că orientarea vectorilor este paralelă cu contul (canalul îngustat ) definit, de asemenea se mai observă că viteza la iesire este de două ori mai mare faţă de intrare.

De asemenea din graficul (vm) prezentat mai jos se poate observa ca factorul de crestere al vitezei la intrare faţă de ieşire este 2 prin analogie dacă canalul ar fi largit proportia sar păstra numai că viteza se va micşora cu o proporţie egală. Bineinteles procesul se petrece în concordanţă cu legea conservări masei şi volumului, principiu ce este incorporat în analiza PDE folosind divergenţa.

Graficul de mai jos de arată distribuţia variaţiei presiuni. Se observă că aceasta diferă foarte puţin la intrare faţă de iesire pe conturul canalului definit.

Un alt lucru remarcat este că viteza (vm) este minimă la colturile unde canalul se îngustează şi maximă la colţurile unde canalul se liniarizează din nou.

Din graficul de mai jos se poate observa că presiunea (p) este minimă la colţurile unde viteza işi atinge vârful.

Codul de culoare indică faptul că acest minim este de 70% din valoare de intrare după cum este definită în figura de mai jos mărită.

Din următorul grafic de asemenea se observa o scădere a presiuni la iesire faţă de intrare dar nu in concordanţă cu lăţimea canalului.

Din ultimele doua grafice se observă că atât divergenţa cât şi vorticitatea dispar .

2.Curgerea viscoasă

2.1. Determinarea regimurilor de curgere a lichidelor în conducte

Mişcarea fluidelor viscoase poate fi laminară sau turbulentă, după cum straturile de fluid alunecă paralele unele cu altele, sau acestea se întrepătrund datorită prezenţei unor componente pulsatorii ale vitezei.

Distincţia dintre regimul laminar şi cel turbulent a fost făcută pentru prima oară de către Reynolds, în 1883. Aceasta a arătat că cele două regimuri de curgere pot fi deosebite cantitativ de mărimea adimensională:

ℜ= v ∙dV

= ρ ∙ v ∙ dη

(2.1)

unde d este o lungime care caracterizează frontiera domeniului mişcării, iar în cazul conductelor circulare este chiar diametrul acesteia. Pentru valori Re < 2300, regimul de curgere este laminar, fiind caracterizat printr-o dependenţă liniară între viteza de mişcare şi gradientul de presiune (legea Hagen-Poisseuille). Pentru valori 2300 < Re < 3000 regimul de curgere este tranzitoriu, iar pentru Re > 3000, regimul de curgere devine turbulent.

Scopul lucrării constă în verificarea experimentală a corelaţiei dintre regimul de mişcare a apei

într-un tub cilindric orizontal şi numărul Reynolds

2.2. Curgerea viscoasă în canale

În această parte vor fi descrise procese în (x,y). În aceste condiţii curl(v) va fi în general diferit de zero.

Mecanica clasică aplicată într-un lichid se reduce la ecuaţia Navier-Stokes. Aceasta ecuaţie exprimând legea de miscare a lui Newton:

ρ0dvdt

=f tot (2.2)

Pentru forta totala (ftot ) pe un element lichid care este transportat impreună cu torentul. Aici, ρ0 este constanta masică de densitate a fluidului. Întrucât viteza într-un element de volum definit este o funcţie de (t,x,y), putem scrie:

dvdt

=∂ v∂ t

+ ∂ v∂ x

∙dxdt

+ ∂ v∂ y

∙dydt

=∂ v∂ t

+v x∂ v∂ x

+v y∂ v∂ y

=∂ v∂t

+(v ∙∇ ) v (2.3)

Folosind această expresie pentru derivare, legea lui Newton capătă forma:

ρ0∂ v∂ t

+ ρ0 ( v ∙∇ ) v−F+∇ p−η∇2 v=0 (2.4)

, unde :-F este o forţă externă ( de exempu gravitaţia),

-∇ p forta datorata presiuni - η∇2 v proportionala cu vascozitatea. Acest vector PDE este cunoscut ca ecutia Navier-Stokes (N-S).

Termenul urmator, ρ0 ¿)v are dimensiune de forţă dar în realitate este o parte a derivatei în functie de timp şi prin urmare este denumita ca o forţă inerţială.

Ultimul termen corespunde forţei determinată de viscozitate în volumul determinatCel mai simplu caz de debit apare la viteze foarte mici astfel încăt forţa de inerţie

nonliniară devine neglijabilă in comparaţie cu forţa de vîscozitate. Raportul dintre forţa inerţială si forta determinată de vascozitatea lichidului este exprimată în forma numărului dimensionat Reynolds, definit de:

ℜ=ρ0 v0 L0

η (2.5)

Unde: -V0 este viteza tipica; -L0 marimea domeniului de soluţii.

La valori suficient de mici le Re, termenul inerţial devine neglijabil în comparţie cu forţa determinată de viscozitate şi problema poate fi tratată liniar în variabile dependente. PDE produce soluţii corespunzatoare unui flux laminar.

Înainte de prima valore critică (Re=1), soluţiile pot rămâne laminare, chiar dacă PDE este neliniar. Până la valori mari, soluţia devine turbulent şi dependentă de timp.

În coordonate carteziene, componentele ecuaţiilor Navier-Stokes pot fi definite astfel:

ρ0 {∂ v x

∂ t∂ v y

∂ t}+ρ0 (v ∙∇ ) v−{Fx

F y}+{∂ p

∂ x∂ p∂ y

}−η {∂2 vv

∂ x2 +∂2 vx

∂ y2

∂2 v y

∂ x2 +∂2 v y

∂ y2} (2.6)

Determinând astfel două soluţii pentru 3 variabile vx vy vz şi p. Conservarea masei la densitate constantă ne conduce la cea de a treia ecuatie:

∇ ∙ ( ρ0 v )=ρ0∇ ∙ v=ρ0( ∂ vx

∂ x+

∂ v y

∂ y )=0 (2.7)

dar din nefericire aceaste formule sunt un PDE de ordinul întâi, pe care FlexPDE nu-l va accepta.Folosind ∇ ∙ v , împreună cu ecuaţia mişcări am putea genera o relaţie conţinând o

derivată de ordin doi în p. Aplicând operatorul divergenţă asupra ecuaţiei N-S vom obţine:

ρ0∂∂t∇ ∙ v+ρ0∇ ∙ [ (v ∙∇ ) v ]−∇ ∙ F+∇2 p−η∇ ∙ (∇2 v )=0 (2.7)

unde primul termen se anulează datorită conservări masei.Acum vom putea elimina ultimul termen utilizand indentitatea:

η∇ ∙ (∇2 v )=η∇2 (∇ ∙ v )=η∇2 (0 )=0 (2.8) De unde rezultă restul ecutiei modificate N-S:2

∇2 p+ ρ0∇ ∙ [ ( v ∙∇ ) v ]−∇ ∙ F=0 (2.9) Dacă volumul fortei F este constant în spaţiu, ultimul termen se va anula.Exprimat în sistemul de coordonate cartezian, acest PDE capătă următoarea formă:

∂2 p∂ x2 + ∂2 p

∂ y2 +ρ0∇ ∙ [ ( v ∙∇ ) v ]−∇ ∙ F=0 (2.10)

De aici se obtin trei funcţii PDE pentru calcularea si vx vy vz şi p. Desi sa derivat ecuaţia pentru (p) folosind conservarea masei, este greşit să deducem că orice soluţie a acestor trei ecuatii va fi necesară pentru a satisface ∇ ∙ v=0. De fapt, una ar putea să fie adevarată doar în cazuri speciale. Deci rezolvarea este prezentă pentru urmatoarea formă:

∂2 p∂ x2 + ∂2 p

∂ y2 +ρ0∇ ∙ [ ( v ∙∇ ) v ]−∇ ∙ F−f ∇∇ ∙ v=0 (2.11)

unde vom putea alege în mod liber factorul în f ∇ comform problemei, pentru a asigura dispariţia divergentei.

Factorul ( f ∇ ) nu poate fi luat ca un număr fix, cu toate acestea din moment ce are o

demensiune fizică, de fapt la fel ca η /L02. Deci ar trebui sa scriem f ∇=C

η

L02 , unde parametrul L0

este marimea domeniului. Numarul (C) este , destul de mare pentru asigura disparitia ∇ ∙ v, dar nu atât de mare astfel încât să afecteze convergenţa în calculele FlexPDE.

2.3 Curgere viscoasă printr-un canal îngustat

Programul care urmeaza este o o continuare a programului de la punctul 1.2 modificat astfel incât sa definească o curgere viscoasă. Pfetele de intrare si iesire sa specificat ∂vx /∂x = 0,

presupunănd că are loc o modificare neglijabila a vx aproape de iesire.Programul sursa este:

TITLE 'Curgere viscoasa printr-un canal ingustat' { fex2.pde }COORDINATES cartesian2SELECT errlim=1e-4 ngrid=1 spectral_colorsVARIABLES vx vy pDEFINITIONS Lx=1.0 Ly=1.0 k=0.5 visc=1e4 delp=100 { presiune } dens=1e3 Re=dens*globalmax( vx)*2*Ly/visc v=vector( vx, vy) vm=magnitude( v) unit_x=vector(1,0) unit_y=vector(0,1) { unitate camp vectorial} nx=normal( unit_x) ny=normal( unit_y) {Directia cosinusurilor }

natp=visc*[ nx*div( grad( vx))+ ny*div( grad( vy))] { Conditii naturale în domeniu pentru p: }EQUATIONS vx: dx( p)- visc*div( grad( vx))=0 vy: dy( p)- visc*div( grad( vy))=0 p: div( grad( p))=0BOUNDARIESregion 'domeniu' start 'exterior' (0,Ly)natural( vx)=0 value( vy)=0 value( p)=delp { intrare }line to (0,-Ly) value( vx)=0 value( vy)=0 natural( p)=natpline to (Lx,-Ly) to (2*Lx,-Ly*k) to (3*Lx,-Ly*k)natural( vx)=0 value( vy)=0 value( p)=0 { iesire }line to (3*Lx,Ly*k) value( vx)=0 value( vy)=0 natural( p)=natpline to (2*Lx,Ly*k) to (Lx,Ly) to closePLOTSelevation( nx, ny) on 'exterior' as 'directia cos'contour( vx) report(Re) contour( vm)vector( v) norm contour( p)contour( div( v)) painted contour( curl( v)) paintedelevation( vx) from (0.5*Lx,-Ly) to (0.5*Lx,Ly)elevation( vx) from (2.5*Lx,-Ly*k) to (2.5*Lx,Ly*k)END

Următoarea diagramă arată cum cos( nx) şi cos( ny) se modifică dea lungul canalului.

Soluţia converge rapid şi diagrama (vm) demonstrează (cum se arată în figura următoare) că viteza dispare la marginea pereţilor. Este surprinzător totuşi că valoarea vitezei de ieşire este mai mică decât valoarea ei de intrare.

Urmatoarea figura arată ca difergenta cu sigurantă nu este 0. De aceea se obtine si o valoarea asemănatoare din graficele de altitudine noi. Valoarea integralei pe linia de jos este evident egala cu volumul de lichid transportat prin seciunea transversală (metru/secunda în z). Cele doua valori ale integralei arată ca fluxul tronsoanelor sunt diferite. Pe scurt solutia nu conservă masa si este gresită.

Cauza acestei discrepanţe este faptul că nu am folosit funcţii suplimentare în PDE pentru a suprima div(v) divergenţa.Factorul 1e4 este adecvat pentru a schimba eroarea.

TITLE 'curgere viascoasa printr-un canal ingustat cu divergenta' { fex3.pde }COORDINATES cartesian2SELECT errlim=1e-4 ngrid=1 spectral_colorsVARIABLES vx vy pDEFINITIONS Lx=1.0 Ly=1.0 k=0.5 visc=1e4 delp=100 { presiune } dens=1e3 Re=dens*globalmax( vx)*2*Ly/visc v=vector( vx, vy) vm=magnitude( v) unit_x=vector(1,0) unit_y=vector(0,1) { unitate camp vectorial} nx=normal( unit_x) ny=normal( unit_y) {Directia cosinusurilor } natp=visc*[ nx*div( grad( vx))+ ny*div( grad( vy))]{ Conditii naturale în domeniu pentru p: }EQUATIONSvx: dx( p)- visc*div( grad( vx))=0vy: dy( p)- visc*div( grad( vy))=0 p: div( grad( p))- 1e4*visc/Ly^2*div(v)=0BOUNDARIESregion 'domeniu' start 'iesire' (0,Ly)natural( vx)=0 value( vy)=0 value( p)=delp { intrare }line to (0,-Ly) value( vx)=0 value( vy)=0 natural( p)=natpline to (Lx,-Ly) to (2*Lx,-Ly*k) to (3*Lx,-Ly*k)natural( vx)=0 value( vy)=0 value( p)=0 { iesire }line to (3*Lx,Ly*k) value( vx)=0 value( vy)=0 natural( p)=natpline to (2*Lx,Ly*k) to (Lx,Ly) to closePLOTSelevation( nx, ny) on 'iesire' as 'directie cos'contour( vx) report(Re) contour( vm)vector( v) norm contour( p)contour( div( v))elevation( natp) on 'iesire'elevation( vx) from (0.5*Lx,-Ly) to (0.5*Lx,Ly)elevation( vx) from (2.5*Lx,-Ly*k) to (2.5*Lx,Ly*k)END

Din urmatorul grafic se observa că Vx viteza maximă la iesire este de aproximativ doua ori fata de intrare.

Funcţia elevation peste o sectiune indică o formă parabolică a profilului vitezei peste canal cum este aratat în figura de mai jos.

În urmatorul grafic se prezintă variaţia natp pe intregul contur. Aceasta defineste o derivata de ordin doi de aceea va aparea ca o funcţie ”de scară” (treaptă).

2.3. Comparatie cu curgerea irotatională

Ar putea fi interesant sa comparăm curgerea viscoasă printr-un canal îngustat cu aceea aparţinând unui potenţial de viteza φ. Pentru a corela condiţiile pe domeniu schimbăm distribuţia vitezei la intrare, astfel încât să se realizeze un profil parabolic al vitezei. Definiţia vx0 în primul program (fex1) trebuie modificată si ar trebui să adaptăm comanda (plot) conform noii situaţii.

TITLE 'Curgere printr-un canal ingustat cu distributie parabolica la intrare a vitezei'COORDINATES cartesian2 { Sistemul de coordonate, 1D,2D,3D, etc }SELECT errlim=1e-5 ngrid=1 spectral_colors { Culori Curcubeu }VARIABLES { Definitre variabile in PDE } phi { Potentialul vitezei }DEFINITIONS { Definire parametri } Lx=1 Ly=1 k=0.5 { Coeficientul de comprimare } vx0=4.0e-4*(Ly^2- y^2)/Ly^2 { Viteza de intrare }

p0=1e5 {Presiunea atmosferica } dens=1e3 { Densitate} vx=dx(phi) vy=dy(phi) { Componentele vitezei } v=vector( vx,vy) vm=sqrt( vx^2+ vy^2) { Viteza } p=p0+ 1/2*dens*(vx0^2-vm^2) { Presiunea } div_v=dx( vx)+ dy( vy) { Divergenta, or div( v) } curl_z=dx( vy)- dy( vx) { Vorticitate, or curl( v) }EQUATIONS dxx( phi)+ dyy( phi)=0 { Sau div( grad( v)) }BOUNDARIES { Definirea domeniului }region 'domeniu' start 'exterior' (0,Ly)natural( phi)=-vx0 line to (0,-Ly) { intrare }natural( phi)=0 line to (Lx,-Ly) to (2*Lx,-Ly*k) to (3*Lx,-Ly*k)value( phi)=0 line to (3*Lx,Ly*k) { iesire }natural( phi)=0 line to (2*Lx,Ly*k) to (Lx,Ly) to closeMonitorsPLOTSelevation( vx0) from (0,-Ly) to (0,Ly)elevation( vx) from (Lx/2,-Ly) to (Lx/2,Ly)elevation( vx) from (3*Lx,-Ly) to (3*Lx,Ly)vector( v) normcontour( vx) painted contour( vm) paintedcontour( div(v)) contour( curl(v))ENDComanda ”elevation plots” ne arată cum distribuţia parabolică de la intrare devine

aproapre dreaptă la iesire.Figura de mai jos unde sunt prezentati vectori relevă acest lucru.

2.4. Injectia la viteza uniforma

Pana acum am injectat fluid într-un canal la presiune uniformă. O alternativă ar fi să impunem o viteză de intrare uniformă, ceea ce ar duce la distribuirea neuniformă a presiunii pe suprafaţa de intrare. Condiţiile pe domeniu pentru presiune sunt derivate naturale, cu excepţia ieşirii unde specificăm valoarea p = 0. Pentru a obţine presiunea totală adăugăm valoarea normală. Rezultatule este urmatorul program PDE:

TITLE 'Injectie la viteza uniforma ' { fex5.pde }SELECT errlim=1e-3 ngrid=1 spectral_colorsVARIABLES vx vy pDEFINITIONS Lx=1.0 Ly=Lx k=0.5 visc=1.0 vx0=1e-5 { Viteza de intrare }

dens=1e3 Re=dens*vx0*2*Ly/visc v=vector( vx, vy) vm=magnitude( v) unit_x=vector(1,0) unit_y=vector(0,1) {campul vectorial} nx=normal( unit_x) ny=normal( unit_y) {Directia cos } natp=visc*[ nx*div( grad( vx))+ ny*div( grad( vy))]EQUATIONS vx: dx( p)- visc*div( grad( vx))=0 vy: dy( p)- visc*div( grad( vy))=0 p: div( grad( p))- 1e4*visc/Ly^2*div( v)=0BOUNDARIESregion 'domeniu' start 'iesire' (0,Ly)value( vx)=vx0 natural( vy)=0 natural( p)=natp { intrare }line to (0,-Ly) value( vx)=0 value( vy)=0 natural( p)=natpline to (Lx,-Ly) to (2*Lx,-Ly*k) to (3*Lx,-Ly*k) { perete }natural( vx)=0 natural( vy)=0 value( p)=0 { iesire }line to (3*Lx,Ly*k) value( vx)=0 value( vy)=0 natural( p)=natpline to (2*Lx,Ly*k) to (Lx,Ly) to close { perete }PLOTSelevation( vx) from (0,-Ly) to (0,Ly)elevation( vx, 0.1*dy( vx)) from (3*Lx,-Ly*k) to (3*Lx,Ly*k)elevation( p) on 'iesire' vector( v) norm report(Re)

contour( vx) contour( vy) contour( vm)contour( p) contour( div( v)) contour( curl( v)) paintedEND

În fugura de mai jos sunt ilustrate schimbarile vitezei pe parcursul traversari canalului.

2.5. Similaritate dinamica

TITLE 'Similaritate dinamica' { fex206a.pde }SELECT errlim=1e-3 ngrid=1 spectral_colorsVARIABLES vx vy pDEFINITIONS Lx=1.0 Ly=Lx k=0.5 visc=1.0 vx0=1e-5 { Viteza de intrare }

dens=1e3 Re=dens*vx0*2*Ly/visc

v=vector( vx, vy) vm=magnitude( v) vxp=vx/vx0 vyp=vy/vx0 vp=vector( vxp, vyp) vpm=magnitude( vp) pp=p/(dens*vx0^2) area=area_integral(1) vpm_mean=area_integral( vpm)/area unit_x=vector(1,0) unit_y=vector(0,1) {campul vectorial} nx=normal( unit_x) ny=normal( unit_y) {Directia cos } natp=visc*[ nx*div( grad( vx))+ ny*div( grad( vy))]EQUATIONS vx: dx( p)- visc*div( grad( vx))=0 vy: dy( p)- visc*div( grad( vy))=0 p: div( grad( p))- 1e4*visc/Ly^2*div( v)=0BOUNDARIESregion 'domeniu' start 'iesire' (0,Ly)value( vx)=vx0 natural( vy)=0 natural( p)=natp { intrare }line to (0,-Ly) value( vx)=0 value( vy)=0 natural( p)=natpline to (Lx,-Ly) to (2*Lx,-Ly*k) to (3*Lx,-Ly*k) { perete }natural( vx)=0 natural( vy)=0 value( p)=0 { iesire }line to (3*Lx,Ly*k) value( vx)=0 value( vy)=0 natural( p)=natpline to (2*Lx,Ly*k) to (Lx,Ly) to close { perete }PLOTSelevation( vx) from (0,-Ly) to (0,Ly)elevation( vx, 0.1*dy( vx)) from (3*Lx,-Ly*k) to (3*Lx,Ly*k)elevation( p) on 'iesire' vector( v) norm report(Re)contour( vx) contour( vy) contour( vm)contour( p) contour( div( v)) contour( curl( v)) paintedvector( vp) norm report(Re) contour( vpm) report(vpm_mean)contour( pp) contour( abs( pp)/area)END