Description 
A simulation using the navierstokes differential equations of the aiflow into a duct at 0.003 m/s (laminar flow). The duct has a small obstruction in the centre that is paralell with the duct walls. The observed spike is mainly due to numerical limitations. This script, which i originally wrote for scilab, but ported to matlab (porting is really really easy, mainly convert comments % > // and change the fprintf and input statements) Matlab was used to generate the image.
%Matlab script to solve a laminar flow
%in a duct problem
%Constants
inVel = 0.003; % Inlet Velocity (m/s)
fluidVisc = 1e5; % Fluid's Viscoisity (Pa.s)
fluidDen = 1.3; %Fluid's Density (kg/m^3)
MAX_RESID = 1e5; %uhh. residual units, yeah...
deltaTime = 1.5; %seconds?
%Kinematic Viscosity
fluidKinVisc = fluidVisc/fluidDen;
%Problem dimensions
ductLen=5; %m
ductWidth=1; %m
%grid resolution
gridPerLen = 50; % m^(1)
gridDelta = 1/gridPerLen;
XVec = 0:gridDelta:ductLengridDelta;
YVec = 0:gridDelta:ductWidthgridDelta;
%Solution grid counts
gridXSize = ductLen*gridPerLen;
gridYSize = ductWidth*gridPerLen;
%Lay grid out with Y increasing down rows
%x decreasing down cols
%so subscripting becomes (y,x) (sorry)
velX= zeros(gridYSize,gridXSize);
velY= zeros(gridYSize,gridXSize);
newVelX= zeros(gridYSize,gridXSize);
newVelY= zeros(gridYSize,gridXSize);
%Set initial condition
for i =2:gridXSize1
for j =2:gridYSize1
velY(j,i)=0;
velX(j,i)=inVel;
end
end
%Set boundary condition on inlet
for i=2:gridYSize1
velX(i,1)=inVel;
end
disp(velY(2:gridYSize1,1));
%Arbitrarily set residual to prevent
%early loop termination
resid=1+MAX_RESID;
simTime=0;
while(deltaTime)
count=0;
while(resid > MAX_RESID && count < 1e2)
count = count +1;
for i=2:gridXSize1
for j=2:gridYSize1
newVelX(j,i) = velX(j,i) + deltaTime*( fluidKinVisc / (gridDelta.^2) * ...
(velX(j,i+1) + velX(j+1,i)  4*velX(j,i) + velX(j1,i) + ...
velX(j,i1))  1/(2*gridDelta) *( velX(j,i) *(velX(j,i+1)  ...
velX(j,i1)) + velY(j,i)*( velX(j+1,i)  velX(j,i+1))));
newVelY(j,i) = velY(j,i) + deltaTime*( fluidKinVisc / (gridDelta.^2) * ...
(velY(j,i+1) + velY(j+1,i)  4*velY(j,i) + velY(j1,i) + ...
velY(j,i1))  1/(2*gridDelta) *( velY(j,i) *(velY(j,i+1)  ...
velY(j,i1)) + velY(j,i)*( velY(j+1,i)  velY(j,i+1))));
end
end
%Copy the data into the front
for i=2:gridXSize  1
for j = 2:gridYSize1
velX(j,i) = newVelX(j,i);
velY(j,i) = newVelY(j,i);
end
end
%Set free boundary condition on inlet (dv_x/dx) = dv_y/dx = 0
for i=1:gridYSize
velX(i,gridXSize)=velX(i,gridXSize1);
velY(i,gridXSize)=velY(i,gridXSize1);
end
%y velocity generating vent
for i=floor(2/6*gridXSize):floor(4/6*gridXSize)
velX(floor(gridYSize/2),i) = 0;
velY(floor(gridYSize/2),i1) = 0;
end
%calculate residual for
%conservation of mass
resid=0;
for i=2:gridXSize1
for j=2:gridYSize1
%mass continuity equation using central difference
%approx to differential
resid = resid + (velX(j,i+ 1)+velY(j+1,i)  ...
(velX(j,i1) + velX(j1,i)))^2;
end
end
resid = resid/(4*(gridDelta.^2))*1/(gridXSize*gridYSize);
fprintf('Time %5.3f \t log10Resid : %5.3f\n',simTime,log10(resid));
simTime = simTime + deltaTime;
end
mesh(XVec,YVec,velX)
deltaTime = input('\nnew delta time:');
end
%Plot the results
mesh(XVec,YVec,velX)
