[Harlotte]

Hello everyone,

I am working on my masters thesis which is about writing a code in Matlab to describe the simulation of a gas exchange in an hydrogen combustion engine.And now I am trying to writing a code( compressible and dimensional).

Allowing air to flow into the cylinder through the inlet valve, at which point the piston is stationary. I used governing equations (below). I assume that the temperature is constant and ignore the energy equation.
I want to simulate the change in velocity, pressure and density of the gas in the cylinder as the air enters the cylinder.

I'm writing my code and here is what I have done:
-have discretized the continuity and momentum equations using FVM and plan to close the system with equation of state
-use SIMPLE update and correct the pressure and velocity

But when I run the code, I cannot see the result, I try to change the values,it doesn't work....

So I hope you could take a look at my code and give me some advice or help me point out/correct errors, I'd really appreciate your help!!!

Here is my code:

clear all
close all
clc

%% Defining the problem domain
n_points = 51; % no_of_points
dom_length = 0.5;% length in m
h = dom_length/(n_points-1);
x = 0:h:dom_length; %X domain span
y = 0:h:dom_length; %Y domain span
T = 300; %Temperature in K
R = 287.05; % air gas constant in J/Kg K
mu = 1.81e-5 % air viscosity in kg/(m·s)
%% Initializing the variables
%Final  variables
u_final(n_points,n_points)=0;
v_final(n_points,n_points)=0;
p_final(1:n_points,1:n_points)=101325; %pressure in Pa
rho(1:n_points,1:n_points) = 1.25; % density in kg/m³
u_final(1,5:10)=-1; %velocity in m/s
u_final(1,16:21)=1;
v_final(1,5:10)=-1;
v_final(1,16:21)=-1;

%Staggered variables
u(n_points+1,n_points)=0;
v(n_points,n_points+1)=0;
p(1:n_points+1,1:n_points+1)=101325;%pressure in Pa
rho(1:n_points+1,1:n_points+1) = 1.25; % density in kg/m³
u(1,5:10)=-2;
u(1,16:21)=2;
v(1,5:10)=-2;
v(1,16:21)=-2;

u_new(n_points+1,n_points)=0;
v_new(n_points,n_points+1)=0;
p_new(1:n_points+1,1:n_points+1)=1; %pressure in Pa
rho_new(1:n_points+1,1:n_points+1) = 1.25; % density in kg/m³
u_new(1,5:10)=-2;
u_new(1,16:21)=2;
v_new(1,5:10)=-2;
v_new(1,16:21)=-2;
%% Solving the governing equations
error = 1;
iterations = 0;
error_req = 1e-7; %final required error residual
figure(1); %for error monitoring

