Forelesning 12. april

Fredrik Eikeland Fossan

Apr 12, 2016

Oppstart av rørstrømning

Dersom vi transformerer diffusjonsligningen $ \frac{\partial u}{\partial t}=\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}+\frac{\partial^2 u}{\partial z^2} $ til sylinder-koordinater og forlanger at $ u $ bare skal være funksjon av tiden $ t $ og radien $ r $, får vi: $$ \begin{equation*} \frac{\partial u}{\partial t}= \frac{\partial^2 u}{\partial r^2}+\frac{1}{r}\frac{\partial u }{\partial r} \end{equation*} $$

Figure 1:

Figuren viser hastighetsprofilet ved strømning i et rør av en inkompressibel fluid ved et gitt tidspunkt. Vi tenker oss at profilet har utviklet seg fra null ved å sette på en konstant trykkgradient $ \frac{dp}{dz} < 0 $, slik at hastighetsprofilet for det stasjonære tilfellet er det velkjente parabolske profilet for Poiseuille-strømning.

Bevegelsesligning: $$ \begin{equation} \label{eq:56010} \frac{\partial U}{\partial \tau}= -\frac{1}{\rho}\frac{dp}{dz}+\nu \left( \frac{\partial^2 U}{\partial R^2}+\frac{1}{R}\frac{\partial U}{\partial R} \right) \end{equation} $$

der hastighetsprofilet $ U=U(R,\ \tau).\ 0\leq R \leq R_0, $ og $ \tau $ er den fysikalske tiden.

Kan vise at problemet kan skrives om til: $$ \begin{equation} \label{eq:56015a} \frac{\partial \omega}{\partial t}= \frac{\partial^2\omega}{\partial r^2}+ \frac{1}{r} \frac{\partial \omega}{\partial r} \end{equation} $$

Randbetingelser: $$ \begin{equation} \label{eq:56015b} \omega(\pm1,t)=0,\ \frac{\partial \omega}{\partial r}(0,t)=0 \end{equation} $$

Startbetingelse: $$ \begin{equation} \label{eq:56015c} \omega(r,0)=u_s=1-r^2 \end{equation} $$

Det opprinnelige problemet er da: $$ \begin{equation} \label{eq:56016} u(r,t)=1-r^2-\omega(r,t) \end{equation} $$

$ \theta $-skjemaet for diffusjonslikningen i sylindriske koordinater kan skrives: $$ \begin{equation} \label{eq:theta_syl} \begin{array}{l} -D\theta\cdot\left(1-\frac{1}{2j}\right)\cdot \omega^{n+1}_{j-1}+(1+2D\theta)\cdot \omega_j^{n+1}-D\theta \cdot\left(1+\frac{1}{2j}\right)\cdot \omega^{n+1}_{j+1}\\ = D(1-\theta)\cdot \left(1-\frac{1}{2j}\right)\cdot \omega^n_{j-1}+[1-2D(1-\theta)]\cdot\omega^n_j +D(1-\theta)\cdot\left(1+\frac{1}{2j}\right)\cdot\omega^n_{j+1} \end{array} \end{equation} $$

men ser at vi får ett problem ved $ r=0 $; Leddet $ \dfrac{1}{r}\dfrac{\partial u}{\partial r} $ må behandles spesielt for $ r=0 $.

Problem 1: L'hopital for r = 0

Versjon 1

L'Hopitals regel: $$ \begin{equation} \label{eq:5606} \lim_{r\to 0}\frac{1}{r}\frac{\partial u}{\partial r} = 1 \frac{\partial^2u}{\partial r^2}\to \frac{\partial u}{\partial t}= 2\frac{\partial^2u}{\partial r^2} \text{ for }r=0 \end{equation} $$

