Die Ausbreitung von Temperaturschwankungen in den Boden hinein wird durch die folgenden beiden Gleichungen beschrieben:

(1)

Numerische Methoden der Umweltphysik

Computer- ¨ Ubungen zur W¨ armeleitungsgleichung 2. Februar 2006

Die Ausbreitung von Temperaturschwankungen in den Boden hinein wird durch die folgenden beiden Gleichungen beschrieben:

F = −K

_h

∂T

∂z und ∂T

∂t = − 1 C

h

∂F

∂z

Die erste Gleichung ist als das Fourier’sche W¨ armeleitungsgesetz bekannt, die zweite Gleichung ist die dazugeh¨ orige Bilanzgleichung. K

_h

[J K

⁻¹

m

⁻¹

s

⁻¹

] bezeichnet die W¨ armeleitf¨ ahigkeit und C

_h

[J m

⁻³

K

⁻¹

] die W¨ armekapazit¨ at pro Volumen.

Kombiniert man die beiden obigen Gleichungen, so ergibt sich die W¨ armeleitungsgleichung:

∂T

∂t = 1 C

_h

∂

∂z

K

_h

∂T

∂z

In dieser ¨ Ubung soll die Ausbreitung der W¨ arme in den Erdboden hinein untersucht werden.

An der Erdoberfl¨ ache wird die Temperatur T(z = 0, t) vorgegeben werden. Der Erdboden selbst besteht aus einer vorgegebenen Zahl von Schichten, wobei in jeder Schicht die Koeffizienten K

h

und C

_h

konstant sind.

1. Installation und Ausf¨ uhrung des Matlab-Programms 1. Holen Sie die Matlab Programme von:

http://www.iac.ethz.ch/staff/schaer/Vorlesungen/NumUmwelt/index.html 2. Starten Sie Matlab.

3. Starten Sie das Programm heat.

Im Startfenster k¨ onnen Sie nun das Modell ausw¨ ahlen und die Parameter eingeben.

Mit dem Start-Button starten Sie die Berechnung, deren Visualisierung in einem neuen Fenster angezeigt wird. Mit dem Button “A” k¨ onnen Sie die analytische L¨ osung ein- und ausblenden; eine solche l¨ asst sich allerdings nur f¨ ur die ersten zwei Modelle berechnen.

2. Experimentelles Absch¨ atzen der Eindringstiefe

Die Eindringtiefe z

_e

einer Temperaturvariation in den Boden wurde in der Vorlesung hergeleitet:

z

e

= r 2κ

ω mit κ = K

h

C

_h

Dabei ist κ die Temperaturleitzahl, K

_h

[J K

⁻¹

m

⁻¹

s

⁻¹

] die W¨ armeleitf¨ ahigkeit und C

_h

[J m

⁻³

K

⁻¹

] die W¨ armekapazit¨ at pro Volumen und ω = 2π/τ die Kreisfrequenz.

1

(2)

Aufgaben:

a. Verifizieren Sie qualitativ die Abh¨ angigkeit der Eindringtiefe von der Periode τ . W¨ ahlen Sie dazu eine Integrationsdauer von 100 Tagen und Perioden von 1, 10 und 100 Tagen.

b. Testen Sie analog zu a) die Abh¨ angigkeit von K

_h

und von C

_h

.

c. Vergleichen Sie das Eindringen in einem Einschichtmodell mit demjenigen in einem Zweischichtmodell. W¨ ahlen Sie zum Beispiel die W¨ armeleitf¨ ahigkeit K

h

der oberen Schicht doppelt so gross wie diejenige der unteren Schicht und die Schichtgrenze bei z = 0.5. Testen Sie auch andere Modelle!

Bemerkung: Die Courant-Zahl, die Sie eingeben, ist definiert durch α = κ ∆t

∆x

²

mit κ = K

_h

C

_h

wobei sich K

_h

und C

_h

auf die oberste Schicht beziehen. W¨ ahlen Sie α so, dass Stabilit¨ at f¨ ur alle Schichten gew¨ ahrleistet ist.

3. Implementation des Crank-Nicolson-Verfahrens

In der Vorlesung wurde das Crank-Nicolson-Verfahren vorgestellt. Die W¨ armeleitungsgleichung im Falle von konstanten Koeffizienten wird durch die implizite Vorschrift

Φ

ⁿ⁺¹_i

− Φ

ⁿ_i

∆t = κ

2 · ( Φ

ⁿ⁺¹_i+1

− 2Φ

ⁿ⁺¹_i

+ Φ

ⁿ⁺¹_i−1

∆x

²

+ Φ

ⁿ_i+1

− 2Φ

ⁿ_i

+ Φ

ⁿ_i−1

∆y

²

)

gegeben. An den Randpunkten m¨ ussen die Werte vorgegeben werden. So resultiert f¨ ur den oberen Rand (Erdoberfl¨ ache) die Vorschrift

Φ

ⁿ⁺¹₁

− Φ

ⁿ₁

∆t = κ

2 · ( Φ

ⁿ⁺¹₂

− 2Φ

ⁿ⁺¹₁

+ Φ

ⁿ⁺¹_L

∆x

²

+ Φ

ⁿ₂

− 2Φ

ⁿ₁

+ Φ

ⁿ_L

∆y

²

)

wobei Φ

ⁿ_L

den vorgebenen Wert am oberen Rand bezeichnet. Ein entsprechender Aus- druck gilt auch f¨ ur den unteren Rand. Wie ¨ ublich f¨ ur implizite Verfahren bei linearen Modellproblemen resultiert ein lineares Gleichungsystem (α = κ∆t/∆x

²

):







