Eksamen MEK1100 2020

Løsningsforslag

Dette l?sningsforslaget bruker b?de programmering og direkte regning til ? l?se de fleste oppgavene. Det er strengt tatt ikke n?dvendig ? bruke programmering p? annet enn plotteoppgavene (2b, 2d, 3b) og beregning av buelengde (3d).

Oppgave 1 - Divergens og curl

Finn divergensen og virvlingen til f?lgende vektorfelt

  • i) $\vec{u} = x \,\mathbf{i} + y \,\mathbf{j} + z \,\mathbf{k}$
  • ii) $\vec{u} = r \cos \theta \,\mathbf{i_{r}} + r \sin \theta \,\mathbf{i}_{\theta} + z \,\mathbf{k}$
  • iii) $\vec{u} = \mathbf{i_r} + \mathbf{i}$

Her er $\mathbf{i}$, $\mathbf{j}$ og $\mathbf{k}$ de Kartesiske enhetsvektorene i henholdsvis $x$, $y$ og $z$-retning, mens $\mathbf{i}_r, \mathbf{i}_{\theta}, \mathbf{k}$ er enhetsvektorene i sylindriske koordinater for $r, \theta$ og $z$-retning.

i)

Divergens

\begin{align*} \nabla \cdot \vec{u} &= \frac{\partial u_x}{\partial x} + \frac{\partial u_y}{\partial y} + \frac{\partial u_z}{\partial z} \\ &= \frac{\partial x}{\partial x} + \frac{\partial y}{\partial y} + \frac{\partial z}{\partial z} \\ &= 3 \end{align*}

Implementering

In [1]:
import sympy as sp
import numpy as np
from sympy.vector import CoordSys3D
from IPython.display import Math

N = CoordSys3D('N')
x, y, z = sp.symbols('x,y,z', real=True)

ue = (x, y, z)
divu = ue[0].diff(x, 1) + ue[1].diff(y, 1) + ue[2].diff(z, 1)
divu
Out[1]:
$\displaystyle 3$
In [2]:
# Eventuelt
us = N.x*N.i + N.y*N.j + N.z*N.k
sp.vector.divergence(us)
Out[2]:
$\displaystyle 3$

Virvling

\begin{align*} \nabla \times \vec{u} &= \begin{vmatrix} \mathbf{i} & \mathbf{j} & \mathbf{k} \\ \frac{\partial }{\partial x} & \frac{\partial }{\partial y} & \frac{\partial }{\partial z} \\ x & y & z \end{vmatrix} \\ &= \mathbf{i}\left(\frac{\partial z}{\partial y} - \frac{\partial y} {\partial z} \right) + \mathbf{j}\left(\frac{\partial x}{\partial z} - \frac{\partial z} {\partial x} \right) + \mathbf{k}\left(\frac{\partial y}{\partial x} - \frac{\partial x} {\partial y} \right) \\ &= \vec{0} \end{align*}

Virvlingen er lik nullvektor.

Implementering

In [3]:
curl = (ue[2].diff(y, 1) - ue[1].diff(z, 1), ue[0].diff(z, 1) - ue[2].diff(x, 1), ue[1].diff(x, 1) - ue[0].diff(y, 1))
curl
Out[3]:
(0, 0, 0)
In [4]:
# Eventuelt
sp.vector.curl(us)
Out[4]:
$\displaystyle \mathbf{\hat{0}}$

ii)

Har n? sylinderkoordinater og

$$ \vec{u} = r \cos \theta \,\mathbf{i_{r}} + r \sin \theta \,\mathbf{i}_{\theta} + z \,\mathbf{k}. $$

Divergens

Finner divergensen

\begin{align*} \nabla \cdot \vec{u} &= \frac{1}{r}\left(\frac{\partial r u_r}{\partial r} + \frac{\partial u_{\theta}}{\partial \theta} + \frac{\partial r u_z}{\partial z} \right) \\ &= \frac{1}{r}\left(\frac{\partial r^2 \cos \theta}{\partial r} + \frac{\partial r \sin \theta}{\partial \theta} + \frac{\partial r z}{\partial z}\right) \\ &= \frac{1}{r}\left( 2r \cos \theta + r \cos \theta + r\right) \\ &= 3 \cos \theta + 1 \end{align*}

Implementering

In [5]:
r, theta, z = sp.symbols('x,y,z', real=True, positive=True)
ue = (r*sp.cos(theta), r*sp.sin(theta), z)
divu = 1/r*((ue[0]*r).diff(r, 1) + ue[1].diff(theta, 1) + (ue[2]*r).diff(z, 1))
divu = sp.simplify(divu)
Math('\\text{div}(u) = ' +sp.latex(divu, symbol_names={r: 'r', theta: '\\theta', z: 'z'}))
Out[5]:
$\displaystyle \text{div}(u) = 3 \cos{\left(\theta \right)} + 1$
In [6]:
# Eventuelt
B = N.create_new('B', transformation='cylindrical') 
us = B.r*sp.cos(B.theta)*B.i + B.r*sp.sin(B.theta)*B.j + B.z*B.k
divu = sp.vector.divergence(us)
Math('\\text{div}(u) = ' + sp.latex(divu).replace('{theta}_{B}', '\\theta'))
Out[6]:
$\displaystyle \text{div}(u) = 3 \cos{\left(\mathbf{\theta} \right)} + 1$

Virvling

\begin{align*} \nabla \times \vec{u} &= \frac{1}{r} \begin{vmatrix} \mathbf{i}_r & r\mathbf{i}_{\theta} & \mathbf{i}_z \\ \frac{\partial }{\partial r} & \frac{\partial}{\partial \theta} & \frac{\partial}{\partial z} \\ r \cos \theta & r^2 \sin \theta & z \end{vmatrix} \\ &= \frac{1}{r}\left( \mathbf{i}_r \left( \frac{\partial z}{\partial \theta} - \frac{\partial r^2 \sin \theta}{\partial z} \right) +r \mathbf{i}_{\theta} \left( \frac{\partial r \cos \theta}{\partial z} - \frac{\partial z}{\partial r}\right) +\mathbf{i}_z \left( \frac{\partial r^2 \sin \theta}{\partial r} - \frac{\partial r \cos \theta}{\partial \theta} \right) \right) \\ &= \frac{1}{r}\mathbf{i}_z\left(2 r \sin \theta + r \sin \theta \right) \\ &= 3 \sin \theta \, \mathbf{i}_z \end{align*}

Implementering

In [7]:
curl = (1/r*(ue[2].diff(theta, 1) - (r*ue[1]).diff(z, 1)),
            (ue[0].diff(z, 1) - ue[2].diff(r, 1)),
            (1/r)*((r*ue[1]).diff(r, 1) - ue[0].diff(theta, 1)))
curl = sp.simplify(curl)
Math('\\text{curl}(u) = ' + sp.latex(curl, symbol_names={theta: '\\theta'}))
Out[7]:
$\displaystyle \text{curl}(u) = \left( 0, \ 0, \ 3 \sin{\left(\theta \right)}\right)$
In [8]:
# Eventuelt
curlu = sp.vector.curl(us)
Math('\\text{curl}(u) = ' + sp.latex(curlu).replace('{theta}_{B}', '\\theta').replace('\\hat{k}_{B}', '\\mathbf{k}' ))
Out[8]:
$\displaystyle \text{curl}(u) = (3 \sin{\left(\mathbf{\theta} \right)})\mathbf{\mathbf{k}}$

iii)

$$ \vec{u} = \mathbf{i_r} + \mathbf{i} $$

Det er flere m?ter ? l?se denne oppgaven p?. Man kan velge ? transformere vektoren til Kartesiske eller sylinder-koordinater, men det er ikke n?dvendig da b?de divergens og virvling er line?re operatorer.

Ved konvertering kan man bruke

$$ \mathbf{i}_r = \cos \theta \,\mathbf{i} + \sin \theta \,\mathbf{j}, $$

og

$$ r = \sqrt{x^2+y^2}, \quad x = r \cos \theta, \quad y = r \sin \theta, $$

og skrive vektoren i Kartesiske koordinater som

$$ \vec{u} = \left(1+\frac{x}{\sqrt{x^2+y^2}} \right) \mathbf{i} + \frac{y}{\sqrt{x^2+y^2}} \,\mathbf{j}. $$

Alternativt, bruk

$$ \mathbf{i} = \cos \theta \mathbf{i}_r - \sin \theta \mathbf{i}_{\theta}, $$

og skriv vektoren som

$$ \vec{u} = \left( 1 + \cos \theta \right) \mathbf{i}_r - \sin \theta \mathbf{i}_{\theta}. $$

Alle fremgangsm?ter skal gi samme svar.

Divergens

Eventuelt

\begin{align*} \nabla \cdot (\mathbf{i}_r + \mathbf{i}) &= \nabla \cdot \mathbf{i}_r + \nabla \cdot \mathbf{i}\\ &= \nabla \cdot \mathbf{i}_r = \frac{1}{r} \frac{\partial r}{\partial r} = \frac{1}{r} \end{align*}

Eller (hopper over litt lang, men triviell utledning.)

\begin{align*} \nabla \cdot \vec{u} &= \frac{\partial }{\partial x}\left(1+\frac{x}{\sqrt{x^2+y^2}} \right) + \frac{\partial}{\partial y}\left(\frac{y}{\sqrt{x^2+y^2}}\right) \\ &= \ldots \\ &= \frac{1}{\sqrt{x^2+y^2}}\\ &= \frac{1}{r} \end{align*}

Til slutt ved rene sylinderkoordinater

\begin{align*} \nabla \cdot \vec{u} &= \frac{1}{r} \frac{\partial r(1+\cos \theta)}{\partial r} - \frac{1}{r}\frac{\partial \sin \theta}{\partial \theta} \\ &= \frac{1+\cos \theta}{r} - \frac{\cos \theta}{r} \\ &= \frac{1}{r} \end{align*}

Implementering

In [9]:
ue = (1+x/sp.sqrt(x**2+y**2), y/sp.sqrt(x**2+y**2), sp.S.Zero)
divu = ue[0].diff(x, 1) + ue[1].diff(y, 1) + ue[2].diff(z, 1)
divu = sp.simplify(divu)
divu
Out[9]:
$\displaystyle \frac{1}{\sqrt{x^{2} + y^{2}}}$
In [10]:
# Eventuelt
us = (1+N.x/sp.sqrt(N.x**2+N.y**2))*N.i + N.y/sp.sqrt(N.x**2+N.y**2)*N.j
divu = sp.vector.divergence(us)
divu = sp.simplify(divu)
divu
Out[10]:
$\displaystyle \frac{1}{\sqrt{\mathbf{{x}_{N}}^{2} + \mathbf{{y}_{N}}^{2}}}$

Her er $x_N$ og $y_N$ coordinatene i coordinatsystem N, alts? $x$ og $y$. Kunne ogs? ha implementert i sylinder-koordinater.

Virvling

Eventuelt

\begin{align*} \nabla \times (\mathbf{i}_r + \mathbf{i}) &= \nabla \times \mathbf{i}_r + \nabla \times \mathbf{i}\\ &= \nabla \times \mathbf{i}_r = \vec{0} \end{align*}

siden kun konstant koeffisient og $h_r=1$.

Eller

\begin{align*} \nabla \times \vec{u} &= \begin{vmatrix} \mathbf{i} & \mathbf{j} & \mathbf{k} \\ \frac{\partial }{\partial x} & \frac{\partial }{\partial y} & \frac{\partial }{\partial z} \\ 1+\frac{x}{r} & \frac{y}{r} & 0 \end{vmatrix} \\ &= \mathbf{i}\left(\frac{\partial 0}{\partial y} - \frac{\partial y/r} {\partial z} \right) + \mathbf{j}\left(\frac{\partial 1+x/r}{\partial z} - \frac{\partial 0} {\partial x} \right) + \mathbf{k}\left(\frac{\partial y/r}{\partial x} - \frac{\partial 1+x/r} {\partial y} \right) \\ &= \mathbf{k} \left( -\frac{xy}{r^3} + \frac{xy}{r^3} \right) \\ &= \vec{0} \end{align*}

Eller

\begin{align*} \nabla \times \vec{u} &= \frac{1}{r}\begin{vmatrix} \mathbf{i}_r & r\mathbf{i}_{\theta} & \mathbf{k} \\ \frac{\partial }{\partial r} & \frac{\partial }{\partial \theta} & \frac{\partial }{\partial z} \\ 1+\cos \theta & \sin \theta & 0 \end{vmatrix} \\ &= \frac{1}{r} \mathbf{i}_z \left(-\sin \theta - (-\sin \theta) \right)\\ &= \vec{0} \end{align*}

Implementering

In [11]:
curl = (ue[2].diff(y, 1) - ue[1].diff(z, 1), ue[0].diff(z, 1) - ue[2].diff(x, 1), ue[1].diff(x, 1) - ue[0].diff(y, 1))
curl
Out[11]:
(0, 0, 0)
In [12]:
# Eventuelt
sp.vector.curl(us)
Out[12]:
$\displaystyle \mathbf{\hat{0}}$

Oppgave 2 - Elliptiske koordinater

Elliptiske koordinater $(u, v)$ er gitt ved posisjonsvektor

$$ \vec{r} = a \cosh u \cos v \mathbf{i} + a \sinh u \sin v \mathbf{j}. $$

2a) Finn enhetsvektorer og skaleringsfaktorer.

  • Finn enhetsvektorene $\mathbf{e}_{u}$, $\mathbf{e}_{v}$ og skaleringsfaktorene til de elliptiske koordinatene. Er enhetsvektorene ortogonale?

Skaleringsfaktorene er gitt ved

$$ h_u = \left|\frac{\partial \vec{r}}{\partial u}\right|, \quad h_v = \left|\frac{\partial \vec{r}}{\partial v}\right|. $$

Enhetsvektorene er gitt ved

$$ \mathbf{e}_u = \frac{1}{h_u}\frac{\partial \vec{r}}{\partial u}, \quad \mathbf{e}_v = \frac{1}{h_v}\frac{\partial \vec{r}}{\partial v}.$$

Finner

\begin{align*} \frac{\partial \vec{r}}{\partial u} &= \frac{\partial \cosh u \cos v}{\partial u} \mathbf{i} + \frac{\partial \sinh u \sin v}{\partial u} \mathbf{j} \\ &= \sinh u \cos v \mathbf{i} + \cosh u \sin v \mathbf{j} \end{align*}

og

\begin{align*} \frac{\partial \vec{r}}{\partial v} &= \frac{\partial \cosh u \cos v}{\partial v} \mathbf{i} + \frac{\partial \sinh u \sin v}{\partial v} \mathbf{j} \\ &= -\cosh u \sin v \mathbf{i} + \sinh u \cos v \mathbf{j} \end{align*}

hvilket gir, hvis vi bruker $\cosh^2 u = 1 + \sinh^2 u$,

\begin{align*} h_u^2 &= (\sinh u \cos v \mathbf{i} + \cosh u \sin v \mathbf{j}) \cdot (\sinh u \cos v \mathbf{i} + \cosh u \sin v \mathbf{j}) \\ &= \sinh^2 u \cos^2 v + \cosh^2 u \sin^2 v \\ &= \sinh^2 u \cos^2 v + (1+\sinh^2 u) \sin^2 v \\ &= \sinh^2 u (\cos^2 v + \sin^2 v) + \sin^2 v \\ &= \sinh^2 u + \sin^2 v \end{align*}

Eventuelt kan vi bruke $\sinh^2 u = \cosh^2 u - 1$ og $\sin^2 v + \cos^2 v = 1$ og f?

\begin{align*} h_u^2 &= \sinh^2 u + \sin^2 v\\ &= \cosh^2 u - 1 + 1 - \cos^2 v \\ &= \cosh^2 u - \cos^2 v \end{align*}

Samme svar f?r vi for $h_v$

\begin{align*} h_v^2 &= (-\cosh u \sin v \mathbf{i} + \sinh u \cos v \mathbf{j}) \cdot (-\cosh u \sin v \mathbf{i} + \sinh u \cos v \mathbf{j}) \\ &= \cosh^2 u \sin^2 v + \sinh^2 u \cos^2 v \\ &= \sinh^2 u + \sin^2 v \\ &= \cosh^2 u - \cos^2 v \end{align*}

Man kan ogs? bruke andre forenklinger. For eksempel

\begin{equation*} \sinh^{2}u = \frac{1}{2}\left(\cosh 2u - 1 \right), \end{equation*}

og da f?r man

\begin{equation*} h_u = h_v = \sqrt{2 \cosh 2u - 2 \cos 2v}. \end{equation*}

S? skaleringsfaktorene er

\begin{equation*} h = h_u = h_v = \sqrt{\sinh^2 u \cos^2 v + \cosh^2 u \sin^2 v}. \end{equation*}

Med forenkling:

\begin{equation*} h = h_u = h_v = \sqrt{\sinh^2 u + \sin^2 v} = \sqrt{\cosh^2 u - \cos^2 v} = \sqrt{2 \cosh 2u - 2 \cos 2v}. \end{equation*}

Alle disse svarene er like gyldige.

Enhetsvektorene er dermed

\begin{align*} \mathbf{e}_u &= \frac{1}{h}\left(\sinh u \cos v \,\mathbf{i} + \cosh u \sin v \,\mathbf{j}\right), \\ \mathbf{e}_v &= \frac{1}{h}\left(-\cosh u \sin v \,\mathbf{i} + \sinh u \cos v \,\mathbf{j}\right). \end{align*}

Implementering

Starter med ? importere funksjonalitet fra sympy, og lager to tupler (Python immutable list) for psi=(u, v) og rv=(sinh(u)*cos(v), sinh(u)*sin(v))

In [13]:
import sympy as sp
import numpy as np
u, v = psi = sp.symbols('u,v', real=True)
rv = (sp.cosh(u)*sp.cos(v), sp.sinh(u)*sp.sin(v))

Bruker noen mer eller mindre generelle funksjoner for a finne basisvektorer, skaleringsfaktorer og enhetsvektorer. Disse er mer eller mindre eksakt som allerede g?tt igjennom p? forelesning for parabolske koordinater:

In [14]:
def basisvektorer(psi, rv):
    """Returner basisvektorer
    
    Parameters
    ----------
    psi : Tuple av nye variable
    rv : Posisjonsvektor
    """
    b = np.zeros((len(psi), len(rv)), dtype=object)
    for i, ui in enumerate(psi):
        for j, rj in enumerate(rv):
            b[i, j] = sp.simplify(rj.diff(ui, 1))
    return b

