Numerical simulation of tracer transport in a 1D channel with the FTBS upwind scheme

1. Aim of the simulation

In this exercise we want to simulate the transport of a tracer cloud along a straight channel.

We imagine releasing a quantity of dye at one point of the channel. The water current carries this concentration downstream with an approximately constant velocity.

The aim is to describe numerically how the spatial distribution of the concentration changes over time.

The main physical phenomenon is convection, that is, the transport of a quantity by an underlying flow.

In this first simulation we consider the simplest possible case:

the channel is one-dimensional;
the flow velocity is constant;
the tracer does not chemically react;
diffusion is neglected;
sources and sinks are neglected;
the concentration is only transported by the flow.

This leads to the one-dimensional linear convection equation.

2. Physical model

We denote by

u(x,t)

the tracer concentration at position x and time t.

The flow velocity is denoted by

c

c > 0

the transport occurs toward the right, that is, toward increasing values of x.

c < 0

the transport occurs toward the left.

In this exercise we consider the case

c > 0

therefore the tracer cloud moves from left to right.

The mathematical model is

\frac{\partial u}{\partial t} + c\frac{\partial u}{\partial x}=0

This is the one-dimensional linear convection equation.

The term

\frac{\partial u}{\partial t}

describes the time variation of the concentration.

The term

c\frac{\partial u}{\partial x}

describes the spatial transport of the concentration with velocity c.

The equation says that concentration is not created and not destroyed: it is simply transported along the channel.

If the initial concentration is a localized cloud, for example a droplet or a rectangular pulse, the ideal solution should be the same cloud translated in time.

If c is positive, after some time the cloud will be farther to the right.

Ideally, without diffusion and without numerical errors, the shape of the cloud should remain unchanged.

In other words:

u(x,t)=u(x-ct,0)

This formula means that the initial profile is transported rigidly with velocity c.

3. Spatial domain

We consider a channel of length L. The spatial domain is

0 \le x \le L

To solve the problem numerically, we divide the channel into a grid of points:

x_i=i\Delta x

where

i=0,1,2,\dots,N-1

and

\Delta x = \frac{L}{N-1}

The numerical value

u_i^n

represents the approximation of u(x_i, t_n), that is, the concentration at the spatial point x_i and time t_n.

4. Time discretization

Time is also discretized:

t_n=n\Delta t

where

n=0,1,2,\dots

and $\Delta t$ is the time step.

Therefore the problem is no longer solved in the continuous sense, but on a space-time grid.

The numerical state of the simulation at a given time is the vector

u^n = \begin{pmatrix} u_0^n \\ u_1^n \\ u_2^n \\ \vdots \\ u_{N-1}^n \end{pmatrix}

Each component represents the concentration at one point of the channel.

5. FTBS upwind scheme

The starting equation is

\frac{\partial u}{\partial t} + c\frac{\partial u}{\partial x}=0

For the time derivative we use a forward difference:

\frac{\partial u}{\partial t} \approx \frac{u_i^{n+1}-u_i^n}{\Delta t}

Since we consider $c > 0$ , information comes from the left. For this reason we use a backward difference in space:

\frac{\partial u}{\partial x} \approx \frac{u_i^n-u_{i-1}^n}{\Delta x}

Substituting these approximations into the convection equation gives

\frac{u_i^{n+1}-u_i^n}{\Delta t} + c \frac{u_i^n-u_{i-1}^n}{\Delta x} =0

Solving for $u_i^{n+1}$ gives

u_i^{n+1} = u_i^n - c\frac{\Delta t}{\Delta x} \left( u_i^n-u_{i-1}^n \right)

This is the FTBS scheme, namely:

Forward Time;
Backward Space.

It is also called an upwind scheme because it uses the value coming from the direction from which information arrives.

6. Courant number and CFL condition

We introduce the Courant number:

C = c\frac{\Delta t}{\Delta x}

The scheme becomes

u_i^{n+1} = u_i^n - C \left( u_i^n-u_{i-1}^n \right)

or equivalently

u_i^{n+1} = (1-C)u_i^n + C u_{i-1}^n

The Courant number measures how far the cloud moves during one time step compared with the size of one spatial cell.

C=1

the cloud moves exactly by one cell at every time step.

C<1

the cloud moves by less than one cell per time step.

C>1

the cloud moves by more than one cell per time step, and the scheme becomes unstable.

For the FTBS upwind scheme applied to the 1D linear convection equation with $c>0$ , the stability condition is

0 \le C \le 1

that is

0 \le c\frac{\Delta t}{\Delta x} \le 1

This condition is called the CFL condition, from Courant, Friedrichs and Lewy.

In practice, it means that the time step cannot be chosen arbitrarily large. It must satisfy

\Delta t \le \frac{\Delta x}{c}

This condition has a clear physical meaning: during a single time step, information must not cross more than one grid cell.

7. First step: defining the physical system

import numpy as np
import matplotlib.pyplot as plt

def cfl_convection(c, dt, dx):
    """Courant number for linear convection."""
    return abs(c) * dt / dx


# Channel length
L = 1.0

# Number of grid points
nx = 201

# Spatial grid
x = np.linspace(0.0, L, nx)

# Spatial step
dx = x[1] - x[0]

# Flow velocity
c = 1.0

# Desired Courant number
C_target = 0.8

# Time step chosen to satisfy the CFL condition
dt = C_target * dx / abs(c)

# Effective Courant number
C = cfl_convection(c, dt, dx)

print("dx =", dx)
print("dt =", dt)
print("Courant number =", C)

