projekt simulationstools und ihre anwendung · 2 2 modellbildung dieses kapitel befasst sich mit...

Fakultät 5Institut für Mechanik

Projekt Simulationstools und ihreAnwendung

«Analytische und numerische Untersuchungen derthermischen Entwicklung eines Planeten mit

Strahlungsrandbedingung»

vorgelegt von:

Julian Gieseler 332702Peter Kienz 341671

Marcus Lauenstein 332886Joel Panepinto 339150

Hochschullehrer: Prof. Dr. rer. nat. habil. W. H. MüllerBetreuer: M.Sc. Felix A. Reich

Technische Universität Berlin, Fakultät 5 – Institut für Mechanik,Fachgebiet für Kontinuumsmechanik und Materialtheorie

Berlin, 22. März 2017

II

Inhaltsverzeichnis

1 Motivation und Zielsetzung 1

2 Modellbildung 22.1 Massen-, Impuls- und Energiebilanz . . . . . . . . . . . . . . . . 32.2 Konstitutive Gleichungen . . . . . . . . . . . . . . . . . . . . . . 82.3 Rand- und Anfangsbedingungen . . . . . . . . . . . . . . . . . . 152.4 Entdimensionalisierung . . . . . . . . . . . . . . . . . . . . . . . . 212.5 Parameter des Erdmodells . . . . . . . . . . . . . . . . . . . . . . 232.6 Quasistationärer Deformationsprozess . . . . . . . . . . . . . . . 25

3 Analytische Lösungen 263.1 Radialsymmetrische, stationäre Wärmeleitung mit ROBIN-Rand

bedingung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 263.2 Radialsymmetrische, instationäre Wärmeleitung mit ROBIN-Rand

bedingung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 283.3 Radialsymmetrische, stationäre Wärmeleitung mit STEFAN-

BOLTZMANN-Randbedingung . . . . . . . . . . . . . . . . . . . . 353.4 Stationäres, thermoelastisches Problem mit STEFAN-BOLTZMANN-

Randbedingung . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

4 Numerische Lösungen 424.1 Schwache Form der Feldgleichungen . . . . . . . . . . . . . . . . 424.2 Stationäre Wärmeleitung mit ROBIN-Randbedingung . . . . . . . 484.3 Instationäre Wärmeleitung mit ROBIN-Randbedingung . . . . . . 494.4 Stationäre Wärmeleitung mit STEFAN-BOLTZMANN-Randbedingung 524.5 Instationäre Wärmeleitung mit STEFAN-BOLTZMANN-Randbeding

ung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 534.6 Stationäres, thermoelastisches Problem mit STEFAN-BOLTZMANN-

Randbedingung . . . . . . . . . . . . . . . . . . . . . . . . . . . . 544.7 Instationäres, thermoelastisches Problem mit STEFAN-

BOLTZMANN-Randbedingung . . . . . . . . . . . . . . . . . . . . 554.8 Erdmodell mit Sonneneinstrahlung . . . . . . . . . . . . . . . . . 58

5 Das Alter der Erde 60

6 Fazit und Ausblick 63

Inhaltsverzeichnis

Wärmeleitung mit Strahlungsrandbedingung III

Literaturverzeichnis III

Inhaltsverzeichnis

Wärmeleitung mit Strahlungsrandbedingung 1

1 Motivation und Zielsetzung

Wärmetransport wird in Wärmeleitung, Konvektion und Strahlungstransporteingeteilt. Bei der Wärmeleitung wird thermische Energie durch Stöße von Atomen untereinander übertragen, ohne dass diese ihre Position im zeitlichen Mitteländern. Im Gegensatz dazu steht die Konvektion, bei der Materie von einem Ortzu einem Anderen getragen wird. Mit dem Materiestrom ist ein Energietransport verbunden, welcher als konvektiver Wärmetransport bezeichnet wird. BeimWärmetransport durch Strahlung führen elektrodynamische Wechselwirkungenauf atomarer Ebene zur Abgabe elektromagnetischer Strahlung. Diese durchthermische Bewegung angeregte Strahlung wird thermische- oder auch Wärmestrahlung genannt. Im thermischen Gleichgewicht wird sie durch die PlanckscheStrahldichte beschrieben.

Von den drei Mechanismen ist nur der Strahlungstransport in der Lage, Wärmedurch Vakuum zu transportieren. Entsprechend kühlen Körper im Weltraumausschließlich über die Abgabe von Wärmestrahlung ab.

Ziel dieses Projektes ist es, die thermische Entwicklung der idealisierten, elastischen Erde zu modellieren und mithilfe der Finiten-Elemente-Methode zuuntersuchen. Auf der Oberfläche wird Wärme über thermische Strahlung abgegeben und empfangen, wie in Abbildung 1.1 skizziert, während innerhalb desPlaneten Energietransport durch Wärmeleitung nach dem Gesetz von Fouriervorliegt. Zusätzlich wird die durch radioaktive Zerfälle zugeführte Energie durcheine Wärmequelle berücksichtigt.

Abb. 1.1: Ein schematisches Segment der Erde, welches einerseits durch Sonneneinstrahlung Energie empfängt, andererseits durch thermische Strahlung Wärme wiederabgibt.

2

2 Modellbildung

Dieses Kapitel befasst sich mit der thermischen Modellbildung der Erde inKugelkoordinaten, wie in Abbildung 2.1 dargestellt ist. Als Ausgangspunktwerden die globale Massen-, Impuls- und Energiebilanz gewählt, welche auflokale Form gebracht werden. Anschließend werden konstitutive Gleichungenfür die Materialfunktionen und Flüsse verwendet, um die zu bestimmendenFeldgleichungen für die Temperatur und das Verschiebungsfeld zu erhalten.

Es wird kein Elektromagnetismus betrachtet. Das erscheint zunächst nicht zielfüh

𝑥

𝑧

𝑦

𝑅

Ω

𝑥

𝜙

𝜗

Abb. 2.1: Die Erdkugel mit kartesischen Koordinatenkreuz. Die Analyse der Feldgleichungen wird in Kugelkoordinaten durchgeführt, dabei ist 𝜗 der Polarwinkel und 𝜙 derAzimuthalwinkel.

Kapitel 2. Modellbildung


rend, da Wärmestrahlung ein elektromagnetisches Phänomen ist. Allerdings kanndiese nicht vollständig mit der klassischen Physik beschrieben werden. Erst durchdie Quantisierung der Energie, die das elektromagnetische Feld transportierenkann, erhält man eine hinreichend genaue Beschreibung der Wärmestrahlung:Das Plancksche Strahlungsgesetz.

Die Oberfläche der Erde besitzt eine Temperatur im Bereich von 300 K. Ausdem Gesetz der Planckschen Strahldichte geht hervor, dass der Großteil derWärmeenergie bei Wellenlängen kleiner als 100 µm abgegeben wird. Diese Skalaist gegenüber dem Durchmesser der Erde mit rund 12000 km sehr klein. Entsprechend wird die Wellennatur der Strahlung vernachlässigt. Stattdessen nutzenwir Methoden der geometrischen Optik, um die Wärmestrahlung in Randbedingungen einzuarbeiten.

2.1 Massen-, Impuls- und Energiebilanz

Als Quelle für dieses Unterkapitel wurde [Dreyer 2013] verwendet. Um dieFeldgleichungen für die Temperatur und die Verschiebung zu erhalten, wertenwir die globale Bilanz der Masse, des Impulses und der Energie aus. Diese werdenauf einem abgeschlossenen, einfach zusammenhängenden Gebiet Ω betrachtet undauf die jeweilige lokale Form gebracht. Anschließend wird aus der Energiebilanzdie Bilanz der inneren Energie abgespalten und mit Hilfe konstitutiver Gesetzedie Feldgleichungen der linearen Thermoelastizität hergeleitet.

Für die Arbeit mit Bilanzgleichungen sind das Reynoldsche Transportheorem,der Integralsatz von Gauss-Ostrogradski und das Lokalisierbarkeitstheoremhilfreich. Das Transportheorem lautet

dd𝑡

ˆ

Ω(𝑡)

𝛷(𝑥,𝑡) d𝑉 =ˆ

Ω(𝑡)

𝜕𝛷

𝜕𝑡(𝑥,𝑡) d𝑉 +

˛

𝜕Ω(𝑡)

𝛷(𝑥, 𝑡)𝑤(𝑥, 𝑡) · 𝑛(𝑥, 𝑡) d𝐴

(2.1)

wobei 𝑤 das Geschwindigkeitsfeld der Randelemente von Ω ist. Das Volumenkann dabei beliebiger Natur sein und das Geschwindigkeitsfeld 𝑤 muss demnachnicht mit dem Geschwindigkeitsfeld 𝑣, dem die eingeschlossene Materie unterliegt,übereinstimmen. Gilt 𝑤 = 𝑣, so spricht man von einem materiellen Volumen. Wirwollen im Folgenden ausschließlich solche Gebiete Ω untersuchen, die materiellsind. Der Integralsatz von Gauss-Ostrogradski lautet

˛

𝜕Ω(𝑡)

𝛷(𝑥, 𝑡) ⊕ 𝑛 d𝐴 =ˆ

Ω(𝑡)

𝛷(𝑥,𝑡) ⊕ ∇ d𝑉 (2.2)

wobei „⊕” ein beliebiges, lineares Produkt kennzeichnet.

Massen-, Impuls- und Energiebilanz

4

Wir können Gleichung (2.2) auf das Reynoldsche Transportheorem anwendenund erhalten

dd𝑡

ˆ

Ω(𝑡)

𝛷(𝑥,𝑡) d𝑉 =ˆ

Ω(𝑡)

(𝜕𝛷

𝜕𝑡(𝑥,𝑡) + (𝛷(𝑥, 𝑡) ⊗ 𝑤(𝑥, 𝑡)) · ∇

)d𝑉 . (2.3)

Zuletzt wird noch das Lokalisierbarkeitstheorem benötigt. Sei 𝐷 ⊂ Ω beliebigund 𝑓 eine Abbildung 𝐷 → R𝑝, dann gilt

ˆ

𝐷

𝑓(𝑥) d𝑉 = 0 , ∀𝐷 ⊂ Ω ⇔ 𝑓 = 0 . (2.4)

Andernfalls ließe sich immer eine Menge 𝐷 finden, auf der das Integral nichtverschwindet.

Im Folgenden sind die Argumente aus Gründen der Übersichtlichkeit unterdrückt.Falls nicht anders angemerkt, sind alle Felder und Parameter eine Funktion vonOrt und Zeit.

Zunächst wird eine globale Bilanz auf einem materiellen Volumen Ω für einphysikalisches Feld 𝛷 mit dem Fluss 𝛹 betrachtet

dd𝑡

ˆ

Ω

𝛷 d𝑣 =˛

𝜕Ω

𝛹 · 𝑛 d𝐴+ˆ

Ω

(𝑍 + 𝛤 ) d𝑉 . (2.5)

Dabei ist 𝑍 die Zufuhr und 𝛤 die Produktion der Größe Φ. Anwenden derGleichungen (2.2) und (2.3) auf die globale Bilanz liefert

ˆ

Ω

(𝜕𝛷

𝜕𝑡+ [𝛷 ⊗ 𝑣 − 𝛹 ] · ∇ − 𝑍 − 𝛤

)d𝑉 = 0 . (2.6)

Durch die Forderung nach einem materiellen Volumen wurde die zeitliche Entwicklung des Gebietes Ω eingeschränkt, nicht aber seine Form. Somit könnenwir Ω auf einem betrachteten Körper beliebig wählen, womit die Voraussetzungfür das Lokalisierbarkeitstheorem erfüllt ist. Dieses angewandt ergibt die lokaleBilanz für Φ

𝜕𝛷

𝜕𝑡+ [𝛷 ⊗ 𝑣 − 𝛹 ] · ∇ − 𝑍 − 𝛤 = 0 . (2.7)

Die globale Massenbilanz für ein materielles Volumen lautet

dd𝑡

ˆ

Ω

𝜌d𝑉 = 0 . (2.8)

In einem materiellen Volumen gibt es keinen Materiefluss über den Rand, entsprechend verschwindet der Flussterm und da Masse eine Erhaltungsgröße ist,



gibt es auch keine Produktion. Eine Massenzufuhr existiert in der klassischenPhysik nicht. Daher ist nach Gleichung (2.7) die lokale Bilanz gegeben durch

𝜕𝜌

𝜕𝑡+ 𝜌𝑣 · ∇ = 0 . (2.9)

Mit der lokalen Bilanz der Masse kann die materielle Zeitableitung des Feldes𝜌E(𝑥, 𝑡)𝛷(𝑥, 𝑡) umgeschrieben werden. Dazu sei 𝜒(𝑋, 𝑡) die Bewegungsfunktion.Sie bildet die Lagrangeschen Koordinaten 𝑋 auf die räumlichen Koordinaten𝑥 ab. Die Funktion 𝛷E(𝑥, 𝑡) wird 𝛷 in Eulerscher Darstellung genannt, analog dazu wird 𝛷L(𝑋, 𝑡) als 𝛷 in Lagrangescher Darstellung bezeichnet. Diematerielle Zeitableitung ist definiert als die zeitliche Änderung eines Feldes beikonstantem 𝑋:

𝛷∙ := 𝜕𝛷E(𝜒(𝑋, 𝑡), 𝑡)𝜕𝑡

𝑋

= 𝜕𝛷L(𝑋, 𝑡)𝜕𝑡

𝑋

. (2.10)

Um hervorzuheben welche partiellen Ableitungen gemeint sind, steht am senkrechten Strich das jeweils konstante Argument. Für ein Feld in LagrangescherDarstellung ist die materielle Zeitableitung identisch mit der partiellen Zeitableitung. In Eulerscher Darstellung 𝛷E(𝑥, 𝑡) gilt hingegen nach der Kettenregel

𝛷∙E = 𝜕𝛷E

𝜕𝑥

𝑡

· 𝜕𝜒

𝜕𝑡

𝑋

+ 𝜕𝛷E𝜕𝑡

𝑥

𝜕𝑡

𝜕𝑡

𝑋

= 𝜕𝛷E𝜕𝑡

+ (𝛷E ⊗ ∇E) · 𝑣

Zusammen mit Gleichung (2.9) gilt dann

𝜌E𝛷∙E = 𝜕𝜌E𝛷E

𝜕𝑡+ (𝛷E ⊗ 𝜌E𝑣) · ∇E (2.11)

wie man durch Ausdifferenzieren bestätigt

𝜕𝜌E𝛷E𝜕𝑡

+ (𝛷E ⊗ 𝜌E𝑣) · ∇E = 𝜕𝜌E𝜕𝑡

𝛷E + 𝜌E𝜕𝛷E𝜕𝑡

+ 𝛷E (𝜌E𝑣 · ∇E) + (𝛷E ⊗ ∇E) · 𝜌E𝑣

= 𝜌E

(𝜕𝛷E𝜕𝑡

+ 𝛷E ⊗ ∇E

)⏟ ⏞

=𝛷∙E

+𝛷E

(𝜕𝜌E𝜕𝑡

+ 𝜌E𝑣 · ∇E

)⏟ ⏞

=0

= 𝜌E𝛷∙E .

Mit Gleichung (2.11) kann die Form der lokalen Bilanz (2.7) weiter umgeformtwerden zu

𝛷∙ − 𝛹 · ∇ − 𝑍 − 𝛤 = 0 . (2.12)

Mit der materiellen Zeitableitung lautet die lokale Bilanz der Masse

𝜌∙ + 𝜌 (𝑣 · ∇) = 0 . (2.13)


6

Hier ist zu beachten, dass ∇ die partiellen Ableitungen nach den räumlichenKoordinaten 𝑥 bezeichnet, 𝑣 hingegen eine Funktion der materiellen Koordinaten𝑋 ist. Um die Ableitung auszuwerten, muss die materielle Geschwindigkeit indie räumliche Darstellung überführt werden.

Die globale Impulsbilanz lautet

dd𝑡

ˆ

Ω

𝜌𝑣 d𝑉 =ˆ

Ω

𝜌𝑏 d𝑉 +˛

𝜕Ω

𝑡 d𝐴 . (2.14)

Das Produkt 𝜌𝑣 wird Impulsdichte genannt. Sie kann nur durch äußere Krafteinwirkung geändert werden. Eingeprägte Kräfte werden in Oberflächenkräfte 𝑡,die man auch Traktion nennt, und in Volumenkräfte 𝜌𝑏 aufgeteilt. Diese Aufteilung geht auf Cauchy zurück, welcher über das sogenannte Tetraederargumentzeigen konnte, dass die Oberflächenkräfte eine lineare Funktion des nach außenorientierten Normalenvektors 𝑛 sind. Da sich jede lineare Abbildung durch einenTensor entsprechender Stufe darstellen lässt, gibt es demnach einen speziellenTensor, der die Abbildung von 𝑛 nach 𝑡 vermittelt. Daher lässt sich die Traktiondarstellen als

𝑡 = 𝑇 · 𝑛 , (2.15)

wobei 𝑇 der Cauchy-Spannungstensor ist. Dies in die globale Bilanz des Impulseseingesetzt ergibt

dd𝑡

ˆ

Ω

𝜌𝑣 d𝑉 =ˆ

Ω

𝜌𝑏 d𝑉 +˛

𝜕Ω

𝑇 · 𝑛 d𝐴 . (2.16)

Man identifiziert 𝑇 als den Impulsfluss und 𝜌𝑏 als die Impulszufuhr. Eine Impulsproduktion existiert nicht, da sie eine Erhaltungsgröße ist. Aus Gleichung (2.7)erhält man die lokale Form der Impulsbilanz,

𝜕𝜌𝑣

𝜕𝑡+ (𝜌𝑣 ⊗ 𝑣 − 𝑇 ) · ∇ − 𝜌𝑏 = 0 , (2.17)

beziehungsweise aus Gleichung (2.12) mit der materiellen Zeitableitung

𝜌𝑣∙ − 𝑇 · ∇ − 𝜌𝑏 = 0 . (2.18)

Die globale Bilanz der Energie für ein materielles Volumen Ω ist

dd𝑡

ˆ

Ω

𝜌

(𝑢+ 𝑣 · 𝑣

2

)d𝑉 =

ˆ

Ω

(𝑟 + 𝜌𝑏 · 𝑣) d𝑉 +˛

𝜕Ω

(𝑡 · 𝑣 − 𝑞 · 𝑛) d𝐴 . (2.19)

Dabei wird 𝑢 als spezifische innere Energie und 𝑣 · 𝑣/2 als spezifische kinetische



Energie bezeichnet. Die Größe 𝑟 ist eine Wärmequelle und 𝑞 ist die Wärmestromdichte. Aus historischen Gründen werden Komponenten von 𝑞 parallel zurFlächennormale 𝑛 positiv gezählt. Da die Flächennormale nach außen gerichtetist, führt ein positiver Wärmestrom 𝑞 zur Abnahme des Energiegehalts desSystems. Daher ist die Änderung der Energie gleich der negativen Projektionvon 𝑞 auf 𝑛.Setzt man erneut den Spannungstensor ein, erhält man

dd𝑡

ˆ

Ω

𝜌

(𝑢+ 𝑣 · 𝑣

2

)d𝑉 =

ˆ

Ω

(𝑟 + 𝜌𝑏 · 𝑣) d𝑉 +˛

𝜕Ω

(𝑣 · 𝑇 − 𝑞) · 𝑛 d𝐴 . (2.20)

