Numerische Simulationsmethoden in der...

Einfuhrung in MATLAB Klassifikation von DGL und RB 1D FD Frequenzbereichs-EM Belegaufgabe MT

Numerische Simulationsmethoden in derGeophysik

Ubung zur gleichnamigen Vorlesung, WS 2014/15

Institut fur Geophysik und Geoinformatik, TU Bergakademie Freiberg

3. Februar 2015


Numerische Simulation physikalischer Felder ist

• wichtig fur das Verstandnis physikalischer Vorgange,

• nutzlich fur die Abschatzung von Messeffekten,

• Grundlage fur die Interpretation von Messungen,

• wichtigster Baustein fur eine Inversion (Ubung imSommersemester) und

• ubergreifend anwendbar in vielen technischen Bereichen.


Inhalt

Inhalt

1. Einfuhrung in MATLAB

2. Klassifikation von Differentialgleichungen undRandbedingungen

3. 1D FD Diffusionsgleichung: Warmeleitung

4. 1D FD Magnetotellurik


MATLAB... Matrix Laboratory

• Metasprache zum wissenschaftlichen (numerischen) Rechnen

• System von freien und kommerziellen Toolboxen

• Funktionen zur Visualisierung von Modellen und Daten

• Erstellung von GUIs

• Compiler zur Erstellung von lauffahigen Programmen


Vektorisierung

Schleifen: Matrixform:

for i=1:10,

for j=2:8,

A(i,j)=B(i,j+1) A(1:10,2:8)=B(1:10,3:9)

end

end

Modularisierung

abgeschlossene Funktionen und Toolboxen


Grundstrukturen

Hilfe

>> help % Liste aller Hilfethemen

>> help topic % Hilfe zu Themen, z.B. help strfun

>> help function % Hilfe zu Funktionen

Kommandozeile

>> a=4*3 % Ausgabe des Ergebnisses

>> b=sqrt(a); % keine Ausgabe


Arbeitsspeicher

>> whos a b % Frage nach Variablen im Arbeitsspeicher

>> clear a % Loschen von Variablen im Arbeitsspeicher

>> whos a b


Vektoren

>> a=[1 2 3] % Zeilenvektor

>> b=[0;2;1] % Spaltenvektor

>> b’ % Transponieren

>> a-b’ % elementweise Subtraktion

>> c=a*b % Skalarprodukt

>> d=b*a % Vektormultiplikation


Matrizen

>> A=[1 2 3; 4 5 6; 0 2 0] % Matrix

>> x=A*b % Matrix-Vektor-

Multiplikation

>> x\A % direkter Gleichungsloser

nach Gauß-Elimination

>> x+b % Addition

>> x.*b % elementweise Multiplikation


Indizierung

Vektoren

>> a=1:10

>> b=0:2:20

>> b(3)

>> b(4:8)

>> b(a)

>> b(8:end)

>> b(6:2:end-1)


Matrizen

>> B=A+1

>> B(2,3)

>> B(2:3,1:2)

>> B(3,:)

>> B(:,2)

>> B(:)


Dateioperationen

>> pwd % aktuelles Verzeichnis

>> cd DIR % wechselt Verzeichnis

>> ls % listet Verzeichnis

>> save FILENAME % speichert Workspace

>> save -ASCII FILE % speichert im ASCII-Format

>> load FILENAME % ladt Workspace

>> load -ASCII a.dat % ladt Variablen aus a.dat


>> edit test.dat % Editieren von Dateien

>> fid = fopen(’test.dat’,’r’) % off. Datei zum Lesen

>> A = fscanf(fid,’%d’,[2 Inf]); % liest Daten ein

>> fclose(fid); % schließt Datei

>> fprintf(’Time=%.2f n=%d\n’,A); % Ausgabe


Zeitmessung

>> tic; % Beginn

>> pause(1.0); % 1s Pause

>> toc % Ende


Verzweigungen und Schleifensteuerung

if-Anweisung

>> if a == 5,

...

else

...

end;

for-Schleife

>> for i = 1:10,

...

end;


while-Schleife

>> while(k < kmax),

...

end;

switch-Anweisung

>> switch msg

case ’ok’

...

case ’cancel’

...

otherwise

...

end;


Nutzliche Befehle

>> size(A) % Große einer Matrix

>> length(a) % Lange eines Vektors

>> max(a), min(a) % Maximum/Minimum

>> find(a), find(a==1) % Finden von Elementen

>> A(find(A==1))=-1; % Ersetzen aller 1 durch -1

