1D LINEAR CONVECTION EQUATION SIMULATION

OBJECTIVE:

To arrive a numerical solution to the wave equation,

1. Let us assume that the domain length is 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. Let us use first order forward differencing for the time derivative.

4. The let us ue the first order backward differencing for the space term.

5. Time step and number of grid points be the variables in the calculation.

6. C = 1.

7. Set time step = 0.01.

EXPECTED RESULTS:

1. For C = 1 and for time = 0.4seconds compare the original and final velocity profiles.

2. Making the comparision for n=20,40,80 and 160.

3. We will get 4 plots in total.

4. To explain the plots.

ABSTRACT:

The PDE for the given problem is dU/dT + C*(dU/dX) = 0;

Solving the first PDE of the course, this is a very simplified PDE yet vastly essential for the bigger problems to be solved. This is also called 1D Linear Convection. Now 1D is because it repsents only one spacial axis, linear is because the constant C is getting multiplied with the velocity. Convective term is because the term to which the constant C is getting multiplied is very similar to the convective term in the navier stokes eqn. 

This will help us understand basic numerical methods, stability criterias and have a strong foundation in CFD.

This eqn can be solved in a time marching procedure or time intergration.

THE CODE:

clear all
close all
clc

% Length of the domain in meters;
L = 1;

% Number of grid points;
n = 20;

% X array or the mesh, the total length of the domain;
x  = linspace(0,L,n);

% Grid spacing = x(2)-x(1);
dx = L/(n-1);

% Time step in seconds;
dt = 0.01;

nt = (0.4/dt) ;

% Linear conviction velocity or the constant part;
c = 1;

% velocity at the starting step in X coordinate;
x_start = 0.1;

% velocity at the ending step in X coordinate;
x_end = 0.3;

% grid point at the start of velocity profile
n_start = (0.1/dx)+1;

% grid point at the end of velocity profile
n_end = (0.3/dx)+1 ;

% Initialising Array and applying boundary condition
u = ones(1,n);
u(n_start:n_end) = 2

u_old = u;
u_initial = u;
tic
for k = 1:nt
    
    for i = 2:n
        u(i) = u_old(i) - (c*dt/dx)*(u_old(i) - u_old(i-1));
    end
    
    % update old values
    u_old = u;
 
end
toc
    plot(x,u_initial)
    hold on
    plot(x,u,'r')
    axis([0 1 0 3])
    xlabel('domain length')
    ylabel('velocity')
    xlabel('Length domain')
    legend('Initial velocity profile','final velocity profile')
    grid on

OUTPUT:

The first time i solved the eqn, I went on with a number of grid points as, n = 20;

The elapsed for the step was 0.038 seconds.

The output graph,

The I pushed the grid size to 50, 100, 150 and upto 200 until the solution became unstable. The graphs looked wiggled and blew up. This happens not only in simple simulations but also in 3D simulation. We can analyse when the solution will go unstable with Von Neumann method of stability analysis.

The time elapsed for grid size of 50 is 0.07 seconds.

The time elapsed for grid size of 100 is 0.14  seconds.

The time elapsed for grid size of 150 is 0.16 seconds.

The time elapsed for grid size of 200 is 0.19 seconds.

The screenshots that support the output are here as follows,

 CONCLUSION: Hence we have arrived at our desired output and results compared.

 The EXPLICIT METHOD is preferred here to the implicit but at the cost of computation time & cost in order to attain a stablility.

However, we can go for the IMPLICIT METHOD as it is UNCONDITIONALLY STABLE as time step effect will not be accounted here.

But the accuracy of the output is to be compromised when the results are compared to the explicit method.