Effect of the time step on the solution of a linear convection equation.

Aim: In this particular challenge, I am going to look at the effect of the time step on the solution of a 1D-linear convection equation.

Governing Equations and discretization:

To solve the equation:

Since the value of C is assumed to be a constant, the above equation is a linear in character.

The time derivative of the above equation is discretized using the forward differencing and the space derivative of the above equation is discretized using the backward differencing scheme. After discretisation the above equation looks something like this.

Since the above equation is explicit in nature, problem of instability comes into place. In the following simulation, number of grid points and the thermal diffusivity are maintained constant. The CFL number is varied by varying the time step.

Problem Setup:

1)Length of the domain, l=1m

2)C=1

3)n=80 %No. of grid points

4)dx=l/(n-1)

5)time steps used= 1e-4, 1e-3, 1e-2, 1e-1.

Expected Solution:

Since there are no diffusive terms in the governing equation, the shape of the initial square pulse is expected to remain same. Because of the advective term in the governing equation, the square pulse is expected to move along the x direction with a speed of 1m/s.

Method of Solving the above problem:

To solve the above partial differential equation we use the finite difference method.

We use the Forward difference method for the time derivative and backward differencing method for spatial derivative. In this problem we assume C to be a constant with a magnitude of 1. The time derivative is taken at a previous time step thereby making the discretization explicit. The initial and boundary conditions are described as shown in the figure.

After discretizing and transposing the PDE shown above we get the following algebraic equation.

Where n is the time index and i is the spatial index.A time marching solution can be obtained using the following algebraic equation.

Expected Solution:

When the initial solution of the above equation is plotted on a x-u axis we get a pulse. Since there is only convective term in the above equation we expect the pulse to move through a specific distance with time without much change in the shape of the pulse.

Coding and Simulation:

clear all
close all
clc


tic 
[u0 x0 init]=convection2(0.1)
toc
soltime1=toc;

tic
[u1 x1]=convection2(0.01)
toc
soltime2=toc;

tic
[u2 x2]=convection2(0.001)
soltime3=toc;

tic
[u3 x3]=convection2(0.0001)
soltime4=toc;

figure(1)
plot(x0,init,'-k','LineWidth',3)
hold on
grid on
plot(x1,u1,'-b','LineWidth',2)
plot(x2,u2,'-k*','LineWidth',1)
plot(x3,u3,'-m','LineWidth',1)
hold off
axis([0 1 0 2.5])
legend('u at t=0','u at t=0.4 and dt=0.01','u at t=0.4 and dt=0.001','u at t=0.4 and dt=0.0001','location','southeast')
xlabel('x-values')
ylabel('u-values')
title('comparision for different time steps')
%saveas(gca,'Comparision_of_all1.jpeg')

%bar graph of time step vs solution time
figure(2)
y=[soltime1 soltime2 soltime3 soltime4]
bar(y,0.3,'k','EdgeColor','r')
set(gca,'XTick', 1:1:4)
set(gca,'xtickLabel',{'dt=0.1','dt=0.01','dt=0.001','dt=0.0001'})
ylabel('Solution time')
ylim([0 600])
legend('Comparision of solution time with time step','location','northeast')
labels={'0.4817', '5.0713', '42.3163' ,'480.0106'}
x=[1,2,3,4]
text(x,y,labels,'HorizontalAlignment','center','VerticalAlignment','baseline')
saveas(gca,'stepsizevssoltime.jpeg')

function file:

function [u,x,init] = convection2(dt)

%defining_parameters
l=1;
n=80
x=linspace(0,1,n)
dx=(l/(n-1));
c=1;
t=0.4
nt=t/dt
cou=(c*(dt)/(dx))

%solution_initialization @t=0
u=1*ones(1,n)
%from x=0.1 to 0.3 ,u=2
xstart=0.1
xend=0.3
y=abs(x-xstart);
z=abs(x-xend);
nstart=find(y==min(y));
nstop=find(z==min(z));
u(nstart:nstop)=2;
init=u
uold=u
for j=1:nt
    for i=2:n
    u(i)=uold(i)+((cou)*(uold(i-1)-uold(i)))
    end
    uold=u     
end

After executing the above codes, the following results were obtained

Discussions:

When dt=0.1 courant number=7.9

for dt=0.01,CFL number=0.79

for dt=0.001, CFL number=0.079

For dt=0.0001, CFL number=0.0079

1) As expected, In the first case, as the courant number crossed the value of 1, the solution became highly unstable. When the courant number employed reached a value close to 1, the result obtained resembled the expected solution. Therefore efforts should me made to adjust the time step and the grid spacing in such a way that the CFL number remains close to 1.

2)After a particular point, further reduction in the value of the time step doesnot alter the solution to a significant extent. This can be clearly seen in the above graphs.