Site map
 Partenov CFD

# Numerical solution of mathematical problems and equations

1. A 1D wave equation with impact is numerically solved using MATLAB script. It uses the Euler method in time and second-order accuracy central finite difference scheme in space. It can apply to different places an impact with desired forces.
```% 1D wave equation after impact
clear
clc

l = 1; % length of road, m
c = 330; % speed of sound in material, m/s

dx = 0.002; % space discretization, m
dt = 0.9*dx/c; % time discretization, s

r = c*dt/dx; % Courant Number
Im = 0.0015; % impact constant

a = 0.5*l; % distance where impact apply
b = l - a;

T = 1200; % amount of time steps
N = floor(l/dx) + 1; % number of space steps

u = zeros(N,1); % fill displacement vector
uminus = zeros(N,1);
uplus = zeros(N,1);

x = zeros(N,1);
tt = zeros(T,1);
uu = zeros(T,1);

% initial condition

x(1) = 0;
for i = 1: N-1
x(i+1) = x(i) + dx;
if x(i) <= a
u(i) = sqrt(Im*b*x(i)*(l^2 - x(i)^2 - b^2)/l)^1.5;
else
u(i) = sqrt(Im*b*(l*(x(i) - a)^3/b + (l^2 - b^2)*x(i) - x(i)^3)/l)^1.5;
end
end

for i = 1: N
u(i) = 0.001*normpdf(x(i),a,0.1);
end

u(1) = 0;
u(N) = 0;

% define uminus
for i = 2:N-1
uminus(i) = u(i) + 0.5*r^2*(u(i+1) - 2*u(i) + u(i-1));
end

figure

t = 0;
for i = 1:T
tt(i) = t;
t = t + dt; % time
for j = 2:N-1
uplus(j) = 2*u(j) - uminus(j) + (r^2)*(u(j+1) - 2*u(j) + u(j-1));
end
uplus(1) = 0; % boundary condition
uplus(N) = 0; % boundary condition
uminus = u;
u = uplus;
uu(i) = u(floor(N/2));

% animation
plot(x, u );
axis([0 l -0.005 0.005]);
xlabel('length of string, m');
ylabel('amplitude of wave, m');
drawnow;

end

figure
plot(tt, uu);
xlabel('time, s');
ylabel('amplitude of wave, m');
```