Diskretiserer \eqref{eq:5606} for $ r=0 $ og utnytter symmetribetingelsen $ \frac{\partial u}{\partial r}(0)=0 $: $$ \begin{equation} \label{eq:5607} (1+4D\theta)\omega^{n+1}_0 -4D\theta \omega_1^{n+1}=\omega_0^n + 4D(1-\theta)(\omega_1^n - \omega_0^n) \end{equation} $$

Hint.

$$ \begin{equation} \left(\begin{array}{ccccccccc} 1 + 4D\theta & -4D\theta & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ -D\theta(1-\frac{0.5}{1})& 1 + 2D \theta & -D\theta(1+\frac{0.5}{1}) & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0& -D\theta(1-\frac{0.5}{2}) & 1 + 2D \theta & -D\theta(1+\frac{0.5}{2}) & 0 & 0 & 0 & 0 & 0 & 0 \\ 0& 0 & -D\theta(1-\frac{0.5}{3}) & 1 + 2D \theta & -D\theta(1+\frac{0.5}{3}) \ & 0 & 0 & 0 & 0 & 0 \\ 0& 0 & 0 & -D\theta(1-\frac{0.5}{4}) & 1 + 2D \theta & -D\theta(1+\frac{0.5}{4}) & 0 & 0 & 0 & 0 \\ 0& 0 & 0 & 0 & -D\theta(1-\frac{0.5}{5}) & 1 + 2D \theta & -D\theta(1+\frac{0.5}{5}) & 0 & 0 & 0 \\ 0& 0 & 0 & 0 & 0 & -D\theta(1-\frac{0.5}{6}) & 1 + 2D \theta & -D\theta(1+\frac{0.5}{6}) \theta & 0 & 0 \\ 0& 0 & 0 & 0 & 0 & 0 & -D\theta(1-\frac{0.5}{7}) & 1 + 2D \theta & -D\theta(1+\frac{0.5}{7}) \theta & 0 \\ 0& 0 & 0 & 0 & 0 & 0 & 0 & -D\theta(1-\frac{0.5}{8}) & 1 + 2D \theta & -D\theta(1+\frac{0.5}{8}) \\ 0& 0 & 0 & 0 & 0 & 0 & 0 & 0 & -D\theta(1-\frac{0.5}{9}) & 1 + 2D \theta & \end{array}\right) \cdot \left(\begin{array}{c} \omega_{0}^{n+1}\\ \omega_{1}^{n+1}\\ \omega_{2}^{n+1}\\ \omega_{3}^{n+1}\\ \omega_{4}^{n+1}\\ \omega_{5}^{n+1}\\ \omega_{6}^{n+1}\\ \omega_{7}^{n+1}\\ \omega_{8}^{n+1}\\ \omega_{9}^{n+1} \end{array}\right) \\ = \left(\begin{array}{c} \omega_0^n + 4D(1-\theta)(\omega_1^n - \omega_0^n)\\ D(1-\theta)(1-\frac{0.5}{1})\omega^n_{0} + \left(1 - 2D\left(1 - \theta \right)\right)\omega^n_{1} + D(1-\theta)(1+\frac{0.5}{1})\omega^n_{2}\\ D(1-\theta)(1-\frac{0.5}{2})\omega^n_{1} + \left(1 - 2D\left(1 - \theta \right)\right)u^n_{2} + D(1-\theta)(1+\frac{0.5}{2})\omega^n_{3} \\ D(1-\theta)(1-\frac{0.5}{3})\omega^n_{2} + \left(1 - 2D\left(1 - \theta \right)\right)u^n_{3} + D(1-\theta)(1+\frac{0.5}{3})\omega^n_{4}\\ D(1-\theta)(1-\frac{0.5}{4})\omega^n_{3} + \left(1 - 2D\left(1 - \theta \right)\right)u^n_{4} + D(1-\theta)(1+\frac{0.5}{4})\omega^n_{5}\\ D(1-\theta)(1-\frac{0.5}{5})\omega^n_{4} + \left(1 - 2D\left(1 - \theta \right)\right)u^n_{5} + D(1-\theta)(1+\frac{0.5}{5})\omega^n_{6}\\ D(1-\theta)(1-\frac{0.5}{6})\omega^n_{5} + \left(1 - 2D\left(1 - \theta \right)\right)u^n_{6} + D(1-\theta)(1+\frac{0.5}{6})\omega^n_{7}\\ D(1-\theta)(1-\frac{0.5}{7})\omega^n_{6} + \left(1 - 2D\left(1 - \theta \right)\right)u^n_{7} + D(1-\theta)(1+\frac{0.5}{7})\omega^n_{8}\\ D(1-\theta)(1-\frac{0.5}{8})\omega^n_{7} + \left(1 - 2D\left(1 - \theta \right)\right)u^n_{8} + D(1-\theta)(1+\frac{0.5}{8})\omega^n_{9}\\ D(1-\theta)(1-\frac{0.5}{9})\omega^n_{8} + \left(1 - 2D\left(1 - \theta \right)\right)u^n_{9} + D(1-\theta)(1+\frac{0.5}{9})\omega^n_{10} \end{array}\right) \label{eq:theta_matrix_cyl} \end{equation} $$

