1 D linear convection

Objective:

Write the program to solve 1 D linear convection problem.

Given:

• Assume that the domain length is L = 1m.
• The initial velocity profile is a step function. It is equal to 2m/s between x= 0.1 and 0.3 and 1m/s everywhere else.
• Constant Velocity, C = 1.
• Set time step = 0.01.
• time = 0.4seconds.
• Make the comparison for n=20,40,80 and 160.

Introduction:

1 D wave equation is given by

∂u/∂t+c∂u/∂x=0

where,

∂u/∂t = time derivative

∂u/∂x = convective terms

Code:

clear all
close all
clc

% define domain length and convection coefficient
L = 1;                        % Length of domain
c= 1;                         % constant Velocity
n = [20,40,80,160];      % no. of grid points

% simulation time
t = 0.4;

% define time step
dt = 0.01;

% define the nodes where velocity is 2m/s
x_start = 0.1;
x_end = 0.3;

% no. of time steps
Nt = t/dt;

% for loop is used to obtain velocity profile for each value of n
for i = 1:length(n)
  
  % discretizing length 'x' w.r.t 'n'
  x = linspace(0,L,n(i));
  dx = L/(n(i)-1);  % grid size
  
  % create an array of nodes at which velocity is 1 m/s
  u = ones(1,n(i));
  
  % evaluating start and end of nodes for different values of 'i'
  n_start(i) = round((x_start/dx)+1);
  n_end(i) = round((x_end/dx)+1);
  
  % define the nodes where velocity iis 2 m/s
  u(n_start(i):n_end(i)) = 2;
  uold = u;
  u_initial = u;
  
  % time marching solution
  for p= 1:Nt         % time loop du/dt
    for j = 2:n(i)      % Space loop  (c*dt/dx)
      u(j) = uold(j) - (c*dt/dx)*(uold(j)-uold(j-1));
    end
    
    % update the old velocity
    uold = u;
    
    % plot
    figure(i)
    plot(x,u,'r')
    hold on
    plot(x,u_initial,'b')

    grid on
    
    % Notation
    xlabel('length of domain')
    ylabel('Velocity')
    title({'1 D linear convection',[ 'No. of grid point =',num2str(n(i))]},'fontsize',16)
    legend('final velocity profile','initial velocity profile')
    axis([0 1  0 2.5])
    pause(dt)
  end
  
end

Output:

1.For grid points = 20

Fig 1 : No. of grid points = 20

2.For grid points = 40

Fig 2 : No. of grid points = 40

3.For grid points = 80

Fig 3 : No. of grid points = 80

4.For grid points = 160

Fig 4 : No. of grid points = 160

Expalnation of code:

• 'for i = 1:length(n)' is used to extract all 4 values in n matrix.
• 'ones' command is used to create unit matrix of 1 coloumn and n rows.
• 'round' command is used to round off the nearby decimal value.

Conclusion:

1.For grid points(n) = 20, initial velocity profile is not square in shape it looks like trapezoidal and as it travel rightwards amplitude is continuosly decreasing and final shape of velocity profile looks like sine wave and length of wave is increased.

2.For grid points(n) = 40, initial velocity profile is not perfect square in shape and amplitude of wave is gradualy decreasing.

3.For grid points(n) = 80, shape of wave is not distorting as it move rightwards and amplitude of wave remains constant hence, it can be said that wave is moving with constant velocity.

4.For grid points(n) = 160, initial velocity profile is perfect square but as wave moves rightward ,it is oscilating continuosly about line of constant velocity of 1m/s thus creating sharp leading and trailing edges of waveforms.

Distortions in the waveforms are occured due to CFL number.

5.From above 4 graphs, it can be said that as the no. of grid points increases waves travels with decreasing amplitude first and then it travels with constant velocity and after optimum no. of grid points extra no. of grid points causes distortion in waveforms due to ignored diffusion terms during numerical calculations.

References:

1. How to solve wave equation in matlab [https://www.youtube.com/watch?v=7xdl8gxq2jM]