Till Eulenspiegels lustige Serie/ Schwingende Objekte

Schwingende Objekte, Saiten, Balken und andere

Dieser Abschnitt handelt von Modellrechnungen, die mechanische Schwingungen und Wellen von drahtförmigen und ähnlich langgezogenen Objekten voraussagen wollen. Speziell sind die Saiten der Musikinstrumente betroffen. Es ist Physikunterricht hier und er stoppt nicht vor der zugehörigen angewandten Mathematik: Differenzial-, Integral-, Variationsrechnung, Wirkungsprinzip, Lagrange-Funktionen. Keine Panik bitte.

Zuerst kommt das einfachste Modell der gespannten Saite, das den Prototyp einer sogenannten Wellengleichung liefert. Wandernde und stehende Wellen und Eigenschwingungen auf einem gespannten Draht werden erklärt. Danach wird berücksichtigt, dass die Saiten wie Stangen und Balken nicht unendlich dünn sind und mit Widerstandskraft gegen das Verbiegen reagieren. Das verfeinerte Modell des Balkens ist schon sehr alt und stammt von Euler und Bernoulli. Es erklärt die inharmonischen Oberschwingungen dicker Saiten, die beim Klavierstimmen für die Spreizung von Oktaven berücksichtigt werden. Auch nichtlineare Effekte können da einziehen und bei großen Amplituden die Tonfrequenz anheben. Schließlich gibt es ein nochmal verfeinertes Modell aus dem zwanzigsten Jahrhundert, die Theorie von Timoshenko und Ehrenfest. Diese räumt ein, dass Elemente elastisch verformter dicker Drähte sich nicht nur senkrecht zur Achse, sondern auch je nach Lage im Querschnitt in Achsrichtung auslenken. Die zusätzlichen Freiheitsgrade rufen dabei eine Materialkonstante namens Schubmodul auf den Plan.

Die Hookeschen Gesetze

Es geht hier um elastische feste Körper. Wirkt eine Kraft an irgendeiner Stelle auf solche Körper ein, dann verformen sie sich proportional zur Größe der Kraft. Verschwindet die Kraft, kommt die ursprüngliche Form zurück. Diese Vorstellung funktioniert gut bis zu einem Grenzwert, ab dem sich das Material dauerhaft verformt, anfängt zu fließen und dann zu zerbrechen. Die schwingenden Körper seien immer im Geltungsbereich der Elastizität betrachtet, also bei hinreichend schwacher Belastung.

Unter anderem gibt es vier Parameter, die an einer elastischen Formveränderung von Klötzen von Metall beispielsweise zu messen sind.

• Erster Fall:

Der Klotz vom Volumen V kommt in eine Druckkammer und wird von allen Seiten gleichmäßig dem Druck p (Einheit 1 Pascal = 1 Newton/m²) ausgesetzt. Die Abnahme des Volumens ist proportional zum Volumen selbst und zum angelegten Druck. In Formeln, die dimensionslose Größe $-\Delta V/V$ ist proportional zur Druckänderung:

$-\Delta V/V=\chi \cdot \Delta p.$

χ ist die Kompressibilität; ihr Kehrwert K= 1/χ heißt der Kompressionsmodul und hat die Dimension Pascal, $(\Delta p=-K\cdot \Delta V/V).$

• Zweiter Fall:

Ein Quader von Länge mal Breite mal Höhe L x B x H wird auf seiner Deckfläche belastet mit einer senkrechten Kraft F, die gleichmäßig auf der Fläche A= L x B verteilt ist. Die Kraft pro Flächeneinheit S = F / A (Einheit Pascal) heißt in diesem Fall eine Spannung, nicht ein Druck. Das deshalb, weil sowohl die Kraft F wie auch die Normale der Bezugsfläche A verschiedene Richtungen haben können. Während ein Druck eine skalare Größe ist, sind zwei Vektoren nötig, um eine Spannung genau zu definieren. Man stelle sich einen Probestempel vor, der mit Sekundenkleber auf einem Stück Oberfläche haftet und dann in alle Richtungen gedrückt oder gezogen werden kann.

Das Hookesche Gesetz sagt aus, dass die Höhe des Quaders abnimmt proportional zur Höhe und zur vertikalen Spannung. Die Parameter Länge und Breite fallen heraus. In Formeln, $S=E\cdot (-\Delta H/H).$ Die Konstante E (Dimension Pascal) ist der Elastizitätsmodul.

Im Arsenal der Physik gibt es also Skalare, Vektoren und Größen, genannt Tensoren, deren Wert von mehreren Vektoren zugleich bestimmt ist. Die allgemeine Spannung ist ein Tensor. Im Inneren eines belasteten komplizierten Objekts, wie dem Kolben eines altmodischen Verbrennungsmotors, herrscht an jeder Stelle eine andere Spannung; es existiert insgesamt ein Tensorfeld.

• Dritter Fall:

Dieser ist nur ein Zusatz zum zweiten Fall. Derselbe Quader mit derselben Art Belastung wird etwas in die Breite gehen, weil nicht einzusehen ist, dass der ganze gerichtete Höhenverlust auch als verlorenes Volumen aufwartet. Länge oder Breite nehmen zu mit der Proportionalität:

$\Delta B/B=\Delta L/L=-\nu \cdot \Delta H/H.$

Der dimensionslose Parameter ν ist die Poissonzahl. Das Phänomen ist eine Querkontraktion, wenn die Spannung als Zugbelastung einwirkt.

• Vierter Fall:

Derselbe Quader, doch nun wird eine Schubspannung S an der Deckfläche mit einer horizontalen Kraft, etwa in x-Richtung ausgeübt, während die Grundfläche des Quaders festgenagelt bleibt. Der Quader wird schief, denn die Deckfläche wird ein Stück $\Delta x$ abwandern. Wieder ist der Effekt proportional zur Höhe und zur angelegten Spannung.

$S=G\cdot (\Delta x/H)=G\cdot \tan(\alpha )$

Die Konstante G (Dimension Pascal) ist der Schubmodul und der Effekt kann auch durch den Neigungswinkel α beschrieben werden.

Ein isotroper Festkörper hat in allen Richtungen die gleichen Eigenschaften. Anisotrope Objekte wie etwa faseriges Holz oder Kristalle haben elastische Parameter, die von den Richtungen abhängen. Für ein isotropes homogenes Material hängen die vier Werte $K,E,\nu ,G$ folgendermaßen zusammen:

$E=2G(1+\nu )=3K(1-2\nu ).$

Nur zwei sind unabhängig. Folglich gilt, wenn K sehr groß ist (inkompressible Materie):

$\nu \simeq 1/2\quad {\text{und}}\quad E\simeq 3G.$

Im allgemeinen Fall hängt eine innere Spannung an jeder Stelle eines kompliziert belasteten Objekts nichttrivial von einer gedachten Flächennormalen ab, sie ist kein bloßes Skalarprodukt mit einer einfachen Kraftrichtung. Sondern der Spannungsvektor $(S_{1},S_{2},S_{2})({\vec {x}})$ hat allgemeinere lineare Beziehungen zu jedem Normalenvektor $(n_{1},n_{2},n_{2}):\;S_{i}({\vec {x}})=\sum _{k}\sigma _{ik}({\vec {x}})\cdot n_{k}.$ Das Objekt σ ist das Cauchy-Tensorfeld, es hat 9 Komponenten (wegen Symmetrie nur sechs unabhängige) und variiert von Ort zu Ort. Dies ist wohl das historisch älteste und namensgebende Beispiel eines Tensors.

Mehr Elastizitäts-Formalismus

Der folgende Stoff wird nur weiter unten zum Timoshenko-Modell gebraucht. Wer das langweilig findet kann beide Abschnitte überspringen.

Cauchy-Spannungstensor σ hat das Indexpaar Flächennormale, Kraftrichtung

Cauchy hat die Kräfte auf einen belasteten Körper beliebiger Form als Summe von Oberflächenkräften beschrieben. An jedem infinitesimalen Testquader im Körper ziehen solche Kräfte auf gegenüberliegenden Flächen in Gegenrichtung, daher sollte am Ende nur ein Oberflächenintegral über den ganzen Körper bleiben. Aus der Impulserhaltung wird bewiesen, dass die Kraft auf beliebig orientierter Fläche eine Linearkombination aus Kräften für drei orthogonale Flächen sind. Genauer, die Kraft ist eine lineare Funktion der Flächennormalen und des Flächeninhalts. Aus der Drehimpulserhaltung an allen Testquadern wird bewiesen, dass die Kraft auf der Fläche A die Form hat

$f_{i}=A\;\sum \nolimits _{j=1}^{3}\sigma _{ji}n_{j};\quad {\vec {f}}=A\;\mathbf {\sigma } \cdot {\vec {n}}$

mit der Flächennormalen ${\vec {n}}$ und dem symmetrischen Cauchy-Spannungstensor $\sigma _{ij}=\sigma _{ji}$. Seine Komponenten haben die Dimension Pa. Seine Symmetrie garantiert, dass er kein Drehmoment hervorruft.

Die unverformten Koordinaten eines Punktes seien ${\vec {x}}=(x,y,z)$. Durch eine elastische Verformung wandert der Punkt an die Position ${\vec {x}}+{\vec {u}}({\vec {x}})$. Die Verformung oder Verzerrung ist die Tabelle der partiellen Ableitungen, genannt die Jacobimatrix, der Abbildung ${\vec {x}}\mapsto {\vec {u}}({\vec {x}}).\;{\vec {u}}({\vec {x}})$ misst, wie weit sich der Punkt ${\vec {x}}$ verschiebt. Ihre Verzerrung dagegen gibt linear angenähert an, wie weit nahe gelegene Punkte sich von einander weg oder aufeinander zu entwickeln.

$F_{ik}=\partial u_{i}/\partial x_{k};\;u_{i}({\vec {x}}+{\vec {d}})\simeq u_{i}({\vec {x}})+\sum \nolimits _{k=1}^{3}F_{ik}({\vec {x}})\;d_{k}$

Das allgemeine Hookesche Gesetz besagt, dass diese Verformung in einer linearen Beziehung zur Cauchy-Spannung steht. Nun stellt jede Abbildung, deren Jacobimatrix antisymmetrisch ist, eine infinitesimale Drehung der benachbarten Punkte um den Auswertungspunkt dar. Ein symmetrischer Cauchy-Tensor entwickelt keine Drehmomente. Reine Drehungen rufen keine elastischen Kräfte auf den Plan. Daher kann nur der Verzerrungstensor, der symmetrische Anteil der Jacobimatrix, zur Hookeschen Regel beitragen. In Koordinaten:

$\epsilon _{ik}:=(F_{ik}+F_{ki})/2;\quad \sigma _{ij}=\sum \nolimits _{k,l}C_{ijkl}\;\epsilon _{kl}$

Zum Glück bleiben von den 81 C-Koeffizienten nur zwei übrig in den homogenen und isotropen Medien, die hier auf dem Programm stehen.

Mathematischer Exkurs. Drehungen konvergieren bei kleinen Winkeln zur Identität plus einer schief- oder anti-symmetrischen Matrix, wegen

${\begin{pmatrix}\cos w&\sin w\\-\sin w&\cos w\end{pmatrix}}\to {\begin{pmatrix}1&w\\-w&1\end{pmatrix}}.$

In allen höheren Dimensionen erfüllen Drehungen als orthogonale Matrizen R die Bedingung $R^{T}R=I$ (Einheitsmatrix).

Auf einer glatten Matrixfunktion $R(x);\;R(0)=I$ liefert das Differenzieren die Antisymmetrie:

$0=(d/dx)[R(x)^{T}R(x)]=(R')^{T}R+R^{T}R'=(R')^{T}+(R')\;{\text{ am Ursprung }}\;x=0.$

Das heißt, die lineare Näherung $R(x)=I+xR'(0)$ einer Rotation ist schiefsymmetrisch.

Die Lehre von den kontinuierlichen Transformationen und von ihren sogenannten infinitesimalen Versionen wurde voll ausgebaut zur mathematischen Theorie der Lie-Gruppen und Lie-Algebren.

In folgenden Formeln sind R,T,U orthogonale Matrizen. Jede symmetrische reelle Matrix $S\;(S^{T}=S\;\Leftrightarrow \;S_{ik}=S_{ki}\;\forall i,k)$ kommt mit orthogonalen Ähnlichkeitstransformationen auf Diagonalform, $RSR^{-1}=RSR^{T}=D$, und die Zahlen auf der Diagonalen sind die reellen Eigenwerte. Dieses wichtige Resultat hat Korollare.

Jede quadratische Matrix kann als Produkt einer orthogonalen und einer symmetrischen Matrix geschrieben werden. Zumindest für umkehrbare Matrizen M folgt das daraus, dass M^TM symmetrisch und positiv definit ist und daher eine umkehrbare Quadratwurzel hat.

$M^{T}M=RDR^{-1}=RE^{2}R^{T}=(RER^{T})(RER^{T})=S^{2};\;E={\sqrt {D}};\;S=RER^{-1}=RER^{T}.$

$U:=MS^{-1}$ ist dann orthogonal, denn auch S^-1 ist symmetrisch:

$(MS^{-1})^{T}(MS^{-1})=S^{-1}(M^{T}M)S^{-1}=S^{-1}SSS^{-1}=I\;\Rightarrow \;M=US.$

Mit einer Diagonalisierung $S=T^{-1}DT$ kann dann die Matrix M weiter in Produkte von drei Matrizen zerlegt werden, $M=U\;S=(UT^{-1})\;D\;T$, so dass zwischen zwei Drehungen eine Diagonalmatrix D steht. Diese D macht die eigentliche Arbeit der Verformung. Sie verwandelt Kugeln in Ellipsoide, welche orthogonale Hauptachsen in Richtung der Koordinatenachsen haben.

Auf dem komplexen Zahlenkörper sehen die privilegierten Matrizen etwas anders aus. Die Eigenschaften {Symmetrisch, Orthogonal} werden zu {Hermitesch, Unitär}.

Abbildung $\mathbf {x} \mapsto \mathbf {u}$ verformt kleine Quadrate zu Parallelogrammen. Partielle Ableitungen messen die relative Verschiebung.

• Spezialfall von σ, ε bei Standardparametern E,G,ν.

Ein Testquader sei mit einer Ecke bei (0,0,0) festgenagelt und habe die Dimension $V=(\Delta x\times \Delta y\times \Delta z)$.

Die Scherspannung $\sigma _{yx}$ setzt auf einer xz-Fläche an (Normale in y-Richtung, erster Index) und wirkt in Richtung x-Achse, anderer Index. Punkt (x,y,z) wird verschoben nach ${\vec {u}}({\vec {x}})=(\delta x,0,0)$, so dass

$G\;\delta x=y\cdot \sigma _{yx}$.

Es gilt $F_{xy}=\delta x/y=\sigma _{yx}/G.$ Laut Postulat soll die Spannung linear mit $\epsilon _{xy}$ zusammenhängen. Weil hier

$F_{yx}=0,\;\Rightarrow \;\epsilon _{xy}=F_{xy}/2;\;\sigma _{yx}=2G\;\epsilon _{xy}.$

Weil es intuitiv nicht so einleuchtet, noch ein Erklärungsversuch. Die Verformung, die zur Definition des Schubmoduls herhielt, zählt nur halb als eine gegenseitige Bewegung der Materie. Die andere (schiefsymmetrische) Hälfte ist eine kleine Drehung ohne Einfluss auf die Form und wird abgezogen. Daher nun: Schubspannung = 2G mal eigentliche, effektive Verzerrung.

Eine Zugspannung $\sigma _{xx}$ auf einer yz-Fläche verlängert die x-Dimensionen mit dem Elastizitätsmodul (1/E) und verkürzt die Dimensionen y,z um den Faktor (ν/E). Also:

${\vec {u}}(x,y,z)=\sigma _{xx}\;(x/E,-\nu \;y/E,-\nu \;z/E)\;\Rightarrow \;(\epsilon _{xx},\epsilon _{yy},\epsilon _{zz})=(1,-\nu ,-\nu )\cdot \sigma _{xx}/E$

Mit Beiträgen aller sigma-Komponenten gilt mit Indizes (xyz)=(123) für die Verzerrung

$\epsilon _{11}=\sigma _{11}/E-\nu (\sigma _{22}+\sigma _{33})/E$

$\epsilon _{12}=\sigma _{12}/(2G)$

Symmetrische Formeln betreffen die anderen Indexpaare des ε-Tensors.

Mit der Drehsymmetrie des isotropen Materials sei nun gezeigt, dass die drei Materialkonstanten E,G,ν nicht unabhängig sind. Gespielt wird in der xy-Ebene mit einer symmetrischen spurfreien Cauchy-Spannung. Unter einer Rotation, dargestellt als die 3x3-Matrix R, transformieren sich folgende Objekte alle gleich: Kräfte,Flächennormalen,Ortsvektoren,freie Vektoren

${\vec {f}},\;{\vec {n}},\;{\vec {x}},\;{\vec {d}}\;\mapsto R{\vec {f}},\;R{\vec {n}},\;R{\vec {x}},\;R{\vec {d}}.$

Die Matrizen für σ, ε verbinden Vektor-Objekte miteinander, etwa

${\vec {f}}=\mathbf {\sigma } \;{\vec {n}}\;{\text{ oder }}\;{\vec {u}}=\mathbf {\epsilon } \;{\vec {d}}.$

Damit solche Beziehungen nach der Rotation zutreffen, müssen die Transformationen gelten

$\mathbf {\sigma } ,\;\mathbf {\epsilon } \;\mapsto \;R\;\mathbf {\sigma } \;R^{-1},\;R\;\mathbf {\epsilon } \;R^{-1}.$

Denn ${\vec {y}}=M{\vec {x}}\;\Rightarrow \;(R{\vec {y}})=(RMR^{-1})(R{\vec {x}})$. Hier meint $R^{-1}$ die inverse Matrix. Bei Rotationen, dargestellt als orthogonale Matrizen, ist sie netterweise identisch mit der transponierten Matrix $R^{T}$.

Wir drehen um 45 Grad im Uhrzeigersinn in der xy-Ebene und bearbeiten symmetrische 2x2-Matrizen.

$R={\begin{pmatrix}\cos \phi &\sin \phi \\-\sin \phi &\cos \phi \end{pmatrix}};\;\phi =-\pi /4;\;R={\begin{pmatrix}1&-1\\1&1\end{pmatrix}}/{\sqrt {2}};\;R^{-1}={\begin{pmatrix}1&1\\-1&1\end{pmatrix}}/{\sqrt {2}}.$

$M={\begin{pmatrix}a&b\\b&c\end{pmatrix}}\quad \mapsto \;RMR^{-1}={\frac {1}{2}}\cdot {\begin{pmatrix}a+2b+c&-a+c\\-a+c&a-2b+c\end{pmatrix}}$

${\begin{array}{lll}{\begin{pmatrix}\sigma _{11}&\sigma _{12}\\\sigma _{21}&\sigma _{22}\end{pmatrix}}&={\begin{pmatrix}a&b\\b&-a\end{pmatrix}}&\mapsto \;{\begin{pmatrix}b&-a\\-a&-b\end{pmatrix}}\\{\begin{pmatrix}\epsilon _{11}&\epsilon _{12}\\\epsilon _{21}&\epsilon _{22}\end{pmatrix}}&={\begin{pmatrix}a(1+\nu )/E&b/(2G)\\b/(2G)&-a(1+\nu )/E\end{pmatrix}}&\mapsto \;{\begin{pmatrix}b/(2G)&-a(1+\nu )/E\\-a(1+\nu /E&b/(2G)\end{pmatrix}}.\end{array}}$

Die allgemein postulierte Form der Epsilon-Matrix kann also nur dann im gedrehten Koordinatensystem zutreffen, wenn gilt, was der vorige Abschnitt unbewiesen erwähnte: $(1+\nu )/E=1/(2G)$. Das zeigt der Vorher-Nachher-Vergleich der Spannungen und Verzerrungen. Anders gesagt, E=2(1+ν)G.

Der Parameter ν kann daher mit E,G eliminiert werden: $\nu =E/(2G)-1$.

$\epsilon _{11}=\sigma _{11}/E+(1/E-1/(2G))\cdot (\sigma _{22}+\sigma _{33})$

$\epsilon _{ik}=\sigma _{ik}/(2G)\quad (i\neq k)$

$\epsilon _{ii}=\sigma _{ii}/(2G)+(1/E-1/(2G))\cdot (\sigma _{11}+\sigma _{22}+\sigma _{33})$

In der Zeile zuletzt taucht die Spur von σ auf, so heißt die Diagonalsumme.

Noch kurz die Relation für den Kompressionsmodul K geprüft. Ein Druck p entspricht der Cauchy-Spannung $\sigma _{xx}=\sigma _{yy}=\sigma _{zz}=-p$ mit Spur -3p. Die relative Änderung des Volumens ist die Summe der diagonalen $\epsilon _{ii}$,

also $(\Delta V/V)=3\cdot [(-p)/(2G)+(-3p)(1/E-1/(2G))]=(3/G-9/E)\;p.$

Vergleich mit der Definition: $p=K\cdot (\Delta V/V):\quad K=EG/(3E-9G).$

Einsetzen von $G=E/(2+2\nu )\;\Rightarrow \quad K=E/[3(1-2\nu )].$

• Die Energiedichte.

Der Testquader sei wieder $V=(\Delta x\times \Delta y\times \Delta z)$ mit kleinen Dimensionen, so dass lineare Näherungen gut zutreffen. Sein Mittelpunkt wird von ${\vec {x}}\;{\text{ nach }}\;{\vec {x}}+{\vec {u}}({\vec {x}})$ verschoben. Der Inhalt wird linear so verformt, dass seine Punkte in Koordinaten ${\vec {d}}$ relativ zu seinem Zentrum, die lineare Transformation erfahren:

$d\mapsto Fd;\;{(Fd)}_{i}=\sum \nolimits _{k}(\partial u_{i}/\partial x_{k})\;d_{k};\;F_{ik}=\partial u_{i}/\partial x_{k}.$

Verallgemeinert nach Hooke bewirkt die Spannung $\sigma _{jl}$ auf der Fläche mit Normalenrichtung (l) in Achsrichtung (j) lineare, windschiefe Verformungen des Quaders in allen 3 Koordinaten (i). Jede Koordinate im Quader $d_{i}$ bekommt Schübe proportional zu den Abständen $d_{k}\;{\text{ von }}\;{\vec {d}}$ zum relativen Ursprung,

$d_{i}\mapsto \sum \nolimits _{k}(S_{ikjl}\sigma _{jl})d_{k}.$

Also pro Spannungselement eine allgemeine lineare Funktion von ${\vec {d}}$. Nach Summierung über alle Paare (j,l) zeigt der Vergleich mit der Formel der Jacobimatrix:

$F_{ik}=\epsilon _{ik}=\sum \nolimits _{j,l}S_{ikjl}\;\sigma _{jl}.$

Hier wird die Symmetrie von σ vorausgesetzt, also verschwindende Drehmomente, und eine 'drehfreie' Matrix $F_{ik}=F_{ki}=\epsilon _{ik}$. Es gilt eine lineare Umkehrformel:

$\sigma _{jl}=\sum \nolimits _{i,k}C_{jlik}\;\epsilon _{ik}.$

Was ist die potenzielle elastische Energie, die im Testquader gespeicht wird, wenn bekannte Spannungen $\{\sigma _{jl}\}$ an den Flächen angelegt werden? Die Mittelpunkte der Flächen k liegen bei ${\vec {d}}_{k}=(\pm \Delta x_{k}/2)\cdot {\vec {n}}_{k}$ mit Normalenvektor ${\vec {n}}_{k}$. Sie wandern nach $(F{\vec {d}}_{k})_{i}=(\pm \Delta x_{k}/2)F_{ik}.$

Mit dem Flächeninhalt $A_{k}$ unterliegen die zwei Flächen senkrecht zur Achse (k) dem Kräftepaar $({\vec {f}}_{k})_{i}=\pm A_{k}\sigma _{ik}$. Nach Hooke ist die Matrix σ eine Funktion der Verzerrungsmatrix [F]. Die Verschiebungswege der Flächenzentren, Angriffspunkte von Kraft, sind ebenfalls lineare Funktionen von [F]. Daher wird die eingespeiste Arbeit, Summe von Kraft mal Weg, wie folgt berechnet. In der Zeit t=0 bis t=1 wird die Verzerrung linear hochgefahren gemäß $t\mapsto tF$. Die Arbeit ist das Zeitintegral

$W=\int _{0}^{1}dt\;{\vec {f}}_{k}(tF)\cdot (d/dt)(tF\;{\vec {d}}_{k})$, Integral von Kraft mal Geschwindigkeit.

Die Abhängigkeit von t ist linear in beiden Faktoren, weshalb das Integral den Faktor $\int _{0}^{1}t\;dt={\frac {1}{2}}$ ergibt. Summiert über alle Flächen und Gegenflächen,

$W={\frac {1}{2}}\sum \nolimits _{i,k}f_{ik}(F{\vec {d}}_{k})_{i}={\frac {1}{2}}\sum \nolimits _{i,k}(A_{k}\;\Delta x_{k})\sigma _{ik}\cdot (\partial u_{i}/\partial x_{k}).$

Der erste Faktor ist das Volumen des Testquaders. Wegen der Symmetrie von σ können die Ableitungen symmetrisiert werden. Es ergibt sich die Energiedichte pro Volumeneinheit,

$(W/V)={\frac {1}{2}}\sum \nolimits _{i,k}\sigma _{ik}\cdot \epsilon _{ik}={\frac {1}{2}}\sum \nolimits _{i,k,j,l}\epsilon _{ik}\;C_{ikjl}\;\epsilon _{jl}.$

Die elastische Energiedichte ist eine quadratische Form der Verzerrungsmatrix, deren Elemente später als dynamische Feldvariablen dienen. In einem beliebigen elastischen Körper variieren $\;\mathbf {\epsilon } ({\vec {x}})\;{\text{ und }}\;\mathbf {\sigma } ({\vec {x}})\;$ von Punkt zu Punkt und die gesamte gespeicherte elastische Energie ist das Volumenintegral der lokalen Energiedichte.

Die gespannte Saite

Eine Saite ist etwa ein Draht aus Stahl oder ein Faden aus Katzendarm, der an zwei Enden festgehalten wird, die Länge L hat und mit einer kontrollierten Kraft F gespannt wird. Hat sie die Schnittfläche A und die Materialdichte ρ, dann ist $\mu =\rho \cdot A$ die Masse der Saite pro Längeneinheit. Es zeigt sich, dass die Schwingungen und Wellen auf der ausgelenkten Saite nur von F und μ abhängen, aber in erster Ordnung gar nicht vom Elastizitätsmodul. Nur bei dicken Saiten werden wir eine Korrektur berechnen.

Erste Herleitung der Wellengleichung. Es geht um die Mechanik eines elastischen Gegenstands. Jeder Massenpunkt bewegt sich nach Newtons Gesetz: Masse mal Beschleunigung gleich Summe aller einwirkenden Kräfte. Eine Kraft wird in Newton gemessen, was dimensionsgleich zu [kg m / s²] ist wegen der Newton-Gleichung. Die Erdanziehung auf ein Gewicht der Masse 1 kg beträgt 9,81 Newton, denn die Erdbeschleunigung beim freien Fall ist 9,81 m/s².

Die Saite sei in x-Richtung gespannt. Jeder Massenpunkt im Intervall x=0 bis L darf sich nur in z-Richtung auslenken. Eine momentane vertikale Auslenkung der Saite ist eine Funktion w(x,t) von zwei Variablen, der Position und der Zeit t. Es bestehen die Randbedingungen w(0,t) = w(L,t) = 0. In Ruhestellung w(x,t)=0 ist die Saite mit der horizontalen Kraft F vorgespannt. Der Draht habe die Dichte ρ und den Radius r; also den Querschnitt A= π r² und die Masse pro Längeneinheit μ = ρ A. Ein kleines Stück Saite von x bis x+ Δx hat die Masse M = μ Δx und die vertikalen Positionen von w(x,t) bis w(x+ Δx,t). Es wird angenommen, dass w() als Funktion von x beliebig aber glatt verläuft und dass die Ableitungen nach der Variablen x ebenfalls glatt verlaufen. Der Winkel, mit dem die Saite am Punkt x momentan schief steht zur Horizontalen ist gegeben durch die Ableitung

$\tan \alpha (x)=w'(x)=\partial w/\partial x.$

Die Saite kann nur deshalb schief verlaufen, weil sie eine elastische Ausdehnung hat und weil eine Kraft ${\tilde {F}}(x)$ genau in der Richtung von w'(x) angreift, die betragsmäßig größer ist als die Horizontalkraft F. Das Massenstück von x bis x+ Δx wird am linken Punkt horizontal von der Kraft -F gezogen und vertikal mit der richtungsmäßig passenden Komponente von ${\tilde {F}}(x),$ nämlich zwangsläufig $-F\tan \alpha =-Fw'(x)$. Am rechten Punkt zieht horizontal die Kraft F nach rechts und es muss vertikal $Fw'(x+\Delta x)$ vorhanden sein. Die horizontalen Kräfte auf dem Drahtstück gleichen sich aus, aber die kleine Differenz der vertikalen Kräfte bewirkt nach Newton eine Beschleunigung, also eine zweite Zeitableitung von w(x,t) am Mittelpunkt:

$\mu \Delta x{\frac {\partial ^{2}w(x+\Delta x/2,t)}{\partial t^{2}}}=F(-w'(x,t)+w'(x+\Delta x,t))$

Wird nun durch Δx geteilt und dessen Grenzwert gegen Null genommen, kommt auf der rechten Seite die zweite Ableitung nach x heraus:

$\mu {\frac {\partial ^{2}w(x,t)}{\partial t^{2}}}=\mu {\ddot {w}}=Fw''=F{\frac {\partial ^{2}w(x,t)}{\partial x^{2}}}.$

Diese Differenzialgleichung ist eine Wellengleichung. Die Notation mit den runden $\partial$ ist üblich, wenn eine Funktion von mehreren Variablen abhängt und nur nach einer bestimmten davon abgeleitet wird, während die anderen Argumente festgezurrt bleiben. Es ist eine partielle Ableitung. Wir erlauben uns manchmal, die partiellen Ableitungen nach der Zeit t mit Punkten und die nach der Position x mit Strichen anzudeuten. Hier gleichbedeutend:

$\{{\dot {w}},w'\}\equiv \{\partial w/\partial t,\partial w/\partial x\}\equiv \{\partial _{0}w,\partial _{1}w\}\equiv \{\partial w/\partial x_{0},\partial w/\partial x_{1}\}.$

Umgeformt gilt einfach: ${\ddot {w}}=(F/\mu )\cdot w''.$ Der Faktor in der Klammer hat die Dimension (Länge durch Zeit) zum Quadrat. Also definiert er eine Geschwindigkeit:

$c:={\sqrt {(F/\mu )}};\;{\ddot {w}}=c^{2}\;w''.$

Hat nun w(x,t) die spezielle Form w(x,t)=u(x-ct) oder auch w(x,t)=v(x+ct), dann wird die Wellengleichung erfüllt. Denn die zweiten Ableitungen sind nach der Kettenregel auszurechnen und beispielsweise

${\ddot {w}}=c^{2}v''.$

Die Funktionen u(x) und v(x) könne beliebige Buckel oder Wellenpakete sein. u(x-ct) beschreibt dann eine wandernde Welle mit Geschwindigkeit c nach rechts, v(x+ct) eine wandernde Welle nach links. Beliebig geformte Wellen können sich in beiden Richtungen auf dem gespannten Draht fortpflanzen. Größere Spannung macht die Wellen schneller, größere Masse der Saite macht sie langsamer.

Nun soll aber der Draht als Saite an beiden Enden festgehalten werden und daher sollen keine Wanderwellen hindurch laufen. Die Lösungen der Gleichung sind dann stehende Wellen, die bei x=0 und x=L die Randbedingungen w(x,t)=0 erfüllen. Dies nennt man eine Randwertaufgabe der Differenzialgleichung. Unsere Gleichung ist linear und homogen; als erste Idee kommt immer ein Separationsansatz auf:

$w(x,t)=f(x)\cdot g(t).$

$f\cdot {\ddot {g}}=c^{2}f''\cdot g.\;$

Es ist lösbar wenn

${\ddot {g}}/g=c^{2}f''/f=\;$ eine Konstante z, unabhängig von (x,t).

Sei also die Gleichung in einer Variablen

${\ddot {g}}(t)=z\cdot g(t).\;$

Einige elementare Funktionen lösen sie:

$g(t)={\exp(\alpha t),\exp(-\alpha t),\sin(\omega t),\cos(\omega t)}.$

Weil etwas mit Schwingungen auf dem Programm steht, sieht es besser aus mit den trigonometrischen Funktionen und deren Ableitungsregeln:

$(d/dt)\sin(\omega t)=\omega \cdot \cos(\omega t);\quad (d/dt)\cos(\omega t)=-\omega \cdot \sin(\omega t).$

Damit folgt $z=-\omega ^{2}.\;\omega$ ist die Kreisfrequenz der zeitperiodischen Funktion $g(t),\;f=\omega /(2\pi )$ ist die eigentliche Frequenz und $T=1/f$ ist die Periode der Schwingung.

Bleibt noch der x-Anteil der Saitenbewegung zu diskutieren:

$c^{2}f''=-\omega ^{2}f$. Die Lösung ist auch eine Sinusfunktion $f(x)=\sin(kx),$

Hier wird f(0)=f(L)=0 verlangt. Das bedeutet, dass kL ein Vielfaches von π sein muss: kL = nπ. Andererseits folgt $-c^{2}k^{2}f=-\omega ^{2}f$,

also $\omega =ck=nc\pi /L$, oder

für die Frequenzen $f_{n}=nc/(2L)$ mit n=1,2,3...

Ergebnis: die Saite kann sich mit Eigenschwingungen bewegen, deren Frequenzen ganze Vielfache der Grundfrequenz sind. Die Oberschwingung der Ordnung n hat n Bäuche längs der Seite und (n-1) Knotenpunkte, die sich nicht bewegen. Die Formel der Frequenzen in Parametern $L,\mu ,F$ lautet demnach $f_{n}=n{\sqrt {(F/\mu )}}\;/(2L)$. Die Frequenz steigt, wenn die Saite kürzer gemacht wird, oder die Spannkraft vergrößert, oder die Masse verkleinert. Wenig zählt, aus welchem Stoff die Saite ist. Die Masse wächst nun wie das Quadrat des Durchmessers D der Saite, die Wurzel der Masse also wie D. Daher folgt die Grundschwingung dem gefühlt seit Jahrtausenden bekannten Gesetz

$f\propto {\sqrt {F}}/(L\cdot D).$

Siehe dazu Pythagoras in der Schmiede.

Andere Herleitung der Saitenschwingung

Es mag erstaunen, dass die Saitenformel gleich ist für Katzendarm, Nylon und Stahl und dass das Hookesche Gesetz mit seiner Materialkonstanten E nicht eingeht. Daher eine Betrachtung dazu, wann und wie E sich wegkürzt. Wird ein Klotz mit Grundfläche A und Höhe h mit Druck oder Zug belastet, hängt die Kraft F(z) mit der Höhenänderung z zusammen: F(z)/A = E (z/h). Die Arbeit, die geleistet wird, wenn die Form von 0 nach z gestreckt/gepresst wird, ist das Integral der Kraft als Funktion des Wegs:

$W(z)=\int dz'\;F(z')={\frac {EA}{h}}\cdot {\frac {z^{2}}{2}}$.

Dies ist die gespeicherte potenzielle Energie, die wegen der perfekten Elastizität wiederverwertbar sein soll.

Die kinetische Energie einer Saite, die an einem gewissen Zeitpunkt t eine Auslenkung mit der Funktion w(x,t) durchläuft, ist die Summe der Energien (mv²/2) aller "Atome" oder infinitesimalen Teilstücke der Länge Δx, nämlich $(\mu \;\Delta x/2)\cdot (\partial w(x,t)/\partial t)^{2}$. Der Grenzwert einer solchen Summe, wenn die Länge der Stücke gegen Null geht, ist das Integral

$T(t)=\int _{0}^{L}(\mu /2)\cdot (\partial w(x,t)/\partial t)^{2}\;dx=\int _{0}^{L}{\frac {\mu {\dot {w}}^{2}}{2}}\;dx$.

Auch die potenzielle Energie soll nun aus ihren Elementen aufsummiert werden.

Die Saite ist in Ruhestellung w(x,t)=0 mit der Kraft F vorgespannt. Das heißt, das Längenstück $\Delta x=h+z=:h\cdot (1+\delta )$ besteht aus der kraftfreien Länge h und der Streckung $z=h\cdot \delta$. Das Hookesche Gesetz besagt $F=EA\cdot \delta$ und die gespeicherte Energie in Stück Δx ist, wegen h=Δx/(1+δ),

$V_{0}\;\Delta x=(EA/h)\cdot (z^{2}/2)={\frac {F}{\delta }}\cdot {\frac {h\delta ^{2}}{2}}={\frac {F\;\Delta x}{\delta (1+\delta )}}\cdot {\frac {\delta ^{2}}{2}}.$

Der Wert von δ hängt nicht vom Intervall Δx ab, sondern er enthält die Kraft und die Eigenschaft des Materials: δ=F/(EA).

Ist nun die Saite beliebig vertikal ausgelenkt mit der Funktion w(x,t), dann steht das Stück über der Horizontalen Δx schräg mit der Steigung $\partial w(x,t)/\partial x$. Die Streckung betrage nun $h\cdot (1+\epsilon )$, wo ε größer ist als δ. Die potenzielle Energie im Segment Δx wird

$V_{1}(x,t)\;\Delta x={\frac {F\;\Delta x}{\delta (1+\delta )}}{\frac {\epsilon ^{2}}{2}}.$

Mit der Steigung der Funktion w(x,t) und dem Satz von Pythagoras:

$(1+\epsilon )^{2}=(1+\delta )^{2}(1+w'(x,t)^{2})\;\Rightarrow \;\epsilon =(1+\delta ){\sqrt {1+w'^{2}}}-1$

Mit der Abkürzung $u(x,t):={\sqrt {1+(\partial w(x,t)/\partial x)^{2}}}$ schreibt sich die potenzielle Energie der Saite so, mit abgezogener Nullpunktsenergie:

$V(x,t)\;\Delta x=(V_{1}(x,t)-V_{0})\;\Delta x={\frac {F\;\Delta x}{2\delta (1+\delta )}}\cdot (\epsilon ^{2}-\delta ^{2}).$

Die letzte Klammer wird:

$(\epsilon +\delta )(\epsilon -\delta )=[(1+\delta )u+(\delta -1)]\cdot [(1+\delta )u-(1+\delta )]$

Es kann einmal (1+δ) gekürzt werden und man bekommt:

$V(x,t):={\frac {F}{2\delta }}(u-1)[(u-1)+\delta (u+1)];\quad u={\sqrt {(1+w'^{2})}}$

Nun sei angenommen, dass die Steigung der Saite so klein ist, dass

• erstens gut angenähert: $u-1={\sqrt {(1+w'^{2})}}-1=w'^{2}/2;\quad u+1=2$
• zweitens höhere Potenzen von (u-1) vernachlässigt werden.

Dann kommt als Dichte der potenziellen Energie des Kleinsignals heraus

$V_{k}(x,t)={\frac {F\;w'(x,t)^{2}}{2}}$

Diese hängt also nicht mehr von materialspezifischer Dehnung δ ab.

Die gesamte potenzielle Energie der Saite ist zur Zeit t das Integral über x

$V(t)=\int _{0}^{L}{\frac {F}{2\delta }}(u(x,t)-1)[(u(x,t)-1)+\delta (u(x,t)+1)]\;dx$

Die Hilfsfunktion u(x,t) misst die die Länge des Graphs der Funktion w(x,t). Denn $u(x,t)\cdot \Delta x$ ist die Länge eines verformten Saitenstücks über der horizontalen Achse, $L_{w}(t)=\int _{0}^{L}u(x,t)\;dx$ ist die Gesamtlänge der Saite zum Zeitpunkt t.

Bis hierhin wird nach der Analyse der kinetischen und der potenziellen Energie klar, dass die Saite nur bei kleinen Auslenkungen mit quadratischen Formen darin beschrieben wird. Bei großen Amplituden tauchen nichtquadratische Terme auf und der Elastizitätsmodul des Materials spielt auch mit.

Aber wo bleibt die Bewegungsgleichung der Saite, also eine Differenzialgleichung für w(x,t) mit partiellen Ableitungen nach allen Variablen? Sie kommt jetzt, mit Hilfe eines der wichtigsten mathematischen Sätze der ganzen Mechanik.

Das Wirkungsprinzip (Variationsprinzip der kleinsten Wirkung)

Die Saite ist ein dynamisches System mit einem Kontinuum von Freiheitsgraden. Ihr Verhalten in Ort und Zeit wird durch eine Funktion w(x,t) beschrieben. Die Funktion w(x,t) der Auslenkung ist der Spezialfall einer Feldfunktion f(x) = f_i(x_k) eines dynamischen Systems. Die Koordinaten x_k sind oft die Vierervektoren $(x_{0},x_{1},x_{2},x_{3})=(t,{\vec {x}})$, welche die Raumzeitpunkte des Systems beschreiben. Die Zahl der Komponenten von f_i() hängt davon ab, ob es Skalare, Vektoren, Tensoren oder andere geometrische Gebilde sind, die an jedem Raumzeitpunkt dynamisch veränderliche Werte annehmen.

Das Wirkungsprinzip der Kontinuumsmechanik behauptet nun folgendes:

• Es gibt ein Wirkungsfunktional auf der Menge der Feldfunktionen f(), also eine Abbildung, die jedem f einen reellen Wert W{f} zuordnet.
• Eine Funktion f() beschreibt genau dann eine physikalisch mögliche Entwicklung des Systems in Raum und Zeit, wenn W{f} einen Extremwert hat bezüglich aller (genügend kleinen) Variationen von f.
• Das Funktional ist lokal, indem es ein Raumzeit-Integral über eine Dichtefunktion ist: $W\{f\}=\int d^{n}x\;L(x,f(x),\partial _{i}f(x),\partial _{i}\partial _{k}f(x)).$

Die Dichtefunktion oder Lagrange-Dichte L am Raumzeit-Punkt x hängt ab:

vom Punkt x und von den Werten f(x)

und von den ersten, zweiten (und höheren) partiellen Ableitungen von f, ausgewertet am selben Punkt x.

• Die Lagrange-Dichte der Newtonschen Mechanik ist die Differenz L = T - V, Dichte der kinetischen minus Dichte der potenziellen Energie.

Beispiel: Die Lagrange-Dichte unserer Saite am Punkt (x,t) ist aus den oben eingeführten Dichten für T(t) und V(t) aufgebaut:

$L(w(x,t),w'(x,t),{\dot {w}}(x,t))={\frac {\mu {\dot {w}}^{2}}{2}}-{\frac {F}{2\delta }}(u(x,t)-1)[(u(x,t)-1)+\delta \cdot (u(x,t)+1)];$

$u(x,t)={\sqrt {1+w'(x,t)^{2}}};\;\delta =F/(EA).$

Die Wirkung ist das Doppelintegral $W\{w\}=\iint dx\;dt\;L(w,w',{\dot {w}})$.

Im Grenzfall kleiner Auslenkungen ist der Integrand L quadratisch:

$L=(\mu /2)\cdot {\dot {w}}(x,t)^{2}-(F/2)\cdot w'(x,t)^{2}$

Das Integral über die Ortskoordinate ergibt dimensionsmäßig eine Energie. Das Doppelintegral über Ort und Zeit hat die Dimension einer Wirkung, nämlich Energie mal Zeit. Daher kommt der Name Wirkungsprinzip.

Das Potenzial V sieht mit seinem Zusatzterm zum quadratischen schöner so aus:

$V={\frac {F}{2\delta }}[\delta \cdot (u^{2}-1)+(u-1)^{2}]={\frac {F}{2}}{\Big [}w'^{2}+{\frac {(u-1)^{2}}{\delta }}{\Big ]}$

Wann ist der Störterm wirklich vernachlässigbar?

Grob gesagt, wenn w'² sehr viel kleiner ist als die bereits winzige Dehnung δ! Denn (u-1)² geht los mit der vierten Potenz von w' .

Zahlenbeispiel.

Kraft F=10 Newton, Querschnitt A=1 mm², Elastizitätsmodul E=200 GPa, Länge 500 mm, Auslenkung ~ 1 mm. δ = 10/(200e9 x 1e-6) = 5 x 10^-5.

Sei die maximale Steigung w' ~ (1/200), ihr Quadrat ~ 2,5 x 10^-5. Sie kommt da schon gefährlich nahe ins nichtlineare Gebiet dieses Modells.

Die Suche nach Extremwerten von Funktionalen ist als Variationsrechnung bekannt und ist eine Verallgemeinerung der Suche nach Maxima und Minima von gewöhnlichen Funktionen. Eine gewöhnliche Funktion y(x) von mehreren Variablen ist extremal am Punkt x, notwendige aber nicht hinreichende Bedingung, wenn alle partiellen Ableitungen an diesem Punkt verschwinden. Ein (glattes) Funktional, das als Integral einer Dichte definiert ist, kandidiert dann als Extremwert für die Funktion f, wenn f() eine partielle Differenzialgleichung erfüllt, die aus der Lagrange-Dichte entsteht. Sie heißt die Euler-Lagrange-Gleichung. Die Bewegungsgleichungen der Mechanik sind die Euler-Lagrange-Gleichungen des Wirkungsfunktionals. Genau genommen, es kommen so viele EL-Gleichungen vor wie reelle Komponenten im Wertebereich des Feldes f.

Die relevanten Feldsysteme der Physik, Mechanik, Elektrodynamik, die Dynamik der Gravitationsfelder, also die allgemeine Relativität, sowie das Standardmodell der Elementarteilchen, sie alle haben Euler-Lagrange-Gleichungen.

Nehmen wir einfachheitshalber f(x) einkomponentig an und definieren: Die Variation W'(f,x) des Funktionals $f\mapsto W\{f\}$ existiert, wenn die Werte von W für alle wenig abweichenden Funktionen $(x\mapsto f(x)+\epsilon (x))$ sich linear mit demselben Integralkern W' approximieren lassen:

$W\{f+\epsilon \}\simeq W\{f\}+\int d^{n}x\;W'(f,x)\cdot \epsilon (x)=:W\{f\}+\int d^{n}x\;{\frac {\delta W}{\delta f(x)}}\cdot \epsilon (x)$

Ist f ist ein Extremwert des Funktionals, dann verschwindet die Variation W'(f,x) an jedem Punkt x des Definitionsbereichs. Denn sonst könnte man mit einer gezielt konzentrierten Beule $(\pm \epsilon (x))$ den Wert des Funktionals in beide Richtungen verändern, Widerspruch. Die Variation ist für Funktionale W{} das Analoge zu den ersten Ableitungen für Funktionen $f:\mathbb {R} ^{n}\to \mathbb {R} ,$ denn solche erzeugen die linearen Näherungen

$f(x+\epsilon )=f(\{x_{i}+\epsilon _{i}\})\simeq f(x)+\sum _{k}\partial _{k}f(x)\cdot \epsilon _{k}$

Die explizite Formel der Variation liefert die gesuchte Differenzialgleichung, wenn man sie zu Null macht. Ohne Beweis und mit der Notation $\partial _{i}f(x)=(\partial f/\partial x_{i})$ hier aufgeschrieben, wird der Integrand für Variationen hergeleitet mit partiellen Integrationen und Annahmen über verschwindende Randwerte. Man beachte das Minuszeichen bei ungerader Zahl der Ableitungen:

$W\{f\}=\int d^{n}x\;L(x,\;f(x),\;\partial _{i}f(x),\;\partial _{i}\partial _{k}f(x),\;...)\quad \Rightarrow$

${\frac {\delta W}{\delta f(x)}}={\frac {\partial L}{\partial f(x)}}-\sum _{i}{\frac {\partial }{\partial x_{i}}}{\Big (}{\frac {\partial L}{\partial (\partial _{i}f(x))}}{\Big )}+\sum _{i,k}{\frac {\partial ^{2}}{\partial x_{i}\partial x_{k}}}{\Big (}{\frac {\partial L}{\partial (\partial _{i}\partial _{k}f(x))}}{\Big )}\;-...$

Unser Beispiel zum Aufwärmen macht nur Terme vom Typ der Einfachsumme,

$L=(\mu /2)\cdot {\dot {w}}(x,t)^{2}-(F/2)\cdot w'(x,t)^{2}$

Variationsgleichung für extremale Funktionen w(x,t):

${\frac {\partial L}{\partial {\dot {w}}}}=(\mu /2)\cdot (2{\dot {w}});\;{\frac {\partial L}{\partial w'}}=(-F/2)\cdot (2w')$

$0=-{\frac {\mu }{2}}\cdot {\frac {\partial }{\partial t}}(2{\dot {w}})+{\frac {F}{2}}\cdot {\frac {\partial }{\partial x}}(2w')$

$0=-\mu {\ddot {w}}(x,t)+Fw''(x,t)\;\Rightarrow \;\mu {\ddot {w}}=Fw''$

Nach all dem Aufwand, die Wellengleichung! Sie ist linear, weil die Lagrange-Dichte nur aus Termen von zweiter Potenz in den Funktionsargumenten besteht, eine quadratische Form also. Die linearen Maxwell-Gleichungen des Elektromagnetismus entstammen genauso einer quadratischen Lagrange-Dichte. Auch die linearen Schrödinger-Gleichungen und Dirac-Gleichungen der Quantenmechanik.

Ein erprobter Vorteil des Lagrange-Formalismus besteht darin, dass man nicht mühsam mit diversen Kräfteparallelogrammen hantieren muss, um an die Bewegungsgleichungen zu kommen. Sondern man kümmert sich nur um die oft einfacher herzuleitenden Energiesummen des Systems. Weiter unten bei den biegesteifen Saiten und Balken zahlt es sich aus.

Rechnen wir aus, was an Stelle des Terms $Fw''(x,t)$ im nichtlinearen Modell als Wellengleichung für große Auslenkungen auftritt. Mit dem definierten

$u={\sqrt {(1+w'^{2})}},$

$-{\frac {\partial L}{\partial w'}}={\frac {\partial V}{\partial w'}}={\frac {\partial V}{\partial u}}{\frac {\partial u}{\partial w'}}={\frac {F}{2\delta }}[(u-1)(1+\delta )+(u-1)+\delta (u+1)]{\Big [}{\frac {w'}{u}}{\Big ]}={\frac {F}{\delta }}[u-1+u\delta ]{\frac {w'}{u}}=$

$={\frac {F}{\delta }}(1+\delta -(1/u))w'=Fw'+{\frac {F}{\delta }}(1-(1/u))w'.$

Von diesem Ausdruck wird die Ableitung nach x gebraucht. Der erste Term ohne 'delta' macht den bekannten linearen Anteil. Im anderen Term ist ein Faktor (w'/u) abzuleiten:

$u'={\frac {w'\cdot w''}{u}}\;\Rightarrow \;{\frac {\partial }{\partial x}}{\frac {w'}{u}}={\frac {(w''u-u'w')}{u^{2}}}={\frac {w''}{u}}\cdot (1-(w'/u)^{2})={\frac {w''}{u^{3}}}$

$(-\partial /\partial x)({\frac {\partial L}{\partial w'}})=Fw''+{\frac {F}{\delta }}w''(1-u^{-3})$

Damit bleibt eine im Argument w nichtlineare Wellengleichung

$\mu {\ddot {w}}=Fw''{\Big [}1+{\big [}1-(1+w'^{2})^{-3/2}{\big ]}/\delta {\Big ]}.$

Sie hat keine einfachen Lösungen vom Typ $w=\sin(\omega t)\cdot g(x)$ mehr.

Mit der Approximation $(1+z)^{-3/2}\simeq (1-3z/2)$ wird die Gleichung zu:

$\mu {\ddot {w}}=Fw''{\big [}1+(3w'^{2})/(2\delta ){\big ]}$

Sobald die Steigungen einer weit ausgelenkten Saite dem Dehnungsfaktor δ zu nahekommen, wächst die effektive Kraft. Die Frequenz einer Schwingung wird im nichtlinearen Bereich höher als bei kleinen Amplituden.

Simulation.

Der nichtlineare Effekt kann numerisch an der vollen Gleichung durchprobiert werden. In etwa 60 Zeilen Python teilt man eine Saite (1 m lang) in 50 Punkte, approximiert die erste und zweite Ortsableitung der Auslenkungen w(x,t) als endliche Differenzen, teilt das Zeitintervall in 1000 Punkte und integriert Schritt für Schritt mit der Runge-Kutta-Formel vierter Ordnung. Die Anfangsbedingung w(x,t=0) ist eine Sinus-Halbwelle mit maximalen Amplituden von 1 mm bis 10 mm. Die Dehnung der ruhend gespannten Saite sei δ=0,001. Die Halbperiode zwischen zwei Nulldurchgängen der Saite wird "gestoppt", mit Parabel-Interpolation um die Minima von Amplitudenquadraten herum. Denn im nichtlinearen Modell behalten die Orte der Saite nicht exakt die gleiche Phase. Bei mehr als 1 mm Auslenkung wächst die Grundfrequenz gewaltig mit der Amplitude, verglichen mit der Frequenz des linearen Modells. Hier als Ergebnis die Faktoren des Frequenzanstiegs, nach wenigen Sekunden Laufzeit auf dem Raspberry Pi. Anhang: Script...

Amplitude	1 mm	2mm	5 mm	10mm
Faktor	1,0014	1,0055	1,0335	1,1235

Die longitudinale Welle im elastischen Stab

Der Stab hat die horizontalen Ortskoordinaten x von 0 bis zur Länge L, die Querschnittsfläche A und eine Materialdichte ρ. Aus der Ruhestellung sei der Querschnitt bei x um den kleinen Betrag w(x) horizontal verschoben. Ein kleines Segment von x bis x+Δx ist also komprimiert oder gestreckt, um dem Wert w(x+Δx)-w(x). Die potenzielle Energie dieser Dimensionsänderung ist $(EA/2)\cdot (w(x+\Delta x)-w(x))^{2}/(\Delta x)$, ganz wie bei der vorgestreckten Saite. Für einen kleinen Abstand ist das einfach $(EA/2)\cdot w'(x)^{2}\;\Delta x$. Die gesamte potenzielle Energie im Stab hat also eine lineare Dichte und beträgt zum Zeitpunkt t:

$V(t)=\int _{0}^{L}(EA/2)\cdot w'(x,t)^{2}\;dx.$

Soll die Verschiebung mechanisch vibrieren, tut sie es mit der kinetische Energie eines Segments $(\rho A\;\Delta x)\cdot {\dot {w}}(x,t)^{2}/2$. Die Lagrange-Dichte lautet daher:

$L(w,w',{\dot {w}})=(A/2)\cdot (\rho {\dot {w}}^{2}-Ew'^{2}).$

Die Euler-Lagrange-Gleichung ergibt wie oben eine Wellengleichung, wo hier der Querschnitt A herausfällt:

${\ddot {w}}=c^{2}w''\quad$ mit der Wellengeschwindigkeit $c={\sqrt {E/\rho }}.$

Angenommen, eine Art Stahl hat $\rho =8\;g/cm^{3}=8000\;kg/m^{3};\quad E=200\;GPa.$

Es folgt $c={\sqrt {200E9/8E3}}=1000\cdot {\sqrt {100/4}}=5\;km/s$, ein Vielfaches der Schallgeschwindigkeit in der Luft, wo der nur einige hundert Meter pro Sekunde entlang kriecht.

Hat ein frei aufgehängter Stab ein Spektrum von Eigenschwingungen? Die Randbedingungen sind nicht wie bei der Saite w(0,t)=w(L,t)=0. Für die Enden sei aber angenommen, dass die letzte Atomschicht keine zeitabhängige Streckung zum Rest erfährt. Diese Randbedingungen besagen, dass die erste Ortsableitung der Auslenkung verschwindet: w'(0,t)=w'(L,0)=0.

Zur Begründung: Innen im Volumen spürt die Atomlage des Querschnitts bei x eine Kraft, weil von rechts mit F_R= EA (w(x+Δx)-w(x))/(Δx) gezogen wird und von links in Gegenrichtung mit F_L= EA (w(x)-w(x-Δx))/(Δx). Die Differenz führt auf die zweite Ortsableitung von w() in der Bewegungsgleichung. An den Enden gibt es keine Ursache für eine Nettokraft, und die einseitige Nullkraft kann nur bedeuten, die erste Ableitung von w() muss dort wegfallen.

Weil die Wellengleichung linear-homogen mit konstanten Koeffizienten ist, erfüllen mit w(x,t) auch w'(x,t) und alle höheren gemischten partiellen Ableitungen dieselbe Gleichung. Denn wie die Analysis lehrt, sind unter milden Bedingungen die partiellen Ableitungen nach verschiedenen Variablen vertauschbar. Also gilt für die Funktion w'() des longitudinal schwingenden Stabes das, was für die Funktion w() der transversal schwingenden Saite herauskam. Es gibt eine Reihe von harmonischen Frequenzen f_n = n c/(2L). Nur sind die Bäuche der Saiten-Auslenkung die Knoten der Streckung des Stabes und umgekehrt. Bei der Grundschwingung ruht der Stab in der Mitte und die Enden zittern in Gegenphase. Der Stab wird periodisch gestreckt und gepresst.

Mit L=1 m und c=5 km/s ist die Grundfrequenz 2,5 kHz im Bereich der Piccoloflöte.

Kurze Stäbe (L= 10 cm) haben die Resonanz auf Ultraschall-Frequenzen.

Gebogene Stäbe und Saiten

Alte Instrumente, Cembalos, Clavichorde, erste Hammerklaviere, hatten dünne Saiten mit mäßigen Zugkräften; sie konnten ganz aus Holz gefertigt werden. Neuere Konzertflügel aber sollten immer lauter werden und ihre Töne länger anhalten. Das ging einher mit dickeren Saiten und höherer Anspannung, entsprechend der Formel für die Schwingungsfrequenz. Im Konzertflügel hält ein solider Metallrahmen die Saitenspannung aus, die insgesamt 20 Tonnen ausmachen kann. Die strenge Physiklehrerin wird uns tadeln, Tonnen bezeichnen eine Masse. Gemeint ist die Zugkraft in Newton, Masse mal Erdbeschleunigung, die ein hängendes Gewicht von 20 Tonnen erzeugt. Bei den Experimenten am historischen Monochord spannen tatsächlich angehängte Gewichte die Saite.

Welches Problem bereitet die Zusatzkraft, die zum Verbiegen in den schwingenden dicken Saiten angreift? Nicht tragisch und leicht zu kompensieren ist, dass der Grundton bei gleicher Zugkraft etwas höher liegt. Jedoch weicht das Spektrum der Obertöne nun etwas ab von den ganzen Vielfachen der Grundfrequenz. Und das prozentual umso mehr, je höher die Ordnungszahl der Harmonischen. Die Pianinos oder aufrechten Klaviere für den Hausgebrauch haben das größere Problem, kürzere und dickere Saiten.

Die nächste Rechenübung besteht darin, dies theoretisch zu begreifen und numerisch auszurechnen.

Die Biegung eines Stabs wird durch einen Krümmungsradius r charakterisiert. Der Stab hat die Querschnittsfläche A und die Länge L. Denkt man sich den elastischen Stab aus parallelen Fasern aufgebaut, dann werden die Fasern auf der Innenseite zusammengedrückt und die auf der Außenseite gedehnt. Mitten durch den Querschnitt geht eine Linie von neutralen Fasern, deren Länge sich nicht ändert. Die gespeicherte potenzielle Energie soll mit Hilfe des Elastizitätsmoduls E als Integral über eine Dichte formuliert werden.

Ein mikroskopisches Stück des Stabes soll die x-Achse zur Tangente haben, und die Koordinaten y,z streichen über den Querschnitt. Mit dem Radius r sei der Stab in der xz-Ebene nach unten gebogen. Die neutralen Fasern in der Mitte machen eine winzige Winkeländerung φ und haben die Bogenlänge Δx=rφ; sie liegen bei (x,y,z)=(0,0,0).

Die Fasern darunter sind gestaucht zu (r+z)φ mit z<0, solche darüber sind gestreckt zu einer Länge (r+z)φ mit positivem z. Die gespeicherte elastische Energie der verformten Faser der Dimension $\Delta x\times \Delta y\times \Delta z$ ist:

$(E\;\Delta y\;\Delta z)/2\cdot {\big [}((r+z)\phi -r\phi )/(r\phi ){\big ]}^{2}\cdot \Delta x,$

Hierbei ergibt die eckige Klammer einfach den Faktor (z/r)².

Die Energie der ganzen Scheibe mit Querschnittsfläche A ist gleich $(EI)/(2r^{2})\cdot \Delta x,$, wo I ein Integral über die Fläche A ist:

$I:=\iint _{A}z^{2}\;dy\;dz.\quad$ Es ist ein Moment dieser Fläche.

Wenn die neutrale Mittelfaser des Stabs nun längs einer Kurve z=w(x) verläuft, dann gehört zum horizontalen Stück Δx ein Stück Länge $\Delta l={\sqrt {(1+w'(x)^{2})}}\;\Delta x$ und die Energie der elastischen Verbiegung $(EI)/(2r\{w(x)\}^{2})\cdot {\sqrt {(1+w'^{2})}}\;\Delta x$. Dabei bezeichnet r{w} den Krümmungsradius, an den sich die Kurve im Punkt (x,w(x)) anschmiegt.

Zur schludrigen Ausrechnung des Radius nehmen wir an, w(x) sei ein Kreisbogen

$(w(x)-w_{0})^{2}+(x-x_{0})^{2}=r^{2},\quad$

Und wir nehmen zwei Ableitungen dieser Formel nach x vor:

$2(w-w_{0})w'+2(x-x_{0})=0;\quad (w-w_{0})w''+w'^{2}+1=0$

Daher:

$\quad (w-w_{0})=-(1+w'^{2})/w'';\quad (x-x_{0})=-w'(w-w_{0})$

Eingesetzt in die erste Gleichung folgt eine Formel für r²:

$r^{2}=[(+w'^{2})/w'']^{2}+w'^{2}[(1+w'^{2})/w'']^{2}=(1+w'^{2})^{3}/w''^{2}.$

Damit bekommen wir die Energie der Verbiegung längs der Kurve w(x) als Integral

$V=\int _{0}^{L}(EI/2)\cdot w''(x)^{2}\cdot (1+w'(x)^{2})^{-5/2}\;dx.$

Der inverse Krümmungsradius wird manchmal als die Krümmung der Kurve y=w(x) definiert:

$\kappa \{w\}(x)=|w''(x)|\cdot (1+w'(x)^{2})^{-3/2}.$

Sei nun der Stab als dicke Saite an beiden Enden festgeklemmt, und zwar nicht nur sei w(x)=w(L)=0, sondern auch die Richtung an beiden Enden sei als horizontal vorgegeben: w'(0)=w'(L)=0. Denn da der Stab sich der Biegung widersetzt, kann er nicht mit einem scharfen Knickwinkel anfangen und enden. Die potenzielle Energie ist die Summe aus der für die Biegungen und für die Streckungen, wenn der Stab mit einer Zugkraft F vorgespannt wird. Mit der schon benutzten Energiedichte der Streckung und dem Dehnungsfaktor δ=F/(EA) ergibt sich:

$V=\int _{0}^{L}\left[{\frac {EI}{2}}\cdot w''(x)^{2}\cdot (1+w'(x)^{2})^{-5/2}+{\frac {F}{2}}[w'(x)^{2}+(u-1)^{2}/\delta ]\right]\;dx;\quad u:={\sqrt {(1+w'(x)^{2})}}.$

In der linearen Näherung wird $w'^{2}\ll \delta \ll 1$ angenommen und es gilt:

$V\simeq \int _{0}^{L}\left[(EI/2)\cdot w''(x)^{2}+(F/2)\cdot w'(x)^{2}\right]\;dx.$

Kurz noch das Moment I für einen runden Draht der Schnittfläche A=π(d/2)² berechnen. Trick: man nutze die Symmetrie aus.

$I=\iint _{A}z^{2}\;dy\;dz={\frac {1}{2}}\iint _{A}(y^{2}+z^{2})\;dy\;dz={\frac {1}{2}}\int _{0}^{2\pi }d\phi \int _{0}^{d/2}(r\;dr)\;r^{2}={\frac {\pi \cdot (d/2)^{4}}{4}}={\frac {A^{2}}{4\pi }}.$

Stehende Wellen auf biegesteifen Saiten

Die Lagrange-Dichte der Saite bei kleiner Auslenkung, mit Krümmungs-Energie:

$L({\dot {w}},w',w'')=(\mu /2){\dot {w}}^{2}-(F/2)w'^{2}-(EI/2)w''^{2}$

Die Bewegungsgleichung wird streng nach Rezept angefertigt. Die Terme mit der ersten Ableitung enthalten mit ihren Minuszeichen die Faktoren

$(-\partial /\partial t){\dot {w}}=-{\ddot {w}}\;{\text{ und }}\;(\partial /\partial x)w'=w''$.

Der Term mit der zweiten Ableitung ergibt eine vierte Ableitung ohne Vorzeichenwechsel; der relevante Faktor wird $-(\partial ^{2}/\partial x^{2})w''=-w''''$. Ergebnis:

$\mu {\ddot {w}}=F\cdot w''-EI\cdot w'''';\;$ Randbedingung w(0)=w(L)=w'(0)=w'(L)=0.

Die Gleichung ist linear für die Funktionen w(x,t). Der Separations-Ansatz $w(x,t)=f(x)\cdot \sin(\omega t)$ macht die Differenzialgleichung:

$\mu \omega ^{2}\;f+F\;f''-EI\;f''''=0.$

Solange der Faktor (EI) klein ist verglichen mit (F), ist dies als eine Verformung der gewöhnliche Wellengleichung anzusehen. Wir erwarten ein diskretes Spektrum von möglichen Frequenzen {ω_k}, die sich womöglich numerisch durch Variation aus den bekannten Werten für die einfache Saitengleichung bestimmen lassen.

Zuerst machen wir die Gleichung dimensionslos mit z:=(x/L); g(z):=f(x). Immer eine gute Idee, bevor irgendwas mit Numerik angegriffen werden soll. Auch legen wir die Saite symmetrisch um den Nullpunkt, so dass die neue Variable z das Intervall [-0,5...0,5] belegt. Mit g'=(dg/dz)=(df/dx)(dx/dz)=Lf' lautet die Gleichung in g:

$\mu \omega ^{2}\;g+(F/L^{2})\;g''-(EI/L^{4})\;g''''=0.$

Mit zwei dimensionslosen Konstanten $k^{2}:=\mu \omega ^{2}L^{2}/F;\;q:=EI/(L^{2}F):$

$k^{2}\;g+g''-q\;g''''=0$

Die Mathematik lehrt, dass jede Lösung einer linear-homogenen gewöhnlichen Differenzialgleichung der Ordnung N eine Linearkombination aus N Basis-Lösungen ist. Für die vierte Ordnung hier wäre folgende Basis aus vier Funktionen $\{\sin(az),\cos(az),\exp(bz),\exp(-bz)\}$ zu probieren, da deren zweite Ableitungen einfach konstante Faktoren absondern. Besser noch, die Kandidaten werden wegen der Symmetrie auf gerade und ungerade Funktionen

$g_{1}(z)=g_{1}(-z)=u\cos(az)+v\cosh(bz)\;\quad g_{2}(z)=-g_{2}(-z)=u\sin(az)+v\sinh(bz)$

beschränkt. Die hyperbolischen Funktionen sind:

$\sinh(x)=(\exp(x)-\exp(-x))/2;\quad (d/dx)\sinh(x)=\cosh(x)$

$\cosh(x)=(\exp(x)+\exp(-x))/2;\quad (d/dx)\cosh(x)=\sinh(x)$

Wegen der Linearität kann u=1 gesetzt werden. Es gibt drei freie Variablen (a,b,v) aber vier Gleichungen, nämlich zwei Koeffizienten aus der Differenzialgleichung und zwei Randbedingungen bei z=0,5. (Die bei z=-0,5 sind redundant wegen der Symmetrie.) Die Rechenaufgabe wird deshalb machbar, weil k² die vierte Unbekannte ist und die gesuchte Eigenfrequenz enthält. Dagegen ist q konstant durch Geometrie und Material der Saite vorgegeben. Daher diese Gleichungen für $(k^{2},a,b,v)$ mit dem Kandidaten g₁:

$0=k^{2}(\cos(az)+v\cosh(bz)+(-a^{2}\cos(az)+vb^{2}\cosh(bz))-q(a^{4}\cos(az)+vb^{4}\cosh(bz))$

$k^{2}-a^{2}-qa^{4}=0;\quad k^{2}+b^{2}-qb^{4}=0$

$\cos(a/2)+v\cosh(b/2)=0;\quad -a\sin(a/2)+vb\sinh(b/2)=0.$

Mit dem Kandidaten g₂ wechseln nur die Randbedingungen:

$\sin(a/2)+v\sinh(b/2)=0;\quad a\cos(a/2)+vb\cosh(b/2)=0.$

Numerik: Wird der Wert k erraten, dann folgen a,b aus den quadratischen Gleichungen (für a²,b²) und v ergibt sich aus der dritten Gleichung. Die vierte und letzte sollte dann Null ergeben; sie wird als Fehlergleichung genutzt. Nach dem Einklammern der Nullstelle wird das gesuchte k so lange interpoliert, bis dieser Fehler unter Epsilon fällt.

Die Anfangswerte auf der Suche nach k kommen aus der Randbedingung der einfachen Wellen, kL=k=nπ mit positiven ganzen Zahlen. Ungerade Zahlen starten die Suche nach Funktionen mit {cos(),cosh()}, die hier auf dem symmetrischen Intervall die geraden Funktionen sind.

Das folgende Skript rechnet ein naives Beispiel durch. Ein Draht mit 1 Millimeter Durchmesser und 1 Meter Länge ist an beiden Enden befestigt und mit 25 Newton gespannt. Elastizitätsmodul 200 GPa. Die Streckgrenze sei über 200 MPa, so dass das Material nicht leidet. Die Eigenfrequenzen N=1 bis 10 werden gezeigt, und zwar geteilt durch N mal die Grundfrequenz von etwa 32 Hz. Die Inharmonizität der biegesteifen Saite spreizt in der Tat die Oktaven. Die Harmonischen N=2, N=4, N=8 sind zu groß um 0,6%, 2,9% und sogar 11,5%. Das Beispiel ist eine Karikatur der Klaviersaite im Bass. Denn in Wahrheit haben solche Saiten dünnere Drahtkerne, umwickelt mit weichen Metall. Daher ist die Spreizung der Oktaven weniger krass, aber auf jeden Fall beim Stimmen zu berücksichtigen.

# -*- coding: utf-8 -*-
from math import exp, log, sqrt, sin, cos, pi, sinh, cosh

def testdraht() : # 1m Stahl, Durchm 1mm, Zug~ Gewicht 2.5 kg. freq 30 Hz.
  length= 1.0; radius=0.5e-3; rho=8000; elast=200e9; force=25 #Newton
  area= pi*radius*radius; iarea= area*radius*radius/4.0; mu= rho*area
  speed= sqrt(force/mu); period= length/speed; fbase= speed/(2*length)
  z= period/length/length; q= (elast*iarea/mu)*z*z
  print('Parameter q='+str(q)+' Frequenz(Hz)='+str(fbase))
  ##   fuers Timoshenko-Modell noch 2 Parameter r,s
  gmodul=10e9 #  uebertrieben, standard 80E9
  kappa=0.9; r=kappa*gmodul*area/force
  z= period/length/length; q= (elast*iarea/mu)*z*z
  s=  area*length*length/iarea # = 4*(length/radius)^2
  return q,r,s, fbase

def fehler(nk,q,sign) : # Euler-Modell. wird nk richtig geraten, muss err=0 sein
  k2=nk*nk; z=1.0/(2*q); p=sqrt(k2/q+z*z) 
  a2=p-z; a=sqrt(a2); b2=p+z; b=sqrt(b2) 
  if sign>0 : v=-cos(a/2)/cosh(b/2); err=-a*sin(a/2)+v*b*sinh(b/2)
  if sign<0 : v=-sin(a/2)/sinh(b/2); err= a*cos(a/2)+v*b*cosh(b/2)
  return err

def drahtspektrum(q,f0, timo=None) : # timo = data for Timoshenko model
  def interpol(xa,xb, ya,yb) : # interpolate x at y=0, given ya*yb < 0
    return xa + ya*(xa-xb)/(yb-ya+0.0)
  def tfehler(k,pars) : z,data=timofehler(k*k,pars); return z # timoshenko
  sign=1; std=True; kval=[0.0]*11; epsi=1e-10; kref=pi # propto ref frequency
  if timo is not None : q,r,s,pars= tuple(timo); std=False; pars[7]=sign
  for n in range(1,11) :
    if n==1 :   k=kref; dk=0.1*k
    elif n==2 : dk=kval[n-1]; k= kval[n-1]+dk; dk=0.1*dk
    else :      dk=kval[n-1]-kval[n-2]; k=kval[n-1]+dk; dk=0.1*dk
    ka=k; bracket=False; ok= False; count=0
    ya=fehler(ka,q,sign) if std else tfehler(ka,pars)
    while (not bracket)and(count<25) : # Klammer um Nullstelle
      kb=k+dk
      yb= fehler(kb,q,sign) if std else tfehler(kb,pars)
      bracket=(yb*ya)<0 
      if not bracket : ka=kb;ya=yb; count+=1
    if bracket : # fast in Euler-Bernoulli, perfect straight line 
      emax=epsi*(abs(ya)+abs(yb)); # print('bracket '+str([count,ya,yb]))
      while (not ok) and (count<50) :
        k=interpol(ka,kb,ya,yb)
        y= fehler(k,q,sign) if std else tfehler(k,pars) 
        if (ya*y)>0 : ya=y;ka=k
        else : yb=y;kb=k
        ok= (abs(y)<emax); count +=1
      if ok :
        kval[n]= k; sign=-sign
        print('N,Iterationen,Freq='+str([n,count,k/n/kval[1]]))
      else : print('exit not ok n='+str(n)+' count='+str(count)); exit()
      if ok and (not std) : pars[7]=sign
    else : print('exit no bracket n='+str(n)+' count='+str(count)); exit()
  print('Grundton mit Biegekraft(Hz)='+ str(f0*kval[1]/kref))

q,r,s,f0= testdraht()
drahtspektrum(q,f0)

Balkenmodell von Timoshenko

Im Euler-Bernoulli-Modell des gebogenen Drahts oder Stabs gibt es den einen Freiheitsgrad w(x) der vertikalen Auslenkung. Die Querschnittsfläche hat an jeder Stelle eine Normale, die in Richtung w'(x) zeigt. Elastische Energie wird durch Ziehen oder Drücken in Drahtrichtung eingespeist. Dafür ist der Elastizitätsmodul zuständig. Timoshenko lässt die Ausrichtung der gebogenen Flächen von der Drahtrichtung abweichen und der Winkelfehler wird eine zweite dynamische Feldvariable. Die Dehnung oder Stauchung einer Faser im Draht agiert mit dem Parameter E. Aber eine Verschiebung benachbarter Fasern zueinander wird jetzt erlaubt und ein Teil der Energiezufuhr fließt da hinein. Der muss vom Schubmodul G abhängen.

Anschaulich wollen die innen liegenden Fasern beim Biegen nicht gern auf die naive Geometrie verkürzt werden, die außen liegenden nicht voll gestreckt. Die Fasern wollen andererseits nicht beliebig gegeneinander verschoben werden. Wegen des Kompromisses neigen sich die Endflächen eines Bogenstücks nicht um den vollen Winkel der Biegung, sondern weniger.

Die Geometrie des verbogenen Drahts wird mit zwei Funktionen beschrieben:

• w(x), die vertikale Auslenkung des Mittelpunkts und
• φ(x), der Winkel des gezerrten und gedrehten Querschnitts zur vertikalen Achse.

Wäre der Verschiebe-Effekt ausgeschaltet, also der Schubmodul Unendlich, dann gälte: $\tan \phi (x)=w'(x)$.

Die Schritte zum zweidimensionalen Modell:

• Die Kinematik der Verformung im Draht $(u_{x},u_{z})={\text{Funktion}}(t,x,y,z),$
• Die Dichte der potenziellen Energie nach dem allgemeinen Hooke-Schema,
• Die Lagrange-Dichte und die Euler-Lagrange-Gleichungen,
• Dynamische (und statische) Lösungen der Gleichungen.

Geometrie der Verformung nach Euler (rot) und Timoshenko (blau)

In der Ruhestellung hat der Draht oder Balken die Koordinaten (x,y,z), wobei y=z=0 die Achse des zentralen Fadens ist. Die z-Achse ist die Vertikale. Die Auslenkung des Punktes (x,y,z) mit den Funktionen w(), φ() ist

$(u_{x},u_{y},u_{z})(x,y,z)=(-z\sin \phi (x),\;0,\;w(x)+z\cos \phi (x)-z)$

Nur die lineare Näherung in Variablen w,φ zieht ins Modell ein:

$(u_{x},u_{z})\simeq (-z\phi (x),w(x))$

Mit vernachlässigter Querkontraktion folgen dann die Epsilon-Sigma-Tensoren und die Energiedichte.

${\begin{array}{ll}\epsilon _{xx}=(\partial u_{x}/\partial x)=-z\phi '(x)&\sigma _{xx}=E\epsilon _{xx}\\\epsilon _{xz}=({\frac {\partial u_{x}}{\partial z}}+{\frac {\partial u_{z}}{\partial x}})/2=(-\phi (x)+w'(x))/2&\sigma _{xz}=(2G)\epsilon _{xz}\\\epsilon _{zz}=(\partial u_{z}/\partial z)=0&\sigma _{zz}=0\end{array}}$

$U({\vec {x}})={\tfrac {1}{2}}\sum \nolimits _{i,k}\epsilon _{ik}({\vec {x}},{\vec {u}}())\;\sigma _{ik}({\vec {x}},{\vec {u}}())$

$U(x,y,z)=(\epsilon _{xx}\sigma _{xx}+\epsilon _{zz}\sigma _{zz})/2+\epsilon _{xz}\sigma _{xz}=E(z\phi ')^{2}/2+G(-\phi +w')^{2}/2$

Daraus macht man eine eindimensionale Energiedichte ${\bar {U}}(x)$, indem über y,z, also den Querschnitt des Balkens, integriert wird. Das Quadrat von z im ersten Term bringt das Moment I wie bei Euler-Bernoulli, für den zweiten kommt die Fläche A.

${\bar {U}}(x)=[EI(\phi ')^{2}+GA(-\phi +w')^{2}]/2$

Die Dichte der kinetischen Energie, mit der Massendichte ρ,

$T(x,y,z)=(\rho /2)({\dot {u}}_{x}^{2}+{\dot {u}}_{z}^{2})=(\rho /2)(z{\dot {\phi }})^{2}+{\dot {w}}^{2})$

wird gleichermaßen zu einer Linien-Dichte integriert:

${\bar {T}}(x)=(\rho /2)(I{\dot {\phi }}^{2}+A{\dot {w}}^{2})$

Die 1D-Lagrangedichte ${\bar {L}}={\bar {T}}-{\bar {U}}$ wird nun mit den zwei dynamischen Feldern w(t,x), φ(t,x) und den zwei unabhängigen Variablen t,x variiert. Nullsetzen der Variation ergibt zwei Bewegungsgleichungen.

${\begin{array}{l}0=(-\partial {\bar {U}}/\partial \phi )+\partial _{x}(\partial {\bar {U}}/\partial \phi ')-\partial _{t}(\partial {\bar {T}}/\partial {\dot {\phi }})\\0=(-\partial {\bar {U}}/\partial w)+\partial _{x}(\partial {\bar {U}}/\partial w')-\partial _{t}(\partial {\bar {T}}/\partial {\dot {w}})\\\rho I\;{\ddot {\phi }}=GA(-\phi +w')+EI\;\phi ''=GA(w'-\phi )+EI\;\phi ''\\\rho A\;{\ddot {w}}=GA\;\partial _{x}(-\phi +w')=GA(w''-\phi ')\end{array}}$

Um manche Approximationen und Nachlässigkeiten zu kompensieren, wurde ein Faktor eingeführt, der den Schubmodul verändert.

$G\;\to \;\kappa \;G;\quad \kappa =5/6\;$ für rechteckigen Querschnitt.

Bei einem runden Querschnitt zitiert Wikipedia (en:Timoshenko_beam_theory) die folgende Formel. In den Timoshenko-Koeffizienten fließt die Querkontraktion über die Poissonzahl ν ein.

$\kappa ={\frac {6+6\nu }{7+6\nu }};\quad \nu =\{0.3,0.5\}\to \kappa =\{0.886,0.9\}.$

Für die Praxis fehlen den Gleichungen noch zwei Terme: eine vertikale Belastung, die von der Position x abhängt, und eine horizontale Vorspannung. Die erste interessiert die Architektin, wenn tragende Balken sich durchbiegen, die zweite den Instrumentenbauer bei dicken Saiten. (Als ob Architekten und Klaviermacherinnen den Schrieb bis hierhin lesen würden.)

Sei q(t,x) die vertikale Kraft pro Länge, in Richtung der positiven z-Achse. Sie kann in Ort und Zeit beliebig variieren. Sie addiert zur Linien-Dichte der potenziellen Energie den Term (-w(t,x)q(t,x)). Eine Vorspannung mit der Kraft F addiert in erster Näherung den Term zur eindimensionalen Energiedichte, der schon bei Euler-Bernoulli benutzt wurde,

${\bar {U}}\to {\bar {U}}+{\tfrac {1}{2}}F\;{w'(t,x)}^{2}.$

Das führt direkt zu diesen Timoshenko-Gleichungen, gut genug für unsere Zwecke.

$\rho A\;{\ddot {w}}(t,x)=\kappa GA(w''(t,x)-\phi '(t,x))+q(t,x)+F\;w''(t,x),$

$\rho I\;{\ddot {\phi }}(t,x)=\kappa GA(w'(t,x)-\phi (t,x))+EI\;\phi ''(t,x).$

Eigenschwingungen dicker Saiten, numerisch

Wieder geht es um Eigenfrequenzen des linear-homogenen Systems. Man errät gerade und ungerade Lösungen mit dem x-Ursprung in der Mitte. Ansätze:

$w(-x)=+w(x):\quad w(t,x)=\sin(\omega t)\;[a\cos(px)+b\cosh(qx)];\quad \phi (t,x)=\sin(\omega t)\;[c\sin(px)+d\sinh(qx)].$

$w(-x)=-w(x):\quad w(t,x)=\sin(\omega t)\;[a\sin(px)+b\sinh(qx)];\quad \phi (t,x)=\sin(\omega t)\;[c\cos(px)+d\cosh(qx)].$

Wegen des Ausdrucks $(w'-\phi )$ wird das Paar mit entgegengesetzter Parität ('Geradheit') angesetzt. Bei Nulldurchgängen, alle w(t,x)=0 für gewisse t, sollte auch φ verschwinden. Daher dieselbe Phase im Zeitverhalten. Die Randbedingungen sagen, dass w(t,x) und φ(t,x) permanent Null sind an den Enden der Saite.

Bemerkung. Die Mathematik weiß ohne Raterei zu dem Ansatz zu kommen. Es gibt für zwei solcher Gleichungen zweiter Ordnung, ebenso für eine von vierter oder für vier von erster Ordnung, eine Basis mit der x-Abhängigkeit $\exp(k_{i}x)$, wo vier komplexe Eigenwerte k_i auflaufen. Die Symmetrien dieses Beispiels machen daraus eine schönere Basis

$\{\sin(fx),\cos(fx),\sinh(gx),\cosh(gx)\}.$

Sei L die Länge der Saite und z=x/L die dimensionslose Koordinate. Vor jedem Ausrechnen geht man über zu dimensionslosen Funktionen, v(t,z)=φ(t,x); u(t,z)=w(t,x)/L. Damit wird wie folgt ersetzt, wo die Striche bei u,v die Ableitungen nach z bezeichnen; $z\in [-{\tfrac {1}{2}},{\tfrac {1}{2}}].$

${\begin{array}{llll}{\ddot {w}}=-L\omega ^{2}u;&w=Lu;&w'=u';&w''=u''/L,\\{\ddot {\phi }}=-\omega ^{2}v;&\phi =v;&\phi '=v'/L;&\phi ''=v''/L^{2}.\end{array}}$

Das ergibt mit dem dimensionslosen 'Frequenzquadrat' $k^{2}:=\rho A\omega ^{2}L^{2}/F,$

${\begin{array}{ll}0=k^{2}u+r(u''-v')+u'';&u(z)=a\cos(fz)+b\cosh(gz);\\0=k^{2}v+sr(w'-v)+sq\;v'';&v(z)=c\sinh(fz)+d\sinh(gz);\\a\cos(f/2)+b\cosh(g/2)=0;&c\sin(f/2)+d\sinh(g/2)=0.\end{array}}$

Hier tauchen drei dimensionslose Parameter q,r,s auf:

$q=EI/(L^{2}F);\quad r=\kappa GA/F;\quad s=AL^{2}/I.$

Der Ansatz funktioniert und findet die Reihe der Eigenschwingungen. Die numerische Auswertung liefert die n-Tupel {a,b,c,d,f,g,k} und benutzt wie im obigen Skript eine Fehlerfunktion zur Nullstellensuche in Eigenfrequenz-Intervallen. Selbst bei übertrieben verschiebungsfreudigem Material, G= 10 Gigapascal, weicht die Frequenzliste nur in den dritten bis vierten Nachkommastellen von der des vorherigen Beispiels ab. Die Tendenz ist richtig: ein schwächerer inharmonischer Effekt, denn der Draht biegt sich etwas leichter.

Folgende Skript-Fetzen enthalten die Initialisierung und die Fehlerfunktion fürs Timoshenko-Beispiel.

def ti_system(x,par) : # nur zur Dokumentation, Timoshenko-Gleichungssystem
  #(1) Av+Bv"+Cu'+kv=0 =  A(cS+dT)+(kcS+kdT)+B(-cffS+cggT)+C(-afS+bgT)
  #(2) Du+Eu"+Fv'+ku=0 =  D(aK+bU)+(kaK+kbU)+E(-affK+bggU)+F(+cfK+dgU)
  # S,T,K,U= sin,sinh,cos,cosh. Paritaet= Vorzeichen letzter Term in y[0],y[2]
  a=0.001; b,c,d,f,g,k= tuple(x); y=[0]*6
  pa,pb,pc,pd,pe,pf,px,sign= tuple(par) # A B C D E F 0.5 1
  cs=cos(f*px); sn=sin(f*px); ch=cosh(g*px); sh=sinh(g*px) #Endpunktwerte
  f2=f*f; g2=g*g # Loesung gesucht, alle y[i]=0.
  y[0]=(pd+k-pe*f2)*a + sign*pf*f*c # coeffs K of (2)
  y[1]=(pd+k+pe*g2)*b + pf*g*d      # coeffs U of (2)
  y[2]=(pa+k-pb*f2)*c - sign*pc*f*a # coeffs S of (1)
  y[3]=(pa+k+pb*g2)*d + pc*g*b      # coeffs T of (1)
  y[4]=a*cs+b*ch # u(z)
  y[5]=c*sn+d*sh # v(z)

def timoshenko_init(q,r,s) :
  # 0 = k u + r(u''-v') +  u''  = ku+ Du+ Eu" + Fv'
  # 0 = k v + sr(u'-v) + sq v'' = kv+ Av+ Bv" + Cu'
  pa= -s*r; pb= s*q; pc=-pa  # fuelle Parameter A B C D E F
  pd= 0; pe= r+1; pf= -r
  a=0.001; px=0.5; sign=1 # a= freie Wellen-Amplitude (linear-homogene Gln.)
  pars=[pa,pb,pc,pd,pe,pf,px,sign]
  return a, pars

def timofehler(k,par) : # berechne f,g, b,c,d aus k. liefert err=0 wenn korrekt
  def quad(a,b,c) : z=b/(2.0*a); return -z+sqrt(z*z-c/a) # axx+bx+c=0
  a=0.001; pa,pb,pc,pd,pe,pf,px,sign= tuple(par) # A...F. a=1e-3, Amplitude 1mm 
  # (pd+k-pe*f2)*a = -sign*pf*f*c # multipliziere y[0]=0 mit y[1]=0
  # (pa+k-pb*f2)*c = sign*pc*f*a; -pc*pf*f2= (pd+k-pe*f2)*(pa+k-pb*f2) 
  # (pd+k)*(pa+k)-f2*(pe*(pa+k)+pb*(pd+k)+pc*pf)+ pe*pb*f2*f2 = 0
  f2=quad(pe*pb, -(pe*(pa+k)+pb*(pd+k))+pc*pf, (pd+k)*(pa+k)); f=sqrt(f2)
  # (pd+k+pe*g2)*b =-pf*g*d  # y[2]=0 mit y[3]=0
  # (pa+k+pb*g2)*d =-pc*g*b;  (pd+k+pe*g2)*(pa+k+pb*g2)=pc*pf*g2 
  g2=quad(pe*pb, pe*(pa+k)+pb*(pd+k)-pc*pf, (pd+k)*(pa+k)); g=sqrt(g2)
  cs=cos(f*px); sn=sin(f*px); ch=cosh(g*px); sh=sinh(g*px)
  c= (pd+k-pe*f2)*a/(-sign*pf*f)     # aus y[0] 
  if sign>0 : b=-a*cs/ch; d=-c*sn/sh # aus Randbedingungen y[4],y[5]
  else :      b=-a*sn/sh; d=-c*cs/ch 
  err= (pa+k+pb*g2)*d+ pc*g*b       # Gleichung y[3]=0, noch ungenutzt 
  return err, [b,c,d,f,g,k]

def draht2() :
  q,r,s,f0= testdraht()
  drahtspektrum(q,f0) # Euler-Bernoulli-Spektrum
  print(''); print('Timoshenko-Modell')
  a,pars= timoshenko_init(q,r,s)
  drahtspektrum(q,f0, timo=[q,r,s,pars])

draht2()

Perioden bei verriegelten Phasen

Wenn ein Klangerzeuger eine Reihe von Frequenzen abgibt, die genau ganzzahlige Vielfache einer Grundfrequenz sind, behalten diese Harmonischen alle dieselbe Phasenbeziehung. Alle Überlagerungen sind streng periodisch mit dem Zeitintervall T, gleich der reziproken Frequenz des Grundtons. Die freien Schwingungen einer steifen Saite nach dem Anschlagen oder Zupfen sind jedoch inharmonisch und der Zeitverlauf hat keine Periode. Die Nulldurchgänge der verschiedenen Obertöne laufen voneinander weg. Abgesehen von Saiten gibt es noch andere inharmonische Schwinger, so wie die zylindrischen oder konischen Rohre der Blasinstrumente, die wegen ihrer Querschnitts-Dimension und wegen der Effekte an den Enden keine streng harmonischen Spektren von Eigenfrequenzen aufweisen.

Ein wichtiges nichtlineares Phänomen kann all diesen Resonatoren doch ein streng ganzzahliges Obertonspektrum aufzwingen. Man nennt es die Verriegelung der Phasen. Die Wechselwirkung geschieht an der Kante, wo Flöten oder Orgelpfeifen angeblasen werden, oder an den Plättchen in den Mundstücken von Oboen, Klarinetten und so weiter. Sie bewirkt, dass einmal pro Grundschwingung das Signal sozusagen vereinheitlicht wird und eine wohldefinierte Periode beibehält. Auch die Streichinstrumente verriegeln durch den ständigen Kontakt des Bogens die Phasen der Teilschwingungen einer Saite, selbst wenn diese etwas dicker und daher inharmonisch ist.

Ein ganz langsamer mechanischer Schwinger, an dem der Effekt in Zeitlupe auftritt, ist die Glocke. Sie besteht aus zwei Pendeln, dem Körper mit einer größenbedingten Pendelfrequenz und dem Klöppel, dessen Eigenfrequenz niedriger liegen muss. Nach jedem Ausschlag der Glocke nach links oder rechts trifft sie auf dem Rückweg den nachzüglerischen Klöppel. Diesem wird ein Impuls versetzt, der ihn vorzeitig zur Umkehr zwingt. Beide Teile der Glocke werden dadurch synchronisiert. Ein Mensch oder ein Motor muss im Rhythmus Energie nachliefern.

Randbemerkung zur Harmonik der Glocken

Im akustischen Bereich ist die Glocke wegen ihrer recht komplizierten Form ein inharmonischer elastischer Schwinger mit vielen Eigenfrequenzen, die sich als stehende Wellen mit mehr oder weniger Knotenlinien auf der Glockenfläche ausbilden. Erst um 1500 herum war die Glockengießerei gut genug entwickelt. Man konnte die Bronzelegierungen und die Geometrie gezielt gestalten, damit ein gewollter Grundton und ein nicht allzu schräges Spektrum von Nebentönen herauskam. Es entstand die Gotische Moll-Oktavglocke. Später wurden im Barock, zuliebe einer geschwungenen Ästhetik (man ist versucht zu spotten einer Walt-Disney-Formgebung), gewisse Fortschritte zeitweise wieder geopfert.

Rückfall auf ein verwaschenes Spektrum.

Der empfundene Ton der Glocke heißt der Schlagton oder Nennton. Er ist nicht unbedingt dominierend in Spektrum; die nächstliegende Eigenfrequenz kann sogar von ihm abweichen. Er ist aber der Residualton, den sich das Ohr praktisch aus allen harmonisch verwandten Teiltönen zusammenbastelt. Die Moll-Oktavglocke hat dank eines ausgeklügelten Profils, auch Rippe genannt, das Spektrum unter Kontrolle. Es besteht aus einem Unterton (der Unteroktave), der Prime gleich dem Schlagton, der Mollterz, Quinte und Oktave. Als schwächere Obertöne kommen dazu noch Dezime, Undezime, Duodezime, Doppeloktave, Tripeloktave.

Der Beginn des gregorianischen Hymnus „Te Deum laudamus“ im C-Schlüssel mit den ersten drei Tönen e – g – a.

Te-Deum-Geläut auf den Tönen e, g und a der Kirche Mater Dolorosa in Berlin-Lankwitz

Ein bescheidenes Geläut aus drei Glocken stimmt zum Beispiel über dem tiefen Nennton seine anderen Glocken ab auf die kleine Terz und die Quarte, gewählt als Anfangsnoten des liturgischen Gesangs Te Deum laudamus.

Schallwellen in der Luft

Vorbemerkungen zur Wellenphysik in Gasen

Die Luft und die Gase sind auch irgendwie elastisch und transportieren Schwingungsformen praktisch unverfälscht über weite Strecken. Wie kommt das? Dieser Abschnitt versucht Erklärungen auf dem Hintergrund der Dynamik von Fluiden. Feste Gegenstände bestehen im naiven Modell aus Kügelchen, die mit Federn aneinander haften. Damit ist das Hookesche Gesetz plausibel. In Wirklichkeit verlangt die Analyse der chemischen Bindungen zwischen Atomen den ganzen Apparat der Quantenmechanik, um die Elastizität herzuleiten. Und was passiert bei Gasen und Flüssigkeiten? Das naive Modell sieht einen Haufen von ungebundenen Kügelchen, die ungeordnet durcheinander fliegen und immer wieder elastisch zusammenstoßen. Bei jedem Stoß gibt es Energieerhaltung.

Ein Stück Luft verfügt also über eine Zustandsgröße, die innere Energie U, Summe aller kinetischen Energien der Moleküle. Sie ändert sich nicht, solange kein Austausch mit der Umgebung erfolgt; auch die unvermeidliche Abstrahlung und der Einfang von Photonen seien mal kurz weggedacht. Die Statistische Mechanik hat es geschafft, die Energie U und andere Zustandsgrößen wie die Entropie S als Funktionen von wenigen makroskopisch messbaren Parametern des Gases herzuleiten. Wir brauchen etwas von der kinetischen Gastheorie, denn diese liegt dem Phänomen der Schallausbreitung zugrunde.

Die wichtigste Erkenntnis ist der Gleichverteilungssatz, er wird pompös auch das Äquipartitions-Theorem genannt. Im statistisch stationären Zustand (unter Vernachlässigung der Quantenmechanik) haben alle mechanischen Freiheitsgrade der Atome oder Moleküle im Mittel dieselbe kinetische Energie. Ein Gasatom hat drei Freiheitsgrade der Translation, ein Molekül dazu noch zwei bis drei relevante für Rotationen. Die mittlere Energie ist eine Zustandsgröße, die konventionell als die Absolute Temperatur T herhalten muss. Aus historisch-praktischen Gründen wird die Temperatur in Kelvin gemessen und die Energie in Joule. Die Konstante der Umrechnung ist nach Ludwig Boltzmann benannt, $\langle E\rangle ={\tfrac {1}{2}}k_{B}\cdot T$. Die statistische Verteilungsfunktion der Energien ist eine Exponentialfunktion. Die relative Häufigkeit, mit der ein System mit kinetischer Gesamtenergie E angetroffen wird, ist proportional zu $W(E)=e^{-\beta E}=\exp {\big (}-E/(k_{B}T){\big )}$ (Maxwell-Boltzmann-Verteilung).

Was ist die Dichte der Energien von Teilchenzuständen, über die aufsummiert wird? Die Energie ist eine Zufallsvariable auf dem Phasenraum mit 6N Koordinaten, den Paaren von (Ort,Impuls) aller Teilchen. Das Maß auf diesem Phasenraum ist $\prod (d^{3}q_{i}\;d^{3}p_{i}\;\exp(-\beta p_{i}^{2}/(2m))$, vor Normalisierung. Zum Erwartungswert der vollen kinetischen Energie $\sum \nolimits _{i}({\vec {p}}_{i}^{2}/2m)$ leistet dann jeder der 3N Freiheitsgrade den Beitrag $1/(2\beta )=k_{B}T/2$. Die Ausrechnung ist langweilige Routine und wir hier weggelassen.

In der Quantenmechanik funktionieren praktisch nur die Translations-Freiheitsgrade mit kontinuierlicher Energie von Null aufwärts. Jedoch die Freiheitsgrade von Rotation und Vibration der Moleküle bekommen diskrete Energieniveaus und sind bei niedriger Temperatur eingefroren. Die Äquipartition gilt nicht und deren Verteilung ist nicht mehr exponentiell. Erst wenn die Temperatur die Schwelle des jeweiligen Quantums übersteigt, taut der Freiheitsgrad auf und trägt bei zur Wärmekapazität. Das heißt, zu der Ableitung der mittleren Energie einer Molekül-Sorte als Funktion der Temperatur.

Adiabatische Prozesse

Ein ideales Gas sei in einem Würfel mit dem Volumen V eingesperrt, dessen elastischen Wände die aufprallenden Atome ohne Energieverlust abschmettern. Jeder Aufprall bedeutet einen Kraftstoß auf die Wand, nämlich Kraft mal Zeit gleich Änderung des Impulses = 2p, wo p den Impuls des ankommenden Atoms senkrecht zur Wand bezeichnet. Wir wollen den Druck P abschätzen, den das Gas auf die Wand ausübt, also die Kraft pro Flächeneinheit im Zeitmittel. Es ist die Summe aller Impulsüberträge, geteilt durch die Zeit der Messung. Wir haben $p^{2}/(2m)=k_{B}T/2;\;p=mv$. Die Dichte des Gases sei ρ Atome pro Volumeneinheit. Die Hälfte der Atome geht mit Geschwindigkeit v Richtung Wand, so dass eine Fläche A in der Zeit t die Anzahl von Stößen $\rho Avt/2$ mitbekommt. Kraft mal Zeit als Summe der Impulse daher

- $Ft=(2mv)\rho Avt/2=At\;\rho mv^{2}$.

- Wegen $(mv^{2}/2)=k_{B}T/2$ folgt $\;P=(F/A)=\rho \cdot k_{B}T$.

- Mit $\rho =(N/V)$ also $\;PV=N\cdot k_{B}T,$ wo N=Gesamtzahl der Atome.

Dies ist ein Zustandsmodell des idealen Gases. In alten Zeiten wurde noch mit Molen, Loschmidt- und Avogadro-Zahlen und irgendwelchen Gaskonstanten R hantiert, aber das lassen wir hier bleiben.

Das Modell der inneren Energie mit f=3 bis 6 Freiheitsgraden sei: $U=(f/2)\;N\cdot k_{B}T$.

Etwas allgemeiner wäre ein Energiemodell U=U(T,V), das an den Freiheitsgraden etwas ändert, wenn die Atome/Moleküle nahe zusammenrücken. Seine Wärmekapazität pro Molekül bei festem Volumen ist definiert als $c_{V}=(\partial U/\partial T);\;=(f/2)\;k_{B}$ hier.

Wenn zwei Objekte ungleicher Temperatur in Kontakt kommen, gibt es einen Energiefluss oder Wärmeaustausch, bis die Temperaturen sich ausgleichen. Denn nur solches ist ein stationärer, stabiler Zustand des Gesamtsystems. Spontane lokale Fluktuationen der Temperatur kommen zwar mikroskopisch und eher selten vor, aber makroskopisch nie! Der Zweite Hauptsatz verbietet sowas.

Der Erste Hauptsatz besagt nichts anderes als die Erhaltung der Energie. Wird zum Beispiel ein Gas mit dem Kolben in einer Luftpumpe komprimiert, geht mechanische Arbeit als Integral von Kraft mal Weg (gleich: Druck mal Änderung von Volumen) in das Gas über. Seine innere Energie und daher seine Temperatur steigen. Man lässt ihm keine Zeit, die Temperatur mit der Umgebung auszugleichen. Die Kompression heißt in diesem Fall ein adiabatischer Prozess.

Was kann man ausrechnen an der Luftpumpe?

Das eingebrachte Element Arbeit $\delta A=-P\;dV$ ist eine Differenzialform auf der PV-Ebene. Das 'δ' soll andeuten, dass die Form nicht exakt ist, nicht die totale Ableitung einer Zustandsfunktion. Aber es gibt die totale Ableitung von U. Mit bekannten Zustandsgleichungen U=U(T,V) und T=T(P,V) ist δA gleich der Energieänderung (Erster Hauptsatz):

$-P\;dV=\delta A=dU=(\partial U/\partial T){\big [}(\partial T/\partial P)\;dP+(\partial T/\partial V)\;dV{\big ]}+(\partial U/\partial V)\;dV.$

Daraus folgt nun eine Differenzialgleichung vom Typ (dP/dV)=f(P,V), die hoffentlich lösbar wird.

$T(P,V)=PV/(Nk_{B});\;(\partial U/\partial T)=(f/2)Nk_{B};\;(\partial U/\partial V)=0\;\Rightarrow$

$-P\;dV=(f/2)V\;dP+(f/2)P\;dV;\quad (f/2)V(dP/dV)=-P(1+(f/2))\;\Rightarrow$

$(f/2)\;d(\log P)=-(1+(f/2))\;d(\log V)\;\Rightarrow \;(P/P_{0})=(V_{0}/V)^{(2+f)/f}$

Das ist die adiabatische Kurve des Paares (P,V) ausgehend von einem bekannten Startpunkt von Druck und Volumen. Bei f=3 beträgt der Exponent (5/3). Er wird auch der Adiabatenkoeffizient (γ) genannt.

Der Schall entsteht, wenn lokal in der Luft eine schnelle Druckänderung auftritt. Eine schwingende Lautsprechermembran beschleunigt beim Ausfahren die auftreffenden und abprallenden Luftmoleküle in Achsrichtung und verlangsamt solche beim Einfahren. Dadurch gibt es im Wechsel Überdruck und Unterdruck gleich vor der Membran. Die ideale Adiabatenkurve beschreibt adäquat, wie der Druck und die Ausdehnung (das Volumen) einer lokalen Luftmenge sich zusammen ändern. Statt mit dem Volumen wird oft mit der Dichte des Gases argumentiert, invers proportional zum gefüllten Raum. Die Druck-Dichte-Gleichung ist also $(P/P_{0})=(\rho /\rho _{0})^{\gamma }$. DieGröße ρ ist entweder die Massendichte oder die Dichte der Molekülzahl, wenn das Gas nur eine Sorte hat.

Die Kompressibilität K misst, welche Änderung des Volumens beziehungsweise der Dichte (mit dem anderen Vorzeichen) eine kleine Druckänderung bewirkt,

$\Delta P=-K\cdot (\Delta V)/V=K\cdot (\Delta \rho )/\rho \;\Rightarrow \;K=\rho (\partial P/\partial \rho )=\gamma P\;$.

Erhaltene Größen und Kontinuitätsgleichung

Für ein fluides dynamisches Medium (Gas, Flüssigkeit) kann man mehrere Dichten einführen, nämlich zusätzlich zur Massendichte $\rho (t,{\vec {x}})$ noch eine Energiedichte $\epsilon (t,{\vec {x}})$ und eine Impulsdichte ${\vec {\pi }}(t,{\vec {x}})$. Ein System von Kontinuitätsgleichungen sorgt dafür, die physikalischen Erhaltungssätze für diese Observablen zu formulieren. Die Dichten sind in Ort und Zeit variabel. Es gibt ein Drift-Geschwindigkeits-Vektorfeld ${\vec {u}}(t,{\vec {x}})$, das ausdrückt, wie schnell sich die Materie am Ort-Zeit-Punkt $(t,{\vec {x}}))$ momentan bewegt. Natürlich ist dieses die mittlere Geschwindigkeit aller Gasmoleküle in einem kleinen Testvolumen, wesentlich langsamer als die thermische Molekülbewegung. Eine Dichte ist ebenso mesoskopisch definiert: über ein kleines Volumen ein gemittelter Wert, der eine quantitative Eigenschaft jedes Atoms ist. So wie etwa seine Masse, seine kinetische Energie, sein Impuls.

Wie variiert eine Dichte irgendeiner Art $\rho (t,{\vec {x}})$ zeitlich in einem winzigen Testwürfel $V=\Delta x\times \Delta y\times \Delta z$ um den Punkt $(t,{\vec {x}})=(t,x,y,z)$ herum?

Es sei vorausgesetzt, dass die Dichte aus zeitlich erhaltenen Werten aufsummiert ist. Die Massendichte ist von dieser Art, da sich die Massen der Moleküle nicht ändern. Die Dichte der kinetischen Energie und des mittleren Impulses ebenfalls, wenn nach jedem elastischen Stoß von Atompaaren die Summe dieser Daten gleich bleibt. Ein Erhaltungssatz sagt dann folgendes:

Die direkte Zunahme $V\cdot (\partial \rho /\partial t)\Delta t.$ der Größe ρ im gegebenen Würfel kommt nur durch die Migration der Quantität mit Driftgeschwindigkeit ${\vec {u}}=(u_{x},u_{y},u_{z})$ aus der oder in die Nachbarschaft zustande. Die saloppe Herleitung der Entwicklung im Testwürfel geht wie folgt, die strengere würde Integralrechnung und einen Stokesschen Satz bemühen. Von links kommt hereintransportiert, mit dem Flächenfaktor brauchbar kodiert:

$(V/\Delta x)\cdot u_{x}(x)\;\rho (x)\;\Delta t.$

Von rechts symmetrisch dazu,

$(V/\Delta x)\cdot (-u_{x}(x+\Delta x))\;\rho (x+\Delta x)\;\Delta t.$

Die Summe beider ist $(V\;\Delta t)$ mal die Ableitung nach $x,\;-(\partial /\partial x)(u_{x}\cdot \rho ).$ Aus den zwei anderen Himmelsrichtungen gibt es gleichartige Beiträge. Die Variation mit der Zeit ist die Summe der hereingedrifteten Mengen und ergibt, geteilt durch $(V\cdot \Delta t),$

$\partial \rho /\partial t=-\partial _{x}(u_{x}\;\rho )-\partial _{y}(u_{y}\;\rho )-\partial _{z}(u_{z}\;\rho )$

Es folgt die Kontinuitätsgleichung mit dem Nabla-Operator.

$\partial \rho /\partial t+\partial _{x}(u_{x}\;\rho )+\partial _{y}(u_{y}\;\rho )+\partial _{z}(u_{z}\;\rho )=\partial _{t}\rho +{\vec {\nabla }}\cdot ({\vec {u}}\;\rho )=0$

Erinnerung: Nabla angewandt auf ein Skalarfeld f produziert ein Vektorfeld ${\vec {\nabla }}f$ und heißt der Gradient von f. Nabla angewandt auf ein Vektorfeld ${\vec {u}}$ in der Form $({\vec {\nabla }}\cdot {\vec {u}})$ produziert ein Skalarfeld und heißt die Divergenz von ${\vec {u}}$.

${\vec {\nabla }}f=(\partial _{1}f,\partial _{2}f,\partial _{3}f);\quad ({\vec {\nabla }}\cdot {\vec {u}})=\partial _{1}u_{1}+\partial _{2}u_{2}+\partial _{3}u_{3}.$

Zu einer Dichte gehört also immer ein Driftstromvektor, $(\rho ,\;{\vec {u}}\rho )$, der angibt, wie schnell die Sache zerfließt, sich räumlich fortbewegt. Die Impulsdichte eines Gases ist natürlich nicht erhalten, wenn eine Schwerkraft von außen angreift. Im kräftefreien Fall aber sei angenommen, dass die gemittelte Impulsdichte einfach als Massendichte mal Driftgeschwindigkeit geschrieben werden kann, ${\vec {\pi }}=\rho \;{\vec {u}}$.

Die Kontinuitätsgleichung der Massenerhaltung für die Massendichte ρ:

$\partial _{t}\rho +{\vec {\nabla }}\cdot ({\vec {u}}\rho )=\partial _{t}\rho +({\vec {u}}\cdot {\vec {\nabla }})\rho +({\vec {\nabla }}\cdot {\vec {u}})\rho =0$

Das Etappenziel sind die Euler-Gleichungen für das Bewegungsverhalten des Gases. Das sind Differenzialgleichungen für folgendes System mit fünf variablen Feldern.

- Das spezifische Volumen $v(t,{\vec {x}}):=1/\rho (t,{\vec {x}});$

- Der lokale Druck des Gases $p(t,{\vec {x}});$

- Die lokale Durchschnitts- oder Driftgeschwindigkeit ${\vec {u}}(t,{\vec {x}}).$

Aus der Kontinuität von ρ folgen zunächst zwei ähnliche Gleichungen für Druck und Volumen bei idealen Gasen. Das Gas soll sich adiabatisch verhalten. Der lokale Druck ist dann eine Funktion der lokalen Dichte der Form $p=C\cdot \rho ^{\gamma }$, mit einer Konstanten C. Das spezifische Volumen wird mit C=1, γ=-1 ein Fall derselben Form. Gesucht sind Kontinuitätsgleichungen für p,v.

$\partial _{t}p+{\vec {\nabla }}({\vec {u}}p)=C\gamma \rho ^{(\gamma -1)}(\partial _{t}\rho +({\vec {u}}\cdot {\vec {\nabla }})\rho )+({\vec {\nabla }}\cdot {\vec {u}})p=(\gamma p/\rho )(-{\vec {\nabla }}\cdot {\vec {u}})\rho +({\vec {\nabla }}\cdot {\vec {u}})p\;\Rightarrow$

$\partial _{t}p+{\vec {\nabla }}({\vec {u}}p)-({\vec {\nabla }}\cdot {\vec {u}})p=\partial _{t}p+({\vec {u}}\cdot {\vec {\nabla }})p=-\gamma ({\vec {\nabla }}\cdot {\vec {u}})p$

Das Paar von Gleichungen für Druck und Volumen ergibt sich daher.

$\partial _{t}p+({\vec {u}}\cdot {\vec {\nabla }})p=-\gamma ({\vec {\nabla }}\cdot {\vec {u}})p$

$\partial _{t}v+({\vec {u}}\cdot {\vec {\nabla }})v=({\vec {\nabla }}\cdot {\vec {u}})v$

Der Differenzialoperator $D_{t}:=[\partial _{t}+({\vec {u}}\cdot {\vec {\nabla }})]$ auf den linken Seiten wird gern interpretiert als die Zeitableitung mit Drift. Er bewertet die zeitliche Änderung der Felder für einen Beobachter, der mit dem Strom mitschwimmt. Wir nennen ihn die konvektive Zeitableitung.

Cauchy-Impulsterm, Euler-Gleichungen

An dieser Stelle zieht Newtons Mechanik in die Dynamik der Fluide ein.

Die Impulsdichte $(\rho {\vec {u}})$ bekommt eine Erhaltungsgleichung, die nicht zu Null aufgeht. Sondern ein Kraftfeld ${\vec {f}}(t,{\vec {x}})=\{f_{i}(t,{\vec {x}}\}$ macht eine Zeitveränderung der Impulse, die keine bloße Bilanz von Migration ist. Das heißt, für jede Kraftkomponente i:

$f_{i}=\partial _{t}(\rho \;u_{i})+{\vec {\nabla }}\cdot ({\vec {u}}\rho u_{i})=(\partial _{t}\rho )u_{i}+\rho \partial _{t}u_{i}+({\vec {\nabla }}\cdot ({\vec {u}}\rho ))\cdot u_{i}+({\vec {u}}\rho )\cdot ({\vec {\nabla }}u_{i})=$

$f_{i}=[\partial _{t}\rho +{\vec {\nabla }}\cdot ({\vec {u}}\rho )]\cdot u_{i}+\rho \cdot (\partial _{t}u_{i}+{\vec {u}}\cdot {\vec {\nabla }}u_{i})=\rho \cdot D_{t}u_{i}.$

Benutzt wird die konvektive oder Drift-Zeitableitung der Geschwindigkeit und die Tatsache, dass die eckige Klammer verschwindet. Denn diese Klammer enthält nur die Kontinuitätsgleichung der Dichte.

Die Zeitableitung der Impulse ist die Summe aller Kräfte. Werden Impulse betrachtet für die Teilchen die im Strom mitschwimmen, zieht die konvektive Ableitung des Geschwindigkeitsfeldes ein. Ganz wie Newton lehrte, Masse mal Zeitveränderung der Geschwindigkeit gleich Kraft.

Eine äußere Kraft ist typisch eine Gravitation $(\rho {\vec {g}})$ mit einem Beschleunigungsvektor. Wichtiger sind die inneren Kräfte zwischen benachbarten Volumen des Gases. Cauchy fand eine allgemeine Beschreibung dafür, mit Hilfe seines berühmten und schon erwähnten Tensors $\{\sigma _{ij}\}$.

Variables Cauchy-Spannungsfeld am Testwürfel: die x-Komponenten der Kraft

Das Spannungs-Tensorelement $\sigma _{ij}(t,{\vec {x}})$ wirkt, zieht, drückt, auf einer Flächennormalen i in Richtung der Achse j. Es übt eine Kraft auf ein Testvolumen aus, die über alle Flächen der Hülle des Volumens zu summieren (integrieren) ist, jeweils proportional zu den Flächenstücken. Hier ist $i,j\in \{x,y,z\}\equiv \{1,2,3\}$. Die Kraft ist proportional zur Fläche, also hat der Spannungstensor die Dimension eines Drucks. Welche Kraft wirkt auf unseren bereits benutzten Testwürfel mit seinem Volumen $V=\Delta x\times \Delta y\times \Delta z$, in einer der drei Richtungen j? Von links auf Achse 1 greift in Richtung j die Kraft

$-(V/\Delta x)\sigma _{1j}(x)$ zu, auf unser Impulselement $(V\rho {\vec {u}}_{j})$.

Von rechts

$(V/\Delta x)\sigma _{1j}(x+\Delta x).$

Die Summe ergibt eine Ableitung im Grenzfall kleiner Δx, nämlich $(V\;\partial _{1}\sigma _{1j})$. Die anderen zwei Himmelsrichtungen wirken entsprechend mit Feldern $\sigma _{2j},\;\sigma _{3j}$, so dass die Kraft die Form hat:

$f_{j}=\partial _{1}\sigma _{1j}+\partial _{2}\sigma _{2j}+\partial _{3}\sigma _{3j}=\sum \nolimits _{i}\partial _{i}\sigma _{ij}.\quad \Rightarrow \;$

Impulsgleichung in Kurzform:

$\rho \;D_{t}{\vec {u}}=(\nabla \cdot \sigma )+\rho {\vec {g}}.$

Das Problem ist nun, an den Spannungstensor heranzukommen. Er enthält allgemein seitlich gerichtete Kräfte, die von der Viskosität des Mediums herrühren. Mit solchen fällt die Navier-Stokes-Gleichung an, die immer noch erfolgreich zur Simulation von Strömungen um Autos und fliegende Vehikel herangezogen wird.

Wir beschränken uns aber auf den Fall, dass {σ} völlig diagonal ist mit gleichen Werten. Der Tensor soll den gewöhnlichen skalaren Druck darstellen.

$\sigma _{xx}()=\sigma _{yy}()=\sigma _{zz}()=-p(t,{\vec {x}}),$ alle anderen Komponenten =0. Das Vorzeichen kommt daher, dass ein positives $\sigma _{xx}$ auf einer, immer nach außen gerichteten, Flächennormalen einen Zug bedeutet. Der Druck auf die Flächen soll nach innen wirken. Die Interpretation der internen Kraft ist hier, dass ein Druckgradient die Teilchen im Mittel beschleunigt. Denn die elastischen Stöße von einer Seite sind heftiger und/oder zahlreicher als die von der anderen; also gibt es einen Netto-übertrag von Impulsen.

Der Spezialfall eines Druckfeldes $p(t,{\vec {x}})$ führt zu

$\rho \;D_{t}{\vec {u}}=(-{\vec {\nabla }}p)+\rho {\vec {g}};\quad D_{t}:=\partial _{t}+({\vec {u}}\cdot {\vec {\nabla }})$.

Wenn das spezifische Volumen v=1/ρ eingesetzt wird, ergibt sich der Satz von Gleichungen fürs adiabatische Gas und sein Feldsystem $\{{\vec {u}},p,v\}$:

${\begin{array}{ll}D_{t}{\vec {u}}&=-v({\vec {\nabla }}p)+{\vec {g}},\\D_{t}p&=-\gamma ({\vec {\nabla }}\cdot {\vec {u}})p,\\D_{t}v&=({\vec {\nabla }}\cdot {\vec {u}})v.\end{array}}$

Leonhard Euler, Spitzenmathematiker und Universalgenie, schrieb die Gleichungen schon um 1755 auf. Nur die Zustandsgleichung mit der Adiabatenkonstante, die p,v verknüpft, wurde von Laplace nachgeliefert.

Das Gleichungssystem ist linear in jeder Zeitableitung der fünf Komponenten und bilinear in allen Paaren, Komponente-A mal Raumableitung-der-Komponente-B. Ein nichtlineares System mit zahlreichen verzwickten Lösungen. Solche liefern brauchbare Modelle für viele Phänomene der Bewegung von Flüssigkeiten und Gasen.

Ein Beispiel ist der Effekt von Bernoulli, ein Zusammenhang zwischen der Strömungsgeschwindigkeit und dem Druck von Gasmassen. Die Gleichung von Bernoulli ist dafür verantwortlich, dass der ganze Flugverkehr überhaupt vom Boden abhebt. Wird eine strömende Luftmasse gezwungen, schneller zu fließen, verliert sie an Druck. Wird sie gebremst, steigt der Druck. Das erste passiert an der Oberseite, das zweite an der Unterseite eines Flügels. Der Druckunterschied bewirkt den Auftrieb.

Eine Variante der Gleichungen mit $({\vec {\nabla }}\cdot {\vec {u}})=0$ modelliert inkompressible Flüssigkeiten. Dies ist äquivalent zu der Forderung, dass die konvektive Zeitableitung der Dichte verschwindet. Druck und Dichte koppeln nicht mehr über die Adiabatengleichung und die Gleichung der Zeitableitung des Drucks geht so nicht mehr durch.

Die Physik der Fluide steht ausgiebiger hier: Mechanik flüssiger und gasförmiger Körper.

Herleitung der Schallgeschwindigkeit

In einer Raumdimension hat das Geschwindigkeitsfeld nur eine Komponente und die Nabla-Operatoren werden einfache x-Ableitungen. Notiert mit Punkt und Strich:

${\begin{array}{ll}{\dot {v}}+uv'&=vu'\\{\dot {u}}+uu'&=-vp'+g\\{\dot {p}}+up'&=-\gamma pu'\end{array}}$

Es sind drei Funktionen von zwei Variablen, in Vektor- und Matrix-Schreibweise:

${\frac {\partial }{\partial t}}{\begin{pmatrix}v\\u\\p\end{pmatrix}}+{\begin{pmatrix}u&-v&0\\0&u&v\\0&\gamma p&u\end{pmatrix}}\cdot {\begin{pmatrix}v'\\u'\\p'\end{pmatrix}}={\begin{pmatrix}0\\g\\0\end{pmatrix}}$

Nun wird eine brutale Näherung für kleine Signale versucht. Die drei Felder seien kleine Variationen um einen Gleichgewichtspunkt herum. Die quadratische Matrix [A] der Gleichung wird durch diese Mittelwerte approximiert, nur in den Ableitungen leben die Fluktuationen. Das ist dann ein lineares System. Reflexartig heißt es: Diagonalisieren! Das geht über die Nullstellen des Polynoms $\det[\mathbf {A} -\lambda \mathbf {U} ]=0$, was ergibt (Determinantenkunde hier nötig):

$(u-\lambda )[(u-\lambda )^{2}-\gamma pv]=0$.

Drei Lösungen: $\lambda =\{u,\;u+u_{1},\;u-u_{1}\}\quad {\text{mit }}u_{1}={\sqrt {(\gamma pv)}}.$

Zu diesen Lambdas gibt es also Moden (Linearkombinationen der Felder) mit der simplen Differenzialgleichung ${\dot {y}}=-\lambda y'$, der Gravitationsterm sei ausgeschaltet. Es sind Wellen vom Typ $y=f(x-\lambda t)$ mit einer beliebigen Signalform f(). Ruht die Luft, ist in der Gleichgewichtslage die mittlere Geschwindigkeit Null (u=0) und es ist genau $\lambda =\pm u_{1}=\pm {\sqrt {(\gamma pv)}}$ die Geschwindigkeit der Wellen.

Unser Gasmodell sagt $pV=k_{B}T$, wo V das Volumen pro Molekül ist. Sei die Molekülmasse m. Die Massendichte ist dann ρ=(1/v)= (m/V). Daher $pv=pV/m=k_{B}T/m$. Schnell in die Quadratwurzel eingebaut:

$u_{1}={\sqrt {(\gamma k_{B}T/m)}}={\sqrt {(\gamma RT/M)}}\quad$ (Schallgeschwindigkeit c)

Zur Probe setzen wir Werte ein für die Luftmoleküle (M ist die Masse für ein Mol):

{ γ=1,4; R=8,3 J/(K mol); T=(273+27=300) K; M= 0,029 kg/mol }.

Als Geschwindigkeit kommt heraus: c = 347 m/s, ganz realistisch für trockene Luft.

Folgerung: die hier ausgebreitete Theorie ist nicht von der Hand zu weisen. Sie sagt die Schallgeschwindigkeit weitgehend richtig voraus, ausgehend von der Masse der Moleküle und den Prinzipien der klassischen Mechanik und Wärmelehre. Mit allen gekoppelten Feldvariablen sollten die Wellen bestehen aus gleichzeitigen Vibrationen der Dichte, des Drucks, der Temperatur, der "mesoskopischen" Geschwindigkeit, der Auslenkung um die Ruhelage. Der Schall wird schneller bei höherer Temperatur und in leichteren Gasen.

Dreidimensionale Schallwellen

Zur groben Linearisierung teilt man wieder alle Felder auf in langsam variierende (Großbuchtstaben) Mittelwerte und schnelle kleine Fluktuationen. Beispielsweise:

$\;p(t,{\vec {x}})\to P(t,{\vec {x}})+p(t,{\vec {x}})$.

Wir zwängen in die Formeln bei Ableitungen nur die Kleinbuchstaben, bei ursprünglichen Feldern die Großbuchstaben. Über gleiche Indizes in demselben Term werde automatisch summiert. Es erfand Einstein diese Konvention, als er sich mit vielen Indizes von Tensoren herumplagte. Damit sehen die Euler-Gleichungen so aus, angenommen ohne äußere Kräfte und in der dritten Spalte bei verschwindendem Driftfeld $\{U_{k}\}$:

${\begin{array}{ll}\partial _{t}u_{i}&=-U_{k}\partial _{k}u_{i}-V\partial _{i}p&=-V\partial _{i}p\\\partial _{t}p&=-U_{k}\partial _{k}p-\gamma P\partial _{k}u_{k}&=-\gamma P\partial _{k}u_{k}\\\partial _{t}v&=-U_{k}\partial _{k}v+V\partial _{k}u_{k}&=V\partial _{k}u_{k}\end{array}}$

Die zweite Zeitableitung des Feldes p wird interessant:

$\partial _{t}^{2}p=-\gamma P\;\partial _{j}(\partial _{t}u_{j})=\gamma PV\;\partial _{j}(\partial _{j}p)=c^{2}\;\partial _{j}\partial _{j}p$

Hier ist c² das Quadrat der bereits diskutierten Schallgeschwindigkeit. Die zweiten Ortsableitungen sind eine Summe

$\partial _{j}\partial _{j}p=(\partial _{x}^{2}+\partial _{y}^{2}+\partial _{z}^{2})p=(\partial _{1}^{2}+\partial _{2}^{2}+\partial _{3}^{2})p=\Delta p$

Dieses Delta ist die Definition des Laplace-Operators und die Gleichung

$\partial _{t}^{2}p={\ddot {p}}=c^{2}\;\Delta p\quad$ ist die dreidimensionale Wellengleichung.

Die Druckfluktuation, genannt der Schalldruck, kann beispielsweise ebene Wellen schlagen der Form

$p(t,{\vec {x}})=f(({\vec {n}}\cdot {\vec {x}})-ct)$,

wo die Amplitude f() jede beliebige glatte Kurve sein darf. Der Vektor ${\vec {n}}$ ist normiert, $|{\vec {n}}|={\sqrt {({\vec {n}}\cdot {\vec {n}})}}=1$. Also ein Einheitsvektor. Er zeigt in die Richtung der Wellenausbreitung.

Für das Schnelle-Feld ${\vec {u}}=\{u_{i}\}$, wir nennen diese oszillierende Drift so denn sie ist was anderes als die Schallgeschwindigkeit, der ebenen Welle folgt

${\dot {u}}_{i}=\partial _{t}u_{i}=-V\;\partial _{i}p=-V\cdot n_{i}\cdot f'({\vec {n}}{\vec {x}}-ct).$

Das heißt, das Beschleunigungs-Vektorfeld zeigt in die Richtung der Wellenbewegung. Das Schnellefeld und das Feld der lokalen Auslenkung folgen aus Zeitintegration der Beschleunigung mit unbekannten ortsabhängigen Konstanten. Da wir annehmen, dass diese Konstanten im langsam variablen Hintergrund $\{U_{k}\}$ stecken, dass der so gut wie Null ist weil kein Wind herrscht, wird gefolgert: Die Vektorfelder der Schnelle und der Auslenkung von Schallwellen sind longitudinal, parallel zur Ausbreitungsrichtung des Schalls. Im Gegensatz zu festen Körpern, in denen sich verschieden schneller transversaler und longitudinaler Schall fortbewegt, gibt es bei Gasen nur longitudinalen Schall. Ein anderer Bereich der Physik, der Elektromagnetismus, macht genau das Gegenteil. Dessen vektoriellen elektrischen und magnetischen Felder erlauben sich ausschließlich transversale Wellen.

Die sinusförmigen ebenen Wellen mit vorgegebener Frequenz f haben die Gleichung $p(t,{\vec {x}})=p_{M}\;\sin(({\vec {k}}\cdot {\vec {x}})-\omega t)$. Darin ist p_M die Spitzenamplitude der Druckvariation, ω=2πf die Kreisfrequenz. Der Wellenvektor ${\vec {k}}$ zeigt in Richtung der Fortbewegung, seine Länge ist die inverse Wellenlänge,

$|{\vec {k}}|=2\pi /\lambda ;\;c=f\cdot \lambda ;\;\omega =c\cdot |{\vec {k}}|.$

Mit der konstanten Geschwindigkeit c sind Wellenvektor und Frequenz streng proportional und es gibt keine Dispersion. Der Schall aller Frequenzen wandert gleich schnell. Die Wellenformen der komplexen musikalischen Signale bleiben unverändert.

Welche Werte haben bei Alltagslärm denn die Amplituden des Schalldrucks, der Schallschnelle, seiner Auslenkung und auch seiner Temperaturvariation? Die Hörschwelle wird beim effektiven Druck 2x10^-5 Pascal angesetzt. Effektivwerte sind etwas kleiner als Spitzenamplituden. Der Druck der Atmosphäre liegt bei 10⁵ Pascal. Ganze 10 Zehnerpotenzen trennen diese Druckwerte. Kein Wunder dass die lineare Näherung der Gleichungen da funktioniert. Die Amplituden werden logarithmisch in Dezibel (dB) verfolgt. Ein Faktor 10 bedeutet einen Zuwachs von 20 dB. Als Nullpunkt (0 dB(A)) der akustischen Dezibel dient die Hörschwelle. Ein Schall von 60 dB(A) = 2x10^-2 Pa ist laut, noch komfortabel. Einer von 120 dB(A) = 20 Pa ist hart an der Schmerzschwelle und zerstört bei Dauerbelastung das Gehör.

Aufgabe. Ein Piepton von f=1000 Hz mit 40 dB(A) = 2 mPa sei in der Luft. Man berechne die Schnelle, die Auslenkung, die Temperaturschwingung.

Mehr zum Schall gibt es bei Grundlegendes zur Akustik. Der Inhalt ergänzt die (arg theoretischen?) Ausführungen hier.

Kugelwellen. Wie können die Lösungen der Wellengleichung aussehen, wenn ein kleiner Rundum-Strahler wie eine Amsel auf dem Baumwipfel die Quelle von Lauten ist? Versuchen wir Funktionen, die nur vom Radius abhängen.

${\bar {p}}(t,{\vec {x}})=p(t,r)=f(r)g(r-ct);\quad {\ddot {p}}=c^{2}fg''\;=?\quad c^{2}\Delta {\bar {p}};\quad r={\sqrt {(x^{2}+y^{2}+z^{2})}}$

Die Funktion g() ist eine beliebige Amplitude, die sich radial ohne Verformung ausbreitet und f() ist ein Faktor, der die Amplitude abklingen läßt. Mit allen bekannten Produkt- und Kettenregeln rechnen wir den Delta-Operator aus für radiale Funktionen.

$\partial _{x}{\bar {p}}=p'(r)\;\partial _{x}r=(x/r)p';\quad \partial _{x}^{2}{\bar {p}}=(x/r)^{2}p''+(p'/r)-(x^{2}p'/r^{3})\;\Rightarrow$

$\Delta {\bar {p}}=p''+(2p'/r)=(fg''+2f'g'+f''g)+2(f'g+fg')/r$

Faktor von g = Null? $\;0=f''+2f'/r\;\Rightarrow \;f''=-2f'/r\;\Rightarrow \;f'=(a/r^{2})\;\Rightarrow \;f=(-a/r)$ als Kandidat.

Faktor von g' =Null? $\;0=2f'+2f/r\;\Rightarrow \;f'=-f/r\;\Rightarrow \;f=b/r$

Beides geht also zugleich und f=1/r ist richtig als Amplitudenfaktor, $\;\Delta {\bar {p}}=fg''$.

Auslaufende Kugelwelle $p(t,r)=g(r-ct)/r;\quad r={\sqrt {x^{2}+y^{2}+z^{2}}}$

Die Amplitude fällt ab wie der Kehrwert des Abstands von der Quelle. Das gilt genauso für die Amplitude der Schallschnelle. Das Quadrat dieser Schnelle zählt bei der kinetischen Energiedichte, die der Schall transportiert, wenn beispielsweise g() ein kurzer Impuls ist. Die Dichte der Leistung fällt mit dem Quadrat des Abstands. Daher bleibt das Integral der Energiedichte über die Kugelfläche eines jeden Abstands konstant. Die Wellengleichung enthält also die Erhaltung der Energie.

Die Energieerhaltung kommt hier davon, dass die Euler-Gleichungen verlustfrei sind. In diesem vereinfachten Modell gibt es keine Reibung an Grenzflächen, keine Viskosität des Mediums, keine Wärmeleitung zum Temperaturausgleich. Nur reversible Prozesse. All die Komplikationen wurden nachträglich in die aufgebohrte Version eingebaut, nämlich die Gleichungen von Navier und Stokes.

Anhang. Saite bei nichtlinearer Auslenkung.

Skript zur numerischen Integration der Schwingung.

# nlsaite.py : Nichtlineare Saiten-Differenzialgleichung
from math import exp, log, sqrt, sin, cos, pi

def div(a,b) : return int(a/b)

def vsum(x,y,f,g) : # Vektorsumme mif Koeffs f,g
  n=len(x); z=[0]*n
  for i in range(n) : z[i]= f*x[i]+g*y[i]
  return z

def vscale(x,f) : # Vektor-Skalierung mit Faktor
  n=len(x); z=[0]*n
  for i in range(n) : z[i]= f*x[i]
  return z

def rungekutta(df,pars,y,t,d) : # schritt d auf Vektor y, df=ableitungsfunktion
  p=1/2.0; q=1/3.0; r=1/6.0; d2=d*p # RK vierter Ordnung
  k1= vscale( df(y,t,pars), d);    z=vsum(y,k1, 1,p)
  k2= vscale( df(z,t+d2,pars), d); z=vsum(y,k2, 1,p)
  k3= vscale( df(z,t+d2,pars), d); z=vsum(y,k3, 1,1)
  k4= vscale( df(z,t+d,pars), d)
  a=vsum(k1,k4, r,r); b=vsum(k2,k3, q,q)
  c=vsum(a,b, 1,1); return vsum(y,c, 1,1)

def integrate(mode, ableit, k, t,dt, params, x) : # Vektor-x Diffgleichung
  # ableit= Ableitungsfunktion, dt Total-Schrittweite, k Zahl der Mikroschritte
  # mode=0 naive Methode, mode=1 Runge-Kutta
  n= len(x); h1= dt/(k+0.0) # n=Zahl der Variablen, h1=Mikro-Schrittweite
  for j in range(k) : 
    if mode==0 :
      dx= ableit(x, t+j*h1, params)
      for i in range(n) : x[i] += h1*dx[i] 
    else : x= rungekutta(ableit,params, x, t+j*h1, h1)
  return x

def ikurven(ableit, pars, n,m, ta,tb, x) : # Kurve x(t) integriert, ta...tb
  yy=[[]]*(n+1); y=x; yy[0]=x; dt= (tb-ta)/(n+0.0) # n steps, m microsteps
  for i in range(n) :
    y= integrate(1, ableit, m, ta+i*dt, dt, pars, y); yy[i+1]=y
  return yy # yy[t=0..n-1][len(x)] loest Dgl in 2D

def saitenmodell(w,r) : # endpunkte x=0,x=1 haben zeitableitung Null
  # model: w"*(1+r*[1-(1+w'^2)^p]), x=0...1, p=-3/2., r=(1/delta)  big
  n=len(w); v=[0.0]*n; d=1.0/(n-1); d2=d*d
  for i in range(1,n-1) :
    w1=(w[i+1]-w[i-1])/2.0/d; w2=(w[i+1]+w[i-1]-2*w[i])/d2
    a=1.0+w1*w1; b=a*sqrt(a); q=r*(1.0-1.0/b); v[i]=w2*(1.0+q)
  return v # geschwindigkeit der Auslenkung w[] 

def stehwelle(dbg=0,r=0.0,wmax=1e-3) : # nonlin Saitenmodell
  # ddot(w)=w" ; L=1; w=sin(pi x); w"=-pi^2 w; w(t)= cos(pi t) Referenzwelle
  def rderiv(x,t,par) : # x ist Paar von Vektoren Ort, geschwindigk.
    r,n= tuple(par); v= saitenmodell(x[:n],r); return x[n:]+v
  n=51 # n=Ortsschritte, nt=Zeitschritte
  x=[0.0]*(2*n); d=1.0/(n-1); z=pi; par=[r,n]; minis=[]
  nt=100; mt=10; t=2.5 # minima near n(nt)=10 30 50, maxima t=1,2...
  for i in range(n) : x[i]= wmax*sin(i*d*z) # Anfangsbedingung
  y= ikurven(rderiv, par, nt,mt, 0,t, x) # gesucht: normen von y[k][0:n]
  if dbg>0 : print('ly='+str(len(y))+' ly0='+str(len(y[0])))
  z=0; ny=len(y); ypeak=[0]*ny
  for k in range(ny) :
    for i in range(n) : z=max(z,abs(y[k][i]))
    if dbg==1 : print(str([k,z])) 
    ypeak[k]=z; z=0
  for k in range(1,ny-1) :
    a=ypeak[k-1]; b=ypeak[k]; c=ypeak[k+1]
    if (b<a)and(b<c) : minis+=[k] 
    if (dbg>0)and(b<a)and(b<c) : print('Min at k='+str([k,a,b,c])) 
  for h in range(len(minis)) :
    k=minis[h]; sqs=[0,0,0] # get arg(min) of squaresum parabola
    for i in range(3) :
      for j in range(n) : z=y[k+i-1][j]; sqs[i]+= z*z # squaresum
    u,v,w= tuple(sqs); kmin= (u-w)/(2.0*(u+w)-4.0*v) # [k-1,k,k+1] args of sqs
    if dbg>0 : print('kmin='+str(k+kmin)) 
    minis[h]=k+kmin
  return minis # Liste der Zeitpunkte bei Minimal-Auslenkung
     
def saite_grosssignal() : # Frequenzanhebung bei grossen Amplituden
  print("Frequenzanhebung. Gross-Signal-Modell Saitenlaenge=1")
  print("  Relative Auslenkung 0.001 bis 0.010, Nichtlinearitaet (1/rho)=1000")
  minis= stehwelle(wmax=1); per=minis[1]-minis[0] # ref periode
  #for wmx in [1e-3,2e-3,5e-3,1e-2] :
  for k in range(1,11) :
    wmx= 1e-3*k # wmx= (Auslenkung / Laenge)
    minis= stehwelle(r=1e3,wmax=wmx); pw= minis[1]-minis[0] # periode
    print('Max Auslenkung='+str(wmx)+' Relative freq='+str(per/pw))

saite_grosssignal() # nichtlineares einfaches Saitenmodell

Hoch zum Seitenanfang