Numerical analysis of 1-D Linear Convection equation for varying grid points (n) using MATLAB
Numerical analysis of 1-D Linear Convection equation for varying grid points (n) using MATLAB
PROBLEM SPECIFICATION:-
Assumptions:
- Domain Length, L=1m.
- 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.
- Using first-order forward differencing for the time derivative.
- Using first-order backward/rearward differencing for the space term.
- The Timestep and the number of grid points should be variable in this calculation.
- Constant velocity, C=1m/s
- Step time step=0.01
Expected Output:
- For C=1 and for time=0.4 second, compare the original and final velocity profiles.
- Make the comparison for n=20, 40, 80, and 160.
- Generate the four plots (each for n=20, 40, 80, and 160).
- 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.
MATLAB code used to solve ‘1-D Linear Convection Equation’:-
- For c=1 and time =0.4sec: Plotting the comparison plot for initial velocity vs. final velocity (n=20, 40, 80, 160)
clear all close all clc %Inputs for the Linear convection problem %Length of domain,L=1m L=1; %No. of grid points,n=[20 40 80 160] n=[20 40 80 160]; %Linear convection veloctiy,c=1m/s c=1; %time step,dt=0.01 dt=1e-2; %no. of time steps, nt=40 nt=40; %grid loop for m=1:length(n) x=linspace(0,L,n(m)); %creating meshgrid dx=L/(n(m)-1); %grid spacing %Initializing arrays and boundary conditions u=1*ones(1,n(m)); %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(m) 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 %plotting the velocities figure(m) plot(x,u,'g') hold on grid on plot(x,u_initial,'r') xlabel('Length of the domain(m)') ylabel('Velocity (m/s)') legend('Final velocity','Initial velocity') axis([0 1 0 2.5]) pause(0.4) end %grid loop ends - For c=1 and time=0.4sec: plotting the velocity trend for different no. of grid points(n=20, 40, 80, 160)
clear all close all clc %Inputs for the Linear convection problem; %Length of domain,L=1m L=1; %No. of grid points,n=[20 40 80 160] n=[20 40 80 160]; %Linear convection veloctiy,c=1m/s c=1; %time step,dt=0.01 dt=1e-2; %no. of time steps, nt=40 nt=40; z=1:length(n); %grid loop for m=1:length(n) x=linspace(0,L,n(m)); %creating meshgrid dx=L/(n(m)-1); %grid spacing %Initializing arrays and boundary conditions u=1*ones(1,n(m)); %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(m) u(i)=uold(i)-(c*dt/dx)*(uold(i)-uold(i-1)); %using explicit approach end %updating old velocity uold=u; %plotting the velocities figure(m) plot(x,u,'b') hold on grid on plot(x,u_initial,'r') xlabel("Length of the domain(m)") ylabel("Velocity (m/s)") axis([0 1 0 2.5]) pause(0.04) legend('Final velocity Profile','Initial velocity Profile') grid on end %time loop ends end %grid loop ends
RESULTS & DISCUSSIONS:-
1. Comparison plot for initial velocity vs. final velocity:
For no. of grid points, n=20
When the number of grid points is 20, we can see from the above graph that, the final velocity loses its sharpness as compared to the initial velocity. Velocity reduces along with the domain due to the artificial diffusion error term; which is higher for this case because of the low number of grid points.
b) For no. of grid points, n=40
As the number of grid points increased to 40, our mesh grid becomes fine and we can notice an improvement in final velocity. It can be inferred from the above graph that, as we increase the number of grid points there’s a decrease in error term (artificial diffusion).
c) For no. of grid points, n=80
When we further increase our number of grid points to 80, our mesh becomes finer and the final velocity approximates to the initial velocity in a good manner; but there are still some error terms (artificial diffusion).
d) For no. of grid points, n=160
As the number of grid points increased to 160, our velocity profile is blown off completely and the solution becomes unstable. It can be concluded from the above graph that, we can increase the number of grid points at a certain level to keep the stability of the solution.
2. Velocity trend plot for different no. of grid points(n=20, 40, 80, 160):
a) For no. of grid points, n=20
From the above graph, it can be observed that; as time increases the initial velocity profile keeps moving towards the right, but as it moves: it is losing sharpness and the peak velocity is decreasing substantially. The decrement in the velocity profile is due to the error generated due to the artificial diffusion. The diffusion term is the effect of the truncation error that we will get while solving the differential equation.
b) For no. of grid points, n=40
As we make our meshing finer there will be less truncation error, less truncation error results in less numerical dissipation or diffusion and we get the more stable solutions than before. This effect can be seen in the above plot.
c) For no. of grid points, n=80
Here we have got the mesh where we can say the solution is near complete stability. The diffusion is comparatively less and the peak velocity is the same as peak velocity for the initial velocity profile.
d) For no. of grid points, n=160
In the above plot, we can see that for the given mesh (n=160) the solution is blown up and the solution gets substantially unstable.
CONCLUSION:-
It can be concluded that the accuracy of the solution increases, as the number of grid points increases or mesh size decreases. Furthermore, we can only increase the number of grid points to a certain threshold; after the threshold, we get an unstable solution.
More the truncation error more will be the numerical dissipation or artificial diffusion; more truncation error leads to the unstable system. It is important to know which differencing scheme is to for a particular system. By using Von Neumann's stability analysis we can figure out when a particular numerical scheme is going to be unstable.