Conductie Termica

Capitolul 0

Transferul termic. Noțiuni fundamentale

0.1 Conducția termică

0.2 Convecția termică

0.3 Radiația termică

0.3.1 Emisia radiației

0.3.2 Absorbția radiației

0.3.3 Coeficientul de transfer termic prin radiație

0.4 Moduri combinate de transfer termic

Ecuațiile conducției termice

Ecuația conducției termice în cazul unidimensional (1D)

Ecuația conductiei termice în cazul tridimensional (3D)

0.6 Condiții la limită

0.6.1 Condiții la limită de tip “temperatură impusă” (Dirichlet)

0.6.2 Condiții la limită de tip “densitatea fluxului termic impusă” (Neumann)

0.6.3 Condiții la limită de tip “convecție” (Robin)

Capitolul 1

Metode de investigare și predicție a realității

Experiența

Calculul teoretic

Calculul numeric (modelare și simulare)

Capitolul 2

Ecuații de conservare și clasificarea ecuațiilor

cu derivate parțiale (EDP)

2.1 Forma generală a ecuațiilor model (EDP)

2.1.1 Principii de conservare

2.1.2 Definiții

Forma integrală sau globală a legii de consevare

2.1.4 Forma locală a legii de conservare

2.1.5 Forma conservativă și neconservativă a ecuațiilor model

2. 2 Clasificarea ecuațiilor diferențiale cu derivate parțiale de ordinul doi

2.3 Coordonate cu simplă și dublă influență

2.3.1 Definiții

2. 4 Ecuații parabolice și eliptice

2.4.1 Definiții

2.4.2 Exemple

Capitolul 3

Principalele metode de discretizare

3.1 Introducere

3.2 Metoda diferențelor finite (MDF)

3.3 Metoda elementelor finite (MEF)

3.4 Metode spectrale (MS)

3.5 Metoda volumelor finite (MVF)

3.5.1 Exemplu: conducția termică 1D staționară

3.5.1 Regulile de bază (rgulile lui Patankar)

4.2.5 Exemple

4.4.1 Forma generală a ecuației discretizate

4.4.2 Exemples

4.5 Conducția termică 3D

4.5.1 Forma generală a ecuației discretizate pentru cazul 3D staționar

4.5.2 Forma generală a ecuației discretizate pentru cazul 3D nestaționar

4.4 Conducția termică 2D nestaționară

4.4.1 Forma generală a ecuației discretizată

4.4.2 Exemple

4.3 Conducția termică 2D staționară

4.3.1 Forma generală a ecuației discretizate

4.3.2 Exemple

4.3.3 Aplicarea algoritmului Thomas la probleme 2D (algoritmul linie cu linie)

4.6 Aplicarea metodei volumelor finite în coordonate cilindrice

4.6.1 Forma generală a ecuației discretizate 1D în regim staționar

4.6.2 Condiții la limită

4.6.3 Aplicație numerică

4.6.4 Validare numerică

4.6.5 Forma generală a ecuației discretizate 2D în regim staționar

4.6.6 Forma generală a ecuației discretizate 2D în regim nestaționar

Capitolul 5

Metoda volumelor finite aplicată la probleme de convecție/difuzie

5.1 Introducere

5.2 Ecuația de convecție/difuzie 1D în regim staționar

Discretizare cu schema centrată (schema cu diferențe centrale)

5.2.2.Discretizare cu schema “upwind” (schema descentrată în amonte)

5.2.3Soluția exactă a problemei de convecție/difuzie 1D staționar

Schema exponențială

Schema hibridă

Schema “Power Law” (Legea Putere)

Recapitularea diferitelor scheme

Formularea generală a diferitelor scheme

Scheme cu diferențe de ordin superior pentru probleme de convecție/difuzie

Capitolul 5

Metoda volumelor finite aplicată la probleme de convecție/difuzie

5.1 Introducere

5.2 Ecuația de convecție/difuzie 1D în regim staționar

5.2.1 Discretizare cu schema centrată (schema cu diferențe centrale)

Discretizare cu schema “upwind” (schema descentrată în amonte)

Soluția exactă a problemei de convecție/difuzie 1D staționar

Schema exponențială

Schema hibridă

Schema “Power Law” (Legea Putere)

Recapitularea diferitelor scheme

5.2.10 Exemple

5.2.11 Proprietățile schemelor de discretizare

5.2.12 Proprietățile schemei cu diferențe centrale

5.2.13 Proprietățile schemei “upwind”

5.2.14 Proprietățile schemei hibride

5.3 Convecție/difuzie 1D nestaționar

5.3.1 Forma generală a ecuației discretizate

5.4 Convecție/difuzie 2D și 3D în regim staționar

5.5 Convecție/difuzie 2D în regim nestaționar

5.5.1 Forma generală a ecuației discretizate

Capitolul 6

Metoda volumelor finite aplicată la tratarea

problemelor de cuplaj viteză/presiune

6.1 Introducere

6.2 Reprezentarea gradientului de presiune în ecuațiile cantității de mișcare discretizate

6.3 Reprezentarea ecuației de continuitate discretizate

6.4 Volume de control decalate (deplasate)

6.5 Ecuațiile de conservare a cantității de mișcare

6.5.1 Integrarea ecuației de conservare a cantității de mișcare

6.6 Algoritmul SIMPLE

6.7 Algoritmul SIMPLER

6.8 Algoritmul SIMPLEC

6.9 Algoritmul PISO

6.10 Condiții la limită pentru ecuația de corecție a presiunii

6.11 Comentarii asupra algoritmilor SIMPLE, SIMPLER, SIMPLEC

și PISO

=== Aplicatia_3 ===

Exemplul 3

Se consideră o bară metalică masivă de lungime infinită și secțiune transversală dreptunghiulară (100/80 mm) parcursă de curentul continuu . Să se determine distribuția în regim staționar a temperaturii în secțiunea transversală a barei.

Sunt cunoscute:

conductivitatea termică, ;

coeficientul global de schimb termic prin convecție și radiație pe suprafața barei aflată în contact cu mediul înconjurător, ;

rezistivitatea electrică a materialului barei la , ;

temperatura mediului ambiant, .

Soluție

Ținând cont de simetria barei, domeniul de analiză se reduce la cel prezentat în fig. 3.1. Pe axele de simetrie se specifică condiții la limită de tip Neumann omogene (flux nul).

Ecuația diferențială ce guvernează transferul termic în regim staționar este:

(3.1)

unde termenul sursă iar densitatea de curent , fiind aria secțiunii transversale a căii de curent () iar este rezistivitatea electrică.. Se obține astfel termenul sursă

(3.2)

Se consideră rețeaua de discretizare din fig. 3.2 și exprimarea matematică corespunzătoare a condițiilor la limită. Ecuația discretizată a ecuației diferențiale (3.1), pentru un nod interior al domeniului de calcul (nodul 16 de exemplu), este

(3.3)

unde

,

.

Termenul sursă fiind constant nu mai este necesară liniarizarea acestuia. Valorile coeficienților din ecuația (3.3), pentru rețeaua de discretizare din figura 3.2, sunt

;

;

.

Ecuațiile discretizate pentru nodurile interioare ale domeniului de calcul (nodurile 8 – 11, 14-17, 20 – 23) sunt următoarele:

Ecuația discretizată pentru un nod interior al frontierei “West”. Aceasta se obține integrând ecuația ( ) pe volumul de control hașurat (jumătate de volum de control) din figura 3.3

(x.x)

Integrând ecuația (mm) se obține:

(x.xx)

Ținând cont de condiția la limită și presupunând o variație liniară a temperaturii între nodurile vecine nodului P, se pot exprima derivatele (gradienții) în punctele și și se obține:

(x.xx)

Regrupând termenii din relația (x.xx) se obține forma generală a ecuației discretizate:

(x.xx)

unde:

; ; ;

; ;

Valorile numerice ale coeficienților din ecuația (x.xx) sunt următoarele:

; ;

Ecuațiile discretizate pentru nodurile frontierei “West” (nodurile 2, 3, 4, 5) sunt:

(x.xx)

Ecuația discretizată pentru un nod interior al frontierei “East”. Aceasta se obține integrând ecuația ( ) pe volumul de control hașurat (jumătate de volum de control) din figura 3.4

(x.x)

Integrând ecuația (xx) se obține:

(x.xx)

Ținând cont că pe această frontieră condiția la limită este

(x.xx)

și exprimând derivatele în punctele , și se obține ecuația:

(x.xx)

Regrupând termenii în ecuația (x.xx) se obține ecuația discretizată sub foma generală:

(x.xx)

unde

,

.

Valorile numerice ale coeficienților sunt următoarele:

, ,

,

.

Ecuațiile discretizate pentru nodurile interioare ale frontierei “East” (nodurile 26, 27, 28 și 29) sunt următoarele:

(x.xx)

Ecuația discretizată pentru un nod interior al frontierei “South”. Aceasta se obține integrând ecuația (xx) pe volumul de control hașurat (jumătate de volum de control) din figura 3.5

(x.x)

Integrând ecuația (xx) se obține:

(x.xx)

Ținând cont că pe această frontieră condiția la limită este

(x.xx)

și exprimând derivatele în punctele , și se obține ecuația:

(x.xx)

Regrupând termenii în ecuația (x.xx) se obține forma generală a ecuației discretizate

(3.3)

unde

,

,

.

Valorile numerice ale coeficienților sunt următoarele

.

Se obțin astfel ecuațiile discretizate pentru nodurile interioare ale frontierei “South” (nodurile 7, 13 și 19):

(x.xx)

Ecuația discretizată pentru un nod interior al frontierei “North”. Aceasta se obține integrând ecuația (xx) pe volumul de control hașurat (jumătate de volum de control) din figura 3.6

(x.x)

(x.xx)

Ținând cont că pe această frontieră condiția la limită este

(x.xx)

și exprimând derivatele în punctele , și se obține ecuația:

(x.xx)

Regrupând termenii în ecuația (x.xx) se obține forma generală a ecuației discretizate:

(x.xx)

unde

,

,

.

Valorile numerice ale coeficienților sunt:

Ecuațiile discretizate pentru nodurile interioare ale frontierei “North” (nodurile 12, 18 și 24) sunt următoarele:

(x.xx)

Ecuația discretizată pentru nodul 1. Se obține integrând ecuația (xx) pe volumul de control hașurat (sfert de volum de control) din figura 3.7

(x.xx)

(x.xx)

Ținând cont de condițiile la limită

și (x.xx)

și exprimând derivatele în punctele și se obține ecuația:

(x.xx)

Regrupând termenii în ecuația (x.xx) se obține forma generală a ecuației discretizate pentru nodul 1:

(x.xx)

unde

,

.

Valorile numerice ale coeficienților sunt următoarele:

Cu valorile calculate, ecuația discretizată pentru nodul 1 este:

(x.xx)

Ecuația discretizată pentru nodul 6. Se obține integrând ecuația (xx) pe volumul de control hașurat (sfert de volum de control) din figura 3.8

(x.xx)

(x.xx)

Ținând cont de condițiile la limită

și (x.xx)

și exprimând derivatele în punctele și se obține ecuația:

(x.xx)

Regrupând termenii în ecuația (x.xx) se obține forma generală a ecuației discretizate pentru nodul 6:

(x.xx)

unde

,

.

Valorile numerice ale coeficienților sunt următoarele:

Cu valorile calculate, ecuația discretizată pentru nodul 6 este:

(x.xx)

Ecuația discretizată pentru nodul 25. Se obține integrând ecuația (xx) pe volumul de control hașurat (sfert de volum de control) din figura 3.9:

(x.xx)

(x.xx)

Ținând cont de condițiile la limită

și (x.xx)

și exprimând derivatele în punctele și se obține ecuația:

(x.xx)

Regrupând termenii în ecuația (x.xx) se obține forma generală a ecuației discretizate pentru nodul 25:

(x.xx)

unde

,

,

Valorile numerice ale coeficienților sunt următoarele:

Cu valorile calculate, ecuația discretizată pentru nodul 25 este:

(x.xx)

Ecuația discretizată pentru nodul 30. Se obține integrând ecuația (xx) pe volumul de control hașurat (sfert de volum de control) din figura 3.10:

(x.xx)

(x.xx)

Ținând cont de condițiile la limită

și (x.xx)

și exprimând derivatele în punctele și se obține ecuația:

(x.xx)

Regrupând termenii în ecuația (x.xx) se obține forma generală a ecuației discretizate pentru nodul 30:

(x.xx)

unde

,

,

Valorile numerice ale coeficienților sunt următoarele:

Cu valorile calculate ale coeficienților, ecuația discretizată pentru nodul 30 este:

(x.xx)

Pentru a aplica algoritmul Thomas după direcția coordonatei “y” formele generale ale ecuațiilor discretizate se scriu sub forma următoare:

(x.xx)

Pe baza ecuațiilor (x.xx) s-a elaboratul programul în Fortran TH2DS3 (a se vedea anexa G). În fig. X.xx se prezintă rezultatele obținute cu acest program

Întrucât gradientul de temperatură în cazul 2D se exprimă cu relația,

modulul gradientului de temperatură într-un nod de coordonate se calculează cu relația

unde componentele modulului gradientului după cele două direcții, “x” și “y” se calculează cu relațiile

Observație: În nodurile de frontieră gradientul de temperatură se calculează folosind diferențe progresive sau regresive.

Folosind componentele gradientului de temperatură se poate calcula, de asemenea, densitatea fluxului termic.

Fluxul termic de convecție pe frontiera “EAST” () a domeniului de calcul, pentru unitatea de lungime a barei, se calculează cu relația

[ W ]

unde:

