@@ -3,8 +3,8 @@ label{ch:advec}
33Wave (Chapter ref{ch:wave}) and diffusion (Chapter ref{ch:diffu})
44equations are solved reliably by finite difference methods. As soon as
55we add a first-order derivative in space, representing *advective*
6- transport, also known as *convective* transport, the numerics gets
7- more complicated, and intuitively attractive methods no longer work
6+ transport ( also known as *convective* transport) , the numerics gets
7+ more complicated and intuitively attractive methods no longer work
88well. We shall show how and why such methods fail and provide
99remedies. The present chapter builds on basic knowledge about finite
1010difference methods for diffusion and wave equations, including the
@@ -44,7 +44,7 @@ label{advec:1D:pde1:I}
4444\end{align}
4545!et
4646In (ref{advec:1D:pde1:u}), $v$ is a given parameter, typically reflecting
47- the velocity of transport of a quantity $u$ with a flow.
47+ the transport velocity of a quantity $u$ with a flow.
4848There is only one boundary condition (ref{advec:1D:pde1:I}) since
4949the spatial derivative is only first order in the PDE (ref{advec:1D:pde1:u}).
5050The information at $x=0$ and the initial condition get
@@ -78,7 +78,7 @@ We see this relation from
7878!bt
7979\begin{align}
8080u(i\Delta x, (n+1)\Delta t) &= I(i\Delta x - v(n+1)\Delta t) \nonumber \\
81- &= I((i-1)\Delta x - vn\Delta t - v\Delta t - \Delta x) \nonumber \\
81+ &= I((i-1)\Delta x - vn\Delta t - v\Delta t + \Delta x) \nonumber \\
8282&= I((i-1)\Delta x - vn\Delta t) \nonumber \\
8383&= u((i-1)\Delta x, n\Delta t), \nonumber
8484\end{align}
@@ -130,6 +130,8 @@ with $C$ as the Courant number
130130
131131=== Implementation ===
132132
133+ idx{`solver_FECS`} idx{Gaussian pulse} idx{cosine pulse! half-truncated}
134+
133135A solver function for our scheme goes as follows.
134136
135137@@@CODE src-advec/advec1D.py fromto: import numpy@def solver\(
@@ -176,9 +178,9 @@ the printout of $u$ values in the `plot` function reveals that
176178$u$ grows very quickly. We may reduce $\Delta t$ and make it
177179very small, yet the solution just grows.
178180Such behavior points to a bug in the code.
179- However, choosing a coarse mesh and performing a time step by
180- hand calculations produce the same numbers as in the code, so
181- it seems that the implementation is correct.
181+ However, choosing a coarse mesh and performing one time step by
182+ hand calculations produces the same numbers as the code, so
183+ the implementation seems to be correct.
182184The hypothesis is therefore that the solution is unstable.
183185
184186===== Analysis of the scheme =====
@@ -204,15 +206,15 @@ simplifies if we introduce the complex Fourier component
204206with
205207
206208!bt
207- \[ \Aex=Be^{-ikv\Delta t} = Be^{iCkx }\tp\]
209+ \[ \Aex=Be^{-ikv\Delta t} = Be^{-iCk\Delta x }\tp\]
208210!et
209211Note that $|\Aex| \leq 1$.
210212
211213It turns out that many schemes also allow a Fourier wave component as
212214solution, and we can use the numerically computed values of $\Aex$
213215(denoted $A$) to learn about the
214- quality of the scheme. Hence, to analyze the difference scheme we just
215- have implemented, we look at how it treats the Fourier component
216+ quality of the scheme. Hence, to analyze the difference scheme we have just
217+ implemented, we look at how it treats the Fourier component
216218
217219!bt
218220\[ u_q^n = A^n e^{ikq\Delta x}\tp\]
@@ -221,13 +223,13 @@ have implemented, we look at how it treats the Fourier component
221223Inserting the numerical component in the scheme,
222224
223225!bt
224- \[ [D_t^+ A e^{ikq\Delta x} + v D_{2x}A e^{ikq\Delta x} = 0]^n_i ,\]
226+ \[ [D_t^+ A e^{ikq\Delta x} + v D_{2x}A e^{ikq\Delta x} = 0]^n_q ,\]
225227!et
226228and making use of (ref{form:exp:fd1c:center})
227229results in
228230
229231!bt
230- \[ [e^{ikq\Delta x} (\frac{A-1}{\Delta t} + v\frac{1}{\Delta x}i\sin (k\Delta x)) = 0]^n_i ,\]
232+ \[ [e^{ikq\Delta x} (\frac{A-1}{\Delta t} + v\frac{1}{\Delta x}i\sin (k\Delta x)) = 0]^n_q ,\]
231233!et
232234which implies
233235
@@ -248,6 +250,8 @@ label{advec:1D:leapfrog}
248250
249251=== Method ===
250252
253+ idx{`advec1D.py`}
254+
251255Another explicit scheme is to do a ``leapfrog'' jump over $2\Delta t$ in
252256time and combine it with central differences in space:
253257
@@ -294,6 +298,8 @@ for n in range(0, Nt):
294298
295299=== Running a test case ===
296300
301+ idx{initial condition! triangular}
302+
297303Let us try a coarse mesh such that the smooth Gaussian initial condition
298304is represented by 1 at mesh node 1 and 0 at all other nodes. This
299305triangular initial condition should then be advected to the right.
@@ -327,6 +333,8 @@ MOVIE: [mov-advec/cosinehat/UP/C08_dt001.ogg] Advection of half a cosine functio
327333
328334=== Analysis ===
329335
336+ idx{stability} idx{CFL condition}
337+
330338We can perform a Fourier analysis again. Inserting the numerical
331339Fourier component in the Leapfrog scheme, we get
332340
@@ -373,12 +381,12 @@ of many schemes with the exact expression.
373381
374382===== Upwind differences in space =====
375383label{advec:1D:FTUP}
376- idx{upwind difference}
384+ idx{upwind difference} idx{diffusion! artificial}
377385
378386Since the PDE reflects transport of information along with a flow in
379387positive $x$ direction, when $v>0$, it could be natural to go (what is called)
380388upstream and not
381- downstream in a spatial derivative to collect information about the
389+ downstream in the spatial derivative to collect information about the
382390change of the function. That is, we approximate
383391
384392!bt
@@ -401,7 +409,7 @@ label{advec:1D:upwind}
401409written out as
402410
403411!bt
404- \[ u^{n+1} = u^n_i - C(u^{n}_{i}-u^{n}_{i-1}),\]
412+ \[ u^{n+1}_i = u^n_i - C(u^{n}_{i}-u^{n}_{i-1}),\]
405413!et
406414gives a generally popular and robust scheme that is stable if $C\leq 1$.
407415As with the Leapfrog scheme, it becomes exact if $C=1$, exactly as shown in
@@ -440,8 +448,8 @@ which means
440448For $C<1$ there is, unfortunately,
441449non-physical damping of discrete Fourier components, giving rise to reduced
442450amplitude of $u^n_i$ as in Figures ref{advec:1D:UP:fig1:C08}
443- and ref{advec:1D:UP:fig2:C08}. The damping
444- is this figure is seen to be quite severe. Stability requires $C\leq 1$.
451+ and ref{advec:1D:UP:fig2:C08}. The damping seen
452+ in this figure is quite severe. Stability requires $C\leq 1$.
445453
446454!bnotice Interpretation of upwind difference as artificial diffusion
447455One can interpret the upwind difference as extra, artificial diffusion
@@ -456,7 +464,7 @@ by a forward difference in time and centered differences in space,
456464!bt
457465\[ D^+_t u + vD_{2x} u = \nu D_xD_x u]^n_i,\]
458466!et
459- gives actually the upwind scheme (ref{advec:1D:upwind}) if
467+ actually gives the upwind scheme (ref{advec:1D:upwind}) if
460468$\nu = v\Delta x/2$. That is, solving the PDE $u_t + vu_x=0$
461469by centered differences in space and forward difference in time is
462470unsuccessful, but by adding some artificial diffusion $\nu u_{xx}$,
@@ -491,8 +499,8 @@ Normally, we can do this in the simple way that `u_1[0]` is updated as
491499
492500In some schemes we may need $u^{n}_{N_x+1}$ and $u^{n}_{-1}$. Periodicity
493501then means that these values are equal to $u^n_1$ and $u^n_{N_x-1}$,
494- respectively. For the upwind scheme it is sufficient to set
495- `u_1[0]=u_1[Nx]` at a new time level before computing `u[1]`, which ensures
502+ respectively. For the upwind scheme, it is sufficient to set
503+ `u_1[0]=u_1[Nx]` at a new time level before computing `u[1]`. This ensures
496504that `u[1]` becomes right and at the next time level `u[0]` at the current
497505time level is correctly updated.
498506For the Leapfrog scheme we must update `u[0]` and `u[Nx]` using the scheme:
@@ -723,6 +731,8 @@ The complete code is found in the file
723731===== A Crank-Nicolson discretization in time and centered differences in space =====
724732label{advec:1D:CN}
725733
734+ idx{dispersion relation}
735+
726736Another obvious candidate for time discretization is the Crank-Nicolson
727737method combined with centered differences in space:
728738
@@ -735,7 +745,7 @@ $\theta$-rule,
735745!bt
736746\[ [D_t u]^n_i + v\theta [D_{2x} u]^{n+1}_i + v(1-\theta)[D_{2x} u]^{n}_i = 0\tp\]
737747!et
738- This gives rise to an *implicit* scheme,
748+ When $\theta$ is different from zero, this gives rise to an *implicit* scheme,
739749
740750!bt
741751\[ u^{n+1}_i + \frac{\theta}{2} C (u^{n+1}_{i+1} - u^{n+1}_{i-1})
@@ -768,7 +778,7 @@ b_{N_x} &= u^n_0
768778\end{align*}
769779!et
770780
771- The dispersion relation follows from inserting $u^n_i = A^ne^{ikx}$
781+ The dispersion relation follows from inserting $u^n_q = A^ne^{ikx}$
772782and using the formula (ref{form:exp:fd1c:center}) for the spatial
773783differences:
774784
@@ -788,7 +798,7 @@ Figure ref{advec:1D:CN:fig:C08} depicts a numerical solution for $C=0.8$
788798and the Crank-Nicolson
789799with severe oscillations behind the main wave. These oscillations are
790800damped as the mesh is refined. Switching to the Backward Euler scheme
791- helps on the oscillations as they are removed , but the amplitude is
801+ removes the oscillations, but the amplitude is
792802significantly reduced. One could expect that the discontinuous derivative
793803in the initial condition of the half a cosine wave would make even
794804stronger demands on producing a smooth profile, but Figure ref{advec:1D:BE:fig:C08} shows that also here, Backward-Euler is capable of producing a
@@ -848,8 +858,7 @@ or written out,
848858This is the explicit Lax-Wendroff scheme.
849859
850860!bnotice Lax-Wendroff works because of artificial viscosity
851- We can immediately from the formulas above
852- see that the Lax-Wendroff method is nothing but
861+ From the formulas above, we notice that the Lax-Wendroff method is nothing but
853862a Forward Euler, central difference in space scheme, which we have shown
854863to be useless because of chronic instability, plus an artificial
855864diffusion term of strength $\half\Delta t v^2$. It means that we can take
@@ -859,8 +868,6 @@ is not sufficiently large to make schemes stable, so then we also add
859868artificial diffusion.
860869!enotice
861870
862- [hpl: Lax-Wendroff does work! It's even unstable for $C=1$.]
863-
864871#FIGURE: [fig-advec/gaussian_LW_C08, width=800 frac=1] Lax-Wendroff scheme, $C=0.8$, $\Delta t = 0.01$ (left) and $\Delta t=0.005$ (right). label{advec:1D:LW:fig:C08}
865872
866873#MOVIE: [fig-advec/
@@ -877,6 +884,8 @@ This means that $|A|=1$ and also that we have an exact solution if $C=1$!
877884===== Analysis of dispersion relations =====
878885label{advec:1D:disprel}
879886
887+ idx{`dispersion_analysis.py`}
888+
880889We have developed expressions for $A(C,p)$ in the exact solution
881890$u_q^n=A^ne^{ikq\Delta x}$ of the discrete equations. These
882891expressions are valuable for investigating the quality of the numerical
@@ -900,12 +909,12 @@ velocity of the wave. The Fourier component
900909\[ D^n e^{ik(x-ct)}\]
901910!et
902911has damping $D$ and wave velocity $c$. Let us express our $A$ in
903- polar form, $A = A_re^{i\phi}$, and insert this expression in
912+ polar form, $A = A_re^{- i\phi}$, and insert this expression in
904913our discrete component $u_q^n = A^ne^{ikq\Delta x} = A^ne^{ikx}$:
905914
906915!bt
907916\[
908- u^n_q = A_r^n e^{i\phi n} e^{ikx} = A_r^n e^{i(kx - n\phi)} =
917+ u^n_q = A_r^n e^{- i\phi n} e^{ikx} = A_r^n e^{i(kx - n\phi)} =
909918A_r^ne^{i(k(x - ct))},\]
910919!et
911920for
@@ -966,12 +975,11 @@ FIGURE: [fig-advec/disprel_C0_5_CN_BE_FE, width=800 frac=1] Dispersion relations
966975The total damping after some time $T=n\Delta t$ is reflected by
967976$A_r(C,p)^n$. Since normally $A_r<1$, the damping goes like
968977$A_r^{1/\Delta t}$ and approaches zero as $\Delta t\rightarrow 0$.
969- The only two ways to reduce the damping are to increase $C$ and the mesh
970- resolution.
978+ The only way to reduce damping is to increase $C$ and/or the mesh resolution.
971979
972980We can learn a lot from the dispersion relation plots. For example,
973981looking at the plots for $C=1$, the schemes LW, UP, and LF has no
974- amplitude reduction, but LF has a wrong phase velocity for the
982+ amplitude reduction, but LF has wrong phase velocity for the
975983shortest wave in the mesh. This wave does not (normally) have enough
976984amplitude to be seen, so for all practical purposes, there is no
977985damping or wrong velocity of the individual waves, so the total shape
@@ -993,15 +1001,17 @@ dampens these waves. The same effect is also present in the Lax-Wendroff
9931001scheme, but the damping of the intermediate waves is hardly present, so
9941002there is visible noise in the total signal.
9951003
996- We realize that there is more understanding of the
997- behavior of the schemes in the dispersion analysis compared with a
998- pure truncation error analysis. The latter just says Lax-Wendroff is
1004+ We realize that, compared to pure truncation error analysis, dispersion
1005+ analysis sheds more light on the behavior of the computational schemes.
1006+ Truncation analysis just says that Lax-Wendroff is
9991007better than upwind, because of the increased order in time, but
1000- most people would say upwind is the better by looking at the plots.
1008+ most people would say upwind is the better one when looking at the plots.
10011009
10021010======= One-dimensional stationary advection-diffusion equation =======
10031011label{advec:1D:stationary}
10041012
1013+ idx{transport phenomena} idx{boundary layer}
1014+
10051015Now we pay attention to a physical process where advection (or convection)
10061016is in balance with diffusion:
10071017
@@ -1030,8 +1040,8 @@ or heat transport. One can also view the equations as a simple model
10301040problem for the Navier-Stokes equations. With the chosen boundary
10311041conditions, the differential equation problem models the phenomenon of
10321042a *boundary layer*, where the solution changes rapidly very close to
1033- the boundary. This is a characteristic of many fluid flow problems and
1034- make strong demands to numerical methods. The fundamental numerical
1043+ the boundary. This is a characteristic of many fluid flow problems, which
1044+ makes strong demands to numerical methods. The fundamental numerical
10351045difficulty is related to non-physical oscillations of the solution
10361046(instability) if the first-derivative spatial term dominates over the
10371047second-derivative term.
@@ -1045,13 +1055,13 @@ number of parameters in the problem. We scale $x$ by $\bar x = x/L$,
10451055and $u$ by
10461056
10471057!bt
1048- \[ \bar u = \frac{u - u_0}{u_L-u_0 }\tp\]
1058+ \[ \bar u = \frac{u - U_0}{U_L-U_0 }\tp\]
10491059!et
10501060Inserted in the governing equation we get
10511061
10521062!bt
1053- \[ \frac{v(u_L-u_0 )}{L}\frac{d\bar u}{d\bar x} =
1054- \frac{\dfc(u_L-u_0 )}{L^2}\frac{d^2\bar u}{d\bar x^2},\quad
1063+ \[ \frac{v(U_L-U_0 )}{L}\frac{d\bar u}{d\bar x} =
1064+ \frac{\dfc(U_L-U_0 )}{L^2}\frac{d^2\bar u}{d\bar x^2},\quad
10551065\bar u(0)=0,\ \bar u(1)=1\tp\]
10561066!et
10571067Dropping the bars is common. We can then simplify to
@@ -1081,12 +1091,14 @@ diffusion. The exact solution reads
10811091!bt
10821092\[ \uex (x) = \frac{e^{x/\epsilon} - 1}{e^{1/\epsilon} - 1}\tp\]
10831093!et
1084- The forthcoming plots illustrates this function for various values of
1094+ The forthcoming plots illustrate this function for various values of
10851095$\epsilon$.
10861096
10871097===== A centered finite difference scheme =====
10881098label{advec:1D:stationary:fdm}
10891099
1100+ idx{viscous boundary layer}
1101+
10901102The most obvious idea to solve (ref{advec:1D:stat:pde1s}) is to apply
10911103centered differences:
10921104
@@ -1158,12 +1170,12 @@ centered difference
11581170we use the one-sided *upwind* difference
11591171
11601172!bt
1161- \[ \frac{du}{dx}_i\approx \frac{u_{i}-u_{i-1}}{2 \Delta x},\]
1173+ \[ \frac{du}{dx}_i\approx \frac{u_{i}-u_{i-1}}{\Delta x},\]
11621174!et
11631175in case $v>0$. For $v<0$ we set
11641176
11651177!bt
1166- \[ \frac{du}{dx}_i\approx \frac{u_{i+1}-u_{i}}{2 \Delta x},\]
1178+ \[ \frac{du}{dx}_i\approx \frac{u_{i+1}-u_{i}}{\Delta x},\]
11671179!et
11681180On compact operator notation form, our upwind scheme can be expressed
11691181as
@@ -1274,7 +1286,7 @@ and centered differences in space everywhere, so the scheme can be expressed as
12741286
12751287!bt
12761288\begin{equation}
1277- [D_t u + vD^-_{2x} u = \dfc \frac{v\Delta x}{2}) D_xD_x u + f]_i^n\tp
1289+ [D_t u + vD^-_{2x} u = \dfc \frac{v\Delta x}{2} D_xD_x u + f]_i^n\tp
12781290\end{equation}
12791291!et
12801292
@@ -1323,9 +1335,11 @@ The mass transport equation for a substance then reads
13231335\end{equation}
13241336!et
13251337
1326- ===== Transport of a heat =====
1338+ ===== Transport of heat in fluids =====
13271339label{advec:app:heat}
13281340
1341+ idx{material derivative} idx{stress}
1342+
13291343The derivation of the heat equation in Section ref{diffu:app:heat} is limited
13301344to heat transport in solid bodies. If we turn the attention to heat transport
13311345in fluids, we get a material derivative of the internal energy in
0 commit comments