% discretize x-momentum using FVM
while error > error_req
    % x-momentum eq. - Interior
      for i = 2:n_points
        for j = 2:n_points - 1
            u_E = 0.5*(u(i,j) + u(i,j+1));
            u_W = 0.5*(u(i,j) + u(i,j-1));
            v_N = 0.5*(v(i-1,j) + v(i-1,j+1));
            v_S = 0.5*(v(i,j) + v(i,j+1));
            
            a_E = -0.5*rho(i,j+1)*u_E*h + mu ;
            a_W = 0.5*rho(i,j)*u_W*h + mu;
            a_N = -0.125*(rho(i,j)+rho(i,j+1)+rho(i-1,j+1)+rho(i-1,j))*v_N*h + mu;
            a_S = 0.125*(rho(i,j)+rho(i,j+1)+rho(i+1,j+1)+rho(i+1,j))*v_S*h +mu;
            
            a_e = 0.5*rho(i,j+1)*u_E*h-0.5*rho(i,j)*u_W*h+0.125*(rho(i,j)+rho(i,j+1)+rho(i-1,j+1)+rho(i-1,j))- 0.125*(rho(i,j)+rho(i,j+1)+rho(i+1,j+1)+rho(i+1,j))+4*mu ;
            
            A_e = -h;
            d_e(i,j) = A_e/a_e;
            
            u_new(i,j) = (a_E*u(i,j+1) + a_W*u(i,j-1) + a_N*u(i-1,j) + a_S*u(i+1,j))/a_e + d_e(i,j)*(p(i,j+1) - p(i,j));
        end
    end
    % x-momentum eq. - Boundary condition
    u_new(1,5:10) = -2-u_new(2,5:10);
    u_new(1,16:21)=2-u_new(2,5:10);
    u_new(n_points + 1,:) = -u_new(n_points,:);
    u_new(2:n_points,1) = 0;
    u_new(2:n_points,n_points) = 0;    
    u_new(1,1:4)=0;
    u_new(1,22:51)=0;
    u_new(1,11:15)=0;
    
    % discretize y-momentum using FVM
    % y-momentum eq. - Interior
    for i = 2:n_points - 1
        for j = 2:n_points
            u_E = 0.5*(u(i,j) + u(i+1,j));
            u_W = 0.5*(u(i,j-1) + u(i+1,j-1));
            v_N = 0.5*(v(i-1,j) + v(i,j));
            v_S = 0.5*(v(i,j) + v(i+1,j));
            
            a_E = -0.125*(rho(i,j)+rho(i,j+1)+rho(i+1,j)+rho(i+1,j+1))*u_E*h + mu ;
            a_W = 0.125*(rho(i,j)+rho(i,j-1)+rho(i+1,j)+rho(i+1,j-1))*u_W*h + mu;
            a_N = -0.5*rho(i,j)*v_N*h + mu;
            a_S = 0.5*rho(i,j+1)*v_S*h +mu;
            
            a_n = 0.125*(rho(i,j)+rho(i,j+1)+rho(i+1,j)+rho(i+1,j+1))*u_E*h-0.125*(rho(i,j)+rho(i,j-1)+rho(i+1,j)+rho(i+1,j-1))*u_W*h+0.5*rho(i,j)*v_N*h - 0.5*rho(i,j+1)*v_S*h+4*mu;
             
            A_n = -h;
            d_n(i,j) = A_n/a_n;
            
            v_new(i,j) = (a_E*v(i,j+1) + a_W*v(i,j-1) + a_N*v(i-1,j) + a_S*v(i+1,j))/a_n + d_n(i,j)*(p(i,j) - p(i+1,j));
        end
    end
    
    % y-momentum eq. - Boundary conditions
    v_new(:,1) = -v_new(:,2);
    v_new(:,n_points + 1) = -v_new(:,n_points);
    v_new(n_points,2:n_points) = 0;
    v_new(1,5:10) = -2-v_new(2,5:10);
    v_new(1,16:21)=-2-v_new(2,16:21);
    v_new(1,1:4)=0;
    v_new(1,11:15)=0;
    v_new(1,22:51)=0; 
    
    % discretize density using FVM
    % Continuity eq. - Interior- problem with the expression of density
    for i = 2:n_points
        for j = 2:n_points
             (rho_new(i,j+1)-rho(i,j-1))*(u_new(i,j)-u(i,j-1))==(v_new(i,j)-v_new(i-1,j))*(rho(i-1,j-rho(i,j+1));
              
        end
    end
    
    % Continuity eq. - Boundary conditions
    rho_new(1,:) = rho_new(2,:);
    rho_new(n_points + 1,:) = rho_new(n_points,:);
    rho_new(:,1) = rho_new(:,2);
    rho_new(:,n_points + 1) = rho_new(:,n_points);
    
    % Calculate Continuity residual as error measure
    error = 0;
    for i = 2:n_points - 1
        for j = 2:n_points - 1
            error = error + abs((u_new(i,j) - u_new(i,j-1) + v_new(i-1,j) - v_new(i,j))/h);
        end
    end
     
    %  Iteration
    u = u_new;
    v = v_new;
    rho = rho_new;
        for i = 1:n_points + 1
        for j = 1:n_points + 1
            p_new(i,j)=rho_new(i,j)*R*T;
        end
        end
    p = p_new;
    iterations = iterations + 1;
end

% After the converged solution, we map the staggered variables to
% collocated variables
for i = 1:n_points
    for j = 1:n_points
        u_final(i,j) = 0.5*(u(i,j) + u(i+1,j));
        v_final(i,j) = 0.5*(v(i,j) + v(i,j+1));
        rho_final(i,j) = 0.25*(rho(i,j) + rho(i,j+1) + rho(i+1,j) + rho(i+1,j+1));
    end
end



%% Contour and vector visuals.
x_dom = ((1:n_points)-1).*h;
y_dom = 1-((1:n_points)-1).*h;
[X,Y] = meshgrid(x_dom,y_dom);
figure(21);
contourf(X,Y,v_final, 21, 'LineStyle', 'none')
c = colorbar;
c.Label.String = 'v final';
colormap('jet')
xlabel('x')
ylabel('y')
title('Velocity in y axis')

figure(22);
contourf(X,Y,u_final, 21, 'LineStyle', 'none')
colorbar
colormap('jet')
xlabel('x')
ylabel('y')
c1 = colorbar;
c1.Label.String = 'u final';
title('Velocity in x axis')

figure(23);
contourf(X,Y,p_final, 21, 'LineStyle', 'none')
colorbar
colormap('jet')
xlabel('x')
ylabel('y')
c2 = colorbar;
c2.Label.String = 'p final';
title('Pressure in Cylinder')

figure(24);
hold on
quiver(X, Y, u_final, v_final, 5, 'k')
axis equal
title('Velocity in Cylinder')
xlabel('u final')
ylabel('v final')