sunt densitatea fluxului termic de convecție corespunzător volumului de control al fiecărui nod de pe frontieră (W/m2) și respectiv aria suprafeței de cedare a căldurii (m2) către mediul ambiant ( și reprezintă numărul de noduri pe direcția “x” și respectiv pe direcția “y”.

In mod asemănător pentru frontiera “SOUTH” se folosesc relațiile:

Folosind relațiile de mai sus, fluxul termic total cedat pe unitatea de lungime (frontierele “East” și “South”) a rezultat de 159.98 W ( și ). În regim staționar, fluxul termic total cedat trebuie să fie egal cu pierderile de energie din bară pe unitatea de timp. Într-adevăr, dacă calculăm pierderile, se obține:

.

Rezultatele prezentate s-au obținut în ipoteza că rezistivitatea electrică este constantă, adică nu variază cu temperatura. În realitate, aceasta variază cu temperatura, în domeniul temperaturilor nu prea mari, conform relației

Aceasta implică și variația termenului sursă. În procesul iterativ de obținere a soluției numerice, termenul sursă se actualizează după fiecare iterație. În această ipoteză se obține o altă distribuție a temperaturii (fig. Xxxx).

Rezultatele au fost validate prin comparație cu rezultatele obținute cu ajutorul programului QuickField, versiunea Student (Student’s QuickField TM 3.40a, Copyright © 1993-1997 TeraAnalysis) care folosește metoda elementelor finite.

Tabelul xx Comparația rezultatelor pe frontiera “West”: Volume finite – QuickField

=== Aplicatia_bara ===

Exemplul 3

Se consideră o cale de curent din cupru având secțiunea transversală dreptunghiulară . La momentul bara se află la temperatura mediului ambiant . La momentul bara este supusă unui regim de încălzire în regim permanent, în curent alternativ, la . Să se determine evoluția în timp (considerând ) și distributia câmpului termic în regim stabilizat de încălzire (se neglijează influența efectului pelicular).

Se cunosc urmatoarele date și proprietati de material ale barei:

coefientul de cedare a căldurii (transmisivitatea termică globală),

conductivitatea termică, ;

densitatea de masă, ;

cáldura specifică, .

Soluție

Ecuația ce guvernează regimul termic tranzitoriu al unei căi de curent este ecuația conducției termice 2D în regim nestaționar

. (2.1)

Profitând de simetrie, se poate calcula distribuția câmpului termic pe un sfert din secțiunea transversală (fig. 2.2).

Domeniul de calcul specificat în fig. 2.2 este discretizat și prezentat în fig. 2.3. Se consideră o rețea de discretizare cu 30 noduri. În programul de calcul însă se va alege un număr de noduri astfel încât, de preferință, să avem (în cazul problemei enunțate, de exemplu 6 noduri pe axa ox si 51 noduri pe axa oy)

Ecuația discretizată pentru un nod interior al domeniului de calcul (nodul 16 de exemplu). Se integrează ecuația (2.1) pe volumul de control hașurat din jurul nodului 16 (fig. 2.3). Utilizând schema total implicitá se obtine:

, (2.6)

unde:

,

,

Termenul sursá (in conditiile neglijarii efectului pelicular) este dat de relatia:

; ; (2.7)

unde – rezistivitatea materialului la temperatura , – densitatea de curent, – sectiunea transversală a barei.

; ;

Termenul sursa inițial va fi:

Valorile numerice ale coeficientilor din ecuatia (5.26) sunt:

;

;

;

;

.

Ecuatia (2.6) se aplică succesiv pentru toate nodurile interioare ale domeniului de calcul (nodurile 8, 9, 10, 11, 14, 15, 16, 17, 20, 21, 22, et 23) si se obtin urmatoatrele ecuatii discretizate:

(2.8)

Ecuația discretizată pentru un nod situat pe frontiera “West” (de exemplu, nodul 4). Se integrează ecuația (2.1) pe volumul de control hașurat (jumătate de volum de control) din fig. 2.4 și se obține:

(2.9)

După integrarea partială (in spațiu și timp a membrului stâng și în spațiu a membrului drept) se obține:

(2.10)

Ținând cont că (condiția la limită pe frontiera “West”), presupunând o variație liniară a gradientului de temperatură și utilizând schema total implicită pentru integrarea în timp, se obține:

,

(2.11)

unde și .

Regrupând termenii în ecuația (2.11), se obține forma generală a ecuației discretizate:

, (2.12)

unde:

;

.

Valorile numerice ale coeficienților din ecuația (2.12) sunt următorii:

;

;

.

Ecuațiile discretizate pentru nodurile interioare ale frontierei “West” (nodurile 2, 3, 4 și 5) sunt urmatoarele:

(2.13)

Ecuația discretizată pentru un nod interior al frontierei “South” (de exemplu nodul 13) Se integrează ecuația (2.1) pe volumul de control hașurat prezentat în figura 2.5 și se obtine:

.

(2.14)

Ținând cont că (condiția la limită pe frontiera “South”), înlocuind gradienții de temperatură din punctele și , și utilizând schema total implicită pentru integrarea în timp, se obține ecuația discretizată sub forma generală urmatoare:

, (2.15)

unde:

;

;

și .

Valorile coeficienților din ecuația (2.15) sunt:

;

;

.

Ecuațiile discretizate pentru nodurile 7, 13 și 19 sunt următoarele:

(2.16)

Ecuația discretizată pentru nodul 1. Se obține integrând ecuația (2.1) pe volumul de control (sfert) hașurat din fig. 2.6. După prima integrare se obține:

(2.17)

În același mod, se obține ecuația discretizată pentru nodul 1 sub forma generală următoare:

, (2.18)

unde:

;

;

et .

Valorile coeficienților din ecuația (2.18) sunt:

;

;

.

Se obține astfel ecuația discretizată pentru nodul 1:

(2.19)

Ecuația discretizată pentru un nod interior al frontierei “East” (de exemplu nodul 28). Se integrează ecuația (2.1) pe volumul de control hașurat din fig. 2.7.

După integrarea în spațiu și timp pentru membrul stâng și respectiv în spațiu pentru membrul drept a ecuației (2.1) se obține:

.

(2.20)

Ținând cont că pe această frontieră condiția la limită este:

, (2.21)

și integrând în timp cu schema total implicită, se obține, înlocuind gradienții de temperatură, ecuația următoare:

.

(2.22)

Regrupând termenii în ecuația (2.22) se obține forma generală discretizată:

, (2.23)

unde:

;

;

, et

Valorile coeficientilor din ecuatia (2.23) sunt:

;

;

;

;

Ecuațiile discretizate pentru nodurile interioare ale frontierei “East” (nodurile 26, 27, 28 și 29) sunt următoarele:

(2.24)

Ecuația discretizată pentru un nod interior al fontierei “North” (de exemplu nodul 18). Se integrează ecuația (2.1) pe volumul de control hașurat din fig. 2.8.

(2.25)

Tinând cont de conditia la limitá pe aceastá frontierá:

. (2.26)

După integrarea ecuației (2.25), se obține forma generală a ecuației discretizate:

, (2.27)

;

și .

Valorile coeficienților din ecuația (2.27) sunt:

;

;

;

;

Ecuațiile ce trebuie rezolvate, pentru nodurile interioare, pe frontiera “North” (nodurile 12, 18 și 24 sunt:

(2.28)

Ecuația discretizată pentru nodul 30. Se integrează ecuația (2.1) pe sfertul de volum de control hașurat din fig. 2.9, adică

Integrând parțial ecuatia (2.29) se obține:

(2.30)

Ținând cont de condițiile la limită pe frontierele “North” și “East” dupa integrare se obține:

. (2.31)

Regrupând termenii in ecuația (2.31), se obține forma generală a ecuației discretizate pentru nodul 30 :

; (2.32)

;

;

et .

Valorile coeficienților din ecuația (2.32) sunt :

;

.

Ecuatia ce trebuie rezolvata pentru nodul 30 este urmatoarea:

(2.33)

Ecuația discretizată pentru nodul 6. Se integrează ecuația (2.1) pe sfertul de volum de control hașurat din fig. 2.10.

. (2.34)

Integrând parțial ecuația (2.34) se obține:

(2.35)

Tinând cont de condițiile la limită pe frontierele “North” și “West”, după integrare, folosind schema total implicită pentru integrarea în timp se obține:

. (2.36)

Regrupând termenii in ecuația (2.36), se obține forma generală a ecuației discretizate pentru nodul 6:

; (2.37)

;

;

și .

Valorile coeficienților din ecuația (2.37) sunt:

;

;

;

Ecuația discretizată ce trebuie rezolvată pentru nodul 6 este următoarea:

(2.38)

Ecuația dscretizată pentru nodul 25. Se integrează ecuația (2.1) pe sfertul de volum de control hașurat din fig. 2.11.

Integrând in acelasi mod ca pentru nodul 6 se obține următoarea formă generală a ecuației discretizate:

; (2.39)

;

;

și .

Valorile coeficienților din ecuația (2.39) sunt:

;

;

.

Ecuația discretizată ce trebuie rezolvată pentru nodul 25 este următoarea:

. (2.40)

Ecuațiile (2.8), (2.13), (2.16), (2.19), (2.24), (2.28), (2.33), (2.38), si (2.40) formează sistemul de ecuații ce trebuie rezolvat (30 ecuații cu 30 necunoscute). Pentru rezolvarea acestui sistem se aplică algoritmul Thomas [2] adaptat pentru probleme 2D. S-a realizat un program de calcul în FORTRAN (codul sursă BARA.for) care permite o discretizare foarte fină a domeniului de calcul în vederea obținerii unei precizii mai bune a rezultatelor. Ținând cont de observațiile de la paragraful 4.3.3 se poate realiza o rețea de discretizare cu 6 noduri pe direcția axei “x” și 51 noduri pe direcția axei “y” astfel încât . Programul de calcul a fost validat [2] cu ajutorul programului comercial QuickField 4.2T (versiunea student) care rezolvă astfel de probleme prin metoda elementelor finite.

2.3 Rezultate

Cu ajutorului programului BARA s-a simulat regimul de incalzire la curentul nominal de 1250 A. Variatia in timp a temperaturii este prezentata in fig. 2.12. Reprezentarea s-a realizat folosind temperatura in centrul barei care este practic egală cu cea de la suprafața barei. In cazul de fata nu se constata diferente între temperatura din centrul barei si temperatura de la suprafata barei (diferenta maxima fiind de 0.01 oC). In cazul in care s-ar tine cont de efectul pelicular este posibil ca temperatura in zona centrala a barei sa scada fata de periferie din cauza scaderii densitatii de curent. Calitativ se poate constata o concordanta intre constanta de timp termică (calculată în capitolul 1) si cea care ar rezulta grafic din fig. 2.12.

În tabelul 3 sunt prezentate temperaturile atinse de bară după 1, 3, 10 și 40 s.

Observație: Se constată că la o temperatură a mediului ambiant de 40 0C bara se încălzește cu numai 9.275 0C. ceea ce înseamnă că după criteriul termic calea de curent este supradimensionată.

Programul de calcul permite pentru o bară de secțiune dreptunghiulară:

simularea regimului permanent de încălzire și determinarea temperaturii de regim staționar pentru o geometrie dată a căii de curent;

dimensionarea căii de curent pentru un raport h/b dat și o temperatură maximă admisă dată;

simularea regimului termic de scurtcircuit având ca date inițiale temperatura de regim staționar a regimului permanent de încălzire.

=== Cap_ 4_A+Br ===

Capitolul 4

Metoda volumelor finite aplicată la probleme

de conducție termică

4.1 Conducția termică 1D staționară

Ecuația diferențială este următoarea:

. (4.1)

Ecuația discretizată pe volumul de control prezentat în figura 4.1 este următoarea (a se vedea relația (3.15) din capitolul 3):

; (4.2)

.

În general , rețeaua de discretizare fiind, în aces caz, neuniformă. Rafinarea rețelei de discretizare se face în zonele cu valori mari ale gradienților.

4.1.1 Determinarea conductivității termice la interfețele volumelor de control

În general , conductivitatea termică variind în funcție de temperatură, , sau chiar în funcție de poziție pentru materialele compozite.

Conductivitatea termică la interfața poate fi determinată prin interpolare liniară între valorile din nodurile și :

. (4.3)

Dacă interfața (fig. 4.2) este situată la jumătatea distanței dintre și , atunci și se obține:

. (4.4)

4.1.2 Conservarea fluxului la interfețele volumelor de control

Dacă se consideră fluxul la interfața (fig. 4.3) se poate scrie:

. (4.5)

Densitatea fluxului termic la interfață poate fi exprimat, de asemenea, astfel:

. (4.6)

Din relația (4.6) rezultă expresia conductivității termice la interfața volumului de control:

(4.7)

Dacă se ține seama de definiția coeficientului , relația (4.7) devine:

. (4.8)

Dacă considerăm cazul particular cu conductivitatea termică la interfața este:

. (4.9)

care reprezintă media armonică a conductivităților termice ale nodurilor vecine ale rețelei de discretizare.

Remarcă

Formula (4.9) este mai recomandată decât formula (4.4):

Dacă (în formula (4.4)), ceea ce semnifică existența unui material izolant în , ceea ce implică ), atunci , și densitatea fluxului termic este dată, la limită, de către relația următoare:

, (4.10)

ceea ce nu convine.

Dacă (în formula (4.9), se obține că ceea ce convine.

Dacă (în formula (4.4), ceea ce implică , , și că temperatura ) atunci ceea ce, de asemenea nu convine.

Dacă (în formula (4.9)), atunci și densitatea fluxului termic la interfața , la limită, devine:

, (4.11)

ceea ce implică , ceea ce convine.

4.1.3 Tratarea neliniarităților

Dacă în ecuația (4.1) și atunci în ecuația discretizată avem

, (4.12)

deci ecuația discetizată (4.2) este o ecuație neliniară.

Rezolvarea sistemului de ecuații algebrice, în acest caz, se face printr-o procedură iterativă. Etapele ce trebuie parcurse, pentru rezolvarea unei probleme neliniare, sunt următoarele:

Se dau valori inițiale arbitrare pentru toate nodurile ale rețelei de discretizare.

Se calculează coeficienții , și .

Se rezolvă sistemul de ecuații algebrice liniare pentru a obține un nou set de valori pentru .

Se iterează (întoarcere la punctul 2.) până când valorile se stabilizează (convergență).

4.1.4 Liniarizarea termenului sursă

Termenul sursă are o dependență neliniară de temperatură. Pentru a obține un sistem de ecuații algebrice liniare, trebuie să liniarizăm termenul sursă în raport cu temperatura T, astfel

, (4.13)

unde este temperatura obținută la iterația precedentă. În același timp trebuie respectată regula Nr. 3, adică trebuie să avem .

La liniarizarea lui ca putem avea situațiile următoare:

care este dat de panta curbei în (dacă această pantă este negativă);

, procesul iterativ diverge (regula Nr. 3 nu este îndeplinită) ;

, există riscul ca procesul iterativ să diveargă (să nu conveargă);

, convergență mai lentă dar asigurată.

În general, există mai multe posibilități de liniarizare a termenului sursă. Acestea sunt ilustrate prin exemplele următoare, [50]:

Exemplul 1

Fie dat, pentru care sunt posibile următoarele liniarizări:

și , este forma optimă recomandată;

și , utilă atunci când expresia lui este mai complicată, ceea ce nu este cazul aici;

și , permite de a avea o pantă negativă mai mare și de a încetini convergența. Acest tip de liniarizare este necesar dacă există multe neliniarități în problemă, întrucât asigură convergența.

Exemplul 2

Fie . În acest caz, liniarizările posibile sunt următoarele:

și , avem divergență întrucât regula Nr. 3 nu este respectată;

și , recomandată atunci când nu putem avea decât ;

și , crearea artificială a pantei negative, , care încetinește convergența.

Exemplul 3

Fie . În acest caz sunt posibile următoarele liniarizări:

și , în acest caz nu se profită de relația cunoscută între și ;

și , liniarizare corectă dar panta nu este destul de mare;

Metoda recomandată. Se caută tangenta la curba în (fig.4.3A)

(4.14)

Astfel, termenul sursă la iterația curentă este:

. (4.15)

Utilizând relația (4.15) pentru termenul sursă considerat ca exemplu, se obține:

, (4.16)

de unde se obține și .

și . Panta este mai mare, ceea ce încetinește convergența.

Dacă considerăm că termenul sursă este dominant în ecuația discretizată se poate aproxima aceasta ecuație astfel

, (4.17)

de unde se obține soluția limită

. (4.18)

Dacă valoarea pantei este mare, soluția tinde către , în timp ce pentru o valoare mai mică a pantei, aceasta implică o mare schimbare a lui , de la valoarea la valoarea . În anumite situații, cazul limită al termenului sursă dominant este utilizat pentru liniarizarea unui termen sursă puternic neliniar și pentru care nu se cunoaște curba analitică (neliniaritatea este dată sub forma unui tabel de date). Dacă presupunem ca pentru valoarea curentă dorim ca la iterația următoare variabila necunoscută să nu depășească valoarea limită , liniarizare se face după cum urmează:

se exprimă derivata în punctul , astfel:

; (4.19)

introducând relația (4.19) în ecuația (4.15) se obține liniarizarea dorită

. (4.20)

4.1.5 Tratarea condițiilor la limită

Se consideră ecuația de conservare a energiei (4.1) în domeniul și rețeaua de discretizare unidimensională prezentată în figura 4.4.

Integrând ecuația (4.1) pe un volum de control (VC), în jurul nodului interior , se obține:

. (4.21)

Exprimând derivatele sub forma discretă, se obține

, (4.22)

unde este valoarea medie a lui pe volumul de control din jurul nodului .

Regrupând termenii, ecuația (4.22) poate fi scrisă astfel:

. (4.23)

Dacă , ecuația (4.23) ia forma următoare:

, (4.24)

unde

.

Pentru punctele (nodurile) și situate la frontierele domeniului de calcul se integrează ecuația (4.1) pe jumătate de volum de control.

Cele trei cazuri tipice de condiții la limită întâlnite, pentru problemele de conducție termică, sunt următoarel:

Temperatură impusă (cunoscută) pe frontieră (condiție de tip Dirichlet);

Densitatea fluxului termic impusă, deci conoscută (condiție de tip Neumann);

Densitatea fluxului termic de convecție specificată printr-un coeficient de schimb termic prin convecție (coeficient de cedare a căldurii ) și o temperatură a fluidului cu care frontiera este în contact () sau printr-un flux de radiație (condiție mixtă sau de tip Fourier)

, (4.24a)

unde este coeficientul de emisie iar este constanta Stefan-Boltzmann.

Condiții la limită de tip Dirichlet

În acest caz, dacă temperatura sau/și este conoscută, nu este necesară scrierea unei ecuații discretizate suplimentară pentru nodul/nodurile sau/și . Când avem o singură condiție la limită de tip Dirichlet (în nodul sau ) numărul ecuațiilor ce trebuie rezolvate simultan este . Dacă avem două condiții la limită de tip Dirichlet (în nodurile de frontieră șii ) atunci numărul ecuațiilor ce trebuie rezolvate simultan este . Pentru din ecuația (4.23) se obține:

. (4.25)

Temperatura fiind cunoscută, termenul trece ca termen sursă, în membrul drept al ecuației, și prima ecuație ce trebuie rezolvată devine:

. (4.26)

Pentru , ecuația (4.26) devient :

. (4.27)

Pentru din ecuația (4.23) se obține:

. (4.28)

Temperatura fiind cunoscută, termenul trece ca un termen sursă și ultima ecuație ce trebuie rezolvată (ecuația pentru nodul ) devine:

. (4.29)

Pentru ecuația (4.29) devine:

. (4.30)

Condiții la limită de tip Neumann

Integrând ecuația (4.1) pe jumătatea de volum de control (VC) ilustrat în figura 4.5 se obține:

, (4.31)

ceea ce înseamnă pentru (în cazul 1D) următoarea integrare:

. (4.32)

După integrare se obține:

. (4.33)

Cum însă densitatea fluxului termic în punctul se exprimă

, (4.34)

atunci ecuația (4.33), exprimând derivata în punctul cu diferențe centrale, devine:

. (4.35)

Regrupând termenii în ecuația (4.35), se obține ecuația discretizată, valabilă pentru nodul de frontieră , pentru condiția la limită de tip Neumann (flux cunoscut):

. (4.36)

Pentru condiția la limită de tip Neumann omogenă avem . Termenul este de fapt media aritmetică între termenul sursă din punctul și termenul sursă din punctul de interfață al jumătății volumului de control.

4.1.6 Algoritmul Thomas sau TDMA (Tri-Diagonal Matrix Algorithm)

Acest algoritm permite calculul soluției unui sistem liniar atunci când matricea este tridiagonală. Acesta este chiar cazul nostru, întrucât ecuațiile discretizate ale sistemului liniar se scriu astfel:

, (4.37)

pentru care rețeaua de discretizare este prezentată în figura 4.6.

Temperatura , din nodul , este exprimată în funcție de temperaturile din nodurile vecine și . Pentru a lua în seamă forma specială a ecuațiilor pentru nodurile dr frontieră și trebuie ca

și (4.38)

Dacă, de exemplu, temperatura este cunoscută avem , , și . Ecuația (4.37) pentru este o relație între , și , dar pentru că este exprimată în funcție de sau este cunoscută, relația între , și , se reduce la o relație între și , adică poate fi exprimată în funcție de . Procesul de substituție continuă până când este exprimată în funcție de care nu joacă nici un rol (), deci se obține, în această etapă, valoarea lui . Se începe în continuare procesul invers în care se determină în funcție de , apoi în funcție de și așa mai departe în funcție de și în funcție de . În procesul de substituție denumit “în înainte” ecuațiile sunt următoarele:

;

;

; (4.39)

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

;

.

În această etapă se caută relații de tipul sub forma,

, (4.40)

dar se pot scrie și relații sub forma următoare:

. (4.41)

Înlocuind ecuația (4.41) în ecuația (4.37) se obține:

. (4.42)

Regrupând termenii în ecuația (4.42) sub forma generală (4.41) se obțin coeficienții și în funcție de coeficienții și :

. (4.43)

Pentru a începe procesul de recurență se constată că pentru ecuația (4.37) este deja sub forma (4.40) și valorile lui și sunt date de relațiile următoare:

și (4.44)

Trebuie precizat că relațiile (4.44) sunt obținute dacă se înlocuiește în relațiile (4.43).

La sfârșitul procesului de recurență se constată că și deci astfel că din ecuația (4.40) se obține

. (4.45)

În acest moment suntem în situația de a începe procesul de substituție “în înapoi” utilizând relația (4.40).

Rezumatul algoritmului

Se calculează și utilizând relațiile (4.44);

Se calculează și , pentru , cu relațiile de recurență (4.43);

Se setează ;

Se utilizează ecuația de la la pentru a obține , , , .

4.1.7 Exemple

Exemplul 1

Se consideră o bară cilindrică (fig. 4.7), fără surse de căldură, având aria transversală și lungimea . Extremitățile, și ale barei sunt menținute la temperaturile constante de de și respectiv .

Să se calculeze distribuția de temperatură de-a lungul barei. Se cunoaște conductivitatea termică .

Soluție

Distribuția temperaturii de-a lungul barei este guvernată de ecuația:

. (4.46)

Considerăm șase puncte (noduri) distribuite uniform de-a lungul barei (fig. 4.8).

Dacă se consideră că și pentru nodurile interioare ale domeniului de calcul (3 și 4 ) ecuația discretizată este

, (4.47)

unde .

Pentru nodurile 2 și 5 se utilizează aceeași ecuație ca pentru un nod interior dar se ține cont că pentru nodurile vecine 1 și 6 temperaturile sunt cunoscute, și respectiv .

Pentru nodul 2 ecuația discretizată este următoarea:

. (4.48)

Pentru nodul 5 ecuația discretizată se scrie astfel:

. (4.49)

Ținând cont că

,

sistemul de ecuații ce trebuie rezolvat este

(4.51)

Regrupând termenii, se obține forma matricială a sistemului de ecuații algeebrice:

ou

Sistemul de ecuații de mai sus comportând un număr redus de ecuații poate fi rezolvat, de exemplu, folosind pachetul software MATHCAD (versiunea MATHCAD 7 Professional, MathSoft Inc., 1997). Soluția ce se obține este

(4.52)

Soluția analitică a ecuației (4.46) reprezintă o distribuție liniară între valorile temperaturii din punctele și , adică:

. (4.53)

În figura 4.9 sunt reprezentate soluțiile analitică și numerică care corespund.

Exemplul 2

Se consideră o placă foarte mare de grosime , având conductivitatea termică constantă și o sursă de căldură uniform distribuită, . Fețele plăcii se află la o temperatură constantă de și respectiv .

Presupunând că dimensiunile plăcii în direcțiile “y” și “z” sunt foarte mari și deci gradientul de temperatură este semnificativ numai pe direcția coordonatei “x”, să se calculeze distribuția de temperatură și să se compare soluția numerică cu soluția analitică.

Soluție

Ecuația diferențială care guvernează distribuția temperaturii este următoarea:

. (4.54)

Domeniul de analiză este discretizat uniform cu șase noduri, ca în exemplul precedent, având . Aria este considerată în planul .

Se integrează ecuația (4.54), pe volumul de control hașurat din figura 4.10, și se obține:

. (4.55)

Prima integrală este evaluată ca la exemplul precedent (exemplul 1). Integrala a doua, care conține termenul sursă, se evaluează considerând o valoare medie a lui pe volumul de control. După integrare din ecuația 4.55) rezultă:

; (4.56)

. (4.57)

Regrupând termenii se obține forma generală a ecuației discretizate, valabilă pentru nodurile 3 și 4, astfel:

, (4.58)

unde:

.

Pentru nodul 2 se utilizează aceeași ecuație discretizată ca pentru un nod interior (nodurile 3 și 4) dar se ține cont că nodul vecin “W” corespunde nodului “1” pentru care temperatura este cunoscută și trece ca un termen sursă suplimentar. Ecuația discretizată este deci:

. (4.59)

Pentru nodul 5 se utilizează aceeași ecuație discretizată ca pentru un nod interior (nodurile 3 și 4) dar se ține cont că nodul vecin “E” corespunde nodului “6” pentru care temperatura este cunoscută și trece ca un termen sursă suplimentar. Ecuația discretizată în acest caz este:

. (4.60)

Ținând cont că și

,

,

lsistemul de ecuații ce trebuie rezolvat este:

Regrupând termenii și înlocuind valorile pentru și se obține sistemul tridiagonal ce trebuie rezolvat:

(4.61)

Soluția numerică a sistemului (4.61) este următoarea:

(4.62)

Soluția analitică se obține integrând ecuația (4.54) și impunând condițiile la limită specificate:

. (4.63)

Compararea rezultatelor numerice, obținute cu metoda volumelor finite, cu soluția analitică este prezentată în tabelul 4.1 și în figura 4.11.

Tabelul 4.1 Comparație: soluție numerică – soluție analitică

Se constată că, în pofida unei rețele grosiere, soluțiile numerică și analitică corespund. Codul sursă în Fortran este prezentat în anexa G.

Exemplul 3

Se consideră o bară cilindrică (fig. 4.12) având aria secțiunii transversale . Una din extremitățile barei este menținută la temperatura constantă de () iar cealaltă extremitate este izolată (fluxul de căldură este nul). Pe toată lungimea barei există un schimb de căldură prin convecție dependent de temperatură. Temperatura mediului înconjurător este .

Să se calculeze distribuția de temperatură și să se compare rezultatele numerice cu soluția analitică. Se cunosc: , .

Soluție

Ecuația diferențială care guvernează transferul termic în acest caz este:

, (4.64)

unde este coeficientul de transfer termic prin convecție, este perimetrul, este conductivitatea termică și este temperatura mediului înconjurător.

Soluția analitică este dată de către relația următoare [48]:

, (4.65)

unde (trebuie precizat că ).

Rețeaua de discretizare utilizată este cea din figura 4.13.

Atunci când , ecuația (4.64) poate fi scrisă astfel:

. (4.66)

Integrând ecuația (4.66) pe volumul de control, din jurul nodului , se obține:

. (4.67)

Prima integrală din ecuația (4.67) va fi tratată ca la exemplele 1 și 2. Cea de a doua integrală, din cauza termenului sursă, se evaluează presupunând că cantitatea ce trebuie integrată este local constantă pe fiecare volum de control și rezultă:

. (4.68)

Pentru nodurile interioare, ecuația discretizată devine:

. (4.69)

Regrupând termenii se obține:

. (4.70)

Pentru un nod interior (nodurile 3, 4 și 5) se obține forma generală a ecuației discretizate:

, (4.71)

unde:

.

Termenul sursă din ecuația de mai sus este identificat ca .

Pentru nodul 2 (pentru nodul 1 nu este necesară scrierea unei ecuații discretizate suplimentare, temperatura fiind cunoscută) este valabilă aceeași ecuație ca pentru un nod interior pentru care nodul vecin “W” corespunde cu nodul “1” pentru care temperatura este cunoscută, . Termenul care conține temperatura nodului “1” este interpretat ca un termen sursă. Astfel, ecuația discretizată, pentru nodul “2”, este:

, (4.72)

unde:

.

Pentru nodul “6” se integrează ecuația (4.66) pe jumătate de volum de control:

; (4.73)

.

Pentru că fluxul în nodul este nul (), se obține ecuația discretizată următoare:

. (4.74)

Regrupând termenii în ecuația (4.74) se obține:

, (4.75)

unde:

.

Ținând cont că

,

sistemul de ecuații algebrice ce trebuie rezolvat este:

Regrupând termenii se obține forma matricială a sistemului de ecuații ce trebuie reyolvat:

(4.76)

Soluția numerică a sistemului de ecuații (4.76) este următoarea:

(4.77)

În tabelul 4.2 se compară soluția numerică cu soluția analitică.

Tabelul 4.2 Comparație: soluție numerică – soluție analitică

Precizia soluției numerice poate fi crescută utilizând o rețea de discretizare mai fină. În tabelul 4.3 și în figura 4.14 sunt prezentate rezultatele numerice și erorile pentru o rețea cu 21 noduri (). Codul sursă în Fortran este prezentat în anexa G.

Tabelul 4.3 Comparație: soluție numerică – soluție analitică

4.2 Conducția termică 1D nestaționară

4.2.1 Forma generală a ecuației discretizate

Ecuația diferențială a conducției termice 1D în regim nestaționar este:

, (4.78)

unde este densitatea de masă (), fiind căldura specifică la presiune constantă ().

Se consideră volumul de control 1D din figura 4.15. Se integrează ecuația (4.78) pe volumul de control și pe un interval de timp de la la astfel:

. (4.79)

Ecuația integrată (4.79) poate fi scrisă astfel:

. (4.80)

unde este aria transversală a volumului de control, – volumul acestuia egal cu .

Dacă temperatura nodului este presupusă aceeași pe volumul de control, partea stângă a ecuației (4.80) poate fi scrisă astfel:

, (4.81)

unde este temperatura la momentul și temperatura la momentul .

Utilizând o schemă cu diferențe centrale, pentru termenii de conducție din membrul drept al ecuației (4.81), se obține:

, (4.82)

Pentru a efectua integralele din membrul drept al ecuației (4.82) trebuie cunoscută variația lui , și în timp. Pentru aceasta există mai multe posibilități, se poate lua temperatura de la momentul , temperatura de la momentul , sau o combinație liniară a temperaturilor de la momentele și . Forma generală a integrării în timp se scrie astfel:

, (4.83)

unde este un factor de ponderare.

Aplicând forma generală de integrare în timp (4.83) punctelor , și se obține:

. (4.84)

Regrupând termenii în ecuația (4.84) se obține:

(4.85)

Dacă se indentifică coeficienții lui și ca și se poate scrie ecuația (4.85) sub forma generală discretizată:

, (4.86)

unde:

;

.

Forma exactă a ecuației discretizate depinde de valoarea factorului . Atunci când se utilizează numai temperaturile , și la momentul , în membrul drept al ecuației (4.86), pentru a calcula la momentul ; o asemenea schemă se numește explicită. Atunci când , se utilizează atât temperaturile la momentul cât și temperaturile la momentul ; schema obținută astfel se numește implicită. Cazul limită când schema se numește total implicită. Dacă schema se numește Crank-Nicolson sau semi-implicită.

4.2.2 Schema explicită

În cazul schemei explicite termenul sursă se liniarizează cu relația . Înlocuind în ecuația (4.86) se obține discretizarea explicite a ecuației conducției termice 1D nestaționar:

, (4.87)

unde:

.

Regula Nr. 2 nu este totdeauna satisfăcută. Coeficientul lui poate fi privit ca fiind coeficientul “vecinului” lui în direcția timp sau un coeficient vecin care face legătura între valorile lui la momentul și cele de la momentul .

Pentru a avea coeficientul lui pozitiv, trebuie ca . În cazul general această condiție devine:

. (4.88)

Dacă și condiția (4.88) devine:

(4.89)

Dacă se notează , se obține numărul lui Fourier care trebuie să fie mai mic sau egal cu :

(4.90)

Relația (4.90) reprezintă criteriul de stabilitate pentru schemele explicite.

Observații

Criteriul de convergență al schemei utilizate, pentru integrarea în timp, rezultă dintr-o considerație fizică (regula Nr. 2).

Dacă se reduce pasul pentru a ameliora precizia în spațiu, trebuie scăzut mult ().

4.2.3 Schema Crank-Nicolson

Înlocuind în ecuația (4.86) se obține discretizarea Crank-Nicolson a ecuației conducției termice 1D nestaționar:

, (4.91)

unde:

.

La momentul mai multe necunoscute sunt prezente în ecuația (4.91), schema este deci implicită și ecuațiile trebuie să fie rezolvate simultan pentru toate nodurile la fiecare pas în timp.

Din punct de vedere matematic schema Crank-Nicolson este necondiționat stabilă, dar din punct de vedere numeric convergența către o soluție acceptabilă din punct de vedere fizic nu este asigurată (de exemplu, soluții oscilante de amplitudine constantă sau descrescătoare).

Regula Nr. 2 este satisfăcută numai atunci când există relația

. (4.92)

Dacă și condiția (4.92) devine:

. (4.93)

Remarcă

Relația (4.93) este mai puțin restrictivă decât relația (4.90) asociată schemei explicite. Precizia schemei Crank-Nicolson este de ordinul doi în timp, deci pentru același pas în timp precizia rezultatelor este mai mare decât în cazul schemei explicite.

4.2.4 Schema total implicită

Atunci când , în ecuația (4.86 ), se obține schema total implicită. Ecuația discretizată este următoarea:

, (4.94)

unde:

.

Regula Nr. 2 este totdeauna verificată, deci schema total implicită (TI) este necondiționat stabilă. Precizia schemei TI este ordinul întâi în timp, deci este necesar ca pasul în timp să fie totuși mic pentru a crește precizia rezultatelor.

4.2.5 Exemple

În continuare se pun în evidență proprietățile schemelor de discretizare explicită și implicită comparând rezultatele numerice, pentru o problemă 1D nestaționară, cu soluția analitică.

Exemplul 1

O placă metalică subțire se găsește inițial la o temperatură uniformă de . La momentul temperatura peretelui “Est” al plăcii este brusc redusă la . Celelalte suprafețe ale plăcii sunt izolate.

Să se utilizeze schema explicită a metodei volumelor finite, pentru un pas de timp adecvat, pentru a calcula distribuția tranzitorie a temperaturii și să se compare rezultatele cu soluția analitică la momentele (i) , (ii) , (iii) ;

Să se calculeze soluția numerică pentru un pas de timp dat de formula (4.90), la momentul și să se compare cu soluția analitică;

Datele problemei sunt: lungimea plăcii , conductivitatea termică și .

Soluție

Ecuația diferențială a conducției termice 1D nestaționar este următoare:

. (4.95)

Condiția inițială este:

Condițiile la limită sunt:

;

.

Soluția analitică este dată de relația următoare [28]:

, (4.96)

unde:

.

Se consideră 6 noduri distribuite uniform pe domeniul de calcul, deci pasul în spațiu este (fig. 4.16).

Pentru un nod interior (nodurile 2, 3 și 4) ecuația discretizată obținută din ecuația (4.87) pentru este următoarea:

, (4.97)

unde:

.

Pentru nodul 1 (nod situat pe frontieră), condiția la limită (temperatura necunoscută) impune obținerea unei ecuații discretizate suplimentare, deci se integrează ecuația (4.95) pe o jumătate de volum de control. Ecuația discretizată astfel obținută, pentru nodul 1, este următoarea:

, (4.98)

unde:

.

Pentru nodul 6 (nod de frontieră), condiția la limită fiind de tip Dirichlet, nu mai este necesară scrierea unei ecuații discretizate suplimentare.

Pentru nodul 5 (nod vecin cu un nod de frontiera pentru care temperatura este cunoscută) se utilizează aceeași ecuație ca pentru un nod interior, dar cum temperatura nodului vecin 6 este cunoscută (C) termenul care conține temperatura nodului 6 trece ca un termen sursă suplimentar. Astfel, se obține ecuația discretizată pentru nodul 5:

, (4.99)

unde:

.

Pasul în timp trebuie să satisfacă condiția de stabilitate (4.90), deci:

.

Pentru că se alege și rezultă:

;

.

După înlocuirea valorilor numerice, în ecuațiile (4.97), (4.98), (4.99) și simplificări, se obține:

(4.100)

În tabelul 4.4 se prezintă un exemplu de calcul, utilizând ecuațiile (4.100) pentru primii doi pași în timp.

Tabelul 4.4

Tabelul 4.5

Rezultatele numerice ale exemplului 1 sunt prezentate în tabelul 4.5. Tabelul 4.6 prezintă rezultatele numerice și analitice la momentele de timp de , și . În figura 4.17 se prezintă o comparație a rezultatelor numerice obținute (6 noduri), pentru diferite valori ale pasului în timp, cu soluția analitică. Comparația rezultatelor numerice cu soluția analitică, la diferite momente de timp este prezentată în figura 4.18 (pentru 21 noduri).

Tableau 4.6

Exemplul 2

Să se rezolve problema de la exemplul 1 utilizând schema total implicită și să se compare rezultatele numerice obținute cu cele obținute cu schema explicită, pentru un pas de timp de .

Soluție

Se utilizează aceeași rețea de discretizare ca în figura 4.16. Ecuația discretizată, folosind schema total implicită, pentru un nod interior al domeniului de calcul (nodurile 2, 3, și 4) este cea descrisă de ecuația (4.94) dar fără termen sursă, adică:

, (4.101)

unde:

;

.

Pentru nodurile 1 și 5 situate pe frontieră și respectiv vecin cu frontiera se impune o tratare specială. Astfel, pentru nodul 1 ecuația discretizată este’următoarea:

, (4.102)

unde:

.

Pentru nodul 5 ecuația discretizată, ținând cont că (cunoscută), este următoarea:

, (4.103)

unde:

.

Chiar dacă schema implicită permite utilizarea unui pas de timp oarecare, în continuare se utilizează totuși o valoare rezonabilă a acestuia , pentru a asigura o bună precizie a rezultatelor. Se obține astfel, după înlocuirea valorilor numerice:

;

.

După înlocuirea valorilor numerice în ecuațiile (4.101), (4.102), (4.103) și după simplificări, se obține:

(4.104)

Ținând cont că forma matricială a sistemului de ecuații algebrice este:

(4.105)

Se constată că ecuația pentru fiecare nod conține temperaturile necunoscute ale nodurilor vecine. Schema implicită implică rezolvarea simultană a sistemului de ecuații (4.105). Valorile temperaturii de la momentul de timp precedent sunt utilizate numai pentru calculul membrului drept al ecuației matriciale (4.105).

Tabelul 4.7 și figura 4.19 prprezintă rezultatele numerice prin comparație cu soluția analitică (pentru o rețea de discretizare cu 6 noduri). Codul sursă în Fortran este prezentat în anexa G (THER1Di2)

În figura 4.20 se prezintă o comparație a rezultatelor numerice la momentul , obținute la utilizarea schemelor explicită și implicită, cu soluția analitică pentru un pas de timp . Se constată că schema explicită, pentru , dă oscilații, în timp ce schema implicită dă rezultate în acord cu soluția analitică. Această constatare arată avantajul schemei implicite care permite utilizarea unui pas de timp mai mare. Trebuie semnalat totodată că o precizie bună se obține utilizând totuși un pas de timp mai mic.

Tabelul 4.7

=== Cap_0r ===

Capitolul 0

Transferul termic. Noțiuni fundamentale

Conceptul de energie este utilizat în termodinamică pentru a preciza starea unui sistem. Este cunoscut faptul că energia nu este nici creată nici distrusă, ci numai transformată dintr-o formă în altă formă. Termodinamica studiază relatia dintre căldură si alte forme de energie în timp ce obiectivul transferului termic (transfer de căldură) este analiza ratei de transfer termic ce are loc într-un sistem. Energia transferată prin transfer de căldură nu este în mod direct măsurabilă dar poate fi apreciată printr-o mărime măsurabilă numită temperatură. S-a constatat experimental că atunci când într-un sistem există o diferență de temperatură, apare un flux de căldură (flux termic) care este orientat dinspre regiunea cu temperatură mai mare către regiunea cu temperatură mai mică. Atunci când într-un sistem există un flux termic, un gradient de temperatură este de asemenea prezent. Cunoașterea distribuției de temperatură într-un sistem este necesară la studiul transferului termic.

Problemele de transfer termic joacă un rol important în aplicatii tehnice, fie că schimburile de căldură trebuie să fie importante și rapide, fie că, dimpotrivă, se caută o izolare termică, fie că se dorește termostatarea cu precizie a unui sistem. În electrotehnică, de exemplu, cunoașterea distribuției de temperatură este necesară pentru a evalua solicitările termice ale unui echipament electric care sunt necesare la dimensionarea diferitelor elemente constitutive (căi de curent, contacte electrice, camera de stingere, etc.)

În studiul transferului termic se disting trei moduri de transmitere a căldurii: conducția, convecția și radiația termică. Conducția termică are loc în solide iar convecția este întâlnită în special în fluide. Transmiterea căldurii prin radiație termică poate avea loc în toate mediile transparente pentru undele electromagnetice. În realitate distribuția de temperatură, într-un mediu, este consecința efectelor cumulate ale celor trei moduri de transmitere a căldurii; este imposibil de a izola un mod de transfer termic de un alt mod. Totuși, pentru simplificarea studiului transferului termic se consideră cele trei moduri de transfer termic în mod separat. De exemplu, se poate studia conducția separat sau conducția cuplată cu convecția, neglijând radiația.

O mărime folosită adesea în studiul transferului termic este densitatea fluxului termic care reprezintă căldura ce traversează unitatea de suprafață în unitatea de timp.

0.1 Conducția termică

Conducția termică este fenomenul de transport al căldurii ce are loc în solide; fenomenul fiind de asemenea prezent în lichidele imobile și într-o măsură mai mică în gaze. Fenomenul microscopic (la scară atomică) ce intervine în cazul conducției termice este propagarea agitației termice a particulelor din zonele mai calde către cele din zonele mai reci. Mecanismul microscopic constă în vibrația moleculară sau atomică (lichide, gaze) și în vibrația cristalină precum și deplasarea electronilor liberi (metale). Conducția termică este deci fenomenul prin care energia este transferată din zonele cu temperatură înaltă către zonele cu temperatură joasă.

Legea lui Fourier (pentru un mediu izotrop, printr-o suprafață izotermă) arată că fluxul termic, prin conducție, într-o direcție dată este proporțional cu aria normală la direcția fluxului termic și cu gradientul de temperatură pe această direcție. Fluxul termic, în direcția , de exemplu, în conformitate cu legea lui Fourier este dat de relația:

(1)

sau dacă se exprimă densitatea fluxului termic:

. (2)

Coeficientul de proporționalitate , numit coeficient de conductivitate termică, depinde de material (natură, structură, temperatură, presiune, densitate, etc.); se măsoară în și este totdeuna pozitiv întrucât căldura (energia termică) se transmite de la zonele calde către zonele reci. Dacă temperatura descrește în direcția pozitivă , atunci este negativ. Fluxul termic și densitatea fluxului termic fiind cantități pozitive în direcția pozitivă , în consecință este necesară introducerea semnului minus în membru drept al expresiilor (1) și (2). Dacă membrul drept al expresiilor (1) și (2) este negativ atunci fluxul termic (și de asemenea densitatea fluxului termic) este orientat în direcția negativă .

În cazul general, în spațiul cu mai multe dimensiuni, legea lui Fourier este dată de relația

, (3)

unde este vectorul normal la aria . Densitatea fluxului termic este:

. (4)

În general conductivitatea termică variază în funcție de temperatură. La temperaturi joase această variație poate fi neglijată.

Exemplul 0.1

Să se determine densitatea fluxului termic și fluxul termic printr-o placă, având aria transversală și grosimea (), când suprafețele sale sunt menținute la temperaturile și respectiv .

Soluție

Dacă gradientul de temperatură este constant atunci distribuția temperaturii în placă este liniară (fig. 0.1).

Dacă se aplică relația (2), densitatea fluxului termic este

.

Densitatea fluxului termic este orientată în direcția pozitivă și valoarea sa este de asemenea pozitivă. Aplicând relația (1) se obține fluxul termic:

.

0.2 Convecția termică

În cazul curgerii unui fluid în contact cu un perete, dacă există o diferență de temperatură între perete și fluid, are loc un transfer termic între fluid și perete ca o consecință a mișcării fluidului în raport cu suprafața peretelui. Acest fenomen de transfer de căldură se numește convecție termică. Propagarea căldurii este realizată prin transport de particule. Sunt cunoscute două tipuri de convecție termică (transport de particule):

naturală (liberă), datorată diferențelor de densitate generate de către gradienții de temperatură (dacă fluidul este izoterm, nu există mișcare de particule);

forțată, datorată unei acțiuni mecanice (ventilator, pompă, etc.).

Pentru a permite acest tip de transport, materia (substanța) trebuie să se deplaseze cu ușurintă.

Fenomenul de convecție liberă se întâlnește atunci când fluidul fiind încălzit în contact cu un corp cald se ridică și este înlocuit de către fluidul rece. Mișcarea fluidului are loc ca urmare a diferenței de densitate. Atunci când fluidul este pus în mișcare în mod mecanic (stimulat de pompe sau ventilatoare) avem o convecție forțată.

Densitatea fluxului termic este dată de legea lui Newton:

, (5)

unde , numit coeficient de convecție sau de cedare a căldurii prin convectie (vezi anexa F), în , depinde de un mare număr de factori și în special de diferența de temperatură dintre solid și fluid.

0.3 Radiația termică

Toate corpurile emit energie din cauza temperaturii lor și această energie emisă se numește radiație termică. Energia radiată de un corp este emisă în spațiu sub formă de unde electromagnetice în conformitate cu teoria lui Maxwell (teoria clasică a undelor electromagnetice) sau sub formă discretă de fotoni în conformitate cu teoria lui Planck. Cele două concepte au fost utilizate pentru studiul transferului termic prin radiație. Exemple de radiație pot fi; radiația solară, radiația unui radiator în infrarosu, radiația unui filament al unei lămpi cu incandescență, radiația arcului electric, etc.

Emisia sau absorbția de energie prin radiație de către un corp este un proces global; adică radiația provenind din interiorul unui corp este emisă prin suprafața sa. Reciproc, radiația incidentă la suprafața unui corp penetrează în profunzimea corpului unde este atenuată pe o distanță foarte mică de la suprafață. Se poate spune deci că radiația este absorbită sau emisă de suprafața unui corp. De exemplu, radiația incidentă pe suprafața unui corp din metal este atenuată pe o distanță de câțiva angströmi de la suprafață; în consecință metalele sunt opace din punct de vedere termic la radiație. Toate corpurile care emit radiație absorb de asemenea radiația primită. Anumite corpuri absorb integral radiația incidentă; acestea sunt numite “corpuri negre” și constituie radiatorul ideal. Alte corpuri, “corpurile gri” absorb doar o fracțiune din radiația incidentă.

0.3.1 Emisia radiației

Densitatea fluxului termic maxim emis de către un corp care are temperatura este dat de legea lui Stefan-Boltzmann

, (6)

unde este temperatura absolută măsurată în , este constanta lui Stefan-Boltzmann () și este numită puterea emisivă a corpului negru. Numai un “ radiator ideal ” sau “ corpul negru ” poate să emită un flux termic conform relatiei (6). Fluxul termic emis de către un radiator real, la temperatura , este mai mic decât fluxul termic dat de relatia (6) și se exprimă astfel

, (7)

unde este emisivitatea care variază între 0 și 1 (pentru toate corpurile reale ).

0.3.2 Absorbția radiației

Dacă un flux termic de radiație este incident la un corp negru atunci acesta este în totalitate absorbit de către corpul negru. Dacă însă fluxul de radiație este incident la corpul real, atunci fluxul absorbit de către corp este dat de relația următoare:

, (8)

unde coeficientul de absorbție variază între 0 și 1 (pentru toate corpurile reale ).

Coeficientul de absorbție al unui corp este în general diferit de coeficientul de emisie . Totuși, în diferite aplicații practice se presupune .

Când două corpuri aflate la diferite temperaturi se află față în față, acestea schimbă căldură între ele prin radiație. Dacă mediul înconjurător este aerul, care este transparent pentru radiație, radiația emisă de un corp traversează mediul înconjurător fără atenuare și atinge celălalt corp, și viceversa. Astfel corpul cald obține un câstig net de căldură ca urmare a schimbului termic prin radiație. Analiza schimbului de căldură prin radiație prin suprafețele corpurilor este în general o problemă complicată.

În figura 0.3 se prezintă o placă caldă având aria suprafeței emisive și emisivitatea . Placa având temperatura este plasată într-un mediu foarte mare de arie () având temperatura . Între placă și mediu se află aer care este transparent pentru radiația termică. Fluxul termic emis de către suprafața este dat de către relația următoare:

. (9)

Mediul având suprafața de arie poate fi aproximat ca un corp negru în raport cu mica suprafață de arie . Atunci, densitatea fluxului termic emis de mediul având suprafața de arie este care este un flux termic incident la suprafata . Fluxul termic absorbit de către suprafata este:

. (10)

Pierderile nete prin radiație (fluxul termic net) prin suprafața reprezintă diferența între fluxul termic emis și fluxul termic absorbit

. (11)

Dacă pentru calculul transferului de căldură între o suprafată mică având aria și mediul înconjurător care are temperatura , se obține relația următoare:

. (12)

Valoarea pozitivă a lui semnifică o pierdere de căldură în placa având suprafata , în timp ce o valoare negativă semnifică un câștig de căldură.

Considerăm de această dată două suprafețe finite având ariile și ilustrate în figura 4. Suprafețele sunt menținute la temperaturile și respectiv , și au emisivitățile și . În această situație numai o parte a radiației emise de suprafața atinge suprafața restul fiind pierdută în mediul înconjurător. Considerații similare pot fi făcute și pentru radiația emisă de suprafața . Analiza transferului termic între cele două suprafețe, pentru un asemenea caz, ar trebui să includă și efectele de orientare a suprafețelor, contribuția radiației mediului precum și reflexia radiației pe suprafețele și . Pentru geometria din figura 0.4 presupunem că fluxul de radiație termică al mediului este neglijabil în raport cu fluxul termic al suprafețelor și . În aceste condiții fluxul net emis prin radiație termică de către suprafața poate fi exprimat sub forma

, (13)

unde este un factor care include efectele orientării suprafețelor și emisivitățile lor. Determinarea acestui factor este o problemă foarte complicată [16].

0.3.3 Coeficientul de transfer termic prin radiație

Pentru a simplifica calculul transferului termic prin radiație, se poate utiliza, în anumite condiții restrictive, coeficientul de transfer termic prin radiație , similar cu coeficientul de ttransfer termic prin convecție:

. (14)

Acest concept poate fi aplicat ecuației (12) care poate fi scrisă astfel

. (15)

Dacă , ecuația (15) poate fi liniarizată astfel

(16)

sau dacă exprimăm densitatea de flux, se obține

. (17)

Comparând ecuațiile (14) și (17) se constată că pentru cazul particular considerat (ecuația (12)), coeficientul de transfer termic prin radiație poate fi definit astfel

. (18)

Exemplul 0.2

Un cilindru având diametrul are o suprafață izolată și cealaltă menținută la temperatura . Suprafața caldă are o emisivitate și este expusă într-un mediu înconjurător având temperatura . Dacă mediul dintre corpul cald și mediu înconjurător este aerul, să se calculeze pierderile prin radiație ale cilindrului către mediul înconjurător.

Soluție

Presupunând că se poate aplica relația (12) și rezultă:

.

0.4 Moduri combinate de transfer termic

Până aici am considerat separat cele trei moduri de transmitere a căldurii, conducția, convecția și radiația. În marea majoritate a situațiilor practice, transferul termic printr-o suprafață are loc simultan prin convecție și prin radiație către mediul ambiant.

În figura 0.5 se prezintă o mică placă, având aria și emisivitatea , care este menținută la temperatura și care schimbă energie prin convecție către un fluid aflat la temperatura cu un coeficient de transfer prin comvecție și prin radiație către mediul înconjurător având temperatura . Densitatea fluxului termic pierdut (transmis) prin placă, prin convecție și radiație, este dată de relația următoare:

. (19)

Dacă al doilea termen din membrul drept al relației (19) poate fi liniarizat și se obține

, (20)

unde este coeficientul de transfer termic prin radiație.

Exemplul 0.3

O placă mică și subtire din metal, având aria , are una din fețe izolată iar cealaltă față este expusă la soare. Placa absoarbe energie solară având densitatea fluxului termic de și pierde această energie prin convecție către aerul înconjurător având cu un coeficient de transfer prin convecție și prin radiație termică către mediul înconjurător presupus un corp negru la . Emisivitatea suprafeței este . Să se determine temperatura de echilibru a plăcii, .

Soluție

Bilanțul de energie prin unitatea de arie a suprafeței expuse poate fi scrisă astfel

sau, introducând valorile numerice, se obține

,

de unde se obține ecuația următoare

.

Soluția ecuației de mai sus este .

Ecuațiile conducției termice

Ecuația conducției termice în cazul unidimensional (1D)

În acest paragraf se demonstrează obținerea ecuației conducției termice 1D nestaționar în coordonate carteziene. Se consideră un perete solid a cărui temperatură este funcție de timp și variază numai după direcția coordonatei . Am văzut în paragraful 0.1 că atunci când temperatura variază de-a lungul coordonatei , există o densitate de flux termic în direcția axei dată de legea lui Fourier sub forma

. (21)

Pentru a generaliza analiza, să presupunem de asemenea o sursă de căldură în placă (de exemplu, pierderi prin efect Joule ) având pierderile specifice . Se consideră un element de volum, având grosimea și aria perpendiculară pe axa , prezentat în figura 0.6.

Ecuația de bilanț pe acest element de volum este următoarea:

(22)

Fie densitatea fluxului termic în punctul în direcția pozitivă a coordonatei pe suprafața a elementului de volum. Atunci fluxul termic care intră în elementul de volum prin suprafața în punctul de abscisă este

. (23)

În mod asemănător, fluxul termic care iese din elementul de volum prin prin suprafața în punctul de abscisă este

. (24)

Câstigul net de flux termic în elementul de volum, prin conducție, este diferența între fluxurile termice de la intrarea si ieșirea din elementul de volum, adică

. (25)

Puterea generată în elementul de volum, având volumul , este

. (26)

Creșterea energiei interne pe unitatea de timp în elementul de volum, rezultând din variația în timp a temperaturii este

, (27)

unde:

= căldura specifică la presiune constantă, ;

= puterea generată pe unitatea de volum, ;

= densitatea fluxului termic în directia , ;

= densitatea materialului, ;

= timpul, .

Dacă termenii din relațiile (25) la (27) sunt introduși în ecuația de bilanț energetic (22) se obține o ecuație de bilanț sub forma următoare

. (28)

Când , primul termen din partea stângă a ecuației (28), prin definiție, devine derivata lui în raport cu și ecuația (28) se scrie astfel

. (29)

Dacă se introduce densitatea fluxului termic dată de relația (21) în ecuația (29) se obține

. (30)

Pentru ecuația (30) nu s-a specificat sistemul de coordonate dar este necesară cunoașterea dependenței ariei de coordonata . Pentru aceasta trebuie să considerăm un anumit sistem de coordonate, carteziene, cilindrice sau sferice.

Coordonate carteziene

Dacă aria nu variază de-a lungul coordonatei atunci ecuația (30) devine

, (31)

care constitue ecuația conducției termice 1D în regim nestaționar în coordonate carteziene.

Coordonate cilindrice

Dacă se înlocuiește în ecuația (30) variabila prin variabila radială și ținem cont că aria este proporțională cu variabila , se obține ecuația conducției termice 1D în regim nestaționar în coordonate cilindrice

. (32)

Coordonate sferice

În mod asemănător dacă se înlocuieste în ecuația (30) variabila prin variabila radială și se ține cont că aria este proporțională cu variabila , se obține ecuația conducției termice 1D, în regim nestaționar, în coordonate sferice

. (33)

Ecuația conductiei termice în cazul tridimensional (3D)

Urmând aceeași procedură ca în paragraful precedent se obține ecuația conducției termice în trei dimensiuni.

Coordonate carteziene

În sistemul de coordonate ecuația conducției termice în regim nestaționar este următoarea

, (34)

dacă mediul este anizotrop. Dacă și ecuația este neliniară (avem o problemă neliniară).

Coordonate cilindrice

În sistemul de coordonate cilindrice ecuația conducției termice în regim nestaționar este următoarea

. (35)

0.6 Condiții la limită

Pentru rezolvarea unei ecuații de transfer termic, în regim nestaționar, sunt necesare condiții inițiale și la limită adecvate. Condițiile inițiale specifică distribuția temperaturii la momentul . Condițiile la limită specifică conditiile termice pe frontiera domeniului de calcul. De exemplu, pe o frontieră se pot specifica distribuția temperaturii, distribuția densității fluxului termic și/sau un transfer termic prin convecție către un fluid având temperatura și coeficientul de transfer termic prin convecție cunoscute.

În continuare se prezintă reprezentarea matematică pentru trei tipuri de condiții la limită denumite temperatură impusă (condiție Dirichlet), densitatea fluxului termic impusă (condiție Neumann) și condiție la limită de tip convecție (condiție Robin).

0.6.1 Condiții la limită de tip “temperatură impusă” (Dirichlet)

În numeroase aplicații practice, temperatura se consideră cunoscută pe frontiera domeniului de calcul. De exemplu, suprafața unei frontiere în contact cu un bloc de gheața care se topește este menținută la temperatură constantă și uniformă sau în alte situații distribuția temperaturii este cunoscută în funcție de timp.

Se consideră o placă având grosimea prezentată în figura 0.7. Presupunem că suprafața frontierei la este menținută la temperatura uniformă iar suprafața frontierei de la este menținută la temperatura uniformă .

Condițiile la limită pot fi scrise astfel

; (36)

. (37)

În anumite cazuri, distribuția temperaturii la frontiere este specificată în funcție de poziție și timp. Când valoarea temperaturii pe frontieră este specificată, se spune că avem condiție la limită de primul tip.

0.6.2 Condiții la limită de tip “densitatea fluxului termic impusă” (Neumann)

În anumite cazuri, la frontiera domeniului de calcul, densitatea fluxului termic este cunoscută. De exemplu, când o suprafată a domeniului de calcul este încălzită cu ajutorul unei instalații electrice densitatea fluxului termic pe această suprafață este cunoscută.

Se consideră o placă de grosime prezentată în figura 0.8. Să presupunem că prin suprafața frontierei domeniului de calcul la , densitatea fluxului termic este iar prin suprafața frontierei la densitatea fluxului termic este .

În continuare se prezintă reprezentarea matematică a acestui tip de condiții la limită. La , sursa de căldură exterioară care intră prin suprafața frontierei trebuie să fie egală cu densitatea fluxului termic de conducție din placă, adică

. (38)

În mod asemănător la densitatea fluxului termic exterior care intră prin suprafața frontierei trebuie să fie egală cu densitatea fluxului termic de conducție din placă, adică

. (39)

În relațiile (38) și (39) o valoare pozitivă pentru și are semificația că sursa de căldură (densitatea fluxului termic) intră în placă, în timp ce o valoare negativă a lui și are semnificația că sursa de căldură (densitatea fluxului termic) iese din placă. Aceste condiții la limită se mai numesc și condiții la limită de al doilea tip.

Dacă este densitatea fluxului de radiație, condiția la limită este

unde este temperatura mediului înconjurător (a fluidului).

0.6.3 Condiții la limită de tip “convecție” (Robin)

În unele aplicații practice este posibil ca transferul termic la frontiera domeniului de analiză să se facă prin convecție, cu un coeficient de transfer termic prin convecție către fluidul cu care este în contact și având temperatura cunoscută. Se consideră o placă de grosime prezentată în figura 0.9. Presupunem că un fluid aflat la temperatura curge pe suprafața plăcii de la , având un coeficient de transfer termic prin convecție . Formularea matematică a acestui tip de condiție la limită este obținută scriind bilanțul energetic pe suprafața de la , astfel

sau

, (40)

care este condiția la limită pe frontiera (condiție la limită de al treilea tip).

Dacă fluidul aflat la temperatura , având coeficientul de transfer termic prin convecție pe suprafață , curge pe suprafața plăcii, bilanțul energetic pe această suprafață la , este

sau în formulare matematică

, (41)

care este condiția la limită pe frontiera de la (condiție la limită de al treilea tip).

Exemplul 0.4

Se consideră conducția termică în regim staționar 1D într-o placă care are conductivitatea termică constantă în domeniul . Fluxul termic generat pe unitatea de volum, în placă, este . Suprafața de la frontiera este izolată isolée iar suprafața frontierei de la cedează un flux termic prin convecție, către mediul înconjurător, cu un coeficient de transfer termic prin convecție .

Să se scrie formularea matematică a acestei probleme de conducție termică.

Soluție

Ecuația conducției termice se obține pornind de la ecuația (30) particularizată pentru regimul staționar, înlocuind termenul sursă și considerând conductivitatea termică constantă. Condiția la limită, la , este obținută din ecuația (38) considerând iar la se obține din ecuația (41). Formularea matematică a ecuației și a condițiilor la limită este următoarea

;

;

.

Exemplul 0.5

Se consideră conducția termică în regim staționar 1D într-un cilindru având conductivitatea termică constantă în domeniul . Fluxul termic generat pe unitatea de volum este , iar suprafețele frontierelor interioară și exterioară cedează (transmite) o densitate de flux termic prin convecție. Coeficienții de transfer termic prin convecție pentru suprafețele interioară și exterioară sunt și respectiv. Temperatura fluidului în contact cu suprafața interioară este și respectiv pentru fluidul în contact cu suprafața exterioară.

Să se scrie formularea matematică a acestei probleme de conducție termică.

Soluție

Ecuația conducției termice se obține plecând de la ecuația (32) particularizată pentru regimul staționar, înlocuind și considerând conductivitatea termică constantă. Condițiile la limită la și se obțin din ecuațiile (40) și respectiv (41) considerând însă coordonatele cilindrice. Formularea matematică a problemei enunțate este următoarea

;

;

.

=== Cap_1r ===

Capitolul 1

Metode de investigare și predicție a realității

Predicția transferului termic și a fenomenelor care au loc la curgerea fluidelor poate fi obținută prin trei metode: experiență, calcul teoretic și calcul numeric.

Experiența

Principiul. Prin experiență, se încearcă izolarea sau reproducerea parțială sau totală a unui fenomen fizic folosind modele la scară normală sau redusă.

Avantaje. Experiența furnizează informația cea mai sigură despre un fenomen fizic (în comparație cu modelele teoretice).

Dezavantaje. Modelele la scară normală sunt adesea foarte costisitoare. Modelele la scară redusă sunt mai puțin costisitoare dar extrapolarea rezultatelor la scară normală este uneori dificilă. Traductoarele și sondele de măsură incluse în model generează perturbații și reprezintă deci surse de erori.

Calculul teoretic

Principiul. Principiul calculului teoretic este ilustrat în figura 1.1

Avantaje

Calculul teoretic:

nu necesită, în general, mijloace puternice de calcul;

nu este costisitor;

furnizează soluții exacte;

furnizează rezultate cu o viteză foarte mare.

Dezavantaje

domeniul de aplicare este extrem de limitat (domeniul simplificat);

condiții de calcul idealizate;

geometrii simple;

fenomene liniare sau slab neliniare;

rar poate fi aplicat pentru probleme 3D.

Calculul numeric (modelare și simulare)

Modelul matematic constituit dintr-o ecuație diferențială cu derivate parțiale (EDP) sau dintr-un sistem de ecuații cu derivate parțiale este transformat, cu ajutorul unei metode de discretizare, într-un sistem de ecuații algebrice.

Principiul

Principiul calculului numeric este ilustrat în figura 1.2. Metodele de discretizare cele mai cunoscute sunt:

metoda diferențelor finite (MDF);

metoda elementelor finite (MEF);

metode spectrale (MS);

metoda volumelor finite (MVF).

Algoritmul de rezolvare numerică implică uneori utilizarea metodei integrării temporale și de decuplaj, ca de exemplu algoritmii SIMPLE, SIMPLER, SIMPLEC și PISO, [50].

Avantaje

Calculul numeric:

permite calculul unei soluții numerice pentru aproape toate problemele practice care au un model matematic;

are un cost foarte scăzut, având tendință de scădere;

prezintă rapiditate (permite modificarea geometriei și a condițiilor la limită, etc.);

prezintă o informație completă despre toate câmpurile de mărimi fizice, în toate punctele domeniului și la orice moment;

are posibilitatea de a simula condiții reale speciale;

are posibilitatea de a simula condiții ideale.

Dezavantaje

totul depinde de modelul matematic inițial;

dificultăți în a selecționa soluția “bună” în cazul unor soluții multiple ale problemei matematice;

uneori este mai costisitoare decât experiența.

Experiență sau simulare ?

Se poate spune că “experiența primează iar simularea urmează” întrucât experiența este:

singura care poate pune în evidență fenomene noi de bază;

indispensabilă pentru validarea codurilor numerice.

Totodată “experiențele și simulările numerice sunt complementare”, întrucât simularea permite reducerea costului proiectelor, prin reducerea numărului de experiențe.

Care sunt componentele unui software de modelare și simulare numerică ?

Aplicarea unei metode numerice de calcul implică elaborarea unui software asociat acesteia. Dacă, cu ocazia primelor încercări de calcul prin metode numerice, se scria, pentru fiecare nouă problemă, un program de calcul diferit ținând cont de geometria sa particulară, de particularitățile fizice și de condițiile la limită, astăzi realizatorii de software se orientează către elaborarea de software general a cărui structură informatică este adaptată la rezolvarea unui mare număr de probleme de același tip. Aceste tipuri de software au, practic toate, trei componente principale: (1) un preprocesor, (2) un procesor de calcul (numit solver în engleză) și (3) un postprocesor. În continuare sunt examinate funcțiile fiecarei componente într-un context general.

Preprocesorul

În această etapă a calculului numeric se realizează următoarele activități:

Descrierea geometriei domeniului de calcul (analiză) – definirea domeniului de calcul;

Generarea rețelei de discretizare – discretizarea domeniului de calcul în elemente finite sau volume de control (volume finite);

Alegerea fenomenului fizic ce trebuie modelat;

Definirea proprietăților fizice;

Specificarea condițiilor la limită adecvate problemei enunțate.

Procesorul de calcul

Acesta poate fi realizat, în general, utilizând patru metode numerice distincte: metoda diferențelor finite, metoda elementelor finite, metode spectrale și metoda volumelor finite.

Orice metodă numerică comportă trei pași distincți:

Aproximarea variabilei necunoscute prin diferite tipuri de funcții simple;

Discretizarea prin substituția aproximării în ecuațiile diferențiale și obținerea unui sistem de ecuații algebrice;

Rezolvarea sistemului de ecuații algebrice.

Principala diferență între cele patru metode numerice constă în tipul de aproximare al variabilei necunoscute precum și în procesul de discretizare.

Metoda diferențelor finite

Variabila necunoscută este descrisă prin mai multe valori în punctele (nodurile) unei rețele de discretizare. Dezvoltarea necunoscutei în serie Taylor trunchiată este utilizată pentru aproximarea derivatelor necunoscutei, în fiecare nod al rețelei de discretizare, prin diferențe finite, utilizând necunoscutele din nodurile vecine. Înlocuind derivatele în ecuațiile diferențiale prin diferențe finite, se obține un sistem de ecuații algebrice pentru valorile necunoscutei din fiecare nod al rețelei de discretizare.

Metoda elementelor finite

În cazul metodei elementelor finite se utilizează funcții liniare sau cuadratice, pe fiecare element, pentru a descrie variația locală a necunoscutei . Ecuațiile diferențiale sunt verificate exact de către soluția exactă. Înlocuind aproximația lui în ecuațiile diferențiale se constată că aceasta nu verifică exact. Pentru a măsura eroarea se definește un rezidu. Acest rezidu este minimizat prin multiplicarea lui cu o funcție de pondere și integrându-l. Rezultatul integrării este un set de ecuații algebrice având ca necunoscute coeficienții funcțiilor de aproximare.

Metode spectrale

Necunoscuta este aproximată cu serii Fourier trunchiate sau cu serii de polinoame Cebâsev. În raport cu metoda diferențelor finite și metoda elementelor finite, aproximația nu este locală dar este validă pe întregul domeniu de calcul. Se utilizează și în acest caz conceptul de rezidu ponderat ca în cazul metodei elementelor finite, impunând condiția ca aproximarea să corespundă cu soluția exactă pentru nodurile rețelei de discretizare.

Metoda volumelor finite

La început, metoda a fost dezvoltată ca o formulare specială a metodei diferențelor finite. Algoritmul numeric are pașii următori:

Domeniul de analiză (de calcul) este divizat în volume finite (generarea rețelei de discretizare);

Integrarea formală a ecuațiilor pe toate volumele de control;

Discretizarea, care implică substituirea termenilor integrați, ce reprezintă diferite procese de curgere cum ar fi convecția, difuzia și termenul sursă, prin diferite aproximări de tip diferențe finite. Rezultatul acestui pas îl reprezintă conversia integralelor într-un sistem de ecuații algebrice;

Rezolvarea sistemului de ecuații algebrice prin utilizarea unei metode iterative.

Primul pas, integrarea pe volumul de control, face distincția între metoda volumelor finite și toate celelalte tehnici (metode) numerice. Rezultatul integrării exprimă conservarea exactă a mărimii fizice, , pe fiecare volum de control. Această relație clară între algoritmul numeric și principiul fizic de conservare determină principala atracție a metodei volumelor finite iar conceptul său devine mai ușor de înțeles decât conceptul metodei elementelor finite sau al metodei spectrale.

Conservarea variabilei generale de curgere , de exemplu o componentă a vitezei sau a entalpiei (temperatura), pe volumul de control poate fi exprimată ca un bilanț între diferitele procese care tind s-o crească sau s-o scadă,

Postprocesorul

La interiorul acestui modul, se prezintă utilizatorului rezultatele sub o formă adaptată percepției sale asupra problemei fizice rezolvate. De exemplu, pentru o problemă de conducție termică, trasarea izotermelor este un rezultat util, mai ales daca software-ul are posibilitatea de vizualizare a curbelor izoterme asociate unor valori preferate ale temperaturii.

Postprocesorul oferă diferite facilități de vizualizare și interpretare a rezultatelor cuprinzând de asemenea:

vizualizarea geometriei și a rețelei de discretizare;

vizualizarea câmpurilor de vectori ai mărimilor calculate;

vizualizarea izoliniilor diferitelor mărimi;

vizualizarea suprafețelor 2D și 3D ale mărimilor calculate;

posibilități de a exporta diferite mărimi sub formă de fișiere;

facilități de animație.

=== Cap_2r ===

Capitolul 2

Ecuații de conservare și clasificarea ecuațiilor

cu derivate parțiale (EDP)

2.1 Forma generală a ecuațiilor model (EDP)

2.1.1 Principii de conservare

Toate ecuațiile cu derivate parțiale (ecuațiile model) rezultă dintr-un principiu (lege) de conservare al unor cantități specifice (masice sau volumice ). Cantitatea specifică este o caracteristică fizică a fluidului numită “variabilă dependentă” unde sunt coordonatele spațiale și temporală numite “variabile independente”. Variabila dependentă poate fi una din mărimile fizice prezentate în tabelul 2.1 pentru diferite ecuații de conservare.

Tabelul 2.1 Mărimile fizice masice și volumice pentru diferite ecuații de conservare.

2.1.2 Definiții

Dacă dimensiunea lui este atunci avem:

, densitatea de suprafață a fluxului variabilei pe unitatea de timp ;

, componenta lui pe direcția ;

, sursa lui pe unitatea de volum: .

Fie “” domeniul din fluid având volumul și suprafața fig.2.1)unde este fluxul care influențează variabila dependentă . Se consideră un volum de control cu dimensiunile , și prezentate în figura 2.2. Fluxul care este componenta lui în direcția intră prin interfața de arie iar fluxul care iese prin interfața opusă este . Fluxul net prin interfața de arie este . Dacă considerăm de asemenea contribuțiile pe direcțiile și , fluxul net pe unitatea de volum este:

