Numerical solution to the wave equation - Explicit Method.

• I am going to write a program in MATLAB which will compare initial and final velocity profile for 1D Linear convection for different value of grid points.

1D linear convection is described as follows:

`ubrace((delu)/(delt) + C(delu)/(delx) = 0)_("Equation - I")`

Problem setup:
1. The domain length is set to ` 1m`
2. The initial velocity profile is a step function. It is equal to `2ms^-1` between `x= 0.1` and `x = 0.3` and `1ms^-1` elsewhere.
3. To discretize this equation in space and time we will use the Forward Difference Scheme (FDS) for time derivative and Backward Difference Scheme (BDS) for space derivative.
4. The speed of the wave '`C`' is set constant to `1 ms^-1`
5. The time step '`Deltat`' is set to `0.01`
6. Total amount of time '`T`' for analysis is set to ` 0.4s`


I will now discritise `eq^n - I`: `obrace((u_i^(n+1) - u_i^n)/(Deltat))^("FDS") + C*obrace((u_i^n - u_(i-1)^n)/(Deltax))^("BDS") = 0`

Transpose the above equation: `u_i^(n+1) = u_i^n - C(Deltat)/(Deltax)*(u_i^n - u_(i-1)^n)`

Here, `i rarr` index of a grid in x and `n rarr` index of a grid in t.

Now, I will write a program in MATLAB to solve this algebraic equation and generate the plots for different value of the grid points used, namely 20, 60, 80, and 160 using marching method.

PROGRAM:

Here is how the function file 'plot_convection.m' looks like:

function plot_convection(n, ax)  %n = grid points, ax = axis for plots 

%%Writing down the assumptions and applying Boundary conditions to initialize the velocity profile
L = 1;  % Length of the domain
C = 1;  % Constant speed for convection
T = 0.4;  % Total time
nt = 40;  % Iterations set to 40 for time step to be 0.01
dt = T/nt;  % Time step, it can be calculated as T/(nt-1)
x = linspace(0,L,n);  % Creating 1D grid in x direction
dx = L/(n-1);  % Difference between two consecutive grid points
u = ones(1,n);  % Creating velocity array with all velocity values = 1 including BC
%Using the formula for nth term of AP, I derived the following two results
n_start = 0.1/dx + 1;  % Starting node of the square wave
n_stop = 0.3/dx + 1;  % Ending node of the square wave
%Using the above results obtained
u(n_start:n_stop) = 2;  % Changing velocity values from 1 to 2 for the nodes where square wave will start and end
u_old = u;  %Previous node velocity
u_initial = u;  % Initial velocity
t = 0;  %Starting time

%%Time loop - Time marching solution
for i = 1:nt
     %Space loop - Space marching
     for j = 2:n
          u(j) = u_old(j) - (C*dt/dx)*(u_old(j) - u_old(j-1));
     end
     %Updating velocity
     u_old = u;
     t = t + dt;
     %Displaying the plot generated
     title(sprintf('For %g grid points and %g Iterations\n Time = %g, dt = 0.01 and dx = %g', (n-1), nt, t, dx),'Fontsize',14)
     grid on
     xlabel('1D Grid','Fontsize',12)
     ylabel('Velocity','Fontsize',12)
     h = plot(x,u_initial,'-bo','linewidth',2,'markersize',6,'markeredgecolor','k','markerfacecolor','y');  % Displaying all iterations
     hold on 
     h1 = plot(x,u,'-r','linewidth',1);  % Plots the initial velocity profile
     axis(ax)
     pause(0.03)
     set(gca,'XTick',(0:0.1:1));       
end    
 
h2 = plot(x,u,'-ko','linewidth',2,'markersize',6,'markeredgecolor','b','markerfacecolor','c');  % Plotting the final velocity profile
H = legend([h h1 h2],{'Initial Velocity Profile','Evolution Of Velocity Profiles','Final Velocity Profile'});
set(H, 'FontSize', 11)

CFL = 0.01/dx;  %CFL = C*dt/dx it can help us knowing the stability.

fprintf('The CFL number for %g grid points is %g. ', (n-1), CFL)

if CFL < 1 
   disp('Thus, stable.')
elseif CFL == 1
   disp('Thus most stable.')
else
   disp('Thus, unstable.')
end  

Here is how the driver script 'explicit_oneD_LC.m' look like:

% Explicit Scheme - 1D Linear Convection

clear all
close all
clc

%%Creating four figure windows to display our plots generated by calling a function - plot_convection
figure(1)  % For n = 21 i.e 20 points, and dx = 0.05
plot_convection(21, [0 1 0.8 2.4])
figure(2)  % For n = 41 i.e 40 points, and dx = 0.025
plot_convection(41, [0 1 0.8 2.4])
figure(3)  % For n = 81 i.e 80 points, and dx = 0.0125
plot_convection(81,[0 1 0.8 2.4])
figure(4)  % For n = 161 i.e 160 points, and dx = 0.00625
plot_convection(161,'tight')
figure(5)  % For n = 101 i.e 100 points, and dx = 0.01 - Optimum stability
plot_convection(101,[0 1 0.8 2.4])

Results:

Plots generated:
Note: I have plotted all steps `i.e` 40 iterations in one graph so that we can observe the evolution of the velocity profile.

Plot 1:

Plot 2:

Plot 3:

Plot 4:

Explaining the plots:

Plot 1: The first order backward differencing scheme for space derivative creates an artificial diffusion. Thus, we can see that the profile is losing its square shape and is going flatter. Another factor for losing the shape is the size of a step size.

Plot 2: The peak is closer to `2` if you compare it with the above plot. The reason is simple, as the value of step size is reduced, the error has reduced too. Thus, we can also see that the profile is getting sharper than the above plot.

Plot 3: The peak is pretty much closer to `2` than the above two plots. As the dx value has decreased, even more, we can observe that the profile has started gaining it's square shape back.

Plot 4: There is extreme disturbance seen in the velocity profile, the shape is lost. This is because on the basis of the assumptions we had made, the value of dx for grid points of 161 `i.e` 0.00625 is very very small and it makes the solution go unstable.

Conclusions:

1. As the value of `C > 0`, the velocity profile is convected to the right with a constant speed of `C = 1ms^-1`.
2. According to the Courant-Friedrichs-Lewy (CFL) stability criterion `i.e` `C(Deltat)/(Deltax) <= 1 <=> C<= (Deltax)/(Deltat)` , the best solution obtained is for `100` grid points used.
   The plot looks like this:
   
   I observed that neither the profile has lost its shape nor the peak has decreased. I could achieve this perfect solution just because of the assumption `C = 1ms^-1` made.
   

The CFL number '`C(Deltat)/(Deltax)`' and stability for respective grid points used above was printed on the command window:

The CFL number for 20 grid points is 0.2. Thus, stable.
The CFL number for 40 grid points is 0.4. Thus, stable.
The CFL number for 80 grid points is 0.8. Thus, stable.
The CFL number for 160 grid points is 1.6. Thus, unstable.
The CFL number for 100 grid points is 1. Thus most stable.