4.1 - 1D First Order Wave Equation

In this project, the 1D Linear wave equation is solved numerically, using FWD and BWD and a variable time step and grid size.

This mathematical approximation was implemented in Matlab as follows.

function FirstOrderWaveEquation(n,dt)
%Inputs for the linear convection problem
% n=21;
% dt=0.01;

L=1;
c=1;
t=0.4;
nt=t/dt;

%Mesh Calculations
x=linspace(0,L,n);
dx=x(2)-x(1); % Or dx=L/(n-1);

%Initialize the square wave and BC
x_start=0.1;
x_end=0.3;
n_start= find(abs(x-x_start)<dx/2);
n_end=find(abs(x-x_end)<dx/2);
u=ones(1,n);
u(n_start:n_end)=2;

%BC
u(1)=1;
u_old= u;
u_t0=u;

%Marching Solution Loops

%Time loop
for k=1:nt
    
    
%Space loop (single time step) BDS  NOT STABLE!! Oscillations keep growing.

    CFL=c*dt/(dx);
    for i=2:n
        u(i)=u_old(i)- CFL*(u_old(i)-u_old(i-1)); 
    end 
    %Update the old velocities
    u_old=u; 
    
    %Plotting the velocities
    title(sprintf(\'Grid Points:%g,Time step:%g s\', (n-1), dt),\'Fontsize\',12)
    graph1= plot(x,u,\'b\');
    axis([0 1 0.8 2.2])
    ylabel(\'Velocity\')
    xlabel(\'X-direction Grid\')
    grid on
    hold on
 
   
end
    graph0= plot(x,u_t0,\'-*r\',\'LineWidth\',2);
    hold on
    graph2= plot(x,u, \'-ok\',\'LineWidth\',2);
    legend ([graph0, graph1, graph2],{\'Initial Velocity\',\'Profile Evolution\',\'Final Velocity\'})

 
end 
    

This function uses as imputs the number of nodes and time step in the model. In order to compare all the conditions required, the following simple main program was written.

clear all
close all
clc

%For n=20 nodes and dt=0.01
figure(1)
FirstOrderWaveEquation(21,0.01)
%For n=40 nodes and dt=0.01
figure(2)
FirstOrderWaveEquation(41,0.01)
%For n=80 nodes and dt=0.01
figure(3)
FirstOrderWaveEquation(81,0.01)
%For n=160 nodes and dt=0.01
figure(4)
FirstOrderWaveEquation(161,0.01)

The output of the code is:

Figure 1) 20 Grid Points

First of all, notice that the initial profile is not a perfect square wave due to the finite number of grid points.

It can be clearly seen that the profile loses its squared shape, and that the velocity profile does not reaches the same max value through time. If only convection occured, there should be only a transport of the bulk property through space. Although the physics of the problem only take into account convection, we have (unintentionally) included what is called as nummerical diffusion. This is, when we approximate the spatial derivative, we are including an error term (in this case first order derivative means error as Deltax^2), that is similar to the diffussion term (second sptial derivative of velocity). This numerical diffusion is the responsible for such a loss in the velocity profile.

Another thing to mention is about the Courant Number. The stability of the mathematical approximation in the code is, as explained above, a function of the Courant Number. In this case, C=0.2, i.e. smaller than 1 -> Stable.

Figure 2) 40 Grid Points

First of all, the initial profile is more similar now to a square wave than before, due to the increase of points in the grid.

In this second case, the loss of shape is reduced, and the maximum value that the final profile reaches is higher. This is because, as it was said before, the numerical diffusion term is proportional to the square of the x-increment in the integration. In this case, using more grid points means that the difussion induced is smaller.

The Courant number here is C=0.4, and therefore the numerical scheme is stable.

Figure 3) 80 Grid Points

For the same reason than before, the profile continues reaching values closer to the initial maximum, and the shape loss decreases. Again, the nummerical diffusion induced keeps decreasing.

The Courant number is now C=0.8, and therefore the numerical scheme is stable.

Figure 4) 160 Grid Points

We clearly see that the trend stopped. Why? If we were increasing the number of grid points, i.e. reducing the x-step in the nummerical integration, why do we suddenly have these results?

Well, the nummerical diffusion might have decreased, but now the mathematical scheme is not stable. This is because of the CFS criteria stated at the beggining. The Courant number now has reached a value C=1.6>Cmax. Therefore the step is too small and the program cant follow the solution. This causes extreme disturbance and oscilations.

Figure 5) Optimum Solution : 100 Grid Points

The optimum solution, this is, the one qhere the minimum numerical diffusion is induced and the scheme is still stable, is when C=Cmax=1. Basics calculations can lead to the optimum size step dx=dt=0.01, this is, 100 Grid points.

In this case, the profile does not lose the squared shape, and the maximum value keeps being 2m/s, as at the beggining. The nummerical diffusion is therefore so small that we cannot appreciate it. Also, notice that the initial profile is almost a square wave now, due to the high number of points in the grid.