Problem 2: 2. ordens foroverdifferanse for r = 0

Versjon 2

kan bruke 2. ordens foroverdifferanser for $ j=0 $ og grensebetingelsen gitt i Eq. \eqref{eq:56015b}: $$ \begin{equation} \label{eq:5609} \frac{\partial u }{\partial r}(0)=0\to \frac{-3u_0^n+4u^n_1-u_2^n}{2\Delta r}\to u^n_0=\frac{1}{3}(4u_1^n-u_2^n) \end{equation} $$

som igjen gir: $$ \begin{equation}\label{eq:56024} \omega^n_0=\frac{1}{3}(4\omega^n_1-\omega_2^n),\ n=0,1,2\dots \end{equation} $$ $$ \begin{align} &-D\theta\cdot\left(1-\frac{1}{2j}\right)\cdot \omega^{n+1}_{j-1}+(1+2D\theta)\cdot \omega_j^{n+1}-D\theta \cdot\left(1+\frac{1}{2j}\right)\cdot \omega^{n+1}_{j+1} \label{_auto1}\\ &= D(1-\theta)\cdot \left(1-\frac{1}{2j}\right)\cdot \omega^n_{j-1}+[1-2D(1-\theta)]\cdot\omega^n_j \label{eq:56025} \\ &+D(1-\theta)\cdot\left(1+\frac{1}{2j}\right)\cdot\omega^n_{j+1} \label{_auto2} \end{align} $$

\eqref{eq:56025} utskrevet for $ j=1 $, innsatt fra \eqref{eq:56024} og sammentrukket: $$ \begin{equation} \label{eq:56026} (1+\frac{4}{3}D\theta)\cdot \omega^{n+1}_1-\frac{4}{3}D\theta \omega_2^{n+1}=[1-\frac{4}{3}D(1-\theta)]\cdot \omega^n_1+\frac{4}{3}D(1-\theta)\omega_2^n \end{equation} $$

Hint 1.

