Prevajanje toplote po tanki palici

Emil Žagar

11. maj 2004

Povzetek

V sestavku obravnavamo prevajanje toplote v tanki palici, torej enodimenzionalnem objektu dolžine L. Enačbo prevajanja toplote prevedemo v brezdimenzionalno obliko in opišemo nekaj znanih metod za numerično reševanje (eksplicitna, implicitna, Crank-Nicholsonova). Teorijo podkrepimo s praktičnim primerom, v katerem numerično izračunamo porazdelitev toplote v jekleni palici ob izbranem času.

Ogledali si bomo naslednji problem: Tanka palica dolžine L ima po celi dolžini začetno temperaturo T_{0}. Konca palice vzdržujemo na temperaturah T_{0} (levi konec) in T_{1} (desni konec). Pri tem privzamemo, da je T_{0} < T_{1} in je palica dovolj tanka, da jo lahko obravnavamo kot enodimenzionalni objekt. Jasno je, da se bo toplota širila od toplejšega konca proti hladnejšemu, zato

[figure: T_{0}; T_{1}; t; x; 0; L; T_{0}; T_{0}; T_{1}]

Slika 1: Palica in robni ter začetna temperatura.

bo temperatura naraščala v obratni smeri. Lahko pričakujemo, da se bo po določenem času temperatura v palici porazdelila v neko stabilno stanje.

Snov  ρ [kg/m^{3}] c [kJ/kg^{∘}C] k [W/m^{∘}C]
aluminij  2707  0.896  204
železo  7897  0.452  73
jeklo  7833  0.465  54
baker  8954  0.383  386
opeka  1600  0.84  0.69
voda  999.8  4.225  0.566

Tabela 1: Tabela konstant za nekatere materiale.

Naj bo U(x, t) temperatura palice na oddaljenosti x od levega roba palice (glejte sliko 1) ob času t. Zaradi krajše pisave bomo včasih pisali tudi U := U(x, t). Za porazdelitev temperature v palici velja enačba

\frac{∂U}{∂t} = κ \frac{∂^{2}U}{∂x^{2}}  (1)

ter robni in začetni pogoji U(0, t) = T_{0}, U(L, t) = T_{1} in U(x, 0) = T_{0}. Konstanta κ se izraža kot κ = k/(ρ c), kjer je ρ gostota, c specifična toplota in k prevodnost palice. Nekaj vrednosti za te konstante je v tabeli 1. Enačbo (1) lahko skaliramo v

U(x, t) = T_{0} + (T_{1} − T_{0})u(ξ, τ),

kjer je x = L ξ in t = (L^{2}/κ)τ. Tako dobimo nove enačbe za brezdimenzijsko funkcijo u, in sicer:

\frac{∂u}{∂τ} = \frac{∂^{2}u}{∂ξ^{2}},  (2)

z robnimi in začetnimi pogoji u(0, τ) = 0, u(1, τ) = 1 ter u(ξ, 0) = 0. Torej je dovolj, če znamo (numerično) rešiti enačbo (2).

Rešujemo jo z metodo končnih diferenc. Časovni odvod aproksimiramo s prvo premo diferenco s korakom Δτ, krajevnega pa z drugo centralno diferenco s korakom Δξ:

\frac{∂u}{∂τ}(ξ, τ) ≈ \frac{u(ξ, τ + Δτ) − u(ξ, τ)}{Δτ},
\frac{∂^{2}u}{∂ξ^{2}}(τ, ξ) ≈ \frac{u(ξ + Δξ, τ) − 2u(ξ, τ) + u(ξ − Δξ, τ)}{(Δξ)^{2}}.

Označimo λ = \frac{Δτ}{(Δξ)^{2}} , zanemarimo napako in prejšnja izraza vstavimo v enačbo (2). Vrednosti na novem časovnem nivoju so

u(ξ, τ + Δτ) = λ u(ξ + Δξ, τ) + (1 − 2λ)u(ξ, τ) + λ u(ξ − Δξ, τ).

Opisani metodi rečemo eksplicitna metoda, saj se vrednosti na novem časovnem nivoju računajo neposredno iz vrednosti na prejšnjem. Metoda je enostavna, vendar ima to pomanjkljivost, da konvergira le za λ < 1/2 (optimalno λ = 1/6). To pomeni, da mora biti časovni korak zelo majhen.

Praktično metodo uporabimo takole: Interval [0, 1] ekvidistantno (ali kako drugače) razdelimo s točkami 0 = ξ_{0} < ξ_{1} < ⋯ < ξ_{n} < ξ_{n+1} = 1. Torej je Δξ = \frac{1}{n+1}, Δτ pa izberemo tako, da je λ = \frac{Δτ}{(Δξ)^{2}} < \frac{1}{2}. Če je u_{i,j} = u(ξ_{i}, τ_{j}), pri čemer je 0 = τ_{0} < τ_{1} < τ_{2} < … ustrezna časovna razdelitev, potem se eksplicitna metoda zapiše (glejte sliko 2):

u_{i,j+1} = λ u_{i−1,j} + (1 − 2λ)u_{i,j} + λ u_{i+1,j}, i = 1, 2, …, n.

[figure: u_{i,j+1}; u_{i−1,j}; u_{i,j}; u_{i+1,j}]

Slika 2: Shema ekplicitne metode.

Pri τ_{0} = 0 je zaradi začetnega pogoja u_{i,0} = 0, na robu ξ_{0} = 0 in ξ_{n+1} = 1 pa zaradi robnih pogojev u_{0,j} = 0 in u_{n+1,j} = 1. Zelo podobno uporabimo implicitno shemo (glejte sliko 3). Pri tej shemi moramo vrednosti na novem nivoju računati kot rešitev tridiagonalnega sistema enačb

−λ u_{i−1,j+1} + (1 + 2λ) u_{i,j+1} − λ u_{i+1,j+1} = u_{i,j}, i = 1, 2, …, n.

Njena prednost pa je, da konvergira za vsak λ.

Obstaja še tretja, Crank-Nicholsonova metoda (glejet sliko 4). Gre za kombinacijo eksplicitne in implicitne sheme. Spet moramo vrednosti na novem nivoju računati kot rešitev tridiagonalnega sistema

−λ u_{i−1,j+1} + 2(1 + λ) u_{i,j+1} − λ u_{i+1,j+1} =λ u_{i−1,j} + 2(1 − λ) u_{i,j} + λ u_{i+1,j},
i = 1, 2, …, n.

Prednost te metode je, da je napaka reda O((Δτ)^{2} + (Δξ)^{2}).

[figure: u_{i,j+1}; u_{i−1,j+1}; u_{i,j}; u_{i+1,j+1}]

Slika 3: Shema implicitne metode.

[figure: u_{i,j+1}; u_{i−1,j+1}; u_{i,j}; u_{i,j}; u_{i,j}; u_{i+1,j+1}]

Slika 4: Shema Crank-Nicholsonove metode.

Primer 1 Izračunati želimo porazdelitev temperature v 0.5 metra dolgi jekleni palici, ki ima na levem robu stalno temperaturo 20^{∘}C, na desnem pa 100^{∘}C. Graf porazdelitve temperature na časovnem intervalu med 0 in 5 minut prikazuje slika 5. Graf porazdelitve temperature v palici po dveh minutah pa 6.

Literatura

[1] N. D. Fowkes in J. J. Mahony, An Introduction to Mathematical Modelling, John Wiley & Sons, New York, 1996.

Slika 5: Slika porazelitve temperature v jekleni palici. 5 Slika 6: Slika porazelitve temperature v jekleni palici po dveh dveh minutah.

6