# 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 Thread Tools Search this Thread Show Printable Version Email this Page Search this Thread: Advanced Search Display Modes Linear Mode Switch to Hybrid Mode Switch to Threaded Mode 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 OffTrackbacks are Off Pingbacks are On Refbacks are On Forum Rules Similar Threads Thread Thread Starter Forum Replies Last Post dowlee OpenFOAM Running, Solving & CFD 11 August 6, 2021 06:40 entitledx Main CFD Forum 8 May 27, 2012 06:22 saeed Main CFD Forum 0 April 24, 2003 02:20 Reza Main CFD Forum 0 August 27, 2002 22:43 Toshiyuki Arima Main CFD Forum 3 May 29, 2001 11:18

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