$$ \begin{equation} \left(\begin{array}{ccccccccc} 1 + \frac{4}{3}D\theta & -\frac{4}{3}D\theta & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ -D\theta(1-\frac{0.5}{2}) & 1 + 2D \theta & -D\theta(1+\frac{0.5}{2}) & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & -D\theta(1-\frac{0.5}{3}) & 1 + 2D \theta & -D\theta(1+\frac{0.5}{3}) \ & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & -D\theta(1-\frac{0.5}{4}) & 1 + 2D \theta & -D\theta(1+\frac{0.5}{4}) & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & -D\theta(1-\frac{0.5}{5}) & 1 + 2D \theta & -D\theta(1+\frac{0.5}{5}) & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & -D\theta(1-\frac{0.5}{6}) & 1 + 2D \theta & -D\theta(1+\frac{0.5}{6}) \theta & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & -D\theta(1-\frac{0.5}{7}) & 1 + 2D \theta & -D\theta(1+\frac{0.5}{7}) \theta & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & -D\theta(1-\frac{0.5}{8}) & 1 + 2D \theta & -D\theta(1+\frac{0.5}{8}) \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & -D\theta(1-\frac{0.5}{9}) & 1 + 2D \theta & \end{array}\right) \cdot \left(\begin{array}{c} \omega_{1}^{n+1}\\ \omega_{2}^{n+1}\\ \omega_{3}^{n+1}\\ \omega_{4}^{n+1}\\ \omega_{5}^{n+1}\\ \omega_{6}^{n+1}\\ \omega_{7}^{n+1}\\ \omega_{8}^{n+1}\\ \omega_{9}^{n+1} \end{array}\right) \\ = \left(\begin{array}{c} (1 - \frac{4}{3}D(1-\theta))\omega^n_{1} + \frac{4}{3}D(1-\theta)\omega^n_{2} \\ D(1-\theta)(1-\frac{0.5}{2})\omega^n_{1} + \left(1 - 2D\left(1 - \theta \right)\right)\omega^n_{2} + D(1-\theta)(1+\frac{0.5}{2})\omega^n_{3} \\ D(1-\theta)(1-\frac{0.5}{3})\omega^n_{2} + \left(1 - 2D\left(1 - \theta \right)\right)\omega^n_{3} + D(1-\theta)(1+\frac{0.5}{3})\omega^n_{4}\\ D(1-\theta)(1-\frac{0.5}{4})\omega^n_{3} + \left(1 - 2D\left(1 - \theta \right)\right)\omega^n_{4} + D(1-\theta)(1+\frac{0.5}{4})\omega^n_{5}\\ D(1-\theta)(1-\frac{0.5}{5})\omega^n_{4} + \left(1 - 2D\left(1 - \theta \right)\right)\omega^n_{5} + D(1-\theta)(1+\frac{0.5}{5})\omega^n_{6}\\ D(1-\theta)(1-\frac{0.5}{6})\omega^n_{5} + \left(1 - 2D\left(1 - \theta \right)\right)\omega^n_{6} + D(1-\theta)(1+\frac{0.5}{6})\omega^n_{7}\\ D(1-\theta)(1-\frac{0.5}{7})\omega^n_{6} + \left(1 - 2D\left(1 - \theta \right)\right)\omega^n_{7} + D(1-\theta)(1+\frac{0.5}{7})\omega^n_{8}\\ D(1-\theta)(1-\frac{0.5}{8})\omega^n_{7} + \left(1 - 2D\left(1 - \theta \right)\right)\omega^n_{8} + D(1-\theta)(1+\frac{0.5}{8})\omega^n_{9}\\ D(1-\theta)(1-\frac{0.5}{9})\omega^n_{8} + \left(1 - 2D\left(1 - \theta \right)\right)\omega^n_{9} + D(1-\theta)(1+\frac{0.5}{9})\omega^n_{10} \end{array}\right) \label{eq:theta_matrix_cyl_V2} \end{equation} $$

Hint 2.

# lec-12-apr/startup.py;Visualization.py @ git@lrhgit/tkt4140/src/src-ch5/Visualization.py;Startupfuncs.py @ git@lrhgit/tkt4140/src/src-ch5/Startupfuncs.py;

import scipy 
import scipy.linalg
import scipy.sparse
import scipy.sparse.linalg

import numpy as np


