Software

2D Heat SIM v.1.0.0

Simulacija 2D nestacionarnog provođenja toplote

Program je besplatan

Glavni prozor aplikacije 2D Heat SIM sa temperaturnim poljem u toku simulacije
Glavni prozor aplikacije — simulacija u toku

Namjena

2D Heat SIM je besplatan edukativan software koji rješava jednačinu nestacionarnog provođenja toplote u pravougaonoj ploči i prikazuje temperaturno polje u toku simulacije. Program je edukativnog karaktera — nije zamjena za komercijalne CFD i FEA pakete, već je alat sa kojim se može vršiti izbor vremenske šeme, gustine mreže i vremenskog koraka.

Zbog toga je koristan svakome ko uči ili predaje prenos toplote i numeričke metode: temperaturno polje, granični uslovi i Fourierov broj su na jednom ekranu, pa se posljedica svake promjene ulaznog podatka vidi odmah, bez pripreme modela i bez čekanja.

Fizički problem

Posmatra se homogena izotropna ploča dimenzija Lx × Ly, bez unutrašnjih izvora toplote, sa uniformnom početnom temperaturom T0. Provođenje toplote opisuje Fourierova (difuziona) jednačina:

∂T/∂t = α (∂2T/∂x2 + ∂2T/∂y2)

gdje je α = k / (ρ cp) termička difuzivnost [m2/s]. Lijeva strana je brzina akumulacije toplotne energije u tački, a desna neto priliv toplote usljed zakrivljenosti temperaturnog polja. Odnos k prema ρcp govori koliko brzo materijal reaguje na toplotni poremećaj — zato bakar i voda, pri istoj geometriji, daju potpuno različitu dinamiku procesa.

Granični uslovi

Granični uslov se zadaje nezavisno za svaku od četiri stranice ploče, što omogućava kombinacije kakve se sreću u praksi — na primjer grijana donja ivica, izolovane bočne ivice i konvekcija na gornjoj površini.

TipUslovFizičko značenje
DirichletT = TwZadana temperatura površine — kontakt sa masivnim grijačem, hladnjakom ili izotermalnim rezervoarom
Izolovano∂T/∂n = 0Adijabatska granica — savršena izolacija ili ravan simetrije
Konvekcija−k ∂T/∂n = h (Ts − T∞)Razmjena sa fluidom poznate temperature, uz koeficijent prelaza toplote h [W/(m2K)]

Polja koja za izabrani tip nemaju smisla program automatski onemogućava i prikazuje sivom bojom — kod izolovane granice nema unosa, kod konvekcije se traže h i Tfluid, a temperatura zida ostaje neaktivna. Time se izbjegava najčešća greška u tumačenju rezultata: unesena vrijednost koja u proračun uopšte ne ulazi.

Numeričke metode

Prostorna diskretizacija je u svim slučajevima metoda konačnih razlika — centralne razlike drugog reda tačnosti na uniformnoj mreži. Drugog reda su i granični uslovi: temperatura graničnog čvora računa se iz jednostrane razlike sa tri tačke, ∂T/∂n ≈ (−3T0 + 4T1 − T2)/(2d), pa red tačnosti ne pada uz sam rub domena. Za vremensku integraciju bira se jedna od tri šeme, koje su matematički objedinjene jednim parametrom θ:

ŠemaθRed tačnostiStabilnostCijena po koraku
Eksplicitna (Forward Euler)0O(Δt)uslovna: Fo ≤ 0,5najniža — bez sistema jednačina
Crank-Nicolson1/2O(Δt2)bezuslovna (A-stabilna)rješavanje rijetkog sistema
Implicitna (Backward Euler)1O(Δt)bezuslovna (L-stabilna)rješavanje rijetkog sistema

Kod implicitnih šema matrica sistema se sastavlja i LU faktoriše samo jednom, na početku simulacije, pa se u svakom koraku rješava samo trougaoni sistem. Praktična posljedica je da veći vremenski korak znači i manje procesorsko vrijeme, a ne samo formalnu stabilnost.

Kod konvektivnih i izolovanih granica temperatura ruba se mijenja iz koraka u korak. Ako se u implicitnom koraku uzme njena vrijednost sa starog vremenskog nivoa, unosi se greška reda O(Δt) koja obara red tačnosti cijele šeme na 1 — bez obzira na θ. Zato se radi predikcija–korekcija: sistem se prvo riješi sa starom granicom, iz te predikcije se odredi granica na novom nivou, pa se riješi ponovo sa ispravno otežanim doprinosom. Kako je LU faktorizacija već gotova, drugo rješavanje je jeftino, a Crank-Nicolson time zadržava svoj drugi red i uz konvekciju.

Fourierov broj i granica stabilnosti

Program u svakom trenutku računa i prikazuje Fourierov broj:

Fo = αΔt/Δx2 + αΔt/Δy2

Za eksplicitnu šemu vrijednost se boji zeleno dok je Fo ≤ 0,5, a crveno čim se granica pređe; prije pokretanja simulacije izdaje se i upozorenje. Za Crank-Nicolson i implicitnu šemu prikazuje se napomena da granica ne važi. Ovo je, uz mogućnost da se ista postavka odmah ponovi drugom šemom, najkraći način da se vidi šta uslovna stabilnost stvarno znači.

Mogućnosti programa

  • Geometrija i mreža — dimenzije ploče Lx i Ly i koraci mreže Δx i Δy zadaju se nezavisno; računska mreža se prikazuje odmah po pokretanju programa, zajedno sa brojem čvorova, i osvježava se dok se koraci mijenjaju
  • Fourierov broj u realnom vremenu — Fo se preračunava čim se promijeni Δx, Δy, Δt ili materijal, pa se prije pokretanja vidi da li je izabrana kombinacija stabilna za eksplicitnu šemu
  • Baza materijala — metali i fluidi (bakar, čelik, vazduh, voda) te građevinski materijali (puna cigla, beton, mineralna vuna, EPS); izborom materijala popunjavaju se ρ, cp i k, a α se računa automatski
  • Vlastiti materijal — opcijom Slobodan izbor unose se ρ, cp i k za bilo koji materijal, a α = k/(ρ·cp) se preračunava dok se kuca. α se namjerno ne unosi zasebno: time je isključena fizički nesaglasna kombinacija u kojoj α ne odgovara unesenim svojstvima, a k se koristi i u konvektivnom graničnom uslovu
  • Vremenska kontrola — vremenski korak Δt i ukupno vrijeme simulacije; simulacija se može zaustaviti u bilo kom trenutku
  • Kriterijum konvergencije — simulacija se automatski prekida kada maksimalna promjena temperature između dva koraka padne ispod zadate vrijednosti ε, čime se dobija stacionarno stanje bez pogađanja potrebnog vremena
  • Provjera unosa — svako polje ima dozvoljeni opseg; neispravan unos se boji crveno, a prije pokretanja proračuna prikazuje se spisak svega što treba ispraviti, uključujući i nesaglasne kombinacije (npr. Δt veći od ukupnog vremena ili premala mreža)
  • Prikaz temperaturnog polja — popunjene konture sa 20 nivoa i preklopljenim izotermama, sa skalom temperature koja se osvježava tokom proračuna
  • Očitavanje vrijednosti mišem — pomjeranjem kursora preko dijagrama prikazuju se koordinate x, y i interpolirana temperatura u toj tački
  • Dijagram temperature u centru ploče — Tcentar u funkciji vremena, kao mjera brzine uspostavljanja ravnoteže
  • Izvoz rezultata — bira se folder za snimanje, a fajlovi dobijaju vremensku oznaku: temperaturno polje i dijagram centra kao PNG (150 DPI) te kompletno polje temperatura kao CSV, za dalju obradu u Excelu ili MATLAB-u
  • Praćenje toka — traka napretka i statusna linija sa razlogom zaustavljanja (isteklo vrijeme ili postignuta konvergencija), temperaturom u centru i dostignutim vremenom

Šta se na programu može pokazati

  • Uticaj termičke difuzivnosti — ista geometrija i isti granični uslovi, a bakar i voda daju vremena uspostavljanja ravnoteže koja se razlikuju za tri reda veličine
  • Granica stabilnosti eksplicitne šeme — postepenim povećavanjem Δt prelazi se Fo = 0,5 i rješenje vidljivo divergira
  • Cijena i korist implicitnih šema — isti zadatak riješen sa 20 umjesto 2000 koraka, uz poređenje tačnosti
  • Uticaj tipa granice — ista ploča sa izolovanom i sa konvektivnom gornjom ivicom, i razlika u stacionarnom polju
  • Konvergencija mreže — usitnjavanjem Δx i Δy prati se kako se rješenje primiče tačnom
  • Prelazak u stacionarno stanje — kada vremenski član iščezne, rješenje zadovoljava Laplaceovu jednačinu
  • Prolaz toplote kroz zid — vidi primjer ispod

Primjer: zid od pune cigle

Zid debljine 0,25 m i visine 0,50 m, unutra 25 °C, napolju vazduh 5 °C. Gornja i donja ivica izolovane (čime se problem svodi na jednodimenzionalan kroz debljinu), lijeva i desna konvektivne, sa koeficijentima prelaza po EN ISO 6946: hi = 7,7 W/m²K iznutra (Rsi = 0,13) i he = 25 W/m²K spolja (Rse = 0,04). Materijal: puna cigla, Crank-Nicolson, Δt = 60 s.

Karakteristično vrijeme takvog zida je L²/α ≈ 37 sati, pa se traži nekoliko desetina sati simuliranog vremena — simulacija se sama zaustavi kad dostigne stacionarno stanje. Rezultat se poklapa sa ručnim proračunom preko toplotnih otpora:

VeličinaRučni proračunSimulacija
Temperatura unutrašnje površine20,07 °C20,07 °C
Temperatura spoljne površine6,52 °C6,52 °C
Toplotni fluks q37,95 W/m²37,95 W/m²

(R = 0,13 + 0,25/0,70 + 0,04 = 0,527 m²K/W → U = 1,90 W/m²K.) Zamjenom cigle mineralnom vunom, uz sve ostalo isto, fluks pada na oko 3 W/m² — k određuje gubitak u ustaljenom stanju, dok α određuje samo koliko brzo se do njega stigne.

Verifikacija

Svaki od tri tipa graničnog uslova provjeren je zasebno — jer program tačan za Dirichleta ne mora biti tačan i za konvekciju.

Dirichletova granica

Ploča sa uniformnom početnom temperaturom kojoj se sve četiri granice naglo postave na zadatu temperaturu. Rješenje se dobija razdvajanjem promjenljivih, kao proizvod dva Fourierova reda. Analiza konvergencije mreže daje empirijski red tačnosti koji odgovara teorijskoj vrijednosti 2:

Δx = Δy [m]Broj čvorovaMaks. odstupanje [°C]Red tačnosti
0,02026 × 260,1664—
0,01051 × 510,04361,93
0,005101 × 1010,01101,98

Konvektivna (Robinova) granica

Ploča sa konvekcijom na sve četiri strane, pri Biotovom broju Bi = 2. Analitičko rješenje je opet proizvod dva 1D rješenja, ali sa svojstvenim vrijednostima iz transcendentne jednačine λ tan λ = Bi. I ovdje se dobija drugi red tačnosti:

Δx = Δy [m]Broj čvorovaMaks. odstupanje [°C]Red tačnosti
0,02026 × 261,4379—
0,01051 × 510,39341,87
0,005101 × 1010,10121,96

Izolovana (Neumannova) granica

Za nju analitičko rješenje nije potrebno — koristi se simetrija. Ploča grijana sa obje strane simetrična je oko svoje sredine, pa polovina domene sa izolovanom granicom na ravni simetrije mora dati isto polje kao puna ploča. Odstupanje iznosi 0,001 °C, što potvrđuje da je uslov ispravno postavljen.

Stabilnost

Ponovljeno sa vremenskim korakom koji odgovara Fo = 1,16 — dakle iznad granice stabilnosti — eksplicitna šema divergira već nakon petnaestak koraka, dok Crank-Nicolson i implicitna šema ostaju stabilne i daju fizički smislen rezultat.

Ograničenja

Da bi rezultati bili ispravno tumačeni, korisno je znati i šta program ne radi:

  • Model je dvodimenzionalan, na pravougaonom domenu i strukturiranoj mreži — nema složene geometrije ni nestrukturiranih mreža
  • Materijal je homogen i izotropan, sa konstantnim svojstvima; nema temperaturne zavisnosti k, ρ i cp, ni višeslojnih konstrukcija
  • Nema unutrašnjih izvora toplote, promjene faze, strujanja fluida ni zračenja — posmatra se čisto provođenje
  • Granični uslovi su konstantni u vremenu i uniformni duž stranice
  • Ugaoni čvorovi se računaju kao srednja vrijednost dvije procjene, po jedna iz svakog smjera; u proračun unutrašnjosti ne ulaze, ali se prikazuju i izvoze

Za probleme izvan ovih okvira potreban je pravi CFD ili FEA paket. U okviru za koji je namijenjen, program daje rezultate koji su validirani analitički.

Tehnologije

Program je pisan u Pythonu, sa otvorenim i široko korištenim bibliotekama:

Python 3.11 · PySide6 · NumPy · SciPy (rijetke matrice) · Matplotlib

Grafički interfejs koristi Qt preko biblioteke PySide6, pod licencom GNU LGPL v3. Qt biblioteke se isporučuju kao zasebni fajlovi, pa ih korisnik može zamijeniti vlastitom verzijom, kako LGPL i zahtijeva. NumPy, SciPy i Matplotlib su pod BSD licencom.

Korisnički interfejs je odvojen i uređuje se u Qt Designeru. Proračun i vizualizacija su nezavisni, pa se numerički dio može čitati i mijenjati nezavisno od interfejsa.

Cijena i sugestije

2D Heat SIM je potpuno besplatan, bez ograničenja u funkcionalnosti i bez naknadnih verzija koje se plaćaju.

Sve sugestije su dobrodošle — predlozi za nove mogućnosti, primjedbe, prijave grešaka.

Kontakt

Za pitanja, prijedloge ili prijave grešaka — direktno na mkozica@outlook.com.

Preuzimanje

Instalacioni program za Windows (64-bitni), veličine 69 MB. Sve potrebne biblioteke su uključene — Program se može instalirati i bez administratorskih prava.

Preuzmi 2D Heat SIM v.1.0.0  ↓

← Nazad na Software