## 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
clc

% 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)
tic
% 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));
end
% update old velocities
u_old = u;   t_exe(1)=toc
end

% 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;
tic
% 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));
end
% update old velocities
u1_old = u1;   t_exe(2)=toc
end

% 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;
tic
% 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));
end
% update old velocities
u2_old = u2;        t_exe(3)=toc
end

% 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;
tic
% 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));
end
% update old velocities
u3_old = u3;          t_exe(4)=toc
end
% plotting the velocities
figure(1)
plot(x,u,x,u1,x,u2,x,u3,'Linewidth',2)
grid on
axis([0 1 0 3])
xlabel('Domain Length')
ylabel('velocity')
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')

figure(2)
z = categorical({'1e-4', '1e-3', '1e-2', '1e-1'})
bar(z,t_exe)
grid on
xlabel('Timestep(dt)')
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.

### Parametric Study of Gate Valve Shravankumar Nagapuri · 2019-10-22 03:14:10

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

### Symmetry vs Wedge vs HP equation Shravankumar Nagapuri · 2019-10-09 22:50:52

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

### Combustion simulation on the combustor model Shravankumar Nagapuri · 2019-10-05 16:23:22

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 challenge Shravankumar Nagapuri · 2019-10-01 17:50:02

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

### Gearbox sloshing effect Shravankumar Nagapuri · 2019-09-28 15:37:58

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

### Simulation of Flow through a pipe in OpenFoam Shravankumar Nagapuri · 2019-09-26 19:10:15

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

### Conjugate Heat Transfer CHT Analysis on a graphics card Shravankumar Nagapuri · 2019-09-20 16:08:56

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 FVM Literature Review Shravankumar Nagapuri · 2019-09-17 13:14:28

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

### Rayleigh Taylor instability Shravankumar Nagapuri · 2019-09-14 17:48:30

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

### BlockMesh Drill down challenge Shravankumar Nagapuri · 2019-09-11 18:17:33

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