CFD analysis of 1D Linear equation using Matlab.

The first step to computing the Navier-Stokes’ equation is a wave propagation equation in 1D. The 1 dimensional linear convection equation:

`(del u)/(del t) + c (del u)/(del x) = 0`

'u' is a quantity that is transported at a constant velocity. The constant velocity makes the equation linear in nature contrary to NS equation which is nonlinear. The physical behaviour of the above equation is that the initial profile at u(x,t=0) = u0(x) is getting transporated depends upon the magnitud of constant velocity c. If the constant velocity 'c' is +ve then initial profile will be tranported rightward and if -ve, it will transport leftward.

To compute this numerically, a space-time discretization is used. In particular, the “Forward difference in time, backward difference in space” numerical scheme is used. This is a finite difference scheme, which takes the form:

`(u_i ^(n+1) - u_i ^n)/(Delta t) +c * (u_i ^n - u_(i-1) ^n)/(Delta x) = 0`

Transposing the above dicretized equation, we can solve for the next time step and it is called time marching method.

`u_i^(n+1) = u_i^n - (c*Delta t)/(Delta x) * (u_i^n - u_(i-1)^(n))`

The initial Profile:

Domain Length - 1m

Initial Condition: u = 2m/sec @0.1-0.3m, else 1m/sec.

Boundary Condition: u =1 m/sec at x =0m & x =1m

Expected output:

Velocity profiles for different timestep and grid size.

Analysis:

I am taking 4 different timsteps and 4 different nodes to check the transported velocity profile for the time period of 0.4 sec. The CFL condition is a important criteria to consider to get the stable solution.

The timsteps chosen are 1e-4, 1e-3, 1e-2 and 1e-1. The number of nodes(Grid Size) chosen is 20,40,80 and 160. The CFL number for all the combinations is,

dx	0.05	0.025	0.0125	6.25E-03
Timestep	n=20	n=60	n = 80	n = 160
1.00E-04	0.002	0.00625	0.008	0.016
1.00E-03	0.02	0.0625	0.08	0.16
1.00E-02	0.2	0.625	0.8	1.6
1.00E-01	2	6.25	8	16

The CFL number more than 1 should give unstable solution and the matlab code to test the claim is,

clear all
close all
clc

% Input for the linear convection problem

n = input('Number of Nodes(21,61,81,161)= ');
dt = input('The size of the timestep(1e-4,1e-3,1e-2,1e-1)= ');

L = 1;      % In meters
c = 1;      % Constant Velocity
t = 0.4;    % Time
nt =t/dt;   % Number of Timestep

x = linspace(0,L,n);    % Domain Size
dx = x(2)-x(1);         % dx = L/(n-1), grid size
Cu = c*dt/dx;           % Courant Number

% Define the start and end position for square wave
 
x_start = x==0.1;
x_end = (x==0.3);
n_start = find(x_start==1);
n_end = find(x_end==1);

% Initialize arrays and apply boundary conditions

u = ones(1,n);
u(n_start:n_end) = 2;

% Boundary Conditions
uold = u;
u_original=u;
figure(1)
plot(x,u_original,'b','lineWidth',2 )
hold on

for k =1:nt
        
    for i = 2:n
            
        u(i) = uold(i)-(Cu*(uold(i) - uold(i-1)));
        
    end
    % update old velocities
    uold =u;
end
    
    %Plotting the Final velocities
    
    plot(x,u,'lineWidth',1)    
    legend('Initial Velocity','Final Velocity','FontSize', 13);
    xlabel('Length of the Domain (m)','FontSize', 15);
    ylabel('Linear Velocity (m/sec)','FontSize', 15);
    ylim([0.5 2.5]);
    title(sprintf('Solution for,n = %d, dt = %f, CFL = %f',n,dt,Cu));

The results for few combinations of timestep and stepsize are,

Discussion:-

1. The first order backward(Upwind scheme) creates false diffusion in the profile because there is no diffusion term in equation.

2. The courant number more than 1 gives unstable solution.

3. The false diffusion seems to have less impact when grid size goes on decreasing or we can say that the error is increasing with increasing grid size.