How to Solve the 1D Advection-Diffusion Equation in Python?

Estimated read time 2 min read

The 1D Advection-Diffusion Equation is a partial differential equation that describes the transport of a quantity (e.g. heat, mass, or momentum) in a fluid. Here is an example of how to solve this equation in Python using the finite difference method:

import numpy as np
import matplotlib.pyplot as plt

# define the parameters of the problem
L = 1.0 # length of the domain
nx = 100 # number of grid points
dx = L / nx # grid spacing
x = np.linspace(0.0, L, nx+1) # grid points
D = 1.0 # diffusion coefficient
u = 1.0 # advection velocity
dt = dx**2 / (4*D) # time step
t = np.arange(0.0, 0.2, dt) # time points

# initialize the solution array
C = np.zeros((len(t), len(x)))
C[0] = np.sin(np.pi*x/L) # initial condition

# solve the equation using the finite difference method
for n in range(1, len(t)):
    # calculate the diffusion term using central differences
    dCdx2 = np.zeros(len(x))
    dCdx2[1:-1] = (C[n-1,:-2] - 2*C[n-1,1:-1] + C[n-1,2:]) / dx**2
    # calculate the advection term using upwind differences
    dCdx = np.zeros(len(x))
    dCdx[1:] = (C[n-1,1:] - C[n-1,:-1]) / dx
    dCdx[0] = (C[n-1,0] - C[n-1,-1]) / dx # periodic boundary condition
    # update the solution using the forward Euler method
    C[n,:] = C[n-1,:] + dt*(-u*dCdx + D*dCdx2)

# plot the solution
fig, ax = plt.subplots()
for i in range(0, len(t), len(t)//10):
    ax.plot(x, C[i], label='t=%.2f' % t[i])
ax.set_xlabel('x')
ax.set_ylabel('C')
ax.legend()
plt.show()

In this solution, we first define the parameters of the problem, including the length of the domain L, the number of grid points nx, the grid spacing dx, the diffusion coefficient D, the advection velocity u, the time step dt, and the time points t.

We then initialize the solution array C with the initial condition, which is a sine wave.

We then solve the equation using the finite difference method. In each time step, we first calculate the diffusion term using central differences and the advection term using upwind differences. We then update the solution using the forward Euler method.

Finally, we plot the solution at 10 equally spaced time points using the plot function from Matplotlib.

To use this code with different parameters, simply modify the values of the parameters at the beginning of the script.

You May Also Like

More From Author

+ There are no comments

Add yours

Leave a Reply