2D Couette Flow ( SIMPLE Algorithm )

MATLAB Script :

%Pressure correction method applied to incompressible flow between two

%parallel plates (i.e. Couette flow 

clc; clear all; close all;

% GEOMETRY CONSTANTS 

% Matlab reads matrices as (rows, columns)

L = 0.5; % length of domain, ft

H = 0.01; % height of domain, ft 

Nxp = 21; % # of points in x direction at which pressure is resolved

Nyp = 11; % # of points in y direction at which pressure is resolved

Nxu = Nxp + 1; % # of points in x direction at which u velocity is resolved

Nyu = Nyp; % # of points in y direction at which u velocity is resolved

Nxv = Nxp + 2; % # of points in x direction at which v velocity is resolved

Nyv = Nyp + 1; % # of points in y direction at which v velocity is resolved

dx = L/(Nxp-1); % in x direction

dy = H/(Nyp-1); % in y direction

dt = 1e-3; % time resolution, seconds

alpha_p = 0.1; % under-relaxation factor for pressure

alpha_u = 0; % under-relaxation factor for u velocity

alpha_v = 0; % under-relaxation factor for v velocity

% CREATE STAGGERED RECTANGULAR GRID

[xp yp] = meshgrid(0:dx:L,0:dy:H);

[xu yu] = meshgrid(-0.5*dx:dx:L+0.5*dx,0:dy:H);

[xv yv] = meshgrid(-1*dx:dx:L+1*dx,-0.5*dy:dy:H+0.5*dy);

plot(xp,yp,'*r'); hold on; plot(xu,yu,'ob'); hold on;

plot(xv,yv,'xm'); hold on;

% FLUID CONSTANTS

rho = 0.002377; % density of air, slug/ft^3

mu = 3.726e-7; % dynamic viscosity of air, lbf s/ft^2

% INITIAL CONDITIONS

pstar = zeros(Nyp,Nxp);

ustar = zeros(Nyu,Nxu); ustar(1,:) = 1;

vstar = zeros(Nyv,Nxv); vstar((Nyv-5)+1,15) = 0.5; % velocity spike at (15,5), in ft/s

rho_ustar = rho.*ustar;

rho_vstar = rho.*vstar;

pprime = zeros(Nyp,Nxp);

pprime1 = zeros(Nyp,Nxp);

uprime = zeros(Nyu,Nxu);

vprime = zeros(Nyv,Nxv);

A = zeros(Nyu,Nxu);

B = zeros(Nyv,Nxv);

% BOUNDARY CONDITIONS

v = vstar;

u = ustar; % horizontal velocity, ft/s

% PRESSURE CORRECTION METHOD

b = -dt/(dx^2);

c = -dt/(dy^2);

a = -2*(b + c);

d = zeros(Nyp,Nxp); % mass flow term in continuity equation

conv = ones(Nyp,Nxp);   

figure; plot(0,d(7,15),'-o'); hold on; grid on;

%%

count = 0;

for m = 1:5

count = count+1;   

loop = [4 20 50 150 300];

for iter = 1 : loop(m),

    for i = 1:Nxp-1, % column, 1 to 20

        for j = 2:Nyu-1, % row, 2 to 10

%      for i = 1:Nxu-2, % column

%         for j = 2:Nyu-1, % row

        k = Nyp-((Nyp-j)+1)+1; % 2 to 10

        kv = Nyv-((Nyv-j))+1; % 3 to 11

        ku = Nyu-((Nyu-j)+1)+1; % 2 to 10

        vbar1 = 0.5*(v(kv-1,i+1)+v(kv-1,i+2));

        vbar2 = 0.5*(v(kv,i+1)+v(kv,i+2));

        A(ku,i+1) = -((rho*(u(ku,i+2)^2)-rho*(u(ku,i)^2))/(2*dx)...

             + ((rho*u(ku-1,i+1)*vbar1)-(rho*u(ku+1,i+1)*vbar2))/(2*dy))...

             + mu*((u(ku,i+2)-2*u(ku,i+1)+u(ku,i))/(dx^2) ...

             + (u(ku-1,i+1)-2*u(ku,i+1)+u(ku+1,i+1))/(dy^2));

        rho_ustar(ku,i+1) = rho_ustar(ku,i+1) + A(ku,i+1)*dt ...

                       - (dt/dx)*(pstar(k,i+1) - pstar(k,i));

        ustar(ku,i+1) = rho_ustar(ku,i+1)/rho;

        u(ku,i+1) = ustar(ku,i+1);

        end

    end

ustar(:,1) = ustar(:,2);

ustar(:,Nxu) = ustar(:,Nxu-1);

u = ustar;

     for i = 1:Nxu-1,     % column, 1 to 21

         for j = 1:Nyu-1, % row, 1 to 10

         k = Nyp-((Nyp-j))+1; % 2 to 11

         kv = Nyv-((Nyv-j))+2; % 3 to 12

         ku = Nyu-((Nyu-j))+1; % 2 to 11

            ubar1 = 0.5*(u(ku,i+1)+u(ku-1,i+1));

            ubar2 = 0.5*(u(ku,i)+u(ku-1,i));

            B(kv-1,i+1) = -(((rho*v(kv-1,i+2)*ubar1)-(rho*v(kv-1,i)*ubar2))/(2*dx)...

               + (rho*(v(kv-2,i+1)^2)-rho*(v(kv,i+1)^2))/(2*dy)) ...

               + mu*((v(kv-1,i+2)-2*v(kv-1,i+1)+v(kv-1,i))/(dx^2) ...

               + (v(kv-2,i+1)-2*v(kv-1,i+1)+v(kv,i+1))/(dy^2));

            rho_vstar(kv-1,i+1) = rho_vstar(kv-1,i+1) + B(kv-1,i+1)*dt...

                           - (dt/dy)*(pstar(k-1,i) - pstar(k,i));

            vstar(kv-1,i+1) = rho_vstar(kv-1,i+1)/rho;

            v(kv-1,i+1) = vstar(kv-1,i+1);

        end

     end

rho_vstar(:,Nxv) = rho_vstar(:,Nxv-1); % v velocity zeroth order interpolation

rho_vstar(:,1) = 0; % vertical velocity at inflow boundary

vstar = rho_vstar./rho;

for n = 1:200,

      for i = 2:Nxp-1, % 2 to 20

        for j = 2:Nyp-1, % 2 to 10

        k = Nyp-((Nyp-j)+1)+1; % 2 to 10

        kv = Nyv-((Nyv-j)+1)+2; % 3 to 11

        ku = Nyu-((Nyu-j)+1)+1; % 2 to 10  

    d(k,i) = (1/dx)*(rho_ustar(ku,i+1)-rho_ustar(ku,i))+(1/dy)*(rho_vstar(kv-1,i+1)-rho_vstar(kv,i+1));

    pprime(:,1) = 0; % inflow

    pprime(:,Nxp) = 0; % outflow

    pprime(1,:) = 0; % top wall

    pprime(Nyp,:) = 0; % bottom wall

    pprime1(k,i) = -(1/a)*(b*pprime(k,i+1)+b*pprime(k,i-1)+c*pprime(k-1,i)+c*pprime(k+1,i)+d(k,i));

    conv(k,i) = pprime1(k,i) - pprime(k,i);

    pprime(k,i) = pprime1(k,i);

        end

      end

end

fprintf('\nCompleted iteration %3d\n',iter);

    pstar = pstar + alpha_p.*pprime;

    u = ustar; %+ alpha_u.*uprime;

    v = vstar; %+ alpha_v.*vprime;

    plot(iter,abs(d(7,15))*1e3,'-o'); hold on; grid on; xlabel('Iteration');

    ylabel('d - (slug/ft^3.s) x 10^3');

end

u_new(:,count) = flipud(u(:,15));

yu_new(:,count)= yu(:,15);

end

%% Velocity plot

figure();

plot((u_new(:,1)), yu_new(:,1),'--or','Markerface','r');

hold on;

plot((u_new(:,2)), yu_new(:,2),'--ob','Markerface','b');

plot((u_new(:,3)), yu_new(:,3),'--ok','Markerface','k');

plot((u_new(:,4)), yu_new(:,4),'--om','Markerface','m');

plot((u_new(:,5)), yu_new(:,5),'--og','Markerface','g');

set(gca,'fontsize',12,'fontweight','bold')

legend('k = 4','k = 20','k = 50','k = 150','k = 300','location','southeast')

xlabel('u, ft/s');

ylabel('y, ft');

axis([0 1 0 H]);

grid on;