2 + 2α −α

−α 2 + 2α −α

.. .

−α 2 + 2α







·







Φ₁ Φ₂ .. .

Φn







n+1

−α·







Φ_L 0 .. . 0 Φ_R







n+1

= · · ·

· · · =







2−2α α

α 2−2α α .. .

α 2−2α







·







Φ₁ Φ2 .. .

Φn







n

+α·







Φ_L 0

.. . 0 Φ_R







n

In dieser Aufgabe sollen Sie das Crank-Nicolson Verfahren implementieren. Die Temper- atur an den R¨ andern sind vorgegeben: Φ

L

= Φ

L

(t), Φ

R

= 0. F¨ ur die Implementation m¨ ussen Sie folgendes wissen:

2

(3)

– Der Kern besteht in einer Routine, die ein tridiagonales lineares Gleichungsystem numerisch l¨ ost. Sie finden im Programm bereits eine solche Routine tridiag vor.

Studieren Sie die Kommentare zu dieser Routine! Im wesentlichen geht es darum, dieser Routine die richtigen Eingabegr¨ ossen zu ¨ ubergeben, und die Ausgabegr¨ ossen richtig zu interpretieren.

– Vervollst¨ andigen Sie dann die Zeilen im Hauptprogramm! Die aktuelle Temperatur steht im Array Tnow(i), die neu berechnete soll in den Array Tnew(i) kommen. Die diffusive Courant-Zahl α steht in cour. Die Invertierungsroutine k¨ onnen Sie mit x

= tridiag (nz,a,b,c,f) aufrufen. Die Temperatur an der Erdoberfl¨ ache f¨ ur den Zeitschritt j finden Sie in Tsur(j).

Aufgaben

a. Vergleichen Sie das Crank-Nicolson Verfahren mit dem “Forward-Time, Centered- Space”-Verfahren, im Bereich der diffusiven Courant-Zahlen, wo beide Schemen stabil sind.

b. Zeigen Sie (durch Simulationen), dass das Crank-Nicolson Verfahren auch f¨ ur diffusive Courant-Zahlen > 1 stabil bleibt. Diskutieren Sie qualitativ den Fehler des Verfahrens mit zunehmender Courant-Zahl!

Anhang: Gestaggerte Gitter bei der Fluss-Methode

F¨ ur die Diskretisierung der Fourier- und der Bilanzgleichung wird ein finites Differenzenver- fahren in Zeit und Ort angewendet. Die Differentialgleichung wird auf dem diskreten Orts- und Zeitgitter

(i∆z, n∆t)

((i + 1/2)∆z, n∆t) i = 0 . . . nz, n = 0 . . . nts

gel¨ ost. Hier ist nz die Anzahl Gitterpunkte und nts die Anzahl Zeitschritte, sowie ∆z die Gitterweite und ∆t die Zeitschrittl¨ ange. Die Temperatur- und W¨ armeflussvariablen werden gestaffelt (gestaggert) in das r¨ aumliche Gitter eingepasst.

Dieses ist in nz Schichten aufgeteilt, in deren Mitten die Temperatur und auf deren Grenzen der W¨ armefluss definiert ist. Schwarze Punkte repr¨ asentieren die gegebenen Anfangs- bzw.

Randbedingungen, weisse die zu berechnenden Gr¨ ossen.

Beim finiten Differenzenverfahren werden die Ableitungen durch Differenzenquotienten angen¨ ahert.

Die Fouriergleichung wird somit mit zentrierter Differenz im Ort zu F

_i+1/2ⁿ

= −K

_i+1/2

T

_i+1ⁿ

− T

_iⁿ

∆z

und die Bilanzgleichung mit dem Forward-Schritt in der Zeit und ebenfalls zentrierter Differenz im Ort zu

T

_iⁿ⁺¹

− T

_iⁿ

∆t = − 1 C

_i

F

_i+1/2ⁿ

− F

_i−1/2ⁿ

∆z

wobei T

_iⁿ

die angen¨ aherte Temperatur am Ort i∆z und zur Zeit n∆t und F

_i+1/2ⁿ

den angen¨ aherten W¨ armefluss bei (i + 1/2)∆z und n∆t bezeichnet. Hierbei wurden die Fehlerterme

3

(4)

F_3/2⁰ T₁⁰

T₂⁰

T_nz⁰

∆t

∆z

F_nz+1/2⁰

F_1/2¹ F_1/2^nts

F_1/2⁰

t= 0 t=time

z= 0

z=z_max

Figure 1: Numerisches Rechengebiet

in O(∆z) und O(∆t) vernachl¨ assigt. Werden die Rand- und Anfangsbedingungen durch T

₀ⁿ

= T

z0

(n∆t)

F

_1/2ⁿ

= F

z0

(n∆t) F

_zn+1/2ⁿ

= F

₁

(n∆t)

T

_i⁰

= T

_t0

(i∆z) ersetzt und wird die Differenzengleichung nach T

_iⁿ⁺¹

aufgel¨ ost,

T

_iⁿ⁺¹

= T

_iⁿ

+ ∆t C

_i

(∆z)

²

K

_i+1/2

(T

_i+1ⁿ

− T

_iⁿ

) − K

_i−1/2

(T

_iⁿ

− T

_i−1ⁿ

)

so k¨ onnen alle Fl¨ usse F

_i+1/2ⁿ

und Temperaturen T

_iⁿ

(i = 1 . . . nz, n = 1 . . . nts) berechnet werden. Die Pfeile in Figur 1 veranschaulichen die Art, wie die Temperaturen und Fl¨ usse auseinander hervorgehen.

4