% plot_2D_duct_flow.m

%

% This MATLAB script file plots the

% analytical solution of the velocity

% profile for laminar, pressure-driven

% flow down a rectangular duct.

%

% K. Beers. MIT ChE. 7/29/2002

% First, set the truncation order of

% the expansion.

N = 10;

resolution = 20;

% Next, set the physical parameters.

G = -1;   % pressure gradient

mu = 1;   % viscosity

Lx = 1;   % duct length in x direction

Ly = 1;   % duct length in y direction

% Next, make a 2-D grid.

x_vect = linspace(0,Lx,resolution);

y_vect = linspace(0,Ly,resolution);

[X,Y] = meshgrid(x_vect,y_vect);

% Now, we calculate the values of the

% velocity profile at each grid point.

VZ = 0*X;

% sum over the orders of the expansion

var1 = -16*G/mu/pi^4;

for i=0:N

    m = 2*i+1;

    for j=0:N

        n = 2*j+1;

        factor1 = var1/m/n/(m^2/Lx^2 + n^2/Ly^2);

        % for every point in the grid

        for p=1:size(X,1)

            for q=1:size(X,2)

                x = X(p,q);

                sin_mx = sin(m*pi*X(p,q)/Lx);

                y = Y(p,q);

                sin_ny = sin(n*pi*Y(p,q)/Ly);

                VZ(p,q) = VZ(p,q) + ...

                    factor1*sin_mx*sin_ny;

            end

        end

    end

end

figure;

contourf(X,Y,VZ,20);

colorbar;

title('Laminar pressure-driven duct flow');

xlabel('x');

ylabel('y');