Numerical analysis of 1-D Linear Convection equation using MATLAB to show the effect of Time Step Size (dt) on the solution

Numerical analysis of 1-D Linear Convection equation using MATLAB to show the effect of Time Step Size (dt) on the solution

PROBLEM SPECIFICATION:-

Assumptions:

1. Domain Length, L=1m.
2. The initial velocity profile is a step function. It is equal to 2m/s between x=0.1 & 0.3, and 1m/s everywhere else.
3. The number of grid points, n=80.
4. Using first-order forward differencing for the time derivative.
5. Using first-order backward/rearward differencing for the space term.
6. The Timestep and the number of grid points should be variable in this calculation.
7. Constant velocity, C=1m/s.
8. Step time step=0.01.

Expected Output:

1. In a single graph compare the numerical solution at t = 0.4 seconds.
2. Make the comparison for, time step, dt=(1e-4, 1e-3, 1e-2, 1e-1).
3. Plot the simulation time as a function of the time step.
4. Write a detailed explanation for your results, plots, or codes.

------------------------------------------------------------------------------------------------------------------------------------------------------------------

THEORY:-

Navier-Stokes x-momentum equation in 2-D:

`(∂u)/(∂t)+u (∂u)/(∂x)+v (∂u)/(∂y)=-(∂p)/(∂x)+μ[(∂^2 u)/(∂x^2 )+(∂^2 u)/(∂y^2 )]`

In CFD most of the flows can be represented by a Navier Stokes Equations. The Navier Stokes Equation is a Combination of convection, pressure drop, diffusion terms, body force, and source terms.

1-D Linear Convection:

This equation is the most accessible in CFD; from the Navier Stokes equation we kept only the accumulation and convection term for the `x` component of the velocity- as we already know, in CFD the variables to be computed are velocities; to make things even simpler, the coefficient of the first derivative of the velocity is constant, making the equation linear.

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

This equation represents the propagation of a wave without change of shape; the speed of the wave is `c` and with the initial condition given by:

`u(x,0)=u_0 x`

The exact solution is:

`u(x,t)=u_0 (x-c.t)`

This means that on an axis with a positive direction to the right, the wave will move from left to right with the speed of `c`.

To discretize this equation in space and time we’ll use the ‘Forward Difference Scheme’ for the time derivative and ‘Backward Difference Scheme’ for space derivative.

The partial derivative can be approximated using the derivative definition, removing the limit operator:-

Forward differencing, `(∂u)/(∂x)≈(u(x+∆x)-u(x))/(∆x)` and Backward differencing, `(∂u)/(∂x)≈(u(x)-u(x-∆x))/(∆x)`

This approach will be applied both for space and time; for time steps we’ll use the superscript `n` and for space steps (grid cells) we’ll use the subscript `i`.

Here, we are using the Explicit approach to get the numerical solution of a one-dimensional linear convection equation. The major difference between explicit and implicit is that, in Implicit approach; when we’re solving, on the right-hand side we take `(n+1)` terms and we get a set of coupled linear equations. Whereas in the Explicit approach, the right-hand side equation contains the term of`(n)` level; and we get uncoupled linear equations.

Explicit form, `(u_i^(n+1)-u_i^n)/(∆t)=-c.(u_i^n-u_(i-1)^n)/(∆x)` and Implicit form, `(u_i^(n+1)-u_i^n)/(∆t)=-c.(u_i^(n+1)-u_(i-1)^(n+1))/(∆x)`

Before writing the code some clarifications are needed for the problem statement; the domain for this problem will be a length of 1m; the problem is solved on the `x`-axis. The initialization is made with a velocity of 1m/s all over the grid point, except those that lie between `x`=0.1 and 0.3; on these points velocity is 2m/s. We’re going to write a function that accepts the time step as an argument and solves the same problem. Here we’ll set our number of grid points to, `n=80` and to show the effect of the time step size on the solution; we’re taking `dt`=1e-4, 1e-3, 1e-2, & 1e-1.

MATLAB code used to solve ‘1-D Linear Convection Equation’:-

1. Creating a function that accepts the time step as an argument:

%Creating a function for time step:-
function [u,x]=time_step_func_1D_linear_convection(dt)
 
%Inputs for the Linear convection problem;
 
%Length of domain,L=1m
L=1;
%No. of grid points,n=[]20 40 80 160]
n=80;
%Linear convection veloctiy,c=1m/s
c=1;
%total time, t=0.4sec
t=0.4;
%time step,dt=[1e-4 1e-3 1e-2 1e-1]
%no. of time steps, nt=t/dt
nt=t/dt;
 
