Linear Convection (Different nodes)

Theory :

To solve a 1D Linear Convection problem with the following assumptions.

• Length of domain = 1 Metre.
• 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 (C) = 1.
• Time Duration (t) = 0.4 seconds
• Time Interval (dt) = 0.01 second.
• Compare output for nodes n=20,40,80,160,320.
• Use FDS for time derivative and BDS for space term
• COmpute the exact node points for change in the velocity of our step function.

***********************Start of Program*******************************

Main Program :

% Program to Compute 1D Linear convection
clear all
close all
clc

% defining inputs for the convection problem

L=1; % Length = 1 metre

c=1; %Constant
n=[20,40,80,160,320]; % Array defining no. of grid points
t=0.4; % time interval = 0.4 seconds
dt=0.01; % time step (dt) = 0.01 second.

%defining the starting and ending point.
x_start=0.1;
x_end=0.3;

nt=t/dt; %defining no. of time steps as nt

for i=1:length(n) %Primary Loop
x=linspace(0,L,n(i)); %% Creating our equidistant mesh
dx= L/(n(i)-1); % computing mesh size dx

u=ones(1,n(i)); %Initialising array 'u' with vale of 1


% Loops to calculate the start and end position of the peak velocity.

for m=1:n(i) % First Loop computes the absolute error between actual values of the array and our required value of start and endpoint.

diff1(m)=abs(x(m)-0.1); %using absolute value to store only magnitudes
diff2(m)=abs(x(m)-0.3);
end
min1=min(diff1); %takes the minimum value from 1st array to calculate start node point
min2=min(diff2); %takes the minimum value from 2nd array to calculate end node point

for m=1:n(i) %Seconf for loop here checks if the duifference is equal to start value and end value, if yes, it assigns that particular value of node as our start and end points.
if min1==diff1(m)
n_start=m;
end
if min2==diff2(m)
n_stop=m;
end
end
u(1)=1; % Initialising 1st node value of array as 1
u(n_start:n_stop)=2; % using ange operator, initialising the value between start and end point as 2.

uold=u; % Copying our original array
u_initial=u; % Copying our original array

for j=1:nt %For loop to compute our time marching solution for each node respectively and displaying the output for the same
for k=2:n(i)
u(k)=uold(k)-(c*(dt/dx)*(uold(k)-uold(k-1)));
end
uold=u;% Copying our array to the old array

figure(i)
grid on
plot(x,u_initial,'b')
hold on
plot(x,u,'r')
xlabel('length of domain')
ylabel('velocity')
title(['Velocity profile for node n = ',sprintf('%d',n(i))]) %Printing the title of the plot to correspond to the no. of nodes it belongs to.
axis([0 1.5 0 2.5]) % Limiting our axis values
end
end

****************************End of Program***************************

Outputs:

• For node n=20;

• For node n=40;

• For node n=80;

• For n=160;

• For n=320;

Explanation and Conclusion:

1. Figure 1 and 2 for n =20 and n=40 respectively shows that the peak amplitude of the function keeps on decreasing as the solution marches forward in time. This implies that fro nodes n=20 and n=40, the peak amplitude is inversely proportional to the length of the domain.
2. Figure 3 for n=80, shows the most accurate graph for our time marching solution. The peak amplitude remains intact for the duration of the plot and the solution seems to be stable. This depicts the optimum result for our function is achieved when n=80.
3. Figure and 4 and 5 for n=160 and n=320 respectively shows the unstable solution. Higher the number of nodes implies that as the diffusion rate is increased for our solution, the time marching solution gets blown up and gives non-reliable results with unwanted peaks all over the place.
4. Thus, the grid spacing has to be optimized for best results as if the spacing is less, the solution blows up and if the spacing is more, the amplitude becomes inversely proportional to the length of the domain and eventually will die off giving non-accurate results.

Key Challenge for the solution was to accurately determine the node at which the velocity shoots up from 1m/s to 2m/s and then comes back from 2m/s to 1m/s.

The validity of the program was dependent on determining the position of these nodes which was successfully developed by introducing the concept of error and minimizing it to find the exact position. For this, extra 2 loops had to be introduced within the program, but it gives the most accurate and verified results for the program.