Numerische Methoden der Umwelphysik Computer- ¨Ubungen zur Wellengleichung 2. Februar 2006

Numerische Methoden der Umwelphysik

Computer- ¨ Ubungen zur Wellengleichung 2. Februar 2006

Die Ausbreitung von Erdbebenwellen (P-Wellen) wird in einfachsten F¨ allen durch die homogene Wellengleichung

∂

Φ

∂t

c

∂

Φ

∂x

= 0

beschrieben, wobei Φ die Auslenkung von der Ruhelage bezeichnet, und c die Ausbreitungs- geschwindigkeit der P-Wellen. Im homogenen Fall ist c ortsunabh¨ angig, im inhomogenen Fall ortsabh¨ angig.

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

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

3. Starten Sie das Programm wave.

Im Startfenster k¨ onnen Sie nun das Problem, die Modellgleichung, die Randbedingun- gen und die Anfangsanomalie ausw¨ ahlen und die entsprechenden Parameter eingeben.

Mit dem Start-Button starten Sie die Berechnung, deren Visualisierung in einem neu- en Fenster angezeigt wird.

2. Homogenes Medium: Verschiedene Arten von Randbedingungen

Wir betrachten die Wellengleichung auf einem eindimensionalen, endlichen Gitter. Damit stellt sich die Frage, wie die beiden R¨ ander behandelt werden sollen.

Φ

Φ

Φ

Φ

Φ

Φ

In dieser Aufgabe geht es darum, verschiedene Arten der Randbehandlung miteinander zu vergleichen. Bereits implementiert sind:

– Periodische Randbedingungen: Der linke und der rechte Randpunkt werden mit- einander identifiziert: Φ

= Φ

. Eine St¨ orung, die am rechten Rand rauswandert, wird dann am linken Rand wieder auftauchen, und umgekehrt.

– “Solid-wall” Randbedingungen: Die Randwerte werden fest vorgegeben, Φ

= Φ

= 0, ohne dass ein ¨ Ubergangsbereich wie bei der Relaxation existiert.

Aufgaben

a. Betrachten Sie die Advektionsgleichung und f¨ uhren Sie einen Test mit den obigen Randbedingungen aus.

b. Wiederholen sie nun die entsprechenden Tests mit der Wellengleichung

1

c. Die allgemeine L¨ osung der Advektionsgleichung ist Φ(x, t) = F(x

ct), wobei F eine (im wesentlichen) beliebige Funktion ist. Die allgemeine L¨ osung der Wellengleichung ist Φ(x, t) = F(x

ct) + G(x + ct), wobei F und G beliebig sind. Was stellen F und G f¨ ur das Advektions- und Wellenproblem dar? Was ist die Beziehung zwischen F und G bei der Wellengleichung mit “solid wall”-Randbedingungen? Warum ist beim Advektionsproblem keine Reflektion m¨ oglich?

d. Implementieren Sie eine Wellen-absorbierende Randbedingung:

Bei diesem Typ Randbedingungen werden die Randwerte so vorgegeben, dass es f¨ ur das numerische Schema scheint, als ob gar kein Rand vorhanden w¨ are. Die theoreti- schen Grundlagen und einige Hinweise zur numerischen Implementierung der Rand- bedingung am rechten Rand sind:

∗

Sie haben oben gesehen, dass die L¨ osung der Wellengleichung aus zwei Termen besteht:

Φ(x, t) = F (x

ct) + G(x + ct)

Der erste Term F stellt die nach rechts propagierende St¨ orung dar, der zweite Term G eine nach links propagierende St¨ orung. Die Reflexion ist auf die nach links laufende St¨ orung zur¨ uckzuf¨ uhren. Wir m¨ ussen also daf¨ ur sorgen, dass dieser Anteil in der L¨ osung nicht auftritt.

∗

Die neuen Werte Φ

f¨ ur 2

i

N

1 werden durch das unter Discretization gew¨ ahlte Verfahren bestimmt. Dabei entspricht die Diskretisierung “Centered time, centered space” dem Schema (7.18) im Skript. Der neue Wert Φ

am rechten Rand soll bestimmt werden aus

Φ

Φ

∆t + c

Φ

Φ

∆x = 0

Dies entspricht einer Diskretisierung der Advektionsgleichung mit Advektions- geschwindigkeit +c. Implementieren Sie dieses Verfahren (innerhalb der Funk- tion LocIntegrate) f¨ ur den rechten Rand. F¨ ur die numerische Implementation m¨ ussen Sie folgendes wissen:

Das Flag irelax bestimmt die Art der lateralen Randbedingung: irelax=3 Absorber am rechten (eastern) Rand, irelax=4 Absorber am linken (western) Rand, irelax=5 Absorber an beiden R¨ andern. Die lokale Wellengeschwindigkeit c steht im Array alpha(i), die zwei Zeitniveaus der Auslenkung sind in unow(i) und unew(i). Der Index i l¨ auft in allen F¨ allen ¨ uber 1:nx+4, wobei i=1,2 und i=nx+3,nx+4 intern reserviert sind. Die relevanten Indices sind also i=3 am linken und i=nx+2 am rechten Rand.

