HOME
PDF (letter size)
PDF (legal size)

Analytical solution to diffusion-convection PDE in 1D

Nasser M. Abbasi

July 7, 2017 compiled on — Friday July 07, 2017 at 08:35 PM

This is a diffusion-convection PDE.

pict

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

pict

Initial conditions are

         {
            0 0 ≤ z < L
F (0,z) =    1   z = L

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

F (t,z) = A (t,z)u (t,z)

Substituting this back in (1) gives

Atu + Aut = k (Azzu + 2Azuz  + Auzz)+  v(Azu + Auz )

Dividing by A  and simplifying

pict

To make (2) pure diffusion PDE, we want

pict

From (4) (2kAz + vA) uz = 0  or 2kAz + vA  = 0  or ∂A-+ -vA = 0
∂z   2k  which has the solution

A (t,z) = C (t) e− v2kz
(5)

From (3) we want   (              )
k  Azz − Akt+  vAz  = 0  . Substituting the result just obtained for A (t,z)  in (3) gives

pict

Hence

           −v2t
C (t) = C1e 4k

For some constant C1   . The constant C1   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          −4v2k t
C (t) = e  . Substituting this into (5) gives the transformation function

           ( 2    )
A (t,z) = e− v4kt + v2zk

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

ut = kuzz
(6)

Converting the original boundary conditions from F  to u  gives

pict

And

pict

And for the initial conditions

pict

Therefore the new PDE to solve is

u = ku
 t     zz

With time varying boundary conditions

pict

And initial conditions

         {  0   0 ≤ z < L
u(0,z) =    vL
           e 2k    z = L

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 (t,z) = ϕ(t,z)+ uE (t,z)
(7)

Where uE (t,z)  is the steady state solution which needs to only satisfy the boundary conditions and ϕ (t,z)  satisfies the PDE but with homogenous boundary conditions. Therefore

pict

And (7) becomes

                    (v2t vL)
u (t,z) = ϕ(t,z)+  ze  4k +2k
                  L

Substituting the above in (6) gives

pict

Or

        2
∂ϕ-= k∂--ϕ + Q (t,z)
∂t     ∂z2
(8)

This is diffusion PDE with homogenous B.C. with source term

Q (t,z) = − d-uE (t,z )
           dt

Now we find ϕ(t,z)  . 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

         ∑∞         ( ∘ ---)
ϕ (t,z) =    bn (t)sin    λnz
         n=1
(8A)

Where eigenvalues are      (n-π)2
λn =   L   for n = 1,2,⋅⋅⋅ and    (√ ---)
sin   λnz are the eigenfunction. Substituting the above in (8) in order to obtain an ODE to solve for bn(t)  gives

∑∞          (∘ ---)     ∞∑              ( ∘ ---)
    b′n (t) sin    λnz  = k    − bn(t)λnsin    λnz  + Q (t,z)
n=1                     n=1
(9)

Expanding Q (t,z )  in terms of eigenfunctions

         ∑∞         ( ∘ ---)
Q (t,z ) =    qn (t) sin    λnz
         n=1

Applying orthogonality

∫  L          (  ---)
    Q (t,z)sin  ∘ λ z  dz = q (t) L
  0               n         n   2
(9A)

But

pict

Hence from (9A) we find

            n 2 (v24kt+v2Lk)
q (t) = (−-1)-v-e√--------
 n          2Lk  λn

Using the above in (9) gives

pict

To solve the above ODE, the integrating factor is μ = ekλnt  , therefore

pict

Using the above in (7) gives

           (      )      ( ∫                            )
        -z  v24kt+v2Lk    ∞∑      t      kλn(τ− t)        − kλnt     (∘ --- )
u(t,z) = L e        +         qn(τ)e       dτ + Cne       sin    λnz
                      n=1   0
(10)

Cn  is now found from initial conditions. At t = 0  the above becomes

                 ∑∞       (∘ --- )
u (0,z)− z-ev2Lk =    Cn sin   λnz
         L       n=1

Applying orthogonality

pict

But ∫L          (√ ---)
 0 u(0,z)sin   λnz  dz = 0  since u(0,z)  is zero everywhere except at the end point. And

  ∫ L         (∘ --- )          n vL
−     ze vL2k-sin    λnz  dz = (−-1√)-e2k
   0  L                         λn

Therefore

pict

And the solution (10) becomes

           (v2t vL)   ∞∑  ( ∫ t                          )    (∘ --- )
u(t,z) =-ze -4k-+2k  +         qn(τ)ekλn(τ− t)dτ + Cne− kλnt  sin    λnz
        L             n=1   0
(11)

But

pict

Hence (11) becomes

                        (                  vL (      tv2   )                     )
           ( 2    )  ∑∞   2 (− 1)n v2e−kλnt+ 2k  ekλnt+-4k-− 1           n  vL           (∘ --- )
u(t,z) = ze  v4tk +v2Lk +   || ---------------------------------- + 2-(− 1√)-e2ke−kλnt|| sin    λnz
         L           n=1(           nπ (4k2λn + v2)            L     λn         )
(12)

We now convert back to F (t,z )  .  Since F (t,z) = A (t,z)u (t,z)  and          − (v2t+ vz)
A (t,z) = e  4   2 then the final solution is

          −(v2t+vz)
F (t,z) = e   4k  2k u(t,z)

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

  

The following is the code used

References