Man identifiziert den Energiefluss als 𝑣 · 𝑇 − 𝑞, die Energiezufuhr durch 𝜌𝑏 · 𝑣und die Energieproduktion mit 𝑟. Energie ist eine Erhaltungsgröße, dennochgibt es in diesem Fall eine Energieproduktion. Diese wurde eingeführt, um diedurch radioaktiven Zerfall zugeführte Energie zu berücksichtigen, ohne ein Mehrkomponentensystem betrachten zu müssen. Eine rigorose Behandlung müssteBilanzen für die Zerfallsprodukte und den ursprünglichen Stoff einführen, wasdie Berechnung erheblich komplizierter macht. Die Energieproduktion wäre danneine Energiezufuhr, welche zu einer Abnahme der Masse im System führt.

Nach Gleichung (2.7) ist die lokale Energiebilanz gegeben durch

𝜕

𝜕𝑡

[𝜌

(𝑢+ 𝑣 · 𝑣

2

)]+[𝜌𝑣

(𝑢+ 𝑣 · 𝑣

2

)+ 𝑞 − 𝑣 · 𝑇

]· ∇ − 𝑟 − 𝜌𝑏 · 𝑣 = 0 (2.21)

und wiederum mit der materiellen Zeitableitung:

𝜌

(𝑢+ 𝑣 · 𝑣

2

)∙+ [𝑞 − 𝑣 · 𝑇 ] · ∇ − 𝑟 − 𝜌𝑏 · 𝑣 = 0 . (2.22)

Es ist möglich, aus der lokalen Energiebilanz die lokale Bilanz der inneren Energieabzuspalten. Dazu wird die Impulsbilanz (2.18) von links mit 𝑣 multipliziert:

𝑣 · 𝜌𝑣∙ − 𝑣 · (𝑇 · ∇) − 𝑣 · 𝜌𝑏 = 0 . (2.23)

Nach der Produktregel gelten

(𝑣 · 𝑣)∙ = 2𝑣 · 𝑣∙ und 𝑣 · (𝑇 · ∇) = (𝑣 · 𝑇 ) · ∇ − (𝑣 ⊗ ∇) ··𝑇 . (2.24)

Eingesetzt ergibt sich die lokale Bilanz der kinetischen Energie:

𝜌

(𝑣 · 𝑣

2

)∙− (𝑣 · 𝑇 ) · ∇ + (𝑣 ⊗ ∇) ··𝑇 − 𝜌𝑣 · 𝑏 = 0 . (2.25)

Zieht man die Bilanz der kinetischen Energie von der Bilanz der Gesamtenergie


8

ab, erhält man die Bilanz der inneren Energie:

𝜌𝑢∙ + 𝑞 · ∇ − 𝑟 − (𝑣 ⊗ ∇) ··𝑇 = 0 . (2.26)

Den Term (𝑣 ⊗ ∇) ··𝑇 nennt man innere Reibung, er wird häufig vernachlässigt.

Nachfolgend sind die grundlegenden Bilanzen noch einmal zusammengefasst:

(Massenbilanz) 𝜌∙ + 𝜌 (𝑣 · ∇) = 0 , (2.27)(Impulsbilanz) 𝜌𝑣∙ − 𝑇 · ∇ − 𝜌𝑏 = 0 , (2.28)

(innere Energiebilanz) 𝜌𝑢∙ + 𝑞 · ∇ − 𝑟 − (𝑣 ⊗ ∇) ··𝑇 = 0 . (2.29)

2.2 Konstitutive Gleichungen

Als Quellen für dieses Unterkapitel wurden [Bertram 2013] und [Truesdell 2013]verwendet. Für den weiteren Bericht werden kleine Deformationen und Drehungen angenommen, sodass Bezugs- und Referenzplatzierung annähernd zusammenfallen und zwischen materiellen und räumlichen Koordinaten nicht weiterunterschieden werden muss. Zudem sei die Dichte als konstant gegeben.

Es werden konstitutive Gleichungen für die Materialfunktionen benötigt, welchemit dem zweiten Hauptsatz der Thermodynamik verträglich sein müssen. Dieserbesagt, dass die Entropieproduktion nie negativ wird. Zu diesem Zweck wird dieglobale Bilanz der Entropie eingeführt

dd𝑡

ˆ

Ω

𝜌𝜂 d𝑉 =ˆ

Ω

(𝑟

𝑇+ 𝜎

)d𝑉 −

ˆ

𝜕Ω

𝑞

𝑇· 𝑛 d𝐴 , (2.30)

wobei 𝜂 die spezifische Entropiedichte, 𝑟/𝑇 die Entropiezufuhr, 𝑞/𝑇 der Entropiestrom und 𝜎 die Entropieproduktion sind. Dann ist die lokale Entropiebilanzgegeben durch

𝜌𝜂∙ +(

𝑞

𝑇

)· ∇ = 𝑟

𝑇+ 𝜎 . (2.31)

Der zweite Hauptsatz wird im Rahmen der rationalen Thermodynamik ausgewertet, weshalb ein geeigneter Zustandsraum gewählt werden muss. Die gesuchtenFelder sind die Dehnungen 𝐸 und die Temperatur 𝑇 . Zudem ist gut belegt, dassder Wärmestrom in Materie verschwindet, falls der Temperaturgradient verschwindet. Daher wird zusätzlich der Temperaturgradient in den Zustandsraumaufgenommen,

Z := {𝑇,𝐸, 𝑇∇} . (2.32)



Der linearisierte Dehnungstensor ist mit den Verschiebungen 𝑢 gegeben durch

𝐸 = 12 (𝑢 ⊗ ∇ + ∇ ⊗ 𝑢) =: def (𝑢) . (2.33)

Des Weiteren gilt

𝐸∙ = def (𝑣) . (2.34)

Alle vom Zustandsraum abhängigen Größen werden Materialfunktionen genannt.In diesem Fall umfasst dies die innere Energie 𝑢, den Wärmestrom 𝑞, die Entropiedichte 𝜂, die freie Energie 𝜓 und den Spannungstensor 𝑇 . Auf die freie Energiewird später weiter eingegangen. Im Rahmen der Thermoelastodynamik nenntman den Zustandsraum auch thermo-kinematischen Zustand, während die Mengeder Materialfunktionen kaloro-dynamischer Zustand genannt wird.

Es werden keine polaren Medien betrachtet, somit gibt es keinen Spin und keineverteilten Momente. In diesem Fall ist die Drehimpulsbilanz genau dann erfüllt,falls der Spannungstensor symmetrisch ist: 𝑇 = 𝑇 T. In das Produkt (𝑣 ⊗ ∇) ··𝑇geht damit nur der symmetrische Anteil von 𝑣 ⊗ ∇ ein und man kann die innereReibung umschreiben zu

(𝑣 ⊗ ∇) ··𝑇 = sym (𝑣 ⊗ ∇) ··𝑇 = 𝐸∙ ··𝑇 . (2.35)

Der zweite Hauptsatz besagt, dass die Entropieproduktion 𝜎 niemals negativ ist.Anders ausgedrückt:

𝜎 ≥ 0 ⇒ 𝜌𝜂∙ +(

𝑞

𝑇

)· ∇ − 𝑟

𝑇≥ 0 , (2.36)

wobei die lokale Bilanz der Entropie nach der Entropieproduktion umgestelltwurde. Führt man die Divergenz in Gleichung (2.36) nach der Produktregel aus,lässt sich die Entropieproduktion weiter umformen zu

𝜌𝜂∙ + 𝑞 · ∇𝑇

− 𝑞 · (𝑇∇)𝑇 2 − 𝑟

𝑇≥ 0 . (2.37)

Die Entropiedichte ist eine Funktion des Zustandsraumes. Führt man die materielle Zeitableitung nach der Kettenregel aus, findet man

𝜌

(𝜕𝜂

𝜕𝑇𝑇 ∙ + 𝜕𝜂

𝜕𝐸··𝐸∙ + 𝜕𝜂

𝜕𝑇∇· (𝑇∇)∙

)+ 𝑞 · ∇

𝑇− 𝑞 · (𝑇∇)

𝑇 2 − 𝑟

𝑇≥ 0 . (2.38)

Für 𝑞 · ∇ wird die lokale Bilanz der inneren Energie eingesetzt, nachdem inihr ebenfalls die materielle Zeitableitung der inneren Energiedichte nach der

Konstitutive Gleichungen

10

Kettenregel ausgewertet wurde:

𝑞 · ∇ = 𝑟 + 𝐸∙ ··𝑇 − 𝜌

(𝜕𝑢

𝜕𝑇𝑇 ∙ + 𝜕𝑢

𝜕𝐸··𝐸∙ + 𝜕𝑢

𝜕𝑇∇· (𝑇∇)∙

). (2.39)

Dann führt Einsetzen in die lokale Entropiebilanz auf

𝜌

(𝜕𝜂

𝜕𝑇𝑇 ∙ + 𝜕𝜂

𝜕𝐸··𝐸∙ + 𝜕𝜂

𝜕𝑇∇· (𝑇∇)∙

)− 𝑞 · (𝑇∇)

𝑇 2 − 𝑟

𝑇+

+ 1𝑇

(𝑟 + 𝐸∙ ··𝑇 − 𝜌

(𝜕𝑢

𝜕𝑇𝑇 ∙ + 𝜕𝑢

𝜕𝐸··𝐸∙ + 𝜕𝑢

𝜕𝑇∇· (𝑇∇)∙

))≥ 0 . (2.40)

In der Entropieproduktion stehen Ableitungen der Zustandsraumvariablen, zumBeispiel 𝑇 ∙, die selbst nicht im Zustandsraum liegen. Diese nennt man höhereAbleitungen. Sortieren von Gleichung (2.40) nach Zustandsraumvariablen undhöheren Ableitungen ergibt:

𝜌

(𝜕𝜂

𝜕𝑇− 1𝑇

𝜕𝑢

𝜕𝑇

)𝑇 ∙ + 𝜌

(𝜕𝜂

𝜕𝐸− 1𝑇

𝜕𝑢

𝜕𝐸+ 𝑇

𝜌𝑇

)··𝐸∙+

+ 𝜌

(𝜕𝜂

𝜕𝑇∇− 1𝑇

𝜕𝑢

𝜕𝑇∇

)· (𝑇∇)∙ − 𝑞 · (𝑇∇)

𝑇 2 ≥ 0 . (2.41)

Mit der freien Energie 𝜓 := 𝑢−𝑇𝜂 lässt sich diese Gleichung kompakter schreiben

− 𝜌

𝑇

(𝜕𝛹

𝜕𝑇+ 𝜂

)𝑇 ∙ − 𝜌

𝑇

(𝜕𝜓

𝜕𝐸− 𝑇

𝜌

)··𝐸∙ − 𝜌

𝑇

(𝜕𝜓

𝜕𝑇∇

)· (𝑇∇)∙ −

− 𝑞 · (𝑇∇)𝑇 2 ≥ 0 . (2.42)

Vor den höheren Ableitungen stehen Faktoren, die aus den Materialfunktionenund den Zustandsraumvariablen aufgebaut sind. Diese Faktoren können dahernicht in Abhängigkeit von den höheren Ableitungen gewählt werden. Es wirdpostuliert, dass jeder thermodynamische Prozess erlaubt sei und die Materialfunktionen die Verträglichkeit mit dem zweiten Hauptsatz unter allen Umständengarantieren. Durch die Beliebigkeit des Prozesses sind die höheren Ableitungenwählbar, womit Gleichung (2.42) nur dann für alle höheren Ableitungen erfülltsein kann, falls die Faktoren vor den höheren Ableitungen separat verschwinden.



Man erhält folgenden Satz von Gleichungen:

(𝑇 ∙) : 𝜕𝜓

𝜕𝑇+ 𝜂 = 0 , (2.43)

(𝐸∙) : 𝜕𝜓

𝜕𝐸− 𝑇

𝜌= 0 , (2.44)

((𝑇∇)∙) : 𝜕𝜓

𝜕𝑇∇= 0 , (2.45)

(Entropieproduktion) : −𝑞 · (𝑇∇)𝑇 2 ≥ 0 . (2.46)

Im Zusammenhang mit dem Verfahren nach Liu werden diese auchLiu-Gleichungen genannt.

Aus Gleichung (2.45) folgt, dass die freie Energie keine Funktion des Temperaturgradienten ist. Weiter erhält man aus Gleichung (2.43), dass die Entropiedichteebenfalls keine Funktion des Temperaturgradienten ist. Stellt man die freieEnergie nach der inneren Energie um, erhält man

𝑢 = 𝛹 + 𝑇𝜂 . (2.47)

Die rechte Seite hängt nicht vom Temperaturgradienten ab und damit ist auchdie innere Energie vom Temperaturgradienten unabhängig.

Aus Gleichung (2.44) folgt, dass die freie Energie ein Potential für die Spannungenist

𝜌𝜕𝜓

𝜕𝐸= 𝑇 . (2.48)

Damit ist auch der Spannungstensor keine Funktion des Temperaturgradienten. Somit sind, den Wärmestrom 𝑞 ausgenommen, alle Materialfunktionenunabhängig vom Temperaturgradienten.

Die Entropieproduktion aus Gleichung (2.46) schränkt die Abhängigkeit des Wärmestroms vom Zustandsraum nicht ein. Sie besagt lediglich, dass Wärmestromund Temperaturgradient einen stumpfen Winkel einschließen. Der einfachsteAnsatz, der den zweiten Hauptsatz erfüllt, ist ein konstitutives Gesetz für denWärmestrom, dass linear im Temperaturgradienten 𝑇∇ ist. Dann ist die Entropieproduktion quadratisch in 𝑇∇ ist und damit stets positiv. Dies führt auf dasGesetz der Fourierschen Wärmeleitung:

𝑞 = −𝜅 · 𝑇∇ . (2.49)

Damit ist die vorläufige Auswertung der Gleichungen (2.43) bis (2.46) abgeschlossen. Als nächstes wird der Spannungstensor betrachtet. Dieser wird in eineTaylorreihe um die Referenzplatzierung entwickelt, welche als spannungs- und


12

dehnungsfrei, sowie mit der Temperatur 𝑇R angenommen wird. Im Rahmen derlinearen Elastizität wird die Entwicklung nach dem ersten Glied abgebrochen:

𝑇 (𝑇,𝐸) = 𝜕𝑇

𝜕𝐸(𝑇R,0) ··𝐸 + 𝜕𝑇

𝜕𝑇(𝑇R,0) (𝑇 − 𝑇R) + 𝒪

(𝐸2, 𝑇 2

). (2.50)

Die Proportionalität zu den Dehnungen ist bereits aus der linearen Elastizitätstheorie bekannt, dort wird der Proportionalitätsfaktor Steifigkeitstetradegenannt, weshalb

𝜕𝑇

𝜕𝐸(𝑇R,0) = 𝜌

𝜕2𝜓

𝜕𝐸2 (𝑇R,0) =: 𝐶 (𝑇R) (2.51)

folgt. Die Proportionalität zur Temperatur ist hingegen weniger bekannt. Mandefiniert

𝜕𝑇

𝜕𝑇(𝑇R,0) = 𝜌

𝜕2𝛹

𝜕𝑇𝜕𝐸(𝑇R,0) = −𝜌 𝜕𝜂

𝜕𝐸(𝑇R,0) := 𝑀(𝑇R) , (2.52)

welchen man „Spannungs-Temperatur-Tensor“nennt. Damit lautet das lineareMaterialgesetz für die Spannungen:

𝑇 (𝑇,𝐸) = 𝐶(𝑇R) ··𝐸 + 𝑀(𝑇R) (𝑇 − 𝑇R) . (2.53)

Bekannter ist ein anderes Maß, dass Spannung und Temperatur koppelt: Der Wärmeausdehnungskoeffizient 𝛼. Er gibt die Proportionalität zwischen einachsigenDehnungen und einer Temperaturdifferenz an. Im dreidimensionalen Fall kanndiese Größe durch einen Tensor verallgemeinert werden, welcher im FolgendenWärmeausdehnungstensor genannt wird. Als nächstes wird der Zusammenhangzwischen Wärmeausdehnungstensor und dem Spannungs-Temperatur-Tensorbetrachtet, wobei anzumerken ist, dass für eine rigorose Untersuchung ein anderer Zustandsraum gewählt werden muss, in dem 𝑇 anstelle von 𝐸 steht. AusGründen der Übersicht wird hier nur eine grobe Abhandlung gegeben. Es wirdeine Abhängigkeit der Gesamtdehnungen 𝐸tot vom Spannungstensor 𝑇 undder Temperatur 𝑇 angenommen. Außerdem wird eine additive Zerlegung derGesamtdehnung 𝐸 in thermische 𝐸th und mechanische 𝐸m Anteile postuliert:

𝐸(𝑇 , 𝑇 ) = 𝐸m(𝑇 ) + 𝐸th(𝑇 ) . (2.54)

Für die Thermodehnungen wird eine Taylorreihe aufgestellt, wobei die Bezugsplatzierung dehnungsfrei sei:

𝐸th(𝑇 ) = 𝜕𝐸th

𝜕𝑇(𝑇R) (𝑇 − 𝑇R) + 𝒪(𝑇 2) . (2.55)

In der linearen Näherung wird die Reihe nach dem ersten Glied abgebrochen



und die Größe

𝜕𝐸th

𝜕𝑇(𝑇R) := 𝛼(𝑇R) (2.56)

als Wärmeausdehnungstensor definiert. Als nächstes wird angenommen, dassdie Spannungen nur von den mechanischen Dehnungen abhängen. Dann gilt imRahmen der linearen Elastizitätstheorie

𝑇 = 𝐶 ··𝐸m = 𝐶 ··𝐸 − 𝐶 ··𝐸th = 𝐶 ··𝐸 − 𝐶 ··𝛼(𝑇R) (𝑇 − 𝑇R) .

Ein Vergleich mit Gleichung (2.53) ergibt

𝐶(𝑇R) ··𝛼(𝑇R) = −𝑀(𝑇R) , (2.57)

den Zusammenhang zwischen dem Wärmeausdehnungtensor und dem Temperatur-Spannungs-Tensor.

Abschließend sei angemerkt, dass der Spannungstensor im vorherigen Abschnittals implizit abhängig von 𝑇 dargestellt wurde. Eine rigorose Rechnung liefertallerdings eine explizite Abhängigkeit. Alternativ kann das Gesetz für den Spannungstensor umgestellt werden, um die Dehnungen in Abhängigkeit der Temperatur und der Spannungen auszudrücken. Sei 𝑆 die inverse Steifigkeitstetrade,welche auch Nachgiebigkeitstetrade genannt wird, so erhält man

𝐸 = 𝑆 ··𝑇 − 𝑆 ··𝑀 (𝑇 − 𝑇R) . (2.58)

Dabei wird

𝛼(𝑇R) = −𝑆(𝑇R) ··𝑀(𝑇R) (2.59)

wiederum als Wärmeausdehnungstensor bezeichnet. Das Ergebnis ist identischzu Gleichung (2.57).

Für den Spannungstensor und den Wärmefluss wurden bereits konstitutiveGesetze gefunden. Es bleibt noch die innere Energie, wobei hier nur das Gesetzfür die zeitliche Ableitung angegeben wird. Die innere Energie ist nicht vomTemperaturgradienten abhängig. Daher gilt

𝑢∙ = 𝜕𝑢

𝜕𝑇𝑇 ∙ + 𝜕𝑢

𝜕𝐸··𝐸∙ . (2.60)

Der Faktor 𝜕𝑢/𝜕𝑇 beschreibt die Wärmekapazität 𝑐E bei konstanten Dehnungen.Für die Ableitung nach den Dehnungen wird eine weitere Identität benötigt. DerSpannungs-Temperatur-Tensor ließ sich darstellen als

𝜕𝑇

𝜕𝑇= 𝑀 . (2.61)


14

Zugleich stellt die freie Energie ein Potential für die Spannungen dar. Mit denGleichungen (2.48) und (2.52) folgt die Beziehung

𝜕𝑇

𝜕𝑇= 𝜌

𝜕2𝛹

𝜕𝐸𝜕𝑇= −𝜌 𝜕𝜂

𝜕𝐸= 𝑀 (2.62)

aus der sich

𝜕𝑢

𝜕𝐸= 𝑇

𝜌+ 𝑇

𝜕𝜂

𝜕𝐸= 𝑇

𝜌− 𝑇

𝜌𝑀 (2.63)

ergibt. Eingesetzt in Gleichung (2.60) erhält man

𝑢∙ = 𝑐E𝑇∙ + 1

𝜌𝑇 ··𝐸∙ − 𝑇

𝜌𝑀 ··𝐸∙ , (2.64)

was wiederum, zusammen mit dem Fourierschen Gesetz, in die lokale Bilanzder inneren Energie,

𝜌𝑢∙ + 𝑞 · ∇ = 𝑟 + 𝐸∙ ··𝑇 , (2.65)

eingesetzt wird, um bei der lokalen Bilanz der inneren Energie der Thermoelastodynamik anzukommen:

𝜌𝑐E𝑇∙ − 𝑇𝑀 ··𝐸∙ − (𝜅 · (𝑇∇)) · ∇ = 𝑟 , (2.66)

oder alternativ mit dem Wärmeausdehnungstensor:

𝜌𝑐E𝑇∙ + 𝑇𝛼 ··𝐶 ··𝐸∙ − (𝜅 · (𝑇∇)) · ∇ = 𝑟 . (2.67)

Dabei wurde die Hauptsymmetrie der Steifigkeitstetrade ausgenutzt. Wird dasSpannungsgesetz aus Gleichung (2.53) in die lokale Impulsbilanz eingesetzt, findetman

𝜌𝑣∙ − [𝐶 ·· (𝐸 − 𝛼 (𝑇 − 𝑇R))] · ∇ − 𝜌𝑏 = 0 . (2.68)

In diesem Projekt werden nur isotrope Materialien betrachtet. Dann gelten

𝛼 ··𝐶 = 𝜆 Sp(𝛼1)1 + 2𝜇𝛼1 = 𝛼(3𝜆+ 2𝜇)1 und 𝑞 = −𝜅𝑇∇ . (2.69)

Dabei sind 𝜇 und 𝜆 die Lamé-Konstanten, 𝛼 der Wärmeausdehnungskoeffizientund 𝜅 die Wärmeleitfähigkeit. Außerdem betrachten wir keine eingeprägtenVolumenkräfte 𝑏, da thermische Deformationen im Vordergrund stehen. Wiesich später zeigen wird, sind die thermischen Deformationen sehr groß, sodasseventuelle Deformationen aufgrund von Eigengravitation gegen die thermischenDeformationen vernachlässigt werden können. Eingesetzt in die Impulsbilanz



und die Bilanz der inneren Energie ergeben sich

𝜌𝑐E𝑇∙ − 𝜅Δ𝑇 + 𝛼(3𝜆+ 2𝜇) (𝑇 − 𝑇R) (𝑢∙ · ∇) = 𝑟 , (2.70)

𝜇Δ𝑢 + (𝜆+ 𝜇) (𝑢 · ∇ ⊗ ∇) − 𝛼 (3𝜆+ 2𝜇) (𝑇 − 𝑇R) ∇ = 𝜌𝑢∙∙ . (2.71)

2.3 Rand- und Anfangsbedingungen

Für eine eindeutige Lösung des thermoelastischen Problems werden Rand- undAnfangsbedingungen für das Verschiebungsfeld und das Temperaturfeld benötigt.

Für die Anfangstemperatur und Verschiebung ist eine Funktion vom Ort 𝑥denkbar:

𝑇 (𝑥, 0) = 𝑇A(𝑥) , 𝑥 ∈ Ω , (2.72)𝑢(𝑥, 0) = 𝑢A(𝑥) , 𝑥 ∈ Ω . (2.73)

Da sich die Erde im Weltraum bewegt und die Materie jenseits der Atmosphäresehr dünn verteilt ist, wird der Rand 𝜕Ω als spannungsfrei angenommen:

𝑇 · 𝑛 = 0 , 𝑥 ∈ 𝜕Ω . (2.74)

Für das Temperaturfeld werden drei verschiedene Randbedingungen betrachtet: Eine Robin-Randbedingung, eine Stefan-Boltzmann-Randbedingungund eine Stefan-Boltzmann-Randbedingung mit Sonneneinstrahlung. DieRobin-Randbedingung hat im Kontext der Erde physikalisch keine Bedeutung, sie dient als Vergleichsproblem. Die Strahlungsrandbedingungen werden mit dem Planckschen Strahldichte Gesetz modelliert, was auf das Stefan-Boltzmann-Gesetz führt. Dazu werden Methoden der geometrischen Optikverwendet, auf die hier nicht weiter eingangen wird. Die Ergebnisse sind allesamt aus [McCluney 2014] und [Howell u. a. 2010] entnommen. Eine tiefereBeschreibung des Planckschen Strahldichte Gesetzes findet sich zum Beispielin [Fitzpatrick 2006].

Die Robin-Randbedingung lautet:

𝑞 · 𝑛 = ℎ (𝑇 − 𝑇ext) , 𝑥 ∈ 𝜕Ω . (2.75)

Dabei ist ℎ der Wärmeübergangskoeffizient und 𝑇ext die Temperatur eines angrenzenden Mediums. Dieses Gesetz erhält man aus dem Newtonschen Abkühlungsgesetz, welches besagt, dass die abgegebene Wärme eines Körpers proportionalzur Temperaturdifferenz des Körpers zur Umgebung ist. Zumeist wird diesesGesetz eingesetzt, um konvektive Wärmeübergange phänomenologisch zu beschreiben. Dabei wird eine Oberfläche von einem Fluid überströmt. Nahe derOberfläche existiert eine Grenzschicht, in der ein steiler Temperaturgradient

Rand- und Anfangsbedingungen

16

𝑛

𝜀

𝑇ext

R

𝜗

Abb. 2.2: Ein Punkt auf einem Ausschnitt der Erdoberfläche unter einer ideal strahlenden Himmelssphäre mit dem Radius 𝑅 und der Temperatur 𝑇ext. Die Erdoberflächestrahlt und empfängt Strahlungsenergie, auch im Fall eines thermischen Gleichgewichts.

vorliegt. Befindet man sich hinreichend weit von der Oberfläche entfernt, ändertsich die Kerntemperatur des Fluides praktisch nicht mehr. Diese Grenzschichtkann sehr klein sein, sodass ein einfacher Modellansatz einen Temperatursprungan der Grenzfläche zwischen Körper und Fluid modelliert, was auf die Robin-Randbedingung führt. Für die Erde ist das offensichtlich kein sinnvollerAnsatz, da sie in kein Medium eingebettet ist, dass Wärme in relevanten Mengenleiten oder konvektiv transportieren kann. Die Robin-Randbedingung dientlediglich als Einstieg in die Wärmeleitung und als Vergleichsfall zur Strahlungsrandbedingung.

Die Stefan-Boltzmann-Randbedingung lautet

𝑞 · 𝑛 = 𝜎𝜀(𝑇 4 − 𝑇 4

ext

), 𝑥 ∈ 𝜕Ω . (2.76)

Dabei ist 𝑇ext die „Umgebungstemperatur“, 𝜎 = 5,67 × 10−8 W m−2 K−4 dieStefan-Boltzmann-Konstante und 𝜀 der Emissionsgrad der Oberfläche . DieseRandbedingung leitet sich aus dem Planckschen Strahldichte Gesetz ab, welchesgegeben ist durch

𝐿p(𝜆, 𝑇 ) = 2ℎ𝑐20

𝜆51

exp(

ℎ𝑐0𝜆𝑘B𝑇

)− 1

. (2.77)

Dabei ist ℎ = 6,626 × 10−34 J s das Plancksche-Wirkungsquantum, 𝑐0 ≈2,998×108 m s−1 die Lichtgeschwindigkeit im Vakuum, 𝜆 die Vakuumwellenlänge,𝑘B = 1,38 × 10−23 J K−1 die Boltzmann-Konstante und 𝑇 die Temperatur derOberfläche. Betrachtet wird das System in Abbildung 2.2. Ein Punkt Q auf derErdoberfläche liegt unter einer Himmelsphäre mit dem Radius 𝑅. Die pro Flächeabgegebene Strahlungsleistung 𝐽 errechnet sich aus der mit dem Polarwinkelgewichteten, abgegebenen Strahldichte 𝜀𝐿p, integriert über alle Raumwinkel und



alle Wellenlängen,

𝐽(𝑇 ) =π/2ˆ

𝜗=0

2πˆ

𝜙=0

∞

𝜆=0

𝜀𝐿p(𝜆, 𝑇 ) cos(𝜗) sin(𝜗) d𝜆 d𝜙 d𝜗 . (2.78)

Dies ist ein Bose-Einstein-Integral, welches gegeben ist durch

𝐽(𝑇 ) = 𝜀8π5𝑘4

B15 (ℎ𝑐0)3𝑇

4 . (2.79)

Der Faktor 8π5𝑘4B/15(ℎ𝑐0)3 =: 𝜎 ist die Stefan-Boltzmann-Konstante und das

Gesetz wird Stefan-Boltzmann-Gesetz genannt. Es wurde 1879 zunächstexperimentell von Josef Stefan und anschließend 1884 analytisch von LudwigBoltzmann gefunden.

Gleichzeitig erhält die Fläche aus der Umgebung einen Strahlungsfluss 𝐼. Nachdem fotometrischen Grundgesetz sind Strahldichte einer sendenden Fläche undBestrahldichte einer empfangenen Fläche identisch. Daher ist der durch dieHimmelssphäre auf den Punkt 𝑃 treffende Strahlungsfluss gleich:

𝐼(𝑇ext) =π/2ˆ

𝜗=0

2πˆ

𝜙=0

∞

𝜆=0

𝐿p(𝜆, 𝑇ext) cos(𝜗) sin(𝜗) d𝜆 d𝜙 d𝜗 = 𝜎𝑇 4ext . (2.80)

Von dieser Leistung wird ein Teil absorbiert, der Rest reflektiert. Der Quotient ausabsorbierter und einfallender Strahlungsleistung wird Absorptionsgrad genannt.Im thermischen Gleichgewicht sind Absorptionsgrad und Emissionsgrad immeridentisch, weshalb die im Punkt 𝑃 absorbierte Strahlungsleistung gleich 𝜀𝐼(𝑇ext)ist.

Die thermische Strahlung ändert die innere Energie im Punkt Q. Daher istder auf den Normalenvektor 𝑛 im Punkt Q projizierte Wärmestrom gleich derDifferenz zwischen abgestrahlter und absorbierter Strahlungsleistung,

𝑞 · 𝑛 = 𝐽(𝑇 ) − 𝐼(𝑇ext) = 𝜀𝜎(𝑇 4 − 𝑇 4

ext

). (2.81)

Im Gegensatz zur Robin-Randbedingung ist die Interpretation der „Umgebung “hier schwieriger. Zur Umgebung gehören alle Körper, mit denen sich die Oberfläche im Strahlungsaustausch befindet. Für ein Oberflächenelement auf derErde gehören damit unzählige Körper zu dieser Umgebung, wie zum Beispieldie Sonne, der Mond, Planeten, ferne Sterne und Galaxien und die kosmischeHintergrundstrahlung. Der Einfluss dieser Körper variiert stark, weshalb hierzunächst die Einflüsse aller anderen Himmelskörper vernachlässigt werden, sodassnur die kosmische Hintergrundstrahlung verbleibt. Der Einfluss der Sonne wird


18

in der nächsten Randbedingung berücksichtigt. Die kosmische Hintergrundstrahlung besitzt ein Plancksches Spektrum mit einer effektiven Temperatur von𝑇ext = 2,7 K.

Die Stefan-Boltzmann-Randbedingung mit Sonneneinstrahlung ist gegebendurch

𝑞 · 𝑛 = 𝜀𝜎𝑇 4 − 𝜀𝑆 sin(𝜗) cos(𝜙) , (𝑥, 𝜙, 𝜗) ∈ 𝜕Ω × [−π/2, π/2] × [0,π] , (2.82)

wobei 𝑆 = 1354 W m−2 die Solarkonstante ist.

Als Ansatz wird die gesamte abgegebene Energie der Sonne betrachtetund der von der Erde abgegebene Strahlungsfluss wieder durch das Stefan-Boltzmann-Gesetz modelliert. Die Sonne gibt nach [Williams 2004b] inguter Näherung ein Plancksches Spektrum mit einer effektiven Temperaturvon 𝑇S = 5772 K ab. Ihr Radius beträgt 𝑅S = 6,957 × 108 m. Sie hat damit eineOberfläche von 𝐴S = 6,0877 × 1018 m2 und strahlt eine Leistung von

𝑃 = 𝜎𝐴S𝑇4S ≈ 3,831 × 1026 W (2.83)

ab. Die Erde umkreist die Sonne in einer mittleren Entfernung von 𝑅AE =1,5 × 1011 m. Die Leistung 𝑃 verteilt sich dabei auf eine Kugeloberfläche 𝐴E mitdem Radius 𝑅AE. Der Strahlungsfluss durch ein Flächenelement dieser Kugel istdamit

𝑆 = 𝑃𝐴S𝐴E

= 𝜎𝑇 4S

(𝑅S𝑅AE

)2≈ 1353,785 W m−2 . (2.84)

Wir nehmen an, dass dieser Strahlungsfluss konstant über die gesamte Erdoberfläche ist. Damit wird vernachlässigt, dass die Distanz zur Sonne über dieOberfläche der Erde variiert, so ist man am Mittag der Sonne 6 000 km näher alsam Abend oder am Morgen.

Als nächstes wird der Strahlungsfluss 𝑃 auf die Normale der Erdoberflächeprojiziert, um so die absorbierte Strahlungsleistung

𝐼 = 𝜀𝑆 cos(𝜗S) (2.85)

und daraus die Randbedingung

𝑞 · 𝑛 = 𝜀𝜎𝑇 4 − 𝜀𝑆 cos(𝜗S) (2.86)

zu erhalten. Für diesen Zweck wird die Geometrie aus Abbildung 2.3 betrachtet.Es wird ein kartesisches Koordinatensystem in die Erde gelegt, sodass die 𝑥-Achseauf die Sonne und die 𝑧-Achse zum Nordpol zeigt. Der Erdradius sei 𝑅, derVektor 𝑅OS zeigt vom Erdmittelpunkt O zum Sonnenmittelpunkt 𝑆, der Vektor𝑅OQ zeigt vom Erdmittelpunkt auf den Punkt Q und der Vektor 𝑅QS zeigt vom



𝑥

𝑧

𝑦

O𝑅QP

Q𝑛Q

𝑅QS

𝑅OS

𝜗S

S

Abb. 2.3: Die Geometrie für die Projektion des Strahlungsflusses auf die Erdoberfläche.Der Winkel 𝜗S ist gesucht.

Punkt Q auf den Sonnenmittelpunkt S. Dann ist der Vektor 𝑅QS gegeben durch

𝑅QS = 𝑅OS − 𝑅OQ . (2.87)

Mit ihm lässt sich der Kosinus von 𝜗S berechnen

cos(𝜗S) = 𝑅QS · 𝑅OQ‖𝑅QS‖ ‖𝑅OQ‖

. (2.88)

Der Normalenvektor 𝑛Q ist parallel zu 𝑅OQ. Daher gilt

𝑛Q = 𝑅OQ‖𝑅OQ‖

. (2.89)

Weiter gelten

𝑅OS = 𝑅AE𝑒𝑥 , (2.90)𝑅OQ = 𝑅𝑒𝑟 = 𝑅 (sin(𝜗) cos(𝜙)𝑒𝑥 + sin(𝜗) sin(𝜙)𝑒𝑦 + cos(𝜗)𝑒𝑧) . (2.91)

Damit erhält man zunächst den Vektor

𝑅QS = (𝑅AE −𝑅 sin(𝜗) cos(𝜙)) 𝑒𝑥 −𝑅 sin(𝜗) sin(𝜙)𝑒𝑦 −𝑅 cos(𝜗)𝑒𝑧 (2.92)

sowie den Normalenvektor in Q:

𝑛Q = sin(𝜗) cos(𝜙)𝑒𝑥 + sin(𝜗) sin(𝜙)𝑒𝑦 + cos(𝜗)𝑒𝑧 = 𝑒⟨𝑟⟩ . (2.93)

Einsetzen von Gleichungen (2.89), (2.92) und (2.93) in Gleichung (2.88) ergibt


20

nach einigen Umformungen:

cos(𝜗S) =sin(𝜗) cos(𝜙) − 𝑅

𝑅AE√1 +

(𝑅

𝑅AE

)2− 𝑅

𝑅AEsin(𝜙) sin(𝜗)

. (2.94)

Falls die Sonne am Horizont steht, gilt cos(𝜗S) = 0, woraus folgt

sin(𝜗) cos(𝜙) = 𝑅

𝑅AE. (2.95)

Aufgrund der endlichen Entfernung 𝑅AE wird nicht die gesamte Halbkugelbeschienen. Die Randbedingung ist daher nur für Winkel gültig, für die

sin(𝜗) cos(𝜙) ≥ 𝑅

𝑅AE(2.96)

gilt. Für das Sonne-Erde-System gilt 𝑅/𝑅AE = 4 × 10−5. In diesem Projekt wirddaher

lim𝑅/𝑅AE → 0

cos(𝜗S) = sin(𝜗) cos(𝜙) (2.97)

betrachtet. Die Bedingung (2.96) zerfällt damit in 𝜙 ∈ [−π/2, π/2] und 𝜗 ∈ [0,π].Gleichung (2.97) eingesetzt in Gleichung (2.86) ergibt die Randbedingung (2.82).

In diesem Projekt steht die thermische Entwicklung der Erde im Vordergrund.Dabei werden Zeitskalen oberhalb von hundert Millionen Jahren betrachtet. DieErdrotation hat eine Periode von 24 Stunden, weshalb sich auf der betrachteten Zeitskala die Erdrotation schwer auflösen lässt. Entsprechend wird dieSonneneinstrahlung über den Azimuthalwinkel 𝜙 gemittelt:

π/2´𝜙=−π/2

𝜀𝑆 sin(𝜗) cos(𝜙) d𝜙

π

𝜙=−π

d𝜙= 𝜀𝑆

πsin(𝜗) . (2.98)

Die gemittelte Randbedingung lautet damit

𝑞 · 𝑛 = 𝜀𝜎𝑇 4 − 𝜀𝑆

πsin(𝜗) , (𝑥, 𝜗) ∈ 𝜕Ω × [0,π] . (2.99)



2.4 Entdimensionalisierung

Die Feldgleichungen (2.70) und (2.71) werden entdimensionalisiert. Dazu werdendie Größen

𝑢 = 𝑢0�� , 𝑇 = 𝑇1𝑇 + 𝑇0 , 𝑡 = 𝜏𝑡 , 𝑥𝑖 = 𝑅��𝑖 und 𝑟(𝑥) = 𝑟A𝛾(��) (2.100)

mit noch zu bestimmenden Normierungsfaktoren eingeführt. Zusammen mitder Temperaturleitfähigkeit 𝑎 := 𝜅/(𝜌𝑐E), der normierten Referenztemperatur𝑇1𝑇R = 𝑇R − 𝑇0 und der normierten Lamé-Konstante �� = 𝜆/𝜇 ergeben sich

