\documentclass[12pt,a4paper]{article}
\usepackage{amssymb,vmargin}
\usepackage[norsk]{babel}
\usepackage{bm}
\usepackage{enumerate}
\usepackage{graphicx}
\usepackage{hyperref}
\usepackage[utf8]{inputenc}
\newcommand{\rmd}{{\mathrm d}}
\newcommand{\oder}[2]{\frac{{\rm d}#1}{{\rm d}#2}}
\newcommand{\pd}[1]{\frac{\partial}{\partial #1}}
\newcommand{\pdt}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\pddt}[2]{\frac{\partial^2 #1}{\partial #2^2}}
\newcommand{\MatL}[1]{\texttt{\bf #1}}

\newcounter{oppg}
\newcounter{deloppg}

\newcommand{\oppgave}[1]{\goodbreak\vspace{2mm}\setcounter{deloppg}{0}\addtocounter{oppg}{1}\vspace{0mm}\noindent {\bf Oppg.~\arabic{oppg} {\em #1}.\quad }}
\newcommand{\eoppgave}{\nobreak\vspace{1mm}\nobreak}

\newcommand{\deloppgave}{\goodbreak\vspace{1mm}\addtocounter{deloppg}{1}\vspace{0mm}\noindent {\bf \alph{deloppg})\ }}

\begin{document}


\title*{\Huge Oblig 2 MEK1100, vår 2016}
\date{}
\section*{Krav til innlevering og godkjenning}

\begin{sloppypar}
  Hvert punkt gir makismalt 10 poeng.  I alt kan du oppnå 100 poeng.  Vi
  krever minimum 70 prosent, eller 70 poeng for å få godkjent.  Dersom dette
  ikke er innfridd ved innlevering, men besvarelsen vurderes som et seriøst
  forsøk, kan det gis anledning til ny innlevering.  For informasjon om
  regler se
  \href{http://www.uio.no/studier/admin/obligatoriske-aktiviteter/mn-math-oblig.html}{\texttt{http://www.uio.no/studier/admin/obligatoriske-aktiviteter/mn-math-oblig.html}}\\
  Tidsfrister er gitt på semestersidene til MEK1100.
\end{sloppypar}

Du kan bruke \LaTeX\ eller levere en håndskrevet besvarelse.  Ta med
autentiske listinger av alle skript og funksjoner i Matlab eller Python, samt
de kjøreeksemplene det bes om.  I utgangspunktet skal alle skript kjøre, men
du kan få noe utteling for ufullstendige forsøk som ikke gir rett resultat.
Kjøreeksempler som ikke kan framkomme som resultat av vedlagte skript, gir
ikke poeng.

Det skal stå filnavn på alle listinger av skript og funksjoner, og
oppgavenummer (med delpunkt) på alle resultatlistinger og plott.

Det er fritt fram for å bruke Matlab eller Python eller andre systemer som
man måtte like bedre (Octave, etc.) så lenge man har til disposisjon de samme
numeriske og grafiske kapasitetene.

\section{Varmeledningslikning og numerisk løsning}
\setcounter{deloppg}{0}

Vi skal studere temperaturen $T(x,t)$ i en stav med lengde $L$ og
varmediffusivitet $\kappa$.  Ved tiden $t=0$ holder staven temperatur $T_0$.
For etterfølgende tider $t>0$ holdes venstre ende av staven stadig på
temperatur $T_0$, mens høyre ende av staven holdes på en høyere temperatur
$T_1>T_0$.  Problemet kan stilles opp slik:
\begin{displaymath}
  \pdt{T}{t} = \kappa \pddt{T}{{x}} \qquad \mbox{for} \quad 0 \leq x \leq L \quad \mbox{og} \quad t>0
\end{displaymath}
\begin{displaymath}
  T(x,0)=T_0 \qquad \mbox{for $0 \leq x \leq L$}
\end{displaymath}
\begin{displaymath}
  T(0,t)=T_0 \qquad \mbox{og} \qquad
  T(L,t)=T_1 \qquad \mbox{for $t>0$}
\end{displaymath}

\deloppgave \textbf{Skalering}

Ved å innføre nye skalerte og dimensjonsløse variabler
\begin{displaymath}
  u^* = \frac{T-T_0}{\alpha} \qquad t^* = \frac{t}{\beta} \qquad x^* = \frac{x}{\gamma}
\end{displaymath}
kan vi omforme problemet til normalisert form som vist nedenfor.
Stjernemerkede variabler angir \emph{dimensjonsløse måltall} (se kapittel 1.6
i Gjevik \& Fagerland).  I det følgende skal vi jobbe videre med denne
normaliserte formen:
\begin{displaymath}
  \pdt{u^*}{t^*} = \pddt{u^*}{{x^*}} \qquad \mbox{for} \quad 0 \leq x^* \leq 1 \quad \mbox{og} \quad t^*>0
\end{displaymath}
\begin{displaymath}
  u^*(x^*,0)=0 \qquad \mbox{for} \quad 0 \leq x^* \leq 1
\end{displaymath}
\begin{displaymath}
  u^*(0,t^*)=0 \qquad \mbox{og} \qquad
  u^*(1,t^*)=1 \qquad \mbox{for} \quad t^*>0
\end{displaymath}
Bestem verdiene til $\alpha$, $\beta$ og $\gamma$.

Hint: Ved hjelp av kjerneregelen for derivasjon har vi for eksempel
\begin{displaymath}
  \pd{t} = \pdt{{t^*}}{t} \pd{{t^*}} = \frac{1}{\beta} \pd{{t^*}}
\end{displaymath}

\framebox{\parbox{\textwidth}{I deloppgaver bcde nedenfor sløyfer vi
    $*$-tegnet, men vi husker at vi arbeider med dimensjonsløse størrelser,
    dermed får vi presentert likningene på pen og enkel form (se nederst side
    21 i Gjevik \& Fagerland).}}

\deloppgave \textbf{Likevektsløsning}

Med en ``likevektsløsning'' mener vi en løsning som ikke avhenger av tid,
$\pdt{u}{t}=0$, denne løsningen skal vi betegne $u = U(x)$.  Finn
likevektsløsningen for det normaliserte problemet.

\deloppgave \textbf{Den generelle løsningen nærmer seg likevektsløsningen}

Definer $\theta = u - U$.  La $\displaystyle I(t) = \int_0^1 \theta^2 \, \rmd
x$.  Vis at $I \geq 0$ og $\oder{I}{t} \leq 0$.

Hint: Gjennomfør en én-dimensjonal utgave av eksempel 8.1 i Matthews.

\bigskip

\framebox{\parbox{\textwidth}{I resten av denne oppgaven skal vi få inntrykk
    av hvordan løsningen nærmer seg likevektsløsningen, vi gjør dette ved å
    løse problemet numerisk.}}

\deloppgave \textbf{Numerisk diskretisering}

Vi skal diskretisere løsningen som $u_{m,n} = u(x_m, t_n)$ hvor $m = 1,
\ldots, M$ og $t_n = n \Delta t$ hvor $n = 0, 1, 2, \ldots$.  I Matlab lar vi
$\MatL{x=linspace(0,1,M)}$ for jevnt fordelte punkter $x_m$ i intervallet
$[0,1]$ med $\Delta x = \frac{1}{M - 1}$.

Førstederivert med hensyn på tid skal vi tilnærme med differensformelen
\begin{equation}
  \label{eq:ut}
  \left.\pdt{u}{t}\right|_{m,n} \approx \frac{u_{m,n+1} - u_{m,n}}{\Delta t}
\end{equation}
Andrederivert med hensyn på rom skal vi tilnærme med differensformelen
\begin{equation}
  \label{eq:uxx}
  \left.\pddt{u}{x}\right|_{m,n} \approx \frac{u_{m+1,n} - 2 u_{m,n} + u_{m-1,n}}{(\Delta x)^2}
\end{equation}

Vis ved Taylor-utvikling at feilen i likning (\ref{eq:ut}) er av
størrelsesorden $\Delta t$.  Finn ut av hvilken størrelsesorden feilen i
likning (\ref{eq:uxx}) er.

\deloppgave \textbf{En numerisk metode}

Ved å sette likningene (\ref{eq:ut}) og (\ref{eq:uxx}) sammen får vi følgende
\emph{eksplisitte} skjema
\begin{displaymath}
  u_{m,n+1} = r u_{m+1,n} + (1 - 2r) u_{m,n} + r u_{m-1,n}
\end{displaymath}
hvor $\displaystyle r = \frac{\Delta t}{(\Delta x)^2}$.  På semestersiden kan
du laste ned en av filene \texttt{varmeledning.m} eller
\texttt{varmeledning.py} som bruker dette skjemaet og som i tillegg sørger
for å oppfylle randkravene.

Lag et skript \texttt{Ttest.m} eller \texttt{Ttest.py} som kaller opp\\
$\MatL{u = varmeledning(ustart, 0, 1, N, r)}$ og som plotter opp
temperaturprofilen for $u(x,t)$ som funksjon av $x$ for noen utvalge tider
$t$, for eksempel $t \in \{0, 0.1, 0.2, 0.3, 0.4, 0.5\}$.  Velg startverdien
$\MatL{u = zeros(1,M)}$.  Velg for eksempel $M = 10$ og $N = 20$.  Legg merke
til at $\Delta t = 0.1/N$ for de nevnte tidene.

Vi må være forsiktig med verdien til parameteren $r$.  Vis at dersom $\Delta
t$ blir så stor i forhold til $\Delta x$ at $r > 0.5$, så kan den numeriske
løsningen bli ustabil og få urealistiske verdier.  Vis også at dersom $\Delta
t$ er så liten at $r < 0.5$, så vil den numeriske løsningen konvergere mot
likevektsløsningen.

\deloppgave \textbf{Karakteristisk tid for konvergens mot likevektsløsning}

Forhåpentligvis har forrige oppgavepunkt tjent til å overbevise om at
temperaturprofilen praktisk talt har konvergert mot likevektsprofilen etter
en dimensjonsløs tid $t^* \approx 0.5$.  Nå skal vi oversette dette til
fysisk virkelighet:

Uttrykk den tilsvarende karakteristiske fysiske tiden $t$ uttrykt ved
tykkelsen $L$ og varmediffusiviteten $\kappa$.

I tabell 10.1 i Gjevik \& Fagerland finner du varmediffusiviteter for fire
forskjellige stoffer.  Bruk disse varmediffusivitetene, anta tykkelsen er
$L=10$cm, og regn ut karakteristiske tider i sekunder.

Merk: Denne utregningen er høyst urealistisk for luft og vann fordi disse to
stoffene er flytende og vil tillate konveksjon som vi ikke har tatt hensyn
til her!

\section{Hastighetsfeltet i en strømvirvel}
\setcounter{deloppg}{0}

Vi kjenner til mange eksempler på strømvirvler: Fyll vasken eller badekaret
med vann, rør litt om, trekk ut proppen, og se på strømningen!  Dersom du
befinner deg i kano eller robåt, studer strømningen som skapes nær kanten av
åra!  Dersom løst løv ligger nært et hushjørne og det blåser litt, kan vinden
få løvet til å gå rundt i sirkelbevegelse.  I alle disse eksemplene er det
lett å innse at den sirkulære hastigheten går mot null vilkårlig langt fra
virvelens sentrum.

En enkel modell for en slik strømvirvel får vi ved å anta at strømningen er
både divergensfri $\nabla \cdot \bm{v} = 0$ og virvelfri
$\nabla \times \bm{v} = \bm{0}$.  Strømningen antas divergensfri fordi
fluidet antas inkompressibelt.  Ettersom strømningen i tillegg er
todimensjonal i det horisontale $xy$-planet, så vet vi at det finnes en
strømfunksjon $\psi$ slik at $\bm{v} = \nabla \times (- \psi \bm{k})$ hvor
$\bm{k}$ er enhetsvektor i vertikal retning.  Strømningen antas virvelfri
fordi dette generelt er en god antakelse for friksjonsfri strømning.  Med
disse antakelsene innser vi at strømfunksjonen må oppfylle Laplace likning
\begin{displaymath}
  \nabla^2 \psi = 0.
\end{displaymath}
Løsningen som representerer en strømvirvel i origo finner vi i Gjevik \&
Fagerland kapittel 9.4.4 i polare koordinater, $\psi = A \ln r$ og
$\bm{v} = v_\theta \bm{i}_\theta$ og
$v_\theta = \partial \psi / \partial r = A/r$.  Her er $x = r \cos \theta$ og
$y = r \sin \theta$, hvor $r$ er radiell avstand fra origo og $\theta$
vinkelen med $x$-aksen.  Den samme løsningen finner vi i Lindstrøm \& Hveberg
eksempel 6.5.8 i kartesiske koordinater for $A = 1$.  Løsningen er
representert grafisk ved den heltrukne kurven i figur 9.7 i Gjevik \&
Fagerland som vi har reprodusert her.  Vi ser at hastighetsfeltet går mot
null vilkårlig langt fra sentrum:

\begin{center}
  \setlength{\unitlength}{1mm}
  \begin{picture}(120,70)
    \put( 10.94,  5.25){\includegraphics[width=0.5\textwidth]{punktvirvel_graf.eps}}
    \put( 45.22, -4.15){\makebox(40,5.0)[ct]{\Large $r$}}
    \put(-31.46, 35.20){\makebox(40,5.0)[rc]{\Large $v_\theta$}}
  \end{picture}
\end{center}

Den enkle modellen ovenfor bryter sammen i origo fordi hastigheten der blir
uendelig stor og peker i alle retninger samtidig!  Naturen oppfører seg
selvfølgelig ikke slik.  Kanskje vi klarer å observere at midt i badekarets
strømvirvel er det ikke vann i det hele tatt?  Det er ikke vanskelig å
observere at løst løv i midten av en strømvirvel nær et hushjørne kun vil
rotere svært langsomt?  Vi skal gripe fatt i dette siste tilfellet: Nær
sentrum av strømvirvelen vil friksjon på grunn av viskositet være viktig,
denne friksjonen vil redusere hastigheten mot null inn mot sentrum, som
antydet ved den stiplede kurven i figuren.  Det som skjer er at antakelsen om
virvelfri bevegelse bryter sammen nær origo!

I et to-dimensjonalt divergensfritt felt i $xy$-planet kan virvlingen til
strømvektoren uttrykkes ved strømfunksjonen
\begin{displaymath}
  \bm{c}(x,y) = \nabla \times \bm {v} = \nabla^{2}\psi \bm{k}
\end{displaymath}
hvor virvlingsvektoren er $\bm{c}(x,y) = c(x,y)\bm{k}$ og $\bm{k}$ stadig er
normalvektoren til planet. Dersom en kjenner virvlingen $c(x,y)$ vil en kunne
finne strømfunksjonen ved å løse likningen
\begin{displaymath}
  \nabla^{2}\psi = c(x,y).
\end{displaymath}
Denne differensiallikningen, som er en generalisering av Laplace-likningen,
kalles Poisson-likningen.  Den dukker opp i mange forbindelser i
hydrodynamikk og elektromagnetisme.

Vi skal betrakte aksesymmetriske felt hvor virvling og strømfunksjon bare
avhenger av avstanden $r$ fra origo.  Bruker vi uttrykket for
Laplace-operatoren i polare koordinater (se kapittel 9.3 i Gjevik \&
Fagerland) finner vi
\begin{equation}
  \label{eq:poisson}
  \nabla^2 \psi =
  \frac{1}{r}\frac{\partial}{\partial r}\Big(r\frac{\partial \psi}{\partial r}\Big) = c(r).
\end{equation}

Først skal vi løse Poisson-likningen for et enkelt tilfelle: La virvlingens
størrelse være
\begin{displaymath}
  c(r) = \left\{
    \begin{array}{ll}
      \tilde{c}  & \mbox{$r \leq a$} \\
      0     & \mbox{$r > a$}
    \end{array}
  \right. 
\end{displaymath}
hvor $\tilde{c}$ er en konstant.  Husk at konstant virvling betyr rotasjon
som fast legeme, altså en ``karusell''.  For $r \leq a$ gjelder da
Poisson-likningen
\begin{displaymath}
  \nabla^2 \psi =
  \frac{1}{r}\frac{\partial}{\partial r}\Big(r\frac{\partial \psi}{\partial r}\Big) = \tilde{c},
\end{displaymath}
og for $r>a$ gjelder Laplace-likningen med høyresiden i likningen ovenfor lik
0.

\deloppgave \textbf{Analytisk løsning av Poisson- og Laplace-likningene}

Vis at den analytiske løsningen av disse likningene er
\begin{displaymath}
  \psi(r) =\left\{
    \begin{array}{ll}
      \frac{1}{4} \tilde{c}r^{2} + \alpha \ln r + \beta & \mbox{$r \leq a$} \\
      A \ln r + B     & \mbox{$r > a$  }
    \end{array}
  \right.
\end{displaymath}
hvor $\alpha$, $\beta$, $A$ og $B$ er konstanter.

Finn strømhastigheten $v_\theta = \partial \psi/\partial r$ og bestem
konstantene slik at både strømfunksjonen og strømhastigheten blir
kontinuerlige funksjoner av $r$ og slik at både strømfunksjonen og
strømhastigheten settes lik null i origo.

Sammenlikn med figur 9.7 i Gjevik \& Fagerland.  Hva har skjedd med det
singulære punktet i origo?

\bigskip

Nå skal vi lage en metode for å løse likning (\ref{eq:poisson}) numerisk.
Vi starter med å skrive
\begin{displaymath}
  \frac{\partial}{\partial r}\Big(r v \Big) = r c(r).
\end{displaymath}
Vi skal diskretisere løsningen som $v_n = v(r_n)$ og virvlingen som
$c_n = c(r_n)$ hvor $n = 1, \ldots, N$.  I Matlab lar vi
$\MatL{r=linspace(0,R,N)}$ for jevnt fordelte punkter $r_n$ i intervallet
$[0,R]$ med $\Delta r = \frac{R}{N - 1}$.

Vi skal tilnærme førstederivert med formelen
\begin{displaymath}
  \pdt{u}{r} \approx \frac{u_{n+1} - u_{n}}{\Delta r}.
\end{displaymath}
Vi tilnærmer høyresiden med middelverdien i de to punktene 
\begin{displaymath}
  r c(r) \approx \frac{r_{n+1} c_{n+1} + r_n c_n}{2}.
\end{displaymath}
Vår numeriske metode blir da
\begin{equation}
  \label{eq:metode1}
  \frac{r_{n+1} v_{n+1} - r_n v_n}{\Delta r} = \frac{r_{n+1} c_{n+1} + r_n c_n}{2}
\end{equation}
sammen med et randkrav i origo, nemlig
\begin{equation}
  \label{eq:metode2}
  v=0 \quad \mbox{for} \quad r = 0.
\end{equation}

\deloppgave \textbf{Vår numeriske metode reproduserer rotasjon som fast
  legeme eksakt}

Vis at denne numeriske metoden er eksakt for konstant virvling
$c = \tilde{c}$, det vil si for en ``karusell''.

Hint: Finn den analytiske løsningen for rotasjon som fast legeme og vis at
den oppfyller likningene (\ref{eq:metode1}) og (\ref{eq:metode2}).

\deloppgave \textbf{En numerisk løsning}

Metoden kan implementeres som en likning på matriseform
\begin{displaymath}
  M V = B
\end{displaymath}
hvor $V$ er en søylevektor med alle $v_n$, $B$ er en søylevektor med
elementer definert av høyresidene i likningene (\ref{eq:metode1}) og
(\ref{eq:metode2}), og $M$ er en matrise definert av venstresidene i
likningene (\ref{eq:metode1}) og (\ref{eq:metode2}).

Dette likningssystemet løser vi på datamaskin ved kommandoen
$V = M\backslash B$, se kapittel 4.5.2 i Lindstrøm \& Hveberg.

Lag et skript \texttt{virvel.m} eller \texttt{virvel.py} som løser likningen
numerisk.  Finn den numeriske løsningen og sammenlikn med den analytiske
løsningen funnet i deloppgave 2a.  Sett for eksempel $R = 10$, $a = 4$,
$\tilde{c} = 1$ og $N = 50$.  Framstill den analytiske og numeriske løsningen
grafisk i samme plott.

\deloppgave \textbf{En annen numerisk løsning}

Nå prøver vi en litt mer komplisert virvling, nemlig
\begin{displaymath}
  c(r) = \left\{
    \begin{array}{ll}
      1  & \mbox{$r \leq a$} \\
      \frac{r-b}{a-b}  & \mbox{$a < r \leq b$} \\
      0 & \mbox{$r > b$}
    \end{array}
  \right. 
\end{displaymath}
Sett for eksempel $R = 10$, $a = 3$, $b = 5$ og $N = 50$.

Får du et resultat som likner på den stiplede kurven i figur 9.7 i Gjevik \&
Fagerland?

\end{document}

%%% Local Variables:
%%% mode: latex
%%% TeX-master: t
%%% End:
