Modeliranje strujanja nestlačivog fluida

Izvor: Hrvatska internetska enciklopedija
Skoči na:orijentacija, traži

Problem opisa i analize strujanja sveprisutan je u tehničkoj praksi, od cjevovoda i ventilacijskih sustava do utjecaja vjetra na građevine i strujanja u mlaznom motoru. S napretkom tehologije rastu i zahtjevi za točnost analiza strujanja fluida. Kako su eksperimentalna mjerenja skupa i dugotrajna poželjno je da se proces optimiranja u što većem dijelu obavi prije samih ispitivanja. Uzimajući u obzir da su danas računala velike procesorske moći lako dostupna, metodama računalne dinamike fluida provode se analize strujanja unutar cijelih gradova, a zahtjevi za veće i detaljnije simulacije neprestano rastu. U nastavku je dan pregled osnovnih matematičkih i fizikalnih modela korištenih za simulaciju strujanja nestlačivog fluida s osvrtom na generalni pristup rješavanju problema strujanja i korištenih numeričkih metoda.

Osnovne jednadžbe za opis strujanja fluida

Strujanje fluida opisano je s pet osnovnih zakona fizike:

U općem slučaju moguće je koristiti svih pet zakona za opisivanje strujanja, no to generalno nije potrebno. Ako se pretpostavi da na materijalni volumen ne djeluju površinski i volumenski momenti to uvjetuje da je tenzor naprezanja simetričan pa je zakon očuvanja momenta količine gibanja implicitno zadovoljen i nije potreban za opisivanje strujanja. Ova pretpostavka je često dovoljno točna pa se ovaj zakon koristi samo u slučajima kada je prisutna rotacija fluida oko neke osi, kao na primjer kod strujanja u turbostrojevima. Drugi zakon termodinamike uvodi jednu novu nepoznanicu u sustav jednadžbi, entropiju, koja se ne pojavljuje u drugim jednadžbama pa se ova jednadžba moze rješavati samostalno nakon što su dobivena polja ostalih fizikalnih veličina i govori o fizikalnosti procesa strujanja.

Iz opće integralne formulacije za promjenu nekog fizikalnog svojstva unutar kontrolnog volumena:

[math]\displaystyle{ \int_{V_{KV}} \frac{\partial \rho\phi}{\partial t}dV + \int_{S_{KV}} \left(\rho v_j \phi - \Gamma\frac{\partial \phi}{\partial x_j} \right) n_j dS = \int_{V_{KV}}S_\phi dV }[/math]

  • [math]\displaystyle{ \rho }[/math] - gustoća fluida
  • [math]\displaystyle{ \phi }[/math] - proizvoljno specifično fizikalno svojstvo
  • [math]\displaystyle{ v_j }[/math] - brzina strujanja fluida
  • [math]\displaystyle{ \Gamma }[/math] - faktor difuzije
  • [math]\displaystyle{ S_{KV} }[/math] - kontrolna površina
  • [math]\displaystyle{ V_{KV} }[/math] - kontrolni volumen

prvi član s lijeve strane jednadžbe predstavlja lokalnu promjenu fizikalnog svojstva u vremenu, a drugi član predstavlja konvektivni i difuzni protok fizikalnog svojstva kroz kontrolnu površinu dok desna strana jednadžbe predstavlja volumenski izvor ili ponor fizikalnog svojstva.

Uz dodatak dopunskih jednadžbi specifičnih za fluid i režim strujanja:

[math]\displaystyle{ \Sigma_{ij} = \mu\left(\frac{\partial v_j}{\partial x_i}+\frac{\partial v_i}{\partial x_j} \right)-\frac{2}{3}\mu\frac{\partial v_k}{\partial x_k}\delta_{ij} }[/math]

  • [math]\displaystyle{ \Sigma_{ij} }[/math] - tenzor viskoznog naprezanja
  • [math]\displaystyle{ \mu }[/math] - dinamička viskoznost
  • [math]\displaystyle{ \delta_{ij} }[/math] - Kroneckerov delta
  • Fourierov zakon o provodenju topline

[math]\displaystyle{ q_i = - \lambda \frac{\partial T}{\partial x_i} }[/math]

  • [math]\displaystyle{ q_i }[/math] - gustoća toplinskog toka
  • [math]\displaystyle{ \lambda }[/math] - koeficijent toplinske provodnosti
  • [math]\displaystyle{ T }[/math] - temperatura

te supstituciju proizvoljnog specifičnog fizikalnog svojstva mogu se dobiti integralni oblici zakona očuvanja mase, količine gibanja i energije:

  • zakon očuvanja mase

