{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Diffusion in heterogeneous media\n", "
\n", "\n", "Diffusion in heterogeneous media normally implies a non-constant\n", "diffusion coefficient $\\alpha = \\alpha (x)$.\n", "A 1D diffusion model with such a variable diffusion coefficient reads" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "\\frac{\\partial u}{\\partial t} =\n", "\\frac{\\partial}{\\partial x}\\left( \\alpha (x) \\frac{\\partial u}{\\partial x}\n", "\\right) + f(x,t), \\quad x\\in (0,L),\\ t\\in (0,T],\n", "\\label{diffu:pde2} \\tag{1}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "u(x,0) = I(x), \\quad x\\in [0,L],\n", "\\label{diffu:pde2:ic:u} \\tag{2}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "u(0,t) = U_0, \\quad t>0,\n", "\\label{diffu:pde2:bc:0} \\tag{3}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "u(L,t) = U_L, \\quad t>0.\n", "\\label{diffu:pde2:bc:L} \\tag{4}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A short form of the diffusion equation with variable coefficients is\n", "$u_t = (\\alpha u_x)_x + f$.\n", "\n", "## Discretization\n", "
\n", "\n", "We can discretize ([1](#diffu:pde2)) by a $\\theta$-rule in time\n", "and centered differences in space:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\lbrack D_t u\\rbrack^{n+\\frac{1}{2}}_i = \\theta\\lbrack D_x(\\overline{\\dfc}^x\n", "D_x u) + f\\rbrack^{n+1}_i +\n", "(1-\\theta)\\lbrack D_x(\\overline{\\dfc}^x\n", "D_x u) + f\\rbrack^{n}_i\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Written out, this becomes" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{align*}\n", "\\frac{u^{n+1}_i-u^{n}_i}{\\Delta t} &=\n", "\\theta\\frac{1}{\\Delta x^2}\n", "(\\dfc_{i+\\frac{1}{2}}(u^{n+1}_{i+1} - u^{n+1}_{i})\n", "- \\dfc_{i-\\frac{1}{2}}(u^{n+1}_i - u^{n+1}_{i-1})) +\\\\\n", "&\\quad (1-\\theta)\\frac{1}{\\Delta x^2}\n", "(\\dfc_{i+\\frac{1}{2}}(u^{n}_{i+1} - u^{n}_{i})\n", "- \\dfc_{i-\\frac{1}{2}}(u^{n}_i - u^{n}_{i-1})) +\\\\\n", "&\\quad \\theta f_i^{n+1} + (1-\\theta)f_i^{n},\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "where, e.g., an arithmetic mean can to be used for $\\dfc_{i+\\frac{1}{2}}$:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\dfc_{i+\\frac{1}{2}} = \\frac{1}{2}(\\dfc_i + \\dfc_{i+1})\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Implementation\n", "
\n", "\n", "Suitable code for solving the discrete equations is very similar to\n", "what we created for a constant $\\dfc$.\n", "Since the Fourier number has no meaning for varying\n", "$\\dfc$, we introduce a related parameter $D=\\Delta t /\\Delta x^2$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def solver_theta(I, a, L, Nx, D, T, theta=0.5, u_L=1, u_R=0,\n", " user_action=None):\n", " x = linspace(0, L, Nx+1) # mesh points in space\n", " dx = x[1] - x[0]\n", " dt = D*dx**2\n", " Nt = int(round(T/float(dt)))\n", " t = linspace(0, T, Nt+1) # mesh points in time\n", "\n", " u = zeros(Nx+1) # solution array at t[n+1]\n", " u_n = zeros(Nx+1) # solution at t[n]\n", "\n", " Dl = 0.5*D*theta\n", " Dr = 0.5*D*(1-theta)\n", "\n", " # Representation of sparse matrix and right-hand side\n", " diagonal = zeros(Nx+1)\n", " lower = zeros(Nx)\n", " upper = zeros(Nx)\n", " b = zeros(Nx+1)\n", "\n", " # Precompute sparse matrix (scipy format)\n", " diagonal[1:-1] = 1 + Dl*(a[2:] + 2*a[1:-1] + a[:-2])\n", " lower[:-1] = -Dl*(a[1:-1] + a[:-2])\n", " upper[1:] = -Dl*(a[2:] + a[1:-1])\n", " # Insert boundary conditions\n", " diagonal[0] = 1\n", " upper[0] = 0\n", " diagonal[Nx] = 1\n", " lower[-1] = 0\n", "\n", " A = scipy.sparse.diags(\n", " diagonals=[diagonal, lower, upper],\n", " offsets=[0, -1, 1],\n", " shape=(Nx+1, Nx+1),\n", " format='csr')\n", "\n", " # Set initial condition\n", " for i in range(0,Nx+1):\n", " u_n[i] = I(x[i])\n", "\n", " if user_action is not None:\n", " user_action(u_n, x, t, 0)\n", "\n", " # Time loop\n", " for n in range(0, Nt):\n", " b[1:-1] = u_n[1:-1] + Dr*(\n", " (a[2:] + a[1:-1])*(u_n[2:] - u_n[1:-1]) -\n", " (a[1:-1] + a[0:-2])*(u_n[1:-1] - u_n[:-2]))\n", " # Boundary conditions\n", " b[0] = u_L(t[n+1])\n", " b[-1] = u_R(t[n+1])\n", " # Solve\n", " u[:] = scipy.sparse.linalg.spsolve(A, b)\n", "\n", " if user_action is not None:\n", " user_action(u, x, t, n+1)\n", "\n", " # Switch variables before next step\n", " u_n, u = u, u_n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The code is found in the file [`diffu1D_vc.py`](${src_diffu}/diffu1D_vc.py).\n", "\n", "## Stationary solution\n", "
\n", "\n", "As $t\\rightarrow\\infty$, the solution of the\n", "problem ([1](#diffu:pde2))-([4](#diffu:pde2:bc:L))\n", "will approach\n", "a stationary limit where $\\partial u/\\partial t=0$. The governing\n", "equation is then" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "\\frac{d}{dx}\\left(\\alpha\\frac{du}{dx}\\right) =0,\n", "\\label{diffu:fd2:pde:st} \\tag{5}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "with boundary conditions $u(0)=U_0$ and $u(L)=U_L$.\n", "mathcal{I}_t is possible to obtain an exact solution of ([5](#diffu:fd2:pde:st))\n", "for any $\\alpha$. Integrating twice and applying the boundary conditions\n", "to determine the integration constants gives" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "u(x) = U_0 + (U_L-U_0)\\frac{\\int_0^x (\\alpha(\\xi))^{-1}d\\xi}{\\int_0^L (\\alpha(\\xi))^{-1}d\\xi}\n", "\\thinspace .\n", "\\label{diffu:fd2:pde:st:sol} \\tag{6}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Piecewise constant medium\n", "
\n", "\n", "Consider a medium built of $M$ layers. The layer boundaries\n", "are denoted $b_0, \\ldots, b_M$,\n", "where $b_0=0$ and $b_M=L$.\n", "If the layers potentially have different material properties, but\n", "these properties are constant within each layer, we can express $\\alpha$ as a\n", "*piecewise constant function* according to" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "\\alpha (x) = \\left\\lbrace\\begin{array}{ll}\n", "\\alpha_0,& b_0 \\leq x < b_1,\\\\\n", "\\vdots &\\\\\n", "\\alpha_i,& b_i \\leq x < b_{i+1},\\\\\n", "\\vdots &\\\\\n", "\\alpha_{M-1},& b_{M-1} \\leq x \\leq b_M.\n", "\\end{array}\\right.\n", "\\end{equation}\n", "\\label{diffu:fd2:pde:st:pc:alpha} \\tag{7}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The exact solution ([6](#diffu:fd2:pde:st:sol)) in case of such a\n", "piecewise constant $\\alpha$ function is easy to derive. Assume that\n", "$x$ is in the $m$-th layer: $x\\in [b_m, b_{m+1}]$. In the integral\n", "$\\int_0^x (a(\\xi))^{-1}d\\xi$ we must integrate through the first\n", "$m-1$ layers and then add the contribution from the remaining part\n", "$x-b_m$ into the $m$-th layer:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "u(x) = U_0 + (U_L-U_0)\n", "\\frac{\\sum_{j=0}^{m-1} (b_{j+1}-b_j)/\\alpha(b_j) + (x-b_m)/\\alpha(b_m)}{\\sum_{j=0}^{M-1} (b_{j+1}-b_j)/\\alpha(b_j)}\n", "\\label{diffu:fd2:pde:st:sol:pc} \\tag{8}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Remark.**\n", "mathcal{I}_t may sound strange to have a discontinuous $\\alpha$ in a differential\n", "equation where one is to differentiate, but a discontinuous $\\alpha$\n", "is compensated by a discontinuous $u_x$ such that $\\alpha u_x$ is\n", "continuous and therefore can be differentiated as $(\\alpha u_x)_x$.\n", "\n", "## Implementation of diffusion in a piecewise constant medium\n", "
\n", "\n", "Programming with piecewise function definitions quickly becomes\n", "cumbersome as the most naive approach is to test for which interval\n", "$x$ lies, and then start evaluating a formula like\n", "([8](#diffu:fd2:pde:st:sol:pc)). In Python, vectorized expressions may\n", "help to speed up the computations.\n", "The convenience classes `PiecewiseConstant` and\n", "`IntegratedPiecewiseConstant` in the [`Heaviside`](${src_diffu}/Heaviside.py)\n", "module were made to simplify programming with\n", "functions like ([7](#diffu:fd2:pde:st:pc:alpha)) and expressions like\n", "([8](#diffu:fd2:pde:st:sol:pc)). These utilities not only represent\n", "piecewise constant functions, but also *smoothed* versions of them\n", "where the discontinuities can be smoothed out in a controlled fashion.\n", "\n", "The `PiecewiseConstant` class is created by sending in the domain as a\n", "2-tuple or 2-list and a `data` object describing the boundaries\n", "$b_0,\\ldots,b_M$ and the corresponding function values\n", "$\\alpha_0,\\ldots,\\alpha_{M-1}$. More precisely, `data` is a nested\n", "list, where `data[i][0]` holds $b_i$ and `data[i][1]` holds the\n", "corresponding value $\\alpha_i$, for $i=0,\\ldots,M-1$. Given $b_i$ and\n", "$\\alpha_i$ in arrays `b` and `a`, it is easy to fill out the nested\n", "list `data`.\n", "\n", "In our application, we want to represent $\\alpha$ and $1/\\alpha$\n", "as piecewise constant functions, in addition to the $u(x)$ function\n", "which involves the integrals of $1/\\alpha$. A class creating the\n", "functions we need and a method for evaluating $u$, can take the\n", "form" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "class SerialLayers:\n", " \"\"\"\n", " b: coordinates of boundaries of layers, b[0] is left boundary\n", " and b[-1] is right boundary of the domain [0,L].\n", " a: values of the functions in each layer (len(a) = len(b)-1).\n", " U_0: u(x) value at left boundary x=0=b[0].\n", " U_L: u(x) value at right boundary x=L=b[0].\n", " \"\"\"\n", "\n", " def __init__(self, a, b, U_0, U_L, eps=0):\n", " self.a, self.b = np.asarray(a), np.asarray(b)\n", " self.eps = eps # smoothing parameter for smoothed a\n", " self.U_0, self.U_L = U_0, U_L\n", "\n", " a_data = [[bi, ai] for bi, ai in zip(self.b, self.a)]\n", " domain = [b[0], b[-1]]\n", " self.a_func = PiecewiseConstant(domain, a_data, eps)\n", "\n", " # inv_a = 1/a is needed in formulas\n", " inv_a_data = [[bi, 1./ai] for bi, ai in zip(self.b, self.a)]\n", " self.inv_a_func = \\\n", " PiecewiseConstant(domain, inv_a_data, eps)\n", " self.integral_of_inv_a_func = \\\n", " IntegratedPiecewiseConstant(domain, inv_a_data, eps)\n", " # Denominator in the exact formula is constant\n", " self.inv_a_0L = self.integral_of_inv_a_func(b[-1])\n", "\n", " def __call__(self, x):\n", " solution = self.U_0 + (self.U_L-self.U_0)*\\\n", " self.integral_of_inv_a_func(x)/self.inv_a_0L\n", " return solution" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A visualization method is also convenient to have. Below we plot $u(x)$\n", "along with $\\alpha (x)$ (which works well as long as $\\max \\alpha(x)$\n", "is of the same size as $\\max u = \\max(U_0,U_L)$)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "\n", "class SerialLayers:\n", " ...\n", "\n", " def plot(self):\n", " x, y_a = self.a_func.plot()\n", " x = np.asarray(x); y_a = np.asarray(y_a)\n", " y_u = self.u_exact(x)\n", " import matplotlib.pyplot as plt\n", " plt.figure()\n", " plt.plot(x, y_u, 'b')\n", " plt.hold('on') # Matlab style\n", " plt.plot(x, y_a, 'r')\n", " ymin = -0.1\n", " ymax = 1.2*max(y_u.max(), y_a.max())\n", " plt.axis([x[0], x[-1], ymin, ymax])\n", " plt.legend(['solution $u$', 'coefficient $a$'], loc='upper left')\n", " if self.eps > 0:\n", " plt.title('Smoothing eps: %s' % self.eps)\n", " plt.savefig('tmp.pdf')\n", " plt.savefig('tmp.png')\n", " plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "[Figure](#diffu:fd2:pde:st:sol:pc:fig1) shows the case where" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "b = [0, 0.25, 0.5, 1] # material boundaries\n", "a = [0.2, 0.4, 4] # material values\n", "U_0 = 0.5; U_L = 5 # boundary conditions" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "
\n", "\n", "

Solution of the stationary diffusion equation corresponding to a piecewise constant diffusion coefficient.

\n", "\n", "\n", "\n", "\n", "\n", "By adding the `eps` parameter to the constructor of the `SerialLayers`\n", "class, we can experiment with smoothed versions of $\\alpha$ and see\n", "the (small) impact on $u$. [Figure](#diffu:fd2:pde:st:sol:pc:fig2)\n", "shows the result.\n", "\n", "\n", "\n", "
\n", "\n", "

Solution of the stationary diffusion equation corresponding to a smoothed piecewise constant diffusion coefficient.

\n", "\n", "\n", "\n", "\n", "\n", "## Axi-symmetric diffusion\n", "
\n", "\n", "\n", "Suppose we have a diffusion process taking place in a straight tube\n", "with radius $R$. We assume axi-symmetry such that $u$ is just a\n", "function of $r$ and $t$, with $r$ being the radial distance from the center\n", "axis of the tube to a point. With such axi-symmetry it is\n", "advantageous to introduce *cylindrical coordinates* $r$, $\\theta$, and\n", "$z$, where $z$ is in the direction of the tube and $(r,\\theta)$ are\n", "polar coordinates in a cross section. Axi-symmetry means that all\n", "quantities are independent of $\\theta$. From the relations $x=\\cos\\theta$,\n", "$y=\\sin\\theta$, and $z=z$, between Cartesian and cylindrical coordinates,\n", "one can (with some effort) derive the diffusion equation in cylindrical\n", "coordinates, which with axi-symmetry takes the form" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\frac{\\partial u}{\\partial t} = \\frac{1}{r}\\frac{\\partial}{\\partial r}\n", "\\left(r\\dfc(r,z)\\frac{\\partial u}{\\partial r}\\right) + \\frac{\\partial}{\\partial z}\n", "\\left(\\alpha(r,z)\\frac{\\partial u}{\\partial z}\\right) + f(r,z,t)\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us assume that $u$ does not change along the tube axis so it\n", "suffices to compute variations in a cross section. Then $\\partial u/\\partial\n", "z = 0$ and we have a 1D diffusion equation in the radial coordinate\n", "$r$ and time $t$. In particular, we shall address the initial-boundary\n", "value problem" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "\\frac{\\partial u}{\\partial t} = \\frac{1}{r}\\frac{\\partial}{\\partial r}\n", "\\left(r\\dfc(r)\\frac{\\partial u}{\\partial r}\\right) + f(t), r\\in (0,R),\\ t\\in (0,T],\n", "\\label{diffu:fd2:radial:PDE} \\tag{9}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "\\frac{\\partial u}{\\partial r}(0,t) = 0, t\\in (0,T],\n", "\\label{diffu:fd2:radial:symmr0} \\tag{10}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "u(R,t) = 0, t\\in (0,T],\n", "\\label{diffu:fd2:radial:uR} \\tag{11}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "u(r,0) = I(r), r\\in [0,R].\n", "\\label{diffu:fd2:radial:initial} \\tag{12}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The condition ([10](#diffu:fd2:radial:symmr0)) is a necessary symmetry condition\n", "at $r=0$, while ([11](#diffu:fd2:radial:uR)) could be any Dirichlet\n", "or Neumann condition (or Robin condition in case of cooling or heating).\n", "\n", "The finite difference approximation will need the discretized version\n", "of the PDE for $r=0$ (just as we use the PDE at the boundary when\n", "implementing Neumann conditions). However, discretizing the PDE at\n", "$r=0$ poses a problem because of the $1/r$ factor. We therefore need\n", "to work out the PDE for discretization at $r=0$ with care.\n", "Let us, for the case of constant $\\dfc$, expand the spatial derivative term to" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\alpha\\frac{\\partial^2 u}{\\partial r^2} + \\alpha\\frac{1}{r}\\frac{\\partial u}{\\partial r}\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The last term faces a difficulty at $r=0$, since it becomes a $0/0$ expression\n", "caused by the symmetry condition at $r=0$.\n", "However, L'Hosptial's rule can be used:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\lim_{r\\rightarrow 0} \\frac{1}{r}\\frac{\\partial u}{\\partial r}\n", "= \\frac{\\partial^2 u}{\\partial r^2}\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The PDE at $r=0$ therefore becomes" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "\\frac{\\partial u}{\\partial t} = 2\\dfc\\frac{\\partial^2 u}{\\partial r^2}\n", "+ f(t)\\thinspace .\n", "\\label{diffu:fd2:radial:eq_PDEr0:aconst} \\tag{13}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For a variable coefficient $\\dfc(r)$ the expanded spatial derivative term reads" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\dfc(r)\\frac{\\partial^2 u}{\\partial r^2} +\n", "\\frac{1}{r}(\\dfc(r) + r\\dfc'(r))\\frac{\\partial u}{\\partial r}\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We are interested in this expression for $r=0$. A necessary condition\n", "for $u$ to be axi-symmetric is that all input data, including $\\alpha$,\n", "must also be axi-symmetric, implying that $\\alpha'(0)=0$ (the second\n", "term vanishes anyway because of $r=0$). The limit of interest is" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\lim_{r\\rightarrow 0}\n", "\\frac{1}{r}\\dfc(r)\\frac{\\partial u}{\\partial r} =\n", "\\dfc(0)\\frac{\\partial^2 u}{\\partial r^2}\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The PDE at $r=0$ now looks like" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "\\frac{\\partial u}{\\partial t} = 2\\dfc(0)\n", "\\frac{\\partial^2 u}{\\partial r^2}\n", "+ f(t),\n", "\\label{diffu:fd2:radial:eq_PDEr0:avar} \\tag{14}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "so there is no essential difference between the constant coefficient\n", "and variable coefficient cases.\n", "\n", "The second-order derivative in ([13](#diffu:fd2:radial:eq_PDEr0:aconst))\n", "and ([14](#diffu:fd2:radial:eq_PDEr0:avar))\n", "is discretized in the usual way." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "2\\dfc\\frac{\\partial^2}{\\partial r^2}u(r_0,t_n) \\approx\n", "[2\\dfc D_rD_r u]^n_0 =\n", "2\\dfc \\frac{u^{n}_{1} - 2u^{n}_0 + u^n_{-1}}{\\Delta r^2}\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The fictitious value $u^n_{-1}$ can be eliminated using the discrete\n", "symmetry condition" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "[D_{2r} u =0]^n_0 \\quad\\Rightarrow\\quad u^n_{-1} = u^n_1,\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "which then gives the modified approximation to the term with the second-order derivative\n", "of $u$ in $r$ at $r=0$:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "4\\dfc \\frac{u^{n}_{1} - u^{n}_0}{\\Delta r^2}\\thinspace .\n", "\\label{_auto1} \\tag{15}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The discretization of the term with the second-order derivative in $r$ at any\n", "internal mesh point is straightforward:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{align*}\n", "\\left[\\frac{1}{r}\\frac{\\partial}{\\partial r}\n", "\\left(r\\dfc\\frac{\\partial u}{\\partial r}\\right)\\right]_i^n\n", "& \\approx [r^{-1} D_r (r \\dfc D_r u)]_i^n\\\\\n", "&= \\frac{1}{r_i}\\frac{1}{\\Delta r^2}\\left(\n", "r_{i+\\frac{1}{2}}\\dfc_{i+\\frac{1}{2}}(u_{i+1}^n - u_i^n) - r_{i-\\frac{1}{2}}\\dfc_{i-\\frac{1}{2}}(u_{i}^n - u_{i-1}^n)\\right)\\thinspace .\n", "\\end{align*}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To complete the discretization, we need a scheme in time, but that can\n", "be done as before and does not interfere with the discretization in space.\n", "\n", "\n", "## Spherically-symmetric diffusion\n", "
\n", "\n", "### Discretization in spherical coordinates\n", "\n", "Let us now pose the problem from the section [Axi-symmetric diffusion](#diffu:fd2:radial)\n", "in spherical coordinates, where $u$ only depends on the radial coordinate\n", "$r$ and time $t$. That is, we have spherical symmetry.\n", "For simplicity we restrict the diffusion coefficient $\\dfc$ to be\n", "a constant. The PDE reads" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "\\frac{\\partial u}{\\partial t} = \\frac{\\dfc}{r^\\gamma}\\frac{\\partial}{\\partial r}\n", "\\left(r^\\gamma\\frac{\\partial u}{\\partial r}\\right) + f(t),\n", "\\label{_auto2} \\tag{16}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "for $r\\in (0,R)$ and $t\\in (0,T]$. The parameter $\\gamma$ is 2 for\n", "spherically-symmetric problems and 1 for axi-symmetric problems.\n", "The boundary and initial conditions\n", "have the same mathematical form as\n", "in ([9](#diffu:fd2:radial:PDE))-([12](#diffu:fd2:radial:initial)).\n", "\n", "Since the PDE in spherical coordinates has the same form as the PDE\n", "in the section [Axi-symmetric diffusion](#diffu:fd2:radial), just with the $\\gamma$ parameter\n", "being different, we can use the same discretization approach.\n", "At the origin $r=0$ we get problems with the term" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\frac{\\gamma}{r}\\frac{\\partial u}{\\partial t},\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "but L'Hosptial's rule shows that this term equals $\\gamma\\partial^2 u/\n", "\\partial r^2$, and the PDE at $r=0$ becomes" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "\\frac{\\partial u}{\\partial t} = (\\gamma+1)\\dfc\\frac{\\partial^2 u}{\\partial r^2}\n", "+ f(t)\\thinspace .\n", "\\label{_auto3} \\tag{17}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The associated discrete form is then" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "[D_t u = \\frac{1}{2} (\\gamma+1)\\dfc D_rD_r \\overline{u}^t + \\overline{f}^t]^{n+\\frac{1}{2}}_i,\n", "\\label{_auto4} \\tag{18}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "for a Crank-Nicolson scheme.\n", "\n", "### Discretization in Cartesian coordinates\n", "\n", "The spherically-symmetric spatial derivative can be transformed to\n", "the Cartesian counterpart by introducing" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "v(r,t) = ru(r,t)\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Inserting $u=v/r$ in" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\frac{1}{r^2}\\frac{\\partial}{\\partial r}\n", "\\left(\\dfc(r)r^2\\frac{\\partial u}{\\partial r}\\right),\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "yields" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "r\\left(\\frac{d \\dfc}{dr}\\frac{\\partial v}{\\partial r} +\n", "\\dfc\\frac{\\partial^2 v}{\\partial r^2}\\right) - \\frac{d \\dfc}{dr}v\n", "\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The two terms in the parenthesis can be combined to" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "r\\frac{\\partial}{\\partial r}\\left( \\dfc\\frac{\\partial v}{\\partial r}\\right)\\thinspace .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The PDE for $v$ takes the form" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "\\frac{\\partial v}{\\partial t} = \\frac{\\partial}{\\partial r}\\left( \\dfc\n", "\\frac{\\partial v}{\\partial r}\\right) - \\frac{1}{r}\\frac{d\\dfc}{dr}v + rf(r,t),\n", "\\quad r\\in (0,R),\\ t\\in (0,T]\\thinspace .\n", "\\label{_auto5} \\tag{19}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For $\\alpha$ constant we immediately realize that we can reuse a\n", "solver in Cartesian coordinates to compute $v$. With variable $\\alpha$,\n", "a \"reaction\" term $v/r$ needs to be added to the Cartesian solver.\n", "The boundary condition $\\partial u/\\partial r=0$ at $r=0$, implied\n", "by symmetry, forces $v(0,t)=0$, because" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\frac{\\partial u}{\\partial r} = \\frac{1}{r^2}\\left(\n", "r\\frac{\\partial v}{\\partial r} - v\\right) = 0,\\quad r=0\\thinspace .\n", "$$" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.2" } }, "nbformat": 4, "nbformat_minor": 4 }