This demo example shows ...



In [3]:

    
%pylab inline









    



Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.zmq.pylab.backend_inline].
For more information, type 'help(pylab)'.

Vi ønsker å løse varmelikningen eksplisitt

Varmelikningen er

$\frac{\delta u}{\delta t} = \kappa \frac{\delta^2u}{\delta x^2}$

Diskretisert explisitt med hensyn på tid får vi

$\frac{1}{\Delta t}\left(u_i^{k+1} - u_i^{k}\right) = \frac{\kappa}{\Delta x^2} \left(u^k_{i-1} - 2u_i^k + u_{i+1}^k \right)$

og vår stensil

$u_i^{k+1} = su^k_{i-1} +(1-2s)u_i^k + su_{i+1}^k , \quad s=\frac{\kappa\Delta t}{\Delta x^2}$

Vi ønsker ikke å anta at vi har kartesiske grid, og vil derfor ha en stensil som beskriver fluxen over en face. Dette kan gjøres som følger:

$ \begin{eqnarray} \Delta u_i^k &=& u_i^{k+1} - u_i^k\\\ &=& su_{i-1}^k - 2su_i^k + s_{i+1}^k\\\ &=& s(u_{i-1}^k - u_i^k) - s(u_i^k - u_{i+1}^k)\\\ &=& F(u)^k_{i-1/2} - F(u)^k_{i+1/2} \end{eqnarray} $

$u_i^{k+1} = u_i^k + F(u)^k_{i-1/2} - F(u)^k_{i+1/2}$

Vi kan nå splitte beregningen opp i to deler:

Beregn flux for hver face
Summer fluxer for hver celle

Løsningen blir ustabil hvis vi ikke oppfyller CFL-kravet

$\frac{1}{2} \lt \frac{\kappa \Delta t}{\Delta x^2}$

eller med andre ord

$\Delta t \lt \frac{\Delta x^2}{2\kappa}$

Her er $u_i^k$ temperaturen i punktet $i\cdot\Delta x$ ved tid $k\cdot\Delta t$.

Vi trenger i tillegg randkrav. I vårt tilfelle benytter vi fixed randkrav, hvor

$u_0^k = r_0, u_n^k = r_1 \forall k$



In [11]:

    
import numpy as np
import time
import matplotlib.pyplot as plt
from IPython.display import clear_output

# Grid spacing
dx = 0.5

# Heat diffusion constant
kappa = 0.3

# Maximum time step size (constrained by CFL)
dt = dx*dx/(2*kappa)

# Number of grid cells
n = 10

# Boundary conditions
r0 = 0.5
r1 = 1.5

# Initial temperatures
u = np.linspace(0.0, 0.0, n)
u[0] = r0
u[n-1] = r1

def F(u_l, u_r):
    return kappa*(u_l - u_r)

fig, ax = plt.subplots()

f = np.linspace(0.0, 0.0, int(n+1)) #n+1 interfaces for n cells 
#Perform one timestep
for k in range(0, 50):
    
    # Loop over all faces and calculate the flux
    for i in range(1, n):
        f[i] = F(u[i-1], u[i]) #f[i] corresponds to the flux between cell i and i-1 (F_{i-1/2})
            
    # For each cell, loop over all faces for that cell, and sum the fluxes
    for i in range(0, n):
        u[i] = u[i] + (f[i] - f[i+1]) * dt / (dx*dx)
        
    # Apply boundary conditions
    u[0] = r0
    u[n-1] = r1
    
    # Plot the solution
    x = np.linspace(0.0, n-1, n) + 0.5
    plot(x, u, 'r.-', label='$u^{k+1}$')
    plot(f, 'b.--', label='$f^{k}$')
    pylab.ylim([r0-0.5,r1+0.5])
    time.sleep(0.2)
    clear_output()
    display(fig)
    ax.cla()
    title("Time = " + str(dt*k))
    legend()
plt.close()



In [ ]: