Skill-Lync Launch Pad – Your Gateway to a Core Engineering Job by 2025! Only 1̶0̶0̶ +50 Seats Available.

03D 23H 28M 01S

Menu

Executive Programs

Workshops

Projects

Blogs

Careers

Placements

Student Reviews

For Business

More

Academic Training

Informative Articles

Find Jobs

We are Hiring!

All Courses

Choose a category

Mechanical

Electrical

Civil

Computer Science

Electronics

Offline Program

All Courses

All Courses

logo

CHOOSE A CATEGORY

Mechanical

Electrical

Civil

Computer Science

Electronics

Offline Program

Top Job Leading Courses

Automotive

CFD

FEA

Design

MBD

Med Tech

Courses by Software

Design

Solver

Automation

Vehicle Dynamics

CFD Solver

Preprocessor

Courses by Semester

First Year

Second Year

Third Year

Fourth Year

Courses by Domain

Automotive

CFD

Design

FEA

Tool-focused Courses

Design

Solver

Automation

Preprocessor

CFD Solver

Vehicle Dynamics

Machine learning

Machine Learning and AI

POPULAR COURSES

coursePost Graduate Program in Hybrid Electric Vehicle Design and Analysis
coursePost Graduate Program in Computational Fluid Dynamics
coursePost Graduate Program in CAD
coursePost Graduate Program in CAE
coursePost Graduate Program in Manufacturing Design
coursePost Graduate Program in Computational Design and Pre-processing
coursePost Graduate Program in Complete Passenger Car Design & Product Development
Executive Programs
Workshops
For Business

Success Stories

Placements

Student Reviews

More

Projects

Blogs

Academic Training

Find Jobs

Informative Articles

We're Hiring!

phone+91 9342691281Log in
  1. Home/
  2. Amith Ganta/
  3. Lid Driven Cavity Flow using SIMPLE Algorithm and Artificial Compressibility method in MATLAB (Steady Flow)

Lid Driven Cavity Flow using SIMPLE Algorithm and Artificial Compressibility method in MATLAB (Steady Flow)

Computational Fluid Dynamics :  The term computational refers to one of the pathways to solve the fluid flow equations. the other possible pathways are experimental fluid dynamics and theoretical fluid dynamics. In experimental fluid dynamics, various experimental techniques are employed to determine the fluid flow…

  • CFD

  • Amith Ganta

    updated on 21 Jun 2021

Computational Fluid Dynamics : 

The term computational refers to one of the pathways to solve the fluid flow equations. the other possible pathways are experimental fluid dynamics and theoretical fluid dynamics. In experimental fluid dynamics, various experimental techniques are employed to determine the fluid flow phenomena. Theoretical fluid dynamics also refers to analytical fluid dynamics where we mathematically solve the governing equations to achieve fluid flow solution. In computational fluid dynamics, various numerical techniques are employed to simulate the flow of fluid using computers. 

Fluid: Any material that cannot resist the flow of shear stress is called fluid. All Liquids and Gases are referred to as Fluids. 

Dynamics: to determine both the kinematics parameters such as velocity, acceleration and the effect these parameters causes such as Force (Lift, Drag, Pressure forces, etc)

Computational fluid dynamics allows estimation of all flow-related phenomena using computers. This helps solve complex fluid flow phenomena and results can be obtained in no time with less cost than experimental technique. 

 

Governing Equations : 

2D incompressible flow:

                                                             Conservation of Mass :

                                                                ∂u∂x+∂v∂y=0∂u∂x+∂v∂y=0

                                                               Conservation of Momentum :

                                 ∂u∂t+u∂u∂x+v∂u∂y=−1ρ∂p∂x+ϑ(∂2u∂x2+∂2u∂y2)+Fx∂u∂t+u∂u∂x+v∂u∂y=-1ρ∂p∂x+ϑ(∂2u∂x2+∂2u∂y2)+Fx

                                 ∂v∂t+u∂v∂x+v∂v∂y=−1ρ∂p∂y+ϑ(∂2v∂x2+∂2v∂y2)+Fy∂v∂t+u∂v∂x+v∂v∂y=-1ρ∂p∂y+ϑ(∂2v∂x2+∂2v∂y2)+Fy

 

In CFD, these differential equations are converted into Algebraic equations so that they can easily be manipulated. Finite differencing methods are one of the ways to convert them. 

 

The fluid flow is governed using 

 