dx = 0.005
dt = 0.004
Courant number = 0.8

8. Initial concentration

# Initial concentration

u0 = np.zeros_like(x)

# Initial tracer cloud
u0[(x >= 0.15) & (x <= 0.30)] = 1.0

plt.figure(figsize=(8, 4))

plt.plot(x, u0, label="Initial concentration")

plt.xlabel("x")
plt.ylabel("u(x,0)")
plt.title("Initial tracer cloud in the channel")
plt.grid(True)
plt.legend()
plt.show()

Initial rectangular tracer cloud in the channel.

9. Second step: creating the numerical simulation method

The following function implements the FTBS upwind formula for positive velocity.

def linear_convection_upwind(u0, c, dx, dt, nt, left_bc=None):
    u = np.asarray(u0, dtype=float).copy()

    if c < 0:
        raise ValueError("This version assumes c > 0; reverse the stencil for c < 0.")

    C = c * dt / dx

    for _ in range(nt):
        un = u.copy()

        u[1:] = un[1:] - C * (un[1:] - un[:-1])

        if left_bc is not None:
            u[0] = left_bc(u, dx)

    return u

10. Third step: initial and boundary conditions

Since the flow is directed from left to right, the left boundary is an inflow boundary. In this example no new tracer enters from the left, so the boundary value is set to zero.

def left_boundary(u, dx):
    return 0.0

t_end = 0.5
nt = int(t_end / dt)

print("nt =", nt)
print("simulated final time =", nt * dt)

nt = 125
simulated final time = 0.5

11. Fourth step: running the simulation

u_final = linear_convection_upwind(
    u0,
    c,
    dx,
    dt,
    nt,
    left_bc=left_boundary
)

plt.figure(figsize=(8, 4))

plt.plot(x, u0, label="Initial concentration")
plt.plot(x, u_final, label="Final concentration")

plt.xlabel("x")
plt.ylabel("u(x,t)")
plt.title("Transport of a tracer cloud with the FTBS upwind scheme")
plt.grid(True)
plt.legend()
plt.show()

Comparison between the initial and final concentration profiles.

12. Evolution in time

To better visualize the transport process, we plot the concentration profile at several different times.

snapshot_times = [0.0, 0.1, 0.2, 0.3, 0.4, 0.5]

plt.figure(figsize=(9, 5))

for t_snap in snapshot_times:
    nt_snap = int(t_snap / dt)

    u_snap = linear_convection_upwind(
        u0,
        c,
        dx,
        dt,
        nt_snap,
        left_bc=left_boundary
    )

    plt.plot(x, u_snap, label=f"t = {nt_snap * dt:.2f}")

plt.xlabel("x")
plt.ylabel("u(x,t)")
plt.title("Evolution of the concentration in the channel")
plt.grid(True)
plt.legend()
plt.show()

Evolution of the rectangular tracer cloud over time.

13. Final remarks

The cloud moves toward the right, as expected from the positive value of the velocity c. The simulation remains stable because the Courant number is $C=0.8$ , which satisfies the CFL condition.

However, the sharp front becomes smoother as time passes. This is not physical diffusion, because diffusion was not included in the model. It is numerical diffusion, introduced by the upwind scheme.

Important observation. The FTBS upwind scheme is stable when the CFL condition is respected, but it tends to smooth sharp profiles. This is why the rectangular pulse becomes less sharp during the simulation.

14. Experiment with a different initial shape

We now repeat the experiment with a different initial tracer cloud. Instead of a rectangular pulse, we use a semicircular profile.

u(x,0)= \sqrt{1-\left(\frac{x-x_c}{R}\right)^2}

inside the interval $|(x-x_c)/R|\le 1$ , and zero outside.

# Initial concentration

u0 = np.zeros_like(x)

# Center and radius of the initial cloud
x_c = 0.225
R = 0.075

# Normalized distance from the center
s = (x - x_c) / R

# Initial cloud with a semicircular shape
mask = np.abs(s) <= 1.0
u0[mask] = np.sqrt(1.0 - s[mask]**2)

plt.figure(figsize=(8, 4))

plt.plot(x, u0, label="Initial concentration")

plt.xlabel("x")
plt.ylabel("u(x,0)")
plt.title("Initial semicircular tracer cloud in the channel")
plt.grid(True)
plt.legend()
plt.show()

Initial semicircular tracer cloud in the channel.

We then run the same FTBS upwind simulation with the new initial condition.

u_final = linear_convection_upwind(
    u0,
    c,
    dx,
    dt,
    nt,
    left_bc=left_boundary
)

snapshot_times = [0.0, 0.1, 0.2, 0.3, 0.4, 0.5]

plt.figure(figsize=(9, 5))

for t_snap in snapshot_times:
    nt_snap = int(t_snap / dt)

    u_snap = linear_convection_upwind(
        u0,
        c,
        dx,
        dt,
        nt_snap,
        left_bc=left_boundary
    )

    plt.plot(x, u_snap, label=f"t = {nt_snap * dt:.2f}")

plt.xlabel("x")
plt.ylabel("u(x,t)")
plt.title("Evolution of the semicircular tracer cloud")
plt.grid(True)
plt.legend()
plt.show()

Evolution of the semicircular tracer cloud over time.

The semicircular profile is also transported toward the right. Since it is smoother than the rectangular pulse, the effect of numerical diffusion is less abrupt, but the profile still becomes slightly smoother during the simulation.

Achtung!

Tracer transport in a 1D channel