>> sort(a) % Sortieren

>> diff(a) % Differenz

>> [ a b c ] % nebeneinander

>> [ a;b;c ] % untereinander

>> zeros(m,n) % erzeugt m x n Matrix aus Nullen

>> ones(m,n) % erzeugt m x n Matrix aus Einsen


Graphik-Befehle

Darstellung von Kurven

>> plot(x,y); % Kurve

>> plot(x,y,’r+:’); % rot gestrichelt mit +

>> plot(x,y1,x,y2); % mehrere Kurven

Beschriftung: xlabel, ylabel, titlelogarithmische Darstellung: semilogx, semilogy, loglog

Flachenhafte Darstellung

>> imagesc(x,y,Z); % x,y..Koordinaten, Z..Matrix

>> contour(x,y,Z) % Isolinien

>> surf(z,z,A) % erhohte Oberflache farbcodiert

>> colorbar % Farbskala


0. Aufgabe: Potential einer Punktquelle

Berechnen Sie das elektrische Potential einer Punktquelle bei~r0 = (0, 2.5 m, 0) uber einem homogenen Halbraum vonρ = 100 Ωm im Gebiet x = −10...10 m, y = −10...10 m,z = 0...10 m und stellen Sie es dar!

>> x=-10:10;

>> y=-10:10; % Erzeugung von x,y,z

>> z=0:10;

>> [X,Y,Z]=ndgrid(x,y,z); % Vektorisierung

>> whos X Y Z


>> R=sqrt(X.^2+Y.^2+Z.^2); % Ortsvektor

>> RA=sqrt(...) % Abstand zwischen Quelle und

Messpunkt

>> PHI=I/(2*pi*sigma*RA); % Formel zur Berechnung

des Potentials einer

Punktquelle

>> imagesc(x,y,PHI(:,:,1)’);% Darstellung

>> title(’...’) % Titel der Graphik

>> xlabel(’...’) % Beschriftung der x-Achse

>> ylabel(’...’) % Beschriftung der y-Achse

>> colormap(bluewhitered)

>> hc=colorbar; % Farbskala

>> set(get(hc,’YLabel’),’string’,’\phi in V’)


Bemerkung

Fur die Darstellung der Potentialwerte in der x-z-Ebene bei y = y1und der y-z-Ebene bei x = x1 mussen aus PHI(:,y1,:) bzw.PHI(x1,:,:) die entsprechenden Matrizen PHIy1(:,:) undPHIx1(:,:) extrahiert werden. Verwenden Sie squeeze, z.B.

>> imagesc(x,z,squeeze(PHI(:,10,:))’);

ZusatzWandeln Sie das programmierte Skript in eine Funktion, die ausder Quellstromstarke I , der Leitfahigkeit σ und der Lage der Quelle(x0, y0, z0) das Potential ϕ an einem beliebigen Ort (x , y , z)berechnet, z.B.

PHI=PTSOURCE(I,SIGMA,X0,Y0,Z0,X,Y,Z)

und kommentieren Sie diese.


Typen von Differentialgleichungen

Man unterscheidet zur Beschreibung der Ausbreitungunterschiedlichster physikalischer Felder im Gebiet Ω drei Typenpartieller Differentialgleichungen (PDE) 2. Ordnung:

• elliptisch: Potentialverfahren (z.B. Geoelektrik),Elektromagnetik (EM) im Frequenzbereich

∇ · (a∇u) + bu = f , ∆u = f , (1)

• parabolisch: Diffusion/Warmeleitung, EM im Zeitbereich

∆u − a∂u

∂t= f , (2)

• hyperbolisch: Wellenausbreitung (z.B. Georadar, Seismik)

∆u − 1

c2

∂2u

∂t2= f . (3)


Randbedingungen

Drei Arten von Randbedingungen (RB) sind auf dem Rand δΩandwendbar:

• Dirichletsche RBu = r , (4)

• Neumannsche RB∂u

∂ne~n = r , (5)

• Gemischte RB

α∂u

∂ne~n + βu = r . (6)


Elliptische partielle Differentialgleichung

... am Beispiel der Ausbreitung ebener diffusiver zeitharmonischerelektromagnetischer Felder

Grundlage zur Beschreibung der Ausbreitung elektromagnetischerFelder sind die Maxwell-Gleichungen. Daraus kann man dieInduktionsgleichungen fur das elektrische Feld ~E und dasMagnetfeld ~H ableiten:

∆~F = µσ~F + µε~F = µσ~F +1

c2~F mit c2 =

1

µε, F = E ,H. (7)


Der Ansatz ebener zeitharmonischer Wellen der Form

~F = ~F0ei(ωt−k~e·~r) (8)

in großer Entfernung von der Quelle mit ~e = (0, 0, ez)T undk~e = (0, 0, kz)T liefert unter Berucksichtigung der quasistationarenNaherung (T ≥ 10−5 s):

∆~F = iωµσ~F = −k2~F , k2 = −iωµσ. (9)


Die Induktionsgleichung fur das elektrisches Feld E:

∂2Ex

∂z2+ k2Ex = 0 bzw. (10)

−∇ · (c∇u) + au = f mit c = 1, a = −k2 = iωµσ, f = 0 (11)


Zusatzlich beschreiben inhomogene Dirichletsche Randbedingungendas Abklingen der Felder an den außeren Modellrandern:

~F = ~F0e−ikz mit F = E ,H und k = −

√−iωµσ. (12)


The 1D MT (magnetotelluric) problem (in frequency domain):

E ′′(z)− iωµ0σ(z)E (z) = 0 (0 < z <∞),

E ′(0) = −iωµ0, E (∞) = 0(13)

(see Haber et al. (2000)). Here,

z ∈ [0,∞) depth,E : [0,∞)→ C electric field,σ(z) conductivity,ω frequency,µ0 vacuum permeability (µ0 = 4π · 10−7).


For general σ the boundary problem (13) has to be solvednumerically:

At E (zf ) = 0, wherezf = 5d

with skin depth d =√

2µ0σω

[m].

We now replace the MT problem (13) by

E ′′(z)− iωµ0σ(z)E (z) = 0 (0 < z < zf ),

E ′(0) = −iωµ0, E (zf ) = 0.(14)

We set E (z) = u1(z) + iu2(z), u1 = Real(E ), u2 = Imag(E ). Then

u′′1 (z) + ωµ0σ(z)u2(z) = 0, u′′2 (z)− ωµ0σ(z)u1(z) = 0,

u′1(0) = 0, u′2(0) = −ωµ0, u1(zf ) = u2(zf ) = 0.(15)


Finally, with u3(z) := u′1(z) and u4(z) := u′2(z), we obtain a firstorder system,

u′1(z)− u3(z) = 0,

u′2(z)− u4(z) = 0,

u′3(z) + ωµ0σ(z)u2(z) = 0,

u′4(z)− ωµ0σ(z)u1(z) = 0

(16)

with ‘initial’ conditions

u4(0) = −ωµ0, u1(zf ) = u2(zf ) = u3(0) = 0.


This can be written asu1(z)u2(z)u3(z)u4(z)

′

−

0 0 1 00 0 0 10 −ωµ0σ(z) 0 0

ωµ0σ(z) 0 0 0

u1(z)u2(z)u3(z)u4(z)

=

0000

with

u1(zf )u2(zf )u3(0)u4(0)

=

000−ωµ0

.We discretize (16) on a mesh

0 = z0 < z1 < · · · < zN = zf

such that the mesh spacing is uniform in t = ln(z) (more precisely,we choose N = 100 equidistant points0 = t1 < t2 < · · · < tN = ln(zf ) and set z0 = 0 , zj = exp(tj)(1 ≤ j ≤ N)).