def thetaScheme_numpy(theta, D, N, wOld):
    """ solve problem by array-slicing and scipy sparse solver
        and using L'hopital at r=0
    """
    superDiag = np.zeros(N - 1)
    subDiag = np.zeros(N - 1)
    mainDiag = np.zeros(N)
    
    RHS = np.zeros(N)
    
    j_array = np.linspace(0, N, N + 1)
    tmp = D*(1. - theta)
    
    superDiag[1:] = -D*theta*(1 + 0.5/j_array[1:-2])
    mainDiag[1:] = np.ones(N - 1)*(1 + 2*D*theta)
    subDiag[:] = -D*theta*(1 - 0.5/j_array[1:-1])
    
    a = tmp*(1 - 1./(2*j_array[1:-1]))*wOld[0:-2]
    b = (1 - 2*tmp)*wOld[1:-1]
    c = tmp*(1 + 1/(2*j_array[1:-1]))*wOld[2:]
    
    RHS[1:] = a + b + c
    
    superDiag[0] = -4*D*theta
    mainDiag[0] = 1 + 4*D*theta
    RHS[0] = (1 - 4*tmp)*wOld[0] + 4*tmp*wOld[1] 
    
    A = scipy.sparse.diags([subDiag, mainDiag, superDiag], [-1, 0, 1], format='csc')
    
    wNew = scipy.sparse.linalg.spsolve(A, RHS)
    wNew = np.append(wNew, 0)

    
    return wNew

def thetaScheme_numpy_V2(theta, D, N, wOld):
    """ solve problem by array-slicing and scipy sparse solver
        and using 2nd order forward difference at r=0
    """
    superDiag = np.zeros(N - ?)
    subDiag = np.zeros(N - ?)
    mainDiag = np.zeros(N - ?)
    
    RHS = np.zeros(N - ?)
    
    j_array = np.linspace(0, N, N + 1)
    tmp = D*(1. - theta)
    
    superDiag[1:] = -D*theta*(1 + 0.5/j_array[?:-?])
    mainDiag[1:] = np.ones(N - ?)*(1 + 2*D*theta)
    subDiag[:] = -D*theta*(1 - 0.5/j_array[?:-?])
    
    a = tmp*(1 - 1./(2*j_array[?:-?]))*wOld[?:-?]
    b = (1 - 2*tmp)*wOld[?:-?]
    c = tmp*(1 + 1/(2*j_array[?:-?]))*wOld[?:?]
    
    RHS[1:] = a + b + c
    
    superDiag[0] = -(4./3)*D*theta
    mainDiag[0] = 1 + (4./3)*D*theta
    RHS[0] = (1 - (4./3)*tmp)*wOld[?] + (4./3)*tmp*wOld[?] 
    
    A = scipy.sparse.diags([subDiag, mainDiag, superDiag], [-1, 0, 1], format='csc')
    
    wNew = scipy.sparse.linalg.spsolve(A, RHS)
    w_0 = (1./3)*(4*wNew[?] - wNew[?])
    
    wNew = np.append(w_0, wNew)
    wNew = np.append(wNew, 0)

    
    return wNew


def thetaScheme_numpy_V3(theta, D, N, wOld):
    """ solve problem by array-slicing and scipy sparse solver
        and using scheme also for r=0. 1/r cancels
    """
    superDiag = np.zeros(N - 1)
    subDiag = np.zeros(N - 1)
    mainDiag = np.zeros(N)
    
    RHS = np.zeros(N)
    
    j_array = np.linspace(0, N, N + 1)
    tmp = D*(1. - theta)
    
    superDiag[1:] = -D*theta*(1 + 0.5/j_array[1:-2])
    mainDiag[1:] = np.ones(N - 1)*(1 + 2*D*theta)
    subDiag[:] = -D*theta*(1 - 0.5/j_array[1:-1])
    
    a = tmp*(1 - 1./(2*j_array[1:-1]))*wOld[0:-2]
    b = (1 - 2*tmp)*wOld[1:-1]
    c = tmp*(1 + 1/(2*j_array[1:-1]))*wOld[2:]
    
    RHS[1:] = a + b + c
    
    superDiag[0] = -2*D*theta
    mainDiag[0] = 1 + 2*D*theta
    RHS[0] = (1 - 2*tmp)*wOld[0] + 2*tmp*wOld[1] 
    
    A = scipy.sparse.diags([subDiag, mainDiag, superDiag], [-1, 0, 1], format='csc')
    
    wNew = scipy.sparse.linalg.spsolve(A, RHS)
    wNew = np.append(wNew, 0)

    
    return wNew






if __name__ == '__main__':
    
    import matplotlib.pylab as plt
    import time as timeModule
    from Visualization import createAnimation
    from Startupfuncs import analyticSolution

    # change some default values to make plots more readable 
    LNWDT=3; FNT=9
    plt.rcParams['lines.linewidth'] = LNWDT; plt.rcParams['font.size'] = FNT
        
    #### main program starts here ####
    
    N = 25 # No. of parts
    h = 1./N # length of r-step
    D = 0.5 # numerical diffusion number
    dt = D*h**2
    nmax = 1500 # No. of time-steps
    T = nmax*dt
    time = np.linspace(0, T, nmax + 1)
    r = np.linspace(0, 1, N + 1)
    eps = 2.2204*10**-16 #tollerance to be used in calculations of analytical solution

    theta = 0
    scheme_name = 'FTCS'
    
    wStart = 1 - r**2 # initial profile is the pouseille profile 
    wOld = wStart
    #initialize solution matrices
    uSolutions = np.zeros((1, len(time), N + 1))
    uAnalytic = np.zeros((len(time), N + 1))
    
    # solve numerically with vectorized solver

    tic = timeModule.time()
    for i in range(nmax):
        
        wNew = thetaScheme_numpy_V2(theta, D, N, wOld)
        uNew = 1 - r**2 - wNew
        
        
        uSolutions[0, i + 1, :] = uNew
        wOld = wNew
    
    toc = timeModule.time()
    print "vectorized scipy sparse solver CPU-time: ", toc - tic
    

    animate = True
    if animate:
        # solve analytically:
        tic = timeModule.time()
        for i in range(nmax):
             
            t = time[i + 1]
            for k, r_variable in enumerate(r):
                uAnalytic[i + 1, k] = 1 - r_variable**2 - analyticSolution(r_variable, t, eps)
     
        toc = timeModule.time()
        print "analytic solution CPU-time: ", toc - tic

        createAnimation(uSolutions, uAnalytic, [scheme_name], r, time, jump = int(round(0.002/dt)))

Avkjøling av kule

$$ \begin{equation} \label{eq:5601} \frac{\partial u}{\partial t}= \frac{\partial^2 u}{\partial r^2}+\frac{\lambda}{r}\frac{\partial u }{\partial r},\ \lambda=0,1,2 \end{equation} $$

Example 3: Avkjøling av kule

Figuren viser en kule som blir avkjølt i vann. Kula har radius $ b=5 $ og har temperatur $ T_k $ før den senkes i vannet. Vannet holder en konstant temperatur $ T_v $ under hele prosessen. Vi ser bort fra varmetap til omgivelsene.

Andre data: $$ \begin{align} &Varmeledningstall: k=0.1W/(cm\cdot^\circ C) \label{eq:56029a} \\ &Varmeovergangstall: \bar h =0.2W/(cm\cdot^\circ C) \label{eq:56029b} \\ &Termisk diffusitivitet: \alpha=0.04 cm^2/s \label{eq:56029c} \end{align} $$

Vi har valgt verdiene slik at forholdet $ \dfrac{\bar h \cdot b}{k}=1 $, noe som fører til en enklere analytisk løsning. Verdiene som er angitt i \eqref{eq:56029a}, \eqref{eq:56029b} og \eqref{eq:56029c}, passer bra for nikkel-legeringer.

Vi skal nå løse følgende problem med $ T=T(r,t) $: $$ \begin{equation} \label{eq:56030a} \frac{\partial T}{\partial t}=\alpha \left(\frac{\partial^2T}{\partial r^2}+\frac{2}{r}\frac{\partial T}{\partial r} \right) \end{equation} $$

Randbetingelser: $$ \begin{equation} \label{eq:56030b} \frac{\partial T}{\partial r}(0,t)=0,\ \text{(symmetribetingelse)} \end{equation} $$

For $ r=b: $ $$ \begin{equation} \label{eq:56030c} k\frac{\partial T}{\partial r}= \bar h \cdot (T_v-T_b) \end{equation} $$

Startbetingelse: $$ \begin{equation} \label{eq:56030d} T(r,0)=T_k \end{equation} $$

I dette eksemplet nøyer vi oss med bruk av FTCS-skjemaet, men beregningen kan lett utvides til $ \theta $-skjemaet som vist i eksempel (Oppstart av rørstrømning). Med $ r_j=\Delta r\cdot j,\ \Delta r= \frac{1}{N+1},\ j=0,1,\dots,N+1 $ og $ D=\alpha \frac{\Delta t}{(\Delta r)^2} $ får vi fra når $ u(r,t)\to T(r,t),\ \theta=0 $ og $ \lambda = 2 $: $$ \begin{equation} \label{eq:56031} T^{n+1}_j=(1-2D)\cdot T^n_j+D[(1-1/j)\cdot T^n_{j-1}+(1+1/j)\cdot T^n_{j+1}],\ j=1,2,\dots \end{equation} $$

Som i eksempel (Oppstart av rørstrømning), bruker vi to versjoner for symmetribetingelsen for $ r=0 $. 1)

Fra \eqref{eq:5607} med $ \lambda = 2 $: $$ \begin{equation} \label{eq:56032} T^{n+1}_0=(1-6D)\cdot T^n_0+6D\cdot T^n_1 \end{equation} $$

2) Fra \eqref{eq:5609}: $$ \begin{equation} \label{eq:56033} T^n_0=\frac{1}{3}(4T^n_1-T_2^n),\ \text{alle } n \end{equation} $$

Randbetingelsen for $ r=b $.

Diskretiserer \eqref{eq:56030b} med bruk av 2. ordens bakoverdifferanser: $$ \begin{equation*} k\cdot \left( \frac{3T^n_{N+1}-4T^n_N+T^n_{N-1}}{2\cdot\Delta r} \right) = \bar h\cdot (T_v-T_b) \text{ som løst m.h.p $T^n_{N+1}$ gir:} \end{equation*} $$ $$ \begin{equation} \label{eq:56034a} T^n_{N+1}=\frac{4T^n_N-T^n_{N-1}+2\delta\cdot T_v}{3+2\delta} \end{equation} $$ $$ \begin{equation} \label{eq:56034b} \text{der } \delta = \frac{\Delta r\cdot \bar h}{k} \end{equation} $$

Problem 4: coded example

Hint.

# src-ch5/sphere.py @ git@lrhgit/tkt4140/src/src-ch5/sphereAnimation.py;


#import matplotlib.pyplot as plt
from matplotlib.pyplot import *
from sphereAnimation import animateSphere
from math import pi, sin, exp

#change some default values to make plots more readable 
LNWDT=3; FNT=11
rcParams['lines.linewidth'] = LNWDT; rcParams['font.size'] = FNT

def analytic(r, t):
    fac1 = 0.5*pi/b
    sign = 1
    criteria = 1
    epsi = 1e-5
    s = 0
    m = -1
    f1, f2 = b*1e-5, b*(1-1e-5)
    if (r >= f1) and (r <= f2):
        while (criteria > epsi):
            m = m + 2
            bm = m*fac1
            bmr = r*bm
            temp = exp(-alfa*bm**2*t)/m**2
            term = sign*sin(bmr)*temp
            s = s + term
            sign = -sign
            criteria = abs(temp/s)
        Tvalue = Tv + 8*b*(Tk - Tv)*s/(pi**2*r)
    
    if r < f1:
        while (criteria > epsi):
            m = m + 2
            bm = m*fac1
            term = sign*exp(-alfa*bm**2*t)/m
            s = s + term
            sign = -sign
            criteria = abs(term/s)
        Tvalue = Tv + 4*(Tk - Tv)*s/pi
    
    if r > f2:
        while (criteria > epsi):
            m = m + 2;  
            bm = m*fac1;  
            term = exp(-alfa*bm**2*t)/m**2
            s = s + term
            criteria = abs(term/s)
            
        Tvalue = Tv + 8*(Tk - Tv)*s/pi**2
    
    return Tvalue


"""script calculationg the temperature variation in a sphere with
    radius b, cooling in water.
    Initial temperature = Tk
    Temperature in water = Tv = constant
    Calculations are done with FTCS scheme, with r - steplength dr
    and numerical Fourier number D given. 
    Using 2. order backward difference in BC for r = b = 5 cm
    rb = 1: Using separate equation
        for r = 0. Stable for D < 1/3
    rb = 2: Using 2. order forward differance
        for r = 0. stable for D < 1/3
    The analytical solution is calculated by method kanalyt. The condition
    H*b/k = 1 has been used, where h is the heat transfer number and k the conductivity.
"""

b = 5. # radius of sphere [cm]
N = 50 # number of parts
r_array = np.linspace(0, b, N+1)
alfa = 0.04 # diffusivity [cm^2/s
dr = b/N
D = 0.4 # numerical fourrier number
dtau = D*dr**2/alfa # time step [s]
tauend = 600 # stop time [s]
ntot = int(tauend/dtau) + 1
K = 0.1 # conductivity [W/cm/C]
H = 0.02 # heat transfer number [W/cm^2/C]
dra = dr*H/K
ta = np.zeros(N + 1) # initialize analytic solution
Tk = 300 # start temperature in the sphere
Tv = 20 # Temperature in the water

Told = Tk*np.ones(N + 1) # initial values
Tnew = Tk*np.ones(N + 1)
Tmatrix = np.zeros((ntot, N + 1)) # initialize solution matrix
Tmatrix[0,:] = Told
Tmatrix_analytic = np.zeros((ntot, N + 1))
Tmatrix_analytic[0,:] = Told
temp_analytic = np.zeros_like(r_array)
rb = 2


for k in range(1,ntot):
    tau = k*dtau
    for j in range(1,N):

        Tnew[j] = (1 - 2*D)*Told[j] + D*((1 - 1./j)*Told[j - 1] + (1 + 1./j)*Told[j + 1]) #FTCS
        
    if rb == 2:
        Tnew[0] = (4*Tnew[1]-Tnew[2])/3.
    elif rb == 1:
        Tnew[0] = (1 - 6*D)*Told[0] + 6*D*Told[1]

    Tnew[N] = (4*Tnew[N - 1] - Tnew[N - 2] + 2*dra*Tv)/(3. + 2*dra)
    
    Told = Tnew.copy()

    for i , r in enumerate(r_array):
        temp_analytic[i]=analytic(r, tau)
        
    Tmatrix_analytic[k:] = temp_analytic
    Tmatrix[k:]= Tnew
    

TsnapshotList = [0, 10, 100, 300, 600]
legendList = []

for Tsnapshot in TsnapshotList:
    Temp = Tmatrix[Tsnapshot,:]
    time = Tsnapshot*dtau
    plot(Temp)
    legendList.append(str(int(time)) + ' s')

for Tsnapshot in TsnapshotList:
    Temp = Tmatrix_analytic[Tsnapshot,:]
    time = Tsnapshot*dtau
    plot( Temp, 'k--')
    #legendList.append(str(int(time)) + ' s')   
 
legend(legendList,loc='best',frameon=False)
title('Temperature in Sphere')
#axis([0, 5, 150, 310])
xlabel('r [cm]')
ylabel('T [$^o C$]')
#savefig('sphere.pdf')

        
showanimation = True

if showanimation:
    animateSphere(r_array, Tmatrix, Tmatrix_analytic, ntot, dtau)
show()