1D Linear Convection in Matlab using first-order discretization schemes and the effect of time-step on the stability of the solution

INTRODUCTION: To solve the 1D linear equation to observe the propagation of wave with time.

Governing Equation

`{partial u}/{partial t} + C{partial u}/{partial x} = 0`

C = Velocity of Propagation

The above equation is discretized using first order forward differencing scheme for time derivative and backward differencing scheme for space derivative. 

The explicit numerical solution for the above equation

 `u_i^(n+1) = u_i^n - {C ∆t}/{∆x}(u_i^n- u_(i-1)^n)`

where 'n' is time step

'i' is node no.

To define a domain by using some inputs

L = 1 m

n = 80

C =1

set time step, dt = [1e-4, 1e-3, 1e-2, 1e-1]

u(t=0)= 2 m/s if `0.1m lexge0.3`

u(t=0)= 2 m/s if otherwise

The Matlab code for linear convection:

The above numerical solution into Matlab, The programme outputs the grid, the initial velocity profile and final velocity profile at t = 0.4 seconds. The code is shown below.

clear all
close all

% Inputs for the linear convection problem
       % Length of the domain in 'm'    
           L = 1;
       % Number of grid points
           n = 80;
       % Linear convection constant
           C = 1
       % Time step in seconds
          dt = [1e-4 1e-3 1e-2 1e-1]
       %  End time
          t_end = 0.4;
% calculations
       % x array along the length of the domain 
          x = linspace(0,L,n);
       % Grid spacing, dx = L/(n-1) or dx =x(2)-x(1)
          dx = x(2) - x(1);
       % time 
             nt1 = t_end/dt(1); nt2 = t_end/dt(2)
             nt3 = t_end/dt(3); nt4 = t_end/dt(4)  
       % Define the start and end positions for the wave
              x_start = 0.1;    x_end = 0.3;
       % n_start = f(x_start)
              n_start = 0.1/dx +1
       % n_end = f(x_end)
              n_end  = 0.3/dx +1
        % Initialize arrays and apply boundary conditions
                u = ones(1,n); u(n_start:n_end) = 2  
        % u(1) = 1;
                u_old = u;    u_initial = u;

        % Time of Execution
                 t_exe = zeros(1,4)
% Time loop(for time step 1e-4, and t_end = 0.4 )
   for k = 1:nt1
         % space loop
             for i = 2:n
                 u(i) = u_old(i) - (C*dt(1)/dx)*(u_old(i)- u_old(i-1));
         % update old velocities
             u_old = u;   t_exe(1)=toc

         % Initialize arrays and apply boundary conditions
              u1 = ones(1,n);    u1(n_start:n_end) = 2  
         % u1(1) = 1;
              u1_old = u1;     u1_initial = u1;
% Time loop(for time step 1e-3, and t_end = 0.4 )
   for k = 1:nt2
        % space loop
        for i = 2:n
            u1(i) = u1_old(i) - (C*dt(2)/dx)*(u1_old(i)- u1_old(i-1));
        % update old velocities
             u1_old = u1;   t_exe(2)=toc

        % Initialize arrays and apply boundary conditions
              u2 = ones(1,n);     u2(n_start:n_end) = 2  
        % u2(1) = 1;
             u2_old = u2;      u2_initial = u2;
% Time loop(for time step 1e-2, and t_end = 0.4 )
   for k = 1:nt3
        % space loop
          for i = 2:n
             u2(i) = u2_old(i) - (C*dt(3)/dx)*(u2_old(i)- u2_old(i-1));
          % update old velocities
             u2_old = u2;        t_exe(3)=toc

          % Initialize arrays and apply boundary conditions
              u3 = ones(1,n);    u3(n_start:n_end) = 2  
          % u3(1) = 1;
              u3_old = u3;    u3_initial = u3;
% Time loop(for time step 1e-1, and t_end = 0.4 )
    for k = 1:nt4
        % space loop
          for i = 2:n
            u3(i) = u3_old(i) - (C*dt(4)/dx)*(u3_old(i)- u3_old(i-1));
         % update old velocities
            u3_old = u3;          t_exe(4)=toc
