Modellering av åpent og isdekt hav med endelig differansemetoder for en ikke-lineær Schrödingerlikning

(1)

endelig differansemetoder for en ikke-lineær Schrödingerlikning

Vilde Hyggen Rosland

Masteroppgave, våren 2019

(2)

studiepoeng.

Forsiden viser et utsnitt av rotsystemet til den eksepsjonelle liegruppen E

8

, projisert

ned i planet. Liegrupper ble oppfunnet av den norske matematikeren Sophus Lie

(1842–1899) for å uttrykke symmetriene til differensiallikninger og spiller i dag en

sentral rolle i flere deler av matematikken.

(3)

Havis dekker store havområder både nord og sør på jordkloden. Med minkende isdekke i Arktis er det viktig å tilegne seg kunnskap om hvordan havet og isen påvirker hverandre. I denne oppgaven blir bølgepropagasjon i åpent hav og i havis studert for utvikling i tid og rom. Resultatene som blir presentert er utregnet ved simulering av ikke-lineær kubisk Schrödingerlikning med 4. ordens Runge-Kutta metode. Det blir lagt frem resultater for tids- og romutvikling av simulering av åpent og isdekt hav. I tilfellet isdekt hav, må det legges til et dempningsledd i Schrödingerlikningen for å få frem tendensene til utviklingen i is.

Abstract

Sea ice covers large areas of the sea in both the nothern and southern parts

of the globe. With decreasing ice cover in the Arctic, it is important to acquire

knowledge about how the sea and the ice affect each other. This thesis conserns

modeling and numerical simulations of wave propagation in open sea and in sea

ice, based on the non-linear cubic Schrödinger equation. A numerical method

is presented based on an explicit finite difference scheme with å fourth-order

Runge-Kutta time integration method. In terms of numerical experiments, results

are presented for time and space development of simulation of open and ice-covered

sea. In the case of ice-covered sea, a damping term must be added in the Schrödinger

equation to capture the features of wave propagation in ice.

(4)

(5)

Først vil jeg sende en takk til min veileder, Karsten Trulsen, for at du delte din kunnskap, gav gode råd og kom med viktige innspill. Jeg vil også takke for at du gjorde arbeide med masteroppgaven til en hyggelig opplevelse med din entusiasme, oppmuntring og generelt avslappede tone.

Flere har tilbudt hjelp underveis, og blant dem vil jeg særlig takke James David Trotter for all den tiden du satte av for å hjelpe meg, den verdifulle innsikten og forståelsen du kom med og inspirasjonen og motivasjonen du gav meg.

Jeg er takknemlig for mine medstudenter som alltid var villige til å besvare spørsmål og gi støtte. Jeg vil da spesielt takke Ingrid Olsen som har vært en fantastisk støttespiller hele veien.

Jeg er takknemlig for de som ellers har hjulpet meg med oppgaven ved å lese igjennom og gi korreksjoner og råd.

Tusen takk.

Vilde Hyggen Rosland

Mai 2019

(6)

Forord iii

Innhold iv

1 Introduksjon 1

2 Metode 5

2.1 Tidsutvikling med ikke-lineær kubisk Schrödingerlikning . . . . 5

2.1.1 Eksplisitt endelig differansemetode . . . . 8

2.1.2 Fysiske løsninger av Schrödingerlikningen . . . . 9

2.2 Romutvikling med ikke-lineær Schrödingerlikning . . . . 11

2.2.1 Analyse av stabiliteten til likningen . . . . 13

2.2.2 Analyse av resultatet fra stabilitetsanalysen . . . . 16

2.2.3 Endelig differansemetode . . . . 17

2.3 Parametere i Schrödingerlikningen . . . . 19

2.3.1 Dispersjonsrelasjonen . . . . 19

2.3.2 Gruppehastigheten og fasefarten . . . . 19

2.3.3 Parameteren β . . . . 20

2.3.4 Parameteren γ . . . . 20

3 Resultater 23 3.1 Rom- og tidsutvikling av overflatehevningen i åpent hav . . . . 23

3.1.1 Måling av karakteristiske bølgeparametere . . . . 26

3.1.2 Overflatehevning . . . . 28

3.1.3 Energispektrum . . . . 30

3.1.4 Energi og varians . . . . 33

3.1.5 Kurtose . . . . 34

3.1.6 Skjevheten . . . . 36

3.1.7 Bølgetall-frekvensspektrum . . . . 36

3.2 Romutvikling av overflatehevningen i isdekt hav . . . . 38

3.2.1 Måling av karakteristiske bølgeparametere . . . . 40

3.2.2 Overflatehevning . . . . 41

3.2.3 Energispektrum . . . . 42

(7)

3.2.4 Energi og varians . . . . 43

3.2.5 Kurtose . . . . 44

3.2.6 Skjevheten . . . . 44

3.2.7 Bølgetall-frekvensspektrum . . . . 45

4 Diskusjon 47 4.1 Sammenligning av resultater for åpent hav . . . . 47

4.2 Sammenligning av resultater for isdekt hav . . . . 48

5 Konklusjon 51

Bibliografi 53

A Python-skript: Stabilitetsanalyse 55

B Python-skript: Løsning av ikke-lineær kubisk Schrödingerlikning på

åpent hav (Runge-Kutta) 59

C Python-skript: Analyse av data fra løser av ikke-lineær kubisk Schrödingerlikning

på åpent hav (Runge-Kutta) 63

D Python-skript: Analyse av data fra eksisterende program [2] 83 E Python-skript: Løsning av ikke-lineær kubisk Schrödingerlikning for

isdekt hav (Runge-Kutta) 103

F Python-skript: Analyse av data fra løser av ikke-lineær kubisk Schrödingerlikning

for isdekt hav (Runge-Kutta) 107

(8)

(9)

KAPITTEL 1

Introduksjon

I denne oppgaven skal jeg undersøke hvordan amplituden til overflatehevningen utvikler seg i tid og rom, for åpent hav og isdekt hav, med en ikke-lineær kubisk Schrödingerlikning.

Det er blitt gjort utallige studier av havbølger, men det er fortsatt mye vi ikke vet. For et land som Norge, med lang kystlinje og med fisk og olje som to av hovedeksportvarene, er havbølger et viktig forskningsområde.

Noen vesentlige deler av havet er dekket av is, og det minkende isdekke i Arktis gir et økt behov for å forstå hvordan havet og havisen oppfører seg i samspill med hverandre.

I denne oppgaven skal jeg studere hvordan bølger forplanter seg i åpent og isdekt hav.

Hvordan bølger forplanter seg i is påvirker hvordan isen brytes opp, og med mer forskning kan man etterhvert føre varsling av havis.

For å beregne bølgepropagasjonen anvender jeg den ikke-lineære kubiske Schrödingerlikningen.

Schrödingerlikningen brukes i flere sammenhenger til å beskrive ikke-lineære bølger, som, for eksempel, propagasjon av laserstråler, havbølger og plasmabølger. Likningen gir en statistisk beskrivelse av envelopedynamikken til et bølgetog Ψe

^i(kx−ωt)

(der steilheten << 1), med en liten og sakte varierende amplitude, Ψ(x, t), bølgetallet k, vinkelfrekvensen ω . Den ikke-lineære kubiske Schrödingerlikningen er definert som

i ∂Ψ

∂t + β ∂

²

Ψ

∂x

²

+ γ|Ψ|

²

Ψ = 0,

der γ er knyttet til utvidelsen i potenser av bølgeintensitet og β er den dobbeltderiverte av vinkelfrekvensen, ω, med hensyn på bølgetallet, k [10].

Schrödingerlikningen er blitt brukt i mange studier av bølgeegenskapene til et åpent hav.

For eksempel, i artikkelen Evolution of a narrow-band spectrum of random surface gravity

waves av Dysthe, Trulsen, Krogstad og Socquet-Juglard [2], undersøker de egenskapene

(10)

til energispekteret til overflatehevningen av havbølger ved å utføre simuleringer basert på en ikke-lineær kubisk Schrödingerlikning og en modifisert Schrödingerlikning. I denne oppgaven ønsker jeg å se om jeg kan produsere liknende resultater, men ved å bruke en eksplisitt endelig differansemetode i stedet for spektral metode som brukes i artikkelen nevnt over. I artikkelen av Dysthe, Trulsen, Krogstad og Socquet-Juglard ser de hovedsakelig på tidsutviklingen av amplituden til overflatehevningen. Jeg vil også studere hvordan amplituden utvikler seg i rom, for å kunne sammenligne med et eksperiment der bølgen propagerer i x. Jeg må da undersøke om det er mulig å skrive om Schrödingerlikningen slik at jeg kan simulere for romutvikling.

