MATLAB Program to solve the 1D linear wave equation

WEEK 4:

(Effect of Grid-Size on output for the solution of 1D linear wave equation)

1. Problem Setup:

Given Partial Differential Equation:

`(delu)/(delt) + c(delu)/(delx) = 0`

Numerical Discretization: 

`u_(i,n+1) = u_(i,n) + (cDeltat)/(Deltax)*(u_(i-1,n+1) - u_(i,n+1))`

where `u_(i,n)` denotes the space index and the time index respectively

Here the domain length is assumed to be L = 1m

Boundary condtions: 1. At x = 0; velocity u = 1m/s
                               2. At x = 1; velocity u = 1m/s
(Drichlet boundary conditons are specified for this problem)

Initial conditions: At t=0; u=1m/s at x = 0

The Initial velocity profile of the problem is shown below for grid size (n=160) for reference:

Discretization Specifications:

1. Time derivative was discretized by first order forward differencing:

∂u/∂t = (u(i,n+1)-u(i,n))/Δt

2. Space derivative was discretized by first order backward differencing:

∂u/∂x = (u(i,n)-u(i-1,n))/Δx

Other specifications:

1. C = 1
2. Time-step of calculation = 0.01s
3. Required time for output = 0.4s

2. MATLAB CODE:

%% MATLAB PROGRAM TO SOLVE THE 1D LINEAR WAVE EQUATION:

% Initial conditions:
L = 1; %Length 
c = 1; 
t = 0.4;
n = input('Enter the grid size: ');
dt = 0.01; %entering the timestep
nt = t/dt; % number of timesteps
x_start = 0.1;
x_stop = 0.3;
error = 1000; %Initializing the value of error
tol = 0.001; %tolerance value
% Initializing variables:
x = linspace(0,L,n+1); %initalizing the x array
dx = x(2) - x(1); %grid spacing
[n_start,n_stop] = func1(x_start, x_stop, n, L);
r = (c*dt/dx);
%Definiing the velocity step function
u = ones(1,n);
u(1,n_start:n_stop) = 2;
u_primitive = u; %Storing the initial velocity profiles
uold = u; % Defining the old value of u initially
% Boundary conditions:
u_initial = 1;
CFL_number = u(1,n_start)*dt/dx;

% Initializing a quantity i.e. a break point to stop the execution
if n == 160
    break_counter = 10000;
end
if n == 80
    break_counter = 5000;
end
if n == 40
    break_counter = 2000;
end
if n == 20
    break_counter = 1000;
end
    
% Solving the velocities by Gauss Siedel method
for i = 1:nt
    num_iter = 1;
    while (error > tol)
        for j=2:n-1
            if num_iter > break_counter
                break
            else
            u(j) = u(j)*(1-r) + r*u(j-1);
            end
        end
        % Defining the error term
        error = max(abs(uold-u));
        % Updating the old values
        uold = u;
        % Incrementing the loop counter
        num_iter = num_iter + 1;
        
        % Plotting the solution dynamically
        subplot(2,1,1);
        figure(1);
        plot(x(2:end),u)
        xlabel('length');
        ylabel('velocity');
        title_text = sprintf('Velocity profile for 1D wave equation \n Number of Iterations: %d \n (Grid-size: %d) (CFL Number: %.3f) (Time-step: %.4f)',num_iter,n,CFL_number,dt);
        title(title_text);
        pause(0.001);
        hold off;
    end
end

% Comparison of the initial and final velocity profiles:
subplot(2,1,2);
hold on;
plot(x(2:end),u_primitive); 
plot(x(2:end),u);
xlabel('length');
ylabel('velocity');
title('Comparison of velocities'); 
legend('Initial velocity profile','Final velocity profile');

3. RESULTS AND DISCUSSIONS:

RUN 1: SETTING UP THE TIME-STEP TO dt = 0.01 s

1. For grid-size (n = 20)
By setting up the grid size to n = 20 we can see that the CFL number being 0.4 gives us a distorted and ill-formed waveform at the final time-step.

2. For grid-size (n=40)

When the grid-size is 40 we can see that the solution completely loses its shape and is close to being unstable because of CFL number being close to 1.

3. For grid-size (n=80)

