Contact
CoCalc Logo Icon
StoreFeaturesDocsShareSupport News Sign UpSign In
| Download
Project: WISB251
Views: 88
Kernel: Python 3 (Ubuntu Linux)

Stromingssimulatie

We kunnen de verspreiding van een concentratie u(t,x)u(t,x) als gevolg van een stroming vv beschrijven als

ut+vux=0,\frac{\partial u}{\partial t} + v\frac{\partial u}{\partial x} = 0,

met x[0,1]x\in [0,1], t[0,T]t\in [0,T] en initiële conditie u(0,x)=f(x)u(0,x) = f(x). We zien meteen dat u(t,x)=f(xvt)u(t,x) = f(x - vt) een oplossing is, en kunnen dus gemakkelijk numerieke methoden uittesten.

Om deze vergelijking op te lossen moeten we zowel in de tijd als in de ruimte discretiseren. Voor de ruimtelijke afgeleiden hebben we 3 keuzes

  • voorwaartse differentie: ux(u(t,x+Δx)u(t,x))/Δx\frac{\partial u}{\partial x} \approx (u(t,x + \Delta x) - u(t,x))/\Delta x,

  • terugwaartse differentie: ux(u(t,x)u(t,xΔx))/Δx\frac{\partial u}{\partial x} \approx (u(t,x) - u(t,x - \Delta x))/\Delta x,

  • centrale differentie: ux(u(t,x+Δx)u(t,xΔx))/(2Δx)\frac{\partial u}{\partial x} \approx (u(t,x+\Delta x) - u(t,x - \Delta x))/(2\Delta x),

Als we nu de functie u(t,x)u(t,x) op nn rooster punten 0,Δx,2Δx,,(n1)Δx0, \Delta x, 2\Delta x, \ldots, (n-1)\Delta x bekijken kunnen kunnen we de voorwaartse differentie bijvoorbeeld schrijven als (uk+1(t)uk(t))/Δx(u_{k+1}(t) - u_{k}(t))/\Delta x, met uk(t)u(kΔx,t)u_k(t) \equiv u(k\Delta x,t). We krijgen nu een stelsel gewonen differentiaalvergelijkingen

u˙(t)+vDu=0,\dot{\mathbf{u}}(t) + v D\mathbf{u} = 0,

waarbij de matrix DD de eerste afgeleide representeert en u(t)=[u0(t),u1(t),,un1(t)]\mathbf{u}(t) = [u_0(t),u_1(t), \ldots, u_{n-1}(t)].

We zullen in eerste instantie EF gebruiken om dit stelsel op te lossen. We gaan kijken welk van deze benaderingen van de afgeleide naar xx het meest geschikt is om de advectievergelijking op te lossen. Hierbij moeten we rekening houden dat vv zowel positief als negatief kan zijn!

Een uitbreiding is het invoeren van een variabele stroomsnelheid v(x)v(x), die dus tegelijk positief als negatief kan zijn. Hoe zouden we hier mee om kunnen gaan?

Een verdere uitbreiding is het bekijken van tweedimensionale advectie

ut+vxux+vyuy=0.\frac{\partial u}{\partial t} + v_x\frac{\partial u}{\partial x} + v_y\frac{\partial u}{\partial y}= 0.

Als vxv_x en vyv_y ruimtelijk varieren moeten ze voldoen aan

vxx+vyy=0.\frac{\partial v_x}{\partial x} + \frac{\partial v_y}{\partial y} = 0.

Vragen voor het verslag

Basis

  • Hoe zien de verschillende matrices DD eruit? Hoe zit het met de nauwkeurigheid van de benadering?

  • Test en analyseer de verschillende ruitemlijke discretisaties in combinatie met EF voor het 1D geval. Kijk naar de orde van de fout in Δt\Delta t en Δx\Delta x en de stabiliteit als functie van Δt\Delta t en Δx\Delta x.

  • Test de verschillende methoden ook numeriek. Gebruik bijvoorbeeld f(x)=exp(100(x.5)2)f(x) = \exp(-100(x-.5)^2).

  • Hoe zit het met behoud van de vorm van initiële conditie? Kun je dit gedrag verklaren?

Extra

  • Als alternatieve methode is er de Lax methode: un+1=AunvΔtDun\mathbf{u}^{n+1} = A\mathbf{u}^{n} - v\Delta t D\mathbf{u}^{n}, waar AA gemiddelde berekend van uu op ieder punt, (u(x+Δx)+u(xΔx))/2(u(x+\Delta x) + u(x-\Delta x))/2 en DD de centrale differentie (u(x+Δx)u(xΔx))/(2Δx)(u(x+\Delta x) -u(x-\Delta x))/(2\Delta x) berekend. Hoe zit het met de stabiliteit en nauwkeurigheid van deze methode?

  • Test je methode ook eens in 2D met vx(x,y)=sin(2πx)sin(2πy)v_x(x,y) = \sin(2\pi x)\sin(2\pi y) en vy(x,y)=cos(2πx)cos(2πy)v_y(x,y) = \cos(2\pi x)\cos(2\pi y).

Code

Om je opweg te helpen staat hieronder alvast wat code

import numpy as np import scipy.sparse as sp import matplotlib.pyplot as plt from scipy.integrate import odeint n = 10 h = 1/(n-1) # voorwaarste differenties Df = (1/h)*sp.diags([-1,1],[0,1],(n,n)) # terugwaarste differenties Db = (1/h)*sp.diags([-1,1],[-1,0],(n,n)) # bekijk hoe de matrices eruit zien plt.spy(Df)
<matplotlib.lines.Line2D at 0x7fca3f27df98>
Image in a Jupyter notebook

We kunnen nu de advectievergelijking oplossen met EF:

v = -1 T = 1 n = 100 h = 1/(n-1) dt = 1e-3 nt = int(T/dt) + 1 # rooster x = np.linspace(0,1,n) t = np.linspace(0,T,int(T/dt)+1) # differentiatiematrix Df = (1/h)*sp.diags([-1,1],[0,1],(n,n)) # beginvoorwaarde u0 = np.exp(-1e3*(x-.5)**2) # los op met EF u = np.zeros((len(t),n)) u[0] = u0 for k in range(len(t)-1): u[k+1] = u[k] - dt*v*Df@u[k] # plot plt.plot(x,u[0,:],x,u[50,:],x,u[100,:]) plt.xlabel('x') plt.legend(('t = 0', 't = 0.05', 't = 0.1'))
<matplotlib.legend.Legend at 0x7fca3efc3860>
Image in a Jupyter notebook

En in 2D kunnen we de methode als volgt implementeren

T = 1 n = 100 h = 1/(n-1) dt = 1e-3 nt = int(T/dt) + 1 # rooster x = np.linspace(0,1,n) t = np.linspace(0,T,int(T/dt)+1) # differentiatiematrix Df = (1/h)*sp.diags([-1,1],[0,1],(n,n)) # vectorveld xx,yy = np.meshgrid(np.linspace(0,1,n),np.linspace(0,1,n)) vx = np.sin(2*np.pi*yy)*np.sin(2*np.pi*xx) vy = np.cos(2*np.pi*yy)*np.cos(2*np.pi*xx) # beginconditie w0 = np.random.randn(n,n); # los op met EF w = np.zeros((n,n,nt)) w[:,:,0] = w0 for k in range(len(t)-1): w[:,:,k+1] = w[:,:,k] - dt*vx*(Df@w[:,:,k]) - dt*vy*(w[:,:,k]@Df.T)