Resolució numèrica de les equacions de Navier-Stokes

Resolució numèrica de les equacions de Navier-Stokes Treball de Final de Grau: Memòria Grau en Enginyeria en Vehicles Aeroespacials Estudiant : Mario Luna López Director: Assensi Oliva Llena Codirector: Carles-David Pérez Segarra Data: 10 de juny de 2018 Resolució numèrica de les equacions de Navier-Stokes 2 Resolució numèrica de les equacions de Navier-Stokes 3 Resum Aquest document és un estudi sobre la resolució computacional de les equacions governants en la dinàmica de fluids mitjançant mètodes numèrics (CFD). Així, es resolen diferents problemes típics que suposen la introducció de nous conceptes per teixir un coneixement sobre la temàtica des del punt de vista acadèmic. Els problemes tractats bàsicament són de conducció de calor, de l’equació de convecciódifusió, i de les equacions de Navier-Stokes. Aquests es resoldran mitjançant la implementació de codis propis en el llenguatge de programació C++. Els codis es verificaran comparant els resultats de les simulacions amb dades de referència de la literatura científica. Resolució numèrica de les equacions de Navier-Stokes 4 Agraïments En primera instància, m’agradaria agrair el suport rebut pels directors d’aquest treball, els doctors Assensi Oliva Llena i Carles-David Pérez Segarra, per la guia que m’han ofert durant la realització d’aquest projecte, així com al Jordi Chiva per l’ajuda rebuda en resoldre els problemes que apareixien en els meus codis. Al llarg de la carrera, han estat molts professors i companys que m’han transmès la seva passió, i de ben segur la seva influència ha estat notable perquè jo m’interessés per la temàtica del CFD; i en últim estadi, el personal del CTTC en els diversos seminaris en què he pogut aprendre encara més. Tanmateix, vull agrair a la meva família i a les amigues i amics que m’han donat el seu suport emocional i de tot tipus quan ha calgut. Resolució numèrica de les equacions de Navier-Stokes 5 Índex RESUM ......................................................................................................................................... 3 AGRAÏMENTS ............................................................................................................................. 4 1 INTRODUCCIÓ ................................................................................................................. 10 1.1 OBJECTIU........................................................................................................................... 10 1.2 ABAST ................................................................................................................................ 10 1.3 REQUERIMENTS ................................................................................................................ 10 1.4 JUSTIFICACIÓ I ESTAT DE L’ART ...................................................................................... 11 2 INTRODUCCIÓ TEÒRICA .............................................................................................. 13 2.1 EQUACIONS I MALLES ....................................................................................................... 13 2.2 CONDICIONS DE CONTORN ............................................................................................... 15 2.3 SOLVERS ............................................................................................................................. 15 2.3.1 MÈTODE GAUSS-SEIDEL ................................................................................................... 16 2.3.2 MÈTODE LINE-BY-LINE .................................................................................................... 16 2.4 ESTRATÈGIA DE RESOLUCIÓ ............................................................................................ 18 3 CAPÍTOL I: EQUACIÓ DE CONDUCCIÓ DE CALOR .............................................. 19 3.1 INTRODUCCIÓ TEÒRICA.................................................................................................... 19 3.2 INTRODUCCIÓ DEL CAS ..................................................................................................... 20 3.3 DISCRETITZACIÓ DEL DOMINI ......................................................................................... 21 3.4 DISCRETITZACIÓ DE LES EQUACIONS .............................................................................. 22 3.5 ESTUDIS PER VERIFICAR EL CODI .................................................................................... 25 3.6 ESTUDIS NUMÈRICS ........................................................................................................... 26 3.7 ESTUDIS FÍSICS .................................................................................................................. 28 4 CAPÍTOL II: EQUACIÓ DE LA CONVECCIÓ-DIFUSIÓ .......................................... 33 4.1 INTRODUCCIÓ TEÒRICA.................................................................................................... 33 4.2 FORMULACIÓ GENERAL ................................................................................................... 34 4.3 EQUACIÓ DE DISCRETITZACIÓ ......................................................................................... 35 4.4 CAS I: FLUX UNIDIMENSIONAL AMB VARIACIÓ DE LA VARIABLE EN LA MATEIXA DIRECCIÓ DEL FLUX ................................................................................................................... 38 4.4.1 INTRODUCCIÓ .................................................................................................................. 38 4.4.2 RESULTATS ...................................................................................................................... 39 4.5 CAS II: FLUX UNIDIMENSIONAL AMB VARIACIÓ DE LA VARIABLE EN DIRECCIÓ PERPENDICULAR AL FLUX ......................................................................................................... 40 4.5.1 INTRODUCCIÓ .................................................................................................................. 40 4.5.2 RESULTATS ...................................................................................................................... 41 4.6 CAS III: FLUX DIAGONAL ................................................................................................. 41 4.6.1 INTRODUCCIÓ .................................................................................................................. 41 4.6.2 RESULTATS I ESTUDI ....................................................................................................... 42 4.7 CAS IV: FLUX SOLENOÏDAL ............................................................................................. 45 4.7.1 INTRODUCCIÓ .................................................................................................................. 45 4.7.2 RESULTATS I ESTUDI ....................................................................................................... 46 4.8 COMPARATIVA DE SOLVERS ............................................................................................. 50 Resolució numèrica de les equacions de Navier-Stokes 6 5 CAPÍTOL III: EQUACIONS DE NAVIER-STOKES ..................................................... 51 5.1 INTRODUCCIÓ TEÒRICA.................................................................................................... 51 5.2 AVALUACIÓ DELS PASSOS ................................................................................................. 53 5.3 INTRODUCCIÓ CAS: DRIVEN-CAVITY ............................................................................... 56 5.4 ESQUEMES DE CONVECCIÓ-DIFUSIÓ ................................................................................ 57 5.5 ESTUDI DE CONVERGÈNCIA .............................................................................................. 57 5.6 CAS A: REYNOLDS 100 ..................................................................................................... 58 5.7 CAS B: REYNOLDS 1000 .................................................................................................... 61 5.8 CAS C: REYNOLDS 3200 ................................................................................................... 62 5.9 CAS D: REYNOLDS 5000 ................................................................................................... 64 5.10 OPTIMITZACIÓ DEL CODI ............................................................................................... 66 6 CONCLUSIONS I LÍNIES FUTURES ............................................................................. 67 7 BIBLIOGRAFIA ................................................................................................................. 68 Resolució numèrica de les equacions de Navier-Stokes 7 Índex de taules TAULA 3.1: PROPIETATS FÍSIQUES DELS MATERIALS. ............................................................ 20 TAULA 3.3: CONDICIONS DE CONTORN DELS COSTATS. .......................................................... 20 TAULA 4.1: FUNCIÓ A(|P|) PER DIFERENTS ESQUEMES. [1] ..................................................... 38 Resolució numèrica de les equacions de Navier-Stokes 8 Índex de figures FIGURA 2.1: ESQUEMA DE LA DISPOSICIÓ DELS NODES. ........................................................................... 14 FIGURA 2.2: PRÀCTICA DE NODES CENTRATS. [1] .................................................................................... 14 FIGURA 2.3: PRÀCTICA DE CARES CENTRADES. [1] .................................................................................. 15 FIGURA 2.4: ESQUEMA LÍNIA LINE-BY-LINE. [1] ...................................................................................... 17 FIGURA 2.5: ORGANIGRAMA DE L’ALGORITME. ...................................................................................... 18 FIGURA 3.1: ESQUEMA DEL CAS DE 4 MATERIALS. ................................................................................... 20 FIGURA 3.2: ESQUEMA DE LA DISCRETITZACIÓ DEL DOMINI. .................................................................. 21 FIGURA 3.3: MAPA DELS NODES. ............................................................................................................... 22 FIGURA 3.4: MAPA DE TEMPERATURES PER T = 5000 S. ........................................................................... 25 FIGURA 3.5: MAPA DE TEMPERATURES PER T = 5000 S DEL GUIÓ. [2] ..................................................... 26 FIGURA 3.6: NÚMERO D’ITERACIONS EN FUNCIÓ DEL FACTOR DE RELAXACIÓ. ..................................... 27 FIGURA 3.7: MAPA DE TEMPERATURES DONATS 24 ºC INCIALS. .............................................................. 28 FIGURA 3.8: MAPA DE TEMPERATURES PER Q = 240.00 W/M.................................................................. 29 FIGURA 3.9: MAPA DE TEMPERATURES PER Α = 27 W/M2K. .................................................................... 29 FIGURA 3.10: MAPA DE TEMPERATURES PER ISOTERMA INFERIOR T = 30 ºC. ........................................ 30 FIGURA 3.11: MAPA DE TEMPERATURES PER ISOTERMA INFERIOR T = 16 ºC. ....................................... 30 FIGURA 3.12: MAPA TEMPERATURES AMB DOBLE PENDENT TEMPS........................................................ 31 FIGURA 3.13: MAPA TEMPERATURES AMB MEITAT PENDENT TEMPS. ..................................................... 31 FIGURA 3.14: MAPES DE TEMPERATURA PER CONDUCTIVITATS DIFERENTS. ......................................... 32 FIGURA 4.1: : ESQUEMA DEL FLUX TOTAL QUE TRAVESSA LA CARA D’UN VOLUM DE CONTROL ENTRE DOS NODES. [1] .................................................................................................................................. 34 FIGURA 4.2: VOLUM DE CONTROL PER LA SITUACIÓ 2-D. [1] .................................................................. 35 FIGURA 4.3: FLUX UNIDIMENSIONAL AMB VARIACIÓ DE LA VARIABLE EN LA DIRECCIÓ DEL FLUX. [4] 38 FIGURA 4.4: RESULTATS PER PECLET 1.................................................................................................... 39 FIGURA 4.5: RESULTATS PER PECLET -10. ............................................................................................... 39 FIGURA 4.6: FLUX UNIDIMENSIONAL AMB VARIACIÓ DE LA VARIABLE EN DIRECCIÓ PERPENDICULAR AL FLUX. [4] ...................................................................................................................................... 40 FIGURA 4.7: RESULTATS DEL 2N CAS. ....................................................................................................... 41 FIGURA 4.8: CAS DE FLUX DIAGONAL. [4] ................................................................................................. 41 FIGURA 4.9: SOLUCIÓ DEL FLUX DIAGONAL PER MALLA DE 100X100 I PECLET D’1. ............................. 42 FIGURA 4.10: SOLUCIÓ DEL FLUX DIAGONAL PER MALLA DE 100X100 I PECLET DE 100. ...................... 43 FIGURA 4.11: SOLUCIÓ DEL FLUX DIAGONAL PER MALLA DE 100X100 I PECLET D’1 M. ....................... 43 FIGURA 4.12: SOLUCIÓ DEL FLUX DIAGONAL PER MALLA DE 500X500 I PECLET D’1M. ........................ 44 FIGURA 4.13: SOLUCIÓ DEL FLUX DIAGONAL PER MALLA DE 1000X1000 I PECLET D’1 M. ................... 44 FIGURA 4.14: ESQUEMA DEL PROBLEMA DE FLUX SOLENOÏDAL. [4] ....................................................... 45 FIGURA 4.15: PERFIL DE LA VARIABLE A L’ENTRADA. ............................................................................. 46 FIGURA 4.16: RESULTATS A LA SORTIDA DEL FLUX SOLENOÏDAL PER 𝝆𝚪=𝟏𝟎. .................................... 46 FIGURA 4.17: MAPA DE LA VARIABLE GENERAL PER 𝝆𝚪=𝟏𝟎. ............................................................... 47 FIGURA 4.18: RESULTATS A LA SORTIDA DEL FLUX SOLENOÏDAL PER 𝝆𝚪=𝟏 𝟎𝟎𝟎. .............................. 47 FIGURA 4.19: MAPA DE LA VARIABLE GENERAL PER 𝝆𝚪=𝟏 𝟎𝟎𝟎. ......................................................... 48 FIGURA 4.20: RESULTATS A LA SORTIDA DEL FLUX SOLENOÏDAL PER 𝝆𝚪=𝟏 𝟎𝟎𝟎 𝟎𝟎𝟎. ...................... 48 FIGURA 4.21: MAPA DE LA VARIABLE GENERAL PER 𝝆𝚪=𝟏 𝟎𝟎𝟎 𝟎𝟎𝟎. ................................................. 49 FIGURA 4.22: LÍNIES DE CORRENT DEL FLUX SOLENOÏDAL. .................................................................... 49 FIGURA 4.23: COMPARATIVA ENTRE ELS SOLVERS DE GAUSS-SEIDEL I LINE-BY-LINE. .......................... 50 FIGURA 5.1: MALLA DESPLAÇADA. [5] ..................................................................................................... 53 FIGURA 5.2: ESQUEMA DEL CAS TRACTAT. [5] ......................................................................................... 56 FIGURA 5.3: ERROR COMPUTAT EN L’ESTUDI DE CONVERGÈNCIA. ......................................................... 57 FIGURA 5.4: TEMPS DE SIMULACIÓ EN L’ESTUDI DE CONVERGÈNCIA. .................................................... 58 FIGURA 5.5: MAPA DE PRESSIONS. ............................................................................................................ 58 FIGURA 5.6: PERFIL DE VELOCITATS U EN LA LÍNIA CENTRAL PER REYNOLDS 100. .............................. 59 FIGURA 5.7: PERFIL DE VELOCITATS V EN LA LÍNIA CENTRAL PER REYNOLDS 100. .............................. 59 FIGURA 5.8: MAPA DEL MÒDUL DE LA VELOCITAT PER REYNOLDS 100. ................................................ 60 FIGURA 5.9: LÍNIES DE CORRENT PER REYNOLDS 100. ............................................................................ 60 FIGURA 5.10: PERFIL DE VELOCITATS U EN LA LÍNIA CENTRAL PER REYNOLDS 1000. .......................... 61 FIGURA 5.11: PERFIL DE VELOCITATS V EN LA LÍNIA CENTRAL PER REYNOLDS 1000. .......................... 61 FIGURA 5.12: MAPA DEL MÒDUL DE LA VELOCITAT PER REYNOLDS 1000. ............................................ 62 Resolució numèrica de les equacions de Navier-Stokes 9 FIGURA 5.13: LÍNIES DE CORRENT PER REYNOLDS 1000. ........................................................................ 62 FIGURA 5.14: PERFIL DE VELOCITATS U EN LA LÍNIA CENTRAL PER REYNOLDS 3200. .......................... 63 FIGURA 5.15: PERFIL DE VELOCITATS V EN LA LÍNIA CENTRAL PER REYNOLDS 3200. .......................... 63 FIGURA 5.16: MAPA DEL MÒDUL DE LA VELOCITAT PER REYNOLDS 3200. ............................................ 64 FIGURA 5.17: LÍNIES DE CORRENT PER REYNOLDS 3200. ........................................................................ 64 FIGURA 5.18: DIFERÈNCIA PER REYNOLDS 5000 I MALLA 50X50. ........................................................... 65 FIGURA 5.19: DIFERÈNCIA PER REYNOLDS 5000 I MALLA 80X80. ........................................................... 65 Resolució numèrica de les equacions de Navier-Stokes 16 2.3.1 Mètode Gauss-Seidel Aquest mètode és un mètode de punt per punt (point-by-point), en què es resol el problema per aquesta via. Cada equació de la forma (2.1) es resol aïllant la variable en el node tractat P: 𝜙𝑃=(𝑎𝐸𝜙𝐸+𝑎𝑊𝜙𝑊+𝑎𝑁𝜙𝑁+𝑎𝑆𝜙𝑆+𝑏𝑃)/𝑎𝑃 (2.6) Així, el valor de la variable en el node P s’obté com a com a combinació lineal dels nodes del voltant, que tenen un valor suposat o l’últim calculat. D’aquesta manera, un cop que s’han visitat tots els nodes, s’ha finalitzat la iteració, i es comprova si la màxima diferència entre el valor suposat i calculat és major a la precisió escollida, és a dir, l’error admissible. Si la solució encara no és suficientment precisa, es torna a iterar amb els valors calculats en l’última iteració, que passen a ser els que se suposen en la iteració actual. 2.3.2 Mètode Line-by-Line TDMA Per explicar aquest mètode, cal explicar prèviament el mètode del TDMA, de l’anglès, Tri-Diagonal Matrix Algorithm, un mètode que requereix d’una matriu de coeficients a tri-diagonal: (𝑎𝑃[1] −𝑎𝐸[1] 0 0 ⋯ −𝑎𝑊[2] 𝑎𝑃[2] −𝑎𝐸[2] 0 … 0⋮−𝑎𝑊[3] ⋮ 𝑎𝑃[3] ⋮ −𝑎𝐸[3] ⋮ ⋱) (2.7) Com es pot comprovar, pel requeriment de la tri-diagonalitat de la matriu, els coeficients que hi apareixen dels nodes dels voltants tan sols són east o west, fet que implica que el mètode del TDMA sigui aplicable només a casos unidimensionals. Aquest mètode no és iteratiu sinó de resolució directa. Es basa en 2 passos: 1. Evaluació de i=1 a N dels coeficients P i R: 𝑃[𝑖]=𝑎𝐸[𝑖] 𝑎𝑃[𝑖]−𝑎𝑊[𝑖]·𝑃[𝑖−1] (2.8) 𝑅[𝑖]=𝑏𝑃[𝑖]+𝑎𝑊[𝑖]·𝑅[𝑖−1] 𝑎𝑃[𝑖]−𝑎𝑊[𝑖]·𝑃[𝑖−1] (2.9) Resolució numèrica de les equacions de Navier-Stokes 17 2. Obtenció dels valors de la variable de i=N a 1: 𝜙[𝑖]=𝑃[𝑖]·𝑇[𝑖+1]+𝑅[𝑖] (2.10) Line-by-Line El mètode de línia per línia és una combinació dels mètodes de Gauss-Seidel i TDMA. La idea és que, per exemple, per un problema bidimensional com els que es tractaran en aquest estudi, es pugui resoldre el problema línia per línia, on els nodes del voltant que no formen part de la línia queden congelats. A partir de l’equació general (2.1), aquesta es pot reescriure com: 𝑎𝑃𝜙𝑃=𝑎𝐸𝜙𝐸+𝑎𝑊𝜙𝑊+𝑏𝑃∗ (2.11) On: 𝑏𝑃∗=𝑏𝑃+𝑎𝑁𝜙𝑁+𝑎𝑆𝜙𝑆 (2.12) És a dir, que els nodes que no formen part de la línia es congelen i passen a formar part del terme independent. En aquest exemple, es resoldria la línia en direcció x ja que hi formen part els nodes east i west, mentre que els south i north estarien congelats. Un cop resolta aquesta línia, es passaria a la següent. Els valors emprats en la resolució, tal i com fa el Gauss-Seidel, serien els suposats o últims calculats. Figura 2.4: Esquema línia Line-by-Line. [1] En la Figura 2.4 anterior, es troba l’esquema d’aquesta metodologia, on els nodes principals de la línia estan marcats per punts, i els nodes dels voltants marcats per creus estan congelats i formen part del nou terme independent. Amb aquesta metodologia, és convenient escombrar en un sentit i després en l’altre; és a dir, per exemple, primer d’esquerra a dreta i després de dalt a baix. Cal comentar que, a fi que la informació es transmeti amb més rapidesa, cal començar escombrant en contorns on la condició de contorn sigui de Dirichlet. Això és perquè el valor en aquells nodes ja és conegut perquè és dada, i així es resol més ràpidament el problema. Resolució numèrica de les equacions de Navier-Stokes 18 2.4 Estratègia de resolució Per resoldre els problemes dels casos que es plantejaran, caldrà implementar un algoritme bàsic, que depenent del cas tindrà més o menys passos però que bàsicament seguiran el mateix organigrama (veure Figura 2.5). Aquest algoritme s’explica a continuació: 1) El primer pas és l’entrada de dades. Les dades físiques que s’introdueixen són la geometria, les propietats dels materials o els fluids, les condicions de contorn i inicials, entre d’altres. Les dades numèriques introduïdes són, per exemple, el número de nodes interns, el paràmetre de l’esquema numèric, la precisió dels càlculs, entre d’altres. 2) Els càlculs previs calculen, per exemple, la posició dels nodes, la posició de les cares dels volums de control, la distància entre nodes, les superfícies, els volums, etc. 3) El mapa inicial assigna el valor de les variables inicials o suposades a les matrius de les variables. 4) El càlcul dels coeficients de discretització es realitza segons el valor dels coeficients que resulten de discretitzar les diverses equacions. 5) El solver emprat és el mètode Gauss-Seidel o Lineby-Line, que calcula les variables amb valors suposats. Un cop calculades, es comprova que les diferències entre les variables calculades i suposades siguin prou petites; és a dir, que es compleixi el criteri de convergència. En cas contrari, en Gauss-Seidel, s’aplica un factor de relaxació per minimitzar les iteracions necessàries per assolir la convergència. 6) Es comprova si cal un nou time-step (no s’hagi arribat al temps de simulació desitjat). En cas de necessitar-ho, es torna a repetir el procés actualitzant el mapa de la variable anterior amb el recentment calculat (salt temporal). Cal notar que en cas de ser transitori, hi ha termes com el bp que es tornaran a calcular després de cada pas de temps ja que depenen de la variables en l’instant anterior. És per això que, al respondre positivament que cal un nou time-step, l’algoritme marca que cal anar entre el pas 4) i 5), ja que només cal recalcular una part del pas 4). Figura 2.5: Organigrama de l’algoritme. Resolució numèrica de les equacions de Navier-Stokes 19 3 Capítol I: Equació de conducció de calor 3.1 Introducció teòrica L’equació de la conducció de calor és l’equació governant en molts processos de transferència de calor. El fenomen es produeix en materials sòlids o alguns fluids sense moviment, en què es transfereix calor a causa d’una diferència de temperatures. La simplicitat d’aquest cas fa que sigui un punt de partida òptim per començar a formular un coneixement sobre les equacions que es tractaran en els següents capítols, ja que contemplar les similituds entre la transferència de calor i la transferència de momentum i concebre la velocitat, en alguns aspectes, de manera anàloga a la temperatura, ajuda a l’enteniment dels fenòmens que es desenvoluparan posteriorment. La calor es transfereix de zones amb major a menor temperatura, és per això que el gradient de temperatures és negatiu. La conducció de calor està determinada per la Llei de Fourier: 𝑞𝑛=−𝜆𝑑𝑇 𝑑𝑛 (3.1) on qn és el flux de calor en direcció n per unitat de superfície (W/m2) i λ és la conductivitat tèrmica. A fi de provar el coneixement adquirit sobre aquest fenomen, s’analitzarà un cas de conducció de calor proposat pel CTTC que es pot trobar en [2]. Resolució numèrica de les equacions de Navier-Stokes 20 3.2 Introducció del cas El cas tracta bàsicament d’una secció rectangular composta per quatre materials diferents (vegeu Figura 3.1). Es pren com a direcció x la direcció horitzontal i com a direcció y la direcció vertical. Les diferents distàncies L1, L2, L3, L4 i L5 conformen les distàncies particulars del problema. Les diferents propietats dels materials es mostren a la Taula 3.1. Cada costat té una condició de contorn diferent, descrites a la Taula 3.2. Figura 3.1: Esquema del cas de 4 materials. Taula 3.1: Propietats físiques dels materials. [2] Taula 3.2: Condicions de contorn dels costats. Costat Condició de contorn Inferior Isoterma a T = 23.00 ºC Superior Flux de calor uniforme Q = 60 .00 W/m Esquerre En contacte amb un fluid a Tg=33 ºC i coef. transf. calor 9.00 W/m2K Dret Temperatura uniforme T=8.00+0.005t ºC (on t és el temps en segons) Resolució numèrica de les equacions de Navier-Stokes 21 3.3 Discretització del domini Per tal d’obtenir el mapa de temperatures del domini per un instant donat de temps, cal discretitzar el domini en volums de control suficientment petits talment les propietats (en aquest cas, la temperatura) romanguin constants en aquell instant de temps i no variïn dins d’aquest. La variable calculada, la temperatura, es captarà col·locant-hi un node en el centre del volum de control, tal i com es comentava anteriorment. Figura 3.2: Esquema de la discretització del domini. A la Figura 3.2 es mostra la pràctica emprada per a construir la malla del problema. Els elements N1, N2 són el número de nodes interns o volums de control en què es divideix el problema en la direcció x, on N1 correspon a 0<x< L1 i N2 a L1<x< L1+L2. A més, a l’extrem esquerre i dret es col·loquen 2 nodes més per captar les condicions de contorn. El mateix passa en la direcció y, on els elements M3, M4 són el número de nodes interns o volums de control en què es divideix el problema en la direcció y, on M3 correspon a 0<y< L3 i M4 a L3<y< L3+L4 +L5. A més, a dalt i avall es col·loquen 2 nodes més per captar les condicions de contorn. D’aquesta manera, es poden diferenciar els nodes interns, els col·locats a les condicions de contorn, i els dels vèrtexs del rectangle. El fet de triar diversos números d’elements en direcció x i y fa possible densificar la malla allà on es creu que caldrà una malla més densa. Si no es fes així i s’emprés una malla uniforme, caldria augmentar la densitat de malla en tot el domini a causa de la necessitat de només una zona, densificant malla on no cal i malbaratant potència computacional i augmentant el temps de càlcul. Aquesta divisió s’ha triat de manera que es pugui densificar la malla en la part inferior del material 2, ja que per les condicions de contorn donades, sembla ser la zona on es forcin uns gradients de temperatura més grans. Aquesta tria també ve donada per la pista que s’ofereix en el guió del cas ([2]), on es mostra el mapa de temperatures a un instant donat, i es pot verificar que es compleix aquesta deducció per la tria de la malla. En relació a la densificació de malla, no obstant, es troba un inconvenient. Per tal de durla a terme, es tria el número de volums de control o nodes en trams de l’eix horitzontal i Resolució numèrica de les equacions de Navier-Stokes 22 vertical. Aquesta pràctica fa que, tot i densificar la malla en la zona desitjada, també fa que es densifiqui més la malla en altres zones. Per exemple, els elements M3 es trien en l’escala vertical, fent que la zona de 0<x< L1 també es vegi afectada. Aquest fenomen, no obstant, no s’optimitzarà doncs fer que només en una zona s’hi densifiqués la malla faria que als voltants hi corresponguessin més d’un volum de control intern per cada volum de control extern a l’hora de veure la calor transferida entre aquests. En altres paraules, la malla serà estructurada i en tot aquest estudi no s’empraran malles no estructurades. En la Figura 3.3 es mostra el mapa dels nodes per N1=10, N2=20, M3=30 i M4=10. Es pot comprovar com l’esquema de la Figura 3.2 intentava representar aquesta malla desigual. Figura 3.3: Mapa dels nodes. 3.4 Discretització de les equacions L’equació governant del problema és l’equació de l’energia per sòlids, transitòria, bidimensional, amb propietats físiques constants i sense fonts internes: 𝜌𝑐𝑃𝜕𝑇 𝜕𝑡=𝜆(𝜕2𝑇 𝜕𝑥2+𝜕2𝑇 𝜕𝑦2) (3.2) El procés de discretització es considera conegut, només caldria integrar-la. Resolució numèrica de les equacions de Navier-Stokes 23 1) Nodes interns L’equació de discretització dels nodes interns és la següent: 𝑎𝑃𝑇𝑃=𝑎𝐸𝑇𝐸+𝑎𝑊𝑇𝑊+𝑎𝑁𝑇𝑁+𝑎𝑆𝑇𝑆+𝑏𝑃 (3.3) On: 𝑎𝐸=𝛽𝜆𝑒𝑆𝑒 𝑑𝑃𝐸 (3.4) 𝑎𝑊=𝛽𝜆𝑤𝑆𝑤 𝑑𝑃𝑊 (3.5) 𝑎𝑁=𝛽𝜆𝑛𝑆𝑛 𝑑𝑃𝑁 (3.6) 𝑎𝑆=𝛽𝜆𝑠𝑆𝑠 𝑑𝑃𝑆 (3.7) 𝑎𝑃=𝑎𝐸+𝑎𝑊+𝑎𝑁+𝑎𝑆+𝜌𝑉𝑐 Δ𝑡 (3.8) 𝑏𝑃=𝑇𝑃0(𝜌𝑉𝑐 Δ𝑡−(1−𝛽)(𝜆𝑤𝑆𝑤 𝑑𝑃𝑊+𝜆𝑒𝑆𝑒 𝑑𝑃𝐸+𝜆𝑠𝑆𝑠 𝑑𝑃𝑆+𝜆𝑛𝑆𝑛 𝑑𝑃𝑁)) +𝑇𝑊 0(1−𝛽)𝜆𝑤𝑆𝑤 𝑑𝑃𝑊+𝑇𝐸0(1−𝛽)𝜆𝑒𝑆𝑒 𝑑𝑃𝐸 +𝑇𝑆0(1−𝛽)𝜆𝑠𝑆𝑠 𝑑𝑃𝑆+𝑇𝑁0(1−𝛽)𝜆𝑛𝑆𝑛 𝑑𝑃𝑁 (3.9) On 𝛽 imposa l’esquema numèric d’integració temporal (endavant s’estudiarà), 𝑆 són les superfícies de les cares dels volums de control, 𝑑𝑃𝐸 és la distància entre els nodes, 𝑉 és el volum del volum de control, Δ𝑡 és el time-step d’integració. Cal notar que les temperatures en l’instant actual no porten cap superíndex, mentre que les de l’instant anterior porten el superíndex 0. Aquesta equació (3.3) relaciona linealment la temperatura del node tractat P amb la temperatura dels nodes que l’envolten mitjançant uns coeficients de discretització i un terme independent, tal i com feia l’equació general (2.1). Quant a les conductivitats tèrmiques, s’emprarà la mitjana harmònica en tots els càlculs ja que es tracta de nodes centrats amb canvis de material: 𝜆𝑒=𝑑𝑃𝐸 𝑑𝑃𝑒 𝜆𝑃+𝑑𝐸𝑒 𝜆𝐸 (3.10) Resolució numèrica de les equacions de Navier-Stokes 24 2) Nodes contorn inferior Els nodes del contorn inferior no tenen una equació de discretització ja que la seva temperatura ja és coneguda: es tracta d’una isoterma de 23 ºC. Altrament, es podria imposar que els coeficients dels nodes dels voltants són nuls, el propi és 1 i el terme independent seria igual a la temperatura donada. 3) Nodes contorn superior L’equació de discretització dels nodes del contorn superior és la següent: 𝑎𝑃𝑇𝑃=𝑎𝑆𝑇𝑆+𝑏𝑃 (3.11) 𝑎𝑆=𝜆𝑠 𝑑𝑃𝑆 (3.12) 𝑎𝑃=𝑎𝑆 (3.13) 𝑏𝑃=𝑞󰇗𝑓𝑙𝑢𝑥 (3.14) Aquesta equació (3.11) relaciona la temperatura dels nodes del contorn superior amb els inferiors interns, i no hi relaciona nodes del contorn entre si ja que entre ells no hi ha transferència de calor. Sorgeix, bàsicament, de fer un balanç de calor entre el flux de calor conegut del contorn i la conducció de calor interior. 4) Nodes contorn esquerre L’equació de discretització dels nodes del contorn esquerre és la següent: 𝑎𝑃𝑇𝑃=𝑎𝐸𝑇𝐸+𝑏𝑃 (3.15) 𝑎𝐸=𝜆𝑒 𝑑𝑃𝐸 (3.16) 𝑎𝑃=𝑎𝐸+𝛼𝑔 (3.17) 𝑏𝑃=𝛼𝑔𝑇𝑔 (3.18) Aquesta equació (3.15) relaciona la temperatura dels nodes del contorn esquerre amb els drets interns, i no hi relaciona nodes del contorn entre si ja que entre ells no hi ha transferència de calor. Resolució numèrica de les equacions de Navier-Stokes 25 Sorgeix, bàsicament, de fer un balanç de calor entre la convecció natural del contorn i la conducció de calor interior. 4) Nodes contorn dret Els nodes del contorn dret no tenen una equació de discretització ja que la seva temperatura ja és coneguda per a un instant donat de temps, segons l’equació donada. Altrament, es podria imposar que els coeficients dels nodes dels voltants són nuls, el propi és 1 i el terme independent seria igual a la temperatura donada en cada instant de temps. 5) Nodes als vèrtexs Els nodes als vèrtexs són irrellevants, ja que no hi ha superfície de transferència de calor, i la temperatura que se’ls hi assignarà serà la mitjana entre els dos nodes més propers a aquests. Aquesta pràctica només es realitzarà de manera que a l’imprimir el mapa de temperatures doni un valor raonable, però no afectarà en la resolució del problema. 3.5 Estudis per verificar el codi A continuació es detallen diversos estudis per verificar que el codi està lliure d’errors. Primer, cal contemplar que es tracta d’un cas de conducció de calor bidimensional en règim transitori. Per tant, es desconeix cap solució analítica per comparar-lo. D’aquesta manera, el guió de la pràctica proporcionava una pista, que és el mapa de temperatures per a l’instant de temps de 5 000 s. Així, es va demanar al programa implementat que imprimís per pantalla els resultats de les temperatures calculades per a aquest mateix instant de temps. El resultat es pot veure a la Figura 3.4, mentre que a la Figura 3.5 es mostra la pista proporcionada al guió per comparar els resultats. Figura 3.4: Mapa de temperatures per t = 5000 s. Resolució numèrica de les equacions de Navier-Stokes 32 Un cop estudiades les variacions de les 4 condicions de contorn possibles, s’arriba a la conclusió que els grans canvis sobre els mapes de temperatures resultants es donen quan es varia el valor d’una isoterma dada. Pel cas de la convecció natural i el flux de calor donat, per variacions no exagerades, el resultat acaba sent similar. Un altre estudi físic que es pot realitzar es variant algunes propietats físiques dels materials. En aquest cas, s’ha decidit variar les conductivitats tèrmiques per veure’n l’efecte i si hi havia més calor transferida. Només s’han variat les conductivitats tèrmiques dels materials de la dreta, que són els que presentaven més canvis de temperatura. El resultat és mostra a la Figura 3.14. Figura 3.14: Mapes de temperatura per conductivitats diferents. Tal i com es pot apreciar, les isotermes s’han eixamplat, ergo no hi ha gradients de temperatura tan grans. Aquest és un resultat correcte ja que dóna a entendre que com la conductivitat ha pujat, es transfereix més calor que redistribueix més les temperatures. Resolució numèrica de les equacions de Navier-Stokes 33 4 Capítol II: Equació de la convecció-difusió 4.1 Introducció teòrica La formulació del tema anterior pot ser emprada per aquest capítol de manera preliminar, només cal afegir el terme de la convecció, fenomen creat pel moviment del fluid. La formulació de l’equació de discretització és anàloga, doncs, a la del cas de conducció de calor, atès que la variable de temperatura 𝑇 i la constant de conductivitat tèrmica 𝜆 són homòlogues a la variable general 𝜙 i el seu coeficient de difusió Γ. D’aquesta manera, l’objectiu de l’estudi de l’equació de convecció-difusió és obtenir la distribució de la variable general 𝜙 donat el camp de velocitats i les propietats del fluid. Aquesta variable general podria representar, per exemple, la temperatura en cada punt. Tot i que l’únic terme nou és el de convecció, la seva formulació no és directa ni banal. Per obtenir resultats amb sentit físic, cal estudiar amb cura com formular la discretització d’aquest terme, el qual està lligat al terme de difusió. És per això que anomenem l’equació següent com equació de la convecció-difusió: 𝜕 𝜕𝑡(𝜌𝜙)+∇·(𝜌𝑣𝜙)=∇·(Γ∇𝜙)+𝑆 (4.1) on el 1r terme és el transitori, el 2n el convectiu, el 3r el difusiu i el 4t el terme font. Per altra banda, cal tenir en compte que el camp fluid donat ha de satisfer l’equació de la continuïtat: 𝜕𝜌 𝜕𝑡+∇·(𝜌𝑣)=0 (4.2) la qual no és més que l’equació de la convecció-difusió amb 𝜙=1 i 𝑆=0. Resolució numèrica de les equacions de Navier-Stokes 34 4.2 Formulació general El complet desenvolupament del terme convectiu i els diferents esquemes de resolució no es durà a terme en aquest document ja que es troben desenvolupats per Patankar [1] i el manuscrit del CTTC [3]. A continuació es desenvoluparan breument conceptes necessaris per entendre bàsicament els paràmetres que intervenen en l’equació de discretització. Considerant flux unidimensional, estacionari i sense fonts, és útil definir el flux 𝐽 total que travessa la cara d’un volum de control, i que contempla el flux convectiu i el flux difusiu, respectivament: 𝐽=𝜌𝑢𝜙−Γ𝑑𝜙 𝑑𝑥 (4.3) Figura 4.1: : Esquema del flux total que travessa la cara d’un volum de control entre dos nodes. [1] Abans de seguir amb aquest desenvolupament cal definir dues variables, la força convectiva 𝐹 i la conductància de la difusivitat 𝐷: 𝐹≡𝜌𝑢 (4.4) 𝐷≡Γ𝛿 (4.5) D’aquestes dues últimes equacions, es pot extreure el número adimensional de Peclet: 𝑃≡F𝐷=𝜌𝑢𝛿 Γ (4.6) El número de Peclet és la relació entre la força convectiva i difusiva. Quan aquest número s’aproxima a 0, tenim un problema amb una difusió quasi bé pura, mentre que quan el Peclet és molt gran, el valor de les variables en un node ve determinat majoritàriament Resolució numèrica de les equacions de Navier-Stokes 35 pels valors que presenten el node aigües amunt, mentre que la influència del node aigües avall és quasi bé nul·la. Tornant al desenvolupament del flux total que travessa la cara d’un volum de control entre dos nodes, es pot redefinir: 𝐽∗=𝐽𝛿 Γ=𝑃𝜙− 𝑑𝜙 𝑑(𝑥/𝛿) (4.7) on el valor de 𝜙 a la cara serà algun tipus de mitjana ponderada entre 𝜙𝑖 i 𝜙𝑖+1, mentre que el gradient serà múltiple de 𝜙𝑖+1-𝜙𝑖. Tenint en compte aquests aspectes i agrupant termes: 𝐽∗=𝐵𝜙𝑖−𝐴𝜙𝑖+1 (4.8) on 𝐴 i 𝐵 són coeficients adimensionals que són funció del número de Peclet P, cadascun associat al node corresponent. És possible escriure tot en termes del coeficient 𝐴 mitjançant relacions amb el coeficient 𝐵, i així es farà d’ara en endavant, ja que ambdós depenen únicament del número de Peclet. 4.3 Equació de discretització Tenint en compte la definició del flux 𝐽 que travessa les cares d’un volum de control donat, és possible discretitzar l’equació de convecció-difusió per dues dimensions i règim transitori amb termes font. Figura 4.2: Volum de control per la situació 2-D. [1] L’equació de discretització obtinguda i els coeficients es presenten a continuació: 𝑎𝑃𝜙𝑃=𝑎𝐸𝜙𝐸+𝑎𝑊𝜙𝑊+𝑎𝑁𝜙𝑁+𝑎𝑆𝜙𝑆+𝑏 (4.9) Resolució numèrica de les equacions de Navier-Stokes 36 on: 𝑎𝐸=𝐷𝑒𝐴(|𝑃𝑒|)+ 𝑚𝑎𝑥 (−𝐹𝑒,0) (4.10) 𝑎𝑊=𝐷𝑤𝐴(|𝑃𝑤|)+ 𝑚𝑎𝑥 (𝐹𝑤,0) (4.11) 𝑎𝑁=𝐷𝑛𝐴(|𝑃𝑛|)+ 𝑚𝑎𝑥 (−𝐹𝑛,0) (4.12) 𝑎𝑆=𝐷𝑠𝐴(|𝑃𝑠|)+ 𝑚𝑎𝑥 (𝐹𝑠,0) (4.13) 𝑎𝑃=𝑎𝐸+𝑎𝑊+𝑎𝑁+𝑎𝑆+𝜌𝑃0∆x∆y ∆𝑡 −𝑆𝑃∆x∆y (4.14) 𝑏=𝑆𝐶∆x∆y+𝑎𝑃0𝜙𝑃0 (4.15) i, on: 𝐷𝑒=Γ𝑒∆y 𝑑𝑃𝐸 (4.16) 𝐷𝑤=Γ𝑤∆y 𝑑𝑃𝑊 (4.17) 𝐷𝑛=Γ𝑛∆x 𝑑𝑃𝑁 (4.18) 𝐷𝑠=Γ𝑠∆x 𝑑𝑃𝑆 (4.19) 𝐹𝑒=(𝜌𝑢)𝑒∆y (4.20) 𝐹𝑤=(𝜌𝑢)𝑤∆y (4.21) 𝐹𝑛=(𝜌𝑣)𝑛∆x (4.22) 𝐹𝑠=(𝜌𝑣)𝑠∆y (4.23) Resolució numèrica de les equacions de Navier-Stokes 37 𝑃𝑒=F𝑒 𝐷𝑒 (4.24) 𝑃𝑤=F𝑤 𝐷𝑤 (4.25) 𝑃𝑛=F𝑛 𝐷𝑛 (4.26) 𝑃𝑠=F𝑠 𝐷𝑠 (4.27) on cal tenir en compte que:  Els subíndexs E i e indiquen al node i a la cara est, respectivament.  L’operador max(a,b) proporciona el valor màxim entre dos valors.  Les variables amb el superíndex 0 indiquen el valor de la variable en l’instant immediatament anterior, mentre que les que no en tenen indiquen que el seu valor és el corresponent a l’instant t tractat.  El terme font està linealitzat de la forma 𝑆𝑃𝜙𝑃+𝑆𝐶. A més, el terme 𝑆𝑃 cal que sigui negatiu per assegurar que el coeficient 𝑎𝑃 sigui del mateix signe que els altres (positius).  La superfície travessada per una cara és de valor “∆𝑥·1” o “∆y·1”, mentre que el volum del volum de control és de valor “∆𝑥∆𝑦·1”, ja que es considera profunditat unitat.  La variable dPE indica la distància entre els nodes P i E, per exemple.  El coeficient A(|P|) és el coeficient del node i+1 abans exposat, i que depèn del Peclet local a la cara corresponent. La tria de la funció és el que determinarà l’esquema que s’emprarà per a resoldre l’esquema. Aquests esquemes es poden trobar a la Taula 4.1. Resolució numèrica de les equacions de Navier-Stokes 38 Taula 4.1: Funció A(|P|) per diferents esquemes. [1] D’acord amb els esquemes que recomana Patankar [1], els esquemes emprats en els següents casos seran el Power-Law i l’Exponential. Aquest últim és la solució exacta, mentre que el primer és l’aproximació quadràtica de la solució exacta. 4.4 Cas I: Flux unidimensional amb variació de la variable en la mateixa direcció del flux 4.4.1 Introducció El cas plantejat es troba representat a la Figura 4.3. Figura 4.3: Flux unidimensional amb variació de la variable en la direcció del flux. [4] El camp de velocitats és el següent: 𝑢(𝑥,𝑦)=𝑈0 (4.28) 𝑣(𝑥,𝑦)=0 (4.29) Per aquest cas estacionari es coneix la solució analítica, que no depèn del valor de la velocitat: ∅−∅0 ∅𝐿−∅0=exp(𝑃𝑥/𝐿)−1 exp(𝑃)−1 (4.30) Resolució numèrica de les equacions de Navier-Stokes 39 4.4.2 Resultats Aquí es presenten els resultats de la simulació d’aquest 1r cas. Figura 4.4: Resultats per Peclet 1. Figura 4.5: Resultats per Peclet -10. Tal i com es pot apreciar, per diferents números de Peclet la simulació numèrica és idèntica a l’analítica. Això és perquè s’ha emprat l’esquema Power-law, que és l’aproximació polinòmica de la solució exacta exponencial. Com es pot observar a la Figura 4.4, per Peclets baixos, la solució en tot el domini va variant entre el valor imposat a l’entrada i a la sortida, ∅0=1 i ∅𝐿=0, respectivament. No obstant, en la Figura 4.5 de Peclet=-10, és un valor alt de signe negatiu, per la qual cosa el valor tendeix a ser el valor de l’extrem dret, ja que la convecció va en el sentit oposat a l’eix x. Resolució numèrica de les equacions de Navier-Stokes 40 4.5 Cas II: Flux unidimensional amb variació de la variable en direcció perpendicular al flux 4.5.1 Introducció El cas plantejat es troba representat a la Figura 4.6 Figura 4.6: Flux unidimensional amb variació de la variable en direcció perpendicular al flux. [4] El camp de velocitats és el següent: 𝑢(𝑥,𝑦)=0 (4.31) 𝑣(𝑥,𝑦)=𝑉0 (4.32) Per aquest cas es coneix la solució analítica, que no depèn del valor de la velocitat ni del Peclet: ∅=∅0+∅𝐿−∅0 𝐿𝑥 (4.33) Resolució numèrica de les equacions de Navier-Stokes 41 4.5.2 Resultats Aquí es presenten els resultats de la simulació d’aquest 2n cas. Figura 4.7: Resultats del 2n cas. Tal i com es pot apreciar, la simulació és idèntica als valors obtinguts amb la solució analítica. 4.6 Cas III: Flux diagonal 4.6.1 Introducció Aquest cas tracta d’un flux contingut en la diagonal del domini rectangular, i es troba representat a la Figura 4.8. Figura 4.8: Cas de flux diagonal. [4] Resolució numèrica de les equacions de Navier-Stokes 48 A la Figura 4.19 es pot observar el mapa de la variable general. Figura 4.19: Mapa de la variable general per 𝝆𝚪 ⁄=𝟏 𝟎𝟎𝟎. Finalment es tractarà el cas de 𝜌Γ ⁄=1 000 000. Figura 4.20: Resultats a la sortida del flux solenoïdal per 𝝆𝚪 ⁄=𝟏 𝟎𝟎𝟎 𝟎𝟎𝟎. En aquest últim cas succeeix igual que en l’anterior: la mida de la malla importa i a mesura que s’augmenta la solució s’apropa cada vegada més a la de referència. A la Figura 4.21 es pot observar el mapa de la variable general. Resolució numèrica de les equacions de Navier-Stokes 49 Figura 4.21: Mapa de la variable general per 𝝆𝚪 ⁄=𝟏 𝟎𝟎𝟎 𝟎𝟎𝟎. Per última instància, es mostrarà les línies de corrent del problema. A [3] es facilita la següent expressió per la funció de corrent: 𝜓=−(1−𝑥2)(1−𝑦2) (4.42) D’aquesta manera, emprant l’anterior equació, s’han calculat les línies de corrent, que no són més que el mapa de iso funcions de corrent. Es representa en la Figura 4.22. Figura 4.22: Línies de corrent del flux solenoïdal. Resolució numèrica de les equacions de Navier-Stokes 50 4.8 Comparativa de solvers A l’acabar el capítol de la convecció-difusió, es va acabar d’implementar el solver del Line-by-Line. D’aquesta manera, a continuació es realitzarà una comparativa entre aquest i el mètode Gauss-Seidel emprat fins ara. Per tal de realitzar aquesta comparativa, s’emprarà l’últim cas, el del flux solenoïdal, amb el cas intermig de 𝜌Γ ⁄=1 000, ja que tant la convecció com la difusió són importants. Figura 4.23: Comparativa entre els solvers de Gauss-Seidel i Line-by-Line. El resultat obtingut d’aquesta comparativa es representa en la Figura 4.23, on clarament es pot observar com el mètode Line-by-Line necessita de menys iteracions que el GaussSeidel per obtenir una solució amb la mateixa precisió. A més, a l’augmentar la malla, la tendència ha necessitar més iteracions pel Gauss-Seidel es dispara mentre que l’increment d’iteracions necessàries del Line-by-Line no és tan gran. Resolució numèrica de les equacions de Navier-Stokes 51 5 Capítol III: Equacions de Navier-Stokes 5.1 Introducció teòrica Les equacions de Navier-Stokes ja es van presentar en la introducció d’aquest estudi. Cal comentar, però, que el cas que es tractarà posteriorment no necessita de l’equació de l’energia, és per això que en aquesta introducció teòrica s’ometrà. Les hipòtesis bàsiques de treball seran les següents:  Convecció forçada  Fluid Newtonià  Propietats termofísiques constants  Dissipació viscosa negligible  Flux incompressible  Flux de viscositat constant Així les equacions de conservació de la massa i de quantitat de moviment es poden reescriure com: ∇·𝑣=0 (5.1) 𝜌𝜕𝑣 𝜕𝑡+𝜌(𝒗·𝛁)𝒗=−𝛁𝑝+𝜇∇2𝒗 (5.2) Quant a la tècnica emprada per resoldre les equacions de Navier-Stokes, serà el Mètode de Pas Fraccionat (FSM, de l’anglès Fractional Step Method). La base teòrica és el Teorema de Helmholtz-Hodge, el qual anuncia que un camp vectorial té una descomposició única en un camp de gradient pur i un vector de divergència nul·la. Això es pot aplicar ja que el flux és incompressible, i per tant el gradient de la velocitat és nul (5.1), mentre que la pressió és un camp de gradient pur. Retornant a les equacions, l’equació (5.2) es pot reescriure com: 𝜌𝜕𝑣 𝜕𝑡=𝑅(𝑣)−𝛁𝑝 (5.3) On: 𝑅(𝑣)=−𝜌(𝒗·𝛁)𝒗+𝜇𝛁2𝒗 (5.4) Resolució numèrica de les equacions de Navier-Stokes 52 Si s’integren les equacions de conservació de massa (5.1) i de quantitat de moviment (5.4) respecte el temps: ∇·𝑣𝑛+1=0 (5.5) 𝜌𝑣𝑛+1−𝑣𝑛 ∆𝑡 =32𝑅(𝑣𝑛)−12𝑅(𝑣𝑛−1)−𝛁𝑝𝑛+1 (5.6) En aquest estadi, si s’aplica el teorema de Helmholtz-Hodge, es pot obtenir la següent descomposició única: 𝑣𝑃=𝑣𝑛+1+∆𝑡 𝜌𝛁𝑝𝑛+1 (5.7) On 𝑣𝑃 és la velocitat predictora, una solució aproximada. Introduint aquesta descomposició de (5.7) a l’equació (5.6) de quantitat de moviment, es pot obtenir una equació de velocitat projectada: 𝑣𝑃=𝑣𝑛+∆𝑡 𝜌[32𝑅(𝑣𝑛)−12𝑅(𝑣𝑛−1)] (5.8) Per altra banda, si s’aplica l’operador divergència sobre la descomposició única (5.7), s’obté una equació per la pressió: l’equació de Poisson. 𝛁𝟐𝑝𝑛+1=𝜌 ∆𝑡∇·𝑣𝑃 (5.9) Finalment, la variable que es vol calcular 𝑣𝑛+1 s’obté reescrivint (5.7): 𝑣𝑛+1=𝑣𝑃−∆𝑡 𝜌𝛁𝑝𝑛+1 (5.10) D’aquesta manera es pot resumir el FSM en 5 passos: 1) Evaluació de 𝑹(𝒗𝒏) 2) Càlcul de la velocitat predictora 𝒗𝑷 3) Resolució del sistema d’equacions del mapa de pressions 4) Obtenció de la velocitat desitjada 𝒗𝒏+𝟏 5) Nou time-step En aquest mètode es poden donar solucions sense sentit físic on el gradient de pressions es desacobli del camp de velocitats (per més informació vegeu [5]). Resolució numèrica de les equacions de Navier-Stokes 53 Per evitar això, s’emprarà una malla diferent, una malla desplaçada (en anglès, staggered mesh), que és bàsicament que els nodes de càlcul de les propietats convectives estan desplaçats respecte els nodes de la malla principal on es calculen totes les variables escalars tals com la pressió. Figura 5.1: Malla desplaçada. [5] A la Figura 5.1 s’esquematitza la malla desplaçada, on la velocitat horitzontal es denomina com u i la velocitat vertical com v. Per altra banda, respecte a la malla, en el codi s’ha implementat una malla no uniforme de concentració hiperbòlica amb la metodologia exposada en [5]. Aleshores, segons el factor d’estretament k que es pot triar, la malla pot ser uniforme o no uniforme amb una forta concentració en els vèrtexs del problema que s’exposarà posteriorment. 5.2 Avaluació dels passos A continuació s’exposarà tant les equacions de discretització, com d’una manera més general, el mètode per calcular i obtenir totes les variables de cadascun dels passos del FSM. En primera instància, pel càlcul de R, s’empra una equació per la velocitat u i una per la velocitat v que són del tot anàlogues: 𝑅(𝑢)Ω𝑥𝑃=−[𝑚󰇗𝑒𝑢𝑒+𝑚󰇗𝑤𝑢𝑤+𝑚󰇗𝑛𝑢𝑛+𝑚󰇗𝑠𝑢𝑠] +[𝜇𝑒𝑢𝐸−𝑢𝑃 𝑑𝑃𝐸 𝐴𝑒−𝜇𝑤𝑢𝑃−𝑢𝑊 𝑑𝑃𝑊 𝐴𝑤 +𝜇𝑛𝑢𝑁−𝑢𝑃 𝑑𝑃𝑁 𝐴𝑛−𝜇𝑠𝑢𝑃−𝑢𝑆 𝑑𝑃𝑆 𝐴𝑠] (5.11) 𝑅(𝑣)Ω𝑦𝑃=−[𝑚󰇗𝑒𝑣𝑒+𝑚󰇗𝑤𝑣𝑤+𝑚󰇗𝑛𝑣𝑛+𝑚󰇗𝑠𝑣𝑠] +[𝜇𝑒𝑣𝐸−𝑣𝑃 𝑑𝑃𝐸 𝐴𝑒−𝜇𝑤𝑣𝑃−𝑣𝑊 𝑑𝑃𝑊 𝐴𝑤 +𝜇𝑛𝑣𝑁−𝑣𝑃 𝑑𝑃𝑁 𝐴𝑛−𝜇𝑠𝑣𝑃−𝑣𝑆 𝑑𝑃𝑆 𝐴𝑠] (5.12) Resolució numèrica de les equacions de Navier-Stokes 54 On Ω𝑃 és el volum del volum de control desplaçat de la velocitat corresponent, i els fluxos 𝑚󰇗 volumètrics s’avaluen amb interpolacions de conservació de massa, mentre que la propietat de transport 𝑢 o 𝑣 s’avalua amb esquemes numèrics convectius. El 2n pas és de càlcul directe de la velocitat predictora 𝑣𝑃 amb l’equació (5.8) un cop s’ha calculat la R. El 3r pas és on es gasta més potència computacional ja que cal resoldre un sistema d’equacions definit per l’equació de Poisson. La discretització porta a la següent equació de discretització: 𝑎𝑃𝑝𝑃𝑛+1=𝑎𝐸𝑝𝐸𝑛+1+𝑎𝑊𝑝𝑊 𝑛+1+𝑎𝑁𝑝𝑁 𝑛+1+𝑎𝑆𝑝𝑆𝑛+1+𝑏𝑃 (5.13) On: 𝑎𝐸=A𝑒 𝑑𝑃𝐸 (5.14) 𝑎𝑊=A𝑤 𝑑𝑃𝑊 (5.15) 𝑎𝑁=A𝑛 𝑑𝑃𝑁 (5.16) 𝑎𝑆=A𝑠 𝑑𝑃𝑆 (5.17) 𝑎𝑃=𝑎𝐸+𝑎𝑊+𝑎𝑁+𝑎𝑆 (5.18) 𝑏𝑃=−1 Δ𝑡[(𝜌𝑢𝑃)𝑒A𝑒−(𝜌𝑢𝑃)𝑤A𝑤+(𝜌𝑣𝑃)𝑛A𝑛−(𝜌𝑣𝑃)𝑠A𝑠] (5.19) On aquesta última equació seria 0 per algoritmes com el SIMPLE on el terme 𝑏𝑃 contindria l’equació de la conservació de la massa , però no es compleix això en aquest cas, ja que les velocitats emprades son les predictores i no les reals. Finalment, en el 4t pas, es pot obtenir el valor desitjat de les velocitats en l’instant de càlcul mitjançant l’equació (5.10). Quant a la nomenclatura, els nodes de la malla principal A i B van en l’ordre de l’abecedari en el sentit creixent de la coordenada x i/o y. 𝑢𝑃𝑛+1=𝑢𝑃𝑃−∆𝑡 𝜌·𝑝𝐵𝑛+1−𝑝𝐴𝑛+1 𝑑𝐴𝐵 (5.20) 𝑣𝑃𝑛+1=𝑣𝑃𝑃−∆𝑡 𝜌·𝑝𝐵𝑛+1−𝑝𝐴𝑛+1 𝑑𝐴𝐵 (5.21) Resolució numèrica de les equacions de Navier-Stokes 55 En última instància cal decidir si cal un nou time-step, una nova integració temporal i un nou salt en la coordenada del temps. El ∆𝑡 emprat serà adaptatiu i es calcularà en cada iteració segons la condició de CourantFriedrich-Levy, on es triarà el mínim (5.24) entre els time-steps en relació a la convecció (5.22) i a la difusió (5.23): Δ𝑡𝑐𝑜𝑛𝑣=min(0.35Δ𝑥 |𝑣|) (5.22) Δ𝑡𝑑𝑖𝑓=min(0.20ρΔ𝑥Δ𝑦 𝜇) (5.23) Δ𝑡=min(Δ𝑡𝑐𝑜𝑛𝑣 ,Δ𝑡𝑑𝑖𝑓) (5.24) Per altra banda, en el cas que s’analitzarà, es voldrà assolir l’estat estacionari. Aquesta condició s’avaluarà mitjançant la següent expressió: max(| 𝑣    𝑛+1− 𝑣    𝑛| Δ𝑡 ) (5.25) On es busca la màxima diferència dels valors de les velocitats en l’instant actual i l’immediatament anterior. Quan l’expressió (5.25) és menor a un criteri de convergència, per exemple, menor a 10−3, es considera que s’ha assolit l’estat estacionari i s’acaba la resolució del problema. Resolució numèrica de les equacions de Navier-Stokes 56 5.3 Introducció cas: Driven-Cavity El cas tractat és el d’una cavitat quadrada on les parets són estàtiques menys la de dalt, que es mou a una velocitat donada horitzontal 𝑢𝑟𝑒𝑓. Figura 5.2: Esquema del cas tractat. [5] D’aquesta manera les condicions de contorn queden definides a tots els voltants com: 𝑣=0 (5.26) 𝜕𝑝 𝜕𝑛=0 (5.27) I en la paret superior: 𝑢=𝑢𝑟𝑒𝑓 (5.28) I per la resta de parets: 𝑢=0 (5.29) Per comparar els resultats de les simulacions, s’agafaran com a referència els de l’article de [7]. Aquests resultats donen alguns valors de la velocitat u per la línia vertical central de l’esquema de la Figura 5.2 i la velocitat v per la línia horitzontal central per diferents números de Reynolds. Resolució numèrica de les equacions de Navier-Stokes 57 5.4 Esquemes de convecció-difusió En aquest capítol s’empraran altres esquemes de convecció-difusió que els emprats en el capítol anterior que eren els recomanats per Patankar [1]. Així, es pot estudiar-ne més. D’aquesta manera, s’empraran altres esquemes trobats a [3], l’UDS (Upwind Direct Scheme) i el CDS (Central Difference Scheme), com també un esquema de més alt ordre trobat a [8], l’SMART. L’UDS és un esquema de 1r ordre, mentre que el CDS és de 2n ordre i l’SMART de 2n4t ordre. 5.5 Estudi de convergència Atès que la potència computacional de l’ordinador en què es duu a terme l’estudi no és elevada, es farà un petit estudi de convergència computant l’error entre la simulació i els valors de referència, a més de veure el temps que triga a fer-se les simulacions. La situació simulada ha estat amb un Reynolds de 100 i l’esquema emprat han estat l’UDS i l’SMART. Figura 5.3: Error computat en l’estudi de convergència. Resolució numèrica de les equacions de Navier-Stokes 64 Figura 5.16: Mapa del mòdul de la velocitat per Reynolds 3200. Figura 5.17: Línies de corrent per Reynolds 3200. 5.9 Cas D: Reynolds 5000 En aquest cas, el Reynolds ja és elevat i el problema està en zones de transició on les inestabilitats del flux aniran conduint cap a la turbulència. D’aquesta manera, és més il·lustratiu mostrar aquests resultats on s’evidencia aquest fet, desmarcant-se dels resultats exposats anteriorment Cal recordar quin era el criteri per considerar que la simulació havia assolit l’estat estacionari. L’equació era la (5.25), que s’anomenarà diferència en les següents figures. Resolució numèrica de les equacions de Navier-Stokes 65 Per una malla de 50x50, la diferència en la simulació es representa en la Figura 5.18. Figura 5.18: Diferència per Reynolds 5000 i malla 50x50. Cal comentar que en les abscisses es troba el temps simulat del fenomen físic, no és el temps que triga l’ordinador a resoldre-ho. A part, el gràfic és semilogarítmic per apreciar millor tot l’espectre. Es pot apreciar com la diferència fluctua molt i li costa arribar al valor de convergència imposat de 10−3. En aquesta anàlisi s’evidencia la presència d’inestabilitats. En la mateixa línia, a posteriori es va simular amb una malla més densa, de 80x80. Per falta de temps, la simulació no va arribar a la convergència. Es representa la diferència en la Figura 5.19. Figura 5.19: Diferència per Reynolds 5000 i malla 80x80. Resolució numèrica de les equacions de Navier-Stokes 66 Tal i com es pot observar, es resistia a baixar cap al valor de 10−3 i oscil·lava. 5.10 Optimització del codi De manera que el codi funcionés més eficientment, es va fer una modificació en aquest que va resultar en què les simulacions, de mitjana, trigaven 3 vegades menys. Quan es crida una funció, l’input pot passar-se a la funció en forma de còpia (l’ordinador es fa una còpia) o per referència. En el 1r cas, el que es fa amb la variable no afecta en el codi principal, mentre que en el 2n cas, per referència, en què no es fa una còpia, si en la funció es modifica la variable, en el codi principal també queda modificada. És per això que molts inputs es passaven amb còpia, per no modificar-los. No obstant, també es poden passar per referència constant, i no es modifica la variable al programa principal, i l’ordinador tampoc fa una còpia. Aquesta optimització del codi feta a finals d’aquest estudi, millora el rendiment dels codis. Resolució numèrica de les equacions de Navier-Stokes 67 6 Conclusions i línies futures Aquest estudi ha permès desenvolupar una base sobre els mètodes bàsics emprats en el camp del CFD. Això és important ja que no sempre es pot realitzar una formació d’aquest tipus pels estudiants d’enginyeria en general. En general, l’abast d’aquest estudi i els requeriments han estat complerts. La limitació computacional de l’ordinador de què es disposa ha jugat el seu paper, però com podria passar-li a qualsevol investigador en problemes més complexos. Els problemes tractats han aportat diversos aspectes clau que ajuden a l’enteniment general del camp del CFD. En trets generals, les solucions obtingudes han estat satisfactòries. Sempre hi ha hagut errors, però s’han pogut minimitzar d’una manera o d’una altra. Un aspecte important és l’experiència que s’obté en verificar que els codis estiguin lliures d’errors, ja que trobar un error en un codi és complicat, i a poc a poc es va millorant en solucionar-los. Els següents passos serien provar casos amb fluxos externs i geometries més complicades o la resolució de problemes tridimensional. Després, una primera aproximació als fenòmens de la turbulència seria un gran pas que es podria anar desenvolupant. Resolució numèrica de les equacions de Navier-Stokes 68 7 Bibliografia [1] Suhas V. Patankar. Numerical Heat Transfer and Fluid Flow. Hemisphere Publishing Corporation, McGraw-Hill Book Company, 1980. [2] Centre Tecnològic de Transferència de Calor (CTTC), Universitat Politècnica de Catalunya. A Two-dimensional Transient Conduction Problem. (PDF). [3] Centre Tecnològic de Transferència de Calor (CTTC), Universitat Politècnica de Catalunya. Convection-difusion equations CTTC manuscript. (PDF). [4] Centre Tecnològic de Transferència de Calor (CTTC), Universitat Politècnica de Catalunya. Convection-difusion exercises. (PDF). [5] Centre Tecnològic de Transferència de Calor (CTTC), Universitat Politècnica de Catalunya. Fractional Step Method: Staggered and Collocated Meshes. (PDF). [6] F.X. Trias, O. Lehmkuhl, A self-adaptive strategy for the time integration of NavierStokes equations, Numerical Heat Transfer, Part B: Fundamentals 60 (2), 116-134, 2011. [7] Ghia et al. High-Re Solutions for Incompressible Flow Using Navier-Stokes Equations and a Multigrid Method. Journal of Computational Physics 48, 387-411 (1982). [8] Darwish, M. S. and Moukalled, F. H. NORMALIZED VARIABLE AND SPACE FORMULATION METHODOLOGY FOR HIGH-RESOLUTION SCHEMES, Numerical Heat Transfer, Part B: Fundamentals,26:1,79 – 96 (1994).