1D Linear convection with different grid and dt value

Objective:

The main objective of this project is to write a code in matlab for solving 1D linear equation. Also to understand the basic concepts of CFD and to perform the square wave equation.

CFL condition:

The Courant–Friedrichs–Lewy or CFL condition is a condition for the stability of unstable numerical methods that model convection or wave phenomena. 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 nodes.

In physical explanation,we can say that for a CFD simulation, the Courant number could be that it tells you something about how fluid is moving through your computational cells. Meaning that If the Courant number is <=1 then the fluid particles move from one cell to another within one time step (at most). While if the Courant number is >1 a fluid particles moves through two or more cells at each time step and this can affect convergence negatively.

CFL Derivative:

We can derive the CFL condition in a simple and intuitive way. Let’s begin by considering the simple linear convection problem of a quantity u:

For the simple problem of linear convection of a quantity u, the first order explicit upwind scheme becomes:

where a is the velocity magnitude, Δt is the timestep and Δx is the length between mesh elements.

By performing a Taylor series expansion on this scheme, this scheme introduces what is called a numerical diffusion, or numerical viscosity, equal to:

As expected from a diffusion coefficient, the numerical diffusion coefficient must be positive, otherwise, the solution will grow indefinitely with respect to time, making the numerical scheme unstable.

The right term inside the parenthesis of the above expression is commonly referred to as the Courant number, which is a dimensionless quantity. Therefore, the Courant number can be stated as follows:

It follows from the numerical diffusion coefficient discussion, that for any explicit simple linear convection problem, the Courant number must be equal or smaller than 1, otherwise, the numerical viscosity would be negative:

The above discussion may only hold true for explicit schemes.

The implicit linear convection is unconditionally stable, and also lacks a CFL condition. This happens because implicit schemes use the entire domain to calculate each timestep. However, implicitly calculating (i.e., solving an implicit linear system) each timestep can be very costly. In addition, time accuracy decreases as the timestep increases. Therefore, it is logical is to use an explicit scheme with a lower timestep.

For more complicated and representative problems, including nonlinear flow, the CFL condition serves as a reference for setting the timestep. Furthermore, when dealing with nonlinear phenomena, even the implicit form can have stability problems if a large Courant number is chosen. Consequently, it might be necessary to apply local stability analysis on a linearized formulation in those cases. Still, the CFL condition is necessary for nonlinear problems.

von Neumann stability analysis:

In Numerical analysis von Neumann stability analysis (also known as Fourier stability analysis) is a procedure used to check the stability of finite difference schemes as applied to linear partial differential equations. The analysis is based on the Fourier decomposition of numerical error. This method is an example of explicit time integration, where the function that defines governing equation is evaluated at the current time.

If the problem of stability analysis can be treated generally for linear equations with constant coefficients and with periodic boundary conditions, as soon as we have to deal with non-constant coefficients and (or) non-linear terms in the basic equations the information on stability becomes very limited.Hence we have to resort to a local stability analysis, with frozen values of the non-linear and non-constant coefficients, to make the formulation linear. In any case, linear stability is a necessary condition for non-linear problems but it is certainly not sufficient.

Assumptions made to solve the Numerical solution to the wave equation

1. 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. First order forward differencing for the time derivative

4. First order backward/rearward differencing for the space term.

5. C = 1 is the constant velocity

6. Set time step = 0.01 and time = 0.4seconds.

7. n=20,40,80 and 160 is the number of grid values

Matlab program:

In below mentioned code, I have used for loop for calculation and also for declaring the start point and end point for the square wave, i have used 'if else' statement. For calculation part the nodes value (n) are declared as array and called from the for loop.

close all 
clear all 
clc

% Input parameters
L = 1;              % length in meter
n = [20,40,80,160]; % no of grid points
c = 1;              % velocity
t = 0.4;            % total time
dt = 0.01;          % time step width
%dt = t/nt; 
nt = t/dt;          % Defining no of time steps with the formula

% Calculation
% creating mesh
for p = 1:length(n)
x = linspace(0,L,n(p));
dx(p)= L/(n(p)-1);
CFL(p) = c*(dt/dx(p)); % CFL value calculation

% Defining start and end positions for the sq wave
% x_start = 0.1
% x_end =0.3

% Initializing array & Boundary condition with the for loop
u = ones(1,n(p));
for j = 1:n(p)
    if x(j)>=0.1 && x(j)<=0.3
        u(j) = 2;
    else
        u(j) = 1;
    end
end
uold = u;                
u_initial = u;
             
                          % du/dt = -cdu/dx  
                          % du/dx = BDS