def skaleringsfaktorer(b):
    """Returner skaleringsfaktorer
    
    Parameters
    ----------
    b : basisvektorer
    """
    h = np.zeros(b.shape[0], dtype=object)
    for i, s in enumerate(np.sum(b**2, axis=1)):
        h[i] = sp.simplify(sp.sqrt(s))
        # Sympy er dessverre ikke smart nok til ? forenkle videre. Hjelper til med denne:
        h[i] = sp.simplify(h[i].replace(sp.cosh(u)**2, 1+sp.sinh(u)**2))
    return h

def enhetsvektorer(psi, rv):
    """Returner enhetsvektorer og skaleringsfaktorer
    
    Parameters
    ----------
    psi : Tuple av nye variable
    rv : Posisjonsvektor
    """
    b = basisvektorer(psi, rv)
    hi = skaleringsfaktorer(b)
    return b / hi[None, :], hi
In [15]:
e, hi = enhetsvektorer(psi, rv)

Merk at vi velger l?sningen med $h = \sqrt{\sinh^2 u + \sin^2 v}$.

In [16]:
hi = sp.sympify(hi)
hi
Out[16]:
$\displaystyle \left[\begin{matrix}\sqrt{\sin^{2}{\left(v \right)} + \sinh^{2}{\left(u \right)}} & \sqrt{\sin^{2}{\left(v \right)} + \sinh^{2}{\left(u \right)}}\end{matrix}\right]$

Skriver ut enhetsvektorene etter f?rst ? forenkle noe s? det ikke blir s? rotete:

In [17]:
h = sp.Symbol('h')
replace = lambda f: f.replace((sp.sin(v)**2+sp.sinh(u)**2)**(-1/2), 1/h)
es = replace(sp.sympify(e))
In [18]:
es
Out[18]:
$\displaystyle \left[\begin{matrix}\frac{\cos{\left(v \right)} \sinh{\left(u \right)}}{h} & \frac{\sin{\left(v \right)} \cosh{\left(u \right)}}{h}\\- \frac{\sin{\left(v \right)} \cosh{\left(u \right)}}{h} & \frac{\cos{\left(v \right)} \sinh{\left(u \right)}}{h}\end{matrix}\right]$

Printer ut sammen med de Kartesiske enhetsvektorene

In [19]:
Math('\\mathbf{e_u} = ' + sp.latex(es[0, 0]) + '\\mathbf{i} +' + sp.latex(es[0, 1]) + '\\mathbf{j} \\\\'
    +'\\mathbf{e_v} = ' + sp.latex(es[1, 0]) + '\\mathbf{i} +' + sp.latex(es[1, 1]) + '\\mathbf{j}') 
Out[19]:
$\displaystyle \mathbf{e_u} = \frac{\cos{\left(v \right)} \sinh{\left(u \right)}}{h}\mathbf{i} +\frac{\sin{\left(v \right)} \cosh{\left(u \right)}}{h}\mathbf{j} \\\mathbf{e_v} = - \frac{\sin{\left(v \right)} \cosh{\left(u \right)}}{h}\mathbf{i} +\frac{\cos{\left(v \right)} \sinh{\left(u \right)}}{h}\mathbf{j}$

Ortogonalitet

Enhetsvektorene er ortogonale siden

\begin{equation*} \mathbf{e}_u \cdot \mathbf{e}_v = 0, \end{equation*}

og

\begin{equation*} \mathbf{e}_u \cdot \mathbf{e}_u = \mathbf{e}_v \cdot \mathbf{e}_v = 1. \end{equation*}

Har i koden at $e[0]=\mathbf{e}_u$ og $e[1]=\mathbf{e}_v$. Bruker dette:

In [20]:
np.sum(e[0]*e[1])
Out[20]:
$\displaystyle 0$
In [21]:
sp.simplify(np.sum(e[0]*e[0]))
Out[21]:
$\displaystyle \frac{\sin^{2}{\left(v \right)} \cosh^{2}{\left(u \right)} + \cos^{2}{\left(v \right)} \sinh^{2}{\left(u \right)}}{\sin^{2}{\left(v \right)} + \sinh^{2}{\left(u \right)}}$

Her ser vi at sympy trenger litt hjelp til ? forenkle

In [22]:
sp.simplify((np.sum(e[0]*e[0])).replace(sp.cosh(u)**2, 1+sp.sinh(u)**2))
Out[22]:
$\displaystyle 1$

Den er med andre ord ortogonal. Samme med $\mathbf{e}_v$:

In [23]:
sp.simplify((np.sum(e[1]*e[1])).replace(sp.cosh(u)**2, 1+sp.sinh(u)**2))
Out[23]:
$\displaystyle 1$

2b) Operatorer

  • Gitt et skalarfelt $f$ og en vektor $\vec{w} = w_u \mathbf{e}_{u} + w_v \mathbf{e}_v$. Gi $\nabla f$, $\nabla \cdot \vec{w}$ og Laplace operatoren $\nabla^2$ i elliptiske koordinater.

Under bruker vi $h = \sqrt{\sinh^2 u + \sin^2 v } = \sqrt{\cosh^2 u - \cos^2 v} = \sqrt{2 \cosh 2u - 2 \cos 2v}$ i formler fra Vector Calculus. Merk at man m? forenkle for $h_u=h_v$ for ? f? full pott.

$$ \nabla f = \frac{\mathbf{e}_u}{h}\frac{\partial f}{\partial u} + \frac{\mathbf{e}_v}{h}\frac{\partial f}{\partial v}, $$$$ \nabla \cdot \vec{w} = \frac{1}{h^2}\left(\frac{\partial w_u h}{\partial u} + \frac{\partial w_v h}{\partial v} \right), $$$$ \nabla^2 = \frac{1}{h^2}\left( \frac{\partial^2}{\partial u^2} + \frac{\partial^2}{\partial v^2}\right). $$

For divergensen kan man ekspandere de deriverte, men uttrykket blir ikke enklere.