We further assume that σ = σj+1/2 is constant in each subinterval[zj , zj+1] (j = 0, 1, . . . ,N − 1). We denote by u`,j ourapproximations to u`(zj) (` = 1, . . . , 4) and set hj = zj − zj−1.

A discrete version of (16) (resulting from the midpoint scheme)now reads as

u1,j − u1,j−1

hj−

u3,j + u3,j−1

2= 0,

u2,j − u2,j−1

hj−

u4,j + u4,j−1

2= 0,

u3,j − u3,j−1

hj+ ωµ0σj−1/2

u2,j + u2,j−1

2= 0,

u4,j − u4,j−1

hj− ωµ0σj−1/2

u1,j + u1,j−1

2= 0 (1 ≤ j ≤ N),

u4,0 = −ωµ0, u1,N = u2,N = u3,0 = 0.


If we write uuu(j) = [u1,j , u2,j , u3,j , u4,j ]T (j = 0, 1, . . . ,N) then

B0uuu(0) =

[0−ωµ0

], where B0 =

[0 0 1 00 0 0 1

],

BNuuu(N) =

[00

], where BN =

[1 0 0 00 1 0 0

],

and, for 1 ≤ j ≤ N,(− 1

hjI4 − Tj(ω)

)uuu(j−1) +

(1

hjI4 − Tj(ω)

)uuu(j) = 000,

where

Tj(ω) =1

2

0 0 1 00 0 0 10 −ωµ0σj−1/2 0 0

ωµ0σj−1/2 0 0 0

.


For a fixed frequency ω = ωk , this results in a linear system

Ak(σ)uuuk = bbbk , where

uuuk =

uuu(0)

uuu(1)

...

uuu(N)

, bbbk = −ωkµ0eee2 ∈ R4(N+1) and

Ak(σ) =

B0

S1 R1

S2 R2

. . .. . .

SN−1 RN−1

SN RN

BN

∈ R4(N+1)×4(N+1)

with Sj = − 1hjI4 − Tj(ωk), Rj = 1

hjI4 − Tj(ωk).


Aufgabenstellung Teil 1

1. Realisieren Sie die Teilmatrizen T , S , R, B0 und BN inMatlab.

2. Erstellen Sie aus diesen funf Matrizen die Systemmatrix A furden Funf-Schicht-Fall ( 1

σ = ρ = [250, 25, 200, 10, 25]Ωm undh = [600, 1400, 3800, 4000, 5000]m) mit einer Frequenzf = 10Hz . Zur Kontrolle laden Sie sich die Matrix AA1

herunter:www.geophysik.tu-freiberg.de/assets/media/Personal/wwilhelms/Modellierung/AA-1.mat. Bilden Sie die DifferenzIhrer Matrix A und AA1. Nutzen Sie ’spy(A-AA1)’ oder’norm(A-AA1)’ um sich von der Richtigkeit IhrerProgrammierung zu uberzeugen.

3. Erweitern Sie die Systemmatrix A, so dass die Simulation furdrei Frequenzen (0.1Hz , 1Hz , 10Hz) moglich wird.


For three different frequencies ω1, ω2, ω3, this results in a(decoupled) linear system Auuu = bbb, where

A(σ) = diag(A1(σ),A2(σ),A3(σ)) ∈ R4·3(N+1)×4·3(N+1)

and bbbT =[bbbT1 ,bbb

T2 ,bbb

T3

]∈ R4·3(N+1).



1. Schreiben Sie eine Funktion MT1dVorwaerts.m, die alsEingabeparameter Schichtmachtigkeiten h, Leitfahigkeiten σund Frequenzen f braucht. Ausgabeparameter sollen dieSystemmatrix A und die rechte Seite bbb sein.

2. Fangen Sie Fehler ab, wenn zu wenige Parameter beimFunktionsaufruf eingegeben werden.

3. Bisher gehen wir davon aus, dass die Leitfahigkeit σ in denKnoten dieselbe ist, wie die in der darunter liegenden Schicht.Das ist jedoch nur beim homogenen Halbraum richtig. Furden allgemeinen Fall muss die Leitfahigkeit σ im Knotengewichtet werden:σ(zi ) = σi+1·hi+1+σi ·hi

hi+1+hifur alle in der Erde befindlichen Knoten

und an der Erdoberflache σ(z0) = σ12 .

4. Die sich dann ergebende Matrix A soll mit AA2 verglichenwerden.



1. Um die Ergebnisse der Vorwartsrechnung richtig bewerten zukonnen, sollen aus den komplexen elektrischen FeldernE = Er + iEi scheinbare spezifische elektrische Widerstande ρaund Phasen φ berechnet werden:ρa = 1

ωµ0|Z |2, φ = tan−1( Zi

Zr) mit Z = E/H

An der Erdoberflache wahlen wir H = 1 Am o.B.d.A., so dass

dann gilt Z = E .

2. Stellen Sie die Ergebnisse in Form einer Sondierungskurvegrafisch dar.



Als letzte Aufgabenstellung soll das Modell fur eine hohereGenauigkeit der Simulationsergebnisse feiner diskretisiert werden.Dazu ist eine Intervallaufteilung dz vorgegeben, die anstelle von him Aufruf der Funktion MT1dVorwaerts.m verwendet werden soll.

1. Laden Sie sich die Intervallaufteilung herunter: http://tu-freiberg.de/geophysik/research/electromagnetics/members/wenke-wilhelms/teaching/simulation/dz-uebung neu.mat.

2. Berechnen Sie mit der Funktion cumsum die tatsachlichenTiefen der jeweiligen Intervalle und vergewissern Sie sich, dassdie Tiefen der Leitfahigkeitsgrenzen des 5-Schicht-Fallsebenfalls vertreten sind.

3. Erstellen Sie nochmals Sondierungskurven.


Belegaufgabe MT

Simulieren Sie MT-Messungen fur elf Frequenzen im Intervall von0.01Hz bis 1Hz uber einem geschichteten Halbraum bestehend auszwei Schichten mit Leitfahigkeiten von 0.01 und 0.1 S

m , beide1000m machtig und einem abschließenden Halbraum mit einerLeitfahigkeit von 1 S

m .Verwenden Sie dazu die Formulierung des Problems wie imGleichungssystem (14) als Grundlage. Im Gegensatz zur Ubung sollnun ein System zweiter Ordnung gelost werden. Das bedeutet, wirnutzen nicht nur erste sondern auch zweite raumliche Ableitungenund bekommen so ein ganz anders strukturiertes Gleichungssystem.Wir werden das elektrische Feld E nicht in Real- und Imaginarteilaufspalten und bekommen so nur eine Gleichung pro Knoten. Einedetaillierte Herleitung der Systemmatrix ist fur die letzte Ubunggeplant.


Belegaufgabe MT Teil 1

Im Theorieteil soll die Herleitung des Problems dargelegt werden.

1. Beginnend vom Gleichungssystem (14) soll dasZustandekommen der Systemmatrix A nachvollzogen werden.

2. Wir splitten die Systemmatrix A mit komplexwertigenEintragen auf in eine Matrix K , die rein reele Eintrage hat,und eine Matrix S , die den Imaginarteil des Problemsubernimmt:A = K − iωµ0S mit S = diag(σ1, σ2, ..., σN).Am Ende soll es moglich sein, anhand der dargelegtenHerleitung das MT-Problem in Matlab-Code zu uberfuhren.



Eine Diskretisierung muss vorgenommen werden. Laden Sie sichdazu folgende Intervallaufteilung dzdzdz herunter: http://tu-freiberg.de/geophysik/research/electromagnetics/members/wenke-wilhelms/teaching/simulation/dz-beleg neu.mat. Stellen Sie sicher,dass Schichtgrenzen des Modells mit Intervallgrenzenubereinstimmen (Matlab-Funktion ’cumsum’ sollte verwendetwerden).

1. Beginnen Sie mit dem homogenen Halbraum, das heißt,σ = 0.01 S

m ist konstant fur alle drei Schichten.

2. Ahnlich wie in der Ubung beginnen Sie mit der Erstellung vonA bzw. K und S fur eine Frequenz.



1. Erweitern Sie im Anschluss das System auf elf Frequenzen mitfunf logarithmisch aquidistant verteilten Frequenzen proDekade (nutzen sie dazu die Matlab-Funktion ’logspace’).

2. Zur Erweiterung der Systemmatrix nutzen Sie moglichst dieMatlab-Funktion ’kron’. Die Matlab-Hilfe erklart dieVerwendung.

3. Erstellen Sie auch die rechte Seite bbb des Gleichungssystems.

4. Losen Sie das Gleichungssystem in Matlab mit Backslash.

5. Fuhren Sie eine Plausibilitatsprufung ihres bisherigen Codesdurch, indem Sie die Ergebnisse in Form einerSondierungskurve fur den scheinbaren spezifischen elektrischenWiderstand ρa fur den homogenen Halbraum darstellen.

6. Diskutieren Sie das Verhalten.



1. Losen Sie nun das Problem fur den geschichteten Halbraumindem sie σ fur dieselbe Intervallaufteilung dzdzdz entsprechendanpassen.

2. Stellen Sie die Ergebnisse der Formulierung erster Ordnung(Aufgabenstellung aus der Ubung) und die zweiter Ordnungfur den geschichteten Halbraum gegenuber. Erstellen Sie dazudrei Sondierungskurven: Real- und Imaginarteil deselektrischen Feldes E an der Erdoberflache, scheinbarenspezifischen elektrischen Widerstand ρa und die Phase φ.

3. Diskutieren Sie das Verhalten und stellen Sie einenZusammenhang her zu dem vorgegebenen Modell.

4. In einer vierten Abbildung stellen Sie den Feldverlauf (Real-und Imaginarteil) mit der Tiefe dar.

5. Diskutieren Sie diesen. Achten Sie dabei besonders auf dasVerhalten an den Schichtgrenzen.

Numerische Simulationsmethoden in der...

Documents

Transcript of Numerische Simulationsmethoden in der...