Numerical solution to the 1D-Linear Convection equation for varying Time steps using MATLAB

OBJECTIVE:

To solve the one-dimensional Linear Convection equation and write a function that accepts the time step as an argument comparing it's effects on numerical solution using MATLAB.

Project Conditions:

1. Assume that the domain length, 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. Grid points, n = 80

6. Constant Velocity, C = 1

7. Set time step, dt = 1e-4, 1e-3, 1e-2, and 1e-1

THEORY:

Linear Convection equation is the most accessible equation in CFD, from the Navier Stokes equation we kept only the accumulation and convection terms for the component of the velocity - as we already know, in CFD the variables to be computed are velocities. To make things even simpler, the coefficient of the first derivative of the velocity is constant, making the equation linear.

where,

`(∂u)/(∂t)` = time derivative

`(∂u)/(∂x)` = convective terms

This equation represents the propagation of a wave without change of shape; the speed of the wave is and with the initial condition given by:

The exact solution is

u(x,t)=u0(x−c⋅t)">u(x,t)=u0(x−c⋅t)u(x,t)=u0(x−c⋅t)

This means that on an axis with the positive direction to the right, the wave will move from left to right with the speed c.

To discretize this equation in space and time we will use the Forward Difference scheme for time derivative and Backward Difference scheme for space derivative.

The partial derivative can be approximated using the derivative definition, removing the limit operator:

`(∂u)/(∂x)≈(u(x+Δx)−u(x))/(Δx)`

This approach will be applied both for space and time; for time steps we will use the superscript n and for space steps (grid cells) we will use the subscript i:

where n+1 and n are two consecutive steps in time and i and i−1 are two neighboring points of the discretized x coordinate (grid points). The only unknown in the equation written in this way is un+1i. We must solve for this unknown to get an equation that allows us to advance in time, as follows:

Here we follow the method of Time Marching, which is well illustated in the image below

PROGRAM:

Program to determine the function for time steps

%Linear Convection - Effect of time step size on the solution
%Program for determining function

function[u, x] = Time_steps(dt)

L = 1;                  %Length of Domain
C = 1;                  %Constant Velocity
n = 80;                 %No. of grid points
t = 0.4;                %Simulation Time

%No. of time steps
nt = t/dt;

%Mesh Generation
x = linspace(0, L, n);
dx = L/(n - 1);


%Define the start and end positions when velocity is 2m/s
x_start = 0.1;
x_end = 0.3;
%Function of start and end values
n_start = (x_start + dx)/dx;
n_end = (x_end + dx)/dx;


%Boundary Condition when velocity is 1m/s
u = ones(1,n);
%Boundary condition when velocity is 2m/s
u(n_start:n_end) = 2;


%Assigning u to u_old
u_old = u;
u_initial = u;

%TIME MARCHING
    %TIME LOOP
    for i = 1:nt 
        
        %SPACE LOOP
        for j = 2:n
            %Difference Equation
            u(j) = u_old(j)-(((C*dt)/dx)*(u_old(j) - u_old(j - 1)));
        end
        
        %Updating Old values
        u_old = u;
    end

end

Explanation

• A function is coded to determine the varying time steps(dt) for outputs (u,x).
• Input parameters are assigned for Length(L), Constant velocity(C), no. of grid points(n), Simulation time(t).
• The formula for no. of time steps is assigned to nt.
• In order for the mesh to generate, x and dx values are inputed.
• Next step we have defined the start and end positions denoting 0.1 & 0.3 at 2m/s and rest of others at 1m/s.
• The values of n_start, n_end, boundary condition for u at 2 & 1m/s are assigned here.
• Since the array u contains new values it can be reused for the time & space loop, so we assign u again with old values (i.e) u_old, u_initial.
• Space Loop is now inputed between Time loop in a function with the difference equation within the know grid value, n = 80
• Old value of u is been reused.
• Function is ended.

Main Program

clear all
close all
clc

%Linear Convection - Effect of time step size on the solution
%Main Program

tic
[u1,x1] = Time_steps(1e-1)
toc
Time_1 = toc

tic
[u2,x2] = Time_steps(1e-2)
toc
Time_2 = toc

tic
[u3,x3] = Time_steps(1e-3)
toc
Time_3 = toc

tic
[u4,x4] = Time_steps(1e-4)
toc
Time_4 = toc

%Plotting 
%1D Linear Convection Velocity profile at various Time steps
figure(1)
plot(x1,u1,'linewidth',1,'color','b')
hold on
plot(x2,u2,'linewidth',1,'color','r')
hold on
plot(x3,u3,'linewidth',1,'color','g')
hold on
plot(x4,u4,'linewidth',1,'color','y')
hold on

xlabel('Domain Length')
ylabel('Velocity')
legend('Time step = 1e-1','Time step = 1e-2','Time step = 1e-3','Time step = 1e-4')
axis([0 1 0.5 2.5])
title('1D Linear Convection Velocity profile at various Time steps')
grid on

%Bar Graph for Simulation Time vs Time steps
figure(2)
Simulation_time = [Time_1,Time_2,Time_3,Time_4];
h = categorical({'1e-1','1e-2','1e-3','1e-4'});
bar(h,Simulation_time,'b')

xlabel('Time steps')
ylabel('Simulation Time')
title('Simulation Time vs Time steps')

Explanation:

• Program is started with commands clear all, close all and clc to establish a fresh start without the manipulation of previous values.
• Time steps is called by function for the given four different values.
• In order to measure the time between the time steps, tic toc command is used
• Plotting is done for both Linear Convection Velocity at various time steps and Simulation time between time steps.

OUTPUTS:

Plot : Comparsion of Velocity profile for various time steps

Observation:

1. Time step 1e-1 shows the most unstable solution, the graph almost blows up in both directions. Reasons may have been due to Explicit approach (where time step should not be greater than stability criterion)
2. Time step 1e-2 has optimium velocity nearly equal to the original curve which makes it the most stable solution and gives less truncation error.
3. Time step 1e-3 & 1e-4 both gives almost same velocity profile which shows that velocity has been defused resulting in less stable solution when compared with 1e-2.

Bar graph : Comparsion of Simulation Time for different time steps

Observation:

1. Time step 1e-1 results blows up after a point and hence it takes large time (This cannot be considered since the result is not accurate)
2. Time step 1e-2 is the stable solution of all, displays the lowest time taken.
3. Time step 1e-3 is similar to 1e-2 in regards of time taken, even though its result is less accurate.
4. Time step 1e-4 takes more simulation time because it is a very small time step which makes it consume time to travel all grids.

CONCLUSION:

From the above plots we can conclude that

• Larger the value of time step, lesser the accuracy of results causing even errors
• Lower value of time step shows positive results but not to low which may increase the simulation time according to CFL number.
• 1e-2 is the most stable and less time consuming time step compared to others.
• 1e-1 is the most unstable of all.