CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > General Forums > Main CFD Forum

2D staggered SIMPLE MATLAB problme

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

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   November 17, 2021, 09:54
Question 2D staggered SIMPLE MATLAB problme
  #1
New Member
 
Gao Shangya
Join Date: Oct 2021
Posts: 18
Rep Power: 3
Harlotte is on a distinguished road
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/(ms)
%% 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
File Type: jpg gas.jpg (98.0 KB, 10 views)

Last edited by Harlotte; November 19, 2021 at 00:07.
Harlotte is offline   Reply With Quote

Reply

Tags
cfd, compressible air, matlab

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Lid Driven Cavity Code using SIMPLE in MATLAB deepmorzaria Main CFD Forum 0 April 23, 2020 16:02
How do I make steady state SIMPLE solver to transient? cfdnewb123 Main CFD Forum 5 April 22, 2020 13:49
SIMPLE Method for Duct Flow obdwinston Main CFD Forum 0 June 11, 2019 11:37
Using the Simple Algorithm for 2D Staggered Grid in Matlab jal121387 Main CFD Forum 7 February 19, 2019 05:38
SIMPLE algorithm Jonathan Castro Main CFD Forum 3 December 10, 1999 05:59


All times are GMT -4. The time now is 00:00.