Week 4.2 - Linear Convection for varying timesteps

QUESTION:

Effect of time step size on the solution

Write a function that accepts the time step as an argument and solves the same problem. You will compare the effect of the time step on the numerical solution.

1. Number of nodes, n = 80
2. Time step, dt = 1e-4, 1e-3, 1e-2, and 1e-1
3. Simulation run time, t = 0.4 seconds.

Expected output:

1. In a single graph compare the initial and final profiles of the numerical solution for all dt's.
2. Explain in detail the CFL criteria that govern the linear convection equation. Do these criteria apply to an implicit solution as well?
3. From the plots obtained: Explain how the change in time step size affects numerical diffusion, instability, the courant number, and accuracy. How does change in time step size affect these phenomena?
4. Plot the simulation time, t as a function of the time step, dt.

ANSWER:

AIM:

• To write a program to solve the linear convection equation with varying time steps.
• To compare the velocity profile of the numerical solution for all timesteps.
• To explain in detail the CFL criterion that governs the linear convection equation.
• To explain how the change in time step size affects the numerical diffusion, instability, the courant number, and accuracy.
• To plot the simulation time, t as a function of the time step, dt.

Assumption and input:

• Assume that the domain length is L = 1m
• 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.
• Use first-order forward differencing for the time derivative.
• Use first-order backward/rearward differencing for the space term.
• Constant convection coefficient, c = 1 m/s
• of nodes = 80
• Time step, dt = 1e-4, 1e-3, 1e-2, and 1e-1
• Simulation run time, t = 0.4 seconds.

Solution:

The one-dimension linear convection equation is written as,

Where u = Velocity and c = Convective coefficient.

Given that we have to use first-order forward differencing for the time derivative term and first-order backward differencing for space term.

Rewriting the above equation, we get,

The above equation is used in MATLAB code to obtain a solution.

INPUTS:

[MATLAB CODE]

%% Clearing
clc
clear all
close all

%% Inputs 

n = 80;                         % Number of nodes
l = 1;                          % Domain length(m)
c = 1;                          % convective coefficient(m/s)
dt = [1e-1,1e-2,1e-3,1e-4];     % Time step
t = 0.4;                        % Simulation run time


%% Solution

tic
[x1,u1] = linear_convection_timestep(l,n,dt(1),c,t);
time1 = toc;

tic
[x2,u2] = linear_convection_timestep(l,n,dt(2),c,t);
time2 = toc;

tic
[x3,u3] = linear_convection_timestep(l,n,dt(3),c,t);
time3 = toc;


tic
[x4,u4] = linear_convection_timestep(l,n,dt(4),c,t);
time4 = toc;

%% Plotting 

% Velocity plot
figure(1)
plot(x1,u1,'r','linewidth',2);
hold on
plot(x2,u2,'g','linewidth',2);
hold on
plot(x3,u3,'b','linewidth',2);
hold on
plot(x4,u4,'black','linewidth',2);
xlabel('Size of Domain');
ylabel('Velocity');
axis([0 1 0 2])
legend('1e-1','1e-2','1e-3','1e-4')
title('Velocity profile at verying time step')
grid on
pause(0.003);

% Simulation time taken at various time step
time_step = categorical({'1e-4','1e-3','1e-2','1e-1'});
simulation_time = [time1 time2 time3 time4];
figure(2)
bar(time_step,simulation_time);
xlabel('Varying time step');
ylabel('Simulation Time');
legend('Simulation time for various time step');
grid on

[MATLAB FUNCTION CODE]

function [x,u] = linear_convection_timestep(l,n,dt,c,t)
    
    x = linspace(0,l,n);        % Length of Domain    
    dx = l/(n-1);               % Grid Spacing

    %% To define the start and endpoint of the square wave

    x_start = 0.1;
    x_stop = 0.3;

    %% Time Marching Solution

    nt = t/dt;

    n_start = round((x_start/dx)+1);
    n_stop = round((x_stop/dx)+1);

    %% Square wave
    u = ones(1,n);
    u(n_start:n_stop) = 2;
    %% Boundary Condition
    u(1) = 1; % Defining the value to the initial node, but in square wave section we have alreaddy defined the value
    uold = u;

    %% 
    for k = 1:nt
        for i = 2:n
            u(i) = uold(i) - (c*dt/dx)*(uold(i) - uold(i-1));
        end

        uold = u;

    end
end

EXPLANATION:

• Initially, the input values are assigned to initialized to the solution.

n = 80; % Number of nodes

l = 1; % Domain length(m)

c = 1; % convective coefficient(m/s)

dt = [1e-1,1e-2,1e-3,1e-4]; % Time step

t = 0.4; % Simulation run time

x = linspace(0,l,n); % Length of Domain

dx = l/(n-1); % Grid Spacing

• Here ‘n’ defines the ‘Number of Nodes’ in which the whole velocity profile is distributed.
• ‘l’ is the length of domain.
• ‘dt’ defines the varying time step value.
• ‘t’ is the total simulation run time.
• ‘dx’ is the grid spacing.
• In order to calculate the simulation for each time step, the command tic and toc are used. Tic command starts the timer and Toc command reads the elapsed time.
• Between the Tic and Toc command, the function linear_convection_timestep is called.
• Now the start and end positions of the square wave are defined. (x_start & x_stop).
• 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.
• Now to find the n_start and n_stop values following expression is used, n_start = round((x_start/dx) +1). Similarly for n_stop.
• Here, x_start/dx+1 will give the nearest point to x_start value. And the round command is used to round off the decimal values.
• The velocity 2m/s is assigned within the points n_start and n_stop, and elsewhere 1m/s is assigned.
• Next, the boundary conditions are assigned to the velocity profile.
• Inside the time loop, the space loop is assigned and is been iterated up to value(n).
• The equation solved in space loop is,

• When the space loop ends, the old velocity value is updated.
• Finally, the time loop comes to an end and the function call is stopped.
• Simulation time is got stored in the variable time1, time2, time3, time4.

RESULTS:

Figure: 1 – Velocity profile at various Time steps.

Explanation:

• Figure 1 shows the final velocity profile for the different time step values.
• For dt = 1e-1 the solution blows up and ripples are seen. No proper velocity profile is seen.
• For dt = 1e-2, the peak value of velocity profile is nearly 2m/s. So, there is only a little artificial diffusion.
• For dt = 1e-3 and dt = 1e-4, the peak value for final velocity profile is less than its initial velocity (2m/s), due to artificial diffusion.

Figure: 2 – Simulation time for Various Time steps.

Courant-Friedrichs-Lewy (CFL) condition:

• It is necessary to condition for PDE numerically. It arises in the numerical analysis of the explicit time integration scheme.
• This CFL number determines the level of instability in the particular problem.
• It is basically the stability criterion for a stable solution.
• Any linear explicit convection problem must be equal to or smaller than 1 to give the stable and appropriate solution.

The CFL number is a dimensionless quantity expressed as,

Effect of change in time step size:

Calculation of Courant number and Artificial numerical diffusion:

• convection coefficient, c = 1 m/s
• dx = 0.0127

• Courant Number –

• Artificial Numerical Diffusion –

• Numerical Diffusion is negative in this case. (Unstable)
• From the above four cases, we can see that how the variation in time step affects the courant number and hence affect the velocity profile.
• In cases 1,2 and 3, the solution is stable whereas in cases 4 the solution becomes unstable.

CONCLUSION:

• Numerical solution to the one-dimensional linear convection equation is derived and the simulation for various time step 1e-1, 1e-2, 1e-3, 1e-4 are performed.
• The velocity profile of the numerical solution is compared with respect to the time step.
• Simulation time is plotted as the function of the time step.
• CFL criteria that govern the linear convection equation are explained.
• The effect of change in time step which affects the numerical diffusion, instability, the courant number, and accuracy is explained.