# 2D staggered SIMPLE MATLAB problme

 Register Blogs Members List Search Today's Posts Mark Forums Read

November 17, 2021, 08:54
2D staggered SIMPLE MATLAB problme
#1
New Member

Gao Shangya
Join Date: Oct 2021
Posts: 18
Rep Power: 4
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:
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')
Attached Images
 gas.jpg (98.0 KB, 13 views)

Last edited by Harlotte; November 18, 2021 at 23:07.

 Tags cfd, compressible air, matlab