Figur 1.1: Fotografi fra Svalbard

Det har blitt gjort få studier av Schrödingerlikningen for å modellere bølgepropagasjon på isdekt hav. I artikkelen Wave Propogation in a Solid Ice Pack av Liu og Mollo-Christensen [8] simulerer de bølgepropagasjon i is ved å utvide parametrene i Schrödingerlikningen slik at de inkluderer egenskapene til is. I Liu og Mollo-Christensen sin artikkel skriver de kun om parametrene i Schrödingerlikningen, selve likningen lar de være uforandret.

Jeg vil undersøke om det er nok å forandre parameterverdien eller om det må legges til et dempningsledd for å få frem tendensene som oppstår i is. Ved å sammenligne simuleringsresultatene med målinger fra kontrollerte forsøk, viser jeg at Schrödingerlikningen er en realistisk modell for bølgepropagasjon i is, dersom det legges til et dempningsledd.

Resten av oppgaven er strukturert som følger: I neste kapittel utledes den ikke-lineære Schrödingerlikningen for både tids- og romutvikling, og den endelige differansemetoden

som blir brukt til å løse likningen blir presentert. I tillegg beskrives likningen til overflatehevningen og parametrene i Schrödingerlikningen. I kapittel 3 blir resultatene fra utregningene av

den ikke-lineære Schrödingerlikningen for utvikling i både tid og rom, og med innsatte

verdier for både åpent hav og isdekt hav lagt frem. Resultatene sammenlignes med

(11)

andre studier. Kapittel 4 inneholder diskusjonen av resultatene. I kapittel 5 presenteres

konklusjonen til oppgaven.

(12)

(13)

KAPITTEL 2

Metode

I dette kapittelet utleder jeg de ikke-lineære kubiske Schrödingerlikningene som blir brukt til å regne ut tidsutviklingen og romutviklingen av amplituden til overflatehevningen.

Jeg utleder likningen for overflatehevningen og presenterer den endelige differansemetoden som blir brukt til numeriske beregninger og simuleringer i kapittel 3. Jeg analyserer stabiliteten til Schrödingerlikningen og beskriver dens parametere for både åpent hav og isdekt hav.

2.1 Tidsutvikling med ikke-lineær kubisk Schrödingerlikning

I denne seksjonen utleder jeg den ikke-lineære kubiske Schrödingerlikningen jeg bruker for å studere hvordan amplituden til overflatehevningen utvikler seg i tid. Jeg utleder et utrykk for overflatehevningen og presenterer den endelige differansemetoden som blir brukt i utregningene.

Følgende beskrivelse tar utgangspunkt i likning (21) og likning (22) fra Liu og Mollo-Christensen sin artikkel Wave Propagation in a Solid Ice Pack [8]. Løsningene for overflatehevingen, η(x, t), og amplituden, A(x, t) er

η(x, t) = Re[A(x, t)e

^i(k^c^x−ω^c^t)

], (2.1)

∂A

∂t + iβ ∂

²

A

∂x

²

+ iγA

²

A

^∗

= 0. (2.2)

Likning (2.1) beskriver overflatehevningen, η(x, t), som er en reel funksjon av en romlig

variabel x og tiden t. A(x, t) er en kompleks, sakte varierende funksjon som representerer

amplituden til bølgen. Konstanten k

_c

representerer det karakteristiske bølgetallet og

konstanten ω

c

representer den karakteristiske vinkelfrekvensen. Sistnevnte beskrives

nærmere i seksjon 2.3 for tilfellene som svarer til åpent og isdekt hav. Likning (2.2)

(14)

kalles den ikke-lineære kubiske Schrödingerlikningen, og er den samme som ble brukt i de numeriske beregningene utført av Liu og Mollo-Christensen. Konstantene β og γ blir også utledet i seksjon 2.3.

I Liu og Mollo-Christensen sin artikkel brukte de et koordinatsystem som beveget seg med gruppehastigheten, C

g

, til bølgen. Dette gjør noen utregninger enklere, men fysisk sett gir det lite mening, da man som regel måler overflatehevningen i et punkt og ikke løpende ved siden av bølgen. Jeg velger derfor å bruke et stasjonært koordinatsystem. For å oppnå dette må likningen (2.2) skrives om. Ved å legge til leddet C

g∂A

∂x

i likning (2.2) går jeg fra et koordinatsystem som beveger seg med gruppehastigheten til et stasjonært koordinatsystem. Den ikke-lineære Schrödingerlikningen blir

∂A

∂t = −C

_g

∂A

∂x − iβ ∂

²

A

∂x

²

− iγ|A|

²

A. (2.3)

I denne oppgaven studerer jeg et initialverdiproblem for likningen på et rektangulært domene, der 0 ≤ x ≤ L og 0 ≤ t ≤ T , med initialbetingelse

A(x, 0) = I(x) , x ∈ [0, L]. (2.4)

I rom velger jeg å bruke periodiske randbetingelser:

A(0, t) = A(L, t) , t ∈ [0, T ] (2.5)

(15)

Figur 2.1: Illustrasjon av hvordan periodiske randbetingelser fungerer. De lilla prikkene representerer punkter i tid, mens de grønne prikkene representerer punkter i x og hvordan verdiene til disse blir byttet ut ved periodiske randbetingelser [5].

Konstantene L og T er sluttverdien i henholdsvis rom og tid. Med periodiske randbetingelser blir startverdien i x satt lik sluttverdien i x slik at jeg kan simulere et helt system ved bare å simulere en liten del. Denne metoden er fremstilt i figur 2.1, der det blir vist hvordan verdiene i x blir byttet ut.

I resten av oppgaven brukes en annen likning for overflatehevningen, der likningen (2.1) blir utvidet ved perturbasjon for å inkludere 2. harmoniske bølger. Jeg får da en mer realistisk overflatehevning, enn hvis jeg kun har med 1. harmoniske bølger.

Denne perturbasjonen ble gjort i Liu og Mollo-Christensen sin artikkel (se likning (29b), (35b) og (37c) i artikkelen). Disse likningene blir skrevet om til å samsvare med et stasjonært koordinatsystem. Den fulle likningen for overflatehevningen med både 1. og 2. harmoniske bølger blir da:

η(x, t) = Re h

A(x, t)e

^i(k^c^x−ω^c^t)

i

+

²

Re h

(A(x, t))

²

η

2

e

^2i(k^c^x−ω^c^t)

i

(2.6)

η

2

er en konstant som blir utledet i seksjon 2.3, og representerer steilheten til bølgen.

(16)

2.1.1 Eksplisitt endelig differansemetode

For å løse likningen (2.3) numerisk diskretiserer jeg området slik at jeg får N

_x

punkter i rom og N

t

punkter i tid. Jeg antar et uniformt rutenett, der

0 = x

0

< x

1

< ... < x

j

< ... < x

N_x−2

< x

N_x−1

= L,

0 = t

⁰

< t

¹

< ... < t

ⁿ

< ... < t

^N^t⁻²

< t

^N^t⁻¹

= T,

og x

j

= j∆x og t

ⁿ

= n∆t, for j = 0, ..., N

x

− 1 og n = 0, ..., N

t

− 1. ∆x og ∆t er steglengden i henholdsvis rom og tid.

For å løse den ikke-lineære kubiske Schrödingerlikningen bruker jeg en eksplisitt endelig differansemetode. En ulempe ved å bruke en eksplisitt endelig differansemetode, i forhold

til implisitt metode, er at den kun er betinget stabil når man simulerer partielle differensialligninger.

Det betyr at det er en grense på hvor stort tidssteget kan være i forhold til steglengden i rom. Hvis tidssteget er større en denne grensen vil metoden bli ustabil og divergere.

Fordelen ved å bruke eksplisitt metode er at den er enkel å implementere.

I Caplan og Carretero-González sin artikkel Numerical stability of explicit Runge–Kutta finite-difference schemes for the nonlinear Schrödinger equation [1] beregnet de hvor stort tidssteget kan være i forhold til steglengden i rom for ulike endelig differansemetoder.

