Numerical simulation of Linear Convection at different grid points using MATLAB

Aim

Study of linear convection using numerical techniques at various grid points.

Governing equation

The general formula for 1D convection is given by

Discretizing the above equation in terms of explicit scheme where the first term is discretized with respect to time and the second term with respect to space gives:

N represents the present time step.

U represents velocity.

Boundary conditions and numerical parameters

1. Length of domain L=1m
2. Initial velocity profile is a step function

V=2m/s, 0.1<x<0.3

V=1m/s, everywhere

3. C=1
4. Time steps =0.01
5. U(1)=1m/s
6. Forward differencing for time derivative and backward differencing for space terms

Algorithm

• Define the constant variables given in the data.
• Define two function call statements which calculates the nodal position where the velocity profile changes to 0.3m/s.
• Define the no of grid points and the boundary conditions
• Define the no of time steps to run the loop (time marching).
• Inside the time marching loop define the space loop which calculates velocity in the next time interval.
• Plot graph between the initial and final velocity profile wrt grid spacing.
• Keeping the time steps to be constant change the no of grid points and run step 5 and 6.

the program code for the above algorithm is given below

Matlab code

main code

clc
clear all
close all
%%%%%solving linear convection problem
%%%%%du/dt=-c*du/dx
%length of the domain
l=1;
c=1; 
%setting time steps
dt=1e-2;
%define the grid points
n=20;% no of grid points is calculated for n=20,40,80,160
%define the mesh
x=linspace(0,l,n);
dx=x(2)-x(1);%grid spacing
%functioncall to calculate the starting node 
n_start=finding_starting_node(x,dx);
%functioncall to calculate the end node
n_end=finding_end_node(x,dx);
%define velocity at the node points
u=ones(1,n);
u(n_start:n_end)=2;
uold=u;
u_initial=u;
%bc
u(1)=1;
%no of time steps
nt=50;
tic
%%%%%%%using forward difference wrt time and backward difference wrt space
for k=1:nt%time marching steps
for i=2:n% space convergence loop
    
u(i)=uold(i)+(-c*dt/dx)*(uold(i)-uold(i-1));
end
%update velocities
uold=u;
figure(1)
plot(x,u,'b')%initial velocity profile
grid on
hold on 
plot(x,u_initial,'color','r','linewidth',2)%final velocity profile
xlabel('grid spacing')
ylabel('velocity')
axis([0 1.5 0 2])
legend('final velocity profile','initial velocity profile')
title('linear convection problem at timestep 1e-4')
end
time=toc;%compute time for simulation

function call for starting node

function output=finding_starting_node(x,dx)
output=1;
for i=1:length(x)
if x(i)<= 0.1
    output=output+1;
end 
end
end

function call for end node

function output=finding_end_node(x,dx)
output=1;
for i=1:length(x)
if x(i)<= 0.3
    output=output+1;
end 
end
end

Results and analysis

For n=20

The velocity profile initially represented as a square velocity profile. Considering the grid points equal to 20, the velocity profile turns to be sinusoidal wave. The peak velocity drops from 2m/s to 1.6m/s as it propagates. The change in velocity profile is highly distorted.

For n=40

The final velocity profile for grid points 40 also turns to be sinusoidal as that of n=20. The peak velocity profile is reduced to 1.8m/s, which is higher than previous final velocity profile(1.6m/s). here the propagation marks more stable solution than former profile.

For n=80

On considering no of grid points to be 80, we see that final profile is same as that obtained from n=20 and n=40. The final peak profile is same as that of its initial profile 2m/s. There is a negligible drop in velocity magnitude. Here the propagation of the solution is more stable than former two.

For n=160

At grid point 160, we see that the entire velocity profile breaks down and discontinuities arises which spoils the entire solution. The velocity profile abruptly changes its magnitude implying no stable solution. Here the defined grid points have exceeded the threshold value of the solution.

Conclusion

On varying the no of grid points, the drop of peak velocity is maximum for lowest no of grid points 20 and minimum for grid points 80. It is seen that at n=80 the solution is most stable and for n=160 it is most unstable. The change in the velocity profile is due to the existence of a numerical diffusion term. This accounts for the relative error at different no of grid points.