% plotting the velocities
grid on
axis([0 1 0 3])
xlabel('Domain Length')
lgd = legend(['dt=', num2str(dt(1))], ['dt=', num2str(dt(2))], ['dt=', num2str(dt(3))],['dt=', num2str(dt(4))]);
title('comparion of outputs for different time steps for T =0.4')

z = categorical({'1e-4', '1e-3', '1e-2', '1e-1'})
grid on
ylabel('Time of Execution')
title('comparion of time taken to solve 1D Linear convective Equation')


Plot 1: 



Plot 2:

Comparison and Conclusion:

1.  From the above plots, it is clear that the solution for dt =0.1 is  unstable, where as for other values of dt the solution is stable. This is because from the numerical formulation of the solution if we take K= `{C ∆t}/{∆x}`, the solution modifies to:

 `u_i^(n+1) = u_i^n - K(u_i^n -u_(i-1)^n)` = `(1-K)u_i^n + u_(i-1)^n`

The (1-K) part on the RHS affects the stability of numerical discretization. If K > 1 the solution diverges to instability. The stability and accuracy of the solution increases as the value of K approaches to unity. For dt = 0.01 , the solution is accurate and smoother. If the solution is unstable, then it can be analyzed by the von Neumann of stable analysis.

2. The Explicit solution is preferred over the implicit, because of the  cost of computation time and cost in order to attain stability. The accuracy of the implicit solution is to be compromised over the explicit solution.

3. The time taken for each value of dt also compared in bar graph, with the decrease of dt the number of time steps required to reach same end time increases, which inturn increase the number of iterations and time taken for solving the equation. 





Projects by Shravankumar Nagapuri

Objectives: To perform a parametric study on the gate valve simulation. To obtain the mass flow rates at the outlet for 5 design points. To show the cut section view for different design points, and also to show the gate disc lift and fluid volume extraction. To s Read more

Objectives: To Write a Matlab program that can generate the computational blockMesh file  automatically for a wedge angle of 3 degrees and compare the results obtained for Symmetry BC and Wedge BC. To Write a Matlab program that takes an angle as input and gene Read more

1Q. Briefly explain about the possible types of combustion simulations in FLUENT. A:  Combustion or burning is a high temperature exothermic chemical reaction between a fuel and an oxidant accompanied by the production of heat, light and unburnt gases in the fo Read more

Cyclone Separator: Principle and Working: A high speed rotating (air)flow is established within a cylindrical or conical container called a cyclone. Air flows in a spiral pattern, beginning at the top (wide end) of the cyclone and ending at the bottom (narrow) end bef Read more

Objective: To analyse the flow pattern of the fluid inside the gearbox, for two different clearances of the same geometry and each geometry flow is analyzed by two different fluids (Water and Oil) as lubricants. Gearbox sloshing effect:  Slosh refers Read more

Objectives: To write a program in Matlab that can generate the computational mesh automatically for any wedge angle and Grading scheme. To calculate, length of the pipe using the entry length formula for laminar flow through a pipe To show that entry length is suf Read more

Objectives: To find out the maximum temperature attained by the processor. To Prove that the simulation has achieved convergence with appropriate images and plots. To Find out the heat transfer coefficient at appropriate areas of the model. To Identify potential h Read more

Finite volume Method: The finite volume method (FVM) is a method for representing and evaluating partial differential equations in the form of algebraic equations. Similar to the finite difference method (FDM) or finite element method (FEM), values are calculated at Read more

Q1. What are some practical CFD models that have been based on the mathematical analysis of Rayleigh Taylor waves? In your own words, explain how these mathematical models have been adapted for CFD calculations  Kelvin–Helmholtz instability: The Kelvin Read more

Objectives: To Generate the BlockMesh file for the given Geometry using OpenFOAM and use the icoFoam solver to simulate the flow through a backward-facing step. To create multiple meshes and will be comparing the results obtained from each mesh. To show the velocit Read more


The End