Numerical Solution To Linear Convection Equation

Aim - Write a matlab code to solve linear convection equation and compare the plots for different number of grid points.

One dimensional linear convection equation is given by -

`(partial u)/(partial t) + c * (partial u)/(partial x) = 0`

First order forward difference for time derivative and first order backward difference for the spatial term has been used for discritizing the partial differential equation. Explicit approach was used. The difference equation is given by -

`(u_i ^(n+1) - u_i ^n )/(Delta t) + c * (u_i ^n - u_(i-1) ^ n)/ (Delta x) = 0`

Since we have considered explicit approach, for a given value of `Delta x` , we are not free to choose any random value for `Delta t` i.e. `Delta t` must be within a certain limit given by stability criteria. To know the value of `Delta t` upto which the solution is stable can be find out by performing von neumann stability analysis. For the solution to be stable `|(epsilon_i ^(n+1))/(epsilon_i ^n)| <= 1` i.e error should not grow up as we march in time. `epsilon_i ^(n+1)` represents error at node `i` at timestep `n+1` and `epsilon_i ^n` represents error at node `i` at timestep `n`. when we apply the above condition on the difference equation we find that `0 < (c.(Delta t)/(Delta x)) <= 1` has to satisfied for stable solution.

For the given problem -

constant velocity c = 1 m/s.

time step `Delta t` = 0.01 sec.

Domain length L = 1 m.

Matlab code for grid points n = 80 is given below.

clear all
close all
clc

L = 1;  % domain length
n = 80; % number of grid points
dx = L/(n-1); 
c = 1; % constant velocity
dt = 0.01; % time step
tEnd = 0.4;
nt = tEnd/dt;
x = linspace(0,L,n);
x_start = 0.1; % starting point of step
x_end = 0.3; % ending point of step
n_start = round(x_start/dx); % finding approximate index to start step 
n_end = round(x_end/dx); % finding approximate index to end step 
u = ones(1,n);
u(n_start:n_end) = 2;

%BC's 
u(1) = 1;

uold = u;
u_initial = u;

% time loop
for t = 1:nt
    clf
    %space loop
    for i = 2:n
        u(i) = uold(i)-(c*dt/dx)*(uold(i)-uold(i-1));
    end
    uold = u;
    plot(x,u,'LineWidth',1.2)
    hold on 
    title(['Courant no : ',num2str(c*dt/dx)])
    xlabel(['Time  = ',num2str(dt*t)])
    axis([0 L 1 2])
    ylabel('Velocity')
    drawnow
end
% plotting initial velocity        
plot(x,u_initial,'LineWidth',1.2)
legend('final velocity','initial velocity')

Plots for different grid points are given below.

Above plot is for n= 20. Maximum final velocity was found to be u = 1.6834 m/s.

Above plot is for n= 40. Maximum final velocity was found to be u = 1.8525 m/s.

Above plot is for n= 80. Maximum final velocity was found to be u = 1.9983 m/s.

Above plot is for n= 160. Since the courant number has exceeded its max value i.e. 1 for the given problem, the solution became unstable.

Conclusion -

1. It has been found that as we reduce the number of grid points the velocity profile starts diffusing.
2. Solution becomes unstable for courant numbers > 1 (note for the above problem).