1 praktikum numerische strömungsmechanik c.-d. munz, s. roller

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,y

uc

yx

ub

x

ua (1) yx2

22

2

2

elliptisch : 04acb

hparabolisc : 04acb

chhyperbolis : 04acb

2

2

2

4

Klassifizierung DGL

1. Elliptische Gleichung

yyxx2

2

2

2

uuy

u

x

u:Δu

Gleichung-Helmholtz : 0f 0, κ

Gleichung-Poisson : 0f 0, κ

Gleichung-Laplace : 0f κ

y)f(x,κuΔu (2)

hungDruckgleic leichung,Potentialg :nAnwendunge

oblemeRandwertpr

5

Klassifizierung DGL

2. Parabolische Gleichung

ung Wärmeleit:Anwendung

oblemRandwertprAnfangs

uu:u

t)y,u(x,u 0ν ,νΔu t

u (3)

yyxx

6

Klassifizierung DGL

3. Hyperbolische Gleichungen

0Δuct

u (4) 2

2

2

Anfangs- Randwertproblem

Anwendung: Wellengleichung

7

:leichungTransportglinearen zur n umschreibe

ngRaumrichtueiner 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(

y

ub

x

ua

t

u

Klassifizierung DGL

8

Numerische Lösung der elliptischen DGL

1.Problemanalyse

a- Problemdarstellung

1J0,...,j ,y jΔyy

1I0,...,i , x iΔx x:teGitterpunk1J

yyΔy ,

1I

xxΔx :tenSchrittwei

JI, :teGitterpunkinneren der Anzahl

orthogonal 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


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

cba

dcba

dcba

dca

ecba

edcba

edcba

edca

ecba

edcba

edcba

edca

ecb

edcb

edcb

edc

A

u

u

u

u

u

u

u

u

u

u

u

u

u

u

u

u

u

Schwach besetzte Matrix

J , ... 1,j , I , ... 1,i fueuducubua (6) ji,1ji,ji,j1,iji,ji,ji,j1,-iji,1ji,ji,

11


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


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


Jacobi Verfahren: Ai = -diag(A) .

Gauß-Seidel Verfahren: Ai = -diag(A) – L

L untere Dreiecksmatrix ohne Diagonale

U 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


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,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

κ2Δy

22Δx

2ji,1ji,1ji,2j1,ij1,i2ji,

1dmit f)u(u

Δy

1)u(u

Δx

1u

d

1

Nach uij aufgelöst:

Iterationsvorschrift:

16


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


Schon bekannt

17


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


J , ... 1,j , I , ... 1,i u~ω u ω-1 u 1)(pji,

(p)ji,

1)(pji,

Δy

Δxβ ,

β1a ,

aω

2

21J

πcos2β

1I

πcosa122

opt

axation Überrel1

laxation Unterre1 sparameterRelaxation

Gauß-Seidel

18


5.LSOR-Verfahren

J , ... 1,j , I , ... 1,i f)u(uΔy

1u

Δx

1uκ

Δy

2

Δx

2u

Δx

1ji,

(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


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

πcosa122

opt

axation Überrel1

laxation Unterre1 sparameterRelaxation

19


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)


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

u

u

uu,TfAuTu2

1(u)F

gilt fAumit * Au)(u2

1Au)(u

2

1Au)(uuf T*T*T*T

** AuT)*(u2

1)uA(uT)*u(u

2

1

AuT)*(u2

1AuT)*(u

2

1AuTu

2

1(u)F

*AuT)*(u2

1*AuT)*(u2

1

**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)(u2

1αf).(Au)(pα.Ap)(p

2

1

)pα(uf)pαA(u)pα(u2

1)pαF(u

kTkk

kTk2k

kTk

kk

kTkk

kTkk

kkk

k

kTk

kTk

kkTk

kkTkk

kk

α Ap)(p

f)(Au)(pα0f)(Au)(pα.Ap)(p)pα(uF


23

TTTTTT rf)(AufAu:F giltu fAuu2

1F(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.


c- Suchrichtungsvektor

kTkk r))F(u(p

kTk

kTk

kTk

kTk

k Ar)(r

r)(r

Ap)(p

f)(Au)(pαistDamit

Steilster Abstieg

24

fuAr Residuum dem und

rAr

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.


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.


CG Verfahren

.sindorthogonalrundrdamit undApundp dass

r)(r

r)(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)(r

r)(rβ

pAαrf)pα(uArResidiumdemund

pβrpngSuchrichtuder

pA)(p

r)(rαteSchrittweidermit

pαuu


27

Numerische Lösung der parabolischen DGL

2

2

2

2

y

u

x

uΔ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 N

tΔt :tweiteZeitschrit

1J

yyΔy ,

1I


N :hritteder Zeitsc Anzahl


]t,[t]y,[y]x,[x :etRechengebi

orthogonal t,äquidistan:Gitter

sj

si

12

sese

21eses

Approximation im Raum: zentrale Differenzen

1.Problemanalyse

0, fuut

instationäres Problem

28


2. Explizites Verfahren 1. Ordnung in der Zeit

ji

nji

nji

nji

nji

nji

nji

nji

nji f

y

uuu

x

uuutuu

,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


3.Implizites Verfahren 1. Ordnung (Euler-Verfahren)

ji

nji

nji

nji

nji

nji

nji

nji

nji f

y

uuu

x

uuu

t

uu,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


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(u2

1u 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(u2

1u 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ν

Δx

uu2uν

Δt

uuji,2

21n

1ji,2

1n

ji,2

1n

1ji,

2

21n

j1,i2

1n

ji,2

1n

j1,in

ji,1n

ji,

ZeitderinnDifferenzeZentrale

)u(u2

1u 1n

ji,n

ji,2

1n

ji,

xi-1 xi xi+1

35


fν

2uΔ

~u

Δt ν

2f~

,Δt ν

2mit κ Gleichung eelliptisch

fν

2uΔ

~u

Δt ν

2u

Δt ν

2uΔ

~

der Zeitin n Differenze zentrale)u(u2

1: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


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ν

Δt

uu (10)

J1,...,jI,1,...,ifür 0Δx

u2uuν

Δt

uu (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

u2uu

2

ν

Δy

u2uu

2

ν

Δt

uu (10)

J1,...,jI,1,...,ifür 0Δx

u2uu

2

ν

Δx

u2uu

2

ν

Δt

uu (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


Die DFL Zahl steht für die dimensionslose Diffusionszahl, die in parabolischen

Gleichungen auftritt:

7.Die DFL Bedingung

.DFLt

x aussich ergibt ungnsausbreitInformatio numerische Die

.beschreibtDiffusion durch ungnsausbreitInformatio die die ,definieren x

alsx Raumgitterein über gkeit"geschwindiDiffusions" einehier sich lässt Es

x

tDFL

max

2

40

2max

max

xDFL

t

:ttweite Zeitschridiefür x festem beidamit und

xDFL

t

x

:n wird vorgegebeDGl die

durch 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 N

tΔt :tweiteZeitschrit

1J

yyΔy ,

1I


N :hritteder Zeitsc Anzahl


]t,[t]y,[y]x,[x :etRechengebi

orthogonal 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 21

21

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


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,x

j1,ij,ij1,ij,ij1,iji,j,i21

j1,iji,j,i21

j,ij,i

buΔy2

ubub

Δy

)u (u b-)u (u b

Δy

hh

auΔx2

uaua

Δx

)u (u a-)u (u a

Δx

gg

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


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)

0dt

dxuu0

dtdx

xu

tu

dtdu/Rt)(x,C xt

durch Identifikation ( 15 und 16 )

(16)

48


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,i

ji,j,ij1,ij,i

j,ij,i

j1,ij,iji,j,i

j,ij,i au

0a,a falls Δx

ua-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

Δtuu

nnnji,

1nji, j1,ij1,i

Idee: Upwind

49

0a falls Δx

u-ua

Δtuu

0a falls Δx

u-ua

Δtuu

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


5.Vollimplizites Verfahren

0x2

gg 1nj1/2,i

1nj1/2,i

Δt

uu nji,

1nji,

Lineares Gleichungssystem

Zentrale Differenz (2. Ordnung) im Raum.Rückwärtsdifferenz (implizit 1.Ordnung) in der Zeit.

51


0gu1/2n

ji,x

1/2n

ji,t

2

uuu

nji,

1nji,1/2n

ji,

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Δx

gg 1/2nj1/2,i

1/2nj1/2,i

Δt

uu nji,

1nji,

52


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

Δx

uunj,i2

Δtnj1,i

nji,2

1n

j,ix2Δtn

j,i

n

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 x

g-gk :Schritt 4.

Δtkuu x

g-gk :Schritt 3.

kuu x

g-gk :Schritt 2.

kuu x

g-gk :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


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


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

u

Monotonic 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


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:

x

uuR,

x

uuL j,ij,1i

ijj,1ij,i

ij

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


Mathematische Theorie für skalare Erhaltungsgleichungen

Erweitert auf Systeme

TVD-Eigenschaft (Total Variation Diminishing)

Hinreichende Bedingung (A. Harten)

iall iall

0i

01i

ni

n1i uuuu

2uu

sx,

uu

sx0

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ürb

0ab,bafüra

b,amodmin

b,kamodmin,kb,amodminmax)a(signb,ask 2k1mit


Limiter: Beispiele

60


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


62


Upwind-Verfahren mit MUSCL

analog h

a falls u

a falls )u (u

a falls u

a u , ugg g

21

21

21

21

21

21

ji,

j,i1/2n

j,)(i

j,i1/2n

j,)(i1/2nj,i2

1

j,i1/2nj,i

j,i1/2n

j,1i1/2nj,i

1/2nj1/2,ij,i

0

0

0

1

1

63


d- MUSCL Prozedur (gesamt)

u , ugg 1/2n1i

1/2ni

1/2n1/2i

Steigung sin

ni

ni

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


Die Neumannsche Stabilitätsanalyse zeigt ,dass die expliziten Verfahren bedingt stabil sind.

1x

ta 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 istt

x

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. , 1x

ta

:bzw. , at

x

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(

)ee

ln(

e

elogq 2

e

e Δx

x2

Δx

x2

2q

Δx

x2

Numerische Untersuchung der Verfahrensordnung

1 praktikum numerische strömungsmechanik c.-d. munz, s. roller

Documents