2c) Skisser koordinatkurvene

  • Presenter en skisse av koordinatkurvene i et Kartesisk koordinatsystem. Det vil si, skisser kurver med konstante verdier for b?de $u$ eller $v$, jamf?r Fig. 6.3 i Vector Calculus.
In [24]:
import matplotlib.pyplot as plt
%matplotlib inline
M = 32
# Velger 3 ganger fler punkter i v retning for ? f? jevne intervall
ui = np.broadcast_to(np.linspace(0, 1, M)[:, None], (M, 3*M))
vi = np.broadcast_to(np.linspace(0, 2*np.pi, 3*M)[None, :], (M, 3*M))

# Varierer c mellom 0 og 1 og veksler mellom ? holde u og v konstant
for c in np.linspace(0, 1, 10):
    plt.plot(np.cosh(c)*np.cos(vi[0]), np.sinh(c)*np.sin(vi[0]), 'b') # Her er c en konstant u
    plt.plot(np.cosh(ui[:, 0])*np.cos(np.pi*c), np.sinh(ui[:, 0])*np.sin(np.pi*c), 'r') # Konstant v (?verste halvdel - 0 til pi)
    plt.plot(np.cosh(ui[:, 0])*np.cos(np.pi*(c+1)), np.sinh(ui[:, 0])*np.sin(np.pi*(c+1)), 'r') # Konstant v (nederste halvdel - pi til 2pi)
plt.xlabel('x')
plt.ylabel('y')
Out[24]:
Text(0, 0.5, 'y')

2d) Kontur- og pilplott

  • Lag konturplott av skalarfeltet $$ f(u, v) = (1-u^2)\cos(2v), $$ i b?de elliptiske og Kartesiske koordinater.

    Finn $\nabla f$ (i elliptiske eller Kartesiske koordinater) og lag et pilplott av denne i samme figur som det Kartesiske konturplottet. Kommenter retningen p? pilene.

Lager skalarfelt $f(u, v) = (1-u^2) \cos(2v)$

In [25]:
f = (1-u**2)*sp.cos(2*v)

Merk at vi bruker $x=2uv$ og $y=u^2-v^2$ evaluert p? et strukturert grid. sp.lambdify er en effektiv (vektorisert) metode ? evaluere en sympy funksjon p?. S? under tilsvarer f(u, v) = sp.lambdify((u, v), f)(ui, vi).

In [26]:
fj = sp.lambdify((u, v), f)(ui, vi)

Hvis vi n? velger ? plotte $f(u, v)$ i det nye koordinatsystemet f?r vi.

In [27]:
plt.contourf(ui, vi, fj)
plt.xlabel('u')
plt.ylabel('v')
Out[27]:
Text(0, 0.5, 'v')

Men det er kanskje mer interessant ? se resultatet i fysiske (Kartesiske) koordinater. Vi trenger derfor ? finne kartesiske x, y fra de gitte u, v. Gj?r dette som f?lger

In [28]:
mesh = []
for rj in rv:
    mesh.append(sp.lambdify((u, v), rj)(ui, vi))
x, y = mesh
plt.contourf(x, y, fj)
plt.xlabel('x')
plt.ylabel('y')
Out[28]:
Text(0, 0.5, 'y')

? plotte gradienten i Kartesiske koordinater er mer involvert siden vi har beregnet gradienten i de nye koordinatene og derfor trenger ? projisere ved ? ta prikk-produktet av gradientvektoren

\begin{align*} \frac{\partial f}{\partial x} &= \nabla f \cdot \mathbf{i},\\ \frac{\partial f}{\partial y} &= \nabla f \cdot \mathbf{j}. \end{align*}

For ? finne gradientvektoren deriverer vi f?rst for ? finne komponentene til $\nabla f$ i nye koordinater

In [29]:
df = np.array((1/hi[0]*f.diff(u, 1), 1/hi[1]*f.diff(v, 1)))

Merk at df n? ikke inneholder enhetsvektorer. S? f?r vi prikker med $\mathbf{i}$ og $\mathbf{j}$ m? vi gange med enhetsvektorene $\mathbf{e_u}$ og $\mathbf{e_v}$ for ? f? $\nabla f$

$$ \nabla f = \frac{\mathbf{e}_u}{h}\frac{\partial f}{\partial u} + \frac{\mathbf{e}_v}{h}\frac{\partial f}{\partial v}, $$
In [30]:
Math('\\nabla f = ' + sp.latex(replace(df[0])) + '\mathbf{e}_u' + sp.latex(replace(df[1])) + '\mathbf{e}_v')
Out[30]:
$\displaystyle \nabla f = - \frac{2 u \cos{\left(2 v \right)}}{h}\mathbf{e}_u- \frac{2 \left(1 - u^{2}\right) \sin{\left(2 v \right)}}{h}\mathbf{e}_v$

Det kan v?re morsomt ? se p? denne ogs? i Kartesiske koordinater selv om vi bare sp?r etter $\nabla f$ i elliptiske eller Kartesiske koordinater. Da m? vi multiplisere med enhetsvektorene

In [31]:
gradf = e[0]*df[0] + e[1]*df[1]
gradf = sp.sympify(gradf)

Merk at vi med denne summen n? har f?tt satt inn for $\mathbf{e_u}$ og $\mathbf{e_v}$, s? vektoren gradf over er gitt ved Kartesiske enhetsvektorer.

In [32]:
Math('\\nabla f = ' + sp.latex(sp.simplify(gradf[0]))+'\mathbf{i} \\\\ ' + sp.latex(sp.simplify(gradf[1]))+'\mathbf{j}')
Out[32]:
$\displaystyle \nabla f = - \frac{2 u \cos{\left(v \right)} \cos{\left(2 v \right)} \sinh{\left(u \right)} + 2 \left(u^{2} - 1\right) \sin{\left(v \right)} \sin{\left(2 v \right)} \cosh{\left(u \right)}}{\sin^{2}{\left(v \right)} + \sinh^{2}{\left(u \right)}}\mathbf{i} \\ - \frac{2 u \sin{\left(v \right)} \cos{\left(2 v \right)} \cosh{\left(u \right)} - 2 \left(u^{2} - 1\right) \sin{\left(2 v \right)} \cos{\left(v \right)} \sinh{\left(u \right)}}{\sin^{2}{\left(v \right)} + \sinh^{2}{\left(u \right)}}\mathbf{j}$

Ved prikking av denne vektoren mot $\mathbf{i}$ er resultatet gradf[0], mens prikking mot $\mathbf{j}$ gir gradf[1]. Derfor hopper vi over direkte regning av prikkproduktet og henter ganske enkelt de Kartesiske vektorkomponentene vi beh?ver til plottet

In [33]:
import warnings
warnings.filterwarnings("ignore")
dfdxi = sp.lambdify((u, v), gradf[0])(ui, vi)
dfdyi = sp.lambdify((u, v), gradf[1])(ui, vi)
plt.contourf(x, y, fj)
plt.quiver(x[::2, ::2], y[::2, ::2], dfdxi[::2, ::2], dfdyi[::2, ::2], pivot='mid')
plt.xlabel('x')
plt.ylabel('y')
Out[33]:
Text(0, 0.5, 'y')

Merk at gradienten peker i retning av ?kende $f$.

Oppgave 3

En Taylor-Green virvel er en to-dimensjonell analytisk modell for periodiske str?mningsvirvler som avtar i styrke over tid. Virvlene er gitt ved

$$ \vec{u}(x,y,t) = \left(\cos x \sin y \,\mathbf{i} - \sin x \cos y \,\mathbf{j} \right) \exp^{-2\nu t}, $$

for tiden $t>=0$ i et domene $\Omega = [-\pi, \pi] \times [-\pi, \pi]$. Parameteren $\nu$ er ? regne som en konstant.

3a) Strømfunksjon og skalarpotensial

  • Vis at str?mfunksjonen er

    $$ \psi(x, y, t) = \cos x \cos y \exp^{-2 \nu t}. $$

    Hvordan kunne vi vite p? forh?nd at dette feltet har en str?mfunksjon?

Vi sjekker at divergensen til $\vec{u}$ er lik 0.

\begin{align*} \nabla \cdot \vec{u} &= \frac{\partial u_x}{\partial x} + \frac{\partial u_y}{\partial y} \\ &= \left(-\sin x \sin y - \sin x (-\sin y) \right) \exp^{-2\nu t} \\ &= 0 \end{align*}

Vi finner $\psi$ enkelt ved ? l?se disse

\begin{align*} -\frac{\partial \psi}{\partial y} &= \cos x \sin y \exp^{-2\nu t} \\ \frac{\partial \psi}{\partial x} &= -\sin x \cos y \exp^{-2\nu t} \end{align*}

Integrerer f?rste og andre likning mhp $y$ og $x$ og f?r

\begin{align*} \psi &= \cos x \cos y \exp^{-2\nu t} + f_1(x) \\ \psi &= \cos x \cos y \exp^{-2\nu t} + f_2(y) \\ \end{align*}

Ser at $f_1(x) = f_2(y) = $ konstant. Vi kan velge konstanten fritt for en str?mfunksjon og velger 0. Da har vi vist at str?mfunksjonen er som gitt i starten av oppgaven.

Implementering

In [34]:
nu, t = sp.symbols('nu t', real=True, positive=True)
us = (sp.cos(N.x)*sp.sin(N.y)*N.i - sp.sin(N.x)*sp.cos(N.y)*N.j)*sp.exp(-2*nu*t)
divu = sp.vector.divergence(us)
divu
Out[34]:
$\displaystyle 0$
  • Har dette vektorfeltet et skalarpotensial? Hvis ja, finn skalarpotensialet.

Vi sjekker om vektorfeltet er virvelfritt. Vi har et 2D vektorfelt, s? virvlingen kan kun ha én komponent, i $z$-retning

\begin{align*} \nabla \times \vec{v} &= \mathbf{k} \left(\frac{\partial u_y}{\partial x} - \frac{\partial u_x}{\partial y} \right)\\ &= \mathbf{k} \left( -\cos x \cos y - \cos x \cos y \right) \exp^{-2 \nu t} \\ &= -2 \mathbf{k} \cos x \cos y \exp^{-2 \nu t} \end{align*}

Implementering

In [35]:
curlu = sp.vector.curl(us)
curlu
Out[35]:
$\displaystyle (- 2 e^{- 2 \nu t} \cos{\left(\mathbf{{x}_{N}} \right)} \cos{\left(\mathbf{{y}_{N}} \right)})\mathbf{\hat{k}_{N}}$

Siden feltet ikke er virvelfritt har det ikke et skalarpotensial.

3b) Pilplott med strømlinjer

  • Lag et pilplott av vektorfeltet $\vec{u}(x, y, 0)$. Plott i samme figur str?mlinjene til $\vec{u}$.
In [36]:
M = 30
xi, yi = np.meshgrid(*(np.linspace(-np.pi, np.pi, M),)*2)
uss = us.subs({t: 0, nu: 1})
usx, usy = uss.dot(N.i), uss.dot(N.j)
ux = sp.lambdify((N.x, N.y), usx)(xi, yi)
uy = sp.lambdify((N.x, N.y), usy)(xi, yi)
plt.quiver(xi, yi, ux, uy)
psi = sp.cos(N.x)*sp.cos(N.y)
psij = sp.lambdify((N.x, N.y), psi)(xi, yi)
plt.contour(xi, yi, psij)
plt.xlabel('x')
plt.ylabel('y')
Out[36]:
Text(0, 0.5, 'y')

3c) Fluks og sirkulasjon

  • Hva blir fluksen

    $$ \oint_{C} \vec{u} \cdot \vec{n} ds, $$

    av vektorfeltet $\vec{u}$ ut av et rektangul?rt omr?de $\Omega = [0, \pi/2] \times [0, \pi/2]$ som omsluttes av kurven $C$? Her er $\vec{n}$ normalvektoren som peker ut fra omr?det og $ds$ er et linjeelement langs kurven $C$. Forklar hvorfor man kan bruke Gauss divergensteorem her selv om det bare er et to-dimensjonalt integral.

Merk at Gauss divergensteorem normalt gjelder en sammenheng mellom et flateintegral og et volumintegral.

Man kan vise at Gauss teorem gjelder i 2D ved ? vise til at det er en generalisering av Greens teorem, som forklart i kap 6.14.3 i Linstr?m og Hveberg "Flervariabel analyse med line?r algebra".

