A simulation using the navier-stokes 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 parallel 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%ConstantsinVel = 0.003; % Inlet Velocity (m/s)fluidVisc = 1e-5; % Fluid's Viscoisity (Pa.s)fluidDen = 1.3; %Fluid's Density (kg/m^3)MAX_RESID = 1e-5; %uhh. residual units, yeah...deltaTime = 1.5; %seconds?%Kinematic ViscosityfluidKinVisc = fluidVisc/fluidDen;%Problem dimensionsductLen=5; %mductWidth=1; %m%grid resolutiongridPerLen = 50; % m^(-1)gridDelta = 1/gridPerLen;XVec = 0:gridDelta:ductLen-gridDelta;YVec = 0:gridDelta:ductWidth-gridDelta; %Solution grid countsgridXSize = 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 conditionfor i =2:gridXSize-1for j =2:gridYSize-1velY(j,i)=0;velX(j,i)=inVel;endend%Set boundary condition on inletfor i=2:gridYSize-1velX(i,1)=inVel;enddisp(velY(2:gridYSize-1,1));%Arbitrarily set residual to prevent%early loop terminationresid=1+MAX_RESID;simTime=0;while(deltaTime) count=0;while(resid > MAX_RESID && count < 1e2) count = count +1;for i=2:gridXSize-1for j=2:gridYSize-1newVelX(j,i) = velX(j,i) + deltaTime*( fluidKinVisc / (gridDelta.^2) * ...(velX(j,i+1) + velX(j+1,i) - 4*velX(j,i) + velX(j-1,i) + ...velX(j,i-1)) - 1/(2*gridDelta) *( velX(j,i) *(velX(j,i+1) - ...velX(j,i-1)) + 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(j-1,i) + ...velY(j,i-1)) - 1/(2*gridDelta) *( velY(j,i) *(velY(j,i+1) - ...velY(j,i-1)) + velY(j,i)*( velY(j+1,i) - velY(j,i+1))));endend%Copy the data into the front for i=2:gridXSize - 1for j = 2:gridYSize-1velX(j,i) = newVelX(j,i);velY(j,i) = newVelY(j,i);endend%Set free boundary condition on inlet (dv_x/dx) = dv_y/dx = 0for i=1:gridYSizevelX(i,gridXSize)=velX(i,gridXSize-1);velY(i,gridXSize)=velY(i,gridXSize-1);    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),i-1) = 0;    end    %calculate residual for %conservation of massresid=0;for i=2:gridXSize-1for j=2:gridYSize-1%mass continuity equation using central difference%approx to differentialresid = resid + (velX(j,i+ 1)+velY(j+1,i) - ...(velX(j,i-1) + velX(j-1,i)))^2;endendresid = resid/(4*(gridDelta.^2))*1/(gridXSize*gridYSize);fprintf('Time %5.3f \t log10Resid : %5.3f\n',simTime,log10(resid));    simTime = simTime + deltaTime;endmesh(XVec,YVec,velX)deltaTime = input('\nnew delta time:');end%Plot the resultsmesh(XVec,YVec,velX)