1. Continuity Equation :

∂u∂x+∂v∂y=0∂u∂x+∂v∂y=0

 

2. Momentum Equation:

 x - momentum = u∂u∂x+v∂u∂y=−∂p∂x+(1Re)(∂2u∂x2+∂2u∂y2)u∂u∂x+v∂u∂y=-∂p∂x+(1Re)(∂2u∂x2+∂2u∂y2)

 y - momentum = u∂v∂x+v∂v∂y=−∂p∂y+(1Re)(∂2v∂x2+∂2v∂y2)u∂v∂x+v∂v∂y=-∂p∂y+(1Re)(∂2v∂x2+∂2v∂y2) 

 To solve the problem in Artificial compressibility method, artificial compressibility is introduced in to the continuity equation. Using that pressure field can be solved. 

 

Steady-state incompressible Navier stokes equations

1δ(∂p∂t)+1δ(∂p∂t)+∂u∂x+∂v∂y=0∂u∂x+∂v∂y=0

x - momentum = ∂u∂t∂u∂t + u∂u∂x+v∂u∂y=−∂p∂x+(1Re)(∂2u∂x2+∂2u∂y2)u∂u∂x+v∂u∂y=-∂p∂x+(1Re)(∂2u∂x2+∂2u∂y2)

 

 

SIMPLE Algorithm : 

clear all
close all
clc

%% Defining the problem domain
n_points = 51; % no_of_points
dom_length = 1;
h = dom_length/(n_points-1);
x = 0:h:dom_length; %X domain span
y = 0:h:dom_length; %Y domain span
Re = 1000; %Reynolds number
nu = 1/Re;
% Under-relaxation factors
alpha = 0.01;
alpha_p = 0.8;
%% Initializing the variables
%Final collocated variables
u_final(n_points,n_points)=0;
v_final(n_points,n_points)=0;
p_final(n_points,n_points)=1;
u_final(1,:) = 1;

%Staggered variables
u(n_points+1,n_points)=0;
u_star(n_points+1,n_points)=0;
d_e(n_points+1,n_points)=0;
v(n_points,n_points+1)=0;
v_star(n_points,n_points+1)=0;
d_n(n_points,n_points+1)=0;
p(n_points+1,n_points+1)=1;
p_star(n_points+1,n_points+1)=1;
pc(n_points+1,n_points+1)=0;
b(n_points+1,n_points+1)=0;
u(1,:)=2;

u_new(n_points+1,n_points)=0;
v_new(n_points,n_points+1)=0;
p_new(n_points+1,n_points+1)=1;
u_new(1,:)=2;

%% Solving the governing equations
error = 1;
iterations = 0;
error_req = 1e-7; %final required error residual
figure(1); %for error monitoring

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*u_E*h + nu;
            a_W = 0.5*u_W*h + nu;
            a_N = -0.5*v_N*h + nu;
            a_S = 0.5*v_S*h + nu;
            
            a_e = 0.5*u_E*h - 0.5*u_W*h + 0.5*v_N*h - 0.5*v_S*h + 4*nu;
            
            A_e = -h;
            d_e(i,j) = A_e/a_e;
            
            u_star_mid = (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));
            u_star(i,j) = (1-alpha)*u(i,j) + alpha*u_star_mid;
        end
    end
    
    % x-momentum eq. - Boundary
    u_star(1,:) = 2 - u_star(2,:);
    u_star(n_points + 1,:) = -u_star(n_points,:);
    u_star(2:n_points,1) = 0;
    u_star(2:n_points,n_points) = 0;
    
    % 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.5*u_E*h + nu;
            a_W = 0.5*u_W*h + nu;
            a_N = -0.5*v_N*h + nu;
            a_S = 0.5*v_S*h + nu;
            
            a_n = 0.5*u_E*h - 0.5*u_W*h + 0.5*v_N*h - 0.5*v_S*h + 4*nu;
            
            A_n = -h;
            d_n(i,j) = A_n/a_n;
            
            v_star_mid = (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));
            v_star(i,j) = (1-alpha)*v(i,j) + alpha*v_star_mid;
        end
    end
    
    % y-momentum eq. - Boundary
    v_star(:,1) = -v_star(:,2);
    v_star(:,n_points + 1) = -v_star(:,n_points);
    v_star(1,2:n_points) = 0;
    v_star(n_points,2:n_points) = 0;
    
    % Zeroing the corrections to begin with
    pc(1:n_points+1,1:n_points+1)=0;
    
    % Continuity equation a.k.a. pressure correction - Interior
    for i = 2:n_points
        for j = 2:n_points
            a_E = -d_e(i,j)*h;
            a_W = -d_e(i,j-1)*h;
            a_N = -d_n(i-1,j)*h;
            a_S = -d_n(i,j)*h;
            a_P = a_E + a_W + a_N + a_S;
            b(i,j) = -(u_star(i,j) - u_star(i,j-1))*h + (v_star(i,j) - v_star(i-1,j))*h;
            
            pc(i,j) = (a_E*pc(i,j+1) + a_W*pc(i,j-1) + a_N*pc(i-1,j) + a_S*pc(i+1,j) + b(i,j))/a_P;
        end
    end
    
    % Correcting the pressure field
    for i = 2:n_points
        for j = 2:n_points
            p_new(i,j) = p(i,j) + alpha_p*pc(i,j);
        end
    end
    
    % Continuity eq. - Boundary
    p_new(1,:) = p_new(2,:);
    p_new(n_points + 1,:) = p_new(n_points,:);
    p_new(:,1) = p_new(:,2);
    p_new(:,n_points + 1) = p_new(:,n_points);
    
    % Correcting the velocities
    for i = 2:n_points
        for j = 2:n_points - 1
            u_new(i,j) = u_star(i,j) + d_e(i,j)*(pc(i,j+1) - pc(i,j));
        end
    end
    
    % x-momentum eq. - Boundary
    u_new(1,:) = 2 - u_new(2,:);
    u_new(n_points + 1,:) = -u_new(n_points,:);
    u_new(2:n_points,1) = 0;
    u_new(2:n_points,n_points) = 0;
    
    for i = 2:n_points - 1
        for j = 2:n_points
            v_new(i,j) = v_star(i,j) + d_n(i,j)*(pc(i,j) - pc(i+1,j));
        end
    end
    
    % y-momentum eq. - Boundary
    v_new(:,1) = -v_new(:,2);
    v_new(:,n_points + 1) = -v_new(:,n_points);
    v_new(1,2:n_points) = 0;
    v_new(n_points,2:n_points) = 0;
            
    
    % Continuity residual as error measure
    error = 0;
    for i = 2:n_points
        for j = 2:n_points
            error = error + abs(b(i,j));
        end
    end
    
    % Error monitoring after every few timesteps
    if(rem(iterations, 1000)) == 0
       figure(1);
       semilogy(iterations, error, '-ko')
       hold on
       xlabel('Iterations')
       ylabel('Residual Error')
    end
    
    % Finishing the iteration
    u = u_new;
    v = v_new;
    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));
        p_final(i,j) = 0.25*(p(i,j) + p(i,j+1) + p(i+1,j) + p(i+1,j+1));
    end
end

%% Centreline u validation - Comparison with benchmark
% figure(11);
% plot(u_final(:,(n_points+1)/2),1-y, 'LineWidth', 1)
% 
% tab_data = xlsread('plot_u_y_Ghia100.csv', 'A2:B18');
% y_ghia = tab_data(:,1);
% u_ghia = tab_data(:,2);
% 
% figure(11); hold on
% plot(u_ghia, y_ghia, 'o', 'LineWidth', 1)
% xlabel('u')
% ylabel('y')
% legend('Numerical', 'Benchmark', 'location', 'southeast')

%% 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,u_final, 21, 'LineStyle', 'none')
colorbar
colormap('jet')
xlabel('x')
ylabel('y')

figure(22);
hold on
quiver(X, Y, u_final, v_final, 5, 'k')

 

 

 

 

 

Artificial Compressibility Method :

clear all
close all
clc

%% Defining the Domain Problem

n_points = 51 ; % Number of points
dom_length = 1; % Domain length
h = dom_length/(n_points - 1); % Grid spacing
x = 0:h:dom_length; % X - coordinates
y = 0:h:dom_length; % Y - coordinates
dt = 0.001; % Time advancement
Re = 100; % Reynolds number
delta = 1.5; % Artificial compressibility


%% Initializing the variables

% Final collocated variables

u_final(n_points,n_points) = 0;
v_final(n_points,n_points) = 0;
p_final(n_points,n_points) = 0;
u_final(1,:) = 1;