Man kan ogs? forestille seg at 2D kvadratet utvides til en kube ved ? ekstrahere i $z$-retningen en lengde $L$ (lengden er ikke viktig, bare at domenet n? blir 3D). Fluksintegralet over denne kuben er

$$ \oint_{A} \vec{u} \cdot \vec{n} \,d\sigma, $$

der $d\sigma$ er et flateelement og $A$ er den omsluttende flaten til kuben. Vi vet at siden $\vec{u}$ ikke er en funksjon av $z$ s? vil, for dette fluks-integralet, de to flateintegralene i $xy$-planet kansellere (flatene ved $z=0$ og $z=L$). S? vi st?r igjen med et flateintegral rundt kurven $C$, ekstrahert en dybde L.

$$ \oint_{A} \vec{u} \cdot \vec{n} d\sigma = L \oint_{C} \vec{u} \cdot \vec{n} \,ds. $$

Venstresiden kan vi bruke vanlig Gauss teorem p?, s?

$$ \int_V \nabla \cdot \vec{u} \, dV = L \oint_{C} \vec{u} \cdot \vec{n} \,ds. $$

Dermed er det klart at vi kan bruke vanlig Gauss teorem til ? finne integralet p? h?yresiden. Vi kan eventuelt forenkle videre med ? bruke $dV = dx dx dz$

\begin{align*} \int_{0}^L \int_{0}^{\pi/2} \int_0^{\pi/2} \nabla \cdot \vec{u} \, dxdydz &= \int_{0}^{L} dz \int_0^{\pi/2} \int_0^{\pi/2} \nabla \cdot \vec{u} \, dxdy, \\ &= L \int_0^{\pi/2} \int_0^{\pi/2} \nabla \cdot \vec{u} \, dxdy, \end{align*}

s? vi blir kvitt $L$'en. Men dette er ikke n?dvendig. Vi vet allerede at $\nabla \cdot \vec{u}=0$, s? vi ender opp med at

$$ \oint_{C} \vec{u} \cdot \vec{n} ds = 0. $$
  • Hva blir sirkulasjonen

    $$ \oint_{C} \vec{u} \cdot d\vec{r}, $$ n?r vi beveger oss mot klokka (alts? fra $(0, 0)$ til $(\pi/2, 0)$ osv.)? Finn resultatet b?de ved direkte regning av kurveintegralet, og ved ? benytte et passende integralteorem.

Ved direkte regning deler vi opp linjeintegralet over de fire siden av kvadratet:

\begin{align*} \oint_{C} \vec{u} \cdot d\vec{r} &= \int_{x=0}^{\pi/2} \vec{u}(x, 0) \cdot \mathbf{i} \,dx \\ &+ \int_{y=0}^{\pi/2} \vec{u}(\pi /2, y) \cdot \mathbf{j} \,dy \\ &+ \int_{x=\pi/2}^{0} \vec{u}(x, \pi/2) \cdot \mathbf{i} \,dx\\ &+ \int_{y=\pi/2}^{0} \vec{u}(0, y) \cdot \mathbf{j} \,dy. \end{align*}

Merk at vi integrerer motsatt vei enten ved ? prikke mot negativ enhetsvektor, eller ved ? snu integrasjonsgrensene, men ikke begge to. Da kanselleres effekten. Integralet er rett frem og vi kan bruke sympy:

In [37]:
(sp.integrate(us.dot(N.i).subs(N.y, 0), (N.x, 0, sp.pi/2))+
sp.integrate(us.dot(N.j).subs(N.x, sp.pi/2), (N.y, 0, sp.pi/2))+
sp.integrate(us.dot(N.i).subs(N.y, sp.pi/2), (N.x, sp.pi/2, 0))+
sp.integrate(us.dot(N.j).subs(N.x, 0), (N.y, sp.pi/2, 0)))
Out[37]:
$\displaystyle - 2 e^{- 2 \nu t}$

Ved ? benytte Stokes' teorem kan vi gjenbruke curlen:

$$ \oint_{C} \vec{u} \cdot d\vec{r} = \int_A (\nabla \times \vec{u}) \cdot \mathbf{k} \, d\sigma. $$
In [38]:
sp.integrate(sp.integrate(curlu.dot(N.k), (N.x, 0, sp.pi/2)), (N.y, 0, sp.pi/2))
Out[38]:
$\displaystyle - 2 e^{- 2 \nu t}$

Siden dette vektorfeltet er i 2D, s? kan man ogs? velge ? bruke Green's teorem, siden Stokes er en generalisering av Green's teorem til 3 dimensjoner.

3d) Ekviskalarflate og buelengde

  • La $\psi(x, y, t=0)=z$ representere h?yde og $\beta(x, y, z) = z - \cos x \cos y = 0$ en ekviskalarflate. Finn flatenormalen.
\begin{align*} \nabla \beta &= \frac{\partial \beta}{\partial x} \mathbf{i} + \frac{\partial \beta}{\partial y} \mathbf{j} + \frac{\partial \beta}{\partial z} \mathbf{k},\\ &= \sin x \cos y \mathbf{i} + \cos x \sin y \mathbf{j} + \mathbf{k}. \end{align*}

Normaliserer (see s 28-29 i Gjevik) slik at flatenormalen $\vec{n}$ er gitt ved

$$ \vec{n} = \frac{\nabla \beta}{|\nabla \beta|}, $$$$ \vec{n} = \frac{\sin x \cos y \mathbf{i} + \cos x \sin y \mathbf{j} + \mathbf{k}}{\sqrt{\sin^2 x \cos^2 y + \cos^2 x \sin^2 y + 1}}. $$

Her kan nevner forenkles videre med ? bruke $\cos^2 x = 1-\sin^2 x$, men det er ikke n?dvendig. Sympy gj?r denne 'forenklingen':