Forma integrală sau globală a legii de consevare

Termenii care intervin în bilanțul lui din legea de conservare sunt următorii:

variația în timp a lui în volumul :

; (2.1)

fluxul lui prin suprafața :

; (2.2)

producția (generarea) de în volumul :

; (2.3)

Ecuația legii de conservare devine:

; (2.4)

2.1.4 Forma locală a legii de conservare

variația în timp a lui în volumul este:

; (2.5)

variația volumică netă a fluxului prin toate fețele volumului elementar este:

(2.6)

producția (generarea) de în volumul elementar este .

În final ecuația legii de conservare este următoarea:

. (2.7)

Pentru că este generat prin convecție și difuzie, se poate scrie:

, (2.8)

unde este coeficientul de difuzie. Forma generală aa legii de conservare devine:

, (2.9)

Pe direcția “” ecuația se scrie astfel:

. (2.10)

În tabelul 2.2 se prezintă variabilele , și pentru diferite ecuații de conservare.

Tabelul 2.2 Variabilele , și pentru diferite ecuații de conservare.

– coeficient de vâscozitate dinamică);

– coeficient de vâscozitate de expansiune volumică;

– conductivitatea termică (coeficient de difuziune termică);

– căldura specifică masică la presiune constantă;

– coeficient de difuzie masică (legea lui Fick) ;