% Time loop               % du/dt = FDS
for k = 1:nt              % defining no.of time steps in for loop
  % Space loop
  for i = 2:n(p)
     u(i) = uold(i) - (c*dt/dx(p))*(uold(i) - uold(i-1)); 
  end
  uold = u;                % updating velocity
  
  % plotting 
  figure(p)
  plot(x,u_initial,'r')
  hold on
  plot(x,u,'b')
  title({'1-D Linear Convection',['time(\itt) = ',num2str(dt*k)],['No.of grid = ',num2str(n(p))]})
  xlabel('Spatial co-ordinate (x) \rightarrow')
  ylabel('Velocity')
  axis([0 1 0 2])
  pause(0.3)

end
end

Plotting:

The plotting is done, which displays time duration and the node values for each plot.

title({'1-D Linear Convection',['time(\itt) = ',num2str(dt*k)],['node = ',num2str(n(p))]})

In the above command title, name for the pltting is given as 1-D Linear convection and the time iteration is converted into string character by 'num2str' command.

Plot:1

In the below plot, we can see here the velocity is getting decreased form the peak velocity 2m/s and diffusion rate is also more.

Plot:2

Here for the same time period by varying the no. of grids we can observa that diffusion rate is minimized compared to the prevoius plot.

Plot:3

In this plot, we achecived the accurate solution by observing that the velocity of the profie remains constant (which is not get decreased) in the peak value(2m/s).Since there is no change in the velocity we can say that there is no diffusion takes place.

Plot:4

The final plot having grid value 160 shows that the profile is blowing up contionously.This means that the velocity profile is in unstable condition.This is because of larger grid size.

Observation and Result:

From the above results we can conclude that if the CFL value is <=1, we will get the stable solution. If it exceeds, then it leads the solution to blowup the velocity profile in an unstable condition.

Program 2:

This program is done by changing the values of dt. That means dt value is defined in the form array [1e-4,1e-3,1e-2,1e-1].

The no of grid n = 80

Plotting is compared for different dt values.

Matlab program:

Function program:

Here in this function, the x value and u value is calculated from the function and stored. Then this function is called in the main program with different dt values and the time duation for each dt values is evaluted.

function [x,u] = linearconvection_func(L,n,c,t,dt,x)
dx = L/(n-1);
nt = t/dt;

% Defining start and end positions for the sq wave
x_start = 0.1;
x_end =0.3;
n_start = round(x_start/dx+1);
n_stop = round(x_end/dx+1);

% Initializing array & Boundary condition
u = ones(1,n);
u(n_start:n_stop) = 2;
uold = u;

% Time loop 
for k = 1:nt
  % du/dt = -cdu/dx
  % du/dx = BDS
  % du/dt = FDS
  % Space loop
  for i = 2:n
     u(i) = uold(i) - (c*dt/dx)*(uold(i) - uold(i-1)); 
  end
  % updating velocity
  uold = u;
  
end
end

Main program:

In the main program input parametres are declared outside the function and then executed with the array of dt values. Then the elapsed time is recorded.

close all 
clear all 
clc

% Input parameters
L = 1;                       % length in meter
n = 80;                      % no of grid points
c = 1;                       % velocity
t = 0.4;
dt = [1e-4,1e-3,1e-2,1e-1];  % time step width
dx = L/(n-1);                % dt = t/nt; 
CFL = abs(c)*dt/dx;          % CFL calculation formula
x = linspace(0,L,n);         % Declaring x values with lispace command 

tic;                         % Time duration recording command 'tic toc' 
[x1,u1] = linearconvection_func(L,n,c,t,1e-4,x);   % calling Function 
ElapsedTime_1 = toc;
tic
[x2,u2] = linearconvection_func(L,n,c,t,1e-3,x);
ElapsedTime_2 = toc;
tic
[x3,u3] = linearconvection_func(L,n,c,t,1e-2,x);
ElapsedTime_3 = toc;
tic
[x4,u4] = linearconvection_func(L,n,c,t,1e-1,x);
ElapsedTime_4 = toc;

% plotting 
 plot(x1,u1,'r')
 hold on 
 plot(x2,u2,'g')
 plot(x3,u3,'b')
 plot(x4,u4,'m')
 grid on
 legend('1e-4','1e-3','1e-2','1e-1')
 title({'1-D Linear Convection',['time-steps(\itt) = ',num2str(dt)]})
 xlabel('Spatial co-ordinate (x) \rightarrow')
 ylabel('Velocity (u)')
 axis([0 1 0 2.5])
  

Plotting:

Plotting is done with differnt colours for each dt value for easy undersatnding. Here we can see all the four velocity profile.

The blue color plot which is for 1e-2 value of dt has the stable solution, maintianing the constant velocity 2m/s.

The most unstable solution is plotted in pink color, in which velocity profile is having blowup. The CFL value for this exceed than 1.

Now , let us see the workspace values

From the above workspace value, we can see that the elapsed time is more for dt =1e-4 and very less for dt = 1e-1.

Thus we can conclude this by saying if the CFL value is less, it will take more elapsed time that means we will get the stable solution with some sort of diffusion and vice versa.