𝑅2

𝜏𝑎

𝜕𝑇

𝜕𝑡− Δ𝑇 + 𝑅𝛼𝑢0𝜇

𝜅𝜏

(3��+ 2

) (𝑇 − 𝑇ref

)(𝜕��

𝜕𝑡· ∇)

= 𝑟A𝑅2

𝜅𝑇1𝛾 , (2.101)

Δ�� +(��+ 1

) (�� · ∇ ⊗ ∇

)− 𝛼𝑇1𝑅

𝑢0

(3��+ 2

) (𝑇 − 𝑇R

)∇ =

= 𝜌𝑅2

𝜇𝜏2𝜕2��

𝜕𝑡2. (2.102)

Mit ∇ ist der Nablaoperator bezüglich der normierten Koordinaten ��𝑖 gemeint.Die charakteristische Zeit 𝜏 wird so gewählt, dass zeitliche und örtliche Änderungdes Temperaturfeldes in der gleichen Größenordnung liegen. Daraus folgt

𝜏 = 𝑅2

𝑎. (2.103)

Des Weiteren werden

𝑢0 = 𝜅𝑅

𝛼𝜇𝑎(3��+ 2

) und 𝑇1 = 𝑟A𝑅2

𝜅(2.104)

gewählt. Zudem wird der Kopplungsfaktor in der Impulsbilanz zusammengefasst:

𝑚 := 𝛼𝑇1𝑅

𝑢0

(3��+ 2

)= 𝛼2𝑅2𝜇𝑎𝑟A

𝜅2

(3��+ 2

)2. (2.105)

Damit sind die Feldgleichungen gegeben durch

𝜕𝑇

𝜕𝑡− Δ𝑇 +

(𝑇 − 𝑇R

)(𝜕��

𝜕𝑡· ∇)

= 𝛤 , (2.106)

Δ�� +(��+ 1

) (�� · ∇ ⊗ ∇

)−𝑚

(𝑇 − 𝑇R

)∇ = 𝜌𝑎2

𝑅2𝜇

𝜕2��

𝜕𝑡2. (2.107)

Entdimensionalisierung

22

Als nächstes werden die Randbedingungen entdimensionalisiert. Mit dem normierten Verzerrungstensor �� gilt für den Verzerrungstensor 𝐸:

𝐸 = sym (𝑢 ⊗ ∇) = 𝑢0𝑅

sym(�� ⊗ ∇

)=: 𝑢0

𝑅�� . (2.108)

Wiederum mit dem normierten Spannungstensor 𝑇 findet man darauf aufbauend

𝑇 = (𝜆 Sp (𝐸) − 𝛼 (3𝜆+ 2𝜇) (𝑇 − 𝑇R)) 1 + 2𝜇𝐸

= 𝜇𝑢0𝑅

((��Sp

(��)

−𝑚(𝑇 − 𝑇R

))1 + 2��

)=: 𝜇𝑢0

𝑅𝑇 . (2.109)

Für die Randbedingung (2.74) folgt damit

𝑇 · 𝑛 = 0 , 𝑥 ∈ 𝜕Ω . (2.110)

Für das Temperaturfeld liegen die drei Randbedingungen (2.75), (2.76) und (2.82)vor. Die Temperaturverschiebung 𝑇0 wird in Abhängigkeit der Randbedingunggesetzt, speziell gilt:

𝑇0 ={𝑇ext , Robin-Randbedingung0 K , Stefan-Boltzmann-Randbedingung

. (2.111)

Für den Wärmestrom findet man

𝑞 = −𝜅𝑇∇ = −𝜅𝑇1𝑅𝑇 ∇ =: 𝜅𝑇1

𝑅�� . (2.112)

Damit folgt für die Robin-Randbedingung

𝜅𝑇1𝑅

�� · 𝑛 = ℎ𝑇1𝑇 ⇒ �� · 𝑛 = ℎ𝑇 . (2.113)

Dabei wurde der normierte Wärmeübergangskoeffizient ℎ = ℎ𝑅/𝜅 definiert.

Für die Stefan-Boltzmann-Randbedingung erhält man

�� · 𝑛 = ��(𝑇 4 − 𝑇 4

ext

), (2.114)

wobei �� = 𝜀𝜎𝑟3A𝑅7/𝜅4 und 𝑇ext = 𝑇ext/𝑇1 gesetzt wurden. Der normierte Strah

lungsparameter �� hängt in hoher Potenz von den physikalischen Eigenschaftendes Systems ab und kann daher schnell sehr groß oder sehr klein werden, wasnumerisch ungünstig sein kann.

Analog ergibt sich für die gemittelte Strahlungsrandbedingung mit Sonneneinstrahlung

�� · 𝑛 = ��𝑇 4 − 𝑆 sin(𝜗) . (2.115)



Der Koeffizient �� ist identisch zur Stefan-Boltzmann-Randbedingung gewählt,die normierte Solarkonstante ist 𝑆 = 𝜀𝑆𝑅/π𝜅𝑇1.

Für die Anfangsbedingungen sind beliebige, hinreichend glatte Funktionen 𝑇Aund 𝑢A vorstellbar. Durch Einsetzen der normierten Temperatur und Verschiebung erhält man die normierte Anfangsbedingung für das Temperatur- undVerschiebungsfeld

𝑇 (��, 0) = 𝑇A(𝑥) − 𝑇0𝑇1

=: 𝑇A(��, 𝑡) , ��(��, 0) = 𝑢A(𝑥)𝑢0

=: ��A(��, 𝑡) . (2.116)

2.5 Parameter des Erdmodells

Für die Feldgleichungen, sowie die Rand- und Anfangsbedingungen werdensinnvolle Parameterwerte benötigt. Da die Erde aus vielen verschiedenen Stoffenbesteht und Daten aus dem Zeitraum der Erdentstehung selten sind, ist dieWahl sinnvoller Material- und Anfangswerte schwierig. Zudem ändern sich Druckund Temperatur von der Erdoberfläche hin zum Kern stark. Der Umstand, dassMaterialeigenschaften im Allgemeinen sowohl eine Funktion des Druckes, alsauch der Temperatur sind, erschwert die Suche nach sinnvollen Parameterwertenweiter. Die elastischen Eigenschaften für typisches Erdgestein wurden [Ji u. a.2010] entnommen, die Wärmekapazität entstammt [Robertson 1988] und dieTemperatur der Erde kurz nach ihrer Entstehung wurde aus [Thomson 1862]genommen. Alle weiteren Parameter sind in [Turcotte u. Schubert 2014] und[Williams 2004a] zu finden.

Die Parameter sind in Tabelle 2.1 aufgeführt, wobei hier nicht auf jeden Wertim Einzelnen eingegangen wird. Die Lamé-Konstanten werden dennoch kurzdiskutiert, da das Auffinden sinnvoller Werte eine besondere Herausforderungdarstellte. Die Erde besteht größtenteils aus Silikatgestein und Eisen, wobeiLetzteres im Erdkern konzentriert ist. Silikatgestein ist ein Sammelbegriff fürsehr viele verschiedene Verbindungen aus Silicium und Sauerstoff, wobei es sichhauptsächlich um Siliciumtrioxidverbindungen handelt. Es unterscheidet sichdemnach vom typischen Sand an Sandstränden, welches eine Siliciumdioxidverbindung ist. Es zeigt sich, dass zumindest in einem Bereich von 20 ∘C bis einigenhundert Grad Celsius, die Lamé-Konstanten in einem Bereich zwischen 10 GPaund 100 GPa liegen, wobei die meisten Gesteine näher an der oberen Grenze zufinden sind. Hier werden

𝜇 = 6,5 × 1010 GPa und 𝜆 = 8 × 1010 GPa (2.117)

gewählt. Diese Werte wurden für Silikatgestein in einer Studie in [Ji u. a. 2010]unter einigen Hundert Kelvin und 1 bar Druck ermittelt. Diese Bedingungen sindweit von den Verhältnissen im Erdinneren entfernt und sind deshalb mit Vorsicht

Parameter des Erdmodells

24

zu behandeln. Für die Wärmekapazität wurde ein Wert von 630 J kg−1K−1 nach[Robertson 1988] gewählt. Für die Starttemperatur wird nach Kelvin einekonstante Temperaturverteilung von 𝑇A = 4 000 K gesetzt.

Tab. 2.1: Die Parameter des Erdmodells.

Bezeichnung Symbol WertErste Lamé-Konstante 𝜆 8 × 1010 GPaZweite Lamé-Konstante 𝜇 6,5 × 1010 GPaWärmekapazität 𝑐𝐸 630 J (kg K)−1

Mittlerer Erdradius 𝑅 6,371 × 106 mErdoberfläche 𝐴 5,1 × 1014 m2

Erdvolumen 𝑉 1,083 × 1021 m3

Erdmasse 𝑀 6 × 1024 kgMittlere Wärmeleitfähigkeit 𝜅 2,2 W (m K)−1

Temperaturgradient 𝑔 3 × 10−2 K m−1

sphärische Albedo 𝐷 0,7Starttemperatur 𝑇A 4000 KUmgebungstemperatur 𝑇ext 2,7 KWärmeausdehnungskoeffizient 𝛼 10−5 m K−1

Solarkonstante 𝑆 1354 W m−2

Dichte 𝜌 = 𝑀/𝑉 5540 kg m−3

Temperaturleitfähigkeit 𝑎 = 𝜅/𝜌𝑐E 6,303 × 10−7 m2 s−1

Emissionsgrad 𝜀 = 1 −𝐷 0,3Wärmeproduktion 𝑄𝜎 = 𝜅𝑔𝐴 3,4 × 1013 WWärmequellenamplitude 𝑟A = 𝑄𝜎/𝑉 3,1 × 10−8 W m−3

Tab. 2.2: Normierte und charakteristische Parameter des Erdmodells.

Bezeichnung Definition WertLamé-Konstante �� = 𝜆/𝜇 1,23Zeitkonstante 𝜏 = 𝑅2/𝑎 6,44 × 1019 sTemperaturkonstante 𝑇1 = 𝑟A𝑅2/𝜅 5,73 × 105 KVerschiebungskonstante 𝑢0 = 𝜅𝑅/(𝛼𝜇𝑎(3��+2)) 6 × 106 mImpulskopplungsfaktor 𝑚 = 𝛼2𝑅2𝜇𝑎𝑟A(3��+2)2

/𝜅2 34,6Hintergrundstrahlung 𝑇ext = 𝑇ext/𝑇1 6,98 × 10−3 KEinstrahlungsparameter 𝑆 = 𝜀𝑆𝑅/(π𝜅𝑇1) 2,05 × 103 W m−2

Strahlungsparameter �� = 𝜀𝜎𝑟3A𝑅7/𝜅4 9,29 × 1016



2.6 Quasistationärer Deformationsprozess

Nachstehend wird ein Argument verwendet, dass an [Dreyer 2013] angelehnt ist.In den normierten Feldgleichungen (2.106) und (2.107) stehen die zwei Faktoren

𝑚 = 𝛼2𝑅2𝜇𝑎𝑟A𝜅2

(3��+ 2

)2und 𝑚1 := 𝜌𝑎2

𝑅2𝜇. (2.118)

Mit den Werten aus Tabelle 2.1 ergeben sich

𝑚 ≈ 34,602 und 𝑚1 = 8,343 × 10−34 . (2.119)

Der Faktor 𝑚1 ist sehr klein. Es erscheint daher sinnvoll den Grenzwert 𝑚1 → 0zu betrachten, woraus die stationäre Feldgleichung der Verschiebungen,

Δ�� +(��+ 1

) (�� · ∇ ⊗ ∇

)−𝑚

(𝑇 − 𝑇R

)∇ = 0 , (2.120)

folgt.

Der Grund für die Größe von 𝑚1 liegt in der Wahl der Zeitskalierung. Diese istan die Änderungsrate des Temperaturfeldes angepasst. Elastische Prozesse laufenin Festkörpern meist sehr viel schneller als thermische Prozesse ab. Daher ist aufder Zeitskala, auf der sich das Temperaturfeld in der Größenordnung 1 ändert,die Änderung der Verschiebungen um viele Größenordnungen höher als die desTemperaturfeldes. Das elastische Feld kann dem Einfluss des Temperaturfeldesnahezu instantan folgen und der Deformationsprozess läuft quasi stationär ab.

Quasistationärer Deformationsprozess

26

3 Analytische Lösungen

Das vollständige, thermoelastische Problem kann nur in speziellen Fällen allgemein gelöst werden. Im Folgenden werden analytische Lösungen für das reinthermische Problem präsentiert. Dazu gehören das radialsymmetrische, stationäre und instationäre Wärmeleitungsproblem mit Robin-Randbedingungund das radialsymmetrische, stationäre Wärmeleitungsproblem mit Stefan-Boltzmann-Randbedingung, wobei die analytischen Lösungen für dasWärmeleitproblem mit Robin zum Verifizieren der Numerik verwendet werden.Den Abschluss bildet eine radialsymmetrische, stationäre Lösung des thermoelastischen Problems mit Stefan-Boltzmann-Randbedingung.

3.1 Radialsymmetrische, stationäre Wärmeleitung mitROBIN-Randbedingung

Wird die innere Reibung vernachlässigt, ist das radialsymmetrische, stationäreProblem mit Robin-Randbedingung gegeben durch

Δ𝑇 (𝑟) + 𝛾(𝑟) = 0 , 𝑥 ∈ Ω , (3.1)�� · 𝑛 = ℎ𝑇 , 𝑥 ∈ 𝜕Ω . (3.2)

Dabei ist 𝑟 =√��2

1 + ��22 + ��2

3 . Ausgeschrieben erhält man

1𝑟2

dd𝑟

(𝑟2 d

d𝑟 𝑇 (𝑟))

+ 𝛾(𝑟) = 0 , (3.3)

d𝑇d𝑟 (1) + ℎ𝑇 (1) = 0 . (3.4)

Die Feldgleichung (3.3) kann unbestimmt integriert werden, um

𝑇 (𝑟) = −ˆ 1𝑟2

ˆ𝑟2𝛾(𝑟) d𝑟 d𝑟 + 𝐶1

𝑟+ 𝐶 (3.5)

zu erhalten. Eine physikalisch sinnvolle Lösung ist in 𝑟 = 0 endlich. Daraus folgt𝐶1 = 0. Die verbleibende Konstante 𝐶 wird aus der Randbedingung gewonnen. Zudiesem Zweck wird die unbestimmte Integration über die normierte Wärmequelle

Kapitel 3. Analytische Lösungen


in einer Funktion 𝛽(𝑟) zusammgefasst:

𝛽(𝑟) := −ˆ 1𝑟2

ˆ𝑟2𝛾(𝑟) d𝑟 d𝑟 . (3.6)

Einsetzen von Gleichung (3.5) in die Randbedingung (3.4) ergibt eine Gleichungfür 𝐶:

d𝛽d𝑟 (1) + ℎ (𝛽(1) + 𝐶) = 0 . (3.7)

Dies ist eine lineare Gleichung in 𝐶, mit genau einer Lösung. Die physikalischeLösung 𝑇 erhält man aus dem entdimensionalisierten Temperaturfeld 𝑇 aus

𝑇 (𝑟) = 𝑇1𝑇 (𝑟) + 𝑇ext . (3.8)

Als Beispiel sei 𝛾 = 1. Dann ist die Lösung des normierten Temperaturfeldesgegeben durch

𝑇 (𝑟) = 𝐶 − 𝑟2

6 . (3.9)

Aus der Randbedingung folgt für 𝐶 die Gleichung

ℎ

(𝐶 − 1

6

)− 1

3 = 0 ⇒ 𝐶 = 13

(12 + 1

ℎ

). (3.10)

Damit liegt sowohl das normierte Temperaturfeld

𝑇 (𝑟) = 16(1 − 𝑟2

)+ 1

3ℎ(3.11)

als auch das physikalische

𝑇 (𝑟) = 𝑇1

(16(1 − 𝑟2

)+ 1

3ℎ

)+ 𝑇ext , (3.12)

mit den normierten Parametern

𝑇1 = 𝑅2𝑟A𝜅

und ℎ = ℎ𝑅

𝜅, (3.13)

fest. Ein exemplarischer Graph der Gleichung (3.11) für den Parameter ℎ = 1 istin Abbildung 3.1 zu finden.

Radialsymmetrische, stationäre Wärmeleitung mit ROBIN-Randbedingung

28

0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1

0,35

0,4

0,45

0,5

Normierter Radius 𝑟/𝑅

Nor

mie

rte

Tem

pera

tur

(𝑇−

𝑇ex

t) /𝑇

1

Abb. 3.1: Beispielhafte Lösung des normierten Temperaturfeldes (3.5) für den Fall 𝛾 = 1und ℎ = 1.

3.2 Radialsymmetrische, instationäre Wärmeleitung mitROBIN-Randbedingung

Die nachfolgende Analyse ist [Hahn u. Özişik 2012] entnommen. Das normierte, radialsymmetrische, instationäre Wärmeleitungsproblem mit Robin-Randbedingung ist gegeben durch

𝜕𝑇

𝜕𝑡(𝑟, 𝑡) − Δ𝑇 (𝑟, 𝑡) = 𝛾(𝑟) , 𝑥 ∈ Ω , (3.14)

�� · 𝑛 = ℎ𝑇 , 𝑥 ∈ 𝜕Ω , (3.15)𝑇 (𝑟, 0) = 𝑇A(𝑟) . (3.16)

Zunächst wird das Temperaturfeld in einen instationären, homogenen Anteil𝑇H(𝑟, 𝑡) und einen stationären, inhomogenen Anteil 𝑇p(𝑟) zerlegt:

𝑇 (𝑟, 𝑡) = 𝑇H(𝑟, 𝑡) + 𝑇p(𝑟) . (3.17)

Der homogene Anteil erfülle die instationäre, homogene Wärmeleitgleichung, derpartikuläre Anteil die stationäre, inhomogene Wärmeleitgleichung. Damit gelten

𝜕𝑇𝐻

𝜕𝑡= Δ𝑇H , (3.18)

Δ𝑇p + 𝛾 = 0 . (3.19)



Für 𝑇H wird ein Separationsansatz verwendet:

𝑇H(𝑟, 𝑡) = 𝑉 (𝑟)𝑍(𝑡) . (3.20)

Einsetzen führt auf

𝑍 ′(𝑡)𝑍(𝑡) = Δ𝑉 (𝑟)

𝑉 (𝑟) , (3.21)

wobei mit 𝑍 ′ die Ableitung von 𝑍 nach 𝑡 gemeint ist. Da beide Seiten Funktionenunabhängiger Koordinaten sind, müssen sowohl 𝑍′(𝑡)/𝑍(𝑡), als auch Δ𝑉 (𝑟)/𝑉 (𝑟)konstant sein. Sei �� ∈ R konstant, dann erhält man ein System entkoppelter,gewöhnlicher Differentialgleichungen für die gesuchten Funktionen 𝑉 und 𝑍:

Δ𝑉 (𝑟) = ��𝑉 (𝑟) , (3.22)𝑍 ′(𝑡) = ��𝑍(𝑡) . (3.23)

Die Differentialgleichung für 𝑉 lautet ausgeschrieben

d2𝑉

d𝑟2 + 2𝑟

d𝑉d𝑟 − ��𝑉 = 0 . (3.24)

Diese Gleichung wird mit

𝑈(𝑟) = 𝑟𝑉 (𝑟) . (3.25)

weiter umgeschrieben. Dazu wird die zweite Ableitung von 𝑈 nach 𝑟 gebildet

d2𝑈

d𝑟2 = 𝑟d2𝑉

d𝑟2 + 2d𝑉d𝑟 . (3.26)

Multiplizieren von Gleichung (3.24) mit 𝑟 und Einsetzen von d2𝑈/d𝑟2 und 𝑈führt auf

d2𝑈

d𝑟2 (𝑟) − ��𝑈(𝑟) = 0 , (3.27)

wodurch sich die Struktur der Differentialgleichung vereinfacht hat. Zusammenfassend gilt es drei gewöhnliche Differentialgleichungen zu lösen:

1𝑟2

dd𝑟

(𝑟2 d𝑇p

d𝑟

)+ 𝛾 = 0, (3.28)

d𝑍d𝑡 − ��𝑍 = 0 , (3.29)

d2𝑈

d𝑟2 − ��𝑈 = 0 . (3.30)

Radialsymmetrische, instationäre Wärmeleitung mit ROBIN-Randbedingung

30

Als nächstes werden Randbedingungen für 𝑈(𝑟) und 𝑇H formuliert. Die Lösungfür 𝑇H wird mit einem Separationsansatz gesucht, weshalb die Randbedingunghomogen sein muss. Einsetzen von Gleichung (3.17) in Gleichung (3.15) ergibt

𝜕𝑇H𝜕𝑟

(1, 𝑡) + ℎ𝑇H(1, 𝑡) = −(

d𝑇pd𝑟 (1) + ℎ𝑇p(1)

). (3.31)

Die Randbedingung für 𝑇H ist genau dann homogen, wenn die partikuläre Lösung

d𝑇pd𝑟 (1) + ℎ𝑇p(1) = 0 (3.32)

erfüllt. Entsprechend gilt dann für die homogene Lösung

𝜕𝑇H𝜕𝑟

(1, 𝑡) + ℎ𝑇H(1, 𝑡) = 0 . (3.33)

Einsetzen des Separationsansatzes aus Gleichung (3.20) liefert( dd𝑟

(𝑈

𝑟

)(1) + ℎ𝑈(1)

)𝑍(𝑡) = 0 . (3.34)

Diese Bedingung muss für alle Zeiten erfüllt sein. Dann ist entweder 𝑍 = 0, wasauf die triviale Lösung führt, oder der Faktor vor 𝑍 muss verschwinden:

dd𝑟

(𝑈

𝑟

)(1) + ℎ𝑈(1) = 0 , (3.35)

beziehungsweise ausdifferenziert:

d𝑈d𝑟 (1) + (ℎ− 1)𝑈(1) = 0 . (3.36)

An dieser Stelle ergibt sich eine weitere Bedingung an 𝑈(𝑟). Es wird verlangt,dass 𝑇 (𝑟, 𝑡) und damit auch 𝑇H(𝑟, 𝑡) auf ganz Ω beschränkt ist. Da 𝑅(𝑟) = 𝑈(𝑟)/𝑟

gesetzt wurde, kann das Temperaturfeld überhaupt nur dann beschränkt sein,falls die notwendige Bedingung 𝑈(𝑟) → 0 für 𝑟 → 0 erfüllt ist. Man erhält einezusätzliche Bedingung

𝑈(0) = 0 . (3.37)

Zunächst wird die homogene Lösung 𝑇H gesucht. Dazu muss Gleichung (3.30)gelöst werden, wobei sich drei Fälle ergeben:

�� = 0 , �� > 0 , �� < 0 , (3.38)

mit den jeweiligen Lösungen 𝑈0, 𝑈+ und 𝑈−.



Zuerst wird �� = 0 untersucht. Die Lösung für 𝑈0 ist

𝑈0(𝑟) = 𝐶1 + 𝐶2𝑟. (3.39)

Aus 𝑈0(0) = 0 folgt 𝐶1 = 0. Einsetzen in die Randbedingung (3.36) ergibt

ℎ𝐶2 = 0 , (3.40)

was für ℎ ∈ R+ nur durch 𝐶2 = 0 gelöst wird, wodurch 𝑈0 die Nullfunktion ist.

Für den Fall �� > 0 wird �� = 𝜆2 mit 𝜆 ∈ R+ gesetzt. Die Lösung ist der Gleichung(3.30) ist

𝑈+(𝑟) = 𝐶1 sinh(𝜆𝑟) + 𝐶2 cosh(𝜆𝑟) . (3.41)

Aus 𝑈+(0) = 0 folgt 𝐶2 = 0. Einsetzen in die Randbedingung (3.36) ergibt

𝐶1(𝜆 cosh(𝜆) + (ℎ− 1) sinh(𝜆)

)= 0 , (3.42)

womit 𝐶1 = 0 oder die charakteristische Gleichung

𝜆 cosh(𝜆) + (ℎ− 1) sinh(𝜆) = 0 (3.43)

folgt. Sie hat auf R+ keine Nullstelle, solange ℎ > 0 erfüllt ist. Für ℎ = 1 erhältman den Fall 𝜆 = 0, welcher ebenfalls ausgeschlossen ist.

Im dritten Fall �� < 0 wird �� = −𝜆2 mit 𝜆 ∈ R+ gesetzt. Die Lösung derGleichung (3.30) lautet

𝑈−𝑟 = 𝐶1 sin(𝜆𝑟) + 𝐶2 cos(𝜆𝑟) . (3.44)

Bedingung (3.37) führt auf 𝐶2 = 0. Aus der Robin Randbedingung (3.36) folgt

𝜆 cos(𝜆) + (ℎ− 1) sin(𝜆) = 0. (3.45)

Für ℎ ∈ R+ besitzt diese unendlich viele Lösungen 𝜆𝑛 ∈ R+.

Als nächstes wird die Lösung für 𝑍(𝑡) bestimmt, wobei als einzig relevanter Fall−𝜆2 = �� < 0 übrig bleibt. Gleichung (3.29) wird durch

𝑍(𝑡) = 𝐶𝑡 exp(−𝜆2𝑡

)(3.46)

für alle 𝜆 ∈ R+ gelöst.

Für die partikuläre Lösung muss Gleichung (3.28) gelöst werden, was im vorheri


32

gen Abschnitt bereits geschehen ist. Die Lösung für 𝑇p ist gegeben durch

𝑇p = 𝛽(𝑟) + 𝐶 , (3.47)d𝛽d𝑟 (1) + ℎ (𝛽(1) + 𝐶) = 0 . (3.48)

Die Gesamtlösung 𝑇 lautet somit

𝑇 (𝑟, 𝑡) = 𝑇H(𝑟, 𝑡) + 𝑇p(𝑟) =∞∑

𝑛=1𝐶𝑛

sin(𝜆𝑛𝑟)𝑟

exp(−𝜆2

𝑛𝑡)

+ 𝑇p(𝑟) . (3.49)

Die Lösung (3.49) soll an die Anfangsbedingung angepasst werden. Es gilt

𝑇 (𝑟, 0) =∞∑

𝑛=1𝐶𝑛

sin(𝜆𝑛𝑟)𝑟

+ 𝑇𝑝(𝑟) = 𝑇A(𝑟) , (3.50)

beziehungsweise∞∑

𝑛=1𝐶𝑛 sin(𝜆𝑛𝑟) =

(𝑇A(𝑟) − 𝑇p(𝑟)

)𝑟. (3.51)

Um die Koeffizienten 𝐶𝑛 zu bestimmen, wird eine Orthogonalitätsrelation dersin(𝜆𝑛𝑟) ausgenutzt. Sei 𝑈𝑛(𝑟) = sin(𝜆𝑛𝑟) mit

𝑈 ′′𝑛(𝑟) + 𝜆2

𝑛𝑈𝑛(𝑟) = 0, 𝑈𝑛(0) = 0 und 𝑈 ′𝑛(1) + (ℎ− 1)𝑈𝑛(1) = 0 .

(3.52)

Dabei bezeichnet (·)′ die Ableitung nach dem Argument von 𝑈 . Dann gilteinerseits

1ˆ

𝑟=0

(𝑈 ′′

𝑛(𝑟)𝑈𝑘(𝑟) − 𝑈 ′′𝑘 (𝑟)𝑈𝑛(𝑟)

)d𝑟 =

1ˆ

𝑟=0

(−𝜆2

𝑛𝑈𝑛(𝑟)𝑈𝑘(𝑟) + 𝜆2𝑘𝑈𝑛𝑈𝑘

)d𝑟

=(𝜆2

𝑘 − 𝜆2𝑛

) 1ˆ

𝑟=0

𝑈𝑛(𝑟)𝑈𝑘(𝑟) d𝑟

und andererseits1ˆ

𝑟=0

(𝑈 ′′

𝑛(𝑟)𝑈𝑘(𝑟) − 𝑈 ′′𝑘 (𝑟)𝑈𝑛(𝑟)

)d𝑟 =

1ˆ

𝑟=0

(𝑈 ′

𝑛(𝑟)𝑈𝑘(𝑟) − 𝑈 ′𝑘(𝑟)𝑈𝑛(𝑟)

)′ d𝑟

=[𝑈 ′

𝑛(𝑟)𝑈𝑘(𝑟) − 𝑈𝑛(𝑟)𝑈 ′𝑘(𝑟)

]1𝑟=0 = (ℎ− 1) (𝑈𝑛(𝑟)𝑈𝑘(𝑟) − 𝑈𝑛(𝑟)𝑈𝑘(𝑟)) = 0



woraus man

(𝜆2

𝑘 − 𝜆2𝑛

) 1ˆ

𝑟=0

𝑈𝑛(𝑟)𝑈𝑘(𝑟) d𝑟 = 0 (3.53)

folgert. Das ist für 𝜆𝑘 = 𝜆𝑛 erfüllt, falls

1ˆ

𝑟=0

𝑈𝑛(𝑟)𝑈𝑘(𝑟) d𝑟 = 0 (3.54)

für alle 𝑛 = 𝑘 gilt. Für 𝑛 = 𝑘 findet sich

1ˆ

𝑟=0

𝑈2𝑛(𝑟) d𝑟 =

(12 − sin(2𝜆𝑛)

4𝜆𝑛

):= 𝑁𝑛. (3.55)

Multiplikation von Gleichung (3.51) mit sin(𝜆𝑘𝑟) und Integration über Ω ergibt

𝐶𝑘𝑁𝑘 =1ˆ

𝑟=0

(𝑇A(𝑟) − 𝑇p(𝑟)

)𝑟 sin(𝜆𝑘𝑟) d𝑟 . (3.56)

An dieser Stelle muss die partikuläre Lösung 𝑇p bekannt sein, um die Koeffizienten𝐶𝑘 explizit bestimmen zu können.

Die Lösung des transienten, radialsymmetrischen Wärmeleitungsproblems mitRobinRandbedingung ist daher gegeben durch

𝑇 (𝑟, 𝑡) = 𝑇1

( ∞∑𝑛=1

𝐶𝑛𝑅

𝑟sin(𝜆𝑛

𝑟

𝑅

)exp

(−𝜆2

𝑛

𝑡

𝜏

)+ 𝑇p (𝑟/𝑅)

)+ 𝑇ext (3.57)

mit

𝜏 = 𝑅2

𝛼, 𝑇1 = 𝑅2𝑟A

𝜅, 𝜆𝑘 cos(𝜆𝑘) +

(𝑅ℎ

𝜅− 1

)sin(𝜆𝑘) = 0 , (3.58)

1𝑟2

dd𝑟

(𝑟2 d𝑇p

d𝑟 (𝑟))

+ 𝛾(𝑟) = 0 und d𝑇pd𝑟 (1) + ℎ𝑇p(1) = 0 . (3.59)

Als Beispiel wird der Fall ℎ = 1, 𝛾 = 1 und 𝑇A ∈ R konstant betrachtet. Diepartikuläre Lösung 𝑇p ist nach Gleichung (3.11) gegeben durch

𝑇p(𝑟) = 12

(1 − 𝑟2

3

). (3.60)


34

Die charakteristische Gleichung lautet

𝜆 cos(𝜆) = 0 (3.61)

woraus die Eigenwerte

𝜆𝑛 = π

(12 + 𝑛

), 𝑛 ∈ N0 (3.62)

folgen. Die Funktionennorm ist

𝑁𝑛 =1ˆ

𝑟=0

sin(π

(12 + 𝑛

)𝑟

)2d𝑟 = 1

2 , (3.63)

woraus die Koeffizienten

𝐶𝑛 = 1𝑁𝑛

1ˆ

𝑟=0

(𝑇A − 1

2

(1 − 𝑟2

3

))𝑟 sin

(π

(12 + 𝑛

)𝑟

)d𝑟

=2 (−1)𝑛

(π2(

12 + 𝑛

)2𝑇A − 1

)π4(

12 + 𝑛

)4 (3.64)

und schließlich das normierte Temperaturfeld

𝑇 (𝑟 , 𝑡) =∞∑

𝑛=0

2 (−1)𝑛(π2(

12 + 𝑛

)2𝑇A − 1

)π4(

12 + 𝑛

)4𝑟

sin(π

(12 + 𝑛

)𝑟

)

exp(

−π2(1

2 + 𝑛

)2𝑡

)+ 1

2

(1 − 𝑟2

3

)(3.65)

folgen. Das entdimensionalisierte Temperaturfeld (3.65) ist für einige exemplarische Zeitpunkte in Abbildung 3.2 dargestellt.



0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1

0,4

0,6

0,8

1


Nor

mie

rte

Tem

pera

tur

(𝑇−

𝑇ex

t) /𝑇

1 𝑡 = 0𝑡 = 0,125𝑡 = 0,25𝑡 = 0,5𝑡 = 1𝑡 = 2stationär

Abb. 3.2: Beispielhafte Lösung des normierten Temperaturfeldes (3.65) für den Fall𝛾 = 1, ℎ = 1 und 𝑇A = 1, wobei die Reihe nach 100 Gliedern abgebrochen wurde. ZumVergleich ist die stationäre Lösung (3.11) mit eingezeichnet.

3.3 Radialsymmetrische, stationäre Wärmeleitung mitSTEFAN-BOLTZMANN-Randbedingung

Das radialsymmetrische, stationäre Wärmeleitungsproblem mit Stefan-Boltzmann-Randbedingung lautet

Δ𝑇 (𝑟) + 𝛾(𝑟) = 0 , 𝑥 ∈ Ω , (3.66)

�� · 𝑛 = ��(𝑇 4 − 𝑇 4

ext

), 𝑥 ∈ 𝜕Ω , (3.67)

beziehungsweise ausgeschrieben:

1𝑟2

dd𝑟

(𝑟2 d𝑇

d𝑟 (𝑟))

+ 𝛾(𝑟) = 0 , 𝑥 ∈ Ω , (3.68)

d𝑇d𝑟 + ��

(𝑇 4 − 𝑇 4

ext

)= 0 , 𝑥 ∈ 𝜕Ω . (3.69)

Nach Gleichung (3.5) ist die Lösung der Feldgleichung gegeben durch

𝑇 (𝑟) = 𝛽(𝑟) + 𝐶 . (3.70)

Radialsymmetrische, stationäre Wärmeleitung mit STEFAN-BOLTZMANN-Randbedingung

36

Dabei ist die Integrationskonstante 𝐶 wieder aus der Randbedingung zu bestimmen. Durch Einsetzen erhält man

(𝛽(1) + 𝐶)4 + d𝛽d𝑟 (1) = 0 . (3.71)

Dies ist eine algebraische Gleichung vierter Ordnung in 𝐶. Sie besitzt damit vier,unter Umständen komplex konjugierte Lösungen. Die Auswahl des physikalischkorrekten 𝐶 erfolgt über die Forderungen

𝐶 ∈ R und 𝑇 (𝑟) > 0 ∀ 𝑟 ∈ [0, 1] . (3.72)

Im Folgenden wird die Lösung exemplarisch für den Fall 𝛾 = 1 berechnet.Unbestimmtes Integrieren von 𝛾 führt auf 𝛽:

𝛽(𝑟) = −ˆ 1𝑟2

ˆ𝑟2𝛾(𝑟) d𝑟 d𝑟 = −𝑟2

6 . (3.73)

Die charakteristische Gleichung für 𝐶 lautet damit

𝐶4 − 23𝐶

3 + 16𝐶

2 − 154𝐶 + 1

1296 − 𝑇 4ext − 1

3�� = 0. (3.74)

Sie besitzt die vier Lösungen

𝐶1,2 = 16 ± i

( 13��(1 + 3𝑇 4

ext��)) 1

4, (3.75)

𝐶3,4 = 16 ±

( 13��(1 + 3𝑇 4

ext��)) 1

4. (3.76)

Die Lösungen 𝐶1,2 sind komplex und damit unphysikalisch. Die Auswahl zwischen𝐶3,4 ist komplizierter. Für eine Entscheidung wird das normierte Temperaturfeldbetrachtet,

𝑇 (𝑟) = 𝛽(𝑟) + 𝐶 = 𝐶 − 𝑟2

6 , (3.77)

welches ein Minimum für 𝑟 = 1 besitzt:

𝑇 (1) = 𝐶 − 16 . (3.78)

Damit die Bedingung 𝑇 > 0 erfüllt ist, muss 𝐶 > 1/6 gelten. Offensichtlich istdas für

𝐶3 = 16 +

( 13��(1 + 3𝑇 4

ext��)) 1

4(3.79)



0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1

1,1

1,15

1,2

1,25


Nor

mie

rte

Tem

pera

tur

𝑇/𝑇

1

Abb. 3.3: Beispielhafte Lösung des normierten Temperaturfeldes (3.81) für den Fall𝛾 = 𝑇ext = �� = 1. Im Gegensatz zum normierten, linearen Problem mit Robin-Randbedingung liegt hier eine dimensionslose, absolute Temperaturskala vor.

immer erfüllt. Hingegen für

𝐶4 = 16 −

( 13��(1 + 3𝑇 4

ext��)) 1

4(3.80)

entspricht das der Forderung nach komplexen Materialparametern. Damit ist 𝐶4immer eine unphysikalische Lösung.

Somit erhält man das normierte,

𝑇 (𝑟) =( 1

3��(1 + 3𝑇 4

ext��)) 1

4+ 1

6(1 − 𝑟2

), (3.81)

und das physikalische Temperaturfeld:

𝑇 (𝑟) = 𝑇1

[( 13��(1 + 3𝑇 4

ext��)) 1

4+ 1

6

(1 −

(𝑟

𝑅

)2)]

. (3.82)

Das normierte Temperaturfeld (3.81) ist in Abbildung 3.3 für den Fall 𝑇ext =�� = 1 dargestellt.

Radialsymmetrische, stationäre Wärmeleitung mit STEFAN-BOLTZMANN-Randbedingung

38

3.4 Stationäres, thermoelastisches Problem mit STEFAN-BOLTZMANN-Randbedingung

Das Problem auf Ω ist gegeben durch

𝜕𝑇

𝜕𝑡− Δ𝑇 = 𝛾 , (3.83)

Δ�� +(��+ 1

) (�� · ∇ ⊗ ∇

)−𝑚

(𝑇 − 𝑇R

)∇ = 0 (3.84)

mit den Randbedingungen auf 𝜕Ω:

�� · 𝑛 = ��(𝑇 4 − 𝑇 4

ext

), (3.85)

𝑇 · 𝑛 = 0 . (3.86)

Zunächst wird die Impulsbilanz, Gleichung (3.84), mit der Identität

�� · ∇ ⊗ ∇ = Δ�� + ∇ × ∇ × �� (3.87)

umgeschrieben:(��+ 2

)Δ�� +

(��+ 1

) (∇ × ∇ × ��

)−𝑚

(𝑇 − 𝑇R

)∇ = 0 . (3.88)

Es wird der radialsymmetrische Fall mit konstanter Referenztemperatur betrachtet und ein semi-inverser Ansatz für die Verschiebungen gewählt:

��(𝑟) = 𝑢𝑟(𝑟)𝑒𝑟 . (3.89)

Damit ist das Verschiebungsfeld wirbelfrei:

∇ × �� = 1𝑟 sin(𝜗)

(𝜕��𝜙 sin(𝜗)

𝜕𝜗− 𝜕��𝜗

𝜕𝜙

)𝑒𝑟+

+ 1𝑟

( 1sin(𝜗)

𝜕��𝑟

𝜕𝜙− 𝜕𝑟��𝜙

𝜕𝑟

)𝑒𝜗 + 1

𝑟

(𝜕𝑟��𝜗

𝜕𝑟− 𝜕��𝑟

𝜕𝜗

)𝑒𝜙 = 0 . (3.90)

Setzt man den Ansatz in die Verschiebungsgleichung ein und streicht den Rotationsterm, erhält man (

��+ 2)

Δ�� −𝑚𝑇 ∇ = 0 . (3.91)

Der vektorielle Laplaceoperator lautet ausgeschrieben

Δ�� =(

Δ��𝑟 − 2��𝑟

𝑟2

)𝑒𝑟 =

(d2��𝑟

d𝑟2 + 2𝑟

d��𝑟

d𝑟 − 2��𝑟

𝑟2

)𝑒𝑟 , (3.92)



woraus die Feldgleichung für das Verschiebungsfeld folgt:((��+ 2

)(d2��𝑟

d𝑟2 + 2𝑟

d��𝑟

d𝑟 − 2��𝑟

𝑟2

)−𝑚

d𝑇d𝑟

)𝑒𝑟 = 0 . (3.93)

Damit ist die Feldgleichung für die radiale Komponente ��𝑟 gegeben durch(d2��𝑟

d𝑟2 + 2𝑟

d��𝑟

d𝑟 − 2��𝑟

𝑟2

)− 𝑚

��+ 2d𝑇d𝑟 = 0 . (3.94)

Im Folgenden wird der Fall 𝛾 = 1 betrachtet. Die Lösung des Wärmeleitproblemsist gegeben durch

𝑇 = 𝐶 − 𝑟2

6 mit 𝐶 = 16 +

( 13�� + 𝑇 4

ext

) 14. (3.95)

Bei bekanntem Temperaturfeld ist die Verschiebungsgleichung eine gewöhnliche,inhomogene Differentialgleichung. Die Lösung ist durch einen partikulären undeinen homogenen Anteil gegeben. Als Ansatz für die partikuläre Lösung wird

��𝑟p = 𝐴𝑟3 . (3.96)

gewählt. Durch Einsetzen des Ansatzes erhält man

𝐴𝑟 (6 + 6 − 2) = 10𝐴𝑟 = − 𝑚𝑟

3(��+ 2

) ⇒ 𝐴 = − 𝑚

30(��+ 2

) . (3.97)

Für die homogene Lösung wird ein Potenzansatz verwendet:

��𝑟H = 𝑟𝑝 . (3.98)

Einsetzen in die homogene Verschiebungsgleichung führt auf

𝑟𝑝−2 (𝑝(𝑝− 1) + 2𝑝− 2) = 0 (3.99)

was für die zwei Werte

𝑝1 = 1 und 𝑝2 = −2 (3.100)

erfüllt ist. Damit erhält man für die homogene Lösung

��𝑟H = 𝐶1𝑟 + 𝐶2𝑟2 . (3.101)

Die Verschiebungen sollen für 𝑟 → 0 endlich sein, woraus 𝐶2 = 0 folgt. Die

Stationäres, thermoelastisches Problem mit STEFAN-BOLTZMANN-Randbedingung

40

Konstante 𝐶1 wird aus der Randbedingung (3.86) gewonnen. Dabei gelten

𝑇 =(��Sp

(��)

−𝑚(𝑇 − 𝑇𝑟

))1 + 2�� (3.102)

und

�� = sym(�� ⊗ ∇

). (3.103)

Für die Verschiebung wurde

��(𝑟) = ��𝑟(𝑟)𝑒𝑟 =

⎛⎝𝐶1𝑟 − 𝑚

30(��+ 2

)𝑟3

⎞⎠ 𝑒𝑟 (3.104)

gefunden. Damit folgt der normierte Verzerrungstensor zu

��(𝑟) =

⎛⎝𝐶1 − 𝑚𝑟2

30(��+ 2

)⎞⎠ 𝑒𝑟 ⊗ 𝑒𝑟+

+

⎛⎝𝐶1 − 𝑚𝑟2

30(��+ 2

)⎞⎠ (𝑒𝜙 ⊗ 𝑒𝜙 + 𝑒𝜗 ⊗ 𝑒𝜗) (3.105)

und der normierte Spannungstensor:

𝑇 (𝑟) =

⎛⎝5𝐶1 −𝑚

⎛⎝𝑇 (𝑟) − 𝑇𝑟 + 11𝑟2

30(��+ 2

)⎞⎠⎞⎠ 𝑒𝑟 ⊗ 𝑒𝑟+

+

⎛⎝5𝐶1 −𝑚

⎛⎝𝑇 (𝑟) − 𝑇𝑟 + 7𝑟2

30(��+ 2

)⎞⎠⎞⎠ (𝑒𝜙 ⊗ 𝑒𝜙 + 𝑒𝜗 ⊗ 𝑒𝜗) . (3.106)

Das Temperaturfeld 𝑇 wird erst zum Schluss eingesetzt. Der Normalenvektorvon 𝜕Ω ist 𝑒𝑟. Damit gilt

𝑇 (1) · 𝑒𝑟 =

⎛⎝5𝐶1 −𝑚

⎛⎝𝑇 (1) − 𝑇𝑟 + 1130(��+ 2

)⎞⎠⎞⎠ 𝑒𝑟 = 0 , (3.107)

beziehungsweise

5𝐶1 −𝑚

⎛⎝𝑇 (1) − 𝑇𝑟 + 1130(��+ 2

)⎞⎠ = 0 , (3.108)



0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1−0,002

0

0,002


Nor

mie

rte

Vers

chie

bung

𝑢𝑟 /

𝑢0

Abb. 3.4: Lösung der normierten, radialen Komponente des Verschiebungsfeldes (3.4)für den Fall 𝛾 = 𝑇ext = �� = 𝑚 = �� = 1 und 𝑇R = 1,15. Das zugehörige Temperaturfeldwurde im vorangegangenen Abschnitt berechnet und ist in Abbildung (3.3) dargestellt.Es treten Bereiche auf, die sowohl unter, als auch oberhalb der Referenztemperaturliegen. Entsprechend ist die Kugel teilweise auf Zug und teilweise auf Druck belastet,wobei die kalten Regionen näher am Rand liegen.

woraus sich die Integrationskonstante ergibt:

𝐶1 = 𝑚

5

⎛⎝𝑇 (1) − 𝑇𝑟 + 1130(��+ 2

)⎞⎠ . (3.109)

Das Verschiebungsfeld ist damit

��𝑟(𝑟) =

⎛⎝𝑚5

⎛⎝𝑇 (1) − 𝑇𝑟 + 1130(��+ 2

)⎞⎠ 𝑟 − 𝑚

30(��+ 2

)𝑟3

⎞⎠ 𝑒𝑟

=

⎛⎝𝑚5

⎛⎝( 13�� + 𝑇 4

ext

) 14

− 𝑇𝑟 + 1130(��+ 2

)⎞⎠ 𝑟 − 𝑚

30(��+ 2

)𝑟3

⎞⎠ 𝑒𝑟 . (3.110)

Die Gleichung (3.4) ist in Abbildung 3.4 für den Fall 𝑇ext = �� = 𝑚 = �� = 1 und𝑇R = 1,15 dargestellt.

Interessanterweise ist die Lösung des Verschiebungsfeldes eindeutig, obwohlein reines Neumann-Problem vorliegt. Der Grund liegt in der Wahl des semi-inversen Ansatzes. Dieser ist so gewählt, dass er radialsymmetrisch ist unddie Verschiebungskomponenten in 𝑒𝜗 oder 𝑒𝜙 Richtung verschwinden, wodurchStarrkörperrotationen und Starrkörpertranslationen ausgeschlossen sind.

Stationäres, thermoelastisches Problem mit STEFAN-BOLTZMANN-Randbedingung

42

4 Numerische Lösungen

In diesem Kapitel werden die Feldgleichungen zuerst für die numerische Analyseaufbereitet und anschließend ausgewählte Ergebnisse präsentiert. Dabei werdendrei Randbedingungen für das thermische und eine Randbedingung für daselastische Problem betrachtet.

Die Feldgleichungen werden in ihre schwache Form überführt, welche mit einerFinite-Elemente-Methode gelöst werden. Die Diskretisierung geschieht mit einer Galerkin-Methode mit Lagrange-Polynomen erster Ordnung. Für dieZeitdomäne wird das implizite Euler-Verfahren verwendet. Der Mesher, sowieder Löser, wurden von der Python Bibliothek FEniCS, [Logg u. a. 2012], bereitgestellt.

Falls nicht anders vermerkt, haben die numerischen Ergebnisse des thermischenProblems 25 203 Freiheitsgrade, während das elastische Problem 75 609 Freiheitsgrade besitzt.

Zudem wird ein Algorithmus vorgestellt, um das gekoppelte, thermoelastischeProblem zu lösen.

4.1 Schwache Form der Feldgleichungen

Ausgangspunkt sind die Feldgleichungen

𝜕𝑇

𝜕𝑡− Δ𝑇 +

(𝑇 − 𝑇R

)(𝜕��

𝜕𝑡· ∇)

− 𝛾 = 0 , (4.1)

Δ�� +(��+ 1

) (�� · ∇ ⊗ ∇

)−𝑚

(𝑇 − 𝑇R

)∇ = 0 . (4.2)

Zunächst wird die Wärmeleitgleichung in die schwache Form überführt. Multiplikation mit einer einmal schwach differenzierbaren Testfunktion 𝛿𝑇 undIntegration über Ω führt auf

ˆ

Ω

(𝜕𝑇

𝜕𝑡− Δ𝑇 +

(𝑇 − 𝑇R

)(𝜕��

𝜕𝑡· ∇)

− 𝛾

)𝛿𝑇 d𝑉 = 0 . (4.3)

Kapitel 4. Numerische Lösungen


Der Laplaceoperator kann mit der Produktregel umgeschrieben werden:(Δ𝑇

)𝛿𝑇 = −

(��𝛿𝑇

)· ∇ −

(𝑇 ∇

)·(𝛿𝑇 ∇

). (4.4)

Zusammen mit dem Integral über Ω lässt sich der Satz von Gauss-Ostrogradskianwenden:

−ˆ

Ω

(��𝛿𝑇

)· ∇ d𝑉 −

ˆ

Ω

(𝑇 ∇

)·(𝛿𝑇 ∇

)d𝑉

= −˛

𝜕Ω

𝛿𝑇 �� · 𝑛 d𝐴−ˆ

Ω

(𝑇 ∇

)·(𝛿𝑇 ∇

)d𝑉 . (4.5)

Einsetzen ergibt

ˆ

Ω

((𝜕𝑇

𝜕𝑡+(𝑇 − 𝑇R

)(𝜕��

𝜕𝑡· ∇)

− 𝛾

)𝛿𝑇 +

(𝑇 ∇

)·(𝛿𝑇 ∇

))d𝑉+

+˛

𝜕Ω

𝛿𝑇 �� · 𝑛 d𝐴 = 0 . (4.6)

Gleichung (4.6) ist die schwache Form der Feldgleichung für die Temperatur.

Nachfolgend werden die schwachen Formen für die drei Randbedingungen präsentiert. Diese sind:

Robin �� · 𝑛 = ℎ𝑇 , (4.7)

Stefan-Boltzmann �� · 𝑛 = ��(𝑇 4 − 𝑇 4

ext

), (4.8)

Sonneneinstrahlung �� · 𝑛 = ��𝑇 4 − 𝑆 sin(𝜗) . (4.9)

Eingesetzt in Gleichung (4.6) erhält man die schwache Form mit Robin-Randbedingung:

ˆ

Ω

((𝜕𝑇


)(𝜕��

𝜕𝑡· ∇)

− 𝛾

)𝛿𝑇 +

(𝑇 ∇

)·(𝛿𝑇 ∇

))d𝑉+

+˛

𝜕Ω

ℎ𝛿𝑇𝑇 d𝐴 = 0 , (4.10)

Schwache Form der Feldgleichungen

44

für die Stefan-Boltzmann-Randbedingung:

ˆ

Ω

((𝜕𝑇


)(𝜕��

𝜕𝑡· ∇)

− 𝛾

)𝛿𝑇 +

(𝑇 ∇

)·(𝛿𝑇 ∇

))d𝑉+

+˛

𝜕Ω

𝛿𝑇 ��(𝑇 4 − 𝑇 4

ext

)d𝐴 = 0 , (4.11)

und für die Stefan-Boltzmann-Randbedingung mit Sonneneinstrahlung:

ˆ

Ω

((𝜕𝑇


)(𝜕��

𝜕𝑡· ∇)

− 𝛾

)𝛿𝑇 +

(𝑇 ∇

)·(𝛿𝑇 ∇

))d𝑉+

+˛

𝜕Ω

𝛿𝑇(��𝑇 4 − 𝑆 sin(𝜗)

)d𝐴 = 0 . (4.12)

Als nächstes wird die Feldgleichung für die Verschiebungen in die schwache Formgebracht. Es ist geschickt, mit der Form

𝑇 · ∇ = 0 (4.13)

zu starten. Multiplikation mit einem einmal schwach differenzierbaren Testfeld𝛿�� und Integration über das Gebiet Ω ergibt

ˆ

Ω

(𝑇 · ∇

)· 𝛿�� d𝑉 = 0 . (4.14)

Abermaliges Anwenden der Produktregel liefert(𝑇 · ∇

)· 𝛿�� =

(𝛿�� · 𝑇

)· ∇ − 𝑇 ··

(𝛿�� ⊗ ∇

). (4.15)

Anwenden des Integralsatzes von Gauss-Ostrogradski auf den ersten Termergibt die gesuchte Identität:

ˆ

Ω

(𝑇 · ∇

)· 𝛿�� d𝑉 =

˛

𝜕Ω

(𝛿�� · 𝑇

)· 𝑛 d𝐴−

ˆ

Ω

𝑇 ··(𝛿�� ⊗ ∇

)d𝑉 , (4.16)

beziehungsweise mit eingesetzter Traktion 𝑇 · 𝑛 =: ��:ˆ

Ω

(𝑇 · ∇

)· 𝛿�� d𝑉 =

˛

𝜕Ω

𝛿�� · �� d𝐴−ˆ

Ω

𝑇 ··(𝛿�� ⊗ ∇

)d𝑉 . (4.17)

Für das elastische Problem wird der Rand 𝜕Ω als spannungsfrei angenommen,



daher verbleibt ˆ

Ω

𝑇 ··(𝛿�� ⊗ ∇

)d𝑉 = 0 . (4.18)

Hier ist zu beachten, dass der Spannungstensor eine Funktion der Verschiebungen und des Temperaturfeldes ist. Das Gesetz für den Spannungstensor ausGleichung (2.4) wird in den Integralkern von Gleichung (4.18) eingesetzt:

𝑇 ··(𝛿�� ⊗ ∇

)=(�� Sp

(��)

−𝑚(𝑇 − 𝑇R

))1 ··

(𝛿�� ⊗ ∇

)+ (4.19)

+2�� ··(𝛿�� ⊗ ∇

). (4.20)

Die Spur des Gradienten ist die Divergenz, daher gelten weiterhin

1 ··(𝛿�� ⊗ ∇

)= Sp

(𝛿�� ⊗ ∇

)= 𝛿�� · ∇ und Sp

(��)

= �� · ∇ . (4.21)

Alles gemeinsam eingesetzt ergibt die schwache Form der Verschiebungsgleichung:ˆ

Ω

((��(�� · ∇

)−𝑚

(𝑇 − 𝑇R

)) (𝛿�� · ∇

)+

+ 2 sym(�� ⊗ ∇

)··(𝛿�� ⊗ ∇

))d𝑉 = 0 . (4.22)

Diese wird noch nach Linear- und Bilinearform sortiertˆ

Ω

(��(�� · ∇

) (𝛿�� · ∇

)+ 2 sym

(�� ⊗ ∇

)··(𝛿�� ⊗ ∇

))d𝑉

=ˆ

Ω

𝑚(𝑇𝑖 − 𝑇R

) (𝛿�� · ∇

)d𝑉 . (4.23)

Für die Zeitdomäne wird das implizite Euler-Verfahren verwendet. Gegeben isteine Differentialgleichung der Form

𝜕𝑢

𝜕𝑡(𝑥, 𝑡) = 𝐹 (𝑥, 𝑡, 𝑢) . (4.24)

Zusammen mit dem Zeitschritt 𝑘 und 𝑢𝑖(𝑥) := 𝑢(𝑥, 𝑖𝑘), sowie 𝐹𝑖(𝑥) :=𝐹 (𝑥, 𝑖𝑘, 𝑢𝑖(𝑥)) lautet das implizite Euler-Verfahren für diese Gleichung

𝑢𝑖 − 𝑢𝑖−1𝑘

= 𝐹𝑖 . (4.25)

Dies angewandt auf die schwache Form des thermischen Problems,


46

Gleichung (4.6), führt auf

ˆ

Ω

((𝑇𝑖 − 𝑇𝑖−1

𝑘+(𝑇𝑖 − 𝑇R

)( ��𝑖 − ��𝑖−1𝑘

· ∇)

− 𝛾

)𝛿𝑇 +

(𝑇𝑖∇

)·(𝛿𝑇 ∇

))d𝑉+

+˛

𝜕Ω

𝛿𝑇 ��𝑖 · 𝑛 d𝐴 = 0 . (4.26)

Setzt man wiederum die drei Randbedingungen ein, erhält man die schwacheForm mit Robin-Randbedingung:

ˆ

Ω

((𝑇𝑖 − 𝑇𝑖−1


)( ��𝑖 − ��𝑖−1𝑘

· ∇)

− 𝛾

)𝛿𝑇 +

(𝑇𝑖∇

)·(𝛿𝑇 ∇

))d𝑉+

+˛

𝜕Ω

𝛿𝑇 ℎ𝑇𝑖 d𝐴 = 0 , (4.27)

mit Stefan-Boltzmann-Randbedingung:

ˆ

Ω

((𝑇𝑖 − 𝑇𝑖−1


)( ��𝑖 − ��𝑖−1𝑘

· ∇)

− 𝛾

)𝛿𝑇 +

(𝑇𝑖∇

)·(𝛿𝑇 ∇

))d𝑉+

+˛

𝜕Ω

𝛿𝑇 ��(𝑇 4

𝑖 − 𝑇 4ext

)d𝐴 = 0 , (4.28)

und mit Stefan-Boltzmann-Randbedingung mit Sonneneinstrahlung

ˆ

Ω

((𝑇𝑖 − 𝑇𝑖−1


)( ��𝑖 − ��𝑖−1𝑘

· ∇)

− 𝛾

)𝛿𝑇 +

(𝑇𝑖∇

)·(𝛿𝑇 ∇

))d𝑉+

+˛

𝜕Ω

𝛿𝑇(��𝑇 4

𝑖 − 𝑆 sin(𝜗))

d𝐴 = 0 . (4.29)

Die schwache Form der Robin-Randbedingung wird weiter nach Bilinear- undLinearform sortiert:ˆ

Ω

((𝑇𝑖

𝑘+ 𝑇𝑖

(��𝑖 − ��𝑖−1

𝑘· ∇))

𝛿𝑇 +(𝑇𝑖∇

)·(𝛿𝑇 ∇

))d𝑉+