∗

Wie lautet der entsprechende Ansatz f¨ ur den linken Rand Φ

? Implementieren Sie auch das und bauen sie beides zusammen f¨ ur den Fall irelax == 5 ein!

∗

F¨ uhren Sie nun einen Advektions- und einen Wellenausbreitungstest aus, indem Sie den rechten Rand auf absobierend, den linken auf reflektierend setzen.

3. Wellenausbreitung in einem inhomogenen Wellenleiter

Eine zeitabh¨ angige St¨ orung am linken Rand f¨ uhrt zu einer Ausbreitung der St¨ orung nach rechts. Am rechten Rand wird eine Wellen-absorbierende Randbedingung verwendet.

2

Die Amplitude der St¨ orung und auch ihre Geschwindigkeit wird massgeblich durch die ortsabh¨ angige Ausbreitungsgeschwindigkeit c(x) (bzw. alpha(i) im Programm) bestimmt.

In dieser Aufgabe geht es darum, qualitativ die ¨ Anderungen zu verstehen.

a. Unter “Define the problem”, w¨ ahlen Sie “Homogenous/inhomogeneous wave con- ductor”. Vergleichen Sie die Ausbreitung einer sinusf¨ ormigen St¨ orung (Men¨ ueintrag activation, Wahl “Sine-like activation”) in einem homogenen Medium mit der Aus- breitung in einem inhomogenen Medium. W¨ ahlen Sie f¨ ur das inhomogene Medium einmal eine positive Anomalie und einmal eine negative Anomalie mit Amplitude 0.3.

b. W¨ ahlen Sie eine transiente St¨ orung, und erh¨ ohen Sie die Anzahl der Gitterpunkte auf 200 oder mehr. Die Amplitude der Anomalie setzen Sie auf -0.4.

4. Inhomogenes Medium: Eigenschwingungen

Es soll die Eigenl¨ osung Φ(x, t) = A(x)

exp(iωt) des Problems

∂

Φ

∂t

c(x)

∂

Φ

∂x

= 0 mit Φ(0, t) = Φ(L, t) = 0

gesucht werden. F¨ ur den Fall eines homogenen Ausbreitungsmediums ist das einfach ana- lytisch herleitbar. Der Ansatz Φ(x, t) = A(x)

exp(iωt) liefert eine Differentialgleichung f¨ ur A(x):

∂

A

∂x

+ ω

c

A = 0 mit A(0) = A(L) = 0

Mit dem Ansatz A(x) = sin(kx) und unter den Randbedingungen A(0) = A(L) = 0 ergeben sich die Eigenfrequenzen zu:

ω

= π

c

L

n mit n = 1, 2, 3, ...

Das l¨ asst sich auf eine Schwingungsdauer τ mittels 2π/τ = ω umrechnen: τ

= 2L/nc.

Aufgaben

a. Im homogenen Fall, c =const, sei c =

und L = 1. Bestimmen Sie ω

f¨ ur n = 1, .., 6 numerisch.

b. Der inhomogene Fall, c = c(x), ist im allgemeinen analytisch nicht mehr l¨ osbar. Um die Eigenfrequenzen ω und die Amplitudenverteilungen der dazugeh¨ origen Eigenl¨ o- sungen zu bestimmen, kann man zum Beispiel die “Shooting”-Methode verwenden.

Die einzelnen Schritte sind:

3

6

x = 0

x = L ω

ω

ω

“Shooting” f¨ ur drei verschiedene Frequenzen ω

, ω

, ω

. Bei der Wahl von ω

wird A(L) = 0, dh. diese Frequenz entspricht einer Eigenfrequenz.

∗

W¨ ahlen Sie eine Eigenfrequenz ω

, die nahe bei der erwarteten liegt. Bei kleinen Inhomogenit¨ aten wird die exakte Eigenfrequenz nur wenig von derjenigen des homogenen Falles abweichen. Sie k¨ onnen also eine solche als erste Sch¨ atzung verwenden.

∗

Jetzt wird die gew¨ ohnliche Differentialgleichung f¨ ur A(x) unter der linken Rand- bedingung A(0) = 0 bis zum rechten Rand x = L integriert.

∗

Ist auch dort A(L) = 0, so hat man eine Eigenl¨ osung und eine Eigenfrequenz gefunden. Geben Sie dazu den zu testenden Eigenwert im Feld Eigen frequency ein. Im Matlab Command Window erhalten Sie dann A(L).

∗

Andernfalls, dh. falls A(L) wesentlich von Null verschieden ist, w¨ ahlt man eine neue Sch¨ atzung f¨ ur die Eigenfrequenz, und f¨ uhrt die obigen Schritte nochmals durch.

c. Verwenden Sie die “Shooting”-Methode, um die Eigenfrequenzen des homogenen Fal- les (aus a.) zu testen.

d. Wiederholen Sie dann den Test f¨ ur ein inhomogenes Ausbreitungsmedium.

4