# SIMPLE Method for Duct Flow

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

 June 11, 2019, 10:37 SIMPLE Method for Duct Flow #1 New Member   Winston Join Date: Feb 2019 Posts: 3 Rep Power: 6 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.

 Tags anderson, grid, pressure, simple, staggered