Post on 03-Feb-2016
description
1
Praktikum
„Numerische Strömungsmechanik“
C.-D. Munz, S. Roller
2
Überblick
I- Klassifizierung Differenzialgleichungen
II- Numerische Lösung der elliptischen Differenzialgleichungen 1-Problemanalyse
a-Problemdarstellungb-Differenzenquotientenc-Aufbau des LGSd-Lösungsverfahren
2-Jacobi-Verfahren3-Gauß-Seidel-Verfahren4-SOR-Verfahren5-LSOR-Verfahren6-Abbruchkriterien
II-Numerische Lösung der parabolischen Differenzialgleichungnen1-Problemanalyse2-Explizit 1.Ordnung3- Implizit 1.Ordnung
4- Explizit 2.Ordnung5-Implizit 2.Ordnung6-Splitting7-DFL Bedingung
III-Numerische Lösung von hyperbolischen Differenzialgleichungen1-Problemanalyse2-Diskretisierung3- Charakteristiken Theorie4- Upwind-Verfahren5-Vollimplizites-Verfahren6-Crank-Nicolson-Verfahren7-Lax-Wendroff-Verfahren8- Runge-Kutta-Verfahren9-MUSCL-Verfahren10-CFL Bedingung
3
Klassifizierung DGL
Lineare partielle Differentialgleichung 2. Ordnung in 2 Dimensionen
)u,uu,y,f(x,yuc
yxub
xua (1) yx2
22
2
2
elliptisch : 04acb
hparabolisc : 04acb
chhyperbolis : 04acb
2
2
2
4
Klassifizierung DGL
1. Elliptische Gleichung
yyxx2
2
2
2
uuyu
xu:Δu
Gleichung-Helmholtz : 0f 0, κGleichung-Poisson : 0f 0, κGleichung-Laplace : 0f κ
y)f(x,κuΔu (2)
hungDruckgleic leichung,Potentialg :nAnwendungeoblemeRandwertpr
5
Klassifizierung DGL
2. Parabolische Gleichung
ung Wärmeleit:AnwendungoblemRandwertprAnfangs
uu:u
t)y,u(x,u 0ν ,νΔu tu (3)
yyxx
6
Klassifizierung DGL
3. Hyperbolische Gleichungen
0Δuctu (4) 22
2
Anfangs- Randwertproblem
Anwendung: Wellengleichung
7
:leichungTransportglinearen zur n umschreibengRaumrichtueiner in ichung Wellenglediesich lässt uq und ucpMit tx
0qcp 0pcq
xt
xt
Als einfachste hyperbolische Gleichung mit 2 Raumrichtungen ergibt sich somit:
0 )5(
yub
xua
tu
Klassifizierung DGL
8
Numerische Lösung der elliptischen DGL
1.Problemanalyse
a- Problemdarstellung
1J0,...,j ,y jΔyy 1I0,...,i , x iΔx x:teGitterpunk
1JyyΔy ,
1IxxΔx :tenSchrittwei
JI, :teGitterpunkinneren der Anzahlorthogonal t,äquidistan :Gitter
y,y x,x :etRechengebi
sj
si
sese
eses
9
Numerische Lösung der elliptischen DGL / Analyse
b- Differenzenquotienten
Zentrales FD-Verfahren 2. Ordnung
i, j+1
i-1, j i, j i+1, j
i, j-1
5 Punktestern
y
x
Pi,j=(xi,yj)
ui,j≈u(xi,yj)
Ersetzen der Ableitungen in Poisson-Gleichung durch Differenzen ergibt:
J , ... 1,j , I , ... 1,i )y,f(xκuΔy
u2uuΔx
u2uu (6) ji2
1ji,ji,1ji,2
j1,iji,j1,i
10
c- Aufbau des Gleichungssystem
Numerische Lösung der elliptischen DGL / Analyse
Mit Sonderbehandlung des Randes ergibt sich:
444444
34343434
24242424
141414
43434343
3333333333
2323232323
13131313
42424242
3232323232
2222222222
12121212
414141
31313131
21212121
111111
44
34
24
14
43
33
23
13
42
32
22
12
41
31
21
11
cbadcba
dcbadca
ecbaedcba
edcbaedca
ecbaedcba
edcbaedca
ecbedcb
edcbedc
A
uuuuuuuuuuuuuuuu
u Schwach besetzte Matrix
J , ... 1,j , I , ... 1,i fueuducubua (6) ji,1ji,ji,j1,iji,ji,ji,j1,-iji,1ji,ji,
11
Numerische Lösung der elliptischen DGL / Analyse
d- LösungsverfahrenGleichungssystem: (7) Au=f
A I·J × I·J-Matrix mit Bandstruktur
u=(u11,u21,…,uI1,u12,…,uIJ)
-Gaußalgorithmus:
Ungünstig, rechnet mit allen Nullen zu hoher Speicheraufwand und Rechenzeit
-Thomasalgorithmus:
Nicht anwendbar wegen den Nullen zwischen den Diagonalen
12
Numerische Lösung der elliptischen DGL / Analyse / Verfahren
sinnvoller :
-Iterationsverfahren:
löst LGS bis zur vorgegebenen Genauigkeit
13
Numerische Lösung der elliptischen DGL / Analyse / Verfahren
Aufspaltung von A: A = -Ai+Ae
Aus (8) erhält man damit die
Iterationsvorschrift
)fu(AAu Δ(p)Δe
1i
1)(pΔ
ΔΔeΔi fuAuA (8)
Iterationsverfahren (Splittingverfahren)
P ist der iterationsindex
(7) Au=f
14
Numerische Lösung der elliptischen DGL / Analyse / Verfahren
Jacobi Verfahren: Ai = -diag(A) .
Gauß-Seidel Verfahren: Ai = -diag(A) – L
L untere Dreiecksmatrix ohne DiagonaleU obere Dreiecksmatrix ohne Diagonale
Iterationsverfahren (Splittingverfahren)
Programmtechnische Umsetzung
Matrizen A, Ai, Ae werden nicht berechnet, Ausgangspunkt ist Gleichung (6)
)fu(AAu Δ(p)Δe
1i
1)(pΔ
15
Numerische Lösung der elliptischen DGL
2.Jakobi-Verfahren
J , ... 1,j , I , ... 1,i f)u(uΔy
1)u(uΔx
1 du ji,(p)
1-ji,(p)
1ji,2(p)
j1,-i(p)
j1,i21)(p
ji,
J , ... 1,j , I , ... 1,i )y,f(xκuΔy
u2uuΔx
u2uu (6) jiji,2
1ji,ji,1ji,2
j1,iji,j1,i
κ2Δy
22Δx
2ji,1ji,1ji,2j1,ij1,i2ji,1dmit f)u(u
Δy1)u(u
Δx1u
d1
Nach uij aufgelöst:
Iterationsvorschrift:
16
Numerische Lösung der elliptischen DGL
3.Gauß-Seidel-Verfahren
J , ... 1,j , I , ... 1,i f)u(uΔy
1)u(uΔx
1 du ji,1)(p1-ji,
(p)1ji,2
1)(pj1,-i
(p)j1,i2
1)(pji,
J , ... 1,j , I , ... 1,i )y,f(xκuΔy
u2uuΔx
u2uu (6) jiji,2
1ji,ji,1ji,2
j1,iji,j1,i
Iterationsvorschrift:
Schon bekannt
17
Numerische Lösung der elliptischen DGL
4.SOR-Verfahren
J , ... 1,j , I , ... 1,i f)u(uΔy
1)u(uΔx
1du~ ji,1)(p1-ji,
(p)1ji,2
1)(pj1,-i
(p)j1,i2
1)(pji,
J , ... 1,j , I , ... 1,i )y,f(xκuΔy
u2uuΔx
u2uu (6) jiji,2
1ji,ji,1ji,2
j1,iji,j1,i
Iterationsvorschrift:
J , ... 1,j , I , ... 1,i u~ω u ω-1 u 1)(pji,
(p)ji,
1)(pji,
ΔyΔxβ ,
β1a ,
aω
2
21J
πcos2β1I
πcosa122opt
axation Überrel1laxation Unterre1
sparameterRelaxation
Gauß-Seidel
18
Numerische Lösung der elliptischen DGL
5.LSOR-Verfahren
J , ... 1,j , I , ... 1,i f)u(uΔy
1uΔx
1uκΔy2
Δx2u
Δx1
ji,(p)
1ji,1)(p1ji,2
1)(pj1,i2
1)(pji,22
1)(pj1,-i2
~~~
J , ... 1,j , I , ... 1,i )y,f(xκuΔy
u2uuΔx
u2uu (6) jiji,2
1ji,ji,1ji,2
j1,iji,j1,i
Iterationsvorschrift:
J , ... 1,j , I , ... 1,i u~ω u ω-1 u 1)(pji,
(p)ji,
1)(pji,
In x-Richtung wird ein tridiagonales Gleichungssystem gelöst
ΔyΔxβ ,
β1a ,
aω
2
21J
πcos2β1I
πcosa122opt
axation Überrel1laxation Unterre1
sparameterRelaxation
19
Numerische Lösung der elliptischen DGL
6.Abbruchkriterien
ε||uu||
und/oderε||fAu||
terienAbbruchkri
Δ(p)Δ
1)(pΔ
ΔΔ1)(p
Δ
Residuum fAu:R Δ1)(p
Δ
Mittel hesquadratisc |||| Maximum ||||
Δ
Δ
Die Verfahren können durch erfüllen der Abbruchkriterien beendet werden:
Genauigkeit des Ergebnis
Genauigkeit der Iteration
20
7.Verfahren der konjugierten Gradienten(CG)
Numerische Lösung der elliptischen DGL
Die Idee der Gradientenverfahren besteht darin,für das Gleichungssystem aus (7) ein Fehlerfunktional zu definieren,um dieses anschließend zu minimieren.
Das Fehlerfunktional:
F(u)=0.5(uTAu) – fTu
hat genau ein globales Minimum in u= u*
Dabei steht u* für die exakte Lösung des Problems aus (7), womit gilt: Au*=f
a- Grundidee
21
b- Mathematische Behandlung
n
n
2
1
R
u
uu
uu,TfAuTu21(u)F
gilt fAumit * Au)(u21Au)(u
21Au)(uuf T*T*T*T
** AuT)*(u21)uA(uT)*u(u
21
AuT)*(u21AuT)*(u
21AuTu
21(u)F
*AuT)*(u21*AuT)*(u
21
**T**T*
*T*
uu0)uA(u)u(u für minimal FcstAu)(u
0)uA(u)u(uwegen
(1)
(2)
(1)+(2)
Numerische Lösung der elliptischen DGL / CG Verfahren
22
Zur Bestimmung des Minimums von F setzen wir
)pαF(umin)pαF(u)F(udass bestimmen, so α wollen und
pαuuk
kk
Rα
kk
k1kk
kk
k1k
k
ungsvektor Suchrichtp,rung von u eine Näheu,Rα kkk
abhängig αnoch von nur )pαF(u minimum dasist p undu festem Bei kk
kkkk
Durch Differenzieren erhalten wir
f)2(Au)(u21αf).(Au)(pα.Ap)(p
21
)pα(uf)pαA(u)pα(u21)pαF(u
kTkk
kTk2k
kTk
kk
kTkk
kTkk
kkk
k
kTk
kTk
kkTk
kkTkk
kk
α Ap)(pf)(Au)(pα0f)(Au)(pα.Ap)(p)pα(uF
Numerische Lösung der elliptischen DGL / CG Verfahren
23
TTTTTT rf)(AufAu:F giltu fAuu21F(u)Mit
Man beginnt nun mit der Suche der Lösung mit einembeliebigen Startvektor und sucht in Richtung des steilsten Abstiegs:
Wir wählen nun als SuchrichtungDiese Wahl scheint geeignet zu sein, da F(u) in negativer Gradientenrichtung am stärksten abfällt.
Numerische Lösung der elliptischen DGL / CG Verfahren
c- Suchrichtungsvektor
kTkk r))F(u(p
kTk
kTk
kTk
kTk
k Ar)(rr)(r
Ap)(pf)(Au)(pαistDamit
Steilster Abstieg
24
fuAr Residuum dem undrAr
rrα teSchrittweider mit
rαuu
kk
kTk
kk
k
kkk1k
Insgesamt ergibt sich so die folgende Iterationsvorschrift:
Diese einfache Suchrichtung konvergiert allerdings nur relativschlecht. Eine deutliche Verbesserung kann durch die
Verwendung von konjugierten Gradienten erzielt werden.
Numerische Lösung der elliptischen DGL / CG Verfahren
Steilster Abstieg
25
Für das Verfahren der konjugierten Gradienten müssen lediglichdie Suchrichtungen so angepasst werden, dass sie A -orthogonalaufeinander stehen. Diese neuen Suchrichtungen werden dann stattdem einfachen negativen Gradienten in obigen Verfahren verwendet.
Numerische Lösung der elliptischen DGL / CG Verfahren
CG Verfahren
.sindorthogonalrundrdamit undApundp dass
r)(rr)(rβ wahldiesichert Dabei
pβrp
tmodifizierRisiduen der Verwendungder unter wirdichtungKorrekturr Die
k1kk1k
kTk
1kT1k
k
kk
1k1k
26
Insgesamt ergibt sich so die folgende Iterationsvorschrift:
kTk
1kT1k
k
kk
kkk
k1k
kk
1k1k
kTk
kTk
k
kk
k1k
r)(rr)(rβ
pAαrf)pα(uArResidiumdemund
pβrpngSuchrichtuder
pA)(pr)(rαteSchrittweidermit
pαuu
Numerische Lösung der elliptischen DGL / CG Verfahren
27
Numerische Lösung der parabolischen DGL
2
2
2
2
yu
xuΔu
Rν , y)f(x,f , t)y,u(x,umit
i, j+1
i-1, j i, j i+1, j
i, j-1y
x
1J0,...,j y,jyy 1I0,...,i x,ix xteGitterpunk
)0(t NtΔt :tweiteZeitschrit
1JyyΔy ,
1IxxΔx :tenSchrittwei
N :hritteder Zeitsc AnzahlJI, :teGitterpunkinneren der Anzahl
]t,[t]y,[y]x,[x :etRechengebiorthogonal t,äquidistan:Gitter
sj
si
12
sese
21eses
Approximation im Raum: zentrale Differenzen
1.Problemanalyse
0, fuut
instationäres Problem
28
Numerische Lösung der parabolischen DGL
2. Explizites Verfahren 1. Ordnung in der Zeit
ji
nji
nji
nji
nji
nji
nji
nji
nji f
yuuu
xuuu
tuu
,21,,1,
2,1,,1,
1, 22
Vorwärts Zentraler Differenzenquotient
tn+1
t n xi-1 xi xi+1
Differenzenstern
Kein LGS
29
Programmtechnische Umsetzung
Die Schrittweite wird über die DFL-Bedingung festgelegt, in dem ein Sicherheitsfaktor eingeführt wird:
)Δy,x(min ν
DFLΔt
41DFL :ngtenbedinguSchrittwei
1-N , ... 0,n , J , ... 1,j , I , ... 1,i
Δtf)u2u(uΔyνΔt)u2u(u
ΔxνΔtuu
22max
max
ji,n
1-ji,n
ji,n
1ji,2n
j1,-in
ji,n
j1,i2n
ji,1n
ji,
Auflösen nach 1
,njiu
Numerische Lösung der parabolischen DGL / Explizit
30
Numerische Lösung der parabolischen DGL
3.Implizites Verfahren 1. Ordnung (Euler-Verfahren)
ji
nji
nji
nji
nji
nji
nji
nji
nji f
yuuu
xuuu
tuu
,2
11,
1,
11,
2
1,1
1,
1,1,
1, 22
Rückwärts Zentraler Differenzenquotient
tn+1
t n
xi-1 xi xi+1
Differenzenstern
31
fν1u
ΔΔν1f
t ν1κ mit Gleichung eelliptisch
fν1u
ΔΔν1u
t ν1u Δ
fuνΔt
u-u
tisierungZeitdiskre nur
n
n1n1n
1nn1n
~,
Numerische Lösung der parabolischen DGL / Implizit
lineares Gleichungssystem
32
Numerische Lösung der parabolischen DGL
4.Explizites Verfahren 2. Ordnung (Du Fort-Frankel)
1,...,J1,...,I,jfür i
i,jf2Δy
n1i,jun
i,ju2n1i,ju
ν2Δx
n1,jiun
i,ju2n1,jiu
νΔt
1ni,ju1n
i,ju
tn+1
tn-1
tn
auflösbarexplizit rungStabilisie implizite
)u(u21u 1n-
ji,1n
ji,n
ji,
t
xxi-1 xi xi-+1
33
Numerische Lösung der parabolischen DGL / Du Fort-Frankel
1-N , ... 0,n , J , ... 1,j , I , ... 1,i
)u(u21u mit
ft 2)uu2(uΔy
tν 2)uu2(uΔx
t2νuu
1-nji,
1nji,ji,
ji,n
1-ji,n
ji,n
1ji,2n
j1,-in
ji,n
j1,i21-nji,
1nji,
Anlaufschritt nötig
34
tn+1
tn
tn+1/2
t
x
5.Implizites Verfahren 2. Ordnung (Crank-Nicolson)
Numerische Lösung der parabolischen DGL / Crank-Nicolson
J1,...,jI,1,...,ifür
fΔy
uu2uν
Δxuu2u
νΔt
uuji,2
21n1ji,
21n
ji,2
1n1ji,
2
21nj1,i
21n
ji,2
1nj1,i
nji,
1nji,
ZeitderinnDifferenzeZentrale
)u(u21u 1n
ji,n
ji,2
1nji,
xi-1 xi xi+1
35
Numerische Lösung der parabolischen DGL
fν2uΔ~u
Δt ν2f
~ ,
Δt ν2mit κ Gleichung eelliptisch
fν2uΔ~u
Δt ν2u
Δt ν2uΔ~
der Zeitin n Differenze zentrale)u(u21:u
fuΔ~νΔt
u-u
ungskretisiernur Zeitdi
nn
nn1n1n
1nn2/1n
2/1nn1n
36
6.Splitting-Verfahren
(Dimensionensplitting, Zwischenschrittmethode, Fractional Step)
Zerlegung:
f u (11)
0 uν - u (10)0uν - u (9)
t
yyt
xxt
fνu -νu - u fuν - u yyxxtt
Numerische Lösung der parabolischen DGL
37
Numerische Lösung der parabolischen DGL / Splitting
In jedem Zeitschritt wird (9), (10), (11) nacheinander 1. Ordnung gelöst
J1,...,jI,1,...,ifür fΔt
uu (11)
J1,...,jI,1,...,ifür 0Δy
u2uuν
Δtuu
(10)
J1,...,jI,1,...,ifür 0Δx
u2uuν
Δtuu
(9)
ji,
**ji,
1nji,
2
**1ji,
**ji,
**1ji,
*ji,
**ji,
2
*j1,i
*ji,
*j1,i
nji,
*ji,
a- Splitting-Methode implizit 1. Ordnung
38
Numerische Lösung der parabolischen DGL / Splitting
(9), (10), (11) werden jeder für sich 2. Ordnung genau gelöst (mit Crank-Nicolson-Verfahren):
Damit das Gesamtverfahren auch 2. Ordnung in der Zeit ist, muß die Reihenfolge der Schritte in jedem Zeitschritt vertauscht werden:
J1,...,jI,1,...,ifür fΔt
uu (11)
J1,...,jI,1,...,ifür 0Δy
u2uu2ν
Δyu2uu
2ν
Δtuu
(10)
J1,...,jI,1,...,ifür 0Δx
u2uu2ν
Δxu2uu
2ν
Δtuu
(9)
ji,
**ji,
1nji,
2
*1ji,
*ji,
*1ji,
2
**1ji,
**ji,
**1ji,
*ji,
**ji,
2
nj1,i
nji,
nj1,i
2
*j1,i
*ji,
*j1,i
nji,
*ji,
... ;Δt
(9) (10), (11), ;Δt
(11) (10), (9), ;Δt
(9) (10), (11), ;Δt
(11) (10), (9),
b- Splitting-Methode implizit 2. Ordnung
39
Numerische Lösung der parabolischen DGL
Die DFL Zahl steht für die dimensionslose Diffusionszahl, die in parabolischen Gleichungen auftritt:
7.Die DFL Bedingung
.DFLtx aussich ergibt ungnsausbreitInformatio numerische Die
.beschreibtDiffusion durch ungnsausbreitInformatio die die ,definieren x
alsx Raumgitterein über gkeit"geschwindiDiffusions" einehier sich lässt Es
xtDFL
max
2
40
2max
max
xDFLt
:ttweite Zeitschridiefür x festem beidamit und
xDFL
tx
:n wird vorgegebeDGl diedurch ist wie groß so mindestensdigkeit gsgeschwinAusbreitun numerische die dass
Bedingung, diesich ergibt ist stabilVerfahren dasDamit
Numerische Lösung der parabolischen DGL / DFL-Bedingung
41
Numerische Lösung der hyperbolischen DGL
f buauu yxt
u)y,f(x,f t),y,b(x,b t),y,a(x,a , t)y,u(x,u :mit
1J0,...,j y,jyy 1I0,...,i x,ix xteGitterpunk
)0(t NtΔt :tweiteZeitschrit
1JyyΔy ,
1IxxΔx :tenSchrittwei
N :hritteder Zeitsc AnzahlJI, :teGitterpunkinneren der Anzahl
]t,[t]y,[y]x,[x :etRechengebiorthogonal t,äquidistan:Gitter
sj
si
12
sese
21eses
1.Problemanalyse
42
Numerische Lösung der hyperbolischen DGL / Problemanalyse
Umformulierung als Erhaltungsgleichung
ubua u)y,f(x, u)y,(x,f
u)y,(x,f (bu)(au)u
yx
yxt
i+1, j
i, j+1
i-1, j i, j
i, j-1
y
x
i+1/2, ji-1/2, ji,j-1/2
i,j+1/2
x(au)(au) j,ij,i 2
12
1
Erhaltungseigenschaft: was aus einer Zelle ausströmt, strömt in die Nachbarzelle ein
Differenz dessen, was links ein und rechts ausströmtFluß über den linken bzw. rechten Rand
43
Numerische Lösung der hyperbolischen DGL / Problemanalyse
Splitting-Methode
u)y,(x,f u (14)
0 (bu)u (13)0 (au)u (12)
t
yt
xt
Verfahren in Erhaltungsform
J1,...,j I,1,...,i fΔtuu (14)
J1,...,j I,1,...,i )h(huu (13)
J1,...,j I,1,...,i )g(guu (12)
**ji,
**ji,
1nji,
*ji,
*ji,Δy
Δt*ji,
**ji,
nj,i
nj,iΔx
Δtnji,
*ji,
21
21
21
21
g, h numerische Flüsse
44
Numerische Lösung der hyperbolischen DGL / Problrmanalyse
0 (au)u(15) xt
Im Weiteren werden Verfahren angegeben, die Gleichungen der Form
lösen, d.h . Verfahren für eine Raumdimension.Treten weitere Dimensionen auf,
so werden sie gemäß des angegebenen Splitting-Vefahrens nacheinander gelöst.
45
Numerische Lösung der hyperbolischen DGL
2. Diskretisierung
)u (u b h)u (u a g
1ji,ji,ji,21
ji,
j1,iji,j,i21
j,i
21
21
21
21
Eingesetzt in das Verfahren ergibt sich für ai+½,,j=ai-½,,j:
ji,y
1ji,ji,1ji,ji,1ji,ji,ji,21
1ji,ji,ji,21
ji,hi,
ji,xj1,ij,ij1,ij,ij1,iji,j,i2
1j1,iji,j,i2
1j,ij,i
buΔy2
ububΔy
)u (u b-)u (u b
Δyhh
auΔx2
uauaΔx
)u (u a-)u (u a
Δxgg
21
21
21
21
21
21
21
21
21
21
21
21
Zentrale Differenz (2. Ordnung) im Raum
Im Raum
46
Numerische Lösung der hyperbolischen DGL / Diskretisierung
Die Ableitungen im Raum werden zu einem gemeinsamen Zeitpunkt gebildet. Für die DGl fehlt nun noch die Ableitung nach der Zeit zu diesem Zeitpunkt. Sie
ermöglicht dann das Fortschreiten in der Zeit. Entscheidend ist dabei der Zeitpunkt, zu dem die DGl angeschrieben wird.
Gebräuchlich für die Diskretisierung der Ableitungen in der Zeit sind:• Zentrale Differenz mit Mittelung (2. Ordnung)• Vorwärts- und Rückwärtsdifferenz (1. Ordnung)• Differenz mit Extrapolation (2. Ordnung)• Runge Kutta Verfahren höherer Ordnung
xb)max(a,
CFLt :sfaktor Sicherheit demmit alsosich ergibt es
genügen Bedingung-CFLder ttweite Zeitschridie mussVerfahren expliziten dieFür
max
In der Zeit
47
Numerische Lösung der hyperbolischen DGL
3.Charakteristiken Theorie0 (au)u xt Die Exakte Lösung von erhält man,Indem man die Kurven(C)
berechnet, auf denen u =const gilt (totale Ableitung=0).
t)a(x,dtdx:C
t
x
C
a konstant C ist eine Gerade
u konstant auf C u(x,t)=u(x-at,0)
0dtdxuu0dt
dxxu
tu
dtdu/Rt)(x,C xt
durch Identifikation ( 15 und 16 )
(16)
48
Numerische Lösung der hyperbolischen DGL
4.Upwind-Verfahren
analog h
0a falls u
0a falls )u (u
0a falls u
a g
21
21
21
21
21
21
ji,
j,ij1,i
j,ij1,iji,21
j,iji,
j,ij,i
Eingesetzt in das Verfahren ergibt sich:
ji,x
j,ij,iji,j,ij1,ij,i
j,ij,ij1,ij,iji,j,i
j,ij,i au 0a,a falls
Δxua-ua
0a,a falls Δx
ua-ua
Δx
gg
21
21
21
21
21
21
21
21
21
21
Linksseitige Differenz für a>0, rechtsseitige Differenz für a<0 (1. Ordnung)
VerfahrenInstabilesΔx2
u-uaΔt
uu nnnji,
1nji, j1,ij1,i
Idee: Upwind
49
0a falls Δx
u-uaΔt
uu
0a falls Δx
u-uaΔt
uu
nnn
ji,1n
ji,
nnnji,
1nji,
ji,j1,i
j1,iji,
Numerische Lösung der hyperbolischen DGL / Upwind
Differenzenbildung in die Richtung, aus der die Information kommt.
Information wird entlang der Charakteristik (PQ) transportiert.
Vorwärtsdifferenz (explizit 1. Ordnung) in der Zeit.Upwind (1.Ordnung) im Raum.
50
Numerische Lösung der hyperbolischen DGL
5.Vollimplizites Verfahren
0x2gg 1n
j1/2,i1n
j1/2,i
Δtuu n
ji,1n
ji,
Lineares Gleichungssystem
Zentrale Differenz (2. Ordnung) im Raum.Rückwärtsdifferenz (implizit 1.Ordnung) in der Zeit.
51
Numerische Lösung der hyperbolischen DGL
0gu 1/2n
ji,x1/2n
ji,t
2uu
un
ji,1n
ji,1/2nji,
6.Crank-Nicolson Verfahren
Durch Mittelung
Implizit 2.Ordnung in Raum und Zeit
Zentrale Differenzen um n+1/2
ergibt
Wie berechnet man die numerischen Flüsse im Zeitpunkt (n+1/2) ?
0Δxgg 1/2n
j1/2,i1/2n
j1/2,i
Δtuu n
ji,1n
ji,
52
Numerische Lösung der hyperbolischen DGL
7.Lax-Wendroff-Verfahren (x-Richtung)
)g(guu 21
21
21
21
nj,i
nj,iΔx
Δtnji,
1nji,
21
21
21
21
21
21
nj,i
nj,i
nj,i uag mit
Explizites Verfahren 2. Ordnung in Raum und Zeit
Δxuun
j,i2Δtn
j1,in
ji,21n
j,ix2Δtn
j,in
j,it2Δtn
j,in
j,i
nji,
nj1,i
21
212
12
121
21
21 a)u(u(au)uuuu
Taylorentwicklung
Prädiktor : berechne Hilfswert un+1/2
53
)k2k2k(kuu xg-g
k :Schritt 4.
Δtkuu xg-g
k :Schritt 3.
kuu xg-g
k :Schritt 2.
kuu xg-g
k :Schritt 1.
12346Δtn
ji,1n
ji,
***j,i
***j,i
4
3**ji,
***ji,
**j,i
**j,i
3
22Δt*
ji,**ji,
*j,i
*j,i
2
12Δtn
ji,*
ji,
nj,i
nj,i
1
21
21
21
21
21
21
21
21
Numerische Lösung der hyperbolischen DGL
8.Runge-Kutta-Verfahren 4. Ordnung (klassische Variante)
54
Dämpfungsterme
)1.0(0
)uu46uu4(uεD :Ordnung 4.
)1.0(0
)u2u(uεD :Ordnung 2.
nj2,-i
nj1,i
nji,
nj1,i
nj2,iee
nj1,i
nji,
nj1,iee
e
e
Numerische Lösung der hyperbolischen DGL / Runge-kutta
55
Numerische Lösung der hyperbolischen DGL
9.MUSCL Verfahren (x-Richtung)
a- Problemdarstellung
Flussformulierung: gi+1/2,j ist eine Approximation an das, was während des
gesamten Zeitintervalls t über den Rand i+1/2,j der Zelle i,j rein oder raus fließt.
Problem: man kennt nur uij, d.h. was zum Zeitpunkt tn insgesamt in der
Zelle i,j ist, aber nicht, wie es verteilt ist oder wie es sich innerhalb des Zeitschritts ändert.
Idee: innerhalb einer Zelle wird eine lineare Verteilung angenommen, so daß man den Fluß am Rand genauer bestimmen kann.
=> Explizites Verfahren 2. Ordnung in Raum und Zeit
56
Numerische Lösung der hyperbolischen DGL / MUSCL
b- Stückweise lineare Rekonstruktion
xi-1 xi xi+1 xi+2 x
uMonotonic Upwind Scheme for Conservation Laws
Statt anzunehmen, dass u konstant ist zwischen xi-1/2 und xi+1/2, nehmen wir
jetzt an, dass u in diesem Bereich linear verteilt ist, d.h. wir bestimmen eine Gerade und werten sie an den Rändern aus, um die Flüsse zu berechnen.
57
Numerische Lösung der hyperbolischen DGL / MUSCL
c- MUSCL Schema
Bestimmung der Geraden: wir benötigen einen Punkt und eine Steigung.Der Punkt ist xij mit dem Funktionswert uij.
Steigung: 2 Möglichkeiten, linksseitige oder rechtsseitige Differenz:
xuu
R,xuu
L j,ij,1iij
j,1ij,iij
Wir müssen eine der beiden oder eine Linearkombination davon auswählen. Dies geschieht mit einem sogenannten Limiter, der bestimmte Bedingungen
erfüllen muß.
58
Numerische Lösung der hyperbolischen DGL / MUSCL
Mathematische Theorie für skalare ErhaltungsgleichungenErweitert auf Systeme
TVD-Eigenschaft (Total Variation Diminishing)
Hinreichende Bedingung (A. Harten)
iall iall
0i
01i
ni
n1i uuuu
2uu
sx,uusx0
i1i
i
1ii
i
Limiter: TVD-Eigenschaft
59
1. Minmod-Funktion
2. Sweby‘s Steigungsberechnung (gewichteter Minmod)
j,ij,ii R,Lmodmins
sonst0
0ab,bafürb0ab,bafüra
b,amodmin
b,kamodmin,kb,amodminmax)a(signb,ask 2k1mit
Numerische Lösung der hyperbolischen DGL / MUSCL
Limiter: Beispiele
60
Numerische Lösung der hyperbolischen DGL / MUSCL
Rekonstruktion im Raum
Steigung sin
ni
ni
ni s
2Δxuu Randwerte zu tn
Rekonstruktion in der Zeit
1/2nn tt ni
ni1/2i
ni
ni
ni
ni
1/2ni uua
x2Δtuufuf
x2Δtuu
Um die 2. Ordnung auch in der Zeit zu bekommen, geht man prinzipiell genauso vor, man extrapoliert vom Zeitpunkt tnden Zeitpunkt tn+1:
Damit kann man von Zellmittelpunkt an den Rand extrapolieren
tni
21n
i u2Δtuu
In der Zeit kann man aber keine Steigung berechnen, da man nur die Werte zu einem Zeitpunkt hat. Man behilft sich, indem man die Zeitableitung durch Raumableitungen ersetzt:
61
Randbehandlung
Jetzt muß von der Zelle auf den Rand umgedacht werden. Der Fluß am Rand, gi+1/2,j ist jetzt:
xi-1 xi xi+1 xi+2 x
u
uij
ui+1,j
2/1nj,iu
2/1nj,1)(iu
u , ugg 1/2n1i
1/2ni
1/2n1/2i
Numerische Lösung der hyperbolischen DGL / MUSCL
62
Numerische Lösung der hyperbolischen DGL / MUSCL
Upwind-Verfahren mit MUSCL
analog h
a falls ua falls )u (ua falls u
a u , ugg g
21
21
21
21
2121
ji,
j,i1/2n
j,)(i
j,i1/2n
j,)(i1/2nj,i2
1j,i
1/2nj,i
j,i1/2n
j,1i1/2nj,i
1/2nj1/2,ij,i
0
0
0
1
1
63
Numerische Lösung der hyperbolischen DGL / MUSCL
d- MUSCL Prozedur (gesamt)
u , ugg 1/2n1i
1/2ni
1/2n1/2i
Steigung sinn
ini
ni s
2Δxuu Randwerte zu tn
1/2nn tt ni
ni
ni
1/2ni ufuf
x2Δtuu
FV-Schema: 1/2n1/2i
1/2n1/2i
ni
1ni gg
ΔxΔtuu
mit
64
10.CFL Bedingung
Numerische Lösung der hyperbolischen DGL
Die Neumannsche Stabilitätsanalyse zeigt ,dass die expliziten Verfahren bedingt stabil sind.
1xta
CFL Bedingung
CFL steht hier für Courant-Friedrichs-Lewy. Die CFL Zahl beschreibt die dimensionslose Konvektionsgeschwindigkeit, die in hyperbolischen Gleichungen auftritt.
mit a als Transportgeschwindigkeit der eindimensionalen linearen Transportgleichung. Die Geschwindigkeit, mit der das Verfahren Information verteilt ist
tx
65
Numerische Lösung der hyperbolischen DGL / CFL
Damit das gewählte Verfahren mit der vorgenommenen Diskretisierung stabil ist, muss die Informationsausbreitung des Verfahrens mindestens so groß sein, wie die der DGl, also bei einer Weitergabe von Information in einem Zeitschritt um ein Raumgitter:
xa
CFLt
:ttweite Zeitschridiefür somit sich ergibt x festesfür
1CFLCFL : bzw. , 1xta
:bzw. , atx
max
max
66
Numerische Lösung auf einem Gitter der Schrittweite x:
Fehler auf einem Gitter der Schrittweite x bzw. 2 x :
Konvergenzordnung q des Verfahrens:
qex ΔxCuu num
qΔxqx2
qΔx
2eΔx)2(Ce
ΔxCe
)2ln(
)eeln(
eelogq 2
ee Δx
x2
Δx
x2
2q
Δx
x2
Numerische Untersuchung der Verfahrensordnung