When the grid size is 80, the CFL number obtained in this case is beyond 1 hence it gives rise to an unstable solution.

4. For grid-size (n=160)

When the grid-size is 160, the CFL Number is 3.2 which is way beyond 1. The solution is completely unstable and does not converge at all.

The next set of observations were made by decreasing the time-step to a very low value so that the CFL Number << 1 and good convergence can be achieved.

RUN 2: SETTING THE TIME-STEP TO dt = 0.0001s

1. For grid-size (n = 20)

2. For grid-size (n=40)

3. For grid-size (n=80)

4. grid-size (n=160)

In the above plots we can see that for low CFL Numbers we get good convergence.
As the grid size increases the graph moves more and more forward in time as the iterations are performed for more number of times.

4. SUMMARY OF RESULTS:

RUN 1 (dt = 0.01)
Grid-size Number of iterations CFL Number Iteration breakpoint Break-point reached

20

52 0.400 1e3 No
40 90 0.800 2e3 No
80 36 1.600 3e3 No
160 102 3.200 4e3 No
RUN 2 (dt = 0.0001)
Grid-size Number of iterations CFL Number iteration breakpoint Break-point reached
20 349 0.004 1e3 No
40 656 0.008 2e3 No
80 1281 0.016 3e3 No
160 2554 0.032 4e3 No

 


Projects by Priyotosh Bairagya

LAMINAR INCOMPRESSIBLE FLOW SIMULATION THROUGH A PIPE IN OPENFOAM: [PART 2/2]: SYMMETRY BOUNDARY CONDITION: OBJECTIVES:1. Creating the Mesh-Script (blockMeshDict file) with symmetry boundary condition.2. Simulation Results and post processing the velocity profile at Read more

LAMINAR INCOMPRESSIBLE FLOW SIMULATION THROUGH A PIPE IN OPENFOAM: [PART 1/2]: WEDGE BOUNDARY CONDITION: OBJECTIVES:1. Calculation of Pre-eliminary quantinites related to the flow.2. Creating the Mesh script (blockMeshDict file) for specifying the geometry and bounda Read more

I. Interpolation Schemes in Finite Volume Method: The approximation of surface and volume integrals may require values of the variable at locations other than the computational nodes of the CV. Values at these locations are obtained using interpolation formulae. Some o Read more

ANALYSIS OF NUMERICAL STABILITY OF VARIOUS ITERATIVE SOLVERS FOR TRANSIENT 2D HEAT CONDUCTION: [Part: 3/3] INTRODUCTION: The criterion of stability of a numerical scheme is determined by the way the errors propagate while the solution moves from one time-step to th Read more

MESH GENERATION AND ANALYSIS USING BLOCKMESH FOR FLOW OVER A BACKWARD FACING STEP: The purpose of the following project is to generate the geometry for a variation of the incompressible cavity flow problem in OpenFOAM. For this purpose we have modified the lid-driven c Read more

ANALYSIS OF VARIOUS ITERATIVE SCHEMES FOR THE SOLUTION OF A 2D HEAT CONDUCTION PROBLEM: [PART: 2/3] In the previous part we had explained the problem statement and the MATLAB Program in detail. In this Part we are going to explain the outputs from the 2D Heat Conduct Read more

NUMERICAL SOLUTION OF 1D SUPERSONIC NOZZLE FLOW SIMULATION BY MAC-CORMACK METHOD PROJECT OBJECTIVES: i. Numerical solution of the governing equations in both conservative and non-conservative forms. ii. Creating user defined functions for calculating the flow quan Read more

Simulation of a 2D Heat Conduction problem in steady and unsteady/transient forms using iterative methods. Project Objectives: 1. Solving the 2 Dimensional Heat conduction equation in the generalized form using various iterative techniques: i. Explicit Solver (for Read more

UNDERSTANDING LINEAR SYSTEMS(ANALYSIS OF VARIOUS ITERATIVE SCHEMES TO SOLVE A SYSTEM OF LINEAR EQUATIONS TO FIND THE EIGEN VALUES AND SPECTRAL RADIUS) (A) PROBLEM STATEMENT: Given coefficient matrix: `A = [[5,1,2],[-3,9,4],[1,2,-7]]` Given Solution Matrix: `X = [[x Read more


Loading...

The End