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

SIMPLE Method for Duct Flow

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

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   June 11, 2019, 10:37
Default SIMPLE Method for Duct Flow
  #1
New Member
 
Winston
Join Date: Feb 2019
Posts: 3
Rep Power: 4
obdwinston is on a distinguished road
I developed a MATLAB code to simulate a simple duct flow, based on the pressure correction technique with staggered grid introduced in J.D. Anderson's textbook on Computational Fluid Dynamics (page 247 - 264).

I managed to get it to work after some trial and error of the space and time steps, relaxation factor, and pressure inflow boundary condition. I am still not exactly sure why it works because I was merely trying my luck with various values. Can someone verify that the code is really doing its job, and not by some fluke that it works.

Thanks in advance!

--- MATLAB CODE ---

clear all, close all, clc

% initialise
m = 25;
n = 25;

p = zeros(m,n);
p_corr = zeros(m,n);
u = zeros(m,n-1);
v = zeros(m-1,n);

% constants
rho = 1.225; % density
nu = 1.5e-5; % kinematic viscosity
dx = 0.0075;
dy = 0.0005;
dt = 0.001;
a_p = .1; % relaxation factor

% boundary conditions
p(1:m,1) = .001; % inflow pressure boundary
p(1:m,2) = .001; % inflow pressure boundary

% SIMPLE algorithm
ctr = 0;
t = 0;
for iter = 1:10000
% x-momentum equation
for i = 2:m-1
for j = 2:n-2
A1 = -((u(i,j+1)*u(i,j+1) - u(i,j-1)*u(i,j-1))/(2*dx) + (u(i-1,j)*.5*(v(i-1,j)+v(i-1,j+1)) - u(i+1,j)*.5*(v(i,j)+v(i,j+1)))/(2*dy));
A2 = nu*((u(i,j+1)-2*u(i,j)+u(i,j-1))/(dx*dx) + (u(i-1,j)-2*u(i,j)+u(i+1,j))/(dy*dy));
A = A1 + A2;
u(i,j) = u(i,j) + A*dt - (dt/dx/rho)*(p(i,j+1)-p(i,j));
end
end

% zeroth-order extrapolation
u(2:m-1,1) = u(2:m-1,2);
u(2:m-1,n-1) = u(2:m-1,n-2);

% y-momentum equation
for i = 2:m-2
for j = 2:n-1
B1 = -((v(i-1,j)*v(i-1,j) - v(i+1,j)*v(i+1,j))/(2*dy) + (v(i,j+1)*.5*(u(i,j)+u(i+1,j)) - v(i,j-1)*.5*(u(i,j-1)+u(i+1,j-1)))/(2*dx));
B2 = nu*((v(i-1,j)-2*v(i,j)+v(i+1,j))/(dy*dy) + (v(i,j+1)-2*v(i,j)+v(i,j-1))/(dx*dx));
B = B1 + B2;
v(i,j) = v(i,j) + B*dt - (dt/dy/rho)*(p(i,j)-p(i+1,j));
end
end

% zeroth-order extrapolation
v(2:m-2,n) = v(2:m-2,n-1);

% pressure correction equation (Poisson equation)
a = 2.*(dt/(dx*dx) + dt/(dy*dy));
b = -dt/(dx*dx);
c = -dt/(dy*dy);

p_corr_diff = ones(m,n);
p_corr_prev = zeros(m,n);
while max(max(abs(p_corr_diff))) > 1e-6
p_corr_prev = p_corr;
for i = 2:m-1
for j = 2:n-1
d = (rho/dx)*(u(i,j) - u(i,j-1)) + (rho/dy)*(v(i-1,j) - v(i,j));
p_corr(i,j) = -(1/a)*(b*p_corr(i,j+1) + b*p_corr(i,j-1) + c*p_corr(i-1,j) + c*p_corr(i+1,j) + d);
end
end
p_corr_diff = p_corr - p_corr_prev;
end

% pressure correction
p = p + a_p*p_corr;
p(1,1:n) = p(2,1:n); % dp/dy = 0 at wall
p(m,1:n) = p(m-1,1:n); % dp/dy = 0 at wall

p(1:m,1) = .1; % inflow pressure boundary
p(1:m,2) = .1; % inflow pressure boundary

ctr = ctr + 1;
t = t + dt;
imagesc(u)
colorbar
str = sprintf('ctr = %d, t = %.3e s', ctr, t);
title(str)
pause(0.001)
end

Last edited by obdwinston; June 12, 2019 at 11:31.
obdwinston is offline   Reply With Quote

Reply

Tags
anderson, grid, pressure, simple, staggered

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
Issues on the simulation of high-speed compressible flow within turbomachinery dowlee OpenFOAM Running, Solving & CFD 9 February 5, 2021 10:15
Sample SIMPLE method CFD code entitledx Main CFD Forum 8 May 27, 2012 06:22
simple method for compressible flow saeed Main CFD Forum 0 April 24, 2003 02:20
SIMPLE method for compressible flow in pipe Reza Main CFD Forum 0 August 27, 2002 22:43
CFD method for High Re flow Toshiyuki Arima Main CFD Forum 3 May 29, 2001 11:18


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