+˛

𝜕Ω

ℎ𝛿𝑇𝑇𝑖 d𝐴 =ˆ

Ω

(𝑇𝑖−1𝑘

+ 𝑇R

(��𝑖 − ��𝑖−1

𝑘· ∇)

+ 𝛾

)𝛿𝑇 d𝑉 . (4.30)



𝑡𝑖 ≤ 𝑇Simu

𝑇𝑖−1 = 𝑇𝑖

𝑢𝑖−1 = 𝑢𝑖

berechne ��(0)𝑖 mit 𝑇𝑖−1

berechne 𝑇 (0)𝑖 mit 𝑇𝑖−1 und ��

(0)𝑖

berechne ��(1)𝑖 mit 𝑇 (0)

𝑖

berechne 𝑇 (1)𝑖 mit 𝑇𝑖−1 und ��

(1)𝑖

𝑒𝑟𝑟th = ||𝑇 (1)𝑖 − 𝑇

(0)𝑖 ||𝐿2

𝑒𝑟𝑟el = ||��(1)𝑖 − ��

(0)𝑖 ||𝐿2

𝑇(0)𝑖 = 𝑇

(1)𝑖

��(0)𝑖 = ��

(1)𝑖

𝑒𝑟𝑟th > 10−8 oder 𝑒𝑟𝑟el > 10−8

𝑇𝑖 = 𝑇(1)𝑖

��𝑖 = ��(1)𝑖

𝑡𝑖 = 𝑡𝑖 + Δ𝑡

Abb. 4.1: Der verwendete Lösungsalgorithmus für das thermoelastische Problem. DieFelder 𝑇𝑖−1 und ��𝑖−1 sind die Lösungen aus dem letzten Zeitschritt. Die Felder 𝑇𝑖 und��𝑖 stellen die aktuellen Lösungen dar und 𝑇 (0/1)

𝑖−1 und ��(0/1)𝑖−1 bezeichnen die Lösungen des

aktuellen, beziehungsweise des letzten Iterationsschrittes zum Bestimmen der aktuellenFelder.

Das vollständige, thermoelastische Problem besitzt einen bilinearen Term inden Feldern �� und 𝑇 , der das Finden einer Lösung des gekoppelten Problemserschwert.

In Abbildung 4.1 ist ein Struktogramm des verwendeten Algorithmus zum Lösender Feldgleichungen dargestellt. Aus dem gegebenen Temperaturfeld 𝑇𝑖−1 desletzten Zeitschrittes wird ein Ansatz für das aktuelle Verschiebungsfeld ��

(0)𝑖

berechnet. Mit diesem Ansatz wird ein Ansatz für das aktuelle Temperaturfeld𝑇

(0)𝑖 berechnet. Anschließend werden aus diesen Startfeldern iterativ Lösungen der

Feldgleichungen errechnet, bis sich die aktuellen Felder zwischen zwei Iterationenin der 𝐿2(Ω)-Norm um nicht mehr als eine Fehlerschranke ändern, welche hier zu10−8 gesetzt wurde. Die letzte Iteration wird als das aktuelle Verschiebungsfeldgesetzt und die Iteration für den nächsten Zeitschritt begonnen. Dies wird solangwiederholt, bis die Simulationszeit 𝑇Simu erreicht wurde. Im ersten Zeitschritt,i=1, wird 𝑇0 = 𝑇A gesetzt.


48

Für das Verschiebungsfeld liegt ein reines Neumann Problem vor, weshalb imAllgemeinen keine eindeutige Lösung existiert. Wenn Starrkörperrotationen undTranslationen herausgerechnet werden, existiert hier trotzdem eine eindeutigeLösung. Wie diese herausgerechnet werden, wird in [Rickert 2016] diskutiert.

4.2 Stationäre Wärmeleitung mit ROBIN-Randbedingung

Zunächst wird die stationäre Lösung des reinen Wärmeleitproblems ohne innereReibung untersucht. Dabei werden die Parameter der Erde verwendet, sowieℎ = 1 W m−2 K−1, 𝑇ext = 1 und 𝛾 = 1 gesetzt. Ein numerisches Ergebnis, zusammen mit der analytischen Lösung aus Gleichung (3.11), sind in Abbildung 4.2dargestellt. Insgesamt ergibt sich eine gute Übereinstimmung der Lösungen.

Für den gleichen Parametersatz wurde zudem eine Konvergenzstudie mit demFehlermaß

𝜀rel =

´Ω

(𝑇num − 𝑇ana

)d𝑉

´Ω𝑇ana d𝑉

(4.31)

durchgeführt. Hierbei bezeichnet 𝑇num die numerische Lösung und 𝑇ana die ana

0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1

0

0,05

0,1

0,15


Nor

mie

rte

Tem

pera

tur

(𝑇−

𝑇ex

t) /𝑇

1 numerischanalytisch

0

0,05

0,1

0,15

Abb. 4.2: Der blaue Graph stellt die numerische Lösung des reinen, entdimensionalisierten Wärmeleitproblems mit Robin-Randbedingung über den normierten Radius dar.Zum Vergleich ist die analytische Lösung als schwarz gestrichelte Linie eingezeichnet. Indas Bild ist ein Schnittbild der Kugel eingesetzt.



103 104 105

10−3

10−2,5

Freiheitsgrade

Rel

ativ

erFe

hler𝜀 r

el

Abb. 4.3: Der relative Fehler des stationären Wärmeleitproblems mit Robin-Randbedingung dargestellt über die Freiheitsgrade im doppelt logarithmischenMaßstab. Mit zunehmend feinerem Mesh nimmt der Fehler ab.

lytische Lösung aus Gleichung (3.11). Das Ergebnis ist in Abbildung 4.3 zu sehen.Die numerische Lösung konvergiert gegen die analytische. Jedoch verlangsamt sichdie Konvergenzordnung mit zunehmender Anzahl der Freiheitsgrade, wodurchdiese nichtlinear ist. Die Gründe hierfür sind unbekannt. Womöglich nimmt dierFehler beim Lösen des Gleichungssystems gegenüber dem Diskretisierungsfehlerder Finiten-Element-Analyse zu größeren Systemen hin zu. Dies wurde allerdingsnicht weiter überprüft.

4.3 Instationäre Wärmeleitung mitROBIN-Randbedingung

Als nächstes wird das instationäre Problem der reinen Wärmeleitung mit Robin-Randbedingung behandelt. Die Implementierung wird mit dem Testfallℎ = 𝛾 = 𝑇A = 1 überprüft, da für diesen eine analytische Lösung vorliegt.Zwei exemplarische Schnittbilder der numerischen Lösung sind in Abbildung 4.4zu sehen. Ein Vergleich zwischen analytischer und numerischer Lösung ist inAbbildung 4.5 dargestellt, wobei für die analytische Lösung 100 Summengliederberechnet wurden.

Die numerische Lösung zeigt eine gute Übereinstimmung mit der analytischenund beide Lösungen konvergieren gegen die stationäre Lösung. Die Ergebnisseder Implementierung können daher als annähernd korrekt betrachtet werden.

In Abbildung 4.6 werden numerische Ergebnisse der instationären Wärmeleitungmit Robin-Randbedingung für die Parameter des Erdmodells präsentiert. Auchhier konvergiert die gefundene Lösung gegen die stationäre Lösung, wobei diese

Instationäre Wärmeleitung mit ROBIN-Randbedingung

50

0

0,05

0,1

0,15

Abb. 4.4: Zwei Schnittbilder der Simulation des instationären, reinen Wärmeleitproblems mit Robin-Randbedingung. Die Oberfläche ist wiederum entsprechend der Temperatur eingefärbt, wobei rote Farben für höhere Temperaturen und blaue Farben fürniedrigere Temperaturen stehen. Das linke Schnittbild zeigt die Kugel zum Zeitpunkt𝑡 = 0,1, das Rechte zeigt in der gleichen Farbskala einen Schnitt zum Zeitpunkt 𝑡 = 0,4.

0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 10,3

0,4

0,5

0,6

0,7

0,8

0,9

1


Nor

mie

rte

Tem

pera

tur

(𝑇−

𝑇ex

t) /𝑇

1

𝑡 = 0𝑡 = 0,12𝑡 = 0,25𝑡 = 0,5𝑡 = 1𝑡 = 2stationär

Abb. 4.5: Dargestellt sind die numerische (gepunktet) und die analytische Lösung(durchgezogen) des instationären, reinen Wärmeleitproblems mit Robin-Randbedingungfür verschiedene Zeitpunkte für den Fall 𝛾 = ℎ = 𝑇A = 1. Die stationäre Lösung istschwarz gestrichelt eingezeichnet. Zum Zeitpunkt =2 liegen alle Lösungen übereinander.



0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1

0

0,05

0,1

0,15


Nor

mie

rte

Tem

pera

tur

(𝑇−

𝑇ex

t) /𝑇

1

𝑡 = 0𝑡 = 0,05𝑡 = 0,1𝑡 = 0,2𝑡 = 0,4stationär

Abb. 4.6: Dargestellt ist die instationäre, numerische Lösung des Ermodells mit Robin-Randbedingung zu verschiedenen Zeitpunkten. Die stationäre Lösung ist schwarzgestrichelt eingetragen. Der stationäre Zustand ist nach 𝑡 = 0,4 bereits weitestgehenderreicht.

bereits nach 𝑡 = 0,4 weitestgehend erreicht ist. Der stationäre Zustand wirddemnach deutlich schneller als im Testfall aus Abbildung 4.5 erreicht.

Instationäre Wärmeleitung mit ROBIN-Randbedingung

52

4.4 Stationäre Wärmeleitung mitSTEFAN-BOLTZMANN-Randbedingung

Ausgehend von der Implementierung des stationären Wärmeleitproblemsmit Robin-Randbedingung wird das stationäre Wärmeleitproblem mit Stefan-Boltzmann-Randbedingung bearbeitet. Es werden ausschließlich die Parameter des Erdmodells verwendet. Ein Vergleich der numerischen Lösung mit deranalytischen aus Gleichung (3.81) ist in Abbildung 4.7 zu sehen. Auch hier zeigtdie numerische Lösung eine gute Übereinstimmung mit der analytischen Lösung.

Wiederum wird mit dem Fehlermaß aus Gleichung (4.31) eine Konvergenzstudiedurchgeführt, welche in Abbildung 4.8 präsentiert wird. Mit zunehmend feineremMesh nimmt der Fehler ab, womit die numerische Lösung konvergiert. DieKonvergenzordnung verhält sich jedoch wieder nicht linear, wobei auch hier diegenauen Gründe unbekannt sind.

0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1

0

0,02

0,04

0,06

0,08

0,1

0,12

0,14

0,16

0,18


Nor

mie

rte

Tem

pera

tur

𝑇/𝑇

1

numerischanalytisch

0

0,05

0,1

0,15

Abb. 4.7: In blau ist die entdimensionalisierte Lösung des Wärmeleitproblems mitStefan-Boltzmann-Randbedingung über den Radius dargestellt. Zum Vergleich istdie analytische Lösung als schwarz gestrichelte Linie mit eingezeichnet. In das Bild istein Schnittbild der Kugel eingesetzt, wobei die Oberfläche entsprechend der Temperatureingefärbt ist.



103 104 105

10−3

10−2,5

Freiheitsgrade

Rel

ativ

erFe

hler𝜀 r

el

Abb. 4.8: Der relative Fehler des stationären Wärmeleitproblems mit Stefan-Boltzmann-Randbedingung dargestellt über die Anzahl der Freiheitsgrade imdoppelt logarithmischen Maßstab. Die Lösung konvergiert, die Kovergenzordnung istwiederum nicht linear.

4.5 Instationäre Wärmeleitung mit STEFAN-BOLTZMANN-Randbedingung

0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1

0

0,05

0,1

0,15


Nor

mie

rte

Tem

pera

tur

𝑇/𝑇

1

𝑡 = 0𝑡 = 0,05𝑡 = 0,1𝑡 = 0,2𝑡 = 0,4stationär

Abb. 4.9: Dargestellt ist die numerische Lösung des instationären Wärmeleitproblemsmit Stefan-Boltzmann-Randbedingung für die Parameter aus Tabelle 2.1, zu verschiedenen Zeitpunkten. Zum Vergleich ist die stationäre Lösung als gestrichtelte, schwarzeLinie mit eingezeichnet. Nach 𝑡 = 0,4 hat das Temperaturfeld weitestgehend seinenstationären Zustand erreicht.

Das instationäre Wärmeleitproblem mit Stefan-Boltzmann-Randbedingung

Instationäre Wärmeleitung mit STEFAN-BOLTZMANN-Randbedingung

54

wird wiederum nur für das Erdmodell untersucht. Zudem liegt keine analytischeLösung vor, sodass kein Vergleich möglich ist, der die Korrektheit der gefundenenLösung bestätigt. Jedoch wurde die Konvergenz gegen die stationäre Lösungüberprüft. In Abbildung 4.9 ist das Ergebnis der Simulation dargestellt. Diestationäre Lösung ist wieder als gestrichelte, schwarze Linie eingezeichnet. Derstationäre Zustand wird nach ungefähr 𝑡 = 0,4 erreicht.

4.6 Stationäres, thermoelastisches Problem mit STEFAN-BOLTZMANN-Randbedingung

Das stationäre, thermoelastische Problem kann ohne weiteren Aufwand gelöstwerden, da die innere Reibung verschwindet. Es wird zunächst das stationäreTemperaturfeld berechnet und mit diesem anschließend das stationäre Verschiebungsfeld gesucht. Für den Fall 𝛾 = 1 liegt eine exakte Lösung aus Gleichung (3.4)vor.

Ein Vergleich der numerischen und analytischen Ergebnisse für die Parameteraus Tabelle 2.1 sind in Abbildung 4.10 und Abbildung 4.11 dargestellt.

Das Temperaturfeld zeigt erneut eine gute Übereinstimmung. Auch die numerische Lösung für die radiale Komponente des Verschiebungsfeldes stimmt mitder analytischen Lösung einigermaßen überein. Anders als bei den Lösungen derstationären Wärmeleitprobleme ist jedoch ein deutlicher Fehler sichtbar.

0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 10

0,05

0,1

0,15


Nor

mie

rte

Tem

pera

tur

𝑇/𝑇

1 numerischanalytisch

Abb. 4.10: Der blaue Graph stellt die numerische Lösung des Temperaturfeldes desthermoelastischen Problems mit Stefan-Boltzmann-Randbedingung dar. Sie ist identisch mit der Lösung aus Abbildung 4.7. Die analytische Lösung ist als gestrichelte,schwarze Linie eingezeichnet.



0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1

0

0,1

0,2

0,3

0,4


Nor

mie

rte

Vers

chie

bung

𝑢/𝑢

0

numerischanalytisch

0

0,1

0,2

0,3

Abb. 4.11: Der blaue Graph stellt die numerische Lösung für die radiale Komponente des Verschiebungsfeldes des stationären thermoelastischen Problems mit Stefan-Boltzmann-Randbedingung dar. Die analytische Lösung ist als gestrichelte, schwarze Linie eingezeichnet. Am unteren, rechten Bildrand ist ein Schnitt durch die Erdegezeigt. Zum Rand hin sind zwischen numerischer und analytischer Lösung deutlicheUnterschiede sichtbar, insgesamt zeigen beide Kurven aber das gleiche Verhalten.

4.7 Instationäres, thermoelastisches Problem mitSTEFAN-BOLTZMANN-Randbedingung

Im Gegensatz zum stationären thermoelastischen Problem ist die Lösung desinstationären thermoelastischen Problems deutlich aufwendiger, da innere Reibung nun eine Rolle spielt. Diese wird zumeist vernachlässigt, weshalb hier derEinfluss der inneren Reibung auf die thermische Entwicklung im Vordergrundsteht. Zum simultanen Ermitteln des Temperatur und des Verschiebungsfeldeswurde der Algorithmus aus Abbildung 4.1 eingesetzt.

Ein numerisches Ergebnis des Temperaturfeldes für das Erdmodell, zusammenmit der analytischen, stationären Lösung ist in Abbildung 4.12 zu verschiedenenZeiten dargestellt. Der stationäre Zustand wird erst nach 𝑡 = 0,8 annähernderreicht. Gegenüber der reinen Wärmeleitung hat die innere Reibung zu einemverzögerten Aufheizen des Körpers geführt.

In Abbildung 4.13 ist die numerische Lösung des normierten Verschiebungsfeldes des Erdmodells zu verschiedenen Zeiten dargestellt. Die gefundene Lösungkonvergiert gegen das stationäre Ergebnis und erreicht diese nach 𝑡 = 0,8 annähernd. Die maximale Verschiebung ist ungefähr 0,4. Angesichts der großen

Instationäres, thermoelastisches Problem mit STEFAN-BOLTZMANN-Randbedingung

56

0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1

0

0,05

0,1

0,15


Nor

mie

rte

Tem

pera

tur

𝑇/𝑇

1 𝑡 = 0𝑡 = 0,1𝑡 = 0,2𝑡 = 0,8stationär

Abb. 4.12: Das normierte Temperaturfeld des instationären, thermoelastischen Problems für das Erdmodell zu verschiedenen Zeiten. Zum Vergleich ist die Lösung desinstationären Problems der reinen Wärmeleitung mit Stefan-Boltzmann-Randbedingung gepunktet eingezeichnet. Die stationäre Lösung ist schwarz gestrichelt dargestellt.

charakteristischen Verschiebung 𝑢0 aus Tabelle 2.2 entspricht dies einer Ausdehnung der Erde um beinahe 0,4𝑅. Dies ist unrealistisch. Der Grund hierfürist unter anderem die große Temperaturdifferenz über den Körper. Über denRadius steigt die Temperatur von beinahe 0 auf 0,16 an, was zusammen mit dercharakteristischen Temperatur von ungefähr 105 K mehreren 10 000 K Temperaturunterschied entspricht. Da sich die Materialparameter weder mit Druck, nochmit der Temperatur ändern führt dies zu einer großen thermischen Ausdehnung.

Zuletzt wird noch der Term der inneren Reibung ausgewertet, um eine Vorstellungdavon zu bekommen, wie lange dieser aktiv ist und welche Werte er annimmt.Die negative Änderung der inneren Energie aufgrund innerer Reibung ist nachGleichung (2.70) gegeben durch

𝑤 := 𝛼 (3𝜆+ 2𝜇) (𝑇 − 𝑇R)(𝜕𝑢

𝜕𝑡· ∇). (4.32)

Diese Leistungsdichte wird mit dem Verschiebungsfeld ausgetauscht und in derelastischen Deformation gespeichert. Setzt man die entdimensionalisierten Größenein, erhält man

�� =(𝑇 − 𝑇R

)(𝜕��

𝜕𝑡· ∇). (4.33)

An dieser Stelle wird die innere Reibung über das gesamte Volumen Ω integriert, wodurch sich die gesamte, zwischen Temperaturfeld und Deformation



0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1

0

0,1

0,2

0,3

0,4


Nor

mie

rte

Vers

chie

bung

𝑢/𝑢

0

𝑡 = 0𝑡 = 0,1𝑡 = 0,2𝑡 = 0,4𝑡 = 0,8stationär

Abb. 4.13: Die normierte Verschiebung des instationären, thermoelastischen Problemszu verschiedenen Zeiten. Die stationäre Lösung (schwarz gestrichelt) wird nach 𝑡 = 0,8annähernd erreicht.

ausgetauschte, normierte Leistung ergibt:

�� :=ˆ

Ω

(𝑇 − 𝑇R

)(𝜕��

𝜕𝑡· ∇)

