Numerical Solution of 1-D Linear Convection Equation using Octave

AIM:-

• To write code to solve the 1D linear convection equation (a PDE) following time marching procedure.
• To compare the final velocity profiles with number of grid points.

THEORY AND ASSUMPTIONS:-

1D Linear Convection Equation is `(del u)/ (del t) + C*(del u)/(del x) = 0`.

1. Domain length is L = 1 m.
2. The initial velocity profile is a step function. It is equal to 2 m/s between x = 0.1 and 0.3 and 1 m/s everywhere else.
3. First order forward differencing scheme is used for the time derivative.
4. First order backward/rearward differencing scheme is used for the space term.
5. Time step and number of grid points are variables.
6. C = 1 (Constant Velocity Term).
7. Time step was set as 0.01 s.
8. Time = 0.4 s.
9. Nodes, n = 20, 40, 80 and 160 are used for comparison of the velocity profiles.

PROGRAM CODE IN OCTAVE:-

clear all
close all
clc

% Linear Convection Equation is do_u/do_t + C*do_u/do_x = 0

% Inputs
% Domain Length = 1 m
L = 1;
% Constant Velocity Term = 1 m/s
C = 1;
% Time step = 0.01 s
dt = 0.01;
% Number of Grid points
n = [20 40 80 160];
length_n = length(n);

% Loop of Grid Size
for k = 1:length_n
    
  % Generating grid points array x
  x = linspace(0,L,n(k));
  % Grid size
  dx = x(2) - x(1);

  % Initial Velocity Profile be a Step Function
  % Everywhere it is 1 m/s except from x = 0.1 to x = 0.3 (2 m/s)
  u_step_1 = 1;
  u_step_2 = 2;
  u = u_step_1*ones(1,n(k));
  x_start = 0.1;
  x_end = 0.3;
  
  % Calculating square wave start and end nodes
  n_start = 1 + round(n(k)*x_start/L);
  n_end = 1 + round(n(k)*x_end/L);
  % Assigning step 2 velocity = 2 m/s
  u(n_start:n_end) = u_step_2;
  
  % Plotting initial velocity profile
  figure(1)
  plot(x,u);
  title_text = sprintf('Linear Convection - Initial and Final Velocity Profiles for n = %d', n(k));
  title(title_text)
  hold on;
  grid on;
  axis([0 1 .5 2.5]);
  xlabel('X Axis (m)')
  ylabel('Velocity (m/s)')

  % Boundary Condition
  u(1) = 1;
  u_old = u;

  % Time = 0.4 s
  t = 0.4;
  % Number of time step
  nt = t/dt;

    % Loop of Time Step
    for i = 1:nt
      % Loop of Nodes
      for j = 2:n(k)
        u(j) = u_old(j) - (C*dt/dx) * (u_old(j) - u_old(j-1));
      endfor
      u_old = u;
    endfor
    
  % Plotting Final Velocity Profile
  plot(x,u,'r');
  legend('Initial Velocity Profile','Final Velocity Profile');
  pause(.3);
  
  % Saving the Plot
  file_text = sprintf('plot_n_%d.png',n(k));
  saveas (gca, file_text);
  delete(figure(1));
  
endfor

PROGRAM CODE EXPLANATION:-

1. Domain length, 1 m is assigned to the variable L.
2. Constant Velocity Term, C is declared and a value of 1 m/s was assigned to it.
3. Time step, 0.02 s is assigned to the variable dt.
4. Number of grid points is defined as an array of 4 elements (20, 40, 80 and 180) and assigned to the variable n.
5. The number of elements of array n is calculated by using the length() function and assigned to the variable length_n.
6. Since we compare the velocity profiles with different number of grid points, a for loop is used for the same. The for loop will execute length_n times (here, 4 times).
7. The grid points array is defined as x and is allocated values by using linspace() function. The starting value is 0 and ending value will be domain length, L.
8. The grid size is uniform as we have used linspace() function. Grid size is equal to the distance between any adjacent nodes. This value is assigned to the variable dx.
9. The initial velocity profile is a step function. As mentioned in the Assumption section, it has two values viz 1 m/s and 2 m/s. These values are assigned to the variables u_step_1 and u_step_2 respectively.
10. Since it is given that the lower value is eveywhere except from 0.1 to 0.3, the velocity variable u is declared and initialized with the lower values, ie u_step_1. The values 0.1 and 0.3 are assigned to the variables x_start and x_end respectively.
11. The higher value, u_step_1 (2 m/s) will start from x_start (0.1) and will end at x_end (0.3). We need to find the corresponding node numbers. The same were calculated by using the equation '1 + round(grid points * x value / Domain Length)'. The function round() is used to get an integer value. The two node numbers were assigned to the variables n_start and n_end respectively.
12. From x = x_start to x = x_end also means that from n = n_start to n = n_end. So the step 2 value (u_step_2 = 2 m/s) is assigned to the velocity variable during this range.
13. Initial Velocity Profile is plotted and suitable title, labels and legends were provided.
14. Boundary condition of velocity at node 1 = 1 m/s is provided separately though it is already given.
15. A variable u_old is declared to store the velocity values at present.
16. Since we follow time marching method, the time variable t is defined and assigned a value of 0.4 s.
17. Since time step is known (dt), the number of time steps is calculated and assinged to the variable nt.
18. Since there are nt time steps, a for loop is used to track it. This for loop will execute nt time.
19. The velocity at each node is calculated for each time step. This is tracked by another for loop inside time step loop.The inside loop will start from 2nd node to the last node. After finishing the inner loop, before going to the next time step, the u_old is renewed as that time step is over.
20. After executing the time step loop completely, the final velocity profile is also plotted in the same plot of initial velocity profile for comparison.
21. The graph is saved with suitable file name.
22. pause() function is used to see the graph once before plotting next one. delete() function is used to delete the existing graph and will help to generate next graph independently.

INITIAL AND FINAL VELOCITY PROFILES:-

• For number of grid points, n = 20:-

• For number of grid points, n = 40:

• For number of grid points, n = 80:

• For number of grid points, n = 160:

CONCLUSIONS:-

• When the grid size is large (number of grid points are less), the diffusion observed is more. The peek velocity reduction is also more.
• As the grid size is getting reduced (number of grid points are increased), the diffusion observed became lesser and lesser. The peek velocity reduction was also becoming lesser. When n = 80, the peek velocity was the same as of initial velocity profile.
• It was also observed that the appearance of initial velocity profile, ie step function was more accurate when small grid size is used.
• In all the cases, convection nature was clearly visible except for n = 160.
• It is true that we will get accurate result when number of grid points is increased, however, when a limit is crossed, there would be some undesirable effects as evident from the last plot (n = 160). In the last plot, we did not get an expected graph nor convection nature was visible. This is due to instability. Hence, choosing number of grid points has atmost importance. If it is very less, the accuracy will be less, at the same time, if it is very large, it will cause inaccurate/undesirable/unexpected results due to instability.