I denne artikkelen ble det vist at en 3. ordens Runge-Kutta metode er den laveste ordenen av Runge-Kutta som er betinget stabil for den ikke-lineære Schrödingerlikningen.

For å unngå numerisk ustabilitet velger jeg å bruke 2. ordens sentral endelig differansemetode i rom med 4. ordens Runge-Kutta metode i tid. I Caplan og Carretero-González sin artikkel viste de at denne metoden er stabil for:

∆t < (∆x)

²

√ 2β (2.7)

Jeg implementerer en 4. ordens Runge-Kutta metode på det tidsderiverte leddet i likning (2.3) og 2. ordens sentral endelig differanse på de x-deriverte leddene i likning (2.3).

∂A

∂t

|{z}

∂A

∂t

= −C

_g

∂A

∂x

|{z}

i)

−iβ ∂

²

A

∂x

²

| {z }

ii)

−iγ |A|

²

A

| {z }

iii)

| {z }

f(t,A)

Jeg deler opp likning (2.3) slik at jeg får tidsderiverte og romderiverte ledd hver for seg:

∂A

∂t = f (t, A) f (t, A) = −C

g

∂A

∂x

|{z}

i)

−iβ ∂

²

A

∂x

²

| {z }

ii)

−iγ |A|

²

A

| {z }

iii)

(2.8)

(17)

Likningen for 4. ordens Runge-Kutta metode finner jeg i kompendiet Numerical Methods for Engineers av Hellevik [6]:

k

₁

(x) = ∆t f(t

ⁿ

, A(x, t

ⁿ

)) k

2

(x) = ∆t f(t

ⁿ

+ 1

2 ∆t, A(x, t

ⁿ

) + 1 2 k

1

(x)) k

3

(x) = ∆t f(t

ⁿ

+ 1

2 ∆t, A(x, t

ⁿ

) + 1 2 k

2

(x)) k

4

(x) = ∆t f(t

ⁿ

+ ∆t, A(x, t

ⁿ

) + k

3

(x)) A(x, t

ⁿ⁺¹

) = A(x, t

ⁿ

) + 1

6 (k

1

(x) + 2k

2

(x) + 2k

3

(x) + k

4

(x))

(2.9)

Jeg implementerer 2. ordens sentral endelig differansemetode på leddene i f (t, A) i likning (2.8):

• Ledd i) i likning (2.8):

∂A

∂x ≈ A

ⁿ_i+1

− A

ⁿ_i−1

2∆x

• Ledd ii) i likning (2.8):

∂

²

A

∂x

²

≈ A

ⁿ_i+1

− 2A

ⁿ_i

+ A

ⁿ_i−1

(∆x)

²

• Ledd iii) i likning (2.8):

|A|

²

A ≈ |A

ⁿ_i

|

²

A

ⁿ_i

Jeg setter leddene sammen:

f(t, A) = −C

g

A

ⁿ_i+1

− A

ⁿ_i−1

2∆x − iβ A

ⁿ_i+1

− 2A

ⁿ_i

+ A

ⁿ_i−1

(∆x)

²

− iγ|A

ⁿ_i

|

²

A

ⁿ_i

(2.10)

2.1.2 Fysiske løsninger av Schrödingerlikningen

Det er kjent at Schrödingerlikningen (2.3) kan ha løsninger som bryter med fysiske lover.

Når Schrödingerlikningen brukes i kvantemekanikken til å beskrive bølgeegenskapene til en partikkel, så kan dette illustreres ved å betrakte en trillende ball som et eksempel.

Ifølge fysiske lover vil en ball som ruller sakte oppover en bakke komme til et stopp og rulle nedover igjen, dersom energien er for lav til å nå toppen. Allikevel forutsier Schrödingerlikningen at det eksisterer en liten sannsynlighet for at ballen vil trille over bakketoppen og til den andre siden.

For å få et resultat som gir fysisk mening, må derfor ikke-fysiske løsninger kuttes vekk.

Metoden jeg bruker for å gjennomføre dette, er at jeg for hvert tidssteg finner Fourier

transformen av amplituden, A(x, t

j

), og alle verdiene til Fourier transformen som er

utenfor et bestemt område setter jeg lik null. Området defineres utifra båndbredden

(18)

til energispekteret, σ. Energispekteret er absoluttverdien av Fourier transformen av overflatehevningen opphøyet i andre (| η| ˆ

²

). Jeg bestemmer at Fourier transformen til amplituden er null utenfor intervallet k ∈ [−(10σ − 0.1), 10σ − 0.1], der k er bølgetallet til overflatehevningen.

Et fullt Python-skript for å løse den ikke-lineære Schrödingerlikningen med 4. ordens

Runge-Kutta metode er lagt ved i vedlegg B

(19)

2.2 Romutvikling med ikke-lineær Schrödingerlikning

I denne seksjonen utleder jeg en annen numerisk metode for den ikke-lineære kubiske Schrödingerlikningen som jeg bruker til å studere hvordan amplituden til overflatehevningen utvikler seg i rom. Jeg analyserer stabiliteten til denne likningen og presenterer den endelige differansemetoden som blir brukt i utregningene.

Jeg ønsker å studere hvordan en bølge utvikler seg som et initialverdiproblem med hensyn på x. For å oppnå dette må jeg ved hjelp av perturbasjon skrive om likning (2.3) fra forrige seksjon:

∂A

∂t = −C

_g

∂A

∂x − iβ ∂

²

A

∂x

²

− iγ|A|

²

A Jeg innfører en perturbasjon ved å anta at løsningen er på formen:

A = A

0

+ A

1

(2.11)

Jeg introduserer en sakte tid og rom skala:

t

₁

= t, x

₁

= x (2.12)

Jeg setter likning (2.11) og (2.12) inn i likning (2.3):

∂A

0

+ A

1

∂

^t¹

= −C

g

∂A

0

+ A

1

∂

^x¹

− iβ ∂

²

A

0

+ A

1

∂(

^x¹

)

²

− iγ(A

0

+ A

1

)

²

(A

0

+ A

1

)

^∗

Jeg får da uttrykket:

∂A

0

∂t

1

+

²

∂A

1

∂t

1

= −C

g

( ∂A

0

∂x

1

+

²

∂A

1

∂x

1

) − iβ(

²

∂

²

A

0

∂x

²₁

+

³

∂

²

A

1

∂x

²₁

)

− iγ(A

²₀

A

^∗₀

+ A

²₀

A

^∗₁

+ 2A

0

A

^∗₀

A

1

+ 2

²

A

0

A

1

A

^∗₁

+

²

A

^∗₀

A

²₁

+

³

A

²₁

A

^∗₁

)

Med perturbasjonsteori finner jeg løsningen A

0

, ved å studere den nulte-ordens tilnærmede likningen av :

⁰

: −iγA

²₀

A

^∗₀

= 0 ⇒ A

0

= 0

Jeg kan da sette alle ledd med A

₀

lik null. Jeg skriver om A

₁

= A, t

₁

= t og x

₁

= x.

Likningen blir da:

²

∂A

∂t = −

²

C

g

∂A

∂x − i

³

β ∂

²

A

∂x

²

− i

³

γ|A|

²

A Jeg deler på

²

og skriver om uttrykket:

∂A

∂x = − 1 C

g

∂A

∂t − iβ C

g

∂

²

A

∂x

²

− iγ C

g

|A|

²

A (2.13)

For å finne et uttrykk for utvikling i rom skriver jeg om leddet

^iβ_C

g

∂²A

∂x²

i likning (2.13),

(20)

slik at jeg kun får tidsderivert ledd på høyresiden. For å oppnå dette deriverer jeg først likning (2.13) med hensyn på x. Jeg antar da at A(x, t) er tilstrekkelig glatt slik at de nødvendige partiellderiverte eksisterer og er kontinuerlige.

∂

²

A

∂x

²

= − 1 C

g

∂

²

A

∂t∂x − iβ C

g

∂

³

A

∂x

³

− iγ C

g

∂(|A|

²

A)

∂x (2.14)

Fra likning (2.13) og (2.14) får jeg da de to likningene:

i) ∂A

∂x = − 1 C

g

∂A

∂t − O() ii) ∂

²

A

∂x

²

= − 1 C

_g

∂

²

A

∂t∂x − O() = − 1 C

_g

∂

∂t

∂A

∂x − O()