2. The numerical solution of a 2D wave equation with impact is shown using MATLAB script. The equation is solved in a circular domain using the polar coordinates. It uses Euler's method in time and second-order accuracy central finite difference scheme in space. It can choose the position and strength of initial impact conditions.
```% 2D wave equation after impact in circular
clear
clc

R = 1; % radius of circle, m
c = 6000; % wave speed in material, m/s

% polar grid
Nr = 40; % number of space steps in radial direction
hr = R/Nr; % space discretization in radial direction, m

Nt = 40; % number of space steps in tangential direction, Nt must equal to Nr, for use function polarplot3d()
ht = 2*pi/Nt; % space discretization in tangential direction, rad

dt = 0.5*hr*ht/(c*sqrt(ht^2 + 1)); % time discretization, s
T = 520; % amount of time steps

tt = zeros(T,1);
uu = zeros(T,1);

rho = zeros(Nr+1,1);
theta = zeros(Nt+1,1);
rhot = zeros(Nt+1,1);

rho0 = 0.0; % position where is apply the initial condition
theta0 = pi/6; % position where is apply the initial condition

rho(1) = 0;
for i = 2 : Nr + 1
rho(i) = rho(i-1) + hr;
end

theta(1) = 0;
for i = 2 : Nt + 1
theta(i) = theta(i-1) + ht;
polar([theta(i-1) theta(i-1)], [0 1]);
hold on
end

% plot polar grid

for j = 1: Nr + 1
for i = 1 : Nt + 1
rhot(i) = rho(j);
end

polar(theta, rhot)
title('Polar grid');
hold on
end

u = zeros(Nr + 1, Nt + 1); % define sorage of wave function
uminus = zeros(Nr + 1, Nt + 1 );
uplus = zeros(Nr + 1, Nt + 1 );
col = zeros(Nr + 1, Nt + 1 );

% set initial condition

for i = 1 : Nr + 1
for j = 1 : Nt
l = sqrt(rho0^2 + rho(i)^2 -2*rho0*rho(i)*cos(theta(j) - theta0));
u(i,j) = -0.005*normpdf(l,0,0.8);
end
end

rl = 0;
for i = 2 : Nr
rl = rl + hr;
for j = 1 : Nt
if j == Nt
term3 = (u(i,1) - 2*u(i,j) + u(i,j-1))/(rl^2*ht^2);
else
if j == 1
term3 = (u(i,j+1) - 2*u(i,j) + u(i,Nt))/(rl^2*ht^2);
else
term3 = (u(i,j+1) - 2*u(i,j) + u(i,j-1))/(rl^2*ht^2);
end
end
term1 = (u(i+1,j) - 2*u(i,j) + u(i-1,j))/hr^2;
term2 = (u(i+1,j) - u(i-1,j))/(rl*hr);
uminus(i,j) = u(i,j) + 0.5*c^2*dt^2*(term1 + term2 + term3);
end
end

sum = 0;
for i = 1:Nt
sum = sum + (u(2,i) - u(1,1));
end

uminus(1,1) = u(1,1) + 2*c^2*dt^2*sum/((Nt+1)*hr^2);
uminus(1,:) = uminus(1,1);

uminus(:, Nt+1) = uminus(:, 1);

figure('color','white');
col(:,:) = 12;
[x, y] = pol2cart(theta,rho);
t = 0;
for k = 1:T
tt(k) = t;
t = t + dt;
sum = 0;
for i = 1:Nt
sum = sum + + (u(2,i) - u(1,1));
end
uplus(1,1) = 2*u(1,1) - uminus(1,1) + 4*c^2*dt^2*sum/((Nt+1)*hr^2);
uplus(1,:) = uplus(1,1);

for i = 2: Nr
rl = rl + hr;
for j = 1 : Nt
if j == Nt
term3 = (u(i,1) - 2*u(i,j) + u(i,j-1))/(rl^2*ht^2);
else
if j == 1
term3 = (u(i,j+1) - 2*u(i,j) + u(i,Nt))/(rl^2*ht^2);
else
term3 = (u(i,j+1) - 2*u(i,j) + u(i,j-1))/(rl^2*ht^2);
end
end
term1 = (u(i+1,j) - 2*u(i,j) + u(i-1,j))/hr^2;
term2 = (u(i+1,j) - u(i-1,j))/(rl*hr);
uplus(i,j) = 2*u(i,j) - uminus(i,j) + c^2*dt^2*(term1 + term2 + term3);
end
end

uplus(Nr+1, :) = 0; % boundary conditions
uplus(:, Nt+1) = uplus(:, 1);

uminus = u;
u = uplus;

uu(k) = u(1,1);

% animation
% if you want to speed up the proccess put next lines under the end of
% loop, so will eliminate animation
%polarplot3d(u );
%axis([-1.4 1.4 -1.4 1.4 -0.005 0.005])
%xlabel('x, m');
%ylabel('y, m');
%zlabel('amplitude of wave, m');
%drawnow;

end

figure
plot(tt, uu);
xlabel('time, s');
ylabel('amplitude of wave, m');
```

3. I present two methods for a numerical solution to Laplace’s equation in the 2D square domain: Jacobi’s iterative and the Gauss-Seidel SOR method. The suggested methods are compared with the analytical solution for accuracy. To discretize the problem is applied finite differences scheme. The number of iterations significantly drops when using the Gauss-Seidel SOR method for the same accuracy.
```% Numerical solution to Laplace's Equation
clear;
clc;

N = 60; % number of space divisions in x and y
L = 1; % length of square

deltax = L/N; % space step in x direction
deltay = deltax; % space step in y direction

% exact solution
n = 210; % number of sumation
u = zeros(N+1, N+1); % initialize the exact solution matrix
u(:,1) = 1; % boundary condition at left side

y = 0;
for i = 2 : N
y = y + deltay;
x = 0;
for j = 2 : N
u(i,j) = 0;
x = x + deltax;
for k = 1 : n
An = (1/sinh(-k*pi))*((-2/(k*pi))*cos(k*pi/L) + 2/(k*pi));
u(i,j) = u(i,j) + An*(sin(k*pi*y/L))*sinh(k*pi*(x - L)/L);
end
end
end

% numerical solution

tol = 1e-04; % tolerance of solution

% Jacobi's iterative method

uj = zeros(N+1, N+1); % initialize the numerical solution matrix, Jacobi method
uj(:,1) = 1; % boundary condition at left side

ujold = uj;
e = 1;
kj = 0;
while e > deltay*tol
kj = kj + 1; % number of iteration
for i = 2 : N
for j = 2 : N
uj(i,j) = 0.25*(ujold(i+1,j) + ujold(i-1,j) + ujold(i,j+1) + ujold(i,j-1)); % solution using Jcobi method
end
end
e = max(max(uj - ujold));

ujold = uj;
end

% Point Gauss-Seidel SOR solver

omega = 1.85; % relaxation parameter

ug = zeros(N+1, N+1); % initialize the numerical solution matrix, Gauss-Seidel SOR solver
ug(:,1) = 1; % boundary condition at left side

l = 0;
for omega = 1.95  % relaxation parameter
l = l + 1;
kg = 0;
ugold = zeros(N+1, N+1); % initial solution
ug(:,1) = 1; % boundary condition
e = 1;
while e > deltay*tol
kg = kg + 1; % number of iteration
for i = 2 : N
for j = 2 : N
ug(i,j) =(1 - omega)*ug(i, j) + omega*0.25*(ugold(i+1,j) + ug(i-1,j) + ugold(i,j+1) + ug(i,j-1)); % solution using Point Gauss-Seidel SOR method
end
end
e = max(max(ug - ugold));

ugold = ug;
end
iteration(l) = kg;
om(l) = omega;
end

X = 0:deltax:1;
Y = 0:deltay:1;

e = abs(ug - u); % error
plot(X(2:N), e(N/2,2:N))
grid on
xlabel('x');
ylabel('Absolute Error')
title('Error in the numerical solution y =0.5');

figure
mesh(X, Y, ug);
xlabel('x');
ylabel('y');
zlabel('value of u');
title('Solution for N = 60');
grid on
```

4. The 2D Poisson equation is solved numerically. The proposed equation describes the stress function in the torsion of a prismatic rod. The solution drove the linear algorithm for Poisson's equation. This method uses the finite difference scheme, in which the average values of the neighboring four points give the solution at a given central point. The suggested algorithm is accelerated using a suitable weighting factor. This factor increases up to two as the number of nodes in the domain becomes infinite. The modified algorithm is much faster than the simple linear algorithm for the same accuracy.
```% Solution of 2D Poisson's equation in square domain
clc
clear
N = 45; % mesh size
H = -0.05; % parameter

toll = 1e-06; % tollerance

h = 1/(N-1); % grid spacing

maxit =  5000; % maximum number of iteration

error_vector = zeros(N^2,1);

b = zeros(N^2,1);
b(1:N^2) =  -H*h^2;
p  = zeros(N^2,1); % vector where is stored the potential function, this is initial guess also
p_old = zeros(N^2,1);

p(1:N)           = 0; % bottom boundary conditionn
p((N-1)*N:N^2)   = 0; % top boundary condition
p(1:N:(N-1)*N+1) = 0; % left boundary condition
p(N:N:N^2)       = 0; % right boundary condition

% main loop of solution, according the equetion (19)
for i = 1: maxit
p_old = p;
for k = 2:N - 1
for j = 2:N - 1
n = j + (k - 1)*N;
p(n) = (p_old(n+1) + p_old(n-1) + p_old(n+N) + p_old(n-N) + b(n))/4;
end
end

for  k = 2:N-1
for j = 2: N-1
n = j + (k - 1)*N;
error_vector(n) = (p(n) - p_old(n))/p(n);
end
end

error = max(abs(error_vector));
if error <= toll
break;
end
end

number_of_iteration_1 = i

x = zeros(N,1);
y = zeros(N,1);
for i = 2: N
x(i) = x(i-1) + h;
y(i) = y(i-1) + h;
end

% prepare to plot

pp = zeros(N,N);
for n  = 1: N*N
k = fix((n-0.1)/N) +1;
j = n-(k-1)*N;
pp(k,j) = p(n);
end
mesh(x,y,pp);

if mod(N,2) == 1 % plot a centerline section
nb=fix(N/2)*N+1;
pc=p(nb:nb+N-1);
figure;
plot(x,pc)
end

% main loop of solution, according the equation (25)
for i = 1: maxit
p_old = p;
for k = 2:N - 1
for j = 2:N - 1
n = j + (k - 1)*N;
p(n) = (p(n+1) + p(n-1) + p(n+N) + p(n-N) + b(n))/4;
end
end

for  k = 2:N-1
for j = 2: N-1
n = j + (k - 1)*N;
error_vector(n) = (p(n) - p_old(n))/p(n);
end
end

error = max(abs(error_vector));
if error <= toll
break;
end
end

number_of_iteration_2 = i

% solution with acceleration scheme

ksi = 1 - 2*(sin(pi/(2*N)))^2;
omega = 2/(1 + sqrt(1 - ksi*ksi)); % weighting factor

v  = zeros(N^2,1);

% main loop of solution, according the equation (28)
for i = 1: maxit
p_old = p;
for k = 2:N - 1
for j = 2:N - 1
n = j + (k - 1)*N;
v(n) = (p(n+1) + p(n-1) + p(n+N) + p(n-N) + b(n))/4;
end
end

p(1:N*N) = omega*v(1:N*N) + (1 - omega)*p_old(1:N*N);

for  k = 2:N-1
for j = 2: N-1
n = j + (k - 1)*N;
error_vector(n) = (p(n) - p_old(n))/p(n);
end
end

error = max(abs(error_vector));
if error <= toll
break;
end
end

number_of_iteration_3 = i

% prepare to plot
figure;
pp = zeros(N,N);
for n  = 1: N*N
k = fix((n-0.1)/N) +1;
j = n-(k-1)*N;
pp(k,j) = p(n);
end
mesh(x,y,pp);

if mod(N,2) == 1 % plot a centerline section
nb=fix(N/2)*N+1;
pc=p(nb:nb+N-1);
figure;
plot(x,pc)
end

% from number_of iteration, it can be view that acceleration is very fast
% compare with another methods
```
```number_of_iteration_1 =
3119
number_of_iteration_2 =
136
number_of_iteration_3 =
2
```

5. The 2D heat equation with heat generation is solved. The domain is quadratic, and the location and geometry of the source can be different. The second-order finite difference method is applied. The right and top boundary conditions are constant temperatures: Dirichlet boundary condition. The left and bottom boundary conditions are Neumann boundary conditions. Validation of code is done by checking the heat flow balance through the walls.

## Contents

```clear;
clc
```

## geometry

outer geometry

```Lx = 1; % length in x direction, [m]
Ly = 1; % length in y direction, [m]

% heat source geometry
Lhx = 0.4;
Lhy = 0.3;

%position of heat source
Lpx = 0.5;
Lpy = 0.2;
```

## initial parameters

```h    = 38; % heat transfer coefficient at left boundary [W/m^2.K]
K    = 8; % thermal conductivity of domain [W/m.K]
Tinf = 8;  % temperature of media [degree]
T1   = 23; % temperature at top boundary of domain [degree]
T2   = 29; % temperature at right boundary of domain [degree]
q    = 36000; % heat source [W/m^3]
```

## grid

```n = 8; % n can accept only positive natural number

dx = 0.1/n; %  mesh size in x direction, [m]
dy = dx; % mesh size in y direction, [m]

Nx = Lx/dx + 1; % number of points in x direction
Ny = Ly/dy + 1; % number of points in y direction

N = Nx*Ny; % total number of points
```

## plot geometry and grid

```figure;

plot([0 Lx Lx 0 0],  [0 0 Ly Ly 0], 'b' ); % plot geometry
hold;
for i = 1:Nx-1
plot([i*dx i*dx], [0 Ly], 'b'); % plot mesh in x direction
end
for i = 1:Ny-1
plot([0 Lx], [i*dy i*dy], 'b'); % plot mesh in y direction
end
plot([Lpx Lpx+Lhx Lpx+Lhx Lpx Lpx], [Lpy Lpy Lpy+Lhy Lpy+Lhy Lpy], 'r'); % plot of heat source geometry
axis equal;
axis([-0.1 1.1 -0.1 1.1]);
grid on;
xlabel('X [m]');
ylabel('Y [m]');
title('Geometry and grid')
```
```Current plot held
```

## identify nodes

```node_type(1:N) = '6'; % inner node
k = 1;
l = 1;
m = Nx;
for i = 1 : N
if i == 1
node_type(i) = '1'; % left bottom corner
else
if i <= Nx-1
node_type(i) = '2'; % bottom boyndary
else
if i == k + Nx && i <= N - 2*Nx + 2
k = k + Nx;
node_type(i) = '5'; % left boundary
else
if i >= l + (Lpy/dy)*Nx + Lpx/dx && i <= l + (Lpy/dy)*Nx + Lpx/dx + Lhx/dx && i < Nx*(Lpy + Lhy)/dy + Lpx/dx + Lhx/dx + 2
node_type(i) = '7'; % source area
if i > Nx*Lpy/dy && i < Nx*Lpy/dy + Nx
if i > l + (Lpy/dy)*Nx + Lpx/dx && i < l + (Lpy/dy)*Nx + Lpx/dx + Lhx/dx
node_type(i) = '9';
end
end
if i > Nx*(Lpy + Lhy)/dy
if i > l + (Lpy/dy)*Nx + Lpx/dx && i < l + (Lpy/dy)*Nx + Lpx/dx + Lhx/dx
node_type(i) = '9';
end
end
if i == l + (Lpy/dy)*Nx + Lpx/dx
if i == 1 + (Lpy/dy)*Nx + Lpx/dx && i < Nx*(Lpy + Lhy)/dy + Lpx/dx + Lhx/dx + 2
node_type(i) = '8';
else
if i <  Nx*(Lpy + Lhy)/dy
node_type(i) = '9';
else
node_type(i) = '8';
end
end
end
if i == l + (Lpy/dy)*Nx + Lpx/dx + Lhx/dx
l = l + Nx;
if l == 1 + Nx
node_type(i) = '8';
else
if i <  Nx*(Lpy + Lhy)/dy
node_type(i) = '9';
else
node_type(i) = '8';
end
end
end
else
if i > N - Nx
node_type(i) = '4'; % top boundary
else
if i == m
m = m + Nx;
node_type(i) = '3'; % right boundary
end
end
end
end
end
end
end
```

## Initialize and populate coefficient matrix A and vector b

```A = zeros(N,N); % initialize matrix A
b = zeros(N,1); % initialize vector b
% populate vector b
for i = 1:N
if node_type(i) == '1'
b(i) = -h*dy*Tinf/K;
else
if node_type(i) == '2'
b(i) = 0;
else
if node_type(i) == '3'
b(i) = T2;
else
if node_type(i) == '4'
b(i) = T1;
else
if node_type(i) == '5'
b(i) = -2*h*dy*Tinf/K;
else
if node_type(i) == '6'
b(i) = 0;
else
if node_type(i) == '7'
b(i) = -q*dx*dy/K;
else
if node_type(i) == '8'
b(i) = -q*dx*dy/(4*K);
else
if node_type(i) == '9'
b(i) = -q*dx*dy/(2*K);
end
end
end
end
end
end
end
end
end
end

a = dy/dx + dx/dy;
c = h*dy/K;
% populate coefficient matrix A
for j = 1:N
if  j == 1
A(j,j) = -(a + c);
A(j,j+1) = dy/dx;
A(j,j+Nx) = dx/dy;
end
if j > 1 && node_type(j) == '2'
A(j,j)    = -2*a;
A(j,j-1)  = dy/dx;
A(j,j+1)  = dy/dx;
A(j,j+Nx) = 2*dx/dy;
else
if node_type(j) == '3' || node_type(j) == '4'
A(j,j) = 1;
else
if node_type(j) == '5'
A(j,j) = -2*(c + a);
A(j,j+1) = 2*dy/dx;
A(j,j-Nx) = dx/dy;
A(j,j+Nx) = dx/dy;
else
if node_type(j) == '6' || node_type(j) == '7' ||  node_type(j) == '8' ||  node_type(j) == '9'
A(j,j) = -2*a;
A(j,j+1) = dy/dx;
A(j,j-1) = dy/dx;
A(j,j+Nx) = dx/dy;
A(j,j-Nx) = dx/dy;
end
end
end
end
end
```

## Solution of system of equation to find T

```T = zeros(Ny,Nx); % initialize solution vector

X = A\b; % solution

% reshape temperature vector into a matrix
o = 0;
for j = 1:Ny
for i = 1:Nx
o = o + 1;
T(j,i) = X(o);
end
end
```

## validation of code

```% heat flow over left boundary
q_left = h*dy*(Tinf - T(1,1))/2;
for j = 2:Ny - 1
q_left = q_left + h*dy*(Tinf - T(j,1));
end

% heat flow over right boundary
q_right = K*dy*(T2 - T(1,Nx-1))/(2*dx);
for j = 2:Ny - 1
q_right = q_right +  K*dy*(T2 - T(j,Nx-1))/dx;
end
q_right = q_right +  K*dy*(T2 - T(Ny,Nx-1))/(2*dx);

% heat flow over top boundary
q_top = K*dx*(T1 - T(Ny-1,1))/(2*dy);
for i = 2: Nx - 1
q_top = q_top + K*dx*(T1 - T(Ny-1,i))/dy;
end
q_top = q_top + K*dx*(T1 - T(Ny-1,Nx))/(2*dy);

error =q*Lhx*Lhy + (q_top+q_left+q_right); % absolute error heat transfer balance
```

## Plot results

```x = 0:dx:Lx;
y = 0:dy:Ly;
[xx, yy] = meshgrid(x, y);

% 3D plot of temperature distribution
figure
surf(xx, yy, T, 'EdgeColor', 'none');
xlabel('X [m]');
ylabel('Y [m]');
zlabel('Temperature [ ^oC]');
colormap(jet)
title('Temperature distribution');

% 2D plot of temperature project over the domain
figure
contour(xx, yy, T, 15);
axis equal;
xlabel('X [m]');
ylabel('Y [m]');
title('Contour of temperature ');
colormap(jet);
colorbar

% plot temperature on edges
figure
subplot(2,1,1);
plot(x, T(1,:)); % plot temperature on bottom boundary
xlabel( 'X [m]');
ylabel( 'Temperature [ ^oC]');
title('Temperature on bottom boundary');
grid on

subplot(2,1,2);
plot(y, T(:,1)); % plot temperature on left boundary
xlabel( 'Y [m]');
ylabel( 'Temperature [ ^oC]');
title('Temperature on left boundary')
grid on
```

6. The quasi-one-dimensional Euler system for nozzle and shock tubes is solved. The Yee-Roe Flux Limited Scheme was successfully implemented and numerically tested. The numerical scheme uses the explicit Euler method to iterate the solution in time and second-order upwind difference in space. The solution examined the influence of the number of nodes, Currant number, and state parameters on the solution. The code changed to simulate other situations as the reflection of shock waves from the side of the tube using the concrete boundary conditions for the nozzle and supersonic wind tunnel. The results showed that using the flux limiter is a significant part of the solver to get the proper solutions. The code performs to test the first-order solution too. The results from this test were with very high oscillation in front of shock waves. The Currant number must be sufficiently small to achieve a stable solution for this case. The solver showed good flexibility to use for a great variety of problems.

## Contents

```% a quasi-one dimensional flow solver applied to nozzles and shock-tubes

clear;
clc;

N = 100; % number of space points
Tf = 0.007; % final time of simulation
```
```% phisical constant
k = 1.40; % heat capacity ratio, of air as ideal gass
```
```% initial data
L = 10; % length of nuzzle or tube, for wind tunnel L = 20m
h = L/(N - 1); % space step
X = 0:h:L;

% geometry of nozzle
S = zeros(1, N);
dS = zeros(1,N);
for i = 1:N
[S(i), dS(i)] = nozzle_geo(X(i), L ,'shock'); % if last variable is 'nozzle', then nozzle geometry is inputed, 'shock' for  shock tube, and  'wind' for wind tunnel
end

ncase = 3; % switch the case 1, case 2 and case 3

switch ncase
case 1
rhoL = 1.1409; % density at left side, kg/m^3
uL = 65.451; % velocity at left side, m/s
pL = 9.7534e+4; % pressure at left side, Pa

rhoR = 1.1008; % density at right side, kg/m^3
uR = 113.060; % velocity at right side, m/s
pR = 9.2772e+4; % pressure at right side, Pa

case 2
rhoL = 1.1308; % density at left side, kg/m^3
uL = 80.054; % velocity at left side, m/s
pL = 9.6328e+4; % pressure at left side, Pa

rhoR = 1.0275; % density at right side, kg/m^3
uR = 150.535; % velocity at right side, m/s
pR = 8.4974e+4; % pressure at right side, Pa

case 3
rhoL = 1.0; % density at left side, kg/m^3
uL = 0; % velocity at left side, m/s
pL = 1.0e+5; % pressure at left side, Pa

rhoR = 0.125; % density at right side, kg/m^3
uR = 0; % velocity at right side, m/s
pR = 1.0e+4; % pressure at right side, Pa

case 4 % wind tunnel parameters
rhoL = 1.0985; % density at left side, kg/m^3
uL = 82.054; % velocity at left side, m/s
pL = 9.406e+4; % pressure at left side, Pa

rhoR = 0.3056; % density at right side, kg/m^3
uR = 0;% velocity at right side, m/s
pR = 1.5741e+4; % pressure at right side, Pa
end

EL = pL/(rhoL*(k - 1)) + 0.5*uL^2; % internal energy of the gas at left side
HL = EL + pL/rhoL; % enthalpy of the at left side
aL = sqrt((k - 1)*(HL - 0.5*uL^2)); % speed of sound at left side, m/s

ER = pR/(rhoR*(k - 1)) +  0.5*uR^2; % internal energy of the gas at right side
HR = ER + pR/rhoR; % enthalpy of the at right side
aR = sqrt((k - 1)*(HR - 0.5*uR^2)); % speed of sound at right side, m/s

CFL = 0.5; % Currant number

l = CFL*h/max(aL, aR); % initial time step
```

## solution

```Q = zeros(3, N); % initialize matrix of solution
Qold = zeros(3, N);
F = zeros(3, N);
deltaF = zeros(3, N);

lambda1 = zeros(3,3);
lambda2 = zeros(3,3);

M1 = zeros(3,1);
M2 = zeros(3,1);
M3 = zeros(3,1);
M4 = zeros(3,1);

N1 = zeros(3,1);
N2 = zeros(3,1);

rho = zeros(1,N);
u = zeros(1,N);
p = zeros(1,N);
H = zeros(1,N);
E = zeros(1,N);
a = zeros(1,N);

rhoav = zeros(1,N-1);
uav   = zeros(1,N-1);
Hav   = zeros(1,N-1);
aav   = zeros(1,N-1);
pav   = zeros(1,N-1);

maxl = zeros(1,N-1);

Res = zeros(1,8000); % residual vector
Time = zeros(1,8000); % time vector
Resx = zeros(1,N);

for i = 1 :N
if X(i) <= L/2
Q(:,i) = [rhoL*S(i);rhoL*uL*S(i);rhoL*EL*S(i)]; % initial value of solution
rho(i) = rhoL; % density
u(i) = uL;  % velocity
p(i) = pL;  % pressure
H(i) = HL;
E(i) = EL;
a(i) = aL;
F(:,i) = [rhoL*uL*S(i); (rhoL*uL^2 + (k - 1)*rhoL*(E(i) - 0.5*uL^2))*S(i); rhoL*uL*H(i)*S(i)];
else
Q(:,i) = [rhoR*S(i);rhoR*uR*S(i);rhoR*ER*S(i)];
rho(i) = rhoR;
u(i) = uR;
p(i) = pR;
H(i) = HR;
E(i) = ER;
a(i) = aR;
F(:,i) = [rhoR*uR*S(i); (rhoR*uR^2 + (k - 1)*rhoR*(E(i) - 0.5*uR^2))*S(i); rhoR*uR*H(i)*S(i)];
end
end

% Initial Roe average state
for i = 1 : N - 1
rhoav(i) = sqrt(rho(i)*rho(i+1));
uav(i)   = (sqrt(rho(i))*u(i) + sqrt(rho(i+1))*u(i+1))/(sqrt(rho(i)) + sqrt(rho(i+1)));
Hav(i)   = (sqrt(rho(i))*H(i) + sqrt(rho(i+1))*H(i+1))/(sqrt(rho(i)) + sqrt(rho(i+1)));
aav(i)   = sqrt((k - 1)*(Hav(i) - 0.5*uav(i)^2));
pav(i)   = rhoav(i)*aav(i)^2/k;
end

% solution time loop
iteration = 0;
for T = 0 :l : Tf
iteration = iteration + 1;
% solution space loop
for i = 2 : N - 1
if i > 2 && i < N - 1
M1 = rright(k, uav(i-2), aav(i-2), rhoav(i-2))*(Q(:,i-1) - Q(:,i-2));
M2 = rright(k, uav(i-1), aav(i-1), rhoav(i-1))*(Q(:,i) - Q(:,i-1));
M3 = rright(k, uav(i), aav(i), rhoav(i))*(Q(:,i+1) - Q(:,i));
M4 = rright(k, uav(i+1), aav(i+1), rhoav(i+1))*(Q(:,i+2) - Q(:,i+1));
end

if i == 2
M2 = rright(k, uav(i-1), aav(i-1), rhoav(i-1))*(Q(:,i) - Q(:,i-1));
M3 = rright(k, uav(i), aav(i), rhoav(i))*(Q(:,i+1) - Q(:,i));
M1 = M3;
M4 = rright(k, uav(i+1), aav(i+1), rhoav(i+1))*(Q(:,i+2) - Q(:,i+1));
end

if i == N - 1
M1 = rright(k, uav(i-2), aav(i-2), rhoav(i-2))*(Q(:,i-1) - Q(:,i-2));
M2 = rright(k, uav(i-1), aav(i-1), rhoav(i-1))*(Q(:,i) - Q(:,i-1));
M3 = rright(k, uav(i), aav(i), rhoav(i))*(Q(:,i+1) - Q(:,i));
M4 = M2;
end

for j = 1: 3
% implement midmod function
if M2(j) > 0 && M3(j) > 0 && M4(j) > 0
minmod = max(0, min([M2(j), M3(j), M4(j)]));
else
minmod = min(0, max([M2(j), M3(j), M4(j)]));
end

%minmod = M2(j); % use when no need to using the limiter
N2(j) = M3(j) - minmod;

% implement midmod function
if M1(j) > 0 && M2(j) > 0 && M3(j) > 0
minmod = max(0, min([M1(j), M2(j), M3(j)]));
else
minmod = min(0, max([M1(j), M2(j), M3(j)]));
end

%minmod = M1(j); % use when no need to using the limiter
N1(j) = M2(j) - minmod;
end

lambda1 = [uav(i-1) 0 0; 0 uav(i-1) + aav(i-1) 0; 0 0 uav(i-1) - aav(i-1)]; % eigenvalue of flux Jacobian
lambda2 = [uav(i) 0 0; 0 uav(i) + aav(i) 0; 0 0 uav(i) - aav(i)]; % eigenvalue of flux Jacobian

maxl(i) = max(max(abs(lambda1)));

t1 = rleft(k, uav(i-1), aav(i-1), rhoav(i-1))*abs(lambda1)*N1;
t2 = rleft(k, uav(i), aav(i), rhoav(i))*abs(lambda2)*N2;

deltaF(:, i) = 0.5*( F(:,i+1) - F(:,i-1) - t2 + t1)/h;

end

Qold(:,:) = Q(:, :);
for i = 2 : N - 1
Q(:, i) = Q(:, i) - l*(deltaF(:,i) - [0; p(i)*dS(i); 0]); % update solution
end

% Q(:,N) = Qold(:, N-1); % use when close shock tube, and for case 4,
% 4.0 c

% use for real boundary conditions, 4.0b
%Q(2,1) = Qold(2,2);
%Q(1,N) = Qold(1,N-1);
%Q(2,N) = Qold(2,N-1);

% update parameters
for i = 1 : N
rho(i) = Q(1,i)/S(i);
u(i)   = Q(2,i)/(rho(i)*S(i));
% u(N)   = 0; % use when close shock tub
E(i)   = Q(3,i)/(rho(i)*S(i));
p(i)   = (k - 1)*rho(i)*(E(i) - 0.5*u(i)^2);
H(i)   = E(i) + p(i)/rho(i);
a(i)   = sqrt((k - 1)*(H(i) - 0.5*u(i)^2));
F(:,i) = [rho(i)*u(i)*S(i); (rho(i)*u(i)^2 + p(i))*S(i); rho(i)*u(i)*H(i)*S(i)];
end

% Update Roe average state
for i = 1 : N - 1
rhoav(i) = sqrt(rho(i)*rho(i+1));
uav(i)   = (sqrt(rho(i))*u(i) + sqrt(rho(i+1))*u(i+1))/(sqrt(rho(i)) + sqrt(rho(i+1)));
Hav(i)   = (sqrt(rho(i))*H(i) + sqrt(rho(i+1))*H(i+1))/(sqrt(rho(i)) + sqrt(rho(i+1)));
aav(i)   = sqrt((k - 1)*(Hav(i) - 0.5*uav(i)^2));
pav(i)   = rhoav(i)*aav(i)^2/k;
end

for i = 1 : N
Resx(i) = abs((Q(2,i)/S(i) - Qold(2,i)/S(i))./(Qold(2,i)/S(i)));
end

if iteration > 8000
break;
end

Time(iteration) = T;
Res(iteration) = max(Resx);

l = CFL*h/max(maxl); % update time step

end

% nozzle, shock tube and supersonic wind tunnel geometry
function [S, dS] = nozzle_geo(x, L, flag)

if isequal(flag, 'nozzle')
if x <= L/2
S = 1 + (3/2)*(1 - 0.2*x)^2;
dS = -(3/5)*(1 - 0.2*x);
else
S = 1 + (1/2)*(1 - 0.2*x)^2;
dS = -0.2*(1 - 0.2*x);
end
else
if isequal(flag, 'shock')
S = 1;
dS = 0;
end
end

if isequal(flag, 'wind')
if x <= L/4
S = 1 + (3/2)*(1 - 0.2*x)^2;
dS = -(3/5)*(1 - 0.2*x);
else
if x> L/4 && x <= L/2
S = 1 + (1/2)*(1 - 0.2*x)^2;
dS = -0.2*(1 - 0.2*x);
else
S = 1.5;
dS = 0;
end
end
end

% left eigenvectors of the Flux Jacobian
function R = rleft(k, u, a, rho)

c1 = 1;
c2 = rho/(sqrt(2)*a);
c3 = c2;
c4 = u;
c5 = c2*(u + a);
c6 = c2*(u - a);
c7 = u^2/2;
c8 = c2*(c7 + a^2/(k - 1) + a*u);
c9 = c2*(c7 + a^2/(k - 1) - a*u);

R = [ c1 c2 c3; c4 c5 c6; c7 c8 c9];

% right eigenvectors of the Flux Jacobian

function R = rright(k, u, a, rho)

c1 = 1 - (k - 1)*u^2/(2*a^2);
c2 = (k - 1)*u/a^2;
c3 = -(k - 1)/a^2;

t = 1/(sqrt(2)*rho*a);

c4 = t*(0.5*(k - 1)*u^2 - a*u);
c5 = t*( a - (k - 1)*u);
c6 = t*(k - 1);
c7 = t*(0.5*(k - 1)*u^2 + a*u);
c8 = -t*(a + (k - 1)*u);
c9 = t*(k - 1);

R = [ c1 c2 c3; c4 c5 c6; c7 c8 c9];

```