Linear Convection for varying number of nodes and Linear Convection for varying timesteps

Aim: Our objective is to solve the PDE 1-dimensional linear convection equation through MatLab code and compare the effect grid size on the solution.

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

Step 1: Inputting the constant values provided in the problem

'L' domain length = 1 meter

t = 0.4 sec

c=1

dt=0.01

and as we need to find the solution for different n values

for this

n given with a provision to enter manually after running based on the given input n value the program executes

Step 2:

U array has been defined as U is 1 at all the X values excepts 0.1 to 0.3 band

So we created a function that predict the nodes at 0.1 and 0.3 and gives you the node number.

So using range operator we defined the value of U in between these nodes as 2

Step 3

The function has been carried out through for loops to move the step function spacely and timely. Where we used first-order forward differencing for time and backward differencing for space

Result

If you run the function It will ask for the number of nodes

at n=20

The solution is stable but the velocity curve goes on reducing this is because even we are solving for the linear convection equation we can see the wave is diffusing as it is moving.Here when we are comming up the first order the term we are ignoring is second order since the error is proportional to d2x/dx2 which is the diffusion term.this is called as artificial diffusion or numerical diffusion.

It can be reduced by using finite grids,since error is proportional to dx.

at n=40

The solution is stable but there is a reduction in the step size while it is still better than the previous.

at n=80

The solution is stable but there is also a minimal reduction in the step size.

at n=160

The solution got blown up because of the excessive number of nodes.

if you try with second order approximation instaed of first order you will the solution got blow up.The scheme to apply mainly depend upon in which direction the information is travelling that depends method of characteristics and the other is stability which is more important.

Grid spacing definitely enhances the quality of the solution to some extent but when it goes beyond solution gets to blow up. This phenomenon is based on courant number which is `c*((dt)/(dx))` if it is greater than 1 then the solution gets to blow up. In the case of n=160, the dx gets is very small which is resulting in higher Courant number more than 1. This is causing the solution blow up.

program :

clear all
close all
clc

%inputs for the linear convection problem

L=1; %1 meter length
n=input('number of nodes='); %the input should be required number of nodes.
c=1;
x=linspace(0,L,n);
%finding dx based on whole length
dx=L/(n-1);
x_start=0.1;
x_end=0.3;
n_start= round((x_start/dx)+1)
n_end=round((x_end/dx)+1)

%calculations
%Creating an array of u as per the given values
u=1*ones(1,n);
u(1)=1;
u(n_start:n_end)=2;
dt=0.01;
uold = u;
u_initial=u;
t=0.4;

%Since we need to find out the solution till 0.4 second we have to carry out t/dt number of iterations
nt=t/dt;
%time loop
for k=1:nt
% space loop
for i=2:n
%first order backward differencing for space term and forward differencing for time
u(i)= uold(i)-(c*dt/dx)*(uold(i)-uold(i-1));
end
uold=u;

plot(x, u,'color', 'g')
hold on
plot(x,u_initial,'color', 'y');
axis([0 1 0 2]);
pause(0.004);
xlabel('length of domain');
ylabel('velocity');
legend('final velocity','initial velocity');
end.

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

Aim: Our objective is to solve the PDE through MatLab code and compare the effect of time step on the solution and compare the simulation time taken for different time steps.

Step 1:Here we are creating a function which gives postion of node on whole length and final velocity value at that node as output by taking dt as input.

Step 2: giving inputs for the linear convection problem

'L' domain length = 1 meter

t = 0.4 sec

c=1

n=80

Step 3:

U array has been defined as U is 1 at all the X values excepts 0.1 to 0.3 band

So we created a function that predict the nodes at 0.1 and 0.3 and gives you the node number.

So using range operator we defined the value of U in between these nodes as 2

Step 4 :

The function has been carried out through two loops one for to move the step function spacely and other timely, where we used first-order forward differencing for time and backwards differencing for space.

step 5:

Then in the main program we are calling this function and calculating the values of final velocity at each node for each time step.then we also calculate the simulation time and compare them at different time steps.

Result:

figure 1

In figure 1 the graph shows the numerical solutions at t = 0.4 seconds for different time steps are compared

figure 2

Figure 2 shows the simulation time as a function of the time step.

Here the small step size took a large simulation time as the number of iterations gets increased nt=t/dt where t is a constant value of 0.4

when

t=0.1 thevelocity profile gets to blow up because of CFL value getting increased more than 1

t=0.01 the velocity profile looks great and no blow up and simulation time is less

t=0.001 as the step time gets lowered the finer the velocity profile becomes and simulation time is more than for 0.01

t=0.0001 it is even more smooth and the simulation time will be very high as the number of iterations is high.

program :

function:

% function which gives postion of node on whole length and final velocity value at that node as output by taking dt as input
function [x,uold] = timestep(dt)

%inputs for the linear convection problem

L=1; %1 meter length
n=80;
c=1;
x=linspace(0,L,n);
%finding dx based on whole length
dx=L/(n-1);
x_start=0.1;
x_end=0.3;
n_start= round((x_start/dx)+1)
n_end=round((x_end/dx)+1)

%calculations
%Creating an array of u as per the given values
u=1*ones(1,n);
u(1)=1;
u(n_start:n_end)=2;

uold = u;
t=0.4;

%Since we need to find out the solution till 0.4 second we have to carry out t/dt number of iterations
nt=t/dt;
%time loop
for k=1:nt
% space loop
for i=2:n
%first order backward differencing for space term and forward differencing for time
u(i)= uold(i)-(c*dt/dx)*(uold(i)-uold(i-1));
end
uold=u;

end
end

main program:

clear all
close all
clc
dt=[1e-4,1e-3,1e-2,1e-1]
% calling the function and finding the simulation time
tic
[x1,u1] = timestep(1e-4)
toc
t1=toc;
tic
[x2,u2] = timestep(1e-3)
toc
t2=toc;

tic
[x3,u3] = timestep(1e-2)
toc
t3=toc;

tic
[x4,u4] = timestep(1e-1)
toc
t4=toc;
figure(1)
% plotting the velocity profile for different time steps
plot(x1, u1,'color', 'g')
hold on
plot(x2,u2,'color', 'y');
hold on
plot(x3, u3,'color', 'r')
hold on
plot(x4,u4,'color', 'b');
axis([0 1 0 2]);
xlabel('length of domain');
ylabel('velocity');
legend('1e-4','1e-3','1e-2','1e-1');

figure(2)
% comparing the simulation time for different time steps on bar graph
simulation_time = [t1 t2 t3 t4]
bar(simulation_time,'g')
set(gca,'xticklabel', {'1e-4','1e-3','1e-2','1e-1'});
xlabel('time steps');
ylabel('simulation time');
title(' simulation time as a function of the time step');
legend('simulation time');