PDF (letter size)

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

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

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\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 diﬀusion 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 )$$ 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

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

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

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

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

• 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.