## MATLAB Program to solve the 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

### Pipe flow simulation in OpenFOAM (Part 2/2: Symmetry boundary conditions) Priyotosh Bairagya · 2018-10-08 11:27:26

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

### Pipe flow simulation in OpenFOAM (Part 1/2: Wedge boundary conditions) Priyotosh Bairagya · 2018-10-07 05:05:54

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

### Flux Limiters and Interpolation Schemes in Finite Volume Method Priyotosh Bairagya · 2018-10-05 12:18:40

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 Solution Stability in a 2D Heat conduction problem Priyotosh Bairagya · 2018-09-20 22:56:03

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

### BlockMesh Analysis of a Backward Facing Step Priyotosh Bairagya · 2018-09-20 17:08:11

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 Steady and Unsteady State solutions of a 2D Heat conduction problem Priyotosh Bairagya · 2018-09-16 10:37:07

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

### Flow simulation of a 1D Super-Sonic nozzle using the Mac-Cormack method Priyotosh Bairagya · 2018-09-02 12:55:50

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 2D heat conduction in steady and unsteady forms Priyotosh Bairagya · 2018-08-31 14:34:23

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

### Iterative solution of a system of linear equations and an analysis of spectral radius of a matrix Priyotosh Bairagya · 2018-08-31 10:13:18

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...