Numerical simulation of a velocity wave in a 1D channel with nonlinear convection

Papadopol Lucian Ioan
l.i.papadopol@gmail.com

1. Goal of the simulation

In this exercise we simulate a phenomenon of nonlinear convection in a one-dimensional straight channel.

In the previous notebook we studied the transport of a tracer cloud in a 1D channel. In that case the current velocity was constant and equal to c.

The model was:

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

This is the one-dimensional linear convection equation.

In that model all parts of the cloud were transported with the same velocity. Ideally, the initial profile moved along the channel without changing its shape.

Now we study a more interesting case: a velocity perturbation in the current that transports itself.

In this new model the transport velocity is no longer an external constant. It depends on the value of the unknown function itself.

The model becomes:

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

This is the nonlinear convection equation.

The goal of the simulation is to observe how an initial perturbation moves and deforms when the higher parts travel faster than the lower parts.

2. Physical system being simulated

We consider a straight channel in which a current perturbation propagates.

The quantity

u(x,t)

represents a normalized local fluid velocity at position x and time t.

We can imagine u(x,t) as a small velocity wave moving along the channel. Where u is larger, the fluid moves faster. Where u is smaller, the fluid moves more slowly.

This is the fundamental difference from the linear case.

In the linear case we had:

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

where c was a constant velocity.

In the nonlinear case we instead have:

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

so the transport velocity is the local value of u itself.

3. Difference between linear and nonlinear convection

In the linear case, the transport velocity is constant. Therefore each point of the profile moves with the same speed. If there are no numerical errors or diffusion, the initial profile is simply translated.

In the nonlinear case, the transport velocity depends on the value of the solution. Each point of the profile moves with a different speed.

linear convection:
all parts of the profile travel with the same velocity c

nonlinear convection:
each point travels with a velocity equal to the local value u

The solution therefore does not simply translate. It changes shape over time.

4. Physical interpretation

Nonlinear convection produces a behavior that is very different from linear convection.

In linear convection the initial profile moves almost rigidly. In nonlinear convection, however, every point of the profile moves with a velocity that depends on its own value.

If part of the profile has a higher value, it tends to move faster. If part of the profile has a lower value, it tends to move more slowly.

As a consequence, the higher parts of the wave can catch up with the lower parts. This produces a progressive deformation of the profile. A typical result is the formation of an increasingly steep front.

high values of u -> travel faster
low values of u  -> travel more slowly

Thus the initial profile is not only transported; it is also deformed.

5. Mathematical model

The mathematical model is:

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

where:

• u(x,t) is the normalized local velocity;
• x is the position along the channel;
• t is time.

The term

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

describes the time variation of the local velocity.

The term

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

describes nonlinear transport.

The equation is nonlinear because the unknown function u multiplies its own spatial derivative.

In the linear case the transport coefficient was a constant:

c

In the nonlinear case the transport coefficient is:

u

This is the main reason why the behavior of the solution changes.

6. Connection with the linear case

In the linear case the numerical scheme contained the term:

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

where c was constant.

In the nonlinear case the velocity is no longer constant. It changes from point to point. Therefore the local coefficient becomes:

u_i^n\frac{\Delta t}{\Delta x}

The central change is therefore:

c \quad \longrightarrow \quad u_i^n

In the linear case all points are transported with the same velocity. In the nonlinear case each point is transported with a velocity that depends on the local value of the solution.

7. Spatial domain

We consider a channel of length L.

The spatial domain is:

0 \le x \le L

We divide the domain into a grid of points:

x_i = i\Delta x

where:

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

The spatial step is:

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

The numerical value

u_i^n

represents the approximation of:

u(x_i,t_n)

that is, the normalized velocity at grid point x_i and time t_n.

8. Time discretization

Time is also discretized. We write:

t_n = n\Delta t

where \Delta t is the time step.

The goal of the numerical method is to compute the profile at the next time level:

u^{n+1}

from the current profile:

u^n

As in the previous notebooks, the simulation advances through successive time steps.

9. Upwind scheme for nonlinear convection

We start from the equation:

\frac{\partial u}{\partial t} + u\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}

Assuming that u is positive, the information moves to the right. Therefore, as in the linear case with c>0, 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 nonlinear convection equation gives:

\frac{u_i^{n+1}-u_i^n}{\Delta t}+u_i^n\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-u_i^n\frac{\Delta t}{\Delta x}\left(u_i^n-u_{i-1}^n\right)