– tensorul eforturilor vâscoase;

– accelerația gravitațională;

– coeficient de dilatare termică;

– entalpia masică;

– fracțiunea masică.

2.1.5 Forma conservativă și neconservativă a ecuațiilor model

Ecuația (2.9) care este forma conservativă generală a ecuațiilor model poate fi explicitată astfel

. (2.11)

Regrupând termenii se obține

, (2.12)

unde termenul de mai jos reprezintă derivata substanțială

, (2.13)

iar termenul

, (2.14)

reprezintă ecuația de continuitate.

În cele din urmă ecuația (2.12) se scrie sub forma neconservativă astfel:

. (2.15)

Remarcă: Metoda volumelor finite se aplică numai formei conservative a ecuațiilor model.

2. 2 Clasificarea ecuațiilor diferențiale cu derivate parțiale de ordinul doi

În general ecuațiile cu derivate parțiale se clasifică în trei categorii, numite eliptice, parabolice, și hiperbolice. Pentru a ilustra această clasificare se consideră forma cea mai generală a ecuației diferențiale de ordinul doi cu două variabile independente et :

. (2.16)

Se presupune că ecuația (2.16) este liniară (această restricție nu este esențială), adică coeficienții , , , și sunt funcții de două variabile independente și acești coeficienți nu depind de variabila dependentă .

Clasificarea se face pe baza coeficienților , și ai derivatelor de ordinul doi, în funcție de valoarea determinantului:

.

Ecuația diferențială este numită:

; (2.17)

; (2.18)

. (2.19)

Exemple

Ecuația conducției termice fără termen sursă,

, (2.20)

este de tip eliptic.

Ecuația conducției termice cu termen sursă,

, (2.21)

este de asemenea de tip eliptic.

O caracteristică a ecuațiilor eliptice este că acestea necesită condiții la limită adecvate pentru toate frontierele domeniului de calcul.

Ecuația conducției termice unidimensională în regim instaționar,

, (2.22)

este de tip parabolic.

Ecuația undei de ordinul doi,

, (2.23)

unde este timpul, este viteza de propagare a undei este coordonata spațială, este de tip hiperbolic.

Criteriul de clasificare prezentat mai sus este un criteriu pur matematic. Dacă considerăm, de exemplu, ecuația (2.20) sau (2.21) se constată că aceste ecuații sunt de ordinul doi atât pentru variabila cât și pentru variabila . Condițiile într-un anumit punct al domeniului de calcul sunt influențate de către schimbările condițiilor de o parte și de alta a acestui punct atât pe direcția variabilei cât și în direcția variabilei . Ecuația conducției termice în regim staționar 2D este eliptică în și în și deci este numită pur și simplu eliptică.

Să considerăm acum ecuația conducției termice 1D în nestaționar (2.22) care este de ordinul doi în variabila și de ordinul întâi în variabila . Ecuația este privită ca o ecuație de tip eliptic în variabila . Totuși, în ceea ce privește variabila , condițiile la un anumit moment de timp sunt influențate numai de către schimbările de la momentul de timp precedent. Ecuația este numită parabolică în timp sau pur și simplu parabolică. Avantajul la rezolvarea unei ecuații parabolice constă în reducerea timpului de calcul și a memoriei folosită pentru stocarea datelor. De exemplu, pentru ecuația conducției termice 2D în regim nestaționar, care este parabolică după timp, trebuie considerată numai o rețea de discretizare 2D pentru câmpul de temperatură Câmpul de temperatură la nu este afectat de câmpul de temperatură de la momentul . Se începe calculul de la un câmp de temperatură inițial cunoscut și apoi se calculează succesiv câmpul de temperatură la diferite momente de timp. În cazul ecuațiilor hiperbolice este imposibil de a avea avantajul care există pentru ecuațiile parabolice

2.3 Coordonate cu simplă și dublă influență

2.3.1 Definiții

O coordonată sau se spune că este est cu dublă influență dacă condițiile într-o poziție dată în direcția coordonatei, sunt influențate de către schimbările condițiilor de o parte sau de alta a acestei poziții.

Exemplu: Într-o placă metalică (fig. 2.3) când temperatura crește, temperatura punctului vecin va crește de asemenea. Dacă crește, va crește de asemenea.

O coordonată sau se spune că este cu simplá influență dacă condițiile într-o poziție dată în direcția coordonatei sunt influențate de către schimbările condițiilor numai în una din vecinătățile acestei poziții.

Exemplu: timpul este totdeauna o variabilă cu simplă influență , adică ceea ce se va petrece în viitor nu influențează niciodată prezentul (ceea ce se va petrece la momentul nu influențează niciodată ce se petrece la momentul ).

Remarcă: În mod normal coordonatele spațiale sunt variabile cu dublă influență, dar sub efectul unei curgeri (convecție) ele pot deveni cu simplă influență.

2. 4 Ecuații parabolice și eliptice

2.4.1 Definiții

Fie o variabilă dependentă care verifică ecuația

,

(2.24)

unde sunt variabile independente.

Dacă toți coeficienții sunt nenuli și de același semn ecuația este eliptică:

.

Dacă toți coeficienții sunt nenuli și de același semn, cu o singură excepție, ecuația este hiperbilică:

.

Dacă unul din coeficienții este nul și dacă coeficientul corespondent este nenul, ecuația este parabolică:

.

În cazul în care verifică ecuația de conservare

, (2.25)

se poate arăta că această ecuație este:

de tip eliptic dacă toate variabilele independente ale ecuației de conservare sunt cu dublă influență;

de tip parabolic dacă ce puțin una din variabilele independente este cu simplă influență.

2.4.2 Exemple

Conducție 1D nestaționară,

, (2.26)

unde este o variabilă cu simplă influență, deci ecuația este parabolică în , în timp ce este o variabilă cu dublă influență, deci ecuația este eliptică în . Ecuația este parabolică ținând cont că , , .

Conducție 1D staționară,

, (2.27)

unde este o variabilă cu dublă influență și , deci ecuația este eliptică.

Conducție 2D staționară,

, (2.28)

unde și sunt variabile cu dublă influență, deci ecuația este eliptică în și eliptică de asemenea în ) și pentru că ecuația este eliptică.

Conducție 2D și convecție 1D, staționare,

, (2.29)

unde și sunt variabile cu dublă influență, deci ecuația este eliptică în și eliptică de asemenea în ) și pentru că și ecuația este eliptică.

=== Cap_3r ===

Capitolul 3

Principalele metode de discretizare

3.1 Introducere

Există două mari familii de metode de discretizare:

Metode de aproximare a ecuațiilor. Potrivit acestor metode, se caută o soluție exactă a ecuațiilor aproximate întrucât operatorii diferențiali sunt discretizați pe o rețea (metoda diferențelor finite și metoda volumelor finite);

Metode de aproximare a soluțiilor. Potrivit acestor metode, se caută o soluție aproximată a unor ecuații exacte. Soluțiile sunt scrise ca serii de funcții trunchiate la ordine de precizie dorite (metodele spectrale și metoda elementelor finite).

3.2 Metoda diferențelor finite (MDF)

Principiul

Domeniul de calcul este discretizat într-un număr finit de puncte pe care se aproximează operatorii de derivare ai ecuațiilor model prin dezvoltări în serii Taylor trunchiate la un ordin de precizie ales.

Exemplu

Fie rețeaua de discretizare uniformă 1D (fig. 3.1):

unde .

Dezvoltarea în serie Taylor, în jurul punctului al rețelei, a variabilei necunoscute, este dată de relațiile următoare:

; (3.1)

. (3.2)

Reținând primii doi termeni din dezvoltarea din relația (3.1) se obține:

. (3.3)

Aceasta semnifică că derivata de ordinul 1, în punctul , este aproximată prin diferențe finite regresive de ordinul 1. Reținând primii doi termeni ai dezvoltării din relația (3.2) se obține:

, (3.4)

ceea ce semnifică că derivata de ordinul 1, în punctul , este aproximată prin diferențe finite progresive de ordinul 1.

Scăzând relația (3.1) din relația (3.2) se obține aproximarea, prin diferențe finite centrale de ordinul 2, următoare:

. (3.5)

En additionnant les relations (3.1) et (3.2) on obtient l’approximation de la dérivée de deuxième ordre par différences finies centrées d’ordre 2 :

. (3.6)

Avantajele metodei

simplicitate în aplicare;

utilizare rezonabilă a memoriei (matrice de tip bandă) și timp de calcul, de asemenea, rezonabil.

Dezavantajele metodei

principiul de conservare nu este asigurat după discretizare;

apariția instabilităților numerice;

dificultăți la tratarea geometriilor mai complexe.

3.3 Metoda elementelor finite (MEF)

Principiul matematic

Principiul matematic se bazează pe următoarele metode:

Metode variaționale (minimizarea unei funcționale);

Metoda rezidurilor ponderate.

Principiul fundamental al metodei elementelor finite constă în decuparea domeniului de analiză în domenii elementare de dimensiuni finite. Pe fiecare domeniu elementar, numit element finit, funcția necunoscută este aproximată printr-un polinom al cărui grad poate varia de la o aplicație la alta, dar rămâne în general mic. Aceste elemente finite, triunghiuri sau patrulatere, rectilinii sau curbilinii, trebuie să realizeze o partiție a domeniului de analiză (elementele finite sunt disjuncte și reuniunea lor acoperă întregul domeniu de analiză). Această partiție denumită decupaj sau discretizarea domeniului trebuie să respecte un anumit număr de reguli care permit asigurarea unei bune desfășurări a calculului.

Etapele de aplicare

Discretizarea domeniului într-un număr finit de elemente;

Alegerea modelului (a funcțiilor) de interpolare (variația variabilei pe element) ;

Scrierea ecuațiilor model sub formă algebrică la nivel local (pe element): determinarea vectorilor și matricelor caracteristice;

Asamblarea vectorilor și matricelor locale într-un vector global și o matrice globală ;

Rezolvarea sistemului de ecuații algebrice (fiind matricea necunoscutelor).

Avantaje

Metoda este bine adaptată pentru geometrii complexe (frontiera domeniului este bine aproximată prin segmente de dreaptă care reprezintă laturile elementelor).

Dezavantaje

Formalismul matematic este mai complicat și mai dificil de aplicat;

Costisitoare în memorie de stocare (matricele sunt pline, sau banda matricei are lățime mare) și în timp de calcul (inversarea matricei) ;

Caracterul conservativ al ecuațiilor nu este neapărat asigurat.

3.4 Metode spectrale (MS)

Principiul

Se înlocuiește, în ecuațiile model, necunoscuta prin dezvoltări trunchiate pe bază de funcții ortogonale (polinoame Cebâșev, Legendre, Fourier) și utilizând proprietatea lor de ortogonalitate sunt transformate în sisteme de ecuații diferențiale ordinare care se rezolvă mai simplu.

Avantaje

Permite obținerea unor soluții cu o foarte mare precizie.

Dezavantaje

Formalismul matematic este complex iar aplicarea este dificilă;

Dificultăți la tratarea geometriilor complexe și a condițiilor la limită nestandard.

3.5 Metoda volumelor finite (MVF)

Metoda volumelor finite a fost descrisă pentru prima dată în 1971 de către Patankar și Spalding și publicată în 1980 de Patankar (Numerical Heat Transfer and Fluid Flow, [50]).

Principiul. Metoda volumelor finite este o tehnică de discretizare care convertește ecuațiile de conservare cu derivate parțiale în ecuații algebrice care pot fi rezolvate numeric. Tehnica volumelor de control constă în integrarea ecuațiilor cu derivate parțiale pe fiecare volum de control al domeniului de analiză pentru a obține ecuațiile discretizate ce conservă toate mărimile fizice pe un volum de control (VC).

Principiul de discretizare poate fi ilustrat considerând ecuația de transport pentru o mărime scalară , valabilă pentru toate ecuațiile de curgere, în regim staționar:

, (3.7)

unde

– densitatea fluidului;

– vectorul viteză ();

– vectorul arie al suprafeței;

– coeficientul de difuzie al mărimii ;

– gradientul lui ( în 2D) ;

– termenul sursă (sursa lui pe unitatea de volum).

Ecuația (3.7) este aplicată pe fiecare volum de control al domeniului de calcul (domeniul de studiu sau de analiză). Din discretizarea ecuației (3.7) se obține:

, (3.8)

unde:

– numărul de fețe (interfețe) ale volumului de control;

– valoarea lui transferată prin convecție prin interfața ;

– fluxul de masă prin interfața ;

– aria interfeței ( în 2D) ;

– componenta normală a lui (normală la interfața );

– volumul volumului de control.

Etapele de aplicare sunt următoarele:

Domeniul de calcul este discretizat într-un număr finit de puncte (numite nodurile rețelei de discretizare), în jurul cărora se definesc volumele elementare (numite volume de control) continue, nejuxtapuse și fără discontinuități la interfețe;

Ecuațiile model, sub formă conservativă, sunt integrate pe fiecare volum de control (VC) al domeniului de calcul;

Integralele pe fiecare volum de control atașat unui nod dat sunt evaluate aproximând variația lui cu profile sau legi de interpolare între nodurile vecine nodului considerat;

Scrierea ecuațiilor algebrice pentru valorile lui în nodurile rețelei de discretizare;

Rezolvarea sistemului de ecuații algebrice liniare obținut.

Avantaje

Prezervarea caracterului conservativ al ecuațiilor pe fiecare volum de control (continuitatea fluxului la interfețe), valabilă pentru orice finețe a rețelei de discretizare;

Aplicare relativ ușoară;

Aplicabilă la geometrii complexe;

Timp de calcul și spatiu de memorie pentru stocare rezonabile (matrice de tip bandă).

Dezavantaj

Este mai puțin precisă decât metodele spectrale.

3.5.1 Exemplu: conducția termică 1D staționară

Preupunem ecuația conducției termice 1D în regim staționar:

. (3.7)

Etapele care trebuie parcurse, pentru obținerea sistemului de ecuații algebrice, sunt următoarele:

Discretizarea domeniului de calcul

Modul de discretizare a domeniului de calcul este prezentat în figura 3.2 (a se vedea și anexa C).

unde:

P – Nodul considerat;

W – Nodul “West” ;

E – Nodul “East” ;

W – interfața (fața) “West” a volumului de control (VC);

e – interfața “East” a volumului de control;

– lățimea volumului de control studiat.

Integrarea ecuației de conducție termică 1D pe volumul de control din jurul nodului .

Integrând ecuația (3.7) pe volumul de control (fig. 3.2) se obține succesiv:

; (3.8)

; (3.9)

, (3.10)

unde este valoarea medie a termenului sursă pe volumul de control. Ținând cont de legea lui Fourier (), fiind densitatea fluxului termic, ecuația (3.10) poate fi scrisă astfel:

, (3.11)

unde și sunt densitățile fluxului termic la interfețele volumului de control.

Alegerea unui profil de temperatură (a unei formule de interpolare) între nodurile vecine nodului .

Există două tipuri de profile pe care le putem prevedea, profilul constant (fig. 3.3 a) și profilul liniar (fig.3.3 b).

În cazul unui profil constant de temperatură (fig. 3.3 a), pe volumul de control, avem o discontinuitate a temperaturii la interfețele și ale volumului de control. În plus, derivata nu este definită și deci acest profil de temperatură nu convine.

În cazul unui profil liniar de temperatură (fig. 3.3 b), între nodurile rețelei de discretizare, discontinuitatea temperaturii nu mai există și derivatele la interfețe sunt definite:

; (3.12)

. (3.13)

Scrierea ecuației conducției termice sub formă algebrică

Ecuația (3.11) se scrie astfel:

, (3.14)

unde este valoarea medie a lui pe volumul de control.

În final, după regruparea termenilor, ecuația algebrică sub forma generală următoare:

, (3.15)

unde:

. (3.16)

Observații

Forma generală a ecuațiilor discretizate este următoarea:

. (3.17)

Derivata ar fi putut să fie evaluată și cu alte funcții de interpolare;

Toate mărimile nu trebuie neapărat să fie evaluate cu aceleași funcții de interpolare;

Pentru o aceeași variabilă, pentru diferiți termeni din ecuația model, se pot folosi funcții de interpolare diferite.

Principii ce trebuie respectate

Chiar pe o rețea de discretizare grosieră, trebuie ca:

Variațiile mărimilor să aibă un comportament fizic realist;

Bilanțul global să fie conservativ.

Tratarea termenului sursă

Dacă termenul sursă are o variație neliniară, , se scrie astfel:

, (3.18)

unde este un termen constant (independent de temperatură). Termenul sursă trebuie liniarizat în ppentru a obține un sistem de ecuații algebrice liniare.

3.5.1 Regulile de bază (rgulile lui Patankar)

Următoarele reguli au fost enunțate de Patankar, [50]

Regula Nr. 1: Consistența fluxului la interfețele volumului de control.

Dacă o interfață este comună la două volume de control, expresia fluxului prin aceste interfețe, în ecuațiile discretizate, trebuie să fie aceeași pentru cele două volume de control vecine considerate.

Regula Nr. 2: Toți coeficienții și trebuie să aibă același semn în ecuația discretizată.

Se poate justifica această regulă printr-un contraexemplu. Să presupunem că în ecuația (3.15) am avea >0 , >0 și <0. Atiunci, dacă crește trebuie ca să scadă, deci ar rezulta un comportament fizic nerealist.

Regula Nr. 3 : Termenul sursă liniarizat trebuie să aibă pantă negativă.

Termenul sursă liniarizat trebuie să aibă panta pentru că altfel putem avea < 0 cu > 0 (contrar cu regula nr. 2).

Regula Nr. 4: Ecuațiile discretizate trebuie să rămână valabile atunci când valoarea unei variabile dependente crește cu o valoare constabntă. Matematic regula poate fi scrisă astfel:

;

(3.19)

.

Demonstrație

Dacă , ecuația model este o ecuație diferențială care conține numai derivatele lui ; și sunt soluții ale ecuației diferențiale pentru problema continuă si respectiv discontinuă. Ecuația discretizată pentru cele două soluții se scrie astfel:

; (3.20)

. (3.21)

Scăzând ecuația (3.20) din ecuația (3.21) se obține:

. (3.22)

ceea ce confirmă regula enunțată mai sus.

=== Cap_4_Br ===

4.2.5 Exemple

În continuare se pun în evidență proprietățile schemelor de discretizare explicită și implicită comparând rezultatele numerice, pentru o problemă 1D nestaționară, cu soluția analitică.

Exemplul 1

O placă metalică subțire se găsește inițial la o temperatură uniformă de . La momentul temperatura peretelui “Est” al plăcii este brusc redusă la . Celelalte suprafețe ale plăcii sunt izolate.

Să se utilizeze schema explicită a metodei volumelor finite, pentru un pas de timp adecvat, pentru a calcula distribuția tranzitorie a temperaturii și să se compare rezultatele cu soluția analitică la momentele (i) , (ii) , (iii) ;

Să se calculeze soluția numerică pentru un pas de timp dat de formula (4.90), la momentul și să se compare cu soluția analitică;

Datele problemei sunt: lungimea plăcii , conductivitatea termică și .

Soluție

Ecuația diferențială a conducției termice 1D nestaționar este următoare:

. (4.95)

Condiția inițială este:

Condițiile la limită sunt:

;

.

Soluția analitică este dată de relația următoare [28]:

, (4.96)

unde:

.

Se consideră 6 noduri distribuite uniform pe domeniul de calcul, deci pasul în spațiu este (fig. 4.16).

Pentru un nod interior (nodurile 2, 3 și 4) ecuația discretizată obținută din ecuația (4.87) pentru este următoarea:

, (4.97)

unde:

.

Pentru nodul 1 (nod situat pe frontieră), condiția la limită (temperatura necunoscută) impune obținerea unei ecuații discretizate suplimentare, deci se integrează ecuația (4.95) pe o jumătate de volum de control. Ecuația discretizată astfel obținută, pentru nodul 1, este următoarea:

, (4.98)

unde:

.

Pentru nodul 6 (nod de frontieră), condiția la limită fiind de tip Dirichlet, nu mai este necesară scrierea unei ecuații discretizate suplimentare.

Pentru nodul 5 (nod vecin cu un nod de frontiera pentru care temperatura este cunoscută) se utilizează aceeași ecuație ca pentru un nod interior, dar cum temperatura nodului vecin 6 este cunoscută (C) termenul care conține temperatura nodului 6 trece ca un termen sursă suplimentar. Astfel, se obține ecuația discretizată pentru nodul 5:

, (4.99)

unde:

.

Pasul în timp trebuie să satisfacă condiția de stabilitate (4.90), deci:

.

Pentru că se alege și rezultă:

;

.

După înlocuirea valorilor numerice, în ecuațiile (4.97), (4.98), (4.99) și simplificări, se obține:

(4.100)

În tabelul 4.4 se prezintă un exemplu de calcul, utilizând ecuațiile (4.100) pentru primii doi pași în timp.

Tabelul 4.4

Rezultatele numerice ale exemplului 1 sunt prezentate în tabelul 4.5.

Tabelul 4.5

Tabelul 4.6 prezintă rezultatele numerice și analitice la momentele de timp de , și . În figura 4.17 se prezintă o comparație a rezultatelor numerice obținute (6 noduri), pentru diferite valori ale pasului în timp, cu soluția analitică. Comparația rezultatelor numerice cu soluția analitică, la diferite momente de timp este prezentată în figura 4.18 (pentru 21 noduri).

Tableau 4.6

Exemplul 2

Să se rezolve problema de la exemplul 1 utilizând schema total implicită și să se compare rezultatele numerice obținute cu cele obținute cu schema explicită, pentru un pas de timp de .

Soluție

Se utilizează aceeași rețea de discretizare ca în figura 4.16. Ecuația discretizată, folosind schema total implicită, pentru un nod interior al domeniului de calcul (nodurile 2, 3, și 4) este cea descrisă de ecuația (4.94) dar fără termen sursă, adică:

, (4.101)

unde:

;

.

Pentru nodurile 1 și 5 situate pe frontieră și respectiv vecin cu frontiera se impune o tratare specială. Astfel, pentru nodul 1 ecuația discretizată este’următoarea:

, (4.102)

unde:

.

Pentru nodul 5 ecuația discretizată, ținând cont că (cunoscută), este următoarea:

, (4.103)

unde:

.

Chiar dacă schema implicită permite utilizarea unui pas de timp oarecare, în continuare se utilizează totuși o valoare rezonabilă a acestuia , pentru a asigura o bună precizie a rezultatelor. Se obține astfel, după înlocuirea valorilor numerice:

;

.

După înlocuirea valorilor numerice în ecuațiile (4.101), (4.102), (4.103) și după simplificări, se obține:

(4.104)

Ținând cont că forma matricială a sistemului de ecuații algebrice este:

(4.105)

Se constată că ecuația pentru fiecare nod conține temperaturile necunoscute ale nodurilor vecine. Schema implicită implică rezolvarea simultană a sistemului de ecuații (4.105). Valorile temperaturii de la momentul de timp precedent sunt utilizate numai pentru calculul membrului drept al ecuației matriciale (4.105).

Tabelul 4.7 și figura 4.19 prprezintă rezultatele numerice prin comparație cu soluția analitică (pentru o rețea de discretizare cu 6 noduri). Codul sursă în Fortran este prezentat în anexa G (THER1Di2)

Tableau 4.7

În figura 4.20 se prezintă o comparație a rezultatelor numerice la momentul , obținute la utilizarea schemelor explicită și implicită, cu soluția analitică pentru un pas de timp . Se constată că schema explicită, pentru , dă oscilații, în timp ce schema implicită dă rezultate în acord cu soluția analitică. Această constatare arată avantajul schemei implicite care permite utilizarea unui pas de timp mai mare. Trebuie semnalat totodată că o precizie bună se obține utilizând totuși un pas de timp mai mic.

=== Cap_4_C_i1r ===

4.4 Conducția termică 2D nestaționară

4.4.1 Forma generală a ecuației discretizate

Ecuația conducției termice 2D în regim nestaționar este următoarea:

(4.151)

L’intégration de l’équation (4.151) sur le volume de contrôle schématisé à la figure 4.35 donne

(4.152)

En remplaçant les gradients de température on obtient

(4.153)

En utilisant le schéma totalement implicite, on obtient

(4.154)

En regroupant les termes dans l’équation (4.154) on obtient la forme générale de l’équation discrétisée

(4.155)

où

Tableau 4.9

Si on peut exprimer (linéarisation du terme source) alors les coefficients , de l’équation discrétisée, sont les suivants

Remarque

L’équation (4.155) est valable pour un noeud intérieur du domaine de calcul. Pour les noeuds situés sur la frontière (dans le cas des conditions aux limites de type flux imposé, de type convection ou de type Neumann) ou voisins avec la frontière (dans le cas des conditions aux limites de type Dirichlet) les équations discrétisées sont obtenues en tenant compte de conditions aux limites. Le traitement des différentes conditions aux limites est illustré dans les exemples qui suivent.

4.4.2 Exemples

Exemple 1

On considère une plaque métallique rectangulaire (). La conductivité thermique du matériel de la plaque est . À l’instant la frontière “North” de la plaque est mise en contact avec une paroi maintenue à la température constante . Les autres frontières de la plaque sont maintenues à la température de . Le terme source est constant et reparti de façon uniforme. . À l’instant la plaque se trouve à la température de .

Calculer la distribution transitoire de la température dans la plaque en utilisant le maillage présenté à la figure 4.31 (). On connaît

Solution

L’équation à résoudre est la suivante

(4.156)

L’équation discrétisée pour un noeud intérieur (le noeud 16 par exemple) est la suivante

(4.157)

où

Le terme source étant constant, il n’est pas nécessaire d’être linéarisé. Les valeurs des coefficients voisins, pour le maillage choisi (), sont

Pour les noeuds situés sur la frontière ( les noeuds 1, 2, 3, 4, 5, 6, 7, 12, 13, 18, 19, 24, 25, 26, 27, 28, 29, 30) la température étant connue, il n’est pas nécessaire d’écrire les équations discrétisées. Alors, le nombre d’équations à résoudre est égal à 12, les équations discrétisées sont les suivantes:

(4.158)

En regroupant les inconnues, le système d’équations à résoudre est le suivant

(4.159)

Sous la forme matricielle le système d’équations à résoudre est le suivant:

(4.160)

Pour la résolution du système (4.160) on applique l’algorithme Thomas, adapté aux problèmes 2D, présenté au paragraphe 4.3.3. Pour cela, l’équation discrétisée doit être mise sous la forme

(4.161)

Aux figures 4.36 et 4.38 sont présentées les solutions numériques obtenues à l’aide du programme THERM2Di (Annexe F), aux différents moments de temps, pour un pas de temps .

À la figure 4.38 est présentée une comparaison entre la solution numérique obtenue l’aide du programme THERM2Di, en utilisant la méthode des volumes finis, et celles obtenue à l’aide du logiciel QuickField (QuickField ™ 4.2T version 4.2.2.2, Copyright © 1993-2001 Tera Analysis) en utilisant la méthode des éléments finis.

La comparaison est présentée la longue de la droite . À cause du fait que la version du logiciel QuickField utilisée est limitée à 200 noeuds les erreurs sont assez grandes. Vers la frontière “North” où les gradients de température sont grands il faut utiliser un maillage très dense. On constate que, pour le maillage , les erreurs dans la zone de forts gradients sont très grandes, la distribution de la température ayant un comportement non réaliste. Au contraire si le maillage devient plus dense dans la zone des forts gradients, mais plus grossière dans le reste du domaine, les erreurs diminuent, mais augmente dans la zone . On constate aussi que pour la méthode des volumes finis la solution est presque la même si on passe d’un maillage à 208 noeuds au maillage à 546 noeuds ce qui nous permet de tirer la conclusion que pour un même nombre de noeuds et pour le même pas dans le temps la méthode des volumes finis est plus précise que la méthode des éléments finis dans le cas des problèmes de conduction thermique instationnaire.

En ce qui concerne le pas de discrétisation dans le temps, même si le schéma utilisé est totalement implicite, il faut noter que pour avoir une bonne précision de la solution, le pas dans le temps doit être très petit. Le schéma de discrétisation est du premier ordre comme précision dans le temps, tandis que dans l’espace est du second ordre.

Dans le tableau ci-dessous on présente l’évolution de la distribution de la température pour différents pas de temps à l’instant .

Tableau 4.10

D’après les résultats présentés au tableau 4.1 on constate qu’une bonne précision est obtenue pour un pas de temps .

=== Cap_4_C_i2r ===

Exemplul 2

Se consideră o bară foarte lungă, având secțiunea transversală rectangulară . La momentul secțiunea transversală a barei are o distribuție uniformă de temperatură . Bara este pusă în contact cu un fluid având temperatura uniformă . Coeficientul de transfer termic prin convecție, pe suprafața de separație între bară și fluid, este .

Să se determine evoluția în timp a distribuției câmpului termic, pe secțiunea transversală a barei, utilizând un pas în timp .

Se cunosc următoarele proprietăți ale materialului barei:

conductivitatea termică, ;

densitatea de masă, ;

căldura specifică, .

Soluție

Profitând de simetrie, se poate calcula distribuția câmpului termic pe un sfert din secțiunea transversală (fig. 4.42).

Domeniul de calcul specificat în fig. 4.42 este discretizat și prezentat în fig. 4.43. Se consideră o rețea de discretizare cu 30 noduri ca în exemplul 1 ().

Ecuația diferențială care guvernează regimul tranzitoriu pe secțiunea barei este următoarea:

. (4.162)

Pentru obținerea ecuației discretizate pentru un nod interior al domeniului de calcul (nodul 16 de exemplu) se integrează ecuația (4.162) pe volumul de control hașurat din jurul nodului 16 (fig. 4.43). Utilizând schema total implicită se obține ecuația discretizată (a se vedea obținerea ecuației (4.155)):

, (4.163)

unde:

,

,

,

termenul sursă fiind nul în cazul acestui exemplu.

Valorile numerice ale coeficiențiilor din ecuația (4.163), în condițiile rețelei uniforme folosite, sunt:

;

;

;

.

Ecuația (4.163) se aplică succesiv pentru toate nodurile interioare ale domeniului de calcul (nodurile 8, 9, 10, 11, 14, 15, 16, 17, 20, 21, 22, și 23) și se obțin ecuațiile algebrice următoare:

(4.164)

Pentru un nod interior pe frontiera “West”, nodul 4 de exemplu, se integrează ecuația (4.162) pe o jumătatae de volum de control hașurat și prezentat în figura 4.44.

(4.165)

După integrarea ecuației (4.165) se obține:

Ținând cont că (condiția la limită pe frontiera “West”), presupunând o variație liniară a gradientului de temperatură și utilizând schema total implicită pentru integrarea în timp, se obține:

, (4.166)

unde:

et .

Regrupând termenii în ecuația (4.166), se bține forma generală a ecuației discretizate:

, (4.167)

unde:

;

.

Valorile numerice ale coeficienților pentru ecuația (4.167) sunt următoarele:

;

;

.

Ecuațiile discretizate pentru nodurile interioare ale frontierei “West” (nodurile 2, 3, 4 și 5) sunt:

(4.168)

Pentru un nod interior pe frontiera “South”, nodul 13 de exemplu, se integrează ecuația (4.162) pe jumătate de volum de control, hașurat și prezentat în figura 4.45, și se obține:

.

(4.169)

Ținând cont că (condiție la limită pe frontiera “South”), înlocuind gradienții de temperatură în punctele și și utilizând schema total implicită pentru integrarea în timp, se obține ecuația discretizată sub forma generală:

, (4.170)

unde:

;

;

et .

Valorile coeficienților din ecuația (4.170) sunt:

;

;

.

Ecuațiile discretizate pentru nodurile 7, 13 și 19 sunt următoarele:

(4.171)

Ecuația discretizată pentru nodul 1 se obține prin integrarea ecuației (4.162) pe sfertul de volum de control hașurat și prezentat în figura 4.46. După o primă integrare se obține:

(4.172)

În mod asemănător ca mai sus, se obține ecuația discretizată sub forma generală pentru nodul 1:

, (4.173)

unde:

;

;

et .

Valorile coeficienților din ecuația (4.173) sunt:

;

;

.

Se obține astfel ecuația discretizată pentru nodul 1:

(4.174)

Pentru un nod interior pe frontiera “East”, nodul 28 de exemplu, se integrează ecuația (4.162) pe jumătate de volum de control hașurat și prezentat în figura 4.47.

După o integrare parțială (în spațiu și timp pentru partea din stânga a ecuației și în spațiu pentru partea dreaptă), se obține:

.

(4.175)

Ținând cont că pe această frontieră condiția la limită este

, (4.176)

și integrând în timp, folosind schema total implicită și înlocuind gradienții de temperatură, se obține ecuația următoare:

(4.177)

Regrupând termenii în ecuația (4.177) se obține forma generală a ecuației discretizate:

, (4.178)

unde:

;

;

, et

Valorile coeficienților din ecuația (4.178) sunt:

;

;

;

;

Ecuațiile discretizate pentru nodurile interioare ale frontierei “East” (nodurile 26, 27, 28 și 29) sunt următoarele:

(4.179)

Pentru un nod situat pe frontiera “North”, nodul 18 de exemplu, se integrează ecuația (4.162) pe jumătate de volum de control, hașurat și ilustrat în figura 4.48.

(4.180)

Ținând cont că pe această frontieră condiția la limită este

, (4.181)

și după integrarea ecuației (4.180), se obține forma generală a ecuației discretizate:

, (4.182)

unde:

;

et .

Valorile coeficienților din ecuația (4.182) sunt:

;

;

;

.

Ecuațiile discretizate, pentru nodurile interioare ale frontierei “North” (nodurile 12, 18 și 24 sunt:

(4.183)

Pentru obținerea ecuației discretizate pentru nodul 30 se integrează ecuația (4.162) pe un sfert din volumul de control hașurat și prezentat în figura 4.49, adică

. (4.184)

Efectuând integrarea parțială în ecuația (4.184) se obține:

(4.185)

Ținând cont de condițiile la limită pe frontierele “North” și “East” după integrare în timp, folosind schema total implicită, se obține:

. (4.186)

Regrupând termenii, în ecuația (4.186), se obține forma generală a ecuației discretizate pentru nodul 30:

; (4.187)

unde:

;

;

et .

Valorile coeficienților din ecuația (4.187) sunt:

;

.

Ecuția discretizată pentru nodul 30 este următoarea:

(4.188)

Pentru obținerea ecuației discretizate pentru nodul 6 se integrează ecuația (4.162) pe un sfert de volum de control, hașurat și prezentat în figura 4.50.

. (4.189)

Integrând parțial ecuația (4.189) se obține:

(4.190)

Ținând cont de condițiile la limită pe frontierele “North” și “West” după integrarea în timp (cu schema total implicită) și înlocuirea gradienților de temperatură se obține:

. (4.191)

Regrupând termenii în ecuația (4.191), se obține forma generală a ecuației discretizate pentru nodul 6:

; (4.192)

unde:

;

;

et .

Valorile coeficienților din ecuația (4.192) sunt:

;

;

;

Ecuația discretizată pentru nodul 6 este următoarea:

. (4.193)

Pentru obținerea ecuației discretizate pentru nodul 25 se integrează ecuația (4.162) pe un sfert de volum de control (hașurat și prezentat în figura 4.51).

Integrând ecuația (4.162) ca în cazul nodului 6, se obține forma generală a ecuației discretizate pentru nodul 25:

; (4.194)

unde:

;

;

et .

Valorile coeficienților din ecuația (4.194) sunt:

;

;

.

Ecuația algebrică pentru nodul 25 este următoarea:

. (4.195)

Ecuațiile (4.164), (4.168), (4.171), (4.174), (4.179), (4.183), (4.188), (4.193), et (4.195) formează sistemul de ecuații algebrice ce trebuie rezolvat. Pentru rezolvarea acestui sistem se aplică algoritmul Thomas, adaptat pentru probleme 2D, ca la exemplul 1.În figura 4.52 se prezintă soluția numerică la momentul iar în figura 4.53 la momentul .

În tabelul 4.11 se prezintă comparația soluției numerice, obținută cu metoda volumelor finite (programul THERM2Di_2) cu soluția obținută cu metoda elementelor finite (cu ajutorul programului QuickField 4.2T).

În figura 4.54 se prezintă comparația soluțiilor numerice (volume finite – elemente finite) în ceea ce privește evoluția în timp a temperaturii punctului de coordonate . Se constată o foarte bună corespondentă (eroarea fiind foarte mică).

Tableau 4.11 Comparația rezultatelor

4.5 Conducția termică 3D

4.5.1 Forma generală a ecuației discretizate pentru cazul 3D staționar

Conducția termică în 3D staționar este guvernată de ecuația:

. (4.196)

Volumul de control ce conține nodul are 6 noduri vecine identificate ca nodurile “West”, “East”, “South”, “North”, “Bottom” și “Top” (W, E, S, N, T, B). Ca și în cazul 2D notațiile și fac referință la fețele (interfețele) “west”, “est”, “sud”, “north”, “top” și “bottom”.

Se integrează ecuația (4.196) pe volumul de control cu trei dimensiuni ), prezentat în figura 4.55, astfel:

(4.197)

Dacă se notează , , după integrarea ecuației (4.197) se obține

(4.198)

Aplicând aceeași procedură ecuației (4.198) ca în cazul 2D se obține ecuația discretizată:

(4.199)

Regrupând termenii în ecuația de mai sus se obține forma generală a ecuației discretizate pentru un nod interior al domeniului de calcul:

, (4.200)

unde:

;

;

Condițiile la limită vor fi tratate ca în cazul problemelor 2D.

4.5.2 Forma generală a ecuației discretizate pentru cazul 3D nestaționar

Conducția termică în 3D nestaționar este guvernată de ecuația:

(4.201)

Ecuația discretizată sub forma generală obținută conform procedurii din paragraful 4.5.1, cu schema total implicită, este:

, (4.202)

unde:

;

;

.

Dacă termenul sursă poate fi liniarizat sub forma atunci coeficienții și , din ecuația discretizată, sunt următorii:

;

.

=== Cap_4_C_ir ===

4.4 Conducția termică 2D nestaționară

4.4.1 Forma generală a ecuației discretizată

Ecuația conducției termice 2D în regim nestaționar este

. (4.151)

Integrând ecuația (4.151) pe intervalul de timp cât și pe volumul de control prezentat în figura 4.35 se obține:

(4.152)

Înlocuind gradienții de temperatură, considerând o variație liniară a temperaturii între nodurile rețelei, se obține:

(4.153)

Utilizând o schemă total implicită pentru integrarea în timp se obține:

(4.154)

Regrupând termenii în ecuația (4.154) se obține forma generală a ecuației discretizate:

, (4.155)

unde:

;

;

.

Tableau 4.9

Dacă termenul sursă se poate exprima ca (liniarizarea termenului sursă) atunci coeficienții și , din ecuația discretizată, vor avea forma următoare:

;

.

Remarcă

Ecuația (4.155) este valabilă pentru un nod interior al domeniului de calcul, adică pentru un nod care are toate nodurile vecine în care variabila dependentă este necunoscută. Pentru nodurile situate pe frontieră (în cazul condițiilor la limită de tip flux impus sau de tip convecție) sau vecine cu frontiera (în cazul condițiilor la limită de tip Dirichlet) ecuațiile discretizate sunt obținute printr-o tratare particulară ținând cont de condițiile la limită. Tratarea diferitelor condiții la limită este ilustrată în exemplele care urmează.

4.4.2 Exemple

Exemplul 1

Se consideră o placă metalică rectangulară (). Conductivitatea termică a materialului plăcii este . La momentul frontiera “North” a plăcii este pusă în contact cu un perete menținut la temperatura constantă . Celelalte frontiere ale plăcii sunt menținute la temperatura de . Termenul sursă este constant și distribuit în mod uniform în placă având valoarea . La momentul placa se găsește la temperatura de .

Să se calculeze distribuția tranzitorie a temperaturii în placă, utilizând rețeaua de discretizare prezentată în figura 4.31 (). Se cunoaște .

Soluție

Ecuația ce trebuie rezolvată în acest caz este:

. (4.156)

Ecuația discretizată pentru un nod din interiorul domeniului de calcul (nodul 16 de exemplu) este următoarea:

, (4.157)

unde:

Termenul sursă fiind constant, nu se mai pune problema liniarizării. Valorile coeficienților din ecuația (4.157), pentru rețeaua aleasă (), sunt:

;

;

;

.

Pentru nodurile situate pe frontieră ( nodurile 1, 2, 3, 4, 5, 6, 7, 12, 13, 18, 19, 24, 25, 26, 27, 28, 29, 30) temperatura fiind cunoscută nu este necesară scrierea ecuațiilor discretizate. În aceste condiții, numărul ecuațiilor ce trebuie rezolvate este egal cu 12, sistemul de ecuații fiind următorul:

(4.158)

Regrupând necunoscutele, sistemul de ecuații este:

(4.159)

Sub formă matricială sistemul de ecuații ce trebuie rezolvat este:

(4.160)

Pentru rezolvarea sistemului (4.160) se aplică algoritmul Thomas, adaptat pentru probleme 2D, prezentat la paragraful 4.3.3. Pentru aceasta, ecuația discretizată trebuie pusă sub forma:

. (4.161)

În figurile 4.37 și 4.39 sunt prezentate soluțiile numerice obținute ont présentées les solutions numériques obtenues cu programul THERM2Di (Anexa G), la diferite momente de timp, pentru un pas de timp .

În figura 4.38 se prezintă o comparație între soluția numerică obținută cu ajutorul programului THERM2Di, utilizând metoda volumelor finite, și soluția obținută cu QuickField (QuickField ™ 4.2T version 4.2.2.2, Copyright © 1993-2001 Tera Analysis) utilizând metoda elementelor finite.

Comparația este prezentată pentru rezultatele situate pe dreapta . Din cauza faptului că pachetul de programe QuickField utilizat nu permite mai mult de 200 noduri, erorile sunt mari. În zona frontierei “North” unde gradienții de temperatură sunt mari trebuie să utilizăm o rețea de discretizare foarte fin (fig. 4.40). Se constată că, pentru rețeaua de discretizare din figura 4.41, erorile în zona gradienților mari sunt foarte mari, distribuția temperaturii fiind nerealistă. Dimpotrivă, dacă rețeaua devine mai fină în zona gradienților mari, dar mai grosieră în restul domeniului (fig. 4.40), erorile scad, dar cresc în zona . Se constată de asemenea că pentru metoda volumelor finite soluția este practic aceeași dacă se trece de la o rețea cu 208 noduri la o rețea cu 546 noduri ceea ce ne permite să tragem concluzia că pentru același număr de noduri și același pas în timp metoda volumelor finite este mai precisă decât metoda elementelor finite în cazul problemelor de conducție termică în regim nestaționar. Această analiză este mai degrabă calitativă întrucât o analiză corectă trebuie făcută pentru același tip de rețea de discretizare, pentru același număr de noduri și pentru aceeași schemă de discretizare în timp.

În ceea ce privește pasul de discretizare în timp, chiar dacă schema utilizată este total implicită, trebuie precizat că pentru a avea o bună precizie a soluției numerice, pasul în timp trebuie să fie foarte mic. Schema de discretizare este de ordinul întâi ca precizie după variabila timp, în timp ce după variabilele spațiale este de ordinul doi ca precizie.

În tabelul 4.10 se prezintă evoluția distribuției de temperatură, pentru diferiți pași în timp, la momentul .

Tabelul 4. 10 Evoluția temperaturii pentru diferiți pași în timp.

După rezultatele prezentate în tabelul 4.1 se constată că o precizie bună se obține pentru un pas în timp .

=== Cap_4_Cr ===

4.3 Conducția termică 2D staționară

4.3.1 Forma generală a ecuației discretizate

Metodologia utilizată pentru discretizarea ecuației diferențiale în cazul unidimensional poate fi cu ușurință aplicătă și în cazul bidimensional (2D). Pentru a ilustra această tehnică, se consideră ecuația conducției termice 2D în regim staționar

. (4.106)

Rețeaua de discretizare folosită în cazul 2D este schematizată în figura 4.21

În plus, față de rețeaua de discretizare 1D, la nodurile vecine, “East” () și “West” () ale nodului (punctului) se adaugă vecinii “North” () și “South” ().

Integrând ecuația (4.106), pe volumul de control hașurat din figura 4.21 se obține:

(4.107)

Dacă se notează și , se obține

(4.108)

Ecuația (4.108) reprezintă bilanțul între sursa lui , în volumul de control, și fluxurile pe fețele volumului de control. Utilizând aceeași aproximație ca în cazul 1D, adică se presupune o variație liniară a gradientului de temperatură între două noduri vecine ale rețelei de discretizare, se pot scrie fluxurile pe fețele volumului de control:

(4.109)

(4.110)

(4.111)

(4.112)

Înlocuind expresiile fluxurilor pe fețele volumului de control, din relațiile (4.109) la (4.112), în ecuația (4.108), se obține

(4.113)

Dacă se ține cont că și și în final regrupând termenii, ecuația (4.113) se scrie astfel:

(4.114)

Ecuația discretizată (4.114) poate fi scrisă sub forma generală, pentru un nod interior, astfel:

(4.115)

unde:

.

Pentru a obține distribuția temperaturii (sau pentru orice altă variabilă dependentă ) pentru o problemă 2D se scrie ecuația discretizată pentru fiecare nod al rețelei de discretizare. Dacă pe o porțiune din frontiera domeniului de analiză temperatura este cunoscută, ecuațiile discretizate pentru nodurile vecine ecuațiile discretizate se modifică corespunzător, ca în cazul 1D, iar pe porțiunile de frontieră pentru care se cunoaște densitatea fluxului termic sau există un transfer termic prin convecție, ecuațiile discretizate sunt modificate corespunzător pentru a ține cont de aceste condiții la limită. Aceste particularități ale tratării nodurilor situate pe frontieră vor fi prezentate în exemplele ce urmează.

4.3.2 Exemple

Exemplul 1

Se consideră o placă dreptunghiulară () având grosimea (fig. 4.22).. Conductivitatea termică a materialului plăcii este . Frontiera “West” a plăcii primește un flux constant iar frontierele “South” și “East” sunt izolate. Frontiera “North” este menținută la temperatura constantă de .

Să se determine distribuția staționară a temperaturii, calculând temperatura în nodurile rețelei de discretizare (1, 2, 3, …etc.) prezentată în figura 4.22 ().

Soluție

Ecuația care guvernează transferul termic, în condițiile enunțate, este următoarea:

(4.116)

Ecuația discretizată pentru un nod interior (nodul 16 de exemplu), conform cu (4.115), este:

(4.117)

unde:

Valorile coeficienților ecuației discretizate pentru nodurile vecine cu nodul , în condițiile unei rețele de discretizare uniforme, sunt:

Ecuațiile discretizate pentru nodurile interioare ale domeniului de calcul (nodurile 8-11, 14-17 și 20-23) sunt

(4.118)

Pentru nodurile situate pe frontiera “West” (nodurile 2, 3, 4 și 5) ecuația discretizată se obține integrând ecuația conducției termice (4.116) pe jumătate de volum de control hașurat și prezentat în figura 4.23.

(4.119)

Prin integrarea ecuației (4.119) se obține

(4.120)

unde și . Înlocuind gradienții de temperatură în punctele , și se obține

unde este fluxul impus (cunoscut) pe frontiera “West”.

Regrupând termenii, se obține ecuația discretizată pentru nodurile interioare ale frontierei “West” sub forma

(4.121)

unde:

Valorile numerice ale coeficienților ecuației (4.121) sunt următoarele:

Ecuațiile discretizate pentru nodurile 2, 3, 4 și 5 sunt următoarele:

(4.122)

Pentru nodurile situate pe frontiera “East” (nodurile 26, 27, 28 și 29) ecuația discretizată se obține integrând ecuația conducției termice (4.116) pe jumătate de volum de control hașurat și prezentat în figura 4.24

(4.123)

(4.124)

unde și . Ținând cont că și presupunând o variație liniară a gradientului de temperatură, se obține

(4.125)

Regrupând termenii ecuației (4.125) se obține ecuația discretizată, pentru un nod interior de pe frontiera “East”, sub forma

(4.126)

unde:

Valorile numerice ale coeficienților ecuației (4.126) sunt următoarele:

Ecuațiile discretizate pentru nodurile 26, 27, 28, și 29 sunt următoarele

(4.127)

Pentru nodurile situate pe frontiera “South” (nodurile 7, 13, și 19) ecuația discretizată se obține integrând ecuația conducției termice (4.116) pe jumătate de volum de control, hașurat și prezentat în figura 4.25.

Prin integrarea ecuației (4.116) pe jumătate de volum de control (fig. 4.25) se obține

, (4.128)

unde și . Aproximând gradienții de temperatură printr-o variație liniară și ținând cont de condiția la limită se obține

. (4.129)

Regrupând termenii în ecuația (4.129) se obține în final forma generală a ecuației discretizate:

(4.130)

unde:

Valorile numerice ale coeficienților ecuației (4.130) sunt următoarele:

Ecuațiile discretizate pentru nodurile 7, 13 și 19 sunt următoarele:

(4.131)

Nodurile 1 și 25 sunt, de asemenea, tratate într-un mod special. Astfel, pentru nodul 1 se integrează ecuația (4.116) pe sfertul de volum de control hașurat (fig. 4.26).

(4.132)

(4.133)

unde și . Ținând cont că și , ecuația (4.133) se reduce la ecuația

(4.134)

Regrupând termenii, se obține forma generală a ecuației discretizate pentru nodul 1.

(4.135)

unde:

Valorile numerice ale coeficienților ecuației (4.135) sunt următoarele:

Ecuația ce trebuie rezolvată este

(4.136)

Pentru a obține ecuația discretizată pentru nodul 25 (nod de colț “East – South”) se integrează ecuația (4.116) pe un sfert de volum de control (fig. 4.27 ).

(4.137)

După integrare se obține ecuația

(4.138)

unde și . Ținând cont că și , se obține

(4.139)

Regrupând termenii se obține ecuația discretizată pentru nodul 25.

(4.140)

unde:

Valorile numerice ale coeficienților ecuației (4.140) sunt următoarele:

Ecuația ce trebuie rezolvată este

(4.141)

În final, sistemul de ecuații ce trebuie rezolvat este următorul:

(4.142)

Soluția numerică a sistemului (4.142) este prezentată în tabelul 4.8. Se constată o bună corespondență, pentru aproape același număr de noduri, între rezultatele obținute cu metoda volumelor finite și cele obținute cu metoda elementelor finite utilizând pachetul software QuickField (Student’s QuickField TM 3.40a, Copyright 1993-97 Tera Analysis). Pentru o mai bună apreciere a preciziei, trebuie comparate rezultatele numerice cu soluția analitică, dacă aceasta se poate obține în anumite condiții. În figura 4.28 sunt prezentate curbele izoterme obținute cu QuickField.

Tabelul 4.8 Soluția numerică a sistemului (4.142) (metoda volumelor finite – MVF, 30 noduri și 546 noduri, și metoda elementelor finite – MEF, 31 noduri și 485 noduri)

Analizând rezultatele din tabelul 4.8 se constată o foarte bună corespondență între rezultatele obținute cu metoda volumelor finite (546 noduri) și cele obținute cu metoda elementelor finite (485 noeuds). Întrucât matricea coeficienților sistemului de ecuații algebrice este o matrice de tip bandă, se poate utiliza algoritmul Thomas adaptat cazului bidimensional.

4.3.3 Aplicarea algoritmului Thomas la probleme 2D (algoritmul linie cu linie)

Algoritmul Thomas poate fi aplicat iterativ pentru a rezolva sistemul de ecuații algebrice în cazul problemelor 2D. Se combină metoda directă a algoritmului Thomas pe direcția unei coordonate și metoda iterativă Gauss-Seidel pe direcția celeilalte coordonate. Fie exemplul rețelei de discretizare din figura 4.29

Ecuațiile discretizate pe linia , de exemplu, sunt următoarele:

(4.143)

Dacă , se rezolvă probleme 1D pe direcția “” cu TDMA, baleind toți indicii . Temperaturile sunt ultimele valori calculate ale temperaturii (iterația precedentă) pe liniile vecine cu linia pentru care se aplică TDMA.

Observații importante

Convergența metodei linie cu linie este mai rapidă decât convergența metodei Gauss-Seidel întrucât informația de pe frontierele (condițiile la limită) liniei pentru care se aplică TDMA este transmisă instantaneu către interiorul domeniului de calcul și aceasta oricare ar fi numărul de noduri pe linie.

Viteza de transmitere a informației pe cealaltă direcție este aceeași ca pentru metoda Gauss-Seidel.

Alternând direcțiile de aplicare a TDMA, se poate accelera propagarea informației conținute pe toate frontierele către interiorul domeniului de calcul.

Geometria și rețeaua de discretizare pot furniza coeficienți mult mai mari într-o direcție față de cealaltă direcție. În acest caz trebuie aleasă ca direcție de aplicare a TDMA direcția cu coeficienții mai mari, întrucât aceasta va conduce la accelerarea convergenței.

Fie exemplul din figura 4.30. Coeficienții ecuației discretizate sunt următorii:

unde .

Aceasta implică și deci trebuie luată ca direcție de aplicarea a TDMA direcția coordonatei “”. Ecuația discretizată este următoarea:

(4.144)

În anumite cazuri, direcția de baleiere pentru iterațiile Gauss-Seidel este impusă de către condițiile la limită sau de către direcția de curgere a fluidului (din amonte către aval)

Exemplul 2

Se consideră o placă metalică rectangulară (). Conductivitatea termică a materialului plăcii este . Toate frontierele plăcii sunt menținute la temperatura de iar termenul sursă, distribuit uniform, este .

Să se determine distribuția staționară a temperaturii în placă, utilizân rețeaua de discretizare din figura 4.31 ().

Soluție

Ecuația ce guvernează distribuția temperaturii în placă este

(4.145)

Ecuația discretizată pentru un nod interior (nodul 16 de exemplu) este următoarea

(4.146)

unde:

Termenul sursă fiind constant, nu se mai pune problema liniarizării. Valorile coeficienților ecuației (4.146), pentru rețeaua de discretizare aleasă, sunt:

Pentru nodurile situate pe frontiera domeniului de calcul ( nodurile 1, 2, 3, 4, 5, 6, 7, 12, 13, 18, 19, 24, 25, 26, 27, 28, 29, 30) temperatura fiind cunoscută nu este necesară scrierea ecuațiilor discretizate. În aceste condiții, numărul de ecuații ale sistemului ce trebuie rezolvat este 12. Ecuațiile discretizate sunt următoarele:

(4.147)

Regrupând necunoscutele în ecuații, sistemul de ecuații ce trebuie rezolvat, sub formă matricială este următorul:

(4.148)

Soluția sistemului de ecuații (4.148) este:

(4.149)

Soluția analitică a problemei pentru domeniul hașurat din figura (fig. 4.34) este [28]:

(4.150)

În figura 4.34 se prezintă soluția analitică obținută cu ajutorul expresiei (4.150) valabilă pentru domeniul hașurat din figura 4.33. În figura 4.34-a se prezintă soluția numerică. Se constată o foarte bună concordanță între soluția numerică și soluția analitică (dar trebuie ținut cont că în acest caz chiar soluția analitică este rezultatul unui calcul numeric de evaluare a unei serii trunchiate). În tabelul de mai jos se prezintă o comparație a rezultatelor obținute cu diferite metode (cu volume finite folosind programul THERM2D2 (anexa G), cu elemente finite folosind QuickField și analitic). Eroarea este calculată în raport cu soluția analitică.

=== Cap_4_Dr ===

4.6 Aplicarea metodei volumelor finite în coordonate cilindrice

4.6.1 Forma generală a ecuației discretizate 1D în regim staționar

Metoda volumelor finite se poate aplica pentru orice tip de rețea de discretizare ortogonală în același mod ca pentru rețeaua de discretizare carteziană. Pentru a ilustra tehnica de aplicare în cazul coordonatelor cilindrice, considerăm ca exemplu ecuația de conservare a energiei (Elenbaas-Heller) valabilă pentru secțiunea transversală a plasmei de arc electric. Ecuația diferențială, în coordonate cilindrice, este următoarea:

, (4.203)

unde este conductivitatea electrică, este intensitatea câmpului electric iar este puterea radiată () de către coloana de arc electric.

Se integrează ecuația (4.203) pe volumul de control asociat nodului “” (fig. 4.56), astfel:

(4.204)

unde .

Înlocuind expresia lui în ecuația (4.204) se obține:

. (4.205)

Integrarea (4.205) poate fi scrisă explicit, după simplificări, astfel:

. (4.206)

După integrarea ecuației (4.206) se obține:

, (4.207)

sau

. (4.208)

Regrupând termenii în ecuația (4.208) se obține forma tridiagonală a ecuației discretizate, valabilă pentru un nod interior:

. (4.209)

4.6.2 Condiții la limită

Condiție la limită de tip Dirichlet în nodul “1”.

Dacă se pune în ecuația (4.209) se obține:

. (4.210)

Ținând cont că temperatura est cunoscută și deci termenul care o conține trece ca termen sursă, se obține:

.

(4.211)

Ecuația (4.211) reprezintă forma generală tridiagonală a ecuației discretizate pentru nodul “” care ține cont că temperatura nodului “”,, este cunoscută.

Condiție la limită de tip Dirichlet în nodul “N”.

Dacă se pune în ecuația (4.209) se obține:

(4.212)