% Staggered variables
u(n_points+1, n_points) = 0;
v(n_points, n_points+1) = 0;
p(n_points+1 , n_points+1) = 1;
u(1,:) = 2;


u_new(n_points+1, n_points) = 0;
v_new(n_points, n_points+1) = 0;
p_new(n_points+1 , n_points+1) = 1;
u_new(1,:) = 2;

%% Solving the Governing Equations
error = 1;
iterations = 0;
error_req = 1e-6; % final required error residual
figure(1); % for error monitoring

while error > error_req
    % x-momentum eq. - Interior
    for i = 2:n_points
        for j = 2:n_points - 1
            pressure = -(p(i,j+1) - p(i,j))/h;
            diffusion = (1/Re)*((u(i+1,j) - 2*u(i,j) + u(i-1,j))/(h*h) + (u(i,j+1) - 2*u(i,j) + u(i,j-1))/(h*h));
            advection_x = ((0.5*(u(i,j)+u(i,j+1)))^2 - (0.5*(u(i,j)+u(i,j-1)))^2)/h;
            advection_y = ((0.25*(u(i,j)+u(i-1,j))*(v(i-1,j)+v(i-1,j+1))) - (0.25*(u(i,j)+u(i+1,j))*(v(i,j)+v(i,j+1))))/h;
            u_new(i,j) = u(i,j) + dt*(diffusion - advection_x - advection_y + pressure);
        end
    end
    
    % x-momentum eq. - Boundary
    u_new(1,:) = 2 - u_new(2,:);
    u_new(n_points + 1,:) = -u_new(n_points,:);
    u_new(2:n_points,1) = 0;
    u_new(2:n_points,n_points) = 0;    
    
    % y-momentum eq. - Interior
    for i = 2:n_points - 1
        for j = 2:n_points
            pressure = -(p(i,j) - p(i+1,j))/h;
            diffusion = (1/Re)*((v(i+1,j) - 2*v(i,j) + v(i-1,j))/(h*h) + (v(i,j+1) - 2*v(i,j) + v(i,j-1))/(h*h));
            advection_y = ((0.5*(v(i,j)+v(i-1,j)))^2 - (0.5*(v(i,j)+v(i+1,j)))^2)/h;
            advection_x = ((0.25*(u(i,j)+u(i+1,j))*(v(i,j)+v(i,j+1))) - (0.25*(u(i,j-1)+u(i+1,j-1))*(v(i,j)+v(i,j-1))))/h;
            v_new(i,j) = v(i,j) + dt*(diffusion - advection_x - advection_y + pressure);
        end
    end
    
    % y-momentum eq. - Boundary
    v_new(:,1) = -v_new(:,2);
    v_new(:,n_points + 1) = -v_new(:,n_points);
    v_new(1,2:n_points) = 0;
    v_new(n_points,2:n_points) = 0;
    
    % Continuity eq. - Interior
    for i = 2:n_points
        for j = 2:n_points
            p_new(i,j) = p(i,j) - delta*dt*(u(i,j) - u(i,j-1) + v(i-1,j) - v(i,j))/h;
        end
    end
    
    % Continuity eq. - Boundary
    p_new(1,:) = p_new(2,:);
    p_new(n_points + 1,:) = p_new(n_points,:);
    p_new(:,1) = p_new(:,2);
    p_new(:,n_points + 1) = p_new(:,n_points);
    
    % 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
    
    % Error monitoring after every few timesteps
    if(rem(iterations, 5000)) == 0
       figure(1);
       semilogy(iterations, error, '-ko')
       hold on
       xlabel('Iterations')
       ylabel('Residual Error')
    end
    
    % Finishing the iteration
    u = u_new;
    v = v_new;
    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));
        p_final(i,j) = 0.25*(p(i,j) + p(i,j+1) + p(i+1,j) + p(i+1,j+1));
    end
end

%% Centreline u validation - Comparison with benchmark
figure(11);
plot(u_final(:,(n_points+1)/2),1-y, 'LineWidth', 1)

tab_data = xlsread('plot_u_y_Ghia100.csv', 'A2:B18');
y_ghia = tab_data(:,1);
u_ghia = tab_data(:,2);

figure(11); hold on
plot(u_ghia, y_ghia, 'o', 'LineWidth', 1)
xlabel('u')
ylabel('y')
legend('Numerical', 'Benchmark', 'location', 'southeast')

%% 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')
colorbar
colormap('jet')
xlabel('x')
ylabel('y')

figure(22);
hold on
quiver(X, Y, u_final, v_final, 5, 'k')

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Leave a comment

Thanks for choosing to leave a comment. Please keep in mind that all the comments are moderated as per our comment policy, and your email will not be published for privacy reasons. Please leave a personal & meaningful conversation.

Please  login to add a comment

Other comments...

No comments yet!
Be the first to add a comment

Read more Projects by Amith Ganta (43)

Quater car modelling using Matlab

Objective:

                                                                                  Quarter Car model …

calendar

25 Sep 2023 03:07 PM IST

    Read more

    Pacjeka combined slip

    Objective:

                                            Pacjeka Magic Tyre model   Pacjeka Magic Tyre model was developed by Prof.Hans Bastiaan Pacejka and it is widely used to represent empirically and interpolate previously measured…

    calendar

    05 Jun 2023 02:04 PM IST

      Read more

      Longitudinal and Lateral brush tyre model using Matlab

      Objective:

                                                                       Tyre Modeling A tyre model is a mathematical or an empirical representation of the mechanism…

      calendar

      16 May 2023 07:14 PM IST

        Read more

        Single cylinder SI engine modelling using GT-Suite

        Objective:

                                                                     Introduction to GT Suite Effective utilization of simulation software with human intelligence is required…

        calendar

        03 Apr 2022 12:50 AM IST

        • CFD
        • GT-POWER
        Read more

        Schedule a counselling session

        Please enter your name
        Please enter a valid email
        Please enter a valid number

        Related Courses

        coursecardcoursetype

        Post Graduate Program in CFD Solver Development

        4.8

        119 Hours of Content

        coursecard

        Introduction to OpenFOAM Development

        4.9

        18 Hours of Content

        coursecardcoursetype

        Post Graduate Program in Battery Technology for Mechanical Engineers

        4.8

        81 Hours of Content

        coursecardcoursetype

        Post Graduate Program in Automation & Pre-Processing for FEA & CFD Analysis

        4.7

        81 Hours of Content

        coursecardcoursetype

        Post Graduate Program in Hybrid Electric Vehicle Design and Analysis

        4.8

        343 Hours of Content

        Schedule a counselling session

        Please enter your name
        Please enter a valid email
        Please enter a valid number

        logo

        Skill-Lync offers industry relevant advanced engineering courses for engineering students by partnering with industry experts.

        https://d27yxarlh48w6q.cloudfront.net/web/v1/images/facebook.svghttps://d27yxarlh48w6q.cloudfront.net/web/v1/images/insta.svghttps://d27yxarlh48w6q.cloudfront.net/web/v1/images/twitter.svghttps://d27yxarlh48w6q.cloudfront.net/web/v1/images/youtube.svghttps://d27yxarlh48w6q.cloudfront.net/web/v1/images/linkedin.svg

        Our Company

        News & EventsBlogCareersGrievance RedressalSkill-Lync ReviewsTermsPrivacy PolicyBecome an Affiliate
        map
        EpowerX Learning Technologies Pvt Ltd.
        4th Floor, BLOCK-B, Velachery - Tambaram Main Rd, Ram Nagar South, Madipakkam, Chennai, Tamil Nadu 600042.
        mail
        info@skill-lync.com
        mail
        ITgrievance@skill-lync.com

        Top Individual Courses

        Computational Combustion Using Python and CanteraIntroduction to Physical Modeling using SimscapeIntroduction to Structural Analysis using ANSYS WorkbenchIntroduction to Structural Analysis using ANSYS Workbench

        Top PG Programs

        Post Graduate Program in Hybrid Electric Vehicle Design and AnalysisPost Graduate Program in Computational Fluid DynamicsPost Graduate Program in CADPost Graduate Program in Electric Vehicle Design & Development

        Skill-Lync Plus

        Executive Program in Electric Vehicle Embedded SoftwareExecutive Program in Electric Vehicle DesignExecutive Program in Cybersecurity

        Trending Blogs

        Heat Transfer Principles in Energy-Efficient Refrigerators and Air Conditioners Advanced Modeling and Result Visualization in Simscape Exploring Simulink and Library Browser in Simscape Advanced Simulink Tools and Libraries in SimscapeExploring Simulink Basics in Simscape

        © 2025 Skill-Lync Inc. All Rights Reserved.

                    Do You Want To Showcase Your Technical Skills?
                    Sign-Up for our projects.