%creating meshgrid
x=linspace(0,L,n);
%grid spacing
dx=L/(n-1); 
%Initializing arrays and boundary conditions
u=1*ones(1,n);
%defining the start and end positions for the square wave
      for j=1:length(x)
          x_start=0.1;
          x_stop=0.3;
          if (x(j)>=x_start && x(j)<=x_stop)
              u(j)=2;
          else
          u(j)=1;
          end  
      end
    uold=u;
    u_initial=u;
  
%time marching solution
%time loop
for k=1:nt
    %space loop
    for i=2:n
        u(i)=uold(i)-(c*dt/dx)*(uold(i)-uold(i-1)); %using explicit approach
    end
%updating old velocity
uold=u;
 
end %time loop ends
 

2. Main program to show the effect of time step size (dt) on the solution:

%Effect of time step size on the solution
 
%Counting the time required to process for different time steps (dt=1e-4, 1e-3, 1e-2, and 1e-1)
 
tic
[u_1,x_1]=time_step_func_1D_linear_convection(1e-4);
toc
time1=toc;
 
tic
[u_2,x_2]=time_step_func_1D_linear_convection(1e-3);
toc
time2=toc;
 
tic
[u_3,x_3]=time_step_func_1D_linear_convection(1e-2);
toc
time3=toc;
 
tic
[u_4,x_4]=time_step_func_1D_linear_convection(1e-1);
toc
time4=toc;
 
%Plotting graph for velocity vs. different time step
 
figure (1)
plot (x_1,u_1,'b','linewidth',2,'linestyle','--')
hold on 
grid on
plot (x_2,u_2,'r','linewidth',2)
hold on
plot (x_3,u_3,'g','linewidth',2)
hold on
plot (x_4,u_4,'c','linewidth',2)
hold on
axis([0 1 0 2.5])
xlabel ('Length of domain, L (m)')
ylabel ('Velocity, u (m/s)')
legend ('time step=1e-4','time step=1e-3','time step=1e-2','time step=1e-1')
title ('Comparison of Velocity profile for various time step')
 
%plotting bar graph for 'simulation time as a function of time-step'
 
figure(2)
simulation_time=[time1,time2,time3,time4];
bar(simulation_time,'g')
set(gca,'XTicklabel',{'1e-4','1e-3','1e-2','1e-1'})
grid on
xlabel('Time Step (dt)')
ylabel('Simulation Time (in sec)')
title('Simulation Time as a function of the Time-Step')
 
%plotting log-plot for 'simulation time as a function of time-step'
 
figure (3)
loglog (1e-4,time1,'marker','*','markersize',20,'linewidth',3)
hold on 
grid on
loglog (1e-3,time2,'marker','+','markersize',20,'linewidth',3)
hold on
loglog (1e-2,time3,'marker','x','markersize',20,'linewidth',3)
hold on
loglog (1e-1,time4,'marker','o','markersize',20,'linewidth',3)
hold on
loglog ([1e-4 1e-3 1e-2 1e-1],[time1 time2 time3 time4])
xlabel('Time Step (dt)')
ylabel('Simulation Time (in sec)')
title('Simulation Time as a function of the Time-Step')

RESULTS & DISCUSSIONS:-

From the above plot, we can observe that the solution for time step ‘1e-3’ and ‘1e-4’ has the optimum velocity nearly equal to 1.8m/s. Here we can see that the artificial diffusion error is greater as compared to time step ‘1e-2’.

Whereas, when we see the velocity profile for the time step ‘1e-2’ it has the optimum velocity approximately equal to 1.9m/s; which is nearly equal to the original velocity. Hence, the least artificial diffusion error can be seen for the time step ‘1e-2’.

Further, if we observe the plot for the time step ‘1e-1’, it can be observed that the solution is completely unstable and blown off. This is occurred due to exceeding the time step limit for an explicit approach. It is always a viable approach to use Neumann stability analysis for the explicit method.

If we observe and compare the time required to calculate the solution for different time steps, we can observe that for the time step ‘1e-2’, simulation time is the smallest and this particular time step gives the most immaculate solution. Also, we can see that for time step ‘1e-4’, the simulation time is the largest; with the fact that it is the smallest amongst all-time steps taken in this problem. It can be inferred from this behavior that for an explicit approach; the smaller the time step, the larger will be the simulation time for a complete solution.

CONCLUSION

We can generalize from our results and discussions that, as the time step increases the artificial diffusion error in the solution decreases. But when the time step increased beyond the stability criterion, the solution becomes unstable.

So, it can be also concluded that making the time step too small makes the simulation run for a long time and as a result; diffusion error increases with each iteration and the final solution has a larger diffusion error. It would be better if we check the stability of finite difference schemes before performing any time-dependent problems because the stability of numerical schemes is closely related to errors.