Effect of grid on 1D Linear Convection equation

Objective:

To write a code solve the 1D linear convection equation for the various grid points 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 grid points(4 values), constant, time dx and dt.
• Create a for loop for the range of length of n and create an array of 'x' by using linspace command and calculate dx value.
• For the given value of n calculate the velocity profile.
• Define the uold and u_init values and create a time loop with the range of 'nt'. Then, calculate the velocity at each grid size and by using backward differenation scheme.
• Then plot the graph between the x value and the updated velocity value and compare it with the original profile
• Plot the graph for the range of that are specified in the problem.

Matlab Script:

% Program to determine the effect of Grid size on 1D Linear Convection Equation.

clear all;
close all;
clc;
% Input for the linear convection problem
L = 1; % Length of the domain(assumation)
n = [20,40,80,160]; % no.of grid points
c = 1; % linear convection velocity 'c'(constant)
time = 0.4;
dt = 1e-2; % time step
x_start = 0.1; % profile should start
x_end = 0.3; % profile should end
ct = 1;
% calculation 
for i = 1:4
x = linspace(0,L,n(i));
dx = L/(n(i)-1);
u = ones(1,n(i));
for j = 1:n(i) % Nested loop 
if (x(j)>=x_start && x(j)<=x_end) 
    u(j) = 2;
else
    u(j) = 1;
end
end

uold = u;
u_init = u;
nt= time/dt;
% Time loop
for j = 1:nt
% Space loop BDS 
% du/dt = -c*du/dx
% du/dx - BD Scheme

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

figure(i)
grid on
plot(x,u_init,'k')
hold on
plot(x,u,'b')
xlabel('X axis')
ylabel('Velocity')
axis([0 1 0.5 2.5])
pause(0.3)

 M(ct)=getframe(gcf);
 ct= ct+1;

end
end
movie(M)
file=VideoWriter('Linear Convection.avi','uncompressed AVI')
open(file)
writeVideo(file,M)
close(file)

Result:

For n =20

From the above graph, it is clear that the profile had been moved as the length of the domain increases.

For n =40

From the plot, It is clear that the velocity profile is been moving as the length of the domain increases but the diffusion rate had been decreased.

For n = 80

From the plot, It is clear that the velocity profile is moved as the length of the domain increases

For n=160

From the graph, It is clear that the results had been blown up i.e it has an unstable solution at 160 grid points.

Animation Link:

The animation attached below shows how the velocity profile changes as the change in the grid size. At 'n=80', the solution is stable compared to other grid sizes.