Ținând cont că temperatura este cunoscută și deci termenul care o conține trece ca termen sursă, se obține:

(4.213)

Ecuația (4.213) reprezintă forma generală tridiagonală a ecuației discretizate pentru nodul “” care ține cont că temperatura nodului “”,, este cunoscută.

Condiție la limită de tip Neumann în nodul “1”.

Se integrează ecuația (4.204) pe o jumătate de volum de control, hașurat și prezentat în figura (4.57), astfel:

. (4.214)

Din integrarea ecuației (4.214) se obține:

. (4.215)

Ținând cont de condiția la limită în nodul “” care se exprimă prin , și prin urmare al doilea termen din ecuația (4.215) se anuleazăest nul, se obține:

. (4.216)

Regrupând termenii în ecuația (4.216) se obține ecuația discretizată, sub formă tridiagonală, pentru nodul “”.

(4.217)

Condiție la limită de tip Neumann în nodul “N”.

Se integrează ecuația (4.204) pe jumătatea de volum de control, hașurat și ilustrat în figura 4.58, astfel:

. (4.218)

După integrarea ecuației (4.218) se obține:

. (4.219)

Ținând cont de condiția la limită în nodul “” exprimată prin , primul termen al ecuației (4.219) se anulează, și se obține ecuația discretizată, sub forma tridiagonală, pentru nodul “”:

. (4.220)

4.6.3 Aplicație numerică

Fie domeniul de analiză 1D având , , și rețeaua de discretizare prezentată în figura 4.59. Condițiile la limită sunt de tip Neumann în nodul și Dirichlet în nodul 5 ().

Ecuația discretizată pentru nodul 1, ținând cont că conductivitatea termică și (a se vedea ecuația (4.217)), devine:

. (4.221)

Ecuațiile discretizate pentru nodurile 2 și 3 devin (a se vedea ecuația (4.209)):

. (4.222)

. (4.223)

Ecuația doscretizată pentru nodul 4, ținând cont că temperatura nodului 5 () este cunoscută, devine:

. (4.224)

Înlocuind valorile numerice în ecuațiile (4.221), (4.222), (4.223) și (4.224) se obține sistemul de ecuații algebrice ce trebuie rezolvat:

(4.225)

sau, sub formă matricială:

(4.226)

Soluția sistemului algebric (4.226) este:

, (4.227)

ceea ce corespunde exact cu valorile calculate (fig. 4.60) cu ajutorul soluției analitice exprimată prin relația următoare:

. (4.228)

În cazul ecuației Elenbaas-Helller, conductivitatea termică este funcție de temperatură, , conductivitatea electrică este de asemenea funcție de temperatură, , și puterea radiată este funcție de temperatură, . Aceste dependențe de temperatură ale lui , și sunt prezentate în figurile 4.61, 4.62 și 4.63 pentru cazul unei plasma de arc în hexafluorură de sulf (SF6), la presiunea de [53]. În aceste condiții ecuația diferențială este puternic neliniară și algoritmul de rezolvare trebuie să cuprindă etapele prezentate la paragraful 4.1.3.

În figura 4.64 sunt prezentate soluțiile numerice, pentru diferite valori ale intensității câmpului electric. Soluția analitică neexistând în cazul neliniar, trebuie verificat dacă codul numeric furnizează soluția corectă sau nu. Pentru aceasta, fie se compară soluția numerică cu date experimentale, fie se apelează la o validare numerică a rezultatelor numerice.

4.6.4 Validare numerică

Dacă în codul numeric se impun valori constante pentru conductivitatea termică și termenul sursă , atunci ecuația (4.204) are o soluție analitică unică dată de relația (4.228). Dacă rezultatele numerice, obținute pentru valori constante ale mărimilor și corespund (făcând abstracție de erorile de calcul) cu soluția analitică, se spune că s-a realizat validarea numerică a codului numeric (programul de calcul) și deci și a rezultatelor pentru cazul neliniar. În figura 4.65 este prezentat rezultatul validării codului numeric pentru două valori ale termenului sursă și pentru o valoare constantă a conductivității termice (). Se constată o corespondență perfectă între soluția numerică și soluția analitică.

Observații

1. Din cauza faptului că ecuația este puternic neliniară, există riscul ca procesul iterativ să devină instabil. Nu trebuie să uităm că la fiecare iterație rezolvarea se face cu coeficienții ecuației calculați folosind soluția iterației precedente. Pentru a evita, în anumite cazuri, instabilitățile numerice, este necesară utilizarea relaxării. Fie, de exemplu, valoarea temperaturii calculate la iterația curentă și valoarea calculată la iterația precedentă. Noua valoare folosită la iterația următoare poate fi exprimată astfel:

, (4.229)

unde este un factor de relaxare.

2. Un proces iterativ converge atunci când iterațiile următoare nu mai produc schimbări semnificative asupra valorii necunoscutei. În practică, procesul iterativ este oprit când criteriul de convergență impus este satisfăcut. Alegerea unui criteriu de convergență adecvat depinde de natura problemei rezolvate precum și de obiectivele calculului numeric. Cel mai adesea, în timpul procesului iterativ, se examinează mărimea cea mai semnificativă a problemei (temperatura, tensiunea de arc, etc.). Se continuă iterațiile până când nu se mai constată o schimbare esențială, a valorii acestei mărimi, între două iterații succesive (de exemplu, cea mai mare schimbare a acestei mărimi, dintre toate nodurile domeniului de calcul, să fie mai mică sau egală cu o valoare impusă).

Când se utilizează subrelaxarea, schimbarea valorii necunoscutei între două iterații succesive este în mod forțat redusă ceea ce creiază iluzia de convergență. Cea mai bună soluție de monitorizare a convergenței este de a examina satisfacerea ecuației discretizate de către valorile curente ale necunoscutei. Astfel pentru fiecare nod al rețelei de discretizare se calculează un rezidu cu relația:

. (4.230)

Valoarea lui devine nulă atunci când ecuația discretizată este satisfăcută. Criteriul de convergență poate fi exprimat, de exemplu, impunând ca valoarea lui să fie mai mică sau egală cu o valoare fixată.

4.6.5 Forma generală a ecuației discretizate 2D în regim staționar

Ecuația conducției termice 2D în regim staționar este:

. (4.231)

Ecuația discretizată se obține integrând ecuația diferențială (4.231) pe volumul de control hașurat și prezentat în figura 4.66.

(4.232)

unde și care înlocuit în ecuația (4.232) se obține:

. (4.233)

Din integrarea de mai sus rezultă:

, (4.234)

sau

.

(4.235)

Înlocuind gradienții de temperatura, în ecuația (4.235), se obține:

,

(4.236)

unde .

Regrupând termenii în ecuația (4.236) se obține forma generală a ecuației discretizate:

, (4.237)

unde:

;

;

;

4.6.6 Forma generală a ecuației discretizate 2D în regim nestaționar

Ecuația diferențială ce guvernează conducția termică 2D nestaționară este:

. (4.238)

Se integrează ecuația (4.238) pe volumul de control din figura 4.47.

, (4.239)

unde . Introducând în ecuația de mai sus expresia lui se obține:

(4.240)

Integrând în spațiu și în timp numai membrul stâng al ecuației (4.240) se obține:

(4.241)

Presupunând o variație liniară a gradienților de temperatură, ecuația (4.241) devine:

(4.242)

Utilizând schema total implicită de integrare în timp se obține ecuația

(4.243)

Regrupând termenii în ecuația (4.243) se obține forma generală a ecuației discretizate:

, (4.244)

unde:

;

;

.

=== Cap_5_Ar ===

Capitolul 5

Metoda volumelor finite aplicată la probleme de convecție/difuzie

5.1 Introducere

Problemele de convecție/difuzie apar în cazul curgerii fluidelor (gaze sau lichide). Problema care se pune este de a găsi câmpul termic (distribuția temperaturii) în prezența unor câmpuri de viteze și densități date. Termenul de difuzie este utilizat aici într-un sens generalizat. Fluxul de difuzie generat de către gradientul variabilei generale este care poate fi un flux de difuzie al unor specii chimice sau un flux termic de exemplu. Pentru problemele în care curgerea fluidului joacă un rol semnificativ trebuie luate în seamă efectele convecției. Ecuația de convecție/difuzie în regim staționar poate fi obținută din ecuația generală de conservare (ecuația de transport) (2.9) neglijând termenul tranzitoriu, adică

. (5.1)

Din integrarea formală pe volumul de control, a ecuației (5.1), se obține:

. (5.2)

Ecuația (5.2) reprezintă bilanțul fluxului într-un volum de control. Membrul stâng al ecuației (5.2) reprezintă fluxul de convecție net, în timp ce primul și al doilea termen al membrului drept reprezintă fluxul de difuzie net și respectiv fluxul de generare a lui (termenul sursă) în volumul de control.

Principala problemă la discretizarea termenului convectiv este calculul variabilei la interfețele volumului de control și a fluxului de convecție prin aceste interfețe (frontiere). În capitolul 4 s-a utilizat schema cu diferențe centrale la obținerea ecuației discretizate pentru termenul de difuzie și termenul sursă. Fenomenul de difuzie afectează distribuția unei mărimi (cantitate) transportate pe direcția gradienților săi în toate direcțiile, în timp ce convecția are influență numai pe direcția de curegere. Prin urmare, schema cu diferențe centrale nu poate fi utilizată pentru tratarea convecției. Pentru a asigura astfel stabilitatea calculului numeric, trebuie să apelăm la alte scheme de discretizare pentru termenul de convecție.

Fie un câmp de curgere care verifică ecuația de continuitate:

. (5.3)

Ecuația diferențială generală sub forma conservativă

, (5.4)

poate fi rescrisă, ținând cont de ecuația de continuitate, astfel:

. (5.5)

Pentru această formă a ecuației (neconservativă) se constată că, pentru distribuții date ale lui , , , și , o anumită soluție a lui verifică ecuația (5.5).

5.2 Ecuația de convecție/difuzie 1D în regim staționar

Pentru a ilustra modul de tratare al termenului convectiv, se consideră ecuația de convecție/difuzie 1D în regim staționar:

. (5.6)

Curgerea trebuie să satisfacă ecuația de continuitate (conservarea masei):

. (5.7)

Se consideră de asemenea volumul de control cu interfețele situate la jumătatea distanței dintre noduri. În continuare se prezintă diferite schema de discretizare.

5.2.1 Discretizare cu schema centrată (schema cu diferențe centrale)

Integrând, pe volumul de control (fig. 5. 1), ecuația (5.6) se obține ecuația:

. (5.8)

După integrarea ecuației (5.8) se obține:

. (5.9)

Aplicând o interpolare liniară și ținând cont că interfețele sunt situate la jumătatea distanței dintre noduri, ecuația (5.9) devine:

, (5.10)

unde:

.

Este convenabil să se definească două variabile și care reprezintă fluxul masic convectiv pe unitatea de suprafață și respectiv conductanța de difuzie, astfel:

, (5.11)

unde:

dacă sau dacă și

Se obține astfel forma generală a ecuației discretizate după schema centrată:

, (5.12)

unde:

;

.

Observații

Ecuația de continuitate implică , deci și .

Interpolarea liniară pentru a obține corespunde diferențelor centrale ale dezvoltărilor în serie Taylor, deci o precizie de ordinul 2.

Pentru că și , și pot fi negativi dacă . În aceste condiții regula No 2 nu este obligatoriu verificată și rezultatele sunt dezastruoase, adică nu se respectă principiul fizic al fenomenului.

Exemplu

Fie și atunci se obțin pentru coeficienți valorile următoare:

.

Cazul 1: Dacă , atunci !

Cazul 2: Dacă , atunci !

Dacă unul din coeficienți (ceilalți coeficienți fiind pozitivi) atunci avem:

