About us    Site map  
       Site map   
Partenov CFD  


   Home   

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');


Published with MATLAB® R2017a
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');

Published with MATLAB® R2017a
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


Published with MATLAB® R2017a
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


Published with MATLAB® R2017a
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


Published with MATLAB® R2017a
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];


Published with MATLAB® R2017a