Jeg substituerer så i) inn i ii):

∂

²

A

∂x

²

= − 1 C

_g

∂

∂t

− 1 C

_g

∂A

∂t − O()

− O() = 1 C

_g²

∂

²

A

∂t

²

− O()

Jeg antar at er liten og setter:

∂

²

A

∂x

²

≈ 1 C

_g²

∂

²

A

∂t

²

Jeg setter dette inn i likning (2.13):

∂A

∂x = − 1 C

g

∂A

∂t − iβ C

g

1 C

_g²

∂

²

A

∂t

²

− iγ C

g

|A|

²

A + O(

²

)

Utvikling i rom med den ikke-lineære Schrödingerlikningen blir da:

∂A

∂x = − 1 C

g

∂A

∂t − iβ C

_g³

∂

²

A

∂t

²

− iγ C

g

|A|

²

A (2.15)

På lik måte som i forrige seksjon, seksjon 2.1 studerer jeg et initialverdiproblem for likningen på et rektangulært domene, med initialverdi:

S(0, t) = I(t) , t ∈ [0, T ] (2.16) Jeg velger å bruke periodiske randbetingelser i tid:

S(x, 0) = S(x, T ) , x ∈ [0, L] (2.17)

(21)

2.2.1 Analyse av stabiliteten til likningen

Stabiliteten til likning (2.15) defineres ut ifra evnen til å motstå en forstyrrelse. Hvis en gitt fysisk tilstand kan motstå en forstyrrelse og returnere til sin opprinnelige tilstand så er den stabil. I denne oppgaven ser jeg blant annet på isdekt hav, så en ustabil fysisk tilstand kan føre til at isen brytes opp. Jeg undersøker stabiliteten til likning (2.15) ved å følge oppskriften fremstilt i kapittel 5 i boka av White, F. M. Viscous fluid flow [11]:

• 1. Finn den grunnleggende løsningen, A

0

.

• 2. Legg til en forstyrrelse, A

⁰

, og substituer A

0

+ A

⁰

inn i den grunnleggende likningen.

• 3. Finn forstyrrelseslikningen ved å trekke ifra den grunnleggende løsningen.

• 4. Lineariser forstyrrelseslikningen ved å anta at forstyrrelsen er liten og fjern alle ledd av andre eller høyre orden.

• 5. Dersom den lineariserte likningen er komplisert kan den forenkles ved å anta at forstyrrelsen har form som en bølge.

1. Jeg finner først den grunnleggende løsning, A

₀

:

Jeg antar at startverdien i x, A(0, t) = I(t) = I for 0 ≤ t ≤ T, er konstant (se likning (2.16)). Siden

^∂A_∂t

= 0 vil flere av leddene i likning (2.15) forsvinne. Jeg står da kun igjen med leddene:

∂A

∂x + iγ

C

_g

A

²

A

^∗

= 0 (2.18)

Jeg antar at løsning, A(x, t), er på formen:

A(x, t) = ae

^iθ(x,t)

Der a er en konstant og θ(x, t) er en funksjon av x og t. Jeg setter A(x, t) inn i likning (2.18) og får:

i

a ∂θ(x, t)

∂x + a

³

γ C

g

e

^iθ(x,t)

= 0

Jeg løser for θ(x, t):

θ(x, t) = −a

²

γ C

g

x + θ

₀

(t)

Jeg får da løsningen:

A(x, t) = ae

^iθ(x,t)

= ae

^i(−a²

γ Cgx+θ₀(t))

= ae

^iθ⁰^(t)

e

^−a²ⁱ

γ Cgx

Ved x = 0 vet jeg at A(0, t) = I. Jeg setter derfor x = 0 for å finne integrasjonskonstanten θ

0

(t):

A(0, t) = ae

^iθ⁰^(t)

e

⁰

= ae

^iθ⁰^(t)

= I

(22)

Den grunnlegende løsningen, A

₀

, blir da:

A

₀

(x) = Ie

^−a²ⁱ^Cg^γ^x

(2.19)

Jeg setter likning (2.19) inn i likning (2.18) og finner den grunnleggende likningen:

h − ia

²

γ C

g

+ a

²

i

A

0

= 0 (2.20)

2. Jeg legger til en liten forstyrrelse (A

⁰

): A ˆ = A

0

+ A

⁰

: Jeg definerer forstyrrelsen som:

A

⁰

= Ie

⁻ⁱ

γ Cga²x

b(x, t) + id(x, t)

= A

0

b(x, t) + id(x, t)

Summen av den grunnleggende løsningen og forstyrrelsen blir da:

A ˆ = Ie

⁻ⁱ

γ Cga²x

1 + b(x, t) + id(x, t)

= A

0

1 + b(x, t) + id(x, t)

Siden A

⁰

er en funksjon av x og t må nå alle leddene i likning (2.15) tas med. Jeg setter A ˆ inn i likning (2.15):

∂ A ˆ

∂x + 1 C

g

∂ A ˆ

∂t + iβ C

_g³

∂

²

A ˆ

∂t

²

+ iγ C

g

A ˆ

²

A ˆ

^∗

= 0 (2.21)

4. Normalt sett lineariserer man likningen etter neste steg, men for å unngå alt for lange uttrykk lineariserer jeg likningen underveis ved å anta at forstyrrelsen er veldig liten (A

0

>> A

⁰

). Alle ikke-lineære A

⁰

ledd kan dermed neglisjeres. Jeg får da likningen:

"

i ∂d

∂x − a

²

γ C

_g

1 + b +

∂b

∂x + a

²

γ C

_g

d

# + 1

C

_g

"

∂b

∂t + i ∂d

∂t

#

+ iβ C

_g³

"

∂

²

b

∂t

²

+ i ∂

²

d

∂t

²

# + iγ

C

_g

h

a

²

1 + 3b + id i

! A

0

= 0

3. Jeg trekker fra den grunnleggende likningen (2.20) og får da den lineariserte forstyrrelseslikningen:

"

i ∂d

∂x − a

²

γ C

_g

b

+

∂b

∂x + a

²

γ C

_g

d

# + 1

C

_g

"

∂b

∂t + i ∂d

∂t

#

+ iβ C

_g³

"

∂

²

b

∂t

²

+ i ∂

²

d

∂t

²

# + iγ

C

_g

h

a

²

3b + id i

!

A

0

= 0

(23)

Jeg deler opp i reell og imaginær del:

Re : ∂b(x, t)

∂x + 1 C

g

∂b(x, t)

∂t − β C

_g³

∂

²

d(x, t)

∂t

²

= 0 Im : ∂d(x, t)

∂x + 1 C

g

∂d(x, t)

∂t + β C

_g³

∂

²

b(x, t)

∂t

²

+ a

²

γ C

g

2b(x, t) = 0

(2.22)

5. Jeg forenkler uttrykket ved å anta at forstyrrelsen har form som en bølge b = ˆ be

^i(Kx−Ωt)

og d = ˆ de

^i(Kx−Ωt)

. Der K er en konstant som representerer bølgetallet og Ω er en konstant som representerer vinkelfrekvensen. ˆ b, d ˆ er konstanter. Jeg setter dette inn i likningen (2.22) og får:

Re : h

ˆ b(iK)e

^i(Kx−Ωt)

i + 1

C

g

h ˆ b(−iΩ)e

^i(Kx−Ωt)

i

− β C

_g³

h d(−iΩ) ˆ

²

e

^i(Kx−Ωt)

i

= 0

Im : h

d(iK)e ˆ

^i(Kx−Ωt)

i + 1

C

g

h d(−iΩ)e ˆ

^i(Kx−Ωt)

i + β

C

_g³

h ˆ b(−iΩ)

²

e

^i(Kx−Ωt)

i

+ a

²

γ C

g

2 h

ˆ be

^i(Kx−Ωt)

i

= 0

Dette kan skrives som likningssystemet:



 i

K −

_C¹

g

Ω

C_g

_β

C_g²

Ω

²

C_g

2a

²

γ −

_C^β2 g

Ω

²

i K −

_C¹

g

Ω



 ·

" ˆ b d ˆ

#

= 0

Hvis dette skal ha en løsning må determinanten til matrisen være 0:

det



 i

K −

_C¹

g

Ω

Cg

_β

C_g²

Ω

²

Cg

2a

²

γ −

_C^β2 g

Ω

²

i

K −

_C¹

g

Ω



 = 0

Determinanten er:

det



 i

K −

_C¹

g

Ω

C_g

β C_g²

Ω

²

C_g

2a

²

γ −

_C^β2 g

Ω

²

i K −

_C¹

g

Ω



 = i K− 1

C

g

Ω i

K− 1 C

g

Ω

− C

g

β C

_g²

Ω

²

C

g

2a

²

γ− β C

_g²

Ω

²

= 0

Jeg løser for K og får at:

K

²

− 2

C

_g

ΩK + Ω

²

C

_g²

1 +

²

C

_g²

2a

²

βγ − β

²

Ω

²

C

_g²

= 0

Jeg løser andregradslikningen for K(Ω):

K(Ω) = Ω C

g

1 ± C

g

s β

²

Ω

²

C

_g²

− 2a

²

βγ

!

(2.23)

Ligningen er stabil dersom innsatt Ω er reell gir at K også er reell. For at likningen skal være stabil må da:

Ω > √ 2aC

g

r γ

β (2.24)

(24)

Dette er i overensstemmelse med tilsvarende resultat av Liu og Mollo-Christensen [8], som viste at likningen 2.2 er stabil dersom:

K > √ 2a

r γ β

Siden de ser på utviklingen i tid har de en betingelse for K i stedet for Ω, og siden de har et koordinatsystem som beveger seg med gruppehastigheten, C

_g

, er denne parameteren ikke med i deres uttrykk.

2.2.2 Analyse av resultatet fra stabilitetsanalysen

Jeg analyserer resultatet fra stabilitetsanalysen ved å plotte vekstraten. Vekstraten angir hvor stor den relative endringen av Im[K(Ω)] over Ω er, se likning (2.23). Jeg setter inn konstantene β , γ, C

g

og for åpent hav og for isdekt hav. Disse konstantene blir utledet i neste seksjon, seksjon 2.3. a setter jeg lik 1. Jeg ønsker å studere området hvor likningen er ustabil fordi det her her veksten er størst. Jeg ser derfor på intervallet Ω ∈ [− √

2C

g

q

_γ

β

, √ 2C

g

q

_γ

β

].

(a) Vekstraten for åpent hav (b) Vekstraten for isdekt hav

Figur 2.2: Forskjellen i vekstraten, den relative endringen av bølgetallet, Im[K], over vinkelfrekvensen, Ω, mellom åpent og isdekt hav.

For åpent hav finner jeg at den maksimale ustabiliteten inntreffer ved Ω

maks

= 8.7 s

⁻¹

med den maksimale vekstraten (ImK)

_max

= 0.8 m

⁻¹

. For isdekt hav inntreffer den maksimale ustabiliteten ved Ω

max

= 15.1 s

⁻¹

. Da er den maksimale vekstraten (ImK)

max

= 7.9 m

⁻¹

. Formasjonen av bølgepakker skjer mye raskere i is en i åpent hav. Dette resultatet stemmer overens med det de finner i artikkelen til Liu og Mollo-Christensen [8].

I Liu og Mollo-Christensen sin artikkel kommer de frem til at den større vekstrate som blir funnet for isdekt hav antyder at formasjonen av bølgepakker i et område med høy trykkspenning kan være med på å bryte opp isen.

Et fullt Python-skript for utregning av vekstraten er lagt ved i vedlegg A.

(25)

2.2.3 Endelig differansemetode

For å løse likningen (2.15) numerisk diskretiserer jeg et rektangulært området på samme måte som i forrige seksjon (seksjon 2.1). Jeg implementerer en 4. ordens Runge-Kutta metode på det x-deriverte leddet i likning (2.15) og 2. ordens sentral endelig differanse på de tidsderiverte leddene i likning (2.15). Fra artikkelen til Caplan og Carretero-González vet jeg at denne metoden er stabil for [1]:

∆x < (∆t)

²

√ 2

_C^β3 g

(2.25)

Jeg deler opp likning (2.15) slik at jeg får tidsderiverte og romderiverte ledd hver for seg:

∂A

∂x

|{z}

∂A

∂x

= − 1 C

g

∂A

∂t

|{z}

i)

− iβ C

_g³

∂

²

A

∂t

²

| {z }

ii)

− iγ C

g

|A|

²

A

| {z }

iii)

| {z }

f(x,A)

∂A

∂x = f (x, A) f (x, A) = − 1

C

g

∂A

∂t

|{z}

i)

− iβ C

_g³

∂

²

A

∂t

²

| {z }

ii)

− iγ C

g

|A|

²

A

| {z }

iii)

(2.26)

Jeg bruker samme likning for 4. ordens Runge-kutta metode som i forrige seksjon.

k

1

(t) = ∆x f (x

j

, A(x

j

, t)) k

2

(t) = ∆x f (x

j

+ 1

2 ∆x, A(x

j

, t) + 1 2 k

1

(t)) k

3

(t) = ∆x f (x

j

+ 1

2 ∆x, A(x

j

, t) + 1 2 k

2

(t)) k

4

(t) = ∆x f (x

j

+ ∆x, A(x

j

, t) + k

3

(t)) A(x

_j+1

, t) = A(x

_j

, t) + 1

6 (k

₁

(t) + 2k

₂

(t) + 2k

₃

(t) + k

₄

(t))

(2.27)

Jeg implementerer 2. ordens sentral endelig differansemetode på leddene i f(x, A) i likning (2.26):

• Ledd i) i likning (2.26):

∂A

∂t ≈ A

ⁿ⁺¹_i

− A

ⁿ⁻¹_i

2∆t

• Ledd ii) i likning (2.26):

∂

²

A

∂t

²

≈ A

ⁿ⁺¹_i

− 2A

ⁿ_i

+ A

ⁿ⁻¹_i

(∆t)

²

• Ledd iii) i likning (2.26):

|A|

²

A ≈ |A

ⁿ_i

|

²

A

ⁿ_i

(26)

Jeg setter leddene sammen:

f (x, A) = − 1 C

g

A

ⁿ⁺¹_i

− A

ⁿ⁻¹_i

2∆t − iβ

C

_g³

A

ⁿ⁺¹_i

− 2A

ⁿ_i

+ A

ⁿ⁻¹_i

(∆t)

²

− iγ

Cg |A

ⁿ_i

|

²

A

ⁿ_i

(2.28)

Et fullt Python-skript for å løse den ikke-lineære Schrödingerlikningen med 4. ordens

Runge-Kutta metode er lagt ved i vedlegg B

(27)

2.3 Parametere i Schrödingerlikningen

I denne seksjonen beskrives parametrene, C

_g

, β og γ, i Schrödingerlikningen for tid- og romutvikling, likning (2.3) og (2.15). Parametrene beskrives for både åpent og isdekt hav.

I artikkelen til Liu og Mollo-Christensens regner de ut parametrene i Schrödingerlikningen [8]. Gruppehastigheten, C

g

, finner de ved å derivere vinkelfrekvensen, ω (kvadratroten av dispersjonsrelasjonen), på bølgetallet, k: C

g

=

^∂ω_∂k

. Parameteren β er den dobbeltderiverte av vinkelfrekvensen med hensyn på k, mens parameteren γ blir utledet ved perturbasjon.

2.3.1 Dispersjonsrelasjonen

Uttrykket for dispersjonsrelasjonen henter jeg fra Liu og Mollo-Christensen sin artikkel (se likning (6) i artikkelen):

ω

²

(k) = gk + Bk

⁵

− Qk

³

1 + kM B = Eh

³

12(1 − s

²

)ρ

H

, Q = P h ρ

H

, M = ρ

_I

h ρ

H

(2.29)

Parameteren B er en konstant som beskriver effektene av bøyning av isen. Konstanten Q angir effektene av isens kompresjon og konstanten M angir effektene av tregheten til isen. k er bølgetallet, E er Young’s modul, s er Poisson’s ratio, ρ

I

er tettheten til isen og ρ

H

er tettheten til havet.

Når jeg modellerer for åpent hav, kan jeg sette tykkelsen til isen lik null, h = 0.

Konstantene B, Q, og M blir da null, og dispersjonsrelasjonen for åpent hav blir:

ω

²

(k) = gk (2.30)

Dette er dispersjonsrelasjonen for dypt vann (kH >> 1, der k er bølgetallet og H er havdypet), der tyngdekraften er mer dominerende enn overflatespenningen for bølgebevegelsen.

