Propagation of a Cloud of Pollutant in River

Propagation of a Cloud of Pollutant in River#

Introduction#

The objective of this workshop is to obtain insight into the time integration and space discretization of the advection equation. In this respect, accuracy and stability play a major role.

As background to this assignment, it is necessary to have some knowledge on the numerical solution to the advection equation. Please read chapter Numerical methods of PDEs of the MUDE book beforehand.

We consider a cloud of pollutant that is propagating in the river. With this exercise we want to explore the numerical properties of spatial discretization to be applied, in particular central differences and the first order upwind scheme. In addition and related to this, we explore stability properties of the forward Euler scheme.

Specifically, we will carry out the following tasks:

  • implement central differences and the first order upwind scheme combined with forward Euler - Python programming

  • evaluate accuracy and stability by varying the time step and grid size - making plots and writing down the results

With respect to programming requirements, you will need to fill in a few missing pieces of the functions, but mostly you will adapt the values of a few Python variables to evaluate different aspects of the problem.

Specifications of the test case#

The assignment involves a river section of 60 km in which one finds a given cloud of pollutant at \(t=0\). Initially, the cloud of pollutant has either a block-shaped profile or a Gaussian bell curve with a maximum concentration of \(c^0 = 1\) kg/m\(^3\). The maximum width of the cloud is 4 km. The initial cloud is located within the river section such that it does not hit the left boundary of the domain, that is, at \(x=0\) the concentration of the pollutant is zero, \(c(0,t) = 0\). This cloud is propagating through the river with a flow speed of 0.5 m/s until it reach almost the right boundary of the domain, that is, the tail of the cloud does not hit that boundary. So, the concentration must remain zero all times at \(x = 60\) km.

Part 1: implementation of spatial discretizations#

Our first task is to complete the code below by implementing both central differences and the first order upwind scheme. In order to do this, we first carry out a pen-and-paper exercise.

Task 1.1

  1. Write down the advection equation for the above problem. Indicate the type of this PDE.

  2. Specify the initial condition. Write down the formulas for both block- and bell-shaped profiles in full.

  3. Determine the duration of the simulation. This duration should be chosen such that the cloud has not left the domain by the last time step.

Task 1.2

The discretization of the advection equation is based on the MOL approach. Explain in your own words what this method involves.

We study two schemes for solving the advection equation, namely, the FTCS scheme and the FTBS scheme. Write down these schemes.

We will now implement both schemes in the code below. The idea is that we will loop over each of the inner grid points in the computational domain, one at a time. So we are not going to set up a system of equations in a matrix form.

import numpy as np
import matplotlib.pylab as plt
%matplotlib inline
from ipywidgets import interact, fixed, widgets

Task 1.3

Use your schemes derived in Task 1.2 to complete the function advection_1D in the code cell below.

Note that other functions are complete. These functions initialize the problem, check your problem and also plot your results.

def initialize_1D(c0, L, Nx, T, Nt, square=True):
    
    dx = L/Nx
    dt = T/Nt
    
    x = np.linspace(0, L, Nx)

    # initialise cloud:
    #  - a square pulse or a Gaussian shape pulse
    if square:
        c_init = np.zeros(Nx)
        c_init[int(2000/dx):int(2000/dx)+int(4000/dx) + 1] = c0
    else:
        s = 4000/6
        c_init = np.exp(-((x-4000)**2)/(2*s*s) )

    c_all = np.zeros((Nt+1, Nx))
    c_all[0] = c_init
    return x, c_all, dx, dt

def advection_1D(c, dx, dt, u, Nx, central=True):

    c_new = np.zeros(Nx)
    for i in range(1, Nx-1):
        if central:
            # central difference scheme
            ### YOUR CODE HERE ###
        else:
            # first order upwind scheme
            ### YOUR CODE HERE ###
    return c_new
    
def run_simulation_1D(c0, u, L, Nx, T, Nt, central=True, square=True):
    
    x, c_all, dx, dt = initialize_1D(c0, L, Nx, T, Nt, square=square)
    
    for t in range(Nt):
        c = advection_1D(c_all[t], dx, dt, u, Nx, central=central)
        c_all[t + 1] = c
        
    return x, c_all
    
def plot_1D(x, c, step=0):
    fig = plt.figure()
    ax = plt.axes(xlim=(0, round(x.max())),
                  ylim=(-0.5, int(np.ceil(c[0].max())) + 0.5))  
    ax.plot(x, c[step], marker='.')
    plt.xlabel('$x$ [m]')
    plt.ylabel('concentration [kg/m$^3$]')
    plt.title('propagation in river')
    plt.show()
    
