URI:
       tdiffusion.tex - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
  HTML git clone git://src.adamsgaard.dk/pism
   DIR Log
   DIR Files
   DIR Refs
   DIR LICENSE
       ---
       tdiffusion.tex (3521B)
       ---
            1 \documentclass{article}
            2 \usepackage[margin=1in]{geometry}
            3 \usepackage{amsmath}
            4 \usepackage{hyperref}
            5 \parindent=0in \parskip=0.5\baselineskip
            6 
            7 \usepackage{fancybox}
            8 \usepackage{xcolor}
            9 \newcommand{\highlight}[1]{{\color{red!80!black} \fbox{$ \displaystyle #1 $} }}
           10 
           11 \begin{document}
           12 
           13 \title{Exact solution for a simple IBVP with Dirichlet and Neumann BC}
           14 \maketitle
           15 
           16 \section{Introduction}
           17 \label{sec:introduction}
           18 
           19 This exact solution can be used to test the implementation of the
           20 enthalpy and temperature column solvers, both in the ice and in the
           21 bedrock.
           22 
           23 This is sufficient to test the bedrock thermal unit code, but only
           24 goes half way towards checking the enthalpy code: in these notes I
           25 consider a simple \emph{diffusion only} problem without advection or
           26 reaction terms.
           27 
           28 \section{The problem}
           29 \label{sec:problem}
           30 
           31 \begin{align}
           32   \label{eq:1}
           33   \text{PDE:} \qquad & u_{t} = \alpha^{2}\, u_{zz},\quad 0 < z < L,\\
           34   \label{eq:2}
           35   \text{BCs:} \qquad &\left\{
           36                        \begin{aligned}
           37                          u(0,t) &= U_{0},\\
           38                          u_{z}(L,t) &= Q_{L},
           39                        \end{aligned}\right.\\
           40   \label{eq:3}
           41   \text{IC:} \qquad & u(z,0) = \phi(z).
           42 \end{align}
           43 
           44 In the PISM context
           45 \begin{equation}
           46   \label{eq:14}
           47   \alpha^{2} = \frac{k}{\rho c_{p}},
           48 \end{equation}
           49 where $k$ is the thermal conductivity, $\rho$ is the density, and
           50 $c_{p}$ is the specific heat capacity of ice.
           51 
           52 To transform this problem with \emph{non-homogeneous} boundary conditions
           53 into one with homogeneous ones I define $v(z,t)$ by subtracting the
           54 steady state solution from $u$:
           55 \begin{equation}
           56   \label{eq:4}
           57   v(z,t) = u(z,t) - \left( U_{0} + Q_{L}\cdot z \right).
           58 \end{equation}
           59 
           60 It is easy to check that $v_{t} = u_{t}$, $v_{zz} = u_{zz}$, and $v$ satisfies
           61 
           62 \begin{align}
           63 \label{eq:15}
           64   \text{PDE:} \qquad & v_{t} = \alpha^{2}\, v_{zz},\quad 0 < z < L,\\
           65   \label{eq:5}
           66   \text{BCs:} \qquad &\left\{
           67                        \begin{aligned}
           68                          v(0,t) &= 0,\\
           69                          v_{z}(L,t) &= 0,
           70                        \end{aligned}\right.\\
           71 \label{eq:16}
           72   \text{IC:} \qquad & v(z,0) = \phi(z) - [U_{0} + Q_{L}z].
           73 \end{align}
           74 
           75 This new problem (equations (\ref{eq:15})--(\ref{eq:16})) can be solved using separation of variables, as follows.
           76 
           77 Assume that $v(z,t)$ can be written as
           78 \begin{equation}
           79   \label{eq:6}
           80   v(z,t) = Z(z)\cdot T(t),
           81 \end{equation}
           82 then the PDE (\ref{eq:15}) becomes
           83 \begin{equation}
           84   \label{eq:7}
           85   ZT' = \alpha^{2}Z''T,
           86 \end{equation}
           87 or
           88 \begin{equation}
           89   \label{eq:8}
           90   \frac{T'}{\alpha^{2}T} = \frac{Z''}{Z} = -\lambda^{2},
           91 \end{equation}
           92 for some constant $\lambda$, and so
           93 \begin{align}
           94   \label{eq:9}
           95   T' + \lambda^{2}\alpha^{2}T &= 0,\\
           96   \label{eq:10}
           97   Z'' + \lambda^{2}Z &= 0.
           98 \end{align}
           99 
          100 Solving equations (\ref{eq:9}) and (\ref{eq:10}) gives
          101 \begin{align}
          102   T(t) &= C_{1}e^{-(\lambda\alpha)^{2}t},\\
          103   Z(z) &= C_{2}\sin(\lambda z) + C_{3}\cos(\lambda z),
          104 \end{align}
          105 or
          106 \begin{equation}
          107   \label{eq:11}
          108   v(z,t) = e^{-(\lambda\alpha)^{2}t}\left( A\sin(\lambda z) + B\cos(\lambda z) \right)
          109 \end{equation}
          110 Combining (\ref{eq:11}) with (\ref{eq:5}) gives
          111 \begin{align}
          112   B &= 0,\\
          113   \cos(\lambda\cdot L) &= 0,
          114 \end{align}
          115 and so fundamental solutions have the form
          116 \begin{align}
          117   \label{eq:12}
          118   v_{n}(z,t) &= \highlight{a_{n}\,e^{-(\lambda_{n}\alpha)^{2}t}\sin\left( \lambda_{n} z \right)},\\
          119   \lambda_{n} &= \highlight{\frac{1}{L} \left( \frac{\pi}{2} + n\,\pi \right),\quad n \in \mathbf{Z}}.
          120 \end{align}
          121 
          122 \end{document}