[math]\displaystyle{ \frac{\partial \rho}{\partial t} = - \frac{\partial \left( \rho v_j \right)}{\partial x_j} }[/math]

  • zakon očuvanja količine gibanja

[math]\displaystyle{ \frac{\partial \left( \rho v_i \right)}{\partial t} = - \frac{\partial \left( \rho v_j v_i \right)}{\partial x_j} + \rho f_i\it - \frac{\partial p}{\partial x_i} + \frac{\partial \Sigma_{ji}}{\partial x_j} }[/math]

  • zakon očuvanja energije

[math]\displaystyle{ \frac{\partial }{\partial t} \left[ \rho \left( \frac{v^{2}}{2} + c_v T \right) \right] = \frac{\partial }{\partial x_j} \left[ \rho v_j \left( \frac{v^{2}}{2} + c_v T \right) \right] + \rho f_i\it v_i - \frac{\partial \left( p v_i \right)}{\partial x_i} + \frac{\partial \left( \Sigma_{ji} v_i \right)}{\partial x_j} + \frac{\partial }{\partial x_j} \left( \lambda \frac{\partial T}{\partial x_i} \right) }[/math]

Za izotermno strujanje nestlačivog fluida jednadžba očuvanja energije svodi se na jednadžbu očuvanja kinetičke energije odnosno jednadbu očuvanja količine gibanja skalarno pomnoženu s brzinom strujanja te kao takva ne pridonosi rješavanju problema strujanja (u slučaju da postoji izmjena topline služi za određivanje temperaturnog polja) pa se kao takva može odbaciti. Jednadžbe očuvanja mase i količine gibanja mogu se dodatno pojednostaviti za nestlačivi fluid, te se grupno nazivaju Navier-Stokes jednadžbe:

[math]\displaystyle{ \frac{\partial v_j}{\partial x_j} = 0 }[/math]

[math]\displaystyle{ \frac{\partial \left( \rho v_i \right)}{\partial t} = - \frac{\partial \left( \rho v_j v_i \right)}{\partial x_j} + \rho f_i\it - \frac{\partial p}{\partial x_i} + \mu \frac{\partial^{2} v_i}{\partial x_j \partial x_j} }[/math]

Jednadžba tlaka

Specifičan problem kod nestlačivog strujanja vidljiv je iz gore danih jednadžbi: ne postoji zasebna jednadžba za polje tlaka! Iz jednadžbe očuvanja količine gibanja moguće je iterativnim metodama izračunati polje brzina no potrebno je poznavati polje tlaka, odnosno gradijent tlaka koji generalno nije poznat. Jedan način rješavanja ovog problema je konstruirati takvo polje tlaka koje će osigurati očuvanje mase. Uzimanjem divergencije gore dane jednadžbe očuvanja količine gibanja te uvrštenjem jednadžbe očuvanja mase dobiva se:

[math]\displaystyle{ \frac{\partial }{\partial x_i} \left( \frac{\partial p}{\partial x_i} \right) = - \frac{\partial }{\partial x_i} \left[ \frac{\partial }{\partial x_j} \left( \rho u_i u_j - \Sigma_{ij} \right) \right] + \frac{\partial \left( \rho f_i\it \right)}{\partial x_i} }[/math]

Za slučaj konstantne viskoznosti i uniformnog polja masenih sila ova jednadžba dodatno se pojednostavljuje:

[math]\displaystyle{ \frac{\partial }{\partial x_i} \left( \frac{\partial p}{\partial x_i} \right) = - \frac{\partial }{\partial x_i} \left[ \frac{\partial \left( \rho u_i u_j \right)}{\partial x_j} \right] }[/math]

Iz ove jednadžbe moguće je izračunati polje tlaka nekom od numeričkih metoda. Bitno je naglasiti da operator gradijenta potiče iz jednadbe očuvanja količine gibanja dok operator divergensa potiče iz jednadžbe očuvanja mase te je nužno očuvati konzistentnost ovih operatora kod diskretizacije ove jednadžbe. Dakle da bi se osiguralo očuvanje mase potrebno je diskretizirati gradijente isto kao i u jednadžbama od kuda su potekli pa se često jednadžba tlaka izvodi iz već prethodno diskretizirane jednadžbe očuvanja količine gibanja.

Vremenska i prostorna diskretizacija

Vidljivo je da su dane jednadžbe nelinearne parcijalne diferencijalne jednadžbe (doduše u integralnom obliku, iako je prevođenje u pravi diferencijalni oblik relativno jednostavno) i u općem slučaju nemaju analitičko rješenje očito je da je potrebno nekim numeričkim postupcima i aproksimacijama doći do zadovoljavajuce točnog rezultata.

Ovdje su dane 3 osnovne metode vremenske diskretizacije primjenom čijih se principa mogu konstruirati i složenije metode višeg reda točnosti.

Eksplicitna Eulerova metoda

Ovo je najjednostavnija metoda prvog reda točnosti pa kao takva ima usko područje primjene.

[math]\displaystyle{ \phi ^{n+1} = \phi ^n + f\it \left( t_n , \phi ^n \right) \Delta t }[/math]

Pozitivno svojstvo ove metode je to što je eksplicitna to jest moguče je izračunati vrijednost veličine [math]\displaystyle{ \phi^{n+1} }[/math] u novom vremenskom trenutku direktno iz već poznatih vrijednosti [math]\displaystyle{ \phi^n }[/math] iz prošlog vremenskog trenutka pa tako ovo zahtijeva manje računalne moći jer nije potrebno rješavati sustav jednadžbi kako bi se izračunao [math]\displaystyle{ \phi^{n+1} }[/math] u novom trenutku. Problem kod ove metode, kao i kod ostalih eksplicitnih metoda je nužno uzimanje relativno malog vremenskog koraka (ovisno o konkretnom problemu koji se rješava) kako bi se osigurala stabilnost metode to jest konvergencija prema nekom konkretnom rješenju. Drugi problem je naravno to što je prvog reda točnosti te se za dva puta manji korak od orginalnog, greška smanjuje samo za dva puta. Ovo je naravno nepovoljno jer je potrebno veće rafiniranje već i ovako malog koraka kako bi se dobilo rješenje zadovoljavajuće točnosti.

Implicitna Eulerova metoda

Ova metoda je poboljšanje prošle i rješava problem ograničenja koraka integracije no i ona je prvog reda točnosti.

[math]\displaystyle{ \phi ^{n+1} = \phi ^n + f\it \left( t_{n+1} , \phi ^{n+1} \right) \Delta t }[/math]

Kako je stabilnost više nije problem moguće je uzimanje proizvoljno velikog vremenskog koraka, ovo je vrlo povoljno posebno za stacionarne probleme gdje je cilj što prije doći do stacionarnog stanja uz korištenje što manje računalne moći. Mana ove metode osim sto je prvog reda točnosti je to što je veličina [math]\displaystyle{ \phi^{n+1} }[/math] u novom vremenskom trenutku definirana preko vrijednosti [math]\displaystyle{ \phi^n }[/math] iz sadašnjeg trenutka i nepoznate (nju se zeli izračunati) vrijednosti [math]\displaystyle{ \phi^{n+1} }[/math] u novom vremenskom trenutku. Dakle da bi se izračunala vrijednost [math]\displaystyle{ \phi^{n+1} }[/math] u novom trenutku potrebno je rješiti sustav jednadžbi što zahtjeva više računalne moći nego direktno računanje iz već poznatih vrijednosti. Usprkos ovim manama ova metoda je vrlo dobra i često se koristi kod računanja stacionarnih problema upravo zbog njene jednostavnosti i stabilnosti.

Crank-Nicholson metoda

Ova metoda moze se gledati kao jednaka mješavina implicitne i eksplicitne Eulerove metode te je drugog reda točnosti.

[math]\displaystyle{ \phi ^{n+1} = \phi ^n + \frac{1}{2} \left[ f\it \left( t_n , \phi ^n \right) + f\it \left( t_{n+1} , \phi ^{n+1} \right) \right] \Delta t }[/math]

Kao i implicitna Eulerova metoda i Crank-Nicholson metoda je stabilna za proizvoljni vremenski korak te ima dodatnu pogodnost, a to je da je drugog reda točnosti to jest za vremenski korak dva puta manji od orginalnog greška se smanjuje četiri puta. Uzimajući u obzir da je metoda drugog reda točnosti, stabilna za proizvoljni vremenski korak te zahtijeva marginalno vise računalne moći od implicitne Eulerove metode i relativno jednostavne implementacije očito je da je vrlo povoljna te se upravo zato vrlo često koristi kako za stacionarne tako i za nestacionarne probleme.

Osim danih, postoji vrlo širok spektar metoda vremenske diskretizacije (npr. skupina prediktor-korektor metoda) od kojih svaka ima neke prednosti i nedostatke pa je prilikom odabira metode potrebno odabrati onu koja zadovoljava zahtjeve određene konkretnim problemom.

Kao i kod vremenske diskretizacije postoje razni načini provođenja prostorne diskretizacije. Često koristen pristup kod formiranja shema prostorne diskretizacije je Taylorov razvoj u red i razne interpolacijske sheme. Ovdje su dani primjeri uzstrujne i centralne sheme diferencije koje su među najjednostavnijima no ostvaruju dobar kompromis izmedu točnosti i jednostavnosti uporabe.

Upwind difference scheme (UDS)

Uzstrujna shema diferencije ili upwind difference scheme - UDS je shema prvog reda točnosti i ima nedostatke tipične za sve sheme prvog reda točnosti to jest potrebna je relativno fina mreža kako bi se dobili zadovoljavajuće točni rezultati i producira numeričku difuziju odnosno u rješenjima se javlja difuzija čak i kad fizikalno difuzija ne postoji.

[math]\displaystyle{ \left( \frac{\partial \phi}{\partial x} \right)_i = \frac{\phi_i - \phi_{i-1}}{x_i - x_{i-1}} }[/math]

Iako se na prvi pogled shema prvog reda točnosti čini odbojnom, ova shema je u širokoj uporabi u večini komercijalnih i istrazivačkih software-a. Razlog ovome je njena vrlo jednostavna implementacija i još bitnije, konvergira na svim mrežama odnosno nikada ne daje oscilatorna rješenja. Ovo je najrelevantnija karakteristika ove sheme i jak argument za njeno korištenje.

Shema centralne diferencije (CDS)

Shema centralne diferencije ili "central difference scheme" - CDS je uvjetno (samo na nekim mrežama) shema drugog reda točnosti.

[math]\displaystyle{ \left( \frac{\partial \phi}{\partial x} \right)_i = \frac{\phi_{i+1} - \phi_{i-1}}{x_{i+1} - x_{i-1}} }[/math]

Razlozi za njeno korištenje su očiti, implementacija je relativno jednostavna i drugog je reda pa je moguće dobiti točnija rješenja na grubljim mrežama u odnosu na sheme prvog reda. Nedostatak je što je uvjetno stabilna odnosno za neke neke probleme moguće je dobiti oscilatorna rješenja, no usprkos tome i ova shema je široko korištena.

Metoda konačnih volumena

Kao što je prije spomenuto u općem slučaju nije moguće analitički rješiti problem strujanja fluida u proizvoljnoj domeni, a razlog ovome je inherentna kompleksnost relevantnih jednadžbi. Iako kompleksne, jednadžbe nisu nerješive, a jedna od metoda rješavanja upravo je metoda konačnih volumena. Metoda se zasniva na na rješavanju integralnih oblika zakona očuvanja mase, količine gibanja i energije tako da se domena strujanja podjeli na velik broj malih konačnih volumena te se jednadžbe postavljaju i rješavaju za svaki od njih. Relevantne veličine računaju se u težištima konačnih volumena te se interpoliraju kako bi se dobile vrijednosti na stranicama konačnog volumena. Volumenski i površinski integrali aproksimiraju se nekom od metoda za diskretizaciju integrala, najjednostavniji pristup je pretpostavljanje da je unutar volumena ili na nekoj od stranica relevantna veličina konstantna i aproksimira se nekom srednjom vrijednosti te se množi sa obujmom ili povrčinom stranice konačnog volumena:

  • aproksimacija površinskog integrala

[math]\displaystyle{ F_e = \int_{S_e} f\it dS \approx f\it_e \Delta S_e }[/math]

  • aproksimacija volumenskog integrala

[math]\displaystyle{ Q_P = \int_V q dV \approx q_P \Delta V }[/math]

Kvaliteta aproksimacije ovisi o finoći mreže, odnosno veličini i broju konačnih volumena i o korištenima shemama diskretizacije relevantnih fizikalnih veličina. Rezultat ove diskretizacije je prevođenje skupa diferencijalnih jednadžbi u velik broj jednostavnih algebarskih jednadžbi s usko povezanim parametrima. Ovakav set algebarskih jednadžbi rješava se nekom od metoda za rješavanje skupa algebarskih jednadžbi, no kako su parametri u jednadbama usko povezani, nelinearni i sheme diferencije generalno implicitne, direktno rješavanje ovakvog skupa jednadžbi bilo bi nepraktično, pa se obično rješavaju nekom iterativnom metodom dok se nepostigne zadovoljavajuće točan rezultat.

Razlozi zbog koji je metoda konačnih volumena generalno najprivlačniji i najprihvačeniji pristup rješavanju problema analize strujanja jest njeno relativno jednostavno prevođenje u računalni kod, inherentno osigurano očuvanje mase (dok god su površinski integrali isti na zajedničkim stranicama susjednih konacnih volumena) i dobra iskoristivost kod kompleksnih geometrija. Ove pozitivne karakteristike su naravno uvjetne i mogu biti narušene za odredene izvedbe metode npr. aproksimacije višeg reda od drugog postaju vrlo komplicirane za korištenje i prevođenje u računalni kod.

Postoje razni pristupi rješavanju jednadbi za opis strujanja i metoda konačnih volumena je samo jedna od njih (iako najčešće korištena), no razne verzije metoda konačnih elemenata i metoda konačnih razlika su također u uporabi.

Rješavanje modela strujanja

Prvi korak u rješavanju modela strujanja fluida je izbor i generacija mreže konačnih volumena. Ovaj korak je među najvažnijima jer nepravilan izbor mreže direktno utječe na kvalitetu dobivenih rezultata. Vrsta i načina generiranja mreža postoji na pretek i nerjetko se koristi zasebni software-ski paket za rješavanje ovog problema. Neki primjeri software-a za generaciju mreza su: Salome, STAR-CD, I-DEAS itd.

Nakon izbora mreže odabiru se sheme diskretizacije za jednadžbu očuvanja količine gibanja i jednadžbu tlaka te se nakon diskretizacije formiraju matrice koeficijenata te se sustav iterativno rješava. Razvijeni su razni algoritmi za rješavanje ovog sustava jednadbi, jedan od najkorištenijih je “Semi-Implicit Method for Pressure Linked Equations” - SIMPLE. Ovdje je dan pregled algoritma:

  1. Pretpostavi vrijednosti polja brzine i tlaka
  2. Rješi diskretizirane jednadžbe količine gibanja i izračunaj polja brzine
  3. Rješi diskretiziranu jednadžbu tlaka i izračunaj korekciju tlaka
  4. Izračunaj korekcije polja brzine
  5. Korigiraj polja brzine sa korekcijama brzine
  6. Korigiraj polje tlaka sa podrelaksiranom (smanjenom) korekcijom tlaka
  7. Ponavljaj 2-6 dok odstupanje ne zadovolji toleranciju
  8. Pomakni se u novi vremenski trenutak i kreni od 1

Postoje razne varijacije ovog algoritma (SIMPLER, SIMPLEC) koji se uglavnom razlikuju kod tretiranja korekcije tlaka.

Turbulencija

Turbulentno strujanje fluida je najčešći oblik strujanja fluida u prirodi (dakle i u tehničkoj praksi) i pojavljuje se pri visokim vrijednostima Reynoldsovog broja. Ovdje su dane glavne karakteristike turbulentnog toka:

  • nestacionarnost; oscilacije tlaka i brzine glavna su karakteristika turbulentnog toka
  • difuzivnost; zbog kaotičnog gubanja čestica proces mješanja je uvelike intenziviran, ovo se zove jos i turbulentna difuzija
  • disipativnost; zbog velikih oscilacija lokalno se pojavljuju veliki gradijenti brzine i tlaka te je disipacija energije također intenzivirana

Kako je turbulentno strujanje zbog po svojoj prirodi kompleksno zbog velikih lokalnih oscilacija fizikalnih veličina, proces opisivanja takvog strujanja nije ni malo jednostavan. Za razliku od modela laminarnog strujanja koji se pretežito razlikuju po načinu diskretizacije osnovnih jednadžbi, modeli turbulencije razlikuju se i po samim jednadžbama koristenim za opis strujanja. Uzimajući u obzir opseg teme turbulentnog strujanja ovdje se neće ulaziti u detaljnije analize fizikalnih i numeričkih modela turbulencije, no spomenut će se tri glavna pristupa:

  • direct numerical simulation - DNS; direktno se rješavaju Navier-Stokes jednadžbe bez dodatnih simplifikacija
  • large eddy simulation - LES; direktno se rješavaju Navier-Stokes jednadžbe za najveće vrtloge dok se manji vrtlozi aproksimiraju
  • Reynolds averaged Navier-Stokes - RANS; rješavaju se osrednjene Navier-Stokes jednadžbe uz dodatno modeliranje oscilatornih komponenti

Literatura i resursi

  • Ferziger, J.H; Peric, M; Computational Methods for Fluid Dynamics
  • Patankar, S.V; Numerical Heat Transfer and Flow
  • Virag, Z; Džijan, I; Računalna Mehanika Fluida - Predavanja