Criteriul lui Scarborough [50] nu se verifică dacă există relația de mai sus. În aceste condiții algoritmul Gauss-Seidel poate diverge (schema centrată este deci limitată la valorile mici ale numărului ().

Dacă atunci rezultă:

,

și sistemul liniar nu este rezolvabil prin metoda Gauss-Seidel și de asemenea prin majoritatea metodelor iterative.

Concluzie: Trebuie găsită o altă schemă de discretizare decât schema centrată.

Discretizare cu schema “upwind” (schema descentrată în amonte)

Schema “upwind” est construită pentru a evita de a spune că dacă convecția este puternică pe o direcție, informația într-un punct poate să vină din aval (cum este cazul la schema centrată).

Formularea schemei “upwind”

Formularea schemei pentru termenul de difuzie este neschimbată, dar termenul de convecție se calculează pe baza următoarelor ipoteze:

dacă ; (5.13)

dacă ; (5.14)

dacă ; (5.15)

dacă ; (5.16)

Dacă se introduce operatorul următor:

, (5.17)

atunci schema descentrată se scrie astfel:

; (5.18)

. (5.19)

Utilizând conceptul descris mai sus, se obține ecuația de convecție difuzie discretizată sub forma:

, (5.20)

unde:

;

;

.

Observații

Coeficienții și nu sunt niciodată negativi;

Soluțiile sunt totdeauna acceptabile din punct de vedere fizic;

Criteriul lui Scarborough [50] este verificat.

Soluția exactă a problemei de convecție/difuzie 1D staționar

Ecuația (5.6) poate fi rezolvată în mod exact dacă coeficientul este presupus constant ( este de asemenea constant în conformitate cu ecuația de continuitate (5.7)). În aceste condiții ecuația ce trebuie rezolvată este:

. (5.21)

Dacă se utilizează domeniul cu condițiile la limită

;

; (5.22)

Soluția analitică a ecuației (5.21) cu condițiile la limită (5.22) este:

, (5.23)

unde este numărul lui Péclet (care reprezintă raportul între forțele de convecție și difuzie) definit de relația:

(5.24)

Observații

Dacă (nu există convecție) avem difuzie pură și se obține un profil liniar pentru pentru .

Dacă avem influență preponderentă a lui .

Dacă avem o influență preponderentă a lui (valori negative ale lui ).

Schema centrată dă , ceea ce este exact pentru și acceptabil pentru valorile mici ale lui .

Schema descentrată, adică,

dacă

dacă

este exactă pentru foarte mare și neexactă dacă este mic.

Atunci când este mare, , deci difuzia este nulă. Însă, în schema “upwind”, se calculează termenul de difuzie folosind un profil liniar pentru între și ,adică la valori mari ale numărului Péclet, se supraestimează difuzia (fenomenul poartă numele de difuzie numerică).

Schema exponențială

Ecuația (5.6) poate fi rescrisă astfel:

sau , (5.25)

unde:

.

Forma discretizată din ecuația (5.25) se scrie:

. (5.26)

Pentru a calcula fluxurile și se folosește soluția exactă între și pentru , și între și pentru , adică

, (5.27)

unde:

și ,

, (5.28)

unde:

și .

Se găsesc astfel expresiile fluxurilor în punctele și :

; (5.29)

. (5.30)

Remarcă: Fluxurile și nu deepind de poziția interfețelor și respectiv .

Înlocuind expresiile (5.29) și (5.30) în ecuația discretizată (5.26) se obține forma generală obișnuită a ecuației discretizate:

, (5.31)

unde: .

Observații

Pentru o problemă 1D în regim staționar această schemă permite obținerea soluției exacte pentru orice valoare a numărului Péclet și pentru orice număr de noduri.

Totuși această schemă este destul de puțin utilizată întrucât:

Exponențialele se calculează numeric mai greu, deci timp de calcul mai mare ;

Schema exponențială nu este exactă pentru probleme 2D și 3D, nestaționare și cu termen sursă.

Schema hibridă

Această schemă are, din punct de vedere calitativ, comportamentul schemei exponențiale dar este mai ieftină din punct de vedere al timpului de calcul. La schema exponențială avem:

. (5.32)

Dacă se înlocuiește se obține forma următoare pentru ecuația (5.32):

. (5.33)

În figura 5.4 se reprezintă variația raportului în funcție de numărul Péclet și se constată că:

Dacă , adică punctul (nodul) al rețelei de discretizare este în aval de nodul, influența nodului scade când valoarea lui crește. La limită când atunci raportul ;

Dacă , adică nodul al rețelei de discretizare este în amonte de nodul , influența nodului crește când valoarea lui crește. La limită când , raportul ;

Dacă , tangenta la curba soluției exacte este .

Concluzii

Schema hibridă este o combinație a acestor trei tangente (a se vedea figura 5.4) astfel:

Pentru ,

(5.34)

Pentru ,

(5.35)

Pentru ,

(5.36)

Sub o formă compactă schema hibridă (relațiile (5.34), (5.35) și (5.36)) se poate scrie astfel:

. (5.37)

Schema hibriddă este identică cu schema centrată pentru :

Schema hibridă este identică cu schema “upwind” pentru , dar cu adică nu se mai supraestimează difuzia;

Schema hibridă este o combinație ameliorată între schema centrată și schema “upwind”.

Ecuația de convecție/difuzie 1D staționar se scrie:

, (5.38)

unde:

.

Schema “Power Law” (Legea Putere)

Se constată pe figura 5.4 că ecartul între soluția exactă (sau schema exponențială) și schema hibridă este maxim când . Pe de altă parte este exagerat de a considera difuzia nulă pentru . Schema “Power Law” realizează o mai bună aproximare a curbei exacte (exponențială) și în plus nu este costisitoare, din punct de vedere al timpului de calcul, în comparație cu schema exponențială.

Pentru schema “Power Law” coeficienții sunt definiți astfel:

Pentru ,

(5.39)

Pentru ,

(5.40)

Pentru ,

(5.41)

Pentru ,

. (5.42)

Sub o formă compactă coeficienții pot fi scriși astfel:

. (5.43)

Observații

Pentru , schema “Power Law” este identică cu schema hibridă;

Schemele “Power Law” și exponențială sunt atât de apropiate încât nu pot fi comparate pe grafic;

Schema “Power Law” este recomandată pentru problemele de convecție/difuzie.

Recapitularea diferitelor scheme

Ecuația discretizată sub forma generală este:

, (5.44)

unde:

.

Coeficienții și , în funcție de schema utilizată, sunt prezentați în tabelul 5.1

Tabloul 5.1 Coeficienții și pentru diferite scheme.

Observații

Schema centrată este valabilă (recomandată) pentru ;

Schema “upwind” este neadecvată pentru valori mici ale lui iar difuzia este supraestimată;

Schema exponențială are un cost (din punct de vedere numeric) ridicat;

Schema hibridă este acceptabilă dar aceasta furnizează erori maxime pentru ;

Schema “Power Law” este cea recomandată.

Formularea generală a diferitelor scheme

Se consideră rețeaua de discretizare din figura 5.5. Fluxul total ce traversează interfața între cele două noduri și poate fi scris, ținând cont de relația (5.23), astfel, definind:

, (5.45)

unde este numărul Péclet, .

Dacă considerăm că fluxul la interfață este o medie ponderată între fluxurile și ,

, (5.46)

și că gradientul lui este proporțional cu ,

(5.47)

atunci fluxul total poate fi scris astfel:

(5.48)

Fluxul total poate fi scris sub forma generală:

, (5.49)

unde și sunt coeficienți fără dimensiune care sunt funcții de .

Proprietățile lui și

Dacă , din relația (5.47), se poate trage concluzia că fluxul de difuzie este nul și deci:

. (5.50)

Din relația (6.46) se obține:

. (5.51)

Combinând relațiile (5.50) și (5.51) se obține:

. (5.52)

Dacă devine , atunci și schimbă rolurile, adică:

sau . (5.53)

Consecințele proprietăților

Din schema exponențială se obține (a se vedea relația (5.29)):

. (5.54)

Fluxul total este deci:

(5.55)

În figura 5.6 sunt reprezentate variațiile exacte ale funcțiilor și .

Funcțiile și sunt cunoscute pentru din moment ce este conoscută pentru . În fapt pentru avem:

(5.56)

Pentru orice valoare a lui , adică se poate scrie:

; (5.57)

. (5.58)

Ținând cont că

, (5.59)

ecuația de convecție/difuzie 1D staționar discretizată (5.26) devine:

(5.60)

Ținând cont de relația (5.49), ecuația (5.60) devine:

. (5.61)

Regrupând termenii, se obține:

, (5.62)

sau

. (5.63)

Coeficienții ecuației (5.63), ținând cont de relațiile (5.57) și (5.58), sunt următorii:

; (5.64)

. (5.65)

Toate schemele pot fi scrise sub forma unică (5.63) cu coeficienții dați de relațiile (5.64) și (5.65) în care coeficienții , ppentru diferite scheme, sunt obținuți prin identificare între relațiile de mai sus și relațiile din tabelul 5.1. Acești coeficienți sunt preyentați în tabelul 5.2.

Tabelul 5.2

În figura 5.7 se preyintă funcția pentru diferite scheme.

Observații

Toate schemele, în afara schemei centrate pentru , dau o soluție acceptabilă din punct de vedere fizic;

Pentru că numărul Péclet decide asupra comportării schemelor, s-ar putea folosi o rețea de discretizare densă ( mic) pentru a avea și a putea astfel alege schema centrată. Totuși pentru majoritatea problemelor fizice aceasta ar necesita un numar mare de volume de control și deci costuri de calcul foarte ridicate.

Scheme cu diferențe de ordin superior pentru probleme de convecție/difuzie

Erorile mai pot fi încă minimizate prin utilizarea schemelor de discretizare de ordin superior. Aceste scheme implică utilizarea mai multor puncte vecine. Schema cu diferențe centrale, care este de ordinul doi ca precizie, poate fi instabilă pentru că această schemă nu ține cont de direcția de curgere. Schemele de ordin superior trebuie să conserve proprietatea schemei și sensibilitatea la direcția de curgere pentru a asigura stabilitatea.

Schema QUICK

Această schemă folosește o interpolare cuadratică ponderată cu trei puncte pentru a calcula valoarea mărimii căutate la iinterfața volumului de control. Valoarea lui la interfață este obținută cu ajutorul unei funcții cuadratice (funcție de gradul doi) care trece prin două noduri ale rețelei de discretizare și un nod situat în amonte de interfață (fig. 5.8).

De exemplu, când și un profil cuadratic care trece prin punctele , și este utilizat pentru a calcula și un alt profil cuadratic care trece prin , și este utilizat pentru a calcula .

Pentru și valorile lui în punctele , și sunt utilizate pentru a calcula în timp ce valorile în punctele , și pentru a calcula . Se poate arăta că, pentru o rețea de discretizare uniformă, (a se vedea Anexa A), valoarea lui la interfața situată între două noduri vecine și este dată de formula următoare:

. (5.66)

Când , nodurile vecine pentru interfața sunt și iar nodul din amonte este nodul (fig.5.8). Valoarea lui la interfață este:

. (5.67)

Când , nodurile vecine pentru interfața sunt și iar nodul din amonte este nodul (fig.5.8). Valoarea lui la interfață este:

. (5.68)

Il est important de noter que pour les termes de diffusion on obtient les mêmes expressions, dans l’équation discrétisée, que pour le schéma aux différences centrales. Si on considère l’équation de convection-diffusion 1D stationnaire discrétisée sous la forme :

. (5.69)

Si et et si on utilise les formules (5.67) et (5.68) pour les termes de convection et les différences centrales pour les termes de diffusion, l’équation (5.69) s’écrit :

.

(5.70)

En regroupant les termes on obtient :

. (5.71)

Sous la forme générale l’équation s’écrit ainsi :

, (5.72)

unde:

;

.

Dans le cas où et les valeurs de aux interfaces et sont données par les expressions :

; (5.73)

. (5.74)

La substitution des expressions (5.73) et (5.74) dans les termes de convection dans l’équation discrétisée (5.69) et l’utilisation des différences centrales pour les termes de diffusion, on obtient la forme générale de l’équation discrétisée :

, (5.75)

unde:

;

.

On peut obtenir des expressions générales, qui sont valables tant pour la direction d’écoulement positive que pour la direction négative, par la combinaison des équations (5.72) et (5.75) ainsi :

, (5.76)

avec le coefficient central :

,

et les coefficients voisins :

;

,

unde: ;

.

=== Cap_5_Ar1 ===

Capitolul 5

Metoda volumelor finite aplicată la probleme de convecție/difuzie

5.1 Introducere

Problemele de convecție/difuzie apar în cazul curgerii fluidelor (gaze sau lichide). Problema care se pune este de a găsi câmpul termic (distribuția temperaturii) în prezența unor câmpuri de viteze și densități date. Termenul de difuzie este utilizat aici într-un sens generalizat. Fluxul de difuzie generat de către gradientul variabilei generale este care poate fi un flux de difuzie al unor specii chimice sau un flux termic de exemplu. Pentru problemele în care curgerea fluidului joacă un rol semnificativ trebuie luate în seamă efectele convecției. Ecuația de convecție/difuzie în regim staționar poate fi obținută din ecuația generală de conservare (ecuația de transport) (2.9) neglijând termenul tranzitoriu, adică

. (5.1)

Din integrarea formală pe volumul de control, a ecuației (5.1), se obține:

. (5.2)

Ecuația (5.2) reprezintă bilanțul fluxului într-un volum de control. Membrul stâng al ecuației (5.2) reprezintă fluxul de convecție net, în timp ce primul și al doilea termen al membrului drept reprezintă fluxul de difuzie net și respectiv fluxul de generare a lui (termenul sursă) în volumul de control.

Principala problemă la discretizarea termenului convectiv este calculul variabilei la interfețele volumului de control și a fluxului de convecție prin aceste interfețe (frontiere). În capitolul 4 s-a utilizat schema cu diferențe centrale la obținerea ecuației discretizate pentru termenul de difuzie și termenul sursă. Fenomenul de difuzie afectează distribuția unei mărimi (cantitate) transportate pe direcția gradienților săi în toate direcțiile, în timp ce convecția are influență numai pe direcția de curegere. Prin urmare, schema cu diferențe centrale nu poate fi utilizată pentru tratarea convecției. Pentru a asigura astfel stabilitatea calculului numeric, trebuie să apelăm la alte scheme de discretizare pentru termenul de convecție.

Fie un câmp de curgere care verifică ecuația de continuitate:

. (5.3)

Ecuația diferențială generală sub forma conservativă

, (5.4)

poate fi rescrisă, ținând cont de ecuația de continuitate, astfel:

. (5.5)

Pentru această formă a ecuației (neconservativă) se constată că, pentru distribuții date ale lui , , , și , o anumită soluție a lui verifică ecuația (5.5).

5.2 Ecuația de convecție/difuzie 1D în regim staționar

Pentru a ilustra modul de tratare al termenului convectiv, se consideră ecuația de convecție/difuzie 1D în regim staționar:

. (5.6)

Curgerea trebuie să satisfacă ecuația de continuitate (conservarea masei):

. (5.7)

Se consideră de asemenea volumul de control cu interfețele situate la jumătatea distanței dintre noduri. În continuare se prezintă diferite schema de discretizare.

5.2.1 Discretizare cu schema centrată (schema cu diferențe centrale)

Integrând, pe volumul de control (fig. 5. 1), ecuația (5.6) se obține ecuația:

. (5.8)

După integrarea ecuației (5.8) se obține:

. (5.9)

Aplicând o interpolare liniară și ținând cont că interfețele sunt situate la jumătatea distanței dintre noduri, ecuația (5.9) devine:

, (5.10)

unde:

.

Este convenabil să se definească două variabile și care reprezintă fluxul masic convectiv pe unitatea de suprafață și respectiv conductanța de difuzie, astfel:

, (5.11)

unde:

dacă sau dacă și

Se obține astfel forma generală a ecuației discretizate după schema centrată:

, (5.12)

unde:

;

.

Observații

Ecuația de continuitate implică , deci și .

Interpolarea liniară pentru a obține corespunde diferențelor centrale ale dezvoltărilor în serie Taylor, deci o precizie de ordinul 2.

Pentru că și , și pot fi negativi dacă . În aceste condiții regula No 2 nu este obligatoriu verificată și rezultatele sunt dezastruoase, adică nu se respectă principiul fizic al fenomenului.

Exemplu

Fie și atunci se obțin pentru coeficienți valorile următoare:

.

Cazul 1: Dacă , atunci !

Cazul 2: Dacă , atunci !

Dacă unul din coeficienți (ceilalți coeficienți fiind pozitivi) atunci avem:

Criteriul lui Scarborough [50] nu se verifică dacă există relația de mai sus. În aceste condiții algoritmul Gauss-Seidel poate diverge (schema centrată este deci limitată la valorile mici ale numărului ().

Dacă atunci rezultă:

,

și sistemul liniar nu este rezolvabil prin metoda Gauss-Seidel și de asemenea prin majoritatea metodelor iterative.

Concluzie: Trebuie găsită o altă schemă de discretizare decât schema centrată.

Discretizare cu schema “upwind” (schema descentrată în amonte)

Schema “upwind” est construită pentru a evita de a spune că dacă convecția este puternică pe o direcție, informația într-un punct poate să vină din aval (cum este cazul la schema centrată).

Formularea schemei “upwind”

Formularea schemei pentru termenul de difuzie este neschimbată, dar termenul de convecție se calculează pe baza următoarelor ipoteze:

dacă ; (5.13)

dacă ; (5.14)

dacă ; (5.15)

dacă ; (5.16)

Dacă se introduce operatorul următor:

, (5.17)

atunci schema descentrată se scrie astfel:

; (5.18)

. (5.19)

Utilizând conceptul descris mai sus, se obține ecuația de convecție difuzie discretizată sub forma:

, (5.20)

unde:

;

;

.

Observații

Coeficienții și nu sunt niciodată negativi;

Soluțiile sunt totdeauna acceptabile din punct de vedere fizic;

Criteriul lui Scarborough [50] este verificat.

Soluția exactă a problemei de convecție/difuzie 1D staționar

Ecuația (5.6) poate fi rezolvată în mod exact dacă coeficientul este presupus constant ( este de asemenea constant în conformitate cu ecuația de continuitate (5.7)). În aceste condiții ecuația ce trebuie rezolvată este:

. (5.21)

Dacă se utilizează domeniul cu condițiile la limită

;

; (5.22)

Soluția analitică a ecuației (5.21) cu condițiile la limită (5.22) este:

, (5.23)

unde este numărul lui Péclet (care reprezintă raportul între forțele de convecție și difuzie) definit de relația:

(5.24)

Observații

Dacă (nu există convecție) avem difuzie pură și se obține un profil liniar pentru pentru .

Dacă avem influență preponderentă a lui .

Dacă avem o influență preponderentă a lui (valori negative ale lui ).

Schema centrată dă , ceea ce este exact pentru și acceptabil pentru valorile mici ale lui .

Schema descentrată, adică,

dacă

dacă

este exactă pentru foarte mare și neexactă dacă este mic.

Atunci când este mare, , deci difuzia este nulă. Însă, în schema “upwind”, se calculează termenul de difuzie folosind un profil liniar pentru între și ,adică la valori mari ale numărului Péclet, se supraestimează difuzia (fenomenul poartă numele de difuzie numerică).

Schema exponențială

Ecuația (5.6) poate fi rescrisă astfel:

sau , (5.25)

unde:

.

Forma discretizată din ecuația (5.25) se scrie:

. (5.26)

Pentru a calcula fluxurile și se folosește soluția exactă între și pentru , și între și pentru , adică

, (5.27)

unde:

și ,

, (5.28)

unde:

și .

Se găsesc astfel expresiile fluxurilor în punctele și :

; (5.29)

. (5.30)

Remarcă: Fluxurile și nu deepind de poziția interfețelor și respectiv .

Înlocuind expresiile (5.29) și (5.30) în ecuația discretizată (5.26) se obține forma generală obișnuită a ecuației discretizate:

, (5.31)

unde: .

Observații

Pentru o problemă 1D în regim staționar această schemă permite obținerea soluției exacte pentru orice valoare a numărului Péclet și pentru orice număr de noduri.

Totuși această schemă este destul de puțin utilizată întrucât:

Exponențialele se calculează numeric mai greu, deci timp de calcul mai mare ;

Schema exponențială nu este exactă pentru probleme 2D și 3D, nestaționare și cu termen sursă.

Schema hibridă

Această schemă are, din punct de vedere calitativ, comportamentul schemei exponențiale dar este mai ieftină din punct de vedere al timpului de calcul. La schema exponențială avem:

. (5.32)

Dacă se înlocuiește se obține forma următoare pentru ecuația (5.32):

. (5.33)

În figura 5.4 se reprezintă variația raportului în funcție de numărul Péclet și se constată că:

Dacă , adică punctul (nodul) al rețelei de discretizare este în aval de nodul, influența nodului scade când valoarea lui crește. La limită când atunci raportul ;

Dacă , adică nodul al rețelei de discretizare este în amonte de nodul , influența nodului crește când valoarea lui crește. La limită când , raportul ;

Dacă , tangenta la curba soluției exacte este .

Concluzii

Schema hibridă este o combinație a acestor trei tangente (a se vedea figura 5.4) astfel:

Pentru ,

(5.34)

Pentru ,

(5.35)

Pentru ,

(5.36)

Sub o formă compactă schema hibridă (relațiile (5.34), (5.35) și (5.36)) se poate scrie astfel:

. (5.37)

Schema hibriddă este identică cu schema centrată pentru :

Schema hibridă este identică cu schema “upwind” pentru , dar cu adică nu se mai supraestimează difuzia;

Schema hibridă este o combinație ameliorată între schema centrată și schema “upwind”.

Ecuația de convecție/difuzie 1D staționar se scrie:

, (5.38)

unde:

.

Schema “Power Law” (Legea Putere)

Se constată pe figura 5.4 că ecartul între soluția exactă (sau schema exponențială) și schema hibridă este maxim când . Pe de altă parte este exagerat de a considera difuzia nulă pentru . Schema “Power Law” realizează o mai bună aproximare a curbei exacte (exponențială) și în plus nu este costisitoare, din punct de vedere al timpului de calcul, în comparație cu schema exponențială.

Pentru schema “Power Law” coeficienții sunt definiți astfel:

Pentru ,

(5.39)

Pentru ,

(5.40)

Pentru ,

(5.41)

Pentru ,

. (5.42)

Sub o formă compactă coeficienții pot fi scriși astfel:

. (5.43)

Observații

Pentru , schema “Power Law” este identică cu schema hibridă;

Schemele “Power Law” și exponențială sunt atât de apropiate încât nu pot fi comparate pe grafic;

Schema “Power Law” este recomandată pentru problemele de convecție/difuzie.

Recapitularea diferitelor scheme

Ecuația discretizată sub forma generală este:

, (5.44)

=== Cap_5_B+C+Dr ===

5.2.10 Exemple

Exemplul 5.1

O anumită proprietate (mărime fizică) , de exemplu temperatura, este transportată prin convecție și difuzie prin intermediul domeniului unidimensional ilustrat în figura 5.9. Condițiile la limită sunt următoarele:

;

.

Utilizând o rețea de discretizare cu 6 noduri și schema cu diferențe centrale, să se calculeze distribuția mărimii în funcție de pentru:

Cazul 1 : și să se compare cu soluția analitică;

Cazul 2 : și să se compare cu soluția analitică;

Cazul 3 : Să se recalculeze soluția numerică pentru utilizând 21 noduri și să se compare cu soluția analitică.

Se cunoaște: , , .

Soluția analitică este:

.

Soluție

Se utilizează rețeaua de discretizare prezentată în figura 5.10 pentru care . Trebuie amintit că , , , .

Ecuația care guvernează transportul mărimii este:

(5.77)

Pentru un nod interior (nodurile 3 și 4) ecuația discretizată poate fi scrisă sub forma următoare:

; (5.78)

.

Pentru nodurile vecine cu nodurile de pe frontiere (nodurile 2 și 5) nu este necesară scrierea unor ecuații discretizate suplimentare. Pentru aceste noduri se utilizează aceeași ecuație ca pentru un nod interior, dar pentru că valorile lui sunt cunoscute în nodurile 1 și 6 (, ) termenii corespunzători acestor noduri trec ca termeni sursă. Astfel, pentru nodul 2 ecuația discretizată este:

. (5.79)

Pentru nodul 5 ecuația discretizată devine:

. (5.80)

Cazul 1

; ;

;

.

Ecuațiile ce trebuie rezolvate sunt următoarele:

Nodul 2 :

Nodul 3 :

Nodul 4 :

Nodul 5 :

Forma matricială a sistemului de ecuații algebrice ce trebuie rezolvat este:

(5.81)

Soluția sistemului (5.81) este:

(5.82)

În figura 5.11 se prezintă comparația între soluția numerică și soluția analitică.

Se constată o foarte bună corespondență între soluția numerică și soluția analitică.

Cazul 2

; ; ;

;

.

Ecuațiile discretizate ce trebuie rezolvate sunt:

Noeud 2 :

Noeud 3 :

Noeud 4 :

Noeud 5 :

Forma matricială a sistemului de ecuații algebrice este:

(5.83)

Soluția sistemului (5.83) este:

(5.84)

În figura 5.12 se prezintă, pentru comparație, soluțiile analitică și numerică. Se constată o instabilitate numerică la schema cu diferențe centrale.

Cazul 3

; ; ;

;

.

Ecuațiile discretizate ce trebuie rezolvate sunt următoarele:

Nodul 2 :

Nodul 3-19 :

Nodul 20 :

În figura 5.13 se prezintă soluțiile numerică și analitică.

În acest caz (cazul 3) se constată o bună corespondență între soluția numerică și soluția analitică (cazul 3). Comparația rezultatelor cazurilor 2 și 3 arată că o rețea de discretizare mai fină (cazul 3) a determinat scăderea numărului Péclet (raportul ) de la 5 la 1,25. Schema cu difrențe centrale furnizează rezultate bune când raportul este mic.

Exemplul 5.2

Să se rezolve problema de la exemplul 1 utilizând schema de discretizare “upwind” pentru cazurile:

;

,

și aceeași rețea de discretizare (cu 6 noduri). Să se compare rezultatele numerice cu soluția nalaitică.

Soluție

Pentru un nod interior (nodurile 3 și 4) ecuația discretizată (5.20) poate fi scrisă sub forma:

, (5.85)

unde:

;

.

Pentru nodurile 2 și 5 sunt valabile ecuațiile (5.79) și (5.80).

Cazul 1

; ; ; ;

; ; .

Sistemul de ecuații algebrice ce trebuie rezolvat este:

Nodul 2:

Nodul 3:

Nodul 4:

Nodul 5:

Forma matricială a sistemului de ecuații algebrice este:

(5.86)

Soluția sistemului (5.86) este:

(5.87)

În tabelul 5.3 se prezintă pentru comparație soluțiile numerică, analitică și erorile.

Tabelul 5.3

Cazul 2

; ; ; ;

; ; .

Sistemul de ecuații algebrice ce trebuie rezolvat este:

Nodul 2:

Nodul 3:

Nodul 4:

Nodul :

Forma matricială a sistemului de ecuații algebrice ce trebuie rezolvat este:

(5.88)

Soluția sistemului (5.88) este:

(5.89)

Tabelul 5.4

Rezultatele numerice, în comparație cu soluția analitică, sunt prezentate în figura 5.15 și în tabelul 5.4.

Remarcă

Se constată că, cu aceeași rețea de discretizare (6 noduri), schema “upwind” dă o soluție realistă în raport cu schema cu diferențe centrale chiar dacă în cazul schemei “upwind” eroarea este destul de mare în vecinătatea frontierei a domeniului de calcul.

Exemplul 5.3

Să se rezolve problema considerată la Cazul 2 de la exemplul 5.1 utilisând schema hibridă pentru . Să se compare soluția unei rețele cu 6 noduri cu soluția unei rețele cu 21 noduri.

Soluție

; ;

;

;

;

;

.

Sistemul de ecuații algebrice ce trebuie rezolvat în acest caz este:

Nodul 2 :

Nodul 3 :

Nodul 4 :

Nodul 5 :

Forma matriceală a sistemului de mai sus este:

(5.90)

Soluția sistemului (5.90) este:

(5.91)

Rezultatele numerice pentru o rețea cu 21 noduri sunt prezentate în tabelul 5.5 și în figura 5.16. Se constată că chiar dacă numărul Péclet este mare (> 2) erorile sunt aproape aceleași ca în cazul schemei “upwind”. Dacă rețeaua devine mai fină, numărul Péclet devine mai mic decât 2, schema hibridă se transformă în schema cu diferențe centrale și soluția este destul de precisă.

Tabelul 5.5

5.2.11 Proprietățile schemelor de discretizare

Până aici s-a prezentat o înlocuire mecanică a derivatelor prin expresii discrete fără a face nici un fel de considerație suplimentară. Dar înainte de a admite reprezentativitatea soluției printr-o schemă cu diferențe, trebuie mai întâi să ne asigurăm că schema reprezintă în mod corect ecuația originală (ecuația diferențială), apoi să fim siguri că erorile, chiar foarte mici, nu se amplifică, și în final trebuie să verificăm că soluția numerică obținută reprezintă soluția problemei.

O bună înțelegere și interpretare a soluției unui algoritm numeric este foarte importantă. În general, există trei concepte matematice [8] utilizate la caracterizarea algoritmilor numerici: consistența, stabilitatea și convergența.

Proprietatea de consistență

Fie o ecuație cu derivate parțiale reprezentată simbolic prin relațiile:

; (5.92)

,

unde este variabila dependentă, fiind domeniul iar frontiera sa. Formele discrete ale acestor relații, adică cele obținute prin înlocuirea derivatelor prin aproximații, vor fi notate astfel:

; (5.93)

.

Definiție

Eroare de trunchiere se definește prin expresia:

. (5.94)

Această eroare este o eroare sistematică care este introdusă de fiecarea dată când se aproximează un operator continu printr-un operator discret.

Definiție

O schemă numerică este consistentă cu o ecuație cu derivate parțiale, dacă eroarea de trunchiere tinde la zero atunci când pasul de discretizare tinde la zero. Ideia de consistență o putem imagina ca un fel de fidelitate locală între ecuația discretă și ecuația cu derivate parțiale. Dacă problema continuă nu este aproximată corect, nu se va putea obține convergența schemei. Calitatea consistenței nu este altceva decât precizia schemei de discretizare.

Proprietatea de stabilitate

Consistența este o garanție locală între discretizare și ecuația cu derivate parțiale. Putem însă să considerăm că ultimul pas, , se face în aceleași condiții ca și primul pas, ?. Diferența este evidentă întrucât la ultimul pas se regăsește acumularea erorilor de discretizare a pașilor precedenți. Această realitate impune o condiție globală, vizând ca la ultimii pași, soluția numerică să rămână aproape de soluția exactă. Altfel spus, pentru ca o soluție să aibă sens, trebuie ca erorile de calcul inevitabile să nu se amplifice. Acest aspect conduce la noțiunea de stabilitate.

Definiție

Un algoritm numeric este stabil dacă soluția problemei discrete rămâne mărginită.

Această definiție implică următoarele observații:

Stabilitatea este un concept care se aplică la algoritmii numerici și nu la ecuațiile diferențiale valabile pe domenii continue.

Definiția este similară noțiunii de stabilitate din mecanică.

Definiția nu impune nici o restricție la oscilații.

Proprietatea de convergență

În pofida existenței noțiunilor de consistență și de stabilitate, nimic nu garantează reprezentativitatea soluției. Aceată caracteristică este dată de către definiția convergenței.

Definiție

Un algoritm numerică este convergent dacă diferența între soluția exactă și soluția numerică tinde către zero când pasul de discretizare tinde la zero.

Evident că diferența nu este calculabilă, totuși teorema următoare ne permite să ieșim din impas.

Teorema lui Lax

Pentru o problemă liniară, consistența și stabilitatea sunt necesare și suficiente pentru a asigura convergența.

Trebuie amintit că aceasta înseamnă că chiar dacă un algoritm numeric convergent este consistent și stabil, contrarul nu este valabil (un algoritm numeric consistent și stabil nu este totdeauna convergent). Totuși, teorema fundamentală a lui Lax asigură că stabilitatea și consistența implică convergența. În cazul problemelor neliniare consistența și stabilitatea sunt condiții necesare pentru convergență dar ele nu sunt suficiente. Pentru astfel de probleme Patankar [50] a formulat anumite reguli (prezentate deja în capitolul 4) care fac din schemele de calcul numeric, cu metoda volumelor finite, scheme numerice de calcul robuste. Metodele de calcul robuste sunt caracterizate de către trei proprietăți fundamentale: proprietatea de conservativitate, proprietatea de mărginire și proprietatea de transportativitate.

Proprietatea de conservativitate

După integrarea ecuației de convecție/difuzie pe un număr finit de volume de control se obține un sistem de ecuații discretizate care implică fluxurile mărimii de transport prin interfețele volumelor de control. Pentru a asigura conservarea lui pe tot domeniul, fluxul lui care iese printr-o anumită interfață trebuie să fie egal cu fluxul lui care intră prin aceeași interfață a volumului de control vecin. Fluxul ce traversează o interfață comună trebuie să fie reprezentat consistent, adică prin aceeași expresie pe cele două volume de control vecine. De exemplu, să considerăm problema de difuzie 1D staționară fără termen sursă (fig. 5.17).

Fluxurile la frontierele domeniului sunt notate cu și . Se consideră 5 noduri cu volumele de control atașate și se aplică schema cu diferențe centrale pentru a calcula fluxul de difuzie ce traversează interfețele. Fluxul care intră prin interfața “” a volumului de control din jurul nodului 3 este iar fluxul care iese prin interfața “” a aceluiași volum de control este . Bilanțul global al fluxului poate fi obținut adunând fluxul net (diferența între fluxul de ieșire și fluxul de intrare) al tuturor volumelor de control și ținând cont de condițiile pe frontierele domeniului, astfel:

(5.95)

Dacă se ține cont că , , , , se constată că fluxurile de difuzie sunt exprimate consistent pentru că numai fluxurile frontierelor rămân în bilanțul global al domeniului. Această formulare, cu diferențe centrale, asigură conservarea fluxului de difuzie al mărimii pe întregul domeniu de calcul.

Utilizarea inconsistentă a formulelor de interpolare poate să genereze un algoritm numeric care nu asigură conservarea variabilei dependente pe întregul domeniu de calcul. De exemplu, se consideră situația în care se utilizează funcția cvadratică 1 prin valorile în nodurile 1, 2 și 3 pentru volumul de control 2 și o altă funcție cvadratică 2 prin valorile în nodurile 2, 3 și 4 pentru volumul de control 3 (fig. 5.18).

Valorile fluxurilor calculate pe interfața “” a volumului de control 2 și pe interfața “” a volumului de control 3 pot fi diferite dacă gradienții funcțiilor, în punctul acestei interfețe, sunt diferite. În aceste condiții conservarea globală a fluxului nu este satisfăcută. Acest exemplu nu sugerează că interpolarea cvadratică ar fi neindicată. Se va arăta în continuare că interpolarea cvadratică numită și schema QUICK este consistentă.

Proprietatea de mărginire

Ecuațiile discretizate pentru fiecare nod al rețelei reprezintă un sistem de ecuații care trebuie ezolvat. De obicei, un asemenea sistem de ecuații se rezolvă utilizând tehnici iterative. Procesul iterativ începe cu o distribuție inițială arbitrară a variabilei și apoi soluția inițială este actualizată după fiecare iterație până la obținerea convergenței. Scarborough (1958) [50] a arătat că o condiție suficientă de convergență pentru o metodă iterativă poate fi exprimată în funcție de valorile coeficienților ecuațiilor discretizate:

, (5.96)

unde este coeficientul efectiv al nodului central (mai exact, ). Dacă schema de discretizare produce coeficienți care satisfac criteriul lui Scarborough, matricea coeficienților este dominant de tip diagonală. Pentru aceasta coeficientul efectiv trebuie să aibă valore mare. Se constată deci, că liniarizarea termenului sursă cu negativ asigură o valoare mare pentru coeficientul pentru că este pozitiv. Dominanța de tip diagonală este o garanție pentru satisfacerea proprietății de mărginire a variabilei .

Pentru problema de conducție termică staționară 1D, fără termen sursă (exemplul 1, paragraful 4.1.7), condițiile la limită, adică temperatura pe frontiere este de și respectiv . Toate valorile temperaturii nodurilor interioare sunt cuprinse între și . O altă condiție pentru ca variabila să fie mărginită este ca toții coeficienții sistemului de ecuații discretizate să aibă același semn (regula nr. 2). Din punct de vedere fizic aceasta inseamnă că o creștere a valorii lui într-un nod va determina o creștere a valorii lui în toate nodurile vecine. Dacă schema de discretizare nu satisface condiția de mărginire pentru , este posibil ca soluția sa nu conveargă pentru toate nodurile, și dacă converge este posibil să aibă vârfuri în soluție (a se vedea fig.5.12). Acest fenomen este bine ilustrat de către rezultatele cazului 2 de la exemplul 5.1.(toți coeficienții sunt negativi, , a se vedea deasemenea matricea coeficienților sistemului (5.83)).

Proprietatea de transportativitate

Proprietatea de transportativitate a variabilei în cazul curgerii poate fi ilustrată considerând o sursă constantă pentru în punctul (fig. 5.19).

Se definește un coeficient fără unitate de măsură ca fiind raportul între forțele de convecție și forțele de difuzie:

. (5.97)

În figura 5.19 sunt prezentate formele generale ale izoliniilor lui pentru diferite valori ale numărului Péclet. Se consideră două cazuri extreme pentru a putea identifica întinderea influenței nodului situat în amonte asupra nodului situat în aval :

nu există convecție, deci există difuzie pură () ;

nu există difuzie, deci există convecție pură ().

În cazul difuziei pure, fluidul nu curge () și izovalorile lui sunt cercuri concentrice cu centrele în punctul . Procesele de difuzie tind să propage pe în toate direcțiile. Condițiile din punctul vor fi influențate atât de fenomenele ce au loc în punctul cât și de fenomenele care au loc în aval de punctul . Odata cu creșterea numărului Péclet, contururile izovalorilor se modifică de la forma circulară până la forma eliptică. Aceste forme ale contururilor sunt deformate în direcția de curgere. Condițiile din punctul vor fi puternic influențate de către condițiile din punctul , în timp ce condițiile din punctul vor fi slab influențate sau chiar deloc de către condițiile din punctul . În cazul unei convecții pure () contururile pentru constant vor fi puternic deplasate în direcția de curgere. Toate proprietățile termenului sursă al lui din punctul vor fi transportate foarte repede în punctul . Dacă nu există difuzie, este egal cu .

5.2.12 Proprietățile schemei cu diferențe centrale

Proprietatea de conservativitate

Schema cu diferențe centrale utilizează expresii consistente pentru evaluarea fluxului convectiv la interfețele volumului de control. Comentariile făcute la paragraful 5.2.10 arată că schema cu diferențe centrale este conservativă.

Proprietatea de mărginire

Coeficienții, pentru nodurile interioare, ai ecuației discretizate (5.12) sunt cei din tabelul de mai jos.

Curgerea staționară unidimensională satisface ecuația de continuitate. Această ecuație discretizată implică relația . În aceste condiții, coeficientul poate fi exprimat de relația . Se constată că coeficienții schemei cu diferențe centrale satisfac criteriul lui Scarborough. Cu contribuția convectivă la coeficientul interfeței “” este negativă; dacă convecția este dominantă atunci coeficientul poate fi negativ. Știind că și (curgerea este unidirecțională), atunci pentru ca să fie pozitiv trebuie ca și să satisfacă condiția următoare:

. (5.98)

Dacă coeficientul va fi negativ și una din exigențele proprietății de mărginire este violată. În exemplul 5.2, cazul 2 avem și deci condiția (5.65) este violată. Consecința asupra soluției numerice poate fi constatată în figura 5.15.

Proprietatea de transportativitate

Schema cu diferențe centrale introduce influențe către toți vecinii nodului atunci când se calculează fluxul de convecție și de difuzie. Schema nu este sensibilă (nu recunoaște) direcția de curgere.

Precizia

Eroarea de trunchiere, pentru dezvoltarea în serie Taylor, în cazul schemei cu diferențe centrale, este de ordinul doi (pentru mai multe detalii, a se vedea anexa A). Exigența de a avea coeficienți pozitivi în cazul schemei cu diferențe centrale este exprimată prin formula (5.65) ceea ce este o garanție că schema va fi stabilă iar precizia va fi, de asemenea acceptabilă. Este important să semnalăm că numărul lui Péclet definit de relația (5.64) este o combinație a proprietăților fluidului ( și ), a proprietății de curgere () și a proprietății rețelei de discretizare (pasul ). Pentru valori date ale lui și , este posibilă satisfacerea relației (5.65) dacă viteza fluidului este mică sau pasul rețelei de discretizare este mic. Ținând cont de aceste limitări, schema cu diferențe centrale nu este recomandată să fie utilizată la calculul curgerilor.

5.2.13 Proprietățile schemei “upwind”

Proprietatea de conservativitate

Schema “upwind” utilizează expresii consistente pentru a calcula fluxul ce traversează interfețele volumelor de control.

Proprietatea de mărginire

Coeficienții ecuației discretizate sunt pozitivi și astfel se satisface exigența de mărginire pentru variabila în cadrul schemei “upwind”. Matricea coeficienților este dominantă de tip diagonal iar “vârfurile” nu mai există în soluție.

Proprietatea de transportativitate

Această proprietate este inclusă în formularea acestei scheme.

Precizia

Schema fiind bazată pe diferențe regresive, precizia este de ordinul unu în ceea ce privește eroarea de trunchiere (a se vedea anexa A). Grație simplicității sale, schema “upwind” a fost aplicată la calcule numerice din dinamica fluidelor. Schema poate fi aplicată ușor pentru probleme 2D și 3D folosind strategia “upwind” pentru toate direcțiile. Un dezavantaj major al schemei “upwind” este că furnizează rezultate cu erori mari când curgerea nu este pe direcția liniilor rețelei de discretizare. Influența erorilor asupra distribuției proprietății de transportativitate este interpretată ca o falsă difuzie (difuzie numerică). Eliminarea difuziei numerice se poate face prin rafinarea rețelei de discretizare dar, în anumite condiții, metoda devine prohibitivă.

5.2.14 Proprietățile schemei hibride

Schema hibridă exploatează proprietățile favorabile ale schemei “upwind” și ale schemei ci diferențe centrale. Schema “upwind” este utilizată acolo unde schema cu diferențe centrale produce erori mari, adică pentru valori mari ale numărului Péclet. Schema este în totalitate conservativă și necondiționat mărginită. Schema satisface exigența proprietății de transportativitate prin utilizarea schemei “upwind” la valori mare ale numărului lui Péclet. Algoritmii numerici bazați pe această schemă sunt foarte stabili. Dezavantajul acestei scheme este că, în funcție de eroarea de trunchiere a dezvoltării în serie Taylor, este numai de ordinul unu.

5.3 Convecție/difuzie 1D nestaționar

5.3.1 Forma generală a ecuației discretizate

Ecuația diferențială generală în cazul 1D nestaționar este următoarea:

, (5.99)

unde este viteza fluidului în direcția .

Dacă se notează:

, (5.100)

unde este fluxul toatal (de convecție și difuzie), se obține ecuația diferențială următoare:

. (5.101)

Se integrează ecuația (5.114) pe volumul de control prezentat în figura 5.20.

(5.102)

Integrând ecuația (5.102) se obține:

, (5.103)

unde

, (5.104)

este fluxul total în direcția prin interfața “” având suprafața (), iar relația

(5.105)

este fluxul total în direcția prin interfața “” având suprafața ().

Se consideră de asemenea ecuația de continuitate care se discretizează folosind aceeași procedură:

. (5.106)

Se integrează ecuația (5.106) pe volumul de control hașurat și prezentat în figura 5.20,

, (5.107)

și se obține după integrarea ecuației (5.107) :

, (5.108)

unde

; (5.109)

. (5.110)

Daca se scade relația (5.108) înmulțită cu din relația (5.103), se obține:

(5.111)

Se evaluează al doilea termen din relația (5.111) în felul următor:

. (5.112)

Ținând cont de relația (5.49) pentru , relația (5.112) devine:

. (5.113)

Ținând cont de relația (5.52), relația (5.113) devine:

. (5.114)

Ținând cont de definiția numărului lui Péclet () relația (5.114) devine:

. (5.115)

Utilizând relația (5.57) pentru se obține succesiv:

; (5.116)

. (5.117)

În mod asemănător, se obține pentru al treilea termen din relația (5.111):

. (5.118)

Ecuația (5.103) se scrie astfel :

.

(5.119)

Regrupând termenii în ecuația (5.119) se obține forma generală a ecuației discretizate:

, (5.120)

cu

;

;

,

unde și sunt valori cunoscute la momentul și

;

.

Funcția , care poate fi aleasă din tabelul 5.2, definește o schemă de discretizare pentru termenul convectiv. Schema “Power Law” este recomandată în acest caz, adică

. (5.121)

5.4 Convecție/difuzie 2D și 3D în regim staționar

Se consideră ecuația diferențială în 3D fără termen sursă:

. (5.122)

Dacă se utilizează, de exemplu, schema hibridă aplicată în același mod ca în cazul 1D și pe celelalte direcții “” și “”, se obține ecuația discretizată care acoperă toate cazurile (1D, 2D și 3D):

, (5.123)

cu coeficientul nodului central :

.

Ceilalți coeficienți sunt prezentați în tabelul 5.6.

Tabelul 5.6 Coeficienții ecuației discretizate pentru schema hibridă

În expresiile de mai sus valorile coeficienților și sunt calculate cu ajutorul formulelor prezentate în tabelul de mai jos:

Condițiile la limită în cazurile 2D și 3D sunt tratate în același mod ca în cazul 1D (a se vedea exemplul 5.3 de la paragraful 5.2.10).

5.5 Convecție/difuzie 2D în regim nestaționar

5.5.1 Forma generală a ecuației discretizate

Ecuația diferențială generală în cazul 2D nestaționar este următoarea:

, (5.124)

unde și sunt componentele vitezei curgerii în direcțiile și .

Dacă se notează:

; (5.125)

, (5.126)

unde și sunt fluxurile totale (de convecție și difuzie), se obține ecuația diferențială sub forma următoare:

. (5.127)

Se integrează ecuația (5.127) pe volumul de control prezentat în figura 5.22

(5.128)

După integrarea ecuației (5.128) se obține:

, (5.129)

unde

;

.

Se consideră, ca în cazul 1D, ecuația de continuitate pentru o problemă 2D:

. (5.130)

Se integrează ecuația (5.130) pe volumul de control:

. (5.131)

După integrarea ecuației (5.131) se obține:

, (5.132)

unde:

;

.

Scăzând relația (5.132) multiplicată cu din relația (5.129), se obține:

. (5.133)

Ca în cazul 1D, se poate arăta că se pot scrie relațiile:

;

;

;

.

Se obține în final forma generală a ecuației discretizate:

, (5.134)

unde:

;

;

;

.

unde:

;

;

;

.

Funcția , care poate fi aleasă din tabelul 5.2, definește o schemă de discretizare pentru termenul convectiv.

=== Cap_6r ===

Capitolul 6

Metoda volumelor finite aplicată la tratarea

problemelor de cuplaj viteză/presiune

6.1 Introducere

Convecția unei variabile scalare depinde atât de valoarea sa cât și de direcția câmpului local al vitezelor. Prezentările din capitolul precedent s-au făcut în ipoteza unui câmp de viteze cunoscut. În general, câmpul de viteze nu este cunoscut și trebuie să fie o parte a soluției globale a variabilelor de curgere.

Se consideră ecuațiile de conservare a cantității de mișcare (a impulsului) în 2D staționar. Câmpul de viteze trebuie să satisfacă de asemenea și ecuația de continuitate.

Ecuația de conservare a cantității de mișcare în direcția :

. (6.1)

Ecuația de conservare a cantității de mișcare în direcția :

. (6.2)

Ecuația de continuitate (de conservare a masei):

. (6.3)

Termenii reprezentând gradienții de presiune în direcțiile și au fost scriși separat pentru a facilita discuțiile care urmează. Rezolvarea sistemului (6.1-6.3) prezintă două noi dificultăți:

prima dificultate – neliniaritatea termenului convectiv din ecuațiile de conservare a cantității de mișcare, de exemplu primul termen din ecuația (6.1), ;

a doua dificultate – în termenul sursă, câmpul de presiune nu este cunoscut și nu există ecuație care să guverneze câmpul de presiune .

Presiunea este specificată indirect prin ecuația de continuitate. Câmpul de presiune corect este câmpul care, atunci când este introdus în ecuația de conservare a cantității de mișcare, pentru un câmp de viteze asociat, satisface ecuația de continuitate.

Cele trei ecuații sunt cuplate pentru ca fiecare componentă a vitezei este prezentă în toate ecuațiile. Dacă gradientul de presiune este cunoscut, obținerea ecuațiilor discretizate este similară cu procedurile prezentate în capitolul 5. Cum câmpul de presiune este o parte a soluției globale a problemei, gradientul de presiune nu este cunoscut. Dacă curgerea este compresibilă, ecuația de continuitate poate fi utilizată ca ecuație de transport pentru densitate și se adaugă la sistemul (6.1-6.3) ecuația de conservare a energiei care este o ecuație de transport pentru temperatură. Presiunea poate fi obținută folosind densitatea și temperatura cu ajutorul ecuației de stare . Totuși, dacă curgerea este incompresibilă atunci densitatea este constantă și nu este legată de presiune. În acest caz cuplajul între presiune și viteză introduce o constrângere asupra soluției câmpului de curgere: dacă soluția corectă a câmpului de presiune este introdusă în ecuațiile de conservare a cantității de mișcare, câmpul de viteze care rezultă verifică ecuația de continuitate.

Toate dificultățile enunțate mai sus vor fi depășite prin utilizarea unei proceduri iterative, cum ar fi algoritmul lui Patankar și Spalding, numit algoritmul SIMPLE (Semi-Implicit Method for Pressure-Linked Equations). În acest algoritm fluxul convectiv pe unitatea de masă , prin interfețele volumului de control, este evaluat folosind un câmp de viteze estimat.

6.2 Reprezentarea gradientului de presiune în ecuațiile cantității de mișcare discretizate

Folosirea metodei volumelor finite începe prin discretizarea domeniului de curgere și a ecuațiilor de transport (6.1–6.3). Mai întâi, trebuie alese punctele în care vor fi stocate vitezele. Ar părea logic să fie stocate în aceleași puncte cu celelalte variabile (presiunea, temperatura, etc.). Totuși, dacă presiunea și viteza sunt definite într-un nod ordinar al volumului de control, un câmp de presiune puternic neuniform poate fi interpretat ca un câmp uniform de către ecuația discretizată a cantității de mișcare. Acest comportament poate fi ilustrat cu ajutorul exemplului prezentat în figura 6.1

Dacă presiunea în punctele ‘’ și ‘’ este obținută prin interpolare liniară, termenul gradientului de presiune din ecuația (6.1) este dat de relația:

. (6.4)

În același mod, termenul gradientului de presiune din ecuația (6.2) este calculat cu relația:

. (6.5)

Se constată că presiunea din nodul central nu apare în expresiile (6.4) și (6.5). Dacă se înlocuiesc valorile presiunii, prezentate în figura 6.1, în relațiile (6.4) și (6.5) se constată că gradienții în punctele rețelei de discretizare sunt nuli chiar dacă câmpul de presiune are oscilații în ambele direcții. În aceste condiții, termenul sursă, în ecuațiile de conservare a cantității de mișcare discretizate, are valoarea zero ca într-un câmp de presiune uniform. Din punct de vedere fizic acest comportament nu este realist. Este clar că dacă viteza este definită în nodurile ordinare ale volumelor de control, influența presiunii nu este în mod corect reprezentată în ecuațiile discretizate ale cantității de mișcare.

6.3 Reprezentarea ecuației de continuitate discretizate

În cazul discretizării ecuației de continuitate apare aceași dificultate. Pentru a ilustra această dificultate se consideră ecuația de continuitate 1D în regim staționar în ipoteza că densitatea este constantă (). Ecuația de continuitate, în aceste condiții, se scrie:

. (6.6)

Dacă se integrează ecuația (6.6) pe volumul de control prezentat în figura 6.2, se obține:

. (6.7)

Considerând interfețele volumului de control situate la jumătatea distanței dintre noduri, ecuația (6.7) devine:

. (6.8)

Consecința ecuației (6.8) este că un câmp de viteze discretizat de tip “dinți de ferăstrău” (nerealist din punct de vedere fizic) poate satisface ecuația de continuitate și poate fi interpretat ca un câmp uniform ().

6.4 Volume de control decalate (deplasate)

Metoda volumelor finite nu impune calcularea diferitelor variabile (, , , , , etc.) în aceleași noduri, adică pe aceeași rețea de volume de control. Ideia este de a calcula variabilele scalare, ca presiunea, densitatea, temperatura, etc. în nodurile ordinare ale volumelor de control, iar componentele vitezei pe volume de control decalate (deplasate), centrate în jurul interfețelor volumelor de control inițiale. În figura 6.3. se prezintă un exemple de astfel de aranjament.

Variabilele scalare, inclusiv presiunea, sunt stocate în nodurile maracate cu (). Componentele vitezei sunt definite ca scalare, între noduri (la interfețe), și sunt indicate prin săgeți. Săgețile orizontale () indică punctele în care este stocată componenta a vitezei, iar săgețile verticale () indică punctele în care este stocată componenta a vitezei. În figura 6.3 s-a introdus un sistem de notație suplimentar bazat pe numerotarea liniilor rețelei și a interfețelor volumelor de control. Acest sistem de notație va fi utilizat mai departe în acest capitol. Pentru moment vom continua să utilizăm notația originală a nodurilor, adică , , , , . Componenta a vitezei este stocată la interfețele ‘’ și ‘’ iar componenta a vitezei este stocată la interfețele ‘’ și ‘’.

În rețeaua decalată, nodurile pentru presiune sunt aceleași cu interfețele volumelor de control pentru componenta a vitezei. Gradientul de presiune este dat de expresia:

, (6.9)

unde este lățimea volumului de control pentru (VC-u). În mod asemănător, gradientul corespunzător volumului de control pentru (VC-v) este dat de relația:

, (6.10)

unde este lățimea lui VC-v. Dacă se calculează gradienții de presiune cu ajutorul relațiilor (6.9) și (6.10) utilizând un câmp de presiune de tip “dinți de ferăstrău” se constată că ei nu mai sunt nuli. O altă consecință importantă a rețelei decalate este că aceasta generează componentele vitezei chiar în punctele în care ele sunt necesare pentru calculul fluxului de convecție ( poate fi calculat fără interpolare la interfețele , , și ).

6.5 Ecuațiile de conservare a cantității de mișcare

La discretizarea ecuațiilor de conservare a cantității de mișcare, se utilizează noul sistem de notare. În figura 6.3 liniile continue sunt numerotate cu litere mari. În direcția numerotarea este …, , , … iar în direcția este …, , , , …etc. Liniile punctate care definesc interfețele sunt numerotate cu litere mici, …, , , , … pentru direcția și …, , , , … pentru direcția . Nodurile scalare (marcate cu ( ) în figura 6.3) sunt localizate la intersecția a două linii ale rețelei de discretizare și sunt identificate prin două litere mari (punctul din figura 6.3 este identificat prin . Componenta a vitezei este stocată în punctele interfețelor ‘’ și ‘’ ale volumelor de control scalare. Aceste puncte sunt localizate la intersecția unei linii ce definește frontiera volumului de control și o linie a rețelei de discretizare. Aceste puncte sunt identificate printr-o combinație de litere mici și litere mari (de exemplu, punctul interfeței ‘’ în jurul punctului este identificat prin . În mod asemănător, componenta a vitezei este identificată printr-o combinație de litere mari – litere mici (de exemplu, punctul interfeței ‘’ din jurul punctului este identificat prin . Se pot utiliza volume de control, pentru calculul vitezei, decalate înainte sau decalate înapoi. Volumele de control pentru componenta (VC-u) și pentru componenta (VC-v) prezentate în figura 6.3 sunt de tipul decalate înapoi.

6.5.1 Integrarea ecuației de conservare a cantității de mișcare

Fie o problemă 2D staționară și volumul de control pentru nodul ‘’, ddecalat înainte, prezentat în figura 6.4. Se integrează ecuația de conservare a cantității de mișcare pentru componenta a vitezei (6.1) pe volumul de control prezentat în figura 6.4.

(6.11)

Integrând (6.11) se obține succesiv:

. (6.12)

(6.13)

Pentru a evalua și se utilizează schemele prezentate în capitolul 5 (upwind, hibridă, QUICK, etc.). Pentru și se alege o interpolare între nodurile vecine, ca de exemplu:

. (6.14)

Relația (6.13), ținând cont de relația (6.14) devine:

(6.15)

În final, după gruparea termenilor, se obține o ecuație discretizată, pentru componenta , sub forma generală următoare:

, (6.16)

sau

, (6.17)

unde

.

Coeficienții și sunt calculați cu una din metodele cu diferențe prezentate la problemele de convecție/difuzie în capitolul 5 (Upwind, Hibridă, Power law, QUICK).

În același mod, se obține ecuația discretiyată pentru componenta , integrând ecuația (6.2), folosind volumul de control din jurul nodului , preyentat în figura 6.5.

Se obține astfel ecuația discretizată următoare:

, (6.18)

unde

.

Dacă se exprimă ecuația (6.17) în funcție de noul sistem de coordonate (din numerotarea nodurilor, fig. 6.3), se obține:

. (6.19)

În noul sistem de numerotare punctele vecine , , și implicate în suma , prezentate în figura 6.6, sunt nodurile , , și . Coeficienții din ecuația (6.19) conțin o combinație a fluxului de convecție și a conductanței de difuzie la interfețele volumelor de control pentru componenta a vitezei. Valorile lui și pentru fiecare interfață , , și a volumului de control sunt date de către relațiile următoare:

(6.20)

(6.21)

(6.22)

; (6.23)

; (6.24)

; (6.25)

; (6.26)

. (6.27)

Formulele (6.20 – 6.27) arată că în punctele în care variabilele scalare sau componentele vitezei nu sunt disponibile, se poate utiliza media a două sau chiar patru puncte vecine pentru care valorile variabilelor sunt disponibile. Pe durata fiecărei iterații valorile componentelor și utilizate pentru calculul coeficienților de mai sus sunt cele de la iterația precedentă. Trebuie remarcat că aceste valori cunoscute, ale lui și , sunt utilizate pentru calculul coeficienților din ecuația (6.19). Aceste valori cunoscute sunt distincte de și care sunt variabilele scalare necunoscute.

În mod asemănător se poate obține ecuația pentru componenta a vitezei:

. (6.28)

Variabilele implicate în suma sunt prezentate în figura 6.7

Coeficienții ecuației (6.28), și sunt constituiți dintr-o combinație a fluxului de convecție și a conductanței de vâscozitate la interfețele volumelor de control pentru componenta a vitezei. Valorile lui și pentru fiecare interfață , , și a volumului de control sunt date de relațiile următoare:

(6.29)

(6.30)

(6.31)

(6.32)

(6.33)

(6.34)

(6.35)

(6.36)

Ecuațiile (6.19) și (6.28) nu pot fi rezolvate decât dacă presiunea este cunoscută sau estimată. Dacă presiunea este cunoscută, câmpul de viteze obținut după rezolvarea sistemului de ecuații algebrice liniare va satisface ecuația de continuitate. Presiunea fiind necunoscută, este necesară o procedură iterativă pentru calculul acesteia.

6.6 Algoritmul SIMPLE

Algoritmul SIMPLE (Semi-Implicit Method for Pressure-Linked Equations) a fost creat de către Patankar și Spalding (1972), [17], și reprezintă o procedură iterativă pentru calculul presiunii folosind o rețea de discretizare deplasată. Procedura iterativă începe cu o presiune estimată. Fie câmpul de presiune estimat. Ecuațiile (6.19) și (6.28) sunt apoi rezolvate pentru a obține câmpul de viteze asociat și :

; (6.37)

. (6.38)

Se definește corecția de presiune ca diferența între presiunea corectă și presiunea estimată :

. (6.39)

În mod asemănător se definesc corecțiile vitezelor și ca diferențe între vitezele corecte și și vitezele estimate și :

; (6.40)

. (6.41)

Substituirea câmpului de presiune corect , în ecuațiile de conservare a cantității de mișcare va furniza câmpul de viteze corect . Ecuațiile discretizate (6.19) și (6.28) leagă câmpul de viteze corect cu câmpul de presiune corect.

Scăzând ecuațiile (6.37) și (6.38) din ecuațiile (6.19) și respectiv (6.28), se obține:

(6.42)

(6.43)

Utilizând formulele de corecție (6.39 – 6.41), ecuațiile (6.42) și (6.43) pot fi rescrise astfel:

; (6.44)

. (6.45)

În acest moment se introduce o aproximare: termenii și sunt neglijați pentru a simplifica ecuațiile (6.44) și (6.45). Neglijarea acestor termeni reprezintă principala aproximare a algoritmului SIMPLE. În aceste condiții se obține:

, (6.46)

, (6.47)

unde

.

Ecuațiile (6.46) și (6.47) descriu corecțiile care trebuie aplicate vitezelor din formulele (6.40) și (6.41). După înlocuiri se obține:

; (6.48)

. (6.49)

Expresii similare se pot scrie și pentru și :

, (6.50)

, (6.51)

unde

.

Până aici s-au considerat numai ecuațiile de conservare a cantității de mișcare, însă câmpul de viteze trebuie să satisfacă în același timp ecuația de continuitate (6.3). Ecuația de continuitate discretizată, obținută prin integrarea ecuației (6.3) pe volumul de control prezentat în figura 6.8, este:

. (6.52)

Substituind ecuațiile corectate (6.48 – 6.51) în ecuația de continuitate discretizată (6.52) se obține:

. (6.53)

Regrupând termenii se obține:

. (6.54)

Identificând coeficienții corecției de presiune ecuația (6.54) poate fi rescrisă sub forma generală următoare:

, (6.55)

unde

.

Ecuația (6.55) reprezintă ecuația de continuitate discretizată ca o ecuație a corecției de presiune . Termenul sursă apare din cauza faptului că se utilizează un câmp de viteze incorect și . Dacă aceasta implică că este nevoie de mai multă corecție de presiune. Prin rezolvarea ecuației (6.55) se obține corecția de presiune pentru toate punctele rețelei de discretizare și atunci presiunea corectă poate fi calculată cu ajutorul formulei (6.39) și componentele vitezei cu ajutorul formulelor de corecție (6.40) și (6.41). Neglijarea termenului nu trebuie să afecteze soluția finală pentru că corecțiile de presiune și de viteză vor fi nule la convergență.

Este posibil ca procesul iterativ să fie divergent. Pentru a remedia acest inconvenient, se poate utiliza subrelaxarea pe durata procesului iterativ:

, (6.56)

unde este un factor de subrelaxare.

Componentele vitezei trebuie de asemenea să fie subrelaxate utilizând relațiile:

; (6.57)

, (6.58)

unde și sunt factorii de subrelaxare pentru componentele vitezei, și fiind componentele corectate fară relaxare iar și reprezintă valorile lor de la iterația precedentă. Se poate demonstra ușor că ecuațiile de conservare a cantității de mișcare, ținând cont de factorul de subrelaxare, pot fi scrise sub forma:

; (6.59)

. (6.60)

Ecuația de corecție a presiunii este de asemenea afectată de către factorul de subrelaxare și se poate arăta că

. (6.61)

Algoritmul SIMPLE este o metodă pentru calculul presiunii și vitezei dar cand alte variabile sunt cuplate la ecuațiile de conservare a cantității de mișcare, de exemplu temperatura, rezolvarea trebuie să fie secvențială. Secvența etapelor algoritmului SIMPLE este prezentată în figura 6.9.

6.7 Algoritmul SIMPLER

Algoritmul SIMPLER (SIMPLE Revised), elaborat de către Patankar (1980), [50], este o versiune ameliorată a algoritmului SIMPLE. Conform acestui algoritm, ecuația de continuitate discretizată (6.52) este utilizată pentru a obține o ecuație discretizată pentru presiune, în locul ecuației de corecție a presiunii din algoritmul SIMPLE. Câmpul de presiune este obținut direct, fără corecție de presiune, însă câmpul de viteze este obținut cu ajutorul corecției folosind ecuațiile (6.48 – 6.51).

Ecuațiile de conservare a cantității de mișcare discretizate (6.19) și (6.28) sunt scrise sub următoarea formă:

; (6.62)

. (6.63)

Conform algoritmului SIMPLER, se definesc pseudo-vitezele și astfel:

; (6.64)

. (6.65)

Ecuațiile (6.62) și (6.63) pot fi scrise astfel:

; (6.66)

. (6.67)

Înlocuind vitezele și date de relațiile (6.66) și (6.67) în ecuația de continuitate discretizată (6.52) se obține:

. (6.68)

Regrupând termenii în ecuația (6.68) se obține ecuația presiunii discretizată:

, (6.69)

unde

.

Se constată că coeficienții din ecuația (6.69) sunt aceeași cu cei din ecuația corecției de presiune (6.55), în afară de termenul sursă care este calculat folosind pseudo-vitezele. Ecuațiile (6.19) și (6.28) sunt rezolvate, utilizând câmpul de presiune obținut cu ajutorul ecuației (6.69), și se obțin componentele vitezelor și . Apoi, ecuațiile de corecție a vitezei (6.48) și (6.49) sunt utilizate pentru a obține vitezele corectate. Totuși, ecuația corecției de presiune (6.55) este folosită pentru corecția vitezei.

Rezumatul algoritmului SIMPLER

Algoritmul constă în rezolvarea ecuației presiunii pentru a obține câmpul de presiune și în rezolvarea ecuației corecției de presiune pentru a corecta vitezele.

Etapele algoritmului sunt:

Se începe cu estimarea (alegerea inițială) a câmpului de viteze, și ,a câmpului de presiune și a variabilei de transport (care în cazul ecuației de conservare a energiei, de exemplu, este temperatura);

Calculul pseudo-vitezelor și , după ce au fost calculați coeficienții ,…,cu ajutorul relațiilor (6.64) și (6.65);

Calculul coeficienților și rezolvarea ecuației presiunii (6.69) pentru a obține câmpul de presiune ;

Inițializarea câmpului de presiune inițial ,cu noul câmp de presiune obținut în etapa 3 () și rezolvarea ecuațiilor de conservare a cantității de mișcare (6.37) și (6.38) pentru a obține și ;

Calculul coeficienților și a termenului sursă și apoi rezolvarea ecuației corecției de presiune (6.55) pentru a obține corecția de presiune ;

Corectarea câmpului de viteze cu ajutorul relațiilor (6.48 – 6.51), fără a corecta presiunea;

Calculul coeficienților și a termenululi sursă și apoi rezolvarea ecuației de transport pentru altă variabilă (temperatura , de exemplu, dacă se utilizează ecuația de conservare a energiei);

Reinițializarea tuturor variabilelor calculate la etapele 3, 6, și 7 (, , , ) și apoi revenire la etapa 2;

Repetarea etapelor 2 la 8 până la obținerea convergenței.

6.8 Algoritmul SIMPLEC

Algoritmul SIMPLEC (SIMPLE Consistent) a fost elaborat de către Van Doormal și Raithby (1984), [35], [81]. Etapele acestui algoritm sunt aproape aceleași ca cele ale algoritmului SIMPLE cu diferența că în ecuațiile de corecție a vitezelor se neglijează termenii nesemnificativi.

Ecuația de corecție pentru componenta este:

, (6.70)

unde

. (6.70a)

În mod asemănător, pentru componenta avem:

, (6.71)

unde

. (6.71a)

Ecuația corecției de presiune este aceeași ca în cazul algoritmului SIMPLE cu excepția termenilor care se calculează cu ajutorul relațiilor (6.70a) și (6.71a). Secvența etapelor este aceeași ca în cazul algoritmului SIMPLE (a se vedea paragraful 6.6).

6.9 Algoritmul PISO

Algoritmul PISO (Pressure Implicit with Splitting of Operators) a fost elaborat de către Issa (1986). Acest algoritm a fost dezvoltat inițial ca o procedură neiterativă pentru calculul curgerilor compresibile nestaționare. Ulterior algoritmul a fost bine adaptat pentru procedura iterativă aplicată la probleme staționare. Algoritmul este o extensie a algoritmului SIMPLE având o etapă suplimentară de corecție.

Ecuațiile de conservare a cantității de mișcare (6.37) și (6.38) sunt rezolvate pornind de la un câmp de presiune pentru a obține componentele și , utilizând aceeași metodă de la algoritmul SIMPLE.

Etapa de corecție 1

Câmpul vitezelor și nu va satisface ecuația de continuitate dacă câmpul de presiune nu este cel corect. Această primă etapă de corecție este introdusă pentru a obține câmpul vitezelor care satisface ecuația de continuitate discretizată. Ecuațiile care rezultă sunt aceleași cu ecuațiile de corecție a vitezelor (6.46) și (6.47) ale algoritmului SIMPLE dar aceasta constitue o primă etapă de corecție în algoritmul PISO, și de aceea se utilizează notații diferite:

;

;

.

Formulele de mai sus sunt folosite pentru a exprima vitezele corectate și :

; (6.72)

. (6.73)

În mod asemănător ca în algoritmul SIMPLE ecuațiile (6.72) și (6.73) sunt substituite în ecuația de continuitate discretizată (6.52) pentru a obține ecuația de corecție a presiunii (6.55). În contextul algoritmului PISO ecuația (6.55) se numește prima ecuație de corecție a presiunii. Rezolvarea acesteia furnizează primul câmp de corecție a presiunii . După ce corecția de presiune este cunoscută, componentele vitezei și pot fi obținute cu ajutorul ecuațiilor (6.72) și (6.73).

Etapa de corecție 2

Ecuațiile de conservare a cantității de mișcare discretizate pentru și sunt:

; (6.74)

. (6.75)

Un al treilea câmp de viteze corectate se obține prin rezolvarea ecuațiilor:

; (6.76)

. (6.77)

Trebuie remarcat că sumele, din ecuațiile de mai sus, sunt evaluate folosind vitezele și calculate în etapa de corecție precedentă.

Scăzând ecuația (6.74) din ecuația (6.76) și ecuația (6.75) din ecuația (6.77) se obține:

; (6.78)

, (6.79)

unde este cea de a doua corecție de presiune, putând fi obținută prin relația următoare:

. (6.80)

Substituind și în ecuația de continuitate discretizată (6.52) se obține a doua ecuație de corecție a presiunii:

, (6.81)

cu și coeficienții vecini dați în tabelul de mai jos:

La obținerea ecuației (6.81) termenul sursă din relația de mai jos

,

este zero dacă componentele și satisfac ecuația de continuitate.

Ecuația (6.81) se rezolvă pentru a obține al doilea câmp de corecție a presiunii și un al treilea câmp de presiune corectat este obținut de către relația:

. (6.82)

În final, se obține un al treilea câmp de viteze corectat cu ajutorul ecuațiilor (6.78) și (6.79).

Rezumatul algoritmului PISO

Etapele algoritmului sunt următoarele:

Estimarea inițială a câmpului de viteze și , a câmpului de presiune și a variabilei de transport ;

Parcurgerea etapelor 1, 2 și 3 ale algoritmului SIMPLE, adică:

Rezolvarea ecuațiilor de conservare a cantității de mișcare;

Rezolvarea ecuației de corecție a presiunii pentru a obține ;

Corectarea presiunii și a vitezelor pentru a actualiza , și .

Calculul coeficienților și a termenului sursă și rezolvarea celei de a doua ecuație de corecție a presiunii (6.81);

Corectarea presiunii și a vitezelor folosind relațiile:

;

;

.

Actualizarea presiunii și a vitezelor, , și ;

Calculul coeficienților și a termenului sursă și apoi rezolvarea ecuației de transport pentru variabila ;

Reinițializarea tuturor variabilelor calculate la etapele 5 și 6 (, , , ) și apoi revenire la etapa 2;

Repetarea etapelor 2 la 7 până la obținerea convergenței.

Acest algoritm implică un considerabil consum de memorie RAM din cauza celei de a doua ecuații de corecție a presiunii. Algoritmul poate fi aplicat cu ușurință la probleme nestaționare.

6.10 Condiții la limită pentru ecuația de corecție a presiunii

Pentru ecuațiile de conservare a cantității de mișcare condițiile la limită sunt tratate ca pentru ecuațiile generale de conservare a variabilei . Pentru ecuația în tratarea condițiilor la limită este specială. Există două tipuri de condiții la limită pentru ecuațiile în :

Fie presiunea este cunoscută pe frontieră: în acest caz, se consideră că pe frontieră și deci pe frontieră;

Fie viteza normală este dată pe frontieră: în acest caz, se alege un volum de control ca cel din figura 6.10, și în ecuația de continuitate nu se înlocuiește prin

,

ci se ia valoarea dată. Nu mai există deci termenul în ecuația de corecție a presiunii, deci nu mai există . Se poate deci considera că în această ecuație (pe frontieră). Concluzie generală: nici o informație despre nnu este necesară pe frontieră.

6.11 Comentarii asupra algoritmilor SIMPLE, SIMPLER, SIMPLEC

și PISO

Neglijarea termenilor de tipul în algoritmul SIMPLE ar fi inacceptabil dacă soluția finală obținută nu ar satisface ecuația de continuitate și ecuația de conservare a cantității de mișcare, dar nu este niciodată cazul.

În cazul general, orice expresie pentru este valabilă atâta timp cât soluția finală (câmpul de viteze) este corectă, în fapt când algoritmul converge și deci .

Rapiditatea convergenței algoritmului SIMPLE depinde de expresia lui . Dacă sunt neglijați prea mulți termeni, algoritmul poate diverge.

“Sursa de masă” în ecuația pentru poate servi ca indicator de convergență a algoritmului SIMPLE: atunci când se tinde către soluție, trebuie să se anuleze pe tot domeniul.

Dacă corecțiile de presiune sunt prea rapide, algoritmul SIMPLE poate diverge și remediul în acest caz este sub-relaxarea.

În algoritmul SIMPLE se pleacă de la un câmp de presiune dat (care este apoi corectat aproximativ) în timp ce în algoritmul SIMPLER, se pornește de la un câmp de viteye dat (verificând ecuația de continuitate) și se deduce câmpul de presiune asociat. Concluzie: avem nevoie de mai puține iterații pentru a ajunge la convergență cu SIMPLER decât cu SIMPLE.

Chiar dacă avem nevoie de mai multe operații elementare în SIMPLER, în final, se poate trage concluzia că SIMPLER converge mai repede decât SIMPLE.

Condițiile la limită pentru ecuația presiunii se tratează în același mod ca la ecuația de corecție a presiunii.

Algoritmii SIMPLEC și PISO sunt la fel de eficienți ca SIMPLER dar nu este clar în ce condiții acești algoritmi sunt mai buni ca SIMPLER. Performanța fiecărui algoritm depinde de condițiile de curgere și de gradul de cuplaj între ecuațiile de conservare a cantității de mișcare și ecuația de transport pentru variabila . În algoritmul PISO se adaugă o corecție suplimentară pentru a ameliora performanța la fiecare iterație. Algoritmul SIMPLER este folosit în diferite software comerciale pentru calculul dinamicii fluidelor (de exemplu FLUENT, http://savas.shef.ac.uk/home/software/fluent )