The Convection Diffusion Equation

Document Sample
The Convection Diffusion Equation Powered By Docstoc
					                                                                          CHAPTER 2. PARABOLIC EQUATIONS

2.3 The Convection Diffusion Equation
                               The Convection­Diffusion Equation: ut + aux = buxx
   If a = 0, we have the heat equation, and if b = 0, we have the wave equation. Define a new function w
such that w(t, x − at) = u(t, x), i.e., w(t, x) = u(t, x + at). Then
                                     bwxx = buxx = ut + aux = wt − awx + awx = wt .
So u is simply the solution to the heat equation translated with speed a. The problem occurs when the
viscosity coefficient b is very small compared to a. Then the obvious numerical methods trick you. We have
a choice of
                         n+1
                        vm − vm   n     v n − vm−1 n                  n    n
                                                             v n − 2vm + vm−1
                                    + a m+1             = b m+1                ,                     (2.2)
                             k                 2h                    h2
                            n+1     n        n       n         n      n    n
                           vm − vm         v      − vm       v    − 2vm + vm−1
                                      + a m+1           = b m+1                ,
                                k               h                    h2
                            n+1
                           vm − vm  n              n
                                           v n − vm−1                 n    n
                                                             v n − 2vm + vm−1
                                      +a m              = b m+1                .                     (2.3)
                                k               h                    h2
We expand in Taylor series to see that the orders of accuracy are
                                           O(k + h2 ),   O(k + h), and O(k + h),
respectively, which speaks strongly in favor of (2.2).
    The heat equation has one nice property. The maximum at a later time is less than the maximum at an
earlier time. We will copy that. Rewrite (2.2) as
                               �          �                       �        �
                     n+1              aλ     n                n         aλ    n
                    vm     =     bµ −       vm+1 + (1 − 2bµ)vm + bµ +        vm−1
                                       2                                 2
                                  �         �                         �         �
                                         aλ     n               n           aλ     n
                           = bµ 1 −            vm+1 + (1 − 2bµ)vm + bµ 1 +        vm−1 .
                                        2bµ                                 2bµ
          aλ        ah                                                           aλ
Let α ≡   2bµ   =   2b .   Now if all the coefficients are positive (i.e., |α| = | 2bµ | < 1), then
                        n+1                        n                     n                     n
                      |vm |      ≤ bµ(1 − α) max |vm | + (1 − 2bµ) max |vm | + bµ(1 + α) max |vm |
                                                  m                     m                      m
                                                                             n
                                 ≤ [bµ(1 − α) + (1 − 2bµ) + bµ(1 + α)] max |vm |
                                                                             m
                                         n
                                 ≤ max |vm |.
                                       m

   The maximum is a decreasing function of time if |α| ≤ 1, i.e., if
                                        |a| h                                         2b
                                           · ≤ 1, which is the same as           h≤       .          (2.4)
                                         b 2                                          |a|
This will, of course, be satisfied eventually as h → 0, but who can wait that long? Say a = 10, b = 10−2 ⇒
h ≈ 10−3 is needed. And remember, for stability we must have 1/2 ≥ bµ = 10−2 · k/(10−3 )2 . This implies
k ≈ 10−4 /2, which is terribly small.
   Now look at (2.3) instead:
                                  n+1  n       n    n           n       n    n
                                 vm − vm + aλ(vm − vm−1 ) = bµ(vm+1 − 2vm + vm−1 ),
which we rewrite as
                               n+1                n                     n       n
                              vm      = (bµ + aλ)vm−1 + (1 − 2bµ − aλ)vm + bµvm+1
                                               aλ n                    aλ
                                      = bµ(1 +    )vm−1 + (1 − 2bµ(1 +               n
                                                                          ))v n + bµvm+1
                                               bµ                      2bµ m
                                                    n                     n        n
                                      = bµ(1 + 2α)vm−1 + (1 − 2bµ(1 + α))vm + bµvm+1 .
2.4. SUMMARY OF SCHEMES FOR THE HEAT EQUATION

Say a > 0, so α > 0. Now the requirement for max­norm stability becomes

                                       2bµ(1 + α) < 1   or   2bµ + aλ < 1,

which is a lot less restrictive than (2.4). We can pick h to be 10−2 instead of 10−3 , i.e., 10 times larger. We
          −3
try k = 105 , h = 10−2 , a = 10, b = 10−2 . Then
                                                             −3
                                                           10
                                    k      k                          2  1
                               2b     2
                                        + a = 2 · 10−2 ·     5
                                                             −2 )2
                                                                   =    + < 1.
                                    h      h             (10         50 5

So we increased the timestep by a factor of 10. But at what price?
   We can rewrite (2.3) as
                       n+1    n       n       n              � n      n    n
                                     vm+1 − vm−1          ah vm+1 − 2vm + vm−1
                                                     �
                      vm − vm
                                 +a                = b+                        .
                           k              2h               2         h2

                                       ah
We introduced an artificial viscosity    2   = bα. This artificial viscosity comes from our numerical method.
   Therefore solving ut + aux = b(1 + α)uxx using (2.2) is equivalent to solving ut + aux = buxx using (2.3).