d𝑉 (4.34)

Die integrierte, innere Reibung ist über die Zeit dargestellt in Abbildung 4.14 zufinden. Es ist zu erkennen, dass sie während des gesamten instationären Prozesseseine Rolle spielt.

Instationäres, thermoelastisches Problem mit STEFAN-BOLTZMANN-Randbedingung

58

0 0,2 0,4 0,6 0,8 1 1,2 1,4 1,6 1,8 20

0,2

0,4

0,6

0,8

Normierte Zeit 𝑡/𝜏

Inne

reR

eibu

ng 𝑊

Abb. 4.14: Die integrierte, entdimensionalisierte, innere Reibung �� dargestellt überdie normierte Zeit 𝑡. Sie nimmt bei 𝑡 = 0,025 ihr Maximum von 0,85 an. Danach fällt sielangsam ab und verschwindet nach ungefähr 𝑡 = 1,5.

4.8 Erdmodell mit Sonneneinstrahlung

Zum Abschluss wird eine numerische Lösung des Erdmodells für das reine Wärmeleitproblem mit Stefan-Boltzmann-Randbedingung und Sonneneinstrahlungpräsentiert.

Das radiale, dimensionsbehaftete Temperaturfeld ist in Abbildung 4.15 dargestellt.Die Temperatur steigt zum Kern hin auf rund 90 000 K an, was weit vomLiteraturwert von 6 000 K entfernt ist.

In Abbildung 4.16 ist das dimensionsbehaftete Temperaturfeld über den Polarwinkel vom Nordpol zum Äquator gezeichnet. Die Temperatur bei 𝜗 = π/2 liegtbei 296,71 K oder auch 23,56 ∘C. Die Temperatur am Nordpol beträgt hingegen88,43 K oder auch −184,72 ∘C. Die mittlere Temperatur der Erde beträgt indiesem Modell 273,65 K, was 0,5 ∘C entspricht. Die Temperatur am Äquator unddie mittlere Oberflächentemperatur stimmen gut mit Beobachtungen überein. DieTemperatur am Pol erscheint etwas niedrig. Allerdings liegt die tiefste, gemesseneTemperatur auf der Erde bei −93,2 ∘C. Zudem liegt im hier vorgestellten Modellder Pol ständig in Dunkelheit, er erhält keine Wärme durch die Sonne. Bewirktwird dies durch die fehlende Achsneigung der Erde gegen ihre Bahnebene um dieSonne. Dies führt dazu, dass die Pole je halbjährlich von der Sonne beschienenwerden. Auch gibt es im vorgestellten Modell keine Atmosphäre, die ausgleichendauf Temperaturdifferenzen wirkt. Daher ist zu erwarten, dass auch die Poltemperatur näher an gemessene Werte rückt, falls diese Mechanismen berücksichtigtwerden.



0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1

20 000

40 000

60 000

80 000

0,997 0,9985 1

200

400

600

800

normierter Radius 𝑟/𝑅

Tem

pera

tur𝑇

inK

𝜗 ∈ [0,π] beliebig𝜗 = 0𝜗 = π/2

Abb. 4.15: Der radiale, dimensionsbehaftete Temperaturverlauf über den normiertenRadius. Auf dieser Skala lässt sich die Variation über den Polarwinkel nicht auflösen. Entsprechend ist ein stark vergrößerter Ausschnit eingefügt, der die Unterschiedeverdeutlicht.

0 0,2 0,4 0,6 0,8 1 1,2 1,4100

150

200

250

300

Polarwinkel 𝜗

Tem

pera

tur𝑇

inK

𝑟/𝑅 = 1Mittelwert

100

150

200

250

Abb. 4.16: Der dimensionsbehaftete Temperaturverlauf für 𝑟/𝑅 = 1, startend amNordpol (𝜗 = 0) und endend am Äquator (𝜗 = π/2). Der Mittelwert ist als schwarze,gestrichelte Linie eingezeichnet. Die Temperatur am Pol beträgt 88 K, am Äquatorhingegen 296,71 K.

Erdmodell mit Sonneneinstrahlung

60

5 Das Alter der Erde

Über die Kenntnis des thermischen Zustands der Erde lässt sich ihr Alterabschätzen. Dazu werden folgende Annahmen getroffen:

1. Homogenes, isotropes Material

2. Wärmeverlust auf der Oberfläche nur durch Strahlung

3. Energietransport im Inneren durch reine Wärmeleitung

4. Die Erde hat ihren stationären Zustand erreicht

Die Punkte 1 bis 3 entsprechen den in der Modellierung der Erde getroffenenAnnahmen, Punkt 4 wird im nachstehenden Abschnitt erläutert.

Die Erde empfängt von der Sonne die Leistung

πˆ

𝜗=0

π/2ˆ

𝜙=−π/2

𝑅2𝜀𝑆 sin(𝜗)2 cos(𝜙) d𝜙d𝜗 ≈ 5,2 × 1016 W . (5.1)

Die Wärmeproduktion innerhalb der Erde beträgt nach Kapitel 2.5 ungefähr3,4 × 1013 W. Das ist im Verhältnis zur eingestrahlten Wärme sehr klein: DieTemperatur der Erdoberfläche wird durch die Sonneneinstrahlung dominiert.Befände sich die Erde nicht im thermischen Gleichgewicht, tauschte sie einenWärmestrom mit der Sonne aus. Speziell erhielte die Erde einen Wärmestrom vonder Sonne, da diese sehr viel wärmer als die Erdoberfläche ist. Dies führt zu einemAnstieg der Oberflächentemperatur. Jedoch gibt es geologische Hinweise darauf,dass sich die Temperatur der Erdoberfläche während der letzten Milliarde Jahrenicht wesentlich geändert hat. Entsprechend befindet sich die Erdoberfläche imthermischen Gleichgewicht, die Erde tauscht mit ihrer Umgebung keine Wärmeaus und befindet sich insgesamt in einem quasi stationären1 Zustand. Daher wirdhier angenommen, dass die Erde ihren stationären Zustand erreicht hat.

Die Simulationen aus Kapitel 4 zeigen, dass im hier vorgestellten Modell derstationäre Zustand nach einer Zeit 0,7𝜏 erreicht wird. Unter den vorgestellten Annahmen kann die Erde demnach nicht jünger als 0,7𝜏 sein, womit aus

1Über sehr große Zeiträume nimmt die Wärmeproduktion durch radioaktiven Zerfall imInneren der Erde ab. Gleichzeitig nimmt die Leuchtkraft der Sonne um ungefähr 10 %pro einer Milliarde Jahre zu, sodass kein stationärer Zustand existiert. Auf einige hundertMillionen Jahre beschränkt ist ein quasi stationärer Zustand dennoch eine gute Näherung.

Kapitel 5. Das Alter der Erde


Tabelle 2.2

0,7𝜏 ≈ 4,51 × 1019 s ≈ 1,43 × 1012 a (5.2)

folgt. Das entspricht ungefähr dem 104-fachen Alter des Universiums von 13,7 ×109 a nach aktuellem Wissensstand. Offensichtlich ist das ein zu großes Alter.

In [Thomson 1862] schätzte Lord Kelvin bereits 1862 das Alter der Erde überihren Wärmeverlust auf 20 bis 400 Millionen Jahre ab. Dazu stellte er die gleichenAnnahmen 1 bis 3 auf, nahm allerdings einen instationären Zustand der Erde an.Bereits damals war bekannt, dass die Temperatur in tiefen Minen um ungefähr3 K pro 100 m zunimmt. Seine Hypothese war, dass diese Wärme aus der Zeit derEntstehung der Erde stammte und sich diese noch immer abkühlt. RadioaktiverZerfall war zu dem Zeitpunkt noch unbekannt.

Unter diesen Annahmen nutzte Kelvin die kurz zuvor veröffentlichten Ergebnissevon Fourier, welcher das Wärmeleitproblem eines halbunendlichen Mediums mitkonstanter Temperatur auf der Oberfläche und im Unendlichen gelöst hatte. DieTemperatur der Erdoberfläche legte Kelvin mit 𝑇O = 0 ∘C fest. Die Temperaturin unendlicher Entfernung setzte er auf die geschätzte Temperatur der Erde kurznach ihrer Entstehung von 𝑇∞ = 4000 K. Das Problem lautet damit

𝜕𝑇

𝜕𝑡= 𝑎

𝜕2𝑇

𝜕𝑥2 , (𝑥, 𝑡) ∈ [0 ,∞)2 , (5.3)

𝑇 (𝑡, 0) = 𝑇O , 𝑡 > 0 , (5.4)lim

𝑥→∞𝑇 (𝑡, 𝑥) = 𝑇∞ , 𝑡 > 0 , (5.5)

𝑇 (0, 𝑥) = 𝑇∞ . (5.6)

Die Lösung dieses Problems erhält man am leichtesten über die Laplace-Transformation, wie zum Beispiel in [Baehr u. Stephan 1994] nachzulesenist. Für das Temperaturfeld erhält man

𝑇 (𝑡, 𝑥) = 𝑇∞ erf(

𝑥

2√𝑎𝑡

). (5.7)

Der Temperaturgradient an der Erdoberfläche ist aus Messungen schon zuKelvins Zeiten bekannt gewesen. Es gilt

𝑔(𝑥, 𝑡) := 𝜕𝑇

𝜕𝑥(𝑥, 𝑡) = 𝑇∞√

π𝑎𝑡exp

(𝑥2

4𝑎𝑡

)(5.8)

beziehungsweise für 𝑥 = 0:

𝑔(0 m, 𝑡) = 𝑇∞√π𝑎𝑡

. (5.9)

62

Dies umgestellt nach der Zeit führt zu

𝑡(𝑔) = 𝑇 2∞

π𝑔2𝑎(5.10)

was mit den Werten aus Tabelle 2.1 das Alter der Erde 𝑡E ergibt:

𝑡E ≈ 8,98 × 1015 ≈ 2,85 × 108 a . (5.11)

Dieses Alter ist viel geringer als das heutzutage akzeptierte Alter der Erde von4 × 109 a und auch Kelvin war sich der Unsicherheiten in den Parameternbewusst, weshalb er einen großen Bereich von 20 bis 400 Millionen Jahre angab.Dieses Ergebnis stieß auf viel Gegenwind von Geologen, die das Alter der Erdeaufgrund von Sedimentationsprozessen auf weit über eine Milliarde Jahre schätzten. Gewicht erlangte Kelvins Argument durch eine unabhängige Schätzungdes Alters der Sonne durch Helmholtz. Dieser schlug gravitative Energie alsEnergiequelle der Sonne vor. Dazu betrachtete er eine Gaswolke mit der Masseder Sonne, die den gesamten Raum homogen ausfüllt. Aufgrund kleiner Inhomogenitäten stürzt diese Gaswolke durch Gravitation in sich zusammen. Dabeiwird potentielle Energie frei, die das Gas aufheizt und es so zum Leuchten bringt.Helmholtz berechnete über den aktuellen Durchmesser der Sonne, dass dieserEnergievorrat für nicht länger als 20 Millionen Jahre gereicht haben könnte,was nah an der Schätzung durch Kelvin lag. Weiterhin argumentierte man,dass die Erde nur so lang Leben beherbergen kann, wie sie Wärme durch dieSonne erhält. Diese unabhängigen, scheinbar übereinstimmenden Entdeckungenverliehen Kelvins Argument derart Gewicht, dass man die Schätzungen derGeologen als falsch erachtete. Auch die Entdeckung des radioaktiven ZerfallsEnde des 19. Jahrhunderts konnte diesen Streit nicht auflösen, denn wie heutebekannt ist, erhält die Sonne ihre Energie aus thermonuklearen Fusionsprozessen.

Beide vorgestellten Schätzungen liegen daneben. Hauptursache dafür ist dieAnnahme eines homogenen Körpers, der im Inneren Wärme nur durch reineWärmeleitung transportiert. Kelvin kommentierte sogar, dass andere Transportmechanismen seine Rechnung obsolet werden ließen. Die Erde besteht aus einerdünnen, thermisch isolierenden Kruste, unterhalb derer ein teilweise flüssigerMantel liegt, der wiederum einen festen Eisenkern umschließt. Mantel und Kernleiten Wärme sehr viel besser als die Kruste, zudem ist im Mantel ein starkerkonvektiver Transport aktiv. Die Temperatur im Mantel und im Erdkern beträgtmehrere tausend Kelvin, wodurch Strahlungstransport innerhalb des Mediumseine weitere Rolle beim Transport innerer Energie spielt, siehe auch [Ayers 1970].

Kapitel 5. Das Alter der Erde


6 Fazit und Ausblick

Dieses Projekt hatte die Modellierung und Berechnung der thermischen Entwicklung der Erde als Ziel. Als Ausgangspunkt wurden die globale Bilanz derMasse, des Impulses und der Energie gewählt, welche in ihre lokale Form überführt wurden. Bei der Modellierung der konstitutiven Gleichungen wurde derzweite Hauptsatz der Thermodynamik in lokaler Form ausgewertet und so diethermodynamische Verträglichkeit der Materialfunktionen gesichert. Dabei ergabsich eine Kopplung zwischen Impuls und innerer Energie, der entsprechendeKopplungsterm wird innere Reibung genannt. Diese ist nur bei instationärenProzessen aktiv, wird aber in gängiger Literatur praktisch immer vernachlässigt. Die Untersuchung dieser Kopplung wurde zu einem sekundären Ziel diesesProjektes.Die Randbedingung der Temperatur wurde durch das Stefan-Boltzmann-Gesetzmodelliert, welche einen Strahlungsaustausch mit einer ideal strahlenden Umgebung darstellt. Diese Umgebung wurde im ersten Modell als die kosmischeHintergrundstrahlung gesetzt, in einem erweiterten Modell wurde Sonneneinstrahlung berücksichtigt. Das so gestellte Problem ist nichtlinear, wodurchanalytische Lösungen schwer zugänglich sind und die Auswahl an Literaturbeschränkt ist. Daher wurde ausgehend von einem bekannten Problem derWärmeleitung mit Robin-Randbedingung die Wärmeleitung mit Strahlungsrandbedingung erarbeitet. Die Robin-Randbedingung diente dabei als robusterTestfall für die eingesetzten analytischen Methoden und den implementiertenCode.Die Feldgleichungen wurden für den stationären, radialsymmetrischen Fall undteilweise auch für den instationären Fall analytisch gelöst. Anschließend wurdendie Feldgleichungen in ihre schwache Form überführt und mit Hilfe einer Finiten-Elemente-Analyse, die ein Galerkin-Verfahren mit Lagrange-Polynomenerster Ordnung umsetzte, gelöst. Der Löser und der Mesher wurden durch diePython Bibliothek FEniCS bereitgestellt. Numerische Ergebnisse wurden, soweitvorhanden, mit analytischen Ergebnissen verglichen, wobei sich größtenteils einegute Übereinstimmung zeigte.Die Oberflächentemperatur der Erde konnte gut abgebildet werden. Die Strukturdes Temperaturfelds im Inneren der Erde weicht jedoch deutlich von Literaturwerten ab. Dies ist auf die Modellierung der Erde als homogenen Festkörperzurückzuführen. Den Abschluss bildet der Vergleich mit Kelvins Argument ausdem 19. Jahrhundert, wobei eine eigene Schätzung aufgestellt wurde.

64

Ausgehend von den hier gewonnenen Ergebnisse eröffnen sich viele Möglichkeiten der Erweiterung des Modells. Ein verbessertes Modell der Erde kann dieSchalenstruktur berücksichtigen und so ein besseres Abbild der thermischen Verhältnisse der Erde im Inneren schaffen. Der Transportmechanismus für Wärmeim Inneren der Erde kann um Konvektion und Strahlung erweitert werden. Auchdie Untersuchung von halbtransparenten Materialien kann das Verständnis umden Strahlungstransport in Materie weiter ausbauen. Die Kopplung zwischenVerschiebungsfeld und innerer Energie kann genutzt werden um den Einflussvon Gezeitenkräften auf das Innere der Erde zu untersuchen. So wird vermutet,dass dieser Mechanismus für den starken Vulkanismus auf dem Jupitermond Ioverantwortlich ist. Auch Europa unterliegt starken Gezeitenkräften, die genugWärme produzieren sollen, dass ein ausgedehnter flüssiger Ozean unter seinemEispanzer existiert.

Kapitel 6. Fazit und Ausblick

Wärmeleitung mit Strahlungsrandbedingung III

Literaturverzeichnis

[Ayers 1970] Ayers, D L.: Transient cooling of a sphere in space. In: J. HeatTransfer 2 (1970), S. 180–182

[Baehr u. Stephan 1994] Baehr, H D. ; Stephan, K: Wärme- und Stoffübertragung. Bd. 8. Springer, 1994

[Bertram 2013] Bertram, A: Festkörpermechanik. Otto-von-Guericke UniversitätMagdeburg, 2013

[Dreyer 2013] Dreyer, W: Kontinuumsphysik thermoelastischer Probleme. TUBerlin, 2013

[Fitzpatrick 2006] Fitzpatrick, R: Thermodynamics and Statistical Mechanics:An intermediate level course. Lulu Enterprises. Inc, 2006

[Hahn u. Özişik 2012] Hahn, D W. ; Özişik, M N.: Heat Conduction. Bd. 3.Wiley Online Library, 2012

[Howell u. a. 2010] Howell, J R. ; Menguc, M P. ; Siegel, R: Thermalradiation heat transfer. Bd. 5. CRC press, 2010

[Ji u. a. 2010] Ji, S ; Sun, S ; Wang, Q ; Marcotte, D: Lamé parametersof common rocks in the Earth’s crust and upper mantle. In: Journal ofGeophysical Research: Solid Earth 115 (2010), Nr. B6

[Logg u. a. 2012] Logg, A ; Mardal, K A. ; Wells, G: Automated solution ofdifferential equations by the finite element method: The FEniCS book. Bd. 84.Springer Science & Business Media, 2012

[McCluney 2014] McCluney, W R.: Introduction to radiometry and photometry.Artech House, 2014

[Rickert 2016] Rickert, W: Kontinuumsphysikalische Untersuchung von Magnetostriktion mittels Methoden der mathematischen Physik sowie gekoppelterFinite-Elemente-und Rand-Elemente-Analyse. TU-Berlin, 2016

[Robertson 1988] Robertson, E C.: Thermal properties of rocks / US GeologicalSurvey,. 1988. – Forschungsbericht

[Thomson 1862] Thomson, W: XV.—On the Secular Cooling of the Earth. In:Transactions of the Royal Society of Edinburgh 23 (1862), Nr. 01, S. 157–169

IV

[Truesdell 2013] Truesdell, C: Linear Theories of Elasticity and Thermoelasticity: Linear and Nonlinear Theories of Rods, Plates, and Shells. Bd. 2.Springer, 2013

[Turcotte u. Schubert 2014] Turcotte, D L. ; Schubert, G: Geodynamics.Bd. 2. Cambridge University Press, 2014

[Williams 2004a] Williams, D R.: Earth fact sheet. In: Structural geology of theEarth’s interior: Proc. Natl. Acad. Sci. NASA (17 Nov 2010) 76 (2004), Nr. 9

[Williams 2004b] Williams, D R.: Sun fact sheet. https://nssdc.gsfc.nasa.gov/planetary/factsheet/sunfact.html, 2004

Literaturverzeichnis

https://nssdc.gsfc.nasa.gov/planetary/factsheet/sunfact.html

https://nssdc.gsfc.nasa.gov/planetary/factsheet/sunfact.html

projekt simulationstools und ihre anwendung · 2 2 modellbildung dieses kapitel befasst sich mit...

Documents