Solving 1D Convection equation using Matlab

In this project, the 1d convection equation was solved and data was plotted comparing the velocities at different number of grid points.

The 1D convection equation is given below.

`(delu)/(delt)+c(delu)/(delx)=0`

Discretizing by forward difference for time derivative and backward difference for space derivative, we get

`((u^(n+1))_i-u_i)/(Deltat)+c(u_j-u_(j-1))/(Deltax) = 0`

Below is the boundary condition used at t=0

`u = {(2 if0.1<=x<=0.3),(1ifx<0.1 andx>0.3):}`

It was found that there was some artificial diffusion (numerical diffusion) involved even though the equation defining the problem only had convection term. This is termed as Courant–Friedrichs–Lewy (CFL) condition . The numerical diffusion (ND) depends on the value of C, `Deltax` and `Deltat`.

`ND= Cdx/2(1-Cdt/dx)`

The solution is stable only when Numerical diffusion is positive. This depends on CFL number.

`CFL = Cdt/dx`

The solution is stable only if this is less than one.

We can see that the numerical diffusion gets minimised as Courant number approaches one. But for Courant number greater than one, it becomes unstable. This can be encountered by chosing timestep and grid spacing wisely.

clear all
close all
clc

%Inputs for the Linear convection problem

L = 1;
n = 161; %Taking 21 instead of 20, because ease of calculation of dx
c = 1;
dt = 0.01;

%Calculating Grid spacing
x = linspace(0,L,n);
dx = L/(n-1);

%Defining start and end points for square wave

%x_start = 0.1;
n_start = (0.1/dx)+1; %+1 for intial term zero
%x_stop = 0.3;
n_stop = (0.3/dx)+1;

%Initialize array and apply boundary condition
u = ones(1,n);
u(n_start:n_stop) = 2;

%Copying array
u_old = u;
u_initial = u;

nt = 500; %number of time steps

%Time Loop

for k = 1:nt

    %Space Loop

    for i = 2:n
      u(i) = u_old(i) - (c*dt/dx)*(u_old(i)-u_old(i-1));
     end
      
     %update old velocities
      u_old = u;
      
      %Plotting the velocities at t = 0.4s
      nt_plot = 0.4/dt;
      
      if(k == nt_plot)
          figure(1)
          plot(x,u,'b',x,u_initial,'r')
          legend('t=0.4s','t=os')
          str = sprintf("Number of Grid points = %d",n); 
          title(str);
          axis([0 1 0 2])
          xlabel('X axis')
          ylabel('Velocity')
       endif

  
  
endfor