In [39]:
beta = N.z - psi
n = sp.vector.gradient(beta)
n = sp.simplify(n / n.magnitude())
n
Out[39]:
$\displaystyle (\frac{\sin{\left(\mathbf{{x}_{N}} \right)} \cos{\left(\mathbf{{y}_{N}} \right)}}{\sqrt{- \frac{\left(1 - 2 \sin^{2}{\left(\mathbf{{x}_{N}} \right)}\right) \left(1 - 2 \sin^{2}{\left(\mathbf{{y}_{N}} \right)}\right)}{2} + \frac{3}{2}}})\mathbf{\hat{i}_{N}} + (\frac{\sin{\left(\mathbf{{y}_{N}} \right)} \cos{\left(\mathbf{{x}_{N}} \right)}}{\sqrt{- \frac{\left(1 - 2 \sin^{2}{\left(\mathbf{{x}_{N}} \right)}\right) \left(1 - 2 \sin^{2}{\left(\mathbf{{y}_{N}} \right)}\right)}{2} + \frac{3}{2}}})\mathbf{\hat{j}_{N}} + (\frac{1}{\sqrt{- \frac{\left(1 - 2 \sin^{2}{\left(\mathbf{{x}_{N}} \right)}\right) \left(1 - 2 \sin^{2}{\left(\mathbf{{y}_{N}} \right)}\right)}{2} + \frac{3}{2}}})\mathbf{\hat{k}_{N}}$
  • Hvis man holder seg p? ekviskalarflaten $\beta$ og g?r en full sirkel med radius 1 ($x^2+y^2=1$) rundt origo, hvor lang er da denne buelengden?

Denne oppgaven er veldig utmattende ? l?se med penn og papir, men kan l?ses enkelt i Geogebra ved ? bruke samme fremgangsm?te som i denne appleten, som ble g?tt igjennom i forelesning. L?ses kanskje enda enklere i sympy. Har parametrisering

$$ \vec{r} = \cos(\theta)\mathbf{i} + \sin(\theta) \mathbf{j} + \cos(\cos(\theta)) \cos(\sin(\theta)) \mathbf{k}, $$

og buelengden er

$$ \int_0^{2\pi} \left|\frac{\partial \vec{r}}{\partial \theta}\right| d\theta. $$

I sympy:

In [40]:
theta = sp.Symbol('theta', real=True, positive=True)
r = sp.cos(theta)*N.i + sp.sin(theta)*N.j + sp.cos(sp.cos(theta))*sp.cos(sp.sin(theta))*N.k
dr = r.diff(theta, 1)
D = sp.Integral(dr.magnitude(), (theta, 0, 2*sp.pi)).evalf(20) # m? evaluere numerisk
print('Buelengden =', D)
Buelengden = 6.2920916462455687220

Vi beregnet buelengde i Python og Sympy i dette scriptet som ble g?tt igjennom i forelesning.

3e) Finn trykket

  • Newton's andre lov for inkompressibel str?mning (inkludert friksjon) gir oss momentumlikningen

    $$ \frac{D \vec{u}}{D t} = -\nabla p + \nu \nabla^2 \vec{u}. \label{eq:NS} $$

    Vis at den partikkelderiverte av feltet $\vec{u}$ er gitt ved

\begin{align*} \frac{D \vec{u}}{ Dt } &= \frac{\partial \vec{u}}{\partial t} + \left(\vec{u} \cdot \nabla \right) \vec{u}, \\ &= -2 \nu \vec{u} - \frac{1}{2}\left( \sin 2 x \mathbf{i} + \sin 2 y \mathbf{j} \right) \exp^{-4 \nu t}. \end{align*}
In [41]:
dudt = us.diff(t, 1)
dudt
Out[41]:
$\displaystyle (- 2 \nu e^{- 2 \nu t} \sin{\left(\mathbf{{y}_{N}} \right)} \cos{\left(\mathbf{{x}_{N}} \right)})\mathbf{\hat{i}_{N}} + (2 \nu e^{- 2 \nu t} \sin{\left(\mathbf{{x}_{N}} \right)} \cos{\left(\mathbf{{y}_{N}} \right)})\mathbf{\hat{j}_{N}}$
In [42]:
conv = sp.simplify(us.dot(sp.vector.Del())(us))
conv
Out[42]:
$\displaystyle (- \frac{e^{- 4 \nu t} \sin{\left(2 \mathbf{{x}_{N}} \right)}}{2})\mathbf{\hat{i}_{N}} + (- \frac{e^{- 4 \nu t} \sin{\left(2 \mathbf{{y}_{N}} \right)}}{2})\mathbf{\hat{j}_{N}}$

Ser at dudt $=-2\nu \vec{u}$ og conv $=- 1/2 \left( \sin 2 x \mathbf{i} + \sin 2 y \mathbf{j} \right) \exp^{-4 \nu t}$.

Bruk dette og videre innsetting i momentumlikningen til ? finne trykket $p$.

Finner f?rst Laplace uttrykket:

In [43]:
nu*sp.vector.Laplacian(us).doit()
Out[43]:
$\displaystyle (- 2 \nu e^{- 2 \nu t} \sin{\left(\mathbf{{y}_{N}} \right)} \cos{\left(\mathbf{{x}_{N}} \right)})\mathbf{\hat{i}_{N}} + (2 \nu e^{- 2 \nu t} \sin{\left(\mathbf{{x}_{N}} \right)} \cos{\left(\mathbf{{y}_{N}} \right)})\mathbf{\hat{j}_{N}}$

Ser at Laplace og tidsderivert kansellerer eksakt:

In [44]:
sp.simplify(nu*sp.vector.Laplacian(us).doit() - dudt)
Out[44]:
$\displaystyle \mathbf{\hat{0}}$

St?r igjen med

$$ (\vec{u} \cdot \nabla) \vec{u} = -\nabla p. $$

Har da

\begin{align*} \frac{\partial p}{\partial x} &= -(\vec{u} \cdot \nabla) \vec{u} \cdot \mathbf{i}, \\ \frac{\partial p}{\partial y} &= -(\vec{u} \cdot \nabla) \vec{u} \cdot \mathbf{j}. \end{align*}

Eller

\begin{align*} \frac{\partial p}{\partial x} &= -\frac{1}{2}\sin 2x \exp^{-4 \nu t}, \\ \frac{\partial p}{\partial y} &= -\frac{1}{2}\sin 2y \exp^{-4 \nu t}. \end{align*}

Som kan l?ses til

$$ p(x, y, t) = -\frac{1}{4}(\cos 2x + \cos 2y) \exp^{-4 \nu t} + p_0, $$

der $p_0$ er en vilk?rlig konstant.

In [ ]: