Vorbemerkung:

Dieses Mathematica- Notebook behandelt nicht die im CA-Rundbrief diskutierte Erwärmung einer Zugbremsscheibe, sondern ein verwandtes Problem, nämlich die Wärmeleitung in einem (dünnen) Stab. Obwohl dieses (Lehrbuch-) Beispiel weniger komplex ist, zeigt es bereits alle Eigenschaften des vollen Problems. Der unten dargestellte Fall (feste Vorgabe der Temperatur an den Stab-Enden) kann auf kompliziertere Randbedingungen erweitert werden.  Die Darstellung erfolgt in der Standard -Mathematica- Notation (nur diese ist programmiertechnisch eindeutig). Alle Erläuterungen sind als Kommentare eingefügt.

Die Umsetzung in die PDF-Datei erfolgte aus dem Mathematica-Menu. Leider ist die Formatierung an einigen Stellen suboptimal.

D. Hackenbracht

Viele Vorgänge werden durch partielle Differentialgleichungen (PDE) beschrieben. Die gesuchte Funktion hängt dabei von mehreren unabhängigen Variablen ab. Die wichtigsten linearen Vertreter dieser PDE sind:

- die Wärmeleitgleichung (allgemeiner: die Diffusionsgleichung)
- die Wellengleichung
- die Potentialgleichung (auch Laplace-Gl. genannt) .

Neben Anfangsbedingungen (Vorgabe der Funktion zur Zeit  t = 0 ) müssen auch Randbedingungen gestellt werden (Vorgabe der Funktion oder ihrer räumlichen Ableitung am Rande des untersuchten Gebietes). Die exakte Behandlung (im Sinne einer einfachen Lösungsformel) ist meistens nicht möglich. Oft benötigt man (unendliche) Reihen, um die Lösung darzustellen. Aber auch dies geht nur in relativ einfachen Fällen.

Im Folgenden wird die Wärmeleitung in einem Stab untersucht. Der Einfachheit halber gehen wir davon aus, dass die Temperatur nur in der Stab- Längsrichtung (die hier die z-Richtung sein soll) variiert. Die Temperatur   T   ist dann eine Funktion von   t  und  z . Man spricht dann auch von der [räumlich] eindimensionalen Wärmeleitgleichung. Die Herleitung der Wärmeleitgleichung findet man in einem Physikbuch, z.B. Alonso/Finn, Fundamental University Physics, Vol. 2 . Außerdem soll der Stab  - außer an seinen Enden -  thermisch isoliert sein.

A) Die Vorbereitung

Remove["Global`*"]

(* Die Materialdaten sind die für reines Kupfer  [Quelle: Deutsches Kupferinstitut] ; 
   eigentlich hängen diese Daten auch von der Temperatur ab, 
   aber dass soll im Folgenden vernachlässigt werden; 

   λ ist die Wärmeleitfähigkeit in W/( K m) , cp die spezifische Wärme in J /(kg K) ,   
   ρ die Dichte in kg/(m^3)  *)

Materialdaten = {λ->395 ,cp->380 ,ρ->8960 }

(* die relevanten geometrischen Daten, d.h. hier linker und rechter Rand des Stabes, Einheit Meter *)

z0=0;

zn=1;

(* die Bestandteile der Wärmeleitgleichung : *)

(* die partielle Ableitung der Temperatur T  nach der Zeit  t : *)

Tpunkt[z_,t_]=D[T[z,t],t]

(* die 2. partielle Ableitung der Temperatur T  nach der Raumkoordinate z :  *)

Tablzz[z_,t_]=D[T[z,t],{z,2}]

(* die Wärmeleitgleichung... *)

dgl= Tpunkt[z,t]== ((λ Tablzz[z,t])/(cp ρ))

(* ...mit den Materialdaten: *)

dgl=dgl/.Materialdaten

(* die zusätzlichen Bedingungen: *)

(* Anfangstemperatur in °C *)

T0=20;

(* 1. 
 als Ausgangslage/Anfangsbedingung zur Zeit t=0
 nehmen wir eine homogene Temperaturverteilung im Stab an: *)

anfangsbedingung=T[z,0]==T0;

(* 2a.
am linken Rand wird die Temperatur zur Zeit t=0
 plötzlich auf 100 °C angehoben und auf diesem Wert festgehalten: *)

randbedingungLinks=T[0,t]== 100;

(* 2b.
 am rechten Rand wird die Temperatur auf dem Wert T0 festgehalten:*)

randbedingungRechts=T[zn,t]== T0;

(* alle Gleichungen zusammen: *)

gl= {dgl,anfangsbedingung, randbedingungLinks ,  randbedingungRechts}

B) Die Berechnung der Lösung

(* alle Gl. werden dem Mathematica- Lösungsbefehl übergeben; 
   DSolve löst Differentialgleichungen exakt! *)

ls = DSolve[ gl,T[z,t],{z,t}]

(* aus den o.g. Gründen kann Mathematica dieses Problen nicht exakt lösen; 
   also ist eine numerische Behandlung erforderlich; 
   NDSolve ermittelt numerisch eine approxiamtive Lösung der Differentialgleichung *)

(* Anfangs- und Endezeit des Simulationsintervalls, in Sekunden: *)

t0=0;

tn=1000;

ls = NDSolve[ gl,T[z,t],{z,z0,zn},{t, t0,tn}]

(* die Aufbereitung des Ergebnisses : *)

lsf=ls//Flatten;

TL[z_,t_]= T[z,t]/.lsf

(* die Gesamtsicht: *)

Plot3D[TL[z,t],{z,z0,zn},{t, t0,tn},AxesLabel->{"Axiale z- Richtung [ m]\t","\tZeit [s]","Temperatur [°C]\t"}]

(* Es ist gar nichts passiert !! 

Die obige Warnung muss man also sehr ernst nehmen; 
 " Boundary and initial conditions are inconsistent " trifft ja auch zu:
 über die Anfangsbedingung geben wir am linken Rand T[0,0] == T0 vor,  
   über die Randbedingung aber T[0,0]== 100 ! 
 Dieser plötzliche Anstieg der Temperatur ist ja auch unphysikalisch. Man muss   
 ihn also anders modellieren, also die Temperatur allmählich hochfahren *)

C) Korrektur des Modells und erneute Berechnung der Lösung

(* also: 
 am linken Rand wird die Temperatur zur Zeit t=0
 allmählich auf 100 °C angehoben und auf diesem Wert festgehalten; 
   dies geschieht hier mit einer geeigneten "Anschalt"-Funktion f(t) : *)

f[t_]:= T0+80*(1-Exp[-t])

Plot[f[t],{t,t0,t0+10},PlotRange->All,AxesLabel->{" t \t"," f \t"}]

(* die e-Funktion klingt so schnell ab, 
   dass die Temperatur nach ca 5s den gewünschten konstanten Wert 100°C hat; 

   korrigierte Randbedingung: *)

randbedingungLinks=T[0,t]==f[t]

(* neue Rechnung: *)

Clear[T]

gl= {dgl,anfangsbedingung, randbedingungLinks ,  randbedingungRechts}

ls = NDSolve[ gl,T[z,t],{z,z0,zn},{t, t0,tn}]

(* die Aufbereitung des Ergebnisses : *)

lsf=ls//Flatten;

TL[z_,t_]= T[z,t]/.lsf

D) Graphische Darstellung der Lösung

(* die Gesamtsicht: *)

Plot3D[TL[z,t],{z,z0,zn},{t, t0,tn},PlotRange->{0,120},AxesLabel->{"Axiale z- Richtung [ m]\t","\tZeit [s]","Temperatur [°C]\t"}]

(* Momentaufnahme: das Temperaturprofil längs des Stabes am Ende der Simulation: *)

Plot[TL[z,tn],{z,z0,zn},PlotRange->{0,120},AxesLabel->{" z \t"," T \t"},PlotLabel-> "räumlicher Temperaturverlauf längs des Stabes am Ende der Simulation"]

(* der zeitliche Temperaturverlauf an einer festen Stelle, hier z= 0.1: *)

Plot[TL[0.1,t],{t,t0,tn},PlotRange->{0,120},AxesLabel->{" z \t"," T \t"},PlotLabel-> "zeitlicher Temperaturverlauf in der Nähe des linken Randes"]

(* daraus ein Detail: *)

Plot[TL[0.1,t],{t,t0,tn/100},AxesLabel->{" z \t"," T \t"},PlotLabel-> "anfänglicher Temperaturverlauf in der Nähe des linken Randes"]

Die unschöne "Delle" im Temperaturverlauf ist, wie man so sagt, ein "numerisches Artefakt" ; 
 die numerische Stabilität der Lösung muss - in jedem Fall - durch eine genauere Rechnung geprüft werden

E) Kontrolle: eine genauere Rechnung

ls = NDSolve[ gl,T[z,t],{z,z0,zn},{t, t0,tn},
                           MaxStepSize->{zn/100,  tn/1000}, 
                           MaxSteps->Infinity, AccuracyGoal->11]

lsf=ls//Flatten;

TLneu[z_,t_]= T[z,t]/.lsf

(* der direkte Vergleich beider Lösungen: *)

Plot3D[TLneu[z,t]-TL[z,t],{z,z0,zn},{t, t0,tn},PlotRange->{-10,10},AxesLabel->{"Axiale z- Richtung [ m]\t","\tZeit [s]","Abweichung [°C]\t"}]

Das sieht gut aus;  insgesamt ist die Lösung also numerisch stabil

(* das kritische Detail: *)

Plot[TLneu[0.1,t],{t,t0,tn/100},AxesLabel->{" z \t"," T \t"},PlotLabel-> "anfänglicher Temperaturverlauf in der Nähe des linken Randes"]

Die Delle ist weg - oder ?

(* zur Kontrolle einige Zahlenwerte: *)

Table[TLneu[0.1,t],{t,0,4,0.1}]

{20.`,20.000000000000966`,20.000000000025647`,19.999999999517808`,19.99999999629978`,19.999999989834823`,19.999999993613375`,20.000000045193175`,20.000000188381946`,20.000000456863607`,20.000000844994958`,20.000001206084377`,20.00000116792606`,20.000000150884556`,19.99999735913098`,19.999992004817617`,19.999983355475646`,19.999970974433463`,19.999955138404804`,19.999937193092485`,19.999919837129394`,19.999908839314052`,19.99991053646658`,19.999933968371465`,19.999993722600973`,20.000102933422667`,20.000283349862844`,20.00055473438746`,20.000944419842586`,20.001481556486766`,20.002198099940426`,20.003134090854655`,20.00432447360984`,20.005810760257454`,20.007640509875024`,20.0098566230452`,20.012507133396706`,20.01564229416877`,20.01930998803875`,20.023559541074608`,20.028440480534698`}

Dieser Verlauf ist akzeptabel,  ansonsten muss man noch genauer rechnen.

(* direkter Vergleich beider Lösungen: *)

Plot[TLneu[0.1,t]-TL[0.1,t],{t,t0,tn/100},AxesLabel->{" z \t"," Abweichung \t"}]