HOME

PDF (letter size)

Analytical solution to diffusion-convection PDE in 1D

Nasser M. Abbasi

September 7, 2023   Compiled on September 7, 2023 at 10:05pm

This is a diffusion-convection PDE.

\begin {align} \frac {\partial F\left ( t,z\right ) }{\partial t} & =k\frac {\partial ^{2}F\left ( t,z\right ) }{\partial z^{2}}+v\frac {\partial F\left ( t,z\right ) }{\partial z}\tag {1}\\ t & >0\nonumber \\ 0 & <z<L\nonumber \end {align}

Where \(k\) is the diffusion constant and \(v\) is the convection speed. Boundary conditions are

\begin {align*} F\left ( t,0\right ) & =0\\ F\left ( t,L\right ) & =1 \end {align*}

Initial conditions are

\[ F\left ( 0,z\right ) =\left \{ \begin {array} [c]{cc}0 & 0\leq z<L\\ 1 & z=L \end {array} \right . \]

The first step is to convert the PDE to pure diffusion PDE using the transformation

\[ F\left ( t,z\right ) =A\left ( t,z\right ) u\left ( t,z\right ) \]

Substituting this back in (1) gives

\[ A_{t}u+Au_{t}=k\left ( A_{zz}u+2A_{z}u_{z}+Au_{zz}\right ) +v\left ( A_{z}u+Au_{z}\right ) \]

Dividing by \(A\) and simplifying

\begin {align} \frac {A_{t}u}{A}+u_{t} & =k\left ( \frac {A_{zz}u+2A_{z}u_{z}}{A}+u_{zz}\right ) +v\left ( \frac {A_{z}u}{A}+u_{z}\right ) \nonumber \\ u_{t} & =ku_{zz}+k\frac {A_{zz}-\frac {A_{t}}{k}+vA_{z}}{A}u+\left ( \frac {2kA_{z}+vA}{A}\right ) u_{z} \tag {2} \end {align}

To make (2) pure diffusion PDE, we want \begin {align} k\frac {A_{zz}-\frac {A_{t}}{k}+vA_{z}}{A}u & =0\tag {3}\\ \left ( \frac {2kA_{z}+vA}{A}\right ) u_{z} & =0 \tag {4} \end {align}

From (4) \(\left ( 2kA_{z}+vA\right ) u_{z}=0\) or \(2kA_{z}+vA=0\) or \(\frac {\partial A}{\partial z}+\frac {v}{2k}A=0\) which has the solution\begin {equation} A\left ( t,z\right ) =C\left ( t\right ) e^{-\frac {v}{2k}z} \tag {5} \end {equation}

From (3) we want \(k\left ( A_{zz}-\frac {A_{t}}{k}+vA_{z}\right ) =0\). Substituting the result just obtained for \(A\left ( t,z\right ) \) in (3) gives

\begin {align*} \frac {v^{2}}{4k^{2}}C\left ( t\right ) e^{-\frac {v}{2k}z}-\frac {dC\left ( t\right ) }{dt}\frac {1}{k}e^{-\frac {v}{2k}z}-\frac {v^{2}}{2k}C\left ( t\right ) e^{-\frac {v}{2k}z} & =0\\ \frac {v^{2}}{4k^{2}}C\left ( t\right ) -\frac {1}{k}C^{\prime }\left ( t\right ) -\frac {v^{2}}{2k}C\left ( t\right ) & =0\\ C^{\prime }\left ( t\right ) +\frac {v^{2}}{4k}C\left ( t\right ) & =0 \end {align*}

Hence \[ C\left ( t\right ) =C_{1}e^{\frac {-v^{2}}{4k}t}\]

For some constant \(C_{1}\). The constant \(C_{1}\) ends up canceling out at the very end. Hence we set it to \(1\) now instead of carrying along in all the derivation in order to simplify notations. Therefore \(C\left ( t\right ) =e^{\frac {-v^{2}}{4k}t}\). Substituting this into (5) gives the transformation function

\[ A\left ( t,z\right ) =e^{-\left ( \frac {v^{2}t}{4k}+\frac {vz}{2k}\right ) }\]

Using this in (2) gives the pure diffusion PDE to solve

\begin {equation} u_{t}=ku_{zz} \tag {6} \end {equation}

Converting the original boundary conditions from \(F\) to \(u\) gives

\begin {align*} F\left ( t,0\right ) & =0\\ A\left ( t,0\right ) u\left ( t,0\right ) & =0\\ e^{-\frac {v^{2}t}{4k}}u\left ( t,0\right ) & =0\\ u\left ( t,0\right ) & =0 \end {align*}

And

\begin {align*} F\left ( t,L\right ) & =1\\ A\left ( t,L\right ) u\left ( t,L\right ) & =1\\ e^{-\left ( \frac {v^{2}t}{4k}+\frac {vL}{2k}\right ) }u\left ( t,L\right ) & =1\\ u\left ( t,L\right ) & =e^{\left ( \frac {v^{2}t}{4k}+\frac {vL}{2k}\right ) } \end {align*}

And for the initial conditions

\begin {align*} F\left ( 0,z\right ) & =\left \{ \begin {array} [c]{cc}0 & 0\leq z<L\\ 1 & z=L \end {array} \right . \\ A\left ( 0,z\right ) u\left ( 0,z\right ) & =\left \{ \begin {array} [c]{cc}0 & 0\leq z<L\\ 1 & z=L \end {array} \right . \\ e^{-\frac {vz}{2k}}u\left ( 0,z\right ) & =\left \{ \begin {array} [c]{cc}0 & 0\leq z<L\\ 1 & z=L \end {array} \right . \\ u\left ( 0,z\right ) & =\left \{ \begin {array} [c]{cc}0 & 0\leq z<L\\ e^{\frac {vL}{2k}} & z=L \end {array} \right . \end {align*}

Therefore the new PDE to solve is

\[ u_{t}=ku_{zz}\]

With time varying boundary conditions

\begin {align*} u\left ( t,0\right ) & =0\\ u\left ( t,L\right ) & =e^{\left ( \frac {v^{2}t}{4k}+\frac {vL}{2k}\right ) } \end {align*}

And initial conditions

\[ u\left ( 0,z\right ) =\left \{ \begin {array} [c]{cc}0 & 0\leq z<L\\ e^{\frac {vL}{2k}} & z=L \end {array} \right . \]

To solve this using separation of variables, the boundary conditions has to be homogenous. Therefore we use standard method to handle this as follows. Let

\begin {equation} u\left ( t,z\right ) =\phi \left ( t,z\right ) +u_{E}\left ( t,z\right ) \tag {7} \end {equation}

Where \(u_{E}\left ( t,z\right ) \) is the steady state solution which needs to only satisfy the boundary conditions and \(\phi \left ( t,z\right ) \) satisfies the PDE but with homogenous boundary conditions. Therefore \begin {align*} u_{E}\left ( t,z\right ) & =u\left ( t,0\right ) +z\left ( \frac {u\left ( t,L\right ) -u\left ( t,0\right ) }{L}\right ) \\ u_{E}\left ( t,z\right ) & =\frac {z}{L}e^{\left ( \frac {v^{2}t}{4k}+\frac {vL}{2k}\right ) } \end {align*}

And (7) becomes

\[ u\left ( t,z\right ) =\phi \left ( t,z\right ) +\frac {z}{L}e^{\left ( \frac {v^{2}t}{4k}+\frac {vL}{2k}\right ) }\]

Substituting the above in (6) gives

\begin {align*} \frac {\partial \phi }{\partial t}+\frac {z}{L}\frac {v^{2}}{4k}e^{\left ( \frac {v^{2}t}{4k}+\frac {vL}{2k}\right ) } & =k\frac {\partial ^{2}\phi }{\partial z^{2}}\\ \frac {\partial \phi }{\partial t} & =k\frac {\partial ^{2}\phi }{\partial z^{2}}-\frac {z}{L}\frac {v^{2}}{4k}e^{\left ( \frac {v^{2}t}{4}+\frac {vL}{2}\right ) } \end {align*}

Or

\begin {equation} \frac {\partial \phi }{\partial t}=k\frac {\partial ^{2}\phi }{\partial z^{2}}+Q\left ( t,z\right ) \tag {8} \end {equation}

This is diffusion PDE with homogenous B.C. with source term \[ Q\left ( t,z\right ) =-\frac {d}{dt}u_{E}\left ( t,z\right ) \] Now we find \(\phi \left ( t,z\right ) \). Since this solution needs to satisfy homogenous boundary conditions, we know the solution to pure diffusion on bounded domain with source present is by given by the following eigenfunction expansion

\begin {equation} \phi \left ( t,z\right ) =\sum _{n=1}^{\infty }b_{n}\left ( t\right ) \sin \left ( \sqrt {\lambda _{n}}z\right ) \tag {8A} \end {equation}

Where eigenvalues are \(\lambda _{n}=\left ( \frac {n\pi }{L}\right ) ^{2}\) for \(n=1,2,\cdots \) and \(\sin \left ( \sqrt {\lambda _{n}}z\right ) \) are the eigenfunction. Substituting the above in (8) in order to obtain an ODE to solve for \(b_{n}\left ( t\right ) \) gives

\begin {equation} \sum _{n=1}^{\infty }b_{n}^{\prime }\left ( t\right ) \sin \left ( \sqrt {\lambda _{n}}z\right ) =k\sum _{n=1}^{\infty }-b_{n}\left ( t\right ) \lambda _{n}\sin \left ( \sqrt {\lambda _{n}}z\right ) +Q\left ( t,z\right ) \tag {9} \end {equation}

Expanding \(Q\left ( t,z\right ) \) in terms of eigenfunctions

\[ Q\left ( t,z\right ) =\sum _{n=1}^{\infty }q_{n}\left ( t\right ) \sin \left ( \sqrt {\lambda _{n}}z\right ) \]

Applying orthogonality

\begin {equation} \int _{0}^{L}Q\left ( t,z\right ) \sin \left ( \sqrt {\lambda _{n}}z\right ) dz=q_{n}\left ( t\right ) \frac {L}{2} \tag {9A} \end {equation}

But \begin {align*} \int _{0}^{a}Q\left ( t,z\right ) \sin \left ( \sqrt {\lambda _{n}}z\right ) dz & =\int _{0}^{L}-\frac {z}{L}\frac {v^{2}}{4k}e^{\left ( \frac {v^{2}t}{4k}+\frac {vL}{2k}\right ) }\sin \left ( \sqrt {\lambda _{n}}z\right ) dz\\ & =-\frac {v^{2}}{4kL}e^{\left ( \frac {v^{2}t}{4k}+\frac {vL}{2k}\right ) }\int _{0}^{L}z\sin \left ( \sqrt {\lambda _{n}}z\right ) dz\\ & =-\frac {v^{2}}{4kL}e^{\left ( \frac {v^{2}t}{4k}+\frac {vL}{2k}\right ) }\frac {\left ( -1\right ) ^{n+1}L^{2}}{n\pi }\\ & =\frac {\left ( -1\right ) ^{n}Lv^{2}e^{\left ( \frac {v^{2}t}{4k}+\frac {vL}{2k}\right ) }}{4kn\pi }\\ & =\frac {\left ( -1\right ) ^{n}v^{2}e^{\left ( \frac {v^{2}t}{4k}+\frac {vL}{2k}\right ) }}{4k\sqrt {\lambda _{n}}} \end {align*}

Hence from (9A) we find \[ q_{n}\left ( t\right ) =\frac {\left ( -1\right ) ^{n}v^{2}e^{\left ( \frac {v^{2}t}{4k}+\frac {vL}{2k}\right ) }}{2Lk\sqrt {\lambda _{n}}}\] Using the above in (9) gives

\begin {align*} \sum _{n=1}^{\infty }b_{n}^{\prime }\left ( t\right ) \sin \left ( \sqrt {\lambda _{n}}z\right ) & =k\sum _{n=1}^{\infty }-b_{n}\left ( t\right ) \lambda _{n}\sin \left ( \sqrt {\lambda _{n}}z\right ) +\sum _{n=1}^{\infty }q_{n}\left ( t\right ) \sin \left ( \sqrt {\lambda _{n}}z\right ) \\ b_{n}^{\prime }\left ( t\right ) +k\lambda _{n}b_{n}\left ( t\right ) & =q_{n}\left ( t\right ) \end {align*}

To solve the above ODE, the integrating factor is \(\mu =e^{k\lambda _{n}t}\), therefore

\begin {align*} b_{n}\left ( t\right ) e^{k\lambda _{n}t} & =\int _{0}^{t}q_{n}\left ( \tau \right ) e^{k\lambda _{n}\tau }d\tau +C_{n}\\ b_{n}\left ( t\right ) & =e^{-k\lambda _{n}t}\int _{0}^{t}q_{n}\left ( \tau \right ) e^{k\lambda _{n}\tau }d\tau +C_{n}e^{-k\lambda _{n}t}\\ & =\int _{0}^{t}q_{n}\left ( \tau \right ) e^{k\lambda _{n}\left ( \tau -t\right ) }d\tau +C_{n}e^{-k\lambda _{n}t} \end {align*}

Using the above in (7) gives

\begin {equation} u\left ( t,z\right ) =\frac {z}{L}e^{\left ( \frac {v^{2}t}{4k}+\frac {vL}{2k}\right ) }+\sum _{n=1}^{\infty }\left ( \int _{0}^{t}q_{n}\left ( \tau \right ) e^{k\lambda _{n}\left ( \tau -t\right ) }d\tau +C_{n}e^{-k\lambda _{n}t}\right ) \sin \left ( \sqrt {\lambda _{n}}z\right ) \tag {10} \end {equation}

\(C_{n}\) is now found from initial conditions. At \(t=0\) the above becomes

\[ u\left ( 0,z\right ) -\frac {z}{L}e^{\frac {vL}{2k}}=\sum _{n=1}^{\infty }C_{n}\sin \left ( \sqrt {\lambda _{n}}z\right ) \]

Applying orthogonality

\begin {align*} \int _{0}^{L}\left ( u\left ( 0,z\right ) -\frac {z}{L}e^{\frac {vL}{2k}}\right ) \sin \left ( \sqrt {\lambda _{m}}z\right ) dz & =\int _{0}^{L}\sum _{n=1}^{\infty }C_{n}\sin \left ( \sqrt {\lambda _{n}}z\right ) \sin \left ( \sqrt {\lambda _{m}}z\right ) dz\\ \int _{0}^{L}u\left ( 0,z\right ) \sin \left ( \sqrt {\lambda _{n}}z\right ) dz-\int _{0}^{L}\frac {z}{L}e^{\frac {vL}{2k}}\sin \left ( \sqrt {\lambda _{n}}z\right ) dz & =C_{n}\frac {L}{2} \end {align*}

But \(\int _{0}^{L}u\left ( 0,z\right ) \sin \left ( \sqrt {\lambda _{n}}z\right ) dz=0\) since \(u\left ( 0,z\right ) \) is zero everywhere except at the end point. And \[ -\int _{0}^{L}\frac {z}{L}e^{\frac {vL}{2k}}\sin \left ( \sqrt {\lambda _{n}}z\right ) dz=\frac {\left ( -1\right ) ^{n}e^{\frac {vL}{2k}}}{\sqrt {\lambda _{n}}}\] Therefore

\begin {align} \frac {\left ( -1\right ) ^{n}e^{\frac {vL}{2k}}}{\sqrt {\lambda _{n}}} & =C_{n}\frac {L}{2}\nonumber \\ C_{n} & =\frac {2}{L}\frac {\left ( -1\right ) ^{n}e^{\frac {vL}{2k}}}{\sqrt {\lambda _{n}}} \tag {10A} \end {align}

And the solution (10) becomes

\begin {equation} u\left ( t,z\right ) =\frac {z}{L}e^{\left ( \frac {v^{2}t}{4k}+\frac {vL}{2k}\right ) }+\sum _{n=1}^{\infty }\left ( \int _{0}^{t}q_{n}\left ( \tau \right ) e^{k\lambda _{n}\left ( \tau -t\right ) }d\tau +C_{n}e^{-k\lambda _{n}t}\right ) \sin \left ( \sqrt {\lambda _{n}}z\right ) \tag {11} \end {equation}

But \begin {align*} \int _{0}^{t}q_{n}\left ( \tau \right ) e^{k\lambda _{n}\left ( \tau -t\right ) }d\tau & =\int _{0}^{t}\frac {\left ( -1\right ) ^{n}v^{2}e^{\left ( \frac {v^{2}\tau }{4k}+\frac {vL}{2k}\right ) }}{2Lk\sqrt {\lambda _{n}}}e^{k\lambda _{n}\left ( \tau -t\right ) }d\tau \\ & =\frac {2\left ( -1\right ) ^{n}v^{2}e^{-k\lambda _{n}t+\frac {vL}{2k}}\left ( e^{k\lambda _{n}t+\frac {tv^{2}}{4k}}-1\right ) }{n\pi \left ( 4k^{2}\lambda _{n}+v^{2}\right ) } \end {align*}

Hence (11) becomes

\begin {equation} u\left ( t,z\right ) =\frac {z}{L}e^{\left ( \frac {v^{2}t}{4k}+\frac {vL}{2k}\right ) }+\sum _{n=1}^{\infty }\left ( \frac {2\left ( -1\right ) ^{n}v^{2}e^{-k\lambda _{n}t+\frac {vL}{2k}}\left ( e^{k\lambda _{n}t+\frac {tv^{2}}{4k}}-1\right ) }{n\pi \left ( 4k^{2}\lambda _{n}+v^{2}\right ) }+\frac {2}{L}\frac {\left ( -1\right ) ^{n}e^{\frac {vL}{2k}}}{\sqrt {\lambda _{n}}}e^{-k\lambda _{n}t}\right ) \sin \left ( \sqrt {\lambda _{n}}z\right ) \tag {12} \end {equation}

We now convert back to \(F\left ( t,z\right ) \).  Since \(F\left ( t,z\right ) =A\left ( t,z\right ) u\left ( t,z\right ) \) and\(\ A\left ( t,z\right ) =e^{-\left ( \frac {v^{2}t}{4}+\frac {vz}{2}\right ) }\) then the final solution is

\[ F\left ( t,z\right ) =e^{-\left ( \frac {v^{2}t}{4k}+\frac {vz}{2k}\right ) }u\left ( t,z\right ) \]

The following is animation of the solution for 30 seconds, side-by-side with numerical solution.

The following is the code used

ClearAll[t, z, L0, k, lam, n, v, f] 
L0 = 10; 
v = 1/2; 
k = 1/2; 
ode = D[f[t, z], t] == k*D[f[t, z], {z, 2}] + v*D[f[t, z], z]; 
ic = f[0, z] == Piecewise[{{1, z == L0}, {0, True}}]; 
bc = {f[t, 0] == 0, f[t, L0] == 1}; 
sol = NDSolve[{ode, bc, ic}, f, {t, 0, 100}, {z, 0, 10}]; 
(*analytical*) 
lam[n_] := (n^2*Pi^2)/L0^2; 
max = 100; 
u[t_, z_] := Module[{}, 
   Sum[((2*(-1)^n*v^2*Exp[(L0*v)/(2*k) - k*lam[n]*t]* 
      (Exp[(t*v^2)/(4*k) + k*lam[n]*t] - 1))/(n* 
      Pi*(4*k^2*lam[n] + v^2)) + 
      (2/L0)*(((-1)^n*Exp[(L0*v)/(2*k)])/Sqrt[lam[n]])* 
      Exp[(-lam[n])*k*t])* 
      Sin[Sqrt[lam[n]]*z], {n, 1, max}] + (z/L0)* 
      Exp[(v^2*t)/(4*k) + (v*L0)/(2*k)]]; 
 
fAnalytical[t_, z_] := Module[{}, Exp[-((v^2*t)/(4*k) + (v*z)/(2*k))]*u[t, z]]; 
 
Manipulate[Grid[{{Row[{"t=", t}], SpanFromLeft}, 
    {Plot[Evaluate[f[t, z] /. sol], {z, 0, 10}, 
     PlotRange -> {{0, 10}, {-0.2, 1}}, 
     ImageSize -> 300, PlotLabel -> "NDSolve solution"], 
     Plot[Evaluate[fAnalytical[t, z]], {z, 0, 10}, 
     PlotRange -> {{0, 10}, {-0.2, 1}}, 
     ImageSize -> 300, PlotLabel -> "Analytical solution", 
     PlotStyle -> Red]}}], 
 
{{t, 0, "t"}, 0, 30, 0.01}, 
TrackedSymbols :> {t}]
 

References