This is the explicit upwind scheme for nonlinear convection in the case u>0.

10. Comparison between the numerical schemes

In the linear case we had:

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

In the nonlinear case we have:

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

The structure of the scheme is very similar. The difference is that the linear case contains the constant velocity c, while the nonlinear case contains the local value u_i^n.

c \quad \text{is replaced by} \quad u_i^n

This small modification changes the behavior of the solution deeply.

11. Courant number in the nonlinear case

In the linear case the Courant number was:

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

In the nonlinear case we no longer have a constant velocity c. The local velocity is given by:

u(x,t)

For this reason the Courant number must be computed using the maximum velocity in the domain:

C = \frac{\max |u|\Delta t}{\Delta x}

The CFL condition becomes:

\frac{\max |u|\Delta t}{\Delta x} \le 1

This means that, during a single time step, even the fastest part of the wave must not cross more than one grid cell.

In the code we use a reusable function:

def cfl_nonlinear_convection(u, dt, dx):
    """Courant number for nonlinear convection."""
    return np.max(np.abs(u)) * dt / dx

12. Choice of the time step

In the linear case we could choose a desired Courant number:

C_{\text{target}}

and obtain the time step as:

\Delta t = C_{\text{target}}\frac{\Delta x}{|c|}

In the nonlinear case we use the same idea, but instead of |c| we use the maximum initial value of |u|:

\Delta t = C_{\text{target}}\frac{\Delta x}{\max |u|}

This choice ensures that the CFL condition is satisfied at least at the beginning of the simulation.

Since the profile can change during the simulation, it is useful to recompute the effective Courant number and check that it remains less than or equal to one.

13. Initial condition

To observe the nonlinear behavior clearly, we use a localized initial perturbation.

In the linear case, an initial profile is mainly transported. In the nonlinear case, the profile is transported and deformed.

The parts with higher value travel faster than the parts with lower value. Therefore we expect the wave front to become progressively steeper.

14. First step — defining the physical system

import numpy as np
import matplotlib.pyplot as plt

def cfl_nonlinear_convection(u, dt, dx):
    """Courant number for nonlinear convection."""
    return np.max(np.abs(u)) * 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]

# Desired Courant number
C_target = 0.8

# Initial condition: velocity wave in the channel
u0 = np.ones_like(x)

x_c = 0.225
R = 0.075

s = (x - x_c) / R
mask = np.abs(s) <= 1.0

u0[mask] = 1.0 + np.sqrt(1.0 - s[mask]**2)

# Only after defining u0 do we compute dt
dt = C_target * dx / np.max(np.abs(u0))

C = cfl_nonlinear_convection(u0, dt, dx)

print("dx =", dx)
print("dt =", dt)
print("Courant number =", C)
print("max u0 =", np.max(np.abs(u0)))

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

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

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

The numerical parameters are:

dx = 0.005
dt = 0.002
Courant number = 0.8
max u0 = 2.0

15. Second step — creating the numerical simulation method

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

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

        u[1:] = un[1:] - un[1:] * dt / dx * (un[1:] - un[:-1])

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

    return u

The update formula implemented by the function is:

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

16. Third step — initial conditions and simulation time

def left_boundary(u, dx):
    return 1.0

t_end = 0.5
nt = int(t_end / dt)

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

With the chosen values, the number of time steps is:

nt = 250
simulated final time = 0.5

17. Fourth step — running the simulation

u_final = nonlinear_convection_upwind(
    u0,
    dx,
    dt,
    nt,
    left_bc=left_boundary
)

18. Plotting the result

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

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

plt.xlabel("x")
plt.ylabel("u(x,t)")
plt.title("Nonlinear convection of a velocity perturbation")
plt.grid(True)
plt.legend()
plt.show()

19. Interpretation of the result

The result shows that, unlike linear convection, the initial profile is not simply transported to the right. Since the local propagation velocity coincides with the value of u, the higher parts of the perturbation move faster than the lower parts.

This causes a progressive deformation of the profile and makes the right front increasingly steep. The upwind scheme keeps the simulation stable, but it also introduces numerical diffusion, visible in the smoothing of the profile and in the reduction of the maximum value.

This behavior is important because it anticipates the study of more advanced models such as Burgers’ equation, shock formation, conservation laws, finite volume methods, and the Lax-Friedrichs and Lax-Wendroff schemes.