HOME
PDF (letter size)
PDF (legal size)

## Analytical solution to diﬀusion-convection PDE in 1D

May 23, 2020   Compiled on May 23, 2020 at 1:46am

This is a diﬀusion-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 diﬀusion 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 ﬁrst step is to convert the PDE to pure diﬀusion 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 diﬀusion 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$$A\left ( t,z\right ) =C\left ( t\right ) e^{-\frac{v}{2k}z} \tag{5}$$

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 diﬀusion PDE to solve

$$u_{t}=ku_{zz} \tag{6}$$

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

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

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 )$$ satisﬁes 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

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

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

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

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

$$\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}$$

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

$$\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}$$

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 ﬁnd $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

$$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}$$

$$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

$$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}$$

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

$$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}$$

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 ﬁnal 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

References

• Paper ”Analytical Solution to the One-Dimensional Advection-Diﬀusion Equation with Temporally Dependent Coeﬃcients”. Dilip Kumar Jaiswal, Atul Kumar, Raja Ram Yadav. Journal of Water Resource and Protection, 2011, 3, 76-84
• Paper ”An analytical solution of the diﬀusion convection equation over a ﬁnite domain”. Mohammad Farrukh N. Mohsen and Mohammed H. Baluch, Appl. Math. Modelling, 1983, Vol. 7, August 285.
• Lecture 20: Heat conduction with time dependent boundary conditions using Eigenfunction Expansions. Introductory lecture notes on Partial Diﬀerential Equations By Anthony Peirce.