FY1005/TFY4165 Termisk fysikk. Institutt for fysikk, NTNU. Våren 2014.

Løsningsforslag til øving 13

a) Med kulesymmetri blir

∇^{2} = \frac{1}{r} \frac{∂^{2}}{∂r^{2}} r,

så vi får

\frac{∂T}{∂t} = a\frac{sin kr}{r} (−D_{T} k^{2})e^{−D_{T} k^{2}t} = −D_{T} k^{2}T,
∇^{2}T = a(−k^{2})\frac{sin kr}{r} e^{−D_{T} k^{2}t} = −k^{2}T.

Følgelig er den gitte T løsning av varmeledningsligningen.

b) Grensebetingelsen T = T_{∞} for r = R innebærer at sin k_{n}R = 0,

eller

k_{n}R = nπ dvs k_{n} = nπ/R,

med n et heltall.

c) Ved å sette inn det gitte uttrykket for T i integralet for a_{m} finner en (t = 0)

a_{m} = \frac{2}{R} ∑_{n} a_{n} ∫_{0}^{R} \frac{sin k_{n}r}{r} r sin k_{m}r dr
= \frac{2}{R} ∑_{n} a_{n} \frac{1}{2} ∫_{0}^{R}\{[cos[(n − m)πr/R] − cos[(n + m)πr/R]\}dr
= \frac{2}{R} ∑_{n} a_{n} \frac{1}{2}Rδ_{nm}
= a_{m},

som skulle vises. Her har vi brukt at

∫_{0}^{R} cos(qπr/R) dr = 0

dersom q er et heltall forskjellig fra null, mens det samme integralet har verdien R dersom q = 0. Med den gitte grensebetingelsen blir koeffisientene

a_{n} = \frac{2}{R}T_{0} ∫_{0}^{R} r sin k_{n}r dr
= \frac{2T_{0}}{R} |_{0}^{R} (−\frac{r cos k_{n}r}{k_{n}} + \frac{sin k_{n}r}{k_{n}^{2}} )
= \frac{2T_{0}}{R} (\frac{−R cos nπ}{k_{n}} + 0) = \frac{2T_{0}}{k_{n}} (−1)^{n−1}
= \frac{2RT_{0}}{nπ} (−1)^{n−1}.

d) Siden

lim_{r→0} \frac{sin k_{n}r}{k_{n}r} = 1,

har en at

T (0, t) − T_{∞} = ∑_{n} a_{n}k_{n}e^{−D_{T} k_{n}^{2}t} →_{t→∞} a_{1}k_{1}e^{−D_{T} k_{1}^{2}t} = 2T_{0}e^{−D_{T} k_{1}^{2}t}.

Ved tida t = τ har en så

T (0, τ) − T_{∞} = 2T_{0}e^{−D_{T} k_{1}^{2}τ} = \frac{1}{10}T_{0}.

Løst mhp τ gir dette

τ = \frac{− ln(1/20)}{D_{T} k_{1}^{2}} ≃ 1.52 h = 5460 s.

Her har vi satt inn k_{1} = π/R, R = 0.05 m, og D_{T} = 0.0005 m^{2}/h.

Litt ekstra puslearbeid for den som liker denslags: Kontroller den oppgitte rekkeutviklingen for T ved å sette inn uttrykk for koeffisientene a_{n} og ”bølgetallene” k_{n} og regne ut

T (r, 0) − T_{∞} = \frac{1}{r} ∑_{n=1}^{∞} a_{n} sin(k_{n}r)

og vise at høyre side her blir lik T_{0}, som oppgitt i punkt c). Tips:

ln(1 + x) = ∑_{n=1}^{∞} \frac{(−1)^{n−1}x^{n}}{n}

og

sin(nα) = \frac{1}{2i}[(e^{iα})^{n} − (e^{−iα})^{n}].

e) Forslag til MATLAB-program som besvarer oppgaven (se utlagt fil ov13.m):

%FY1005/TFY4165 Termisk fysikk, oving 13, oppgave e, 2014.
%Avkjoling av vannkule med perfekt kobling til omgivende
%varmereservoar.
clear;
format long;
%Tabell med tidspunkter (i sekunder)
t=0:60:10800;
%Tabell med posisjoner (i meter), fra 1 mm til 50 mm
r=0.001:0.001:0.050;
%Kulas radius (m):
R=0.050;
%Temperaturforskjell mellom kule og omgivelser ved t=0:
T0=80;
%Temperatur i omgivende reservoar (K):
T1=293;
%Termisk diffusivitet for vann (m^2/s):
D=(25E-7)/18;
a=D*pi*pi/(R*R);
%Antall ledd som tas med i losningen:
N=1000;
%Starttabell for T(r) med like mange elementer som
%tabellen for posisjonen r:
T = linspace(T1,T1+T0,length(r));
%Neste linje setter EraseMode til xor, bra for "smooth animation", se
%http://nd.edu/~dtl/cheg258/notes/doc/tec2.5.html
p=plot(r,T,'-','EraseMode','xor');
%Sett aksegrenser [xmin xmax ymin ymax]:
axis([0 0.05 T1-3 T1+1.25*T0]);
%Nodvendig med "hold on" for fortsatt bruk av samme figur:
hold on;
xlabel('Posisjon (m)');
ylabel('Temperatur (K)');
title('Avkjoling av vannkule, animasjon');
for i = 1 : length(t),
 for k = 1 : length(r),
  T(k) = T1;
  for j = 1 : N,
  T(k)=T(k)+2*T0*(-1)^(j+1)*exp(-j*j*a*t(i))*sin((j*pi/R)*r(k))/((j*pi/R)*r(k));
  end
  %Tsentrum = T(t) i sentrum av kula:
  if k == 1,
      Tsentrum(i) = T(k);
  end
  %Tmidti = T(t) i avstand R/2 fra sentrum:
  if k == int8(length(r)*1/2),
      Tmidti(i) = T(k);
  end
  %Tlangtut = T(t) i avstand 2R/3 fra sentrum:

if k == int8(length(r)*2/3),
    Tlangtut(i) = T(k);
end
%Tytterst = T(t) paa kulas overflate:
if k == length(r),
    Tytterst(i) = T(k);
end
%T60 = T(r) ved t = 60 s:
if i == 2,
    T60 = T;
end
%T1200 = T(r) ved t = 1200 s:
if i == 21,
    T1200 = T;
end
%T3600 = T(r) ved t = 3600 s:
if i == 61,
    T3600 = T;
end
%T10800 = T(r) ved t = 10800 s:
if i == 181,
    T10800 = T;
end
end
%Plott T(r) for aktuelt tidspunkt:
set(p,'XData',r,'YData',T)
%Oppdaterer grafen i figuren:
drawnow
%Forsinker framvisningen i 0.2 s (juster etter behov):
pause(0.2);
end
hold off;
%Ber om ny figur for aa plotte T(t) for valgte avstander fra sentrum:
figure;
plot(t,Tsentrum,t,Tmidti,t,Tlangtut,t,Tytterst);
legend('r=0','r=R/2','r=2R/3','r=R');
xlabel('Tid (s)');
ylabel('Temperatur (K)');
title('Avkjoling av vannkule');
%Ber om ny figur for aa plotte T(r) for valgte tidspunkter:
figure;
plot(r,T60,r,T1200,r,T3600,r,T10800);
legend('t=1min','t=20min','t=60min','t=180min');
axis([0 0.05 T1-3 T1+1.5*T0]);
xlabel('Posisjon (m)');
ylabel('Temperatur (K)');
title('Avkjoling av vannkule');

Figur som viser T (t) for utvalgte verdier av posisjonen r: Figur som viser T (r) for utvalgte tidspunkter t: