Effect of Time step on 1D Linear Convection.
Objective:
To write a code solve the 1D linear convection equation for the various time step and to generate the plot for the velocity profile using the Matlab.
Convection: The flow that combines diffusion and the advection is called convection.
Advection: The bulk transport of mass, heat or momentum of the molecules.
Diffusion: Momentum of the molecules from the area of higher concentration to the lowe concentration is called Diffusion.
2D Navier Stroke Equation:
`∂u/∂t+u*∂u/∂x+v*∂u/∂y= -∂P/∂x+u*((∂^2 u)/(∂x^2 )+(∂^2 u)/(∂y^2 ))`
Here `∂P/∂x` - Pressure drop term.
`u*∂u/∂x+v*∂u/∂y`- Convective term.
`((∂^2 u)/(∂x^2 )+(∂^2 u)/(∂y^2 ))`- Diffusive term.
Linear Convection Term:
`∂u/∂t+c*∂u/∂x` - Linear convection term is solved by the time marching solution.
Here, We are using the forward differentiation scheme(FDS) for the time loop and backward differentiation scheme(BDS) for space loop.
The formula for the BDS is given by
`(u_i^(n+1)-u_i^n)/∆t= -c*((u_i^n-u_(i-1)^n)/∆x)`.
Basically, there are three solver methods for discretization they are
- Implicit method
- Explicit method
- Crank Nichalsion method (Implicit + Explicit method).
Code:
- Define inputs that are required to solve the 1D linear convection suck as the length of the domain, the array of the time steps(4 values), constant, time dx and grid size.
- Define a function called 'timestep' to calculate the velocity for the range of 'x'. Then calculate the velocity for the range of 'dt' and plot the graphs
- The simulation time taken for each time step is calculated by using 'tic-toc' command.
- Plot a bar graph of the simulation time and label it.
Matlab Script:
clear all;
close all;
clc;
% Input for the linear convection problem
L = 1; % Length of the domain(assumation)
n = 80; % no.of grid points
x= linspace(0,L,n);
dx = L/(n-1);
c = 1; % linear convection velocity 'c'(constant)
time = 0.4;
x_start = 0.1; % profile should start
x_end = 0.3; % profile should end
dt = [1e-4,1e-3,1e-2,1e-1];
for i=1:4
tic;
[x,u] = timestep(dt(i),L,n,x,dx,x_start,x_end,time,c)
t(i) =toc
% plotting
subplot(2,2,i)
plot(x,u)
grid on
hold on
xlabel('Length of Domain')
ylabel('Velocity(u)')
end
figure(2)
bar(t)
set(gca,'XTickLabel',{'1e-4','1e-3','1e-2','1e-1'})
title('Simulation Time at different time step')
Function Script:
function [x,u] = timestep(dt,L,n,x,dx,x_start,x_end,time,c)
n_start = (x_start/dx)+1;
n_end = (x_end/dx)+1;
u = ones(1,n);
u(n_start:n_end) = 2;
% Boundary Condition
uold = u;
u_initial = u;
%nt = 100; % total time we need to perform integration
nt = time/dt;
% time loop
for j = 1:nt
% Space loop BDS
% du/dt = -c*du/dx
% du/dx - BDS FDS Skewed schemes
for i = 2:n
u(i) = uold(i)-(c*(dt/dx)*(uold(i)-uold(i-1)));
end
uold = u; % update the old velocity
end
end
Result:
The plots represented as 1e-4(2,2,1),1e-4(2,2,2),1e-3(2,2,3), 1e-1(2,2,4)
From the graph, It is clear that as the time step decreases (right to left) the solution becoming unstable. This is because an explicit scheme is sensitive to the time step. In order to avoid this, it is preferred to use the implicit method, Implicit methods are not to the step size. To determine the accurate time step we have to consider the stability criteria like CFL(Courant–Friedrichs–Lewy ), Von Neumann stability condition and Eigen stability condition.
The CFL condition expresses that the distance that any information travels during the timestep length within the mesh must be lower than the distance between mesh elements. In other words, information from a given cell or mesh element must propagate only to its immediate neighbors.
The CFL stability condition for the given linear convection term is that 'C<=1'
Plot for the Simulation:
From the graph, it is clear that as the time step increases the simulation time also increased. This phenomenon is due to an increase in the number of iterations.