def plot_1D_all():
    play = widgets.Play(min=0, max=Nt-1, step=1, value=0,
                        interval=100, disabled=False)
    slider = widgets.IntSlider(min=0, max=Nt-1, step=1, value=0)
    widgets.jslink((play, 'value'), (slider, 'value'))
    
    interact(plot_1D, x=fixed(x), c=fixed(c_all), step=play)

    return widgets.HBox([slider])
    
def check_variables_1D():
    dx = L/Nx
    dt = T/Nt
    print('Current variables values:')
    print(f'  c0 [---]: {c0:0.2f}')
    print(f'  u  [m/s]: {u:0.2f}')
    print(f'  L  [ m ]: {L:0.1f}')
    print(f'  Nx [---]: {Nx:4d}')
    print(f'  T  [ s ]: {T:0.1f}')
    print(f'  Nt [---]: {Nt:4d}')
    print(f'  dx [ m ]: {dx:0.2e}')
    print(f'  dt [ s ]: {dt:0.2e}')
    print(f'Using central differences?: {central}')
    print(f'Using square init. cond.? : {square}')
    print(f'CFL: {u*dt/dx:.2e}')

Part 2: numerical experiments#

The following variables will be defined to set up the problem:

c0 = initial value of our "pulse" representing the maximum of c(x,t) [-]
u = flow velocity field [m/s]
L = length of the domain [m]
Nx = number of grid cells in the direction x
T = duration of the simulation [s]
Nt = number of time steps

There are also two flag variables:

  1. central allows you to switch between central differences and upwind discretization schemes, and

  2. square changes the pulse from a box to a smooth bell curve.

Task 2.1

Below, some variables have been set. Note that \(\Delta x = 100\) m. Fill in the missing variables by appropriate numbers. (Hint: use the CFL condition to estimate the time step \(\Delta t\).)

You should use the function provided check_variables_1D, prior to running a simulation to make sure you that your numerical parameters are correct.

c0 = 1.0
u = 0.5

L = 60000.
Nx = 600

T = ### YOUR CODE HERE ###
Nt = ### YOUR CODE HERE ###

central = True
square = True

check_variables_1D()

Run the cell below to carry out the simulation.

You may ignore any warning.

x, c_all = run_simulation_1D(c0, u, L, Nx, T, Nt, central, square)

Use the plotting function to check your initial values. It should look like a box shape right at the start of the domain.

plot_1D(x, c_all)

Task 2.2

  1. Run the cell below to see what happens to the box-shaped cloud when moving from the left to right. Explain your findings based on the theory.

  2. Adapt the shape of the cloud to make it more smooth and run again (see second cell below). What do you observe?

  3. Repeat the two runs again but this time with the first order upwind scheme. (See third cell below.)

plot_1D_all()

square=False
x, c_all = run_simulation_1D(c0, u, L, Nx, T, Nt, central, square)
plot_1D_all()

# re-define key variables for upwind and re-run with run_simulation_1D() and plot_1D_all()
### YOUR CODE HERE ###

Part 3: stability and numerical diffusion of FTBS#

With this exercise we explore the stability of the FTBS scheme. From theory, we found its stability condition as given by:

\[ \sigma = \frac{u \, \Delta t}{\Delta x} \leq 1 \]

In addition, the FTBS scheme produces some amount of numerical diffusion and is given by

\[ \kappa_a = -\frac 12 \, u \left ( u\,\Delta t - \Delta x \right ) \]

This amount consists of two contributions:

  • the forward Euler scheme produces some amount equal to \(-u^2\,\Delta t/2\) and

  • the first order upwind scheme produces some other amount, namely, \(u\,\Delta x/2\).

We observe that forward Euler decreases(!) the total amount of numerical diffusion while the upwind scheme increases that amount.

Task 3.1

  1. Derive above \(\kappa_a\) for the FTBS scheme using the Taylor series expansion. (Hint: write the result of substituting the expansions in the FTBS scheme as an advection-diffusion equation.) Finally, calculate \(\kappa_a\) based on your current numerical parameters.

  2. Observe how the results of the upwind scheme changes with different time steps. For example, what happens when the time step is doubled? Compute the Courant number.

  3. Follow up to previous question, how would you adapt the mesh size such that the solution becomes stable again? Compare the results with the ones obtained in Task 2.2. Look at the movement of the center of gravity and also at the shape. Explain your findings based on the notion of numerical diffusion.

  4. What happens when the time step is halved? (Use again \(\Delta x = 100\) m.) Compute again \(\kappa_a\). Compare with the previous results. What do you observe? Explain your findings.

  5. Compute the time step \(\Delta t\) such that the Courant number \(\sigma\) is precisely 1. Rerun the simulation. What do you observe? What is \(\kappa_a\)? Can you explain your observations?

# run the simulation with different time steps
### YOUR CODE HERE ###

By Marcel Zijlema, Delft University of Technology. CC BY 4.0, more info on the Credits page of Workbook.

Contents