1D Linear Convection

Objective:

Solve 1D Linear Convection equation using Matlab.

1) Based on varying number of grid points

The equation is given by

A time marching procedure is used to numerically solve the equation by employing the Backward-differencing method for the space terms.

Numerical Solution

Assumptions:

a) The domain length is L=1m.

b) The intial velocity profile is a step function. It is equal to 2m/s between x=0.1 and x=0.3 and 1m/s everywhere else.

c) First order forward differencing for the time derivative. First order backward differencing for space terms.

d) C = 1 and time step = 0.01.Total time = 0.4 seconds.

The comparison is based on n=20,40,80 and 160. Hence, the accuracy depends on dx (grid spacing).

The above graph represents the graph of the initial and final velocity profiles when the number of grid points is 21. The accuracy of the final velocity profile depends on the number of grid points since 'dt' is constant. We see that the profile loses its square wave type shape and does not reach the peak velocity. This is due to the artificial diffusion that takes place over time which is nothing but the error term.

This graph represents the velocity profile when number of grid points is increased to 41. Here, there is an improvement in the shape of the profile. This is due to the increase in number of grid points. Also the peak velocity is closer to the actual value.

For the number of grid points equal to 81, the profile is more accurate. The peak velocity is very close to the actual value. But there is a loss of shape as it moves downstream. The aaccuracy of the solution is increased and the diffusion is reduced. The solution is stable.

From this graph, we see that further increase in the number of grid points causes the solution to blow up. The soultion is not stable and it produces oscillations even though the diffusion is reduced. The peak velocity is also not visible.

Snippet of the program:

x_start = 0.1; % length from which velocity becomes 2 m/s
x_stop = 0.3; % upto which velocity is 2 m/s

n_start = round((x_start/dx)+1); 
n_stop = round((x_stop/dx)+1);

u(n_start:n_stop) = 2; % square wave

u_old = u;
u_initial = u;

% Time marching

  % 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))); % Backward differencing
             
         end 
         
         % update velocity 
         u_old = u;
         plot(x,u_initial,'r');
         hold on;
         plot(x,u,'b');
         axis([0 1 0.5 2.5])
         xlabel('Length, m');
         ylabel('Velocity, m/s');
         pause(0.1);
         grid on;
         
     end 

2) Based on varying time steps

Assumptions:

a) Grid points n = 80

b) Time steps used are 1e-4,1e-3,1e-2 and 1e-1

c) Numerical solution compared at t = 0.4 seconds.

The above graph represents the velocity profiles at different time steps. We see for the value of 'dt = 1e-1', the solution blows up and this due to the unstability. For a value of 'dt = 1e-2', it is very close to the initial profile, but overtime it loses shape. The profiles for '1e-3' and '1e-4' time steps are similar. The peak velocity is not reached.

The computaional time for various time steps is represented in a bar graph. We see that the smaller the time step used that is 'dt = 1e-4', longer the computational time. Thus an optimum time step must be chosen to get the required efficiency without compromising on the accuracy.

Snippet of the program:

function Onedimensional_linear_convection_time_step(dt,colour)

L = 1; % Length = 1 metre
C = 1; % coefficient
n = 81; % number of grid points
t = 0.4; % total time
nt = t/dt; % number of time steps

% grid or mesh
x = linspace(0,L,n); % creates a mesh of equal spacing
dx = x(2) - x(1); % grid spacing

% Boundary conditions
u = ones(1,n);

x_start = 0.1; % length from which velocity becomes 2 m/s
x_stop = 0.3; % upto which velocity is 2 m/s

n_start = round((x_start/dx)+1); 
n_stop = round((x_stop/dx)+1);

u(n_start:n_stop) = 2; % square wave

u_old = u;

% Time marching

  % 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))); % Backwarddifferencing
             
         end 
         
         % update velocity 
         u_old = u;
         
     end 
     
     plot(x,u,colour);