Linear Convection - Solving 1D linear convection equation for different values of total number of grid points (n)

Objective:-

1. Assume that the domain length is L = 1m

2. 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

3. Use first order forward differencing for the time derivative

4. Use first order backward/rearward differencing for the space term.

5. Time step and number of grid points should be variables in your calculation

6. C = 1

7. Set time step = 0.01

8. For C = 1 and for time = 0.4seconds compare the original and final velocity profiles.

9. Make the comparison for n=20,40,80 and 160.

Software used:- MATLAB

Solution:-

Matlab code for solving 1D linear convection equation

clear all
close all
clc

L = 1 % Domain lenght in meter
n = 20 % no. of grid points
c = 1 %linear convection velocity
x = linspace(0,L,n) % mesh
dx = x(2) - x(1) 
dt = 0.01 % time step
time = 0.4 % total time in seconds
nt = time/dt % total no. of time steps

%start point for square wave
x_start = 0.1 
n_start = round((x_start/dx) + 1) % finding the integer value of the node for x_start

%stop point for square wave
x_stop = 0.3
n_stop = round((x_stop/dx) + 1) % finding the integer value of the node for x_stop

u = ones(L,n)
u(:,n_start:n_stop) = 2
u_initial = u
uold = u
ct = 1
%time loop
for k = 1:nt
    % space loop
   for i = 2:n
       u(i) = uold(i) - (c*(dt/dx))*(uold(i) - uold(i-1)); % solving linear convection equation
   end
   %update old velocities
   uold = u

%plotting the velocities
   plot(x,u_initial,'r','linewidth',1) % original velocity profile
   hold on
   plot(x,u,'b') % final velocity profile
   xlabel('x aixs')
   ylabel('velocity')
   pause(0.03)
   
   M(ct) = getframe(gcf); % function used to capture all the generate frames in the loop
   ct = ct+1;

end
% Creating movie to visualize convection
movie(M) 
convection_video = VideoWriter('convection_n=20.avi', 'Uncompressed AVI');
open(convection_video)
writeVideo(convection_video,M)
close(convection_video)

Report:

• First we define all the input data, like Domain length (D), no. of grid points (n), linear convection velocity (c), mesh, dx, dt, no of time steps(nt).
• Then we find the node at which x=0.1 and x=0.3. We have to find this because we have to plot a square wave form x=0.1 to x=0.3
• then we define the square wave for velocity "u"
• after that we copy the value of "x" in "u_initial", for plotting the original square wave plot. And then we define "uold=u" for solving the linear convection equation inside the for loop.
• We discretize the linear convection equation to find the numerical solution.
• Now we use "for loop" for finding U(i) from 2 to "n" number of grid points.
• this for loop is called space loop and this is defined inside the time loop for finding the time marching solution from i=1 to total number of time steps (nt).
• Inside the time loop we also define use "getframe(gcf)" command to capture all the frame generated by the programe.
• These frame are stiched together to form a video of linear convection.

The above programe was executed at four different values of number of grid points, n=20, 40,80,160. This generated four graphs for linear convection and four videos were also generated.

• The above graph is generated at n=20. The velocity profile is moving from left to right.
• The red colour profile represents the original velocity profile and the blue colour profile is the final velocity profile
• As the profile moves from left to right, it gets smoother and the peak value of velocity is continuously reducing.

• The above graph is generated at n=40. The velocity profile is moving from left to right.
• The red colour profile represents the original velocity profile and the blue colour profile is the final velocity profile
• As compared to the first graph, this graph is less smoother and its peak velocity did not decreased drasticaly as it did in the first graph.

• The above graph is generated at n=80. The velocity profile is moving from left to right.
• The red colour profile represents the original velocity profile and the blue colour profile is the final velocity profile.
• In this graph we observe that the sharpness of the square wave is maintained for a longer time as compared to the first two graphs. The paek value of the velocity is almost maintained till the final profile.

• The above graph is generated at n=160. The velocity profile is moving from left to right.
• No square wave is generated in this graph. Here, we see a lot of wiggles and oscilations. The solution is ustable. 

Conclusion:-

• Four graphs were generated at different number of grid points. At first the velocity profile baceme more stable and more squre wave with the increase in number of grid points from n=20 to n=80. But when we generated the profile at n=160 the solution was ustable, square wave profile was completely lost. This means there is a limit on how much finner grid we can use to find the solution. In general the finner the grid the better the solution will be, but there is a limit for finner grid after which the solution will be ustable.

Animation of linear convection:-

for n=20

https://youtu.be/xMNOgHfbwSc

for n=40

https://youtu.be/nf0BN-LzDVk

for n=80

https://youtu.be/jKDtyNPPhVI

for n=160

https://youtu.be/3zX4D6KYRxY