2.3.2 Gruppehastigheten og fasefarten

Gruppehastigheten er definert som den k-deriverte av vinkelfrekvensen, ω (kvadratroten av dispersjonsrelasjonen). Denne utregningen blir gjort i Liu og Mollo-Christensen sin artikkel, se likning (11) i artikkelen.

C

g

= ∂ω

∂k = g + Bk

⁴

(5 + 4kM) − Qk

²

(3 + 2kM )

2ω(1 + kM )

²

(2.31)

(28)

Fasefarten er definert som vinkelfrekvensen, ω delt på bølgetallet k:

C

_p

= ω k =

s

gk + Bk

⁵

− Qk

³

k

²

(1 + kM) (2.32)

Gruppehastighet og fasefarten for åpent hav, der h = 0, blir da:

C

g

= ∂ω

∂k = g 2ω = 2

r g

k (2.33)

C

p

= ω k =

r g k = 1

2 C

g

(2.34)

2.3.3 Parameteren β

Parameteren β er definert som den dobbeltderiverte av vinkelfrekvensen, ω, med hensyn på bølgetallet, k. Denne parameteren er utledet i likning (41) i Liu og Mollo-Christensen sin artikkel.

β = − 1 2

∂

²

ω

∂k

²

= g

²

(1 + 4kM) − 6gBk

⁴

(5 + 8kM + 4k

²

M

²

) + 2gQk

²

(3 + 2kM + 2k

²

M

²

) 8ω

³

(1 + kM)

⁴

−3Q

²

k

⁴

− B

²

k

⁸

(15 − 20kM + 8k

²

M

²

) 8ω

³

(1 + kM)

⁴

(2.35)

For åpent hav blir konstantene B, Q og M satt lik 0 og jeg får:

β = 1 8

ω

k

²

(2.36)

2.3.4 Parameteren γ

Parameteren γ blir funnet ved perturbasjon av hastighetspotensialet φ, overflatehevningen

η og fasefarten C

_p

. Parameteren blir utledet i likning (39) og (40) i Liu og Mollo-Christensen

(29)

sin artikkel [8]:

γ = C

2

k Der:

C

2

= k

²

C

0

(g + 127Bk

⁴

− 10Q

0

k

²

− 3Qk

²

) − 36kη

2

C

0

(10Bk

⁴

− Q

0

k

²

) − 2kφ

2

(7g + 187Bk

⁴

− 25Q

0

k

²

) 2g

0

(2.37) Konstantene i likning (2.37) er utledet i likning (33), (37b) og (37c) i artikkelen til Liu og Mollo-Christensen [8].

Q

₀

= Q + M C

²

g

₀

= g + Bk

⁴

− Q

₀

k

²

C

0

= r g

0

k

η

2

= − g

0

k(−g + 45Bk

⁴

− 9Q

0

k

²

)

2(g + 15Bk

⁴

− 3Q

₀

k

²

)(g + 16Bk

⁴

− 4Q

₀

k

²

) (2.38) φ

2

= kC

0

(15Bk

⁴

− 3Q

0

k

²

)

g + 15Bk

⁴

− 3Q

₀

k

²

Likning (2.38) for η

2

blir brukt i den fulle likningen for overflatehevningen (se likning (2.6) i seksjon 2.1) som inkluderer 1. og 2. harmoniske bølger:

η(x, t) = Re h

A(x, t)e

^i(k^c^x−ω^c^t)

i

+

²

Re h

(A(x, t))

²

η

2

e

^2i(k^c^x−ω^c^t)

i

I åpent hav blir konstantene B, Q og M satt lik 0 og jeg får:

γ = 1

2 ωk

²

(2.39)

η

2

= k

2 (2.40)

(30)

(31)

KAPITTEL 3

Resultater

I dette kapittelet presenterer jeg resultatene fra utregningen av den ikke-lineære Schrödingerlikningen for utvikling i både tid og rom, og med innsatte verdier for både åpent hav og isdekt

hav. Blant resultatene som presenteres er overflatehevningen fra simuleringer og ulike ulike statistiske egenskaper. Seksjon 3.1 inneholder resultater for åpent hav og seksjon 3.2 for isdekt hav.

3.1 Rom- og tidsutvikling av overflatehevningen i åpent hav

I denne seksjonen legger jeg frem resultatene fra utregningen av den ikke-lineære Schrödingerlikningen for utvikling i både tid og rom, og med innsatte verdier for åpent hav. Jeg sammenligner

mine data med tidligere studier. Resultatene som blir presentert er utviklingen av karakteristiske bølgeparametere, overflatehevningen, energispekteret, variansen, kurtosen, skjevheten og bølgetall-frekvensspekteret.

For å undersøke om mitt program produserer realistiske resultater, tester jeg det imot eksisterende studie utført i artikkelen Evolution of a narrow-band spectrum of random surface gravity waves av Dysthe, Trulsen, Krogstad og Socquet-Juglard [2]. Jeg får tilgang på programmet som ble brukt for å produsere resultatene i denne artikkelen (NLS-programmet), og som gjør det lettere å sammenligne. I utregningene, som ble gjort i artikkelen, brukte de samme Schrödingerlikning som jeg utledet i forrige kapittel, likning (2.3) i seksjon 2.1 (se likning (1) i artikkelen): Ikke-lineær kubisk Schrödingerlikning for utvikling i tid

∂A

∂t = −C

g

∂A

∂x − iβ ∂

²

A

∂x

²

− iγ|A|

²

A

For utvikling i rom bruker jeg likningen (2.15) som ble utledet tidligere i seksjon 2.2:

∂A

∂x = − 1 C

g

∂A

∂t − iβ C

_g³

∂

²

A

∂t

²

− iγ C

g

|A|

²

A

(32)

Her bruker jeg parametrene, C

_g

, β og γ, for åpent hav, som ble utledet i forrige kapittel, i seksjon 2.3.

I artikkelen av Dysthe, Trulsen, Krogstad og Socquet-Juglard har de definert overflatehevningen som (se likning (2) i artikkelen):

η(x, t) = 1 2

A(x, t)e

^i(k^c^x−ω^c^t)

+ A

2

(x, t)e

^2i(k^c^x−ω^c^t)

Der utrykket for den 2. harmoniske bølgen er:

A

₂

(x, t) = (A(x, t))

²

2 − iA(x, t) 2

∂A(x, t)

∂x

Dette er en litt annen definisjon enn den jeg bruker i mine utregninger (se likning (2.6) i seksjon 2.1). For å tilnærme meg deres resultater deler jeg funksjonen A(x, t) i likning (2.6) på i når jeg utvikler i tid og på 3.5 når jeg utvikler i rom.

η(x, t) = Re

"

A(x, t) eller 3.5

e

^i(k^c^x−ω^c^t)

# +

²

Re

"

(A(x, t))

²

eller 3.5

η

₂

e

^2i(k^c^x−ω^c^t)

)

# (3.1)

Jeg legger inn de samme konstantene i begge programmene og er påpasselig med at

steglengden i tid ikke overskrider det kravet som ble utledet i likning (2.7) i seksjon

2.1.1, når jeg ser på utviklingen i tid. For utvikling i rom passer jeg på at steglengden i

x ikke overskrider kravet utledet i likning (2.25) i seksjon 2.2.3 i forrige kapittel.

(33)

Tabell 3.1: Innsatte verdier for utvikling i tid og rom på åpent hav

Parameter Verdi i tidsutvikling Verdi i romutvikling

Antall kjøringer N 50 50

Antall steg i tid N

t

1000 2000

Antall steg i rom N

x

2000 1000

Sluttverdi i tid T 315 315

Sluttverdi i rom L 315 315

Steglengde i t ∆t 0.158 < 0.56 0.315

Steglengde i x ∆x 0.315 0.158 < 2.8

Gravitasjons akselerasjonen g 1.0 1.0

Karakteristisk bølgetall k

c

1.0 1.0

Karakteristisk vinkelfrekvens ω

_c

1.0 1.0

Karakteristisk amplitude a

c

1.0 1.0

Steilheten (a

_c

k

_c

) 0.1 0.1

Gruppehastigheten C

g

0.5 0.5

Parameter i Schrödingerlikningen β 0.125 0.125

Parameter i Schrödingerlikningen γ 0.5 0.5

Antall perioder i startverdien N

P

50.1 50.1

Antall bølgelengder i startverdien N

λ

50.1 50.1

Diskretisering av frekvensplanet ∆N

P

0.02 0.02

Diskretisering av bølgetallplanet ∆N

P

0.02 0.02

Båndbredden σ 0.1 0.1

Jeg ønsker at energispekteret skal ha form som en Gaussisk bjelle. Jeg velger derfor at startverdien i initialverdiproblemet (se likning (2.4) og likning (2.16) i henholdsvis seksjon 2.1 og seksjon 2.2) skal være kvadratroten av den inverse Fourier transformen av en Gaussisk bjelle. I artikkelen av Dysthe, Trulsen, Krogstad og Socquet-Juglard har de satt startverdien til å være (se likning (7) i artikkelen):

I(x) = F

⁻¹

s ∆N

_λ

√ 2πσ e

⁻^k

2 4σ2

e

^iθ

Der θ er en stokastisk variabel som er uniformt fordelt på intervallet [0, 2π], og som har like mange punkter som x i tidsutvikling og som t i romutvikling. For hver gang programmet kjøres vil jeg dermed få litt ulikt resultat. Ved å utføre flere kjøringer (her N = 50) kan jeg finne gjennomsnittet av verdiene, som gir et sikrere resultat. For romutvikling byttes k og x ut med ω og t

I mitt program bruker jeg en litt annen startverdi for å tilnærme meg resultatene i artikkelen. For utvikling i tid setter jeg startverdien i t = 0 lik:

I(x) = F

⁻¹

1

²

s

∆N

λ

√ 2πσ e

⁻^k

2 4σ2

e

^iθ

(34)

For utvikling i rom setter jeg startverdien i x = 0 lik:

I(t) = F

⁻¹

s

N

_P²

√ 2πσ e

⁻^ω

2 4σ2

e

^iθ

3.1.1 Måling av karakteristiske bølgeparametere

For å undersøke egenskapene til bølgefeltet og hvordan disse egenskapene utvikler seg i tid og rom måler jeg ulike bølgeparametere på starten av simuleringen og på slutten av simuleringen, se figur 3.1.

(a) Tidsutvikling (b) Romutviling

Figur 3.1: Illustrasjon av hvordan en bølge vil utvikle seg fra innsatt initialverdi til slutten av simuleringen i tidsutvikling og i romutvikling.

Parametrene gruppehastigheten, C

g

, og fasefarten, C

p

, finner jeg med likningene (2.33) og (2.34) fra seksjon 2.3. Jeg finner antall bølgetopper, N

topp

og den gjennomsnittlige amplituden, a, til disse bølgetoppene med Python-funksjonen find_peaks. Fra disse parametrene kan jeg regne ut: bølgelengden, λ, perioden, P , bølgetallet, k, frekvensen, f og vinkelfrekvensen, ω.

λ = L N

topp

, P = λ C

p

, k = 2π

λ , f = 1

P , ω = 2πf

I tabell 3.2 og tabell 3.3 er det satt inn gjennomsnittsverdien av parametrene etter

N = 50 kjøringer av begge programmene.

(35)

Tabell 3.2: Utviklingen av karakteristiske bølgeparametere i tid. Kolonne 1 viser parametrene, der 0 indikerer starten på simuleringen og T slutten på simuleringen.

Kolonne 2 viser verdiene som ble produsert med mitt program, og kolonne 3 viser utviklingen av parameteren fra starten av simuleringen til slutten av simuleringen i prosent (k

T

/k

0

) . Kolonne 4 og 5 har samme oppsett som 2 og 3 bare for verdiene fra NLS-programmet.

Parameter Mine data Utvikling Data fra NLS-programmet Utvikling

C

_g

0.507 0.507

C

p

1.013 1.013

k

₀

0.995 0.992

k

T

0.956 - 3.9 % 0.948 - 4.4 %

ω

₀

1.009 1.005

ω

T

0.969 - 4.0 % 0.961 - 4.4 %

f

₀

0.161 0.160

f

T

0.154 - 4.3 % 0.153 - 4.4 %

a

₀

0.096 0.096

a

T

0.092 - 4.2 % 0.091 - 5.2 %

0

0.096 0.095

T

0.088 - 8.3 % 0.087 - 8.4 %

λ

0

6.313 6.335

λ

T

6.573 + 4.1 % 6.626 + 4.6 %

P

0

6.229 6.252

P

_T

6.486 + 4.1 % 6.538 + 4.6 %

Resultatene fra de to forskjellige programmene samsvarer godt, og det er kun små forskjeller mellom dem. Generelt synker parameterverdiene ettersom simuleringen utvikler seg fra t = 0 til t = T , med unntak av bølgelengden, λ og perioden, P som øker.

Parametrene regnet ut med NLS-programmet synker/øker jevnt over mer enn mine data,

men forskjellen er innenfor 1%.

(36)

Tabell 3.3: Utviklingen av karakteristiske bølgeparametere i rom. Kolonne 1 viser parametrene, der 0 indikerer starten på simuleringen og L slutten på simuleringen.

Kolonne 2 viser verdiene som ble produsert med mitt program, og kolonne 3 viser utviklingen av parameteren fra starten av simuleringen til slutten av simuleringen i prosent (k

L

/k

0

) . Kolonne 4 og 5 har samme oppsett som 2 og 3 bare for verdiene fra NLS-programmet.

Parameter Mine data Utvikling Data fra NLS-programmet Utvikling

C

_g

0.502 0.502

C

p

1.003 1.003

k

₀

1.002 0.996

k

L

0.977 - 2.5 % 0.968 - 2.8 %

ω

₀

0.999 0.993

ω

L

0.973 - 2.6 % 0.965 - 2.8 %

f

₀

0.159 0.158

f

L

0.155 - 2.5 % 0.154 - 2.5 %

a

₀

0.097 0.096

a

L

0.092 - 5.2 % 0.094 - 2.1 %

0

0.097 0.096

L

0.090 - 7.2 % 0.091 - 5.2 %

λ

0

6.269 6.309

λ

L

6.434 + 2.6 % 6.489 + 2.9 %

P

0

6.290 6.330

P

_L

6.455 + 2.6 % 6.511 + 2.9 %

Som for utvikling i tid er det god overensstemmelse mellom dataene fra de to programmene.

Man kan se at parameterverdiene generelt øker/stiger mindre enn de gjør i tidsutviklingen.

Det eneste unntaket er amplituden, a, som har en høyre nedgang i mine utregninger for romutvikling enn for tidsutvikling.

3.1.2 Overflatehevning

Figur 3.2 og figur 3.3 viser overflatehevningen ved starten av simuleringen og ved slutten av en simulering, utregnet med begge programmene i henholdsvis tid- og romutvikling.

Lengden og tiden er dimensjonsløs, da jeg har multiplisert med det karakteristiske bølgetallet, k

c

, i tidsutvikling og den karakteristiske vinkelfrekvensen, ω

c

i romutvikling.

Her vises kun resultatene etter én kjøring (N = 1).

(37)

(a) Data fra NLS-programmet, t = 0 (b) Mine data, t = 0

(c) Data fra NLS-programmet, t = T (d) Mine data, t = T

Figur 3.2: Sammenligning av tidsutviklingen av overflatehevning over x fra t = 0 til t = T mellom de to datasettene

(a) Data fra NLS-programmet, x = 0 (b) Mine data, x=0

(c) Data fra NLS-programmet, x = L (d) Mine data, x=L

Figur 3.3: Sammenligning av romutviklingen av overflatehevning over t fra x = 0 til x = L mellom de to datasettene

Resultatene fra de to programmene samsvarer godt både i tidsutvikling og i romutvikling.

Man ser at amplituden, a, er omtrent lik for mine data og datane fra NLS-programmet.

Dette kan man også lese utifra tabellene 3.2 og 3.3. Man kan også se at antall bølgetopper

i hvert bølgetog er omtrent lik for de to programmene. Ved slutten av simuleringen er

bølgegruppene mindre enn ved starten av simuleringen, som tilsier at båndbredden, σ,

(38)

har økt.

Man kan også legge merke til at ved slutten av simuleringen er overflaten mer irregulær enn ved starten av simuleringen, se figur 3.4.

(a) Regulær ekvivalentsjø, (b) Irregulær langkammet sjø

(c) Kortkammet sjø

Figur 3.4: Sammenlikning av ulike havoverflater, for utvikling i x og y [3].

3.1.3 Energispektrum

Figur 3.5 og figur 3.6 viser bølgetallspektrumet ved fire ulike tider og frekvensspektrumet ved fire ulike punkter i x. Bølgetallspektrum og frekvensspektrum er navn på energispektrum over henholdsvis bølgetallet, k, og vinkelfrekvensen, ω. Energispekteret finner jeg ved å ta absoluttverdien av Fourier transformen av overflatehevningen opphøyd i andre (|ˆ η|

²

).

Her har jeg plottet det ensidige spektrumet, som kun tar med de reelle verdiene. De negative bølgetallene og de negative vinkelfrekvensene er blitt kuttet vekk. Siden de negative verdiene tilsvarer den kompleks konjugerte av de positive verdiene, er det ikke nødvendig å ha dem med. I et fullt spektrum ville spekteret vært speilet om origo.

Bølgetallet og vinkelfrekvensen er dimensjonsløs, da jeg har delt på det karakteristiske

bølgetallet k

c

i bølgetallspektrumet og den karakteristiske vinkelfrekvensen ω

c

i frekvensspekteret.

Figurene under viser ensemble gjennomsnittet av de 50 ulike energispektrene funnet etter

50 kjøringer (N = 50). Som i artikkelen av Dysthe, Trulsen, Krogstad og Socquet-Juglard,

har gjennomsnitts-spekteret deretter blitt glattet med Python-funksjonen convolve.

(39)

(a) Data fra NLS-programmet (b) Mine data

Figur 3.5: Tidsutvikling av glattet gjennomsnittlig ensidig bølgetallspektrum ved fire tidspunkter. Sammenligning mellom de to datasettene

(a) Data fra NLS-programmet (b) Mine data

Figur 3.6: Romutvikling av glattet gjennomsnittlig ensidig frekvensspekter ved fire punkter i rom. Sammenligning mellom de to datasettene

I bølgetallspekteret (figur 3.5) ser man tydelig at båndbredden, σ, øker ettersom tiden går. Båndbredden går fra omtrent σ

₀

= 0.15, til omtrent σ

_T

= 0.2. Dette stemmer overens med det man ser i utviklingen av overflatehevningen (se figur 3.2). Ettersom tiden går synker spekteret, men siden båndbredden øker vil energien holde seg tilnærmet konstant.

Veileder Karsten Trulsen informerer om at det kan være noen feil i romutviklingen til NLS-programmet. Jeg velger derfor å fokusere på mine resultater. Som i bølgetallspekteret ser man i frekvensspekteret at båndbredden, σ, øker ettersom x øker. Båndbredden går fra omtrent σ

0

= 0.15, til omtrent σ

L

= 0.28.

Man kan også legge merke til at spektrene har samme maksverdien ved starten av simuleringen, 1.7 · 10

⁻⁴

[J Hz

⁻¹

].

Figur 3.7 og figur 3.8 viser det semilogaritmiske bølgetallspektrumet ved fire ulike tider

og frekvensspektrumet ved fire ulike punkter i x. Ved bruk av logaritmisk skala blir

spekteret forstørret og variasjonene blir fremhevet.

(40)

(a) Data fra NLS-programmet (b) Mine data

Figur 3.7: Tidsutvikling av glattet gjennomsnittlig semilogaritmisk ensidig bølgetallspektrum ved fire tidspunkter. Sammenligning mellom de to datasettene

(a) Data fra NLS-programmet (b) Mine data

Figur 3.8: Romutvikling av glattet gjennomsnittlig semilogaritmisk ensidig frekvensspektrum ved fire punkter i rom. Sammenligning mellom de to datasettene

I figurene 3.7 og 3.8 kan man tydelig se at energispektrene har to topper. Dette kom ikke tydelig frem i energispektrene med lineær skala (se figur 3.5 og figur 3.6). De to toppene oppstår fordi overflatehevningen inneholder 2. ordens harmoniske bølger.

Toppene i bølgetallspekteret og frekvensspekteret befinner seg ved henholdsvis k = 1 og k = 2, og ω = 1 og ω = 2, som stemmer godt med at den 2. harmoniske bølgen har dobbelt så stor fase (2i(k

c

x − ω

c

t)) som den 1. harmoniske bølgen (i(k

c

x − ω

c

t)), se likning (2.6) i seksjon 2.1.

I figur 3.8b kan man se en skarp endring av verdiene for energispektrene x = L/2 og

x = L ved ω/ω

c

= 2. Dette skyldes at metoden jeg bruker for å kutte vekk ikke-fysiske

verdier fra løsninger av Schrödingerlikningen som forklart i seksjon 2.1, kutter nettopp

her.

(41)

3.1.4 Energi og varians

Figur 3.9 og figur 3.10 viser hvordan energien og variansen utvikler seg i tid og rom.

Energien finner jeg ved å ta integralet av energispekteret. Variansen, som gir et mål på hvor langt fra gjennomsnittet målingen er, finner jeg med Python-funksjonen var av overflatehevningen.

(a) Data fra NLS-programmet (b) Mine data

Figur 3.9: Tidsutvikling av gjennomsnittlig energi ( R

S(k)dk) og gjennomsnittlig varians (V ar(η(x, t

ⁿ

)))

(a) Data fra NLS-programmet (b) Mine data

Figur 3.10: Romutvikling av gjennomsnittlig energi ( R

S(ω)dω) og gjennomsnittlig varians (V ar(η(x

_j

, t)))

Man kan se at energien og variansen er lik, fordi de to plottene i figur 3.9 og figur 3.10 er identiske. Det betyr at likningen er normalisert riktig ( R

S(ω)dω = V ar(η(x

j

, t))).

Differansen av energien ved starten av simuleringen til slutten av simuleringen er omtrent

2.5 · 10

⁻⁵

for begge programmene i tidsutvikling, og omtrent −6.22 · 10

⁻⁵

i utregningen

av mine data og 1.6 · 10

⁻⁵

i dataene fra NLS-programmet. Forandringen er såpass liten

at energien er tilnærmet konstant og kan da sies å være bevart. Dette stemmer overens

med det man ser i figurene 3.5 og 3.6

(42)

3.1.5 Kurtose

Figur 3.11 og figur 3.12 viser hvordan kurtosen til overflatehevningen utvikler seg i tid og rom. Kurtosen er et mål på fordelingen av data rundt gjennomsnittet. Kurtosen ble funnet med Python-funksjonen kurtosis.

(a) Data fra NLS-programmet (b) Mine data

Figur 3.11: Tidsutvikling av den gjennomsnittlige kurtosen av overflatehevningen

(a) Data fra NLS-programmet (b) Mine data

Figur 3.12: Romutvikling av gjennomsnittlig kurtosen av overflatehevningen En normaltfordelt eller Gaussisk overflatehevning, η(x, t), har kurtoseverdi på 3. I figurene 3.11 og 3.12 starter kurtoseverdien på omtrent 3 og ettersom simuleringen går, øker kurtoseverdien. De høyeste verdiene er oppe i 6-6.5. Når kurtoseverdien er høyere enn 3, er det større sannsynlighet for at ekstreme bølger kan forekommer. Kurtoseverdien er noe lavere i romutviklingen enn i tidsutviklingen

Kurtose er også et mål på spissheten"og på tykkelsen på halen til fordelingen. I figurene 3.7 og 3.8 øker tykkelse til halen ettersom tiden går og x øker, og derfor øker også kurtosen.

I artikkelen Influence of crest and group length on the occurrence of freak waves av

Gramstad og Trulsen finner de kurtosen som en funksjon av BFI (Benjamin–Feir index,

lengden på en bølgegruppe) og kamlengden, L

_c

, se figur 3.13 [4].

(43)

Figur 3.13: Kurtosen som en funksjon av BFI og kamlengden, med Gaussisk initial spektrum.

BF I = σ

ω

/ω

c

, L

c

= k

c

σ

k

Der σ

ω

og σ

k

er båndbredden i henholdsvis frekvensspekteret og bølgetallspekteret.

Kurtoseverdiene jeg får ligger langs x-aksen, 1/L

c

= 0, i figur 3.13. Det betyr at bølgene er ideelt langkammede (se figur 3.4). I utviklingen av overflatehevningen, se figur 3.2 og figur 3.3, ser man at overflatehevningen blir mer irregulær mot slutten av simuleringen.