Image:Airflow-Obstructed-Duct.png

From Wikipedia, the free encyclopedia

Image
File history
File links

Size of this preview: 800 × 571 pixel
Image in higher resolution (1270 × 907 pixel, file size: 85 KB, MIME type: image/png)

[edit] Summary

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 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 = 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 Viscosity
fluidKinVisc = fluidVisc/fluidDen;

%Problem dimensions
ductLen=5; %m
ductWidth=1; %m

%grid resolution
gridPerLen = 50; % m^(-1)
gridDelta = 1/gridPerLen;
XVec = 0:gridDelta:ductLen-gridDelta;
YVec = 0:gridDelta:ductWidth-gridDelta; 

%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:gridXSize-1
        for j =2:gridYSize-1
                velY(j,i)=0;
                velX(j,i)=inVel;
        end
end


%Set boundary condition on inlet
for i=2:gridYSize-1
        velX(i,1)=inVel;
end

disp(velY(2:gridYSize-1,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:gridXSize-1
                for j=2:gridYSize-1
                        newVelX(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))));
                end
        end
        
        %Copy the data into the front 
        for i=2:gridXSize - 1
                for j = 2:gridYSize-1
                        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,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 mass
        resid=0;
        for i=2:gridXSize-1
                for j=2:gridYSize-1
                        %mass continuity equation using central difference
                        %approx to differential
                        resid = resid + (velX(j,i+ 1)+velY(j+1,i) - ...
                                (velX(j,i-1) + velX(j-1,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)

[edit] Licensing

I, the creator of this work, hereby release it into the public domain. This applies worldwide.
In case this is not legally possible,
I grant any entity the right to use this work for any purpose, without any conditions, unless such conditions are required by law.

File history

Legend: (cur) = this is the current file, (del) = delete this old version, (rev) = revert to this old version.
Click on date to download the file or see the image uploaded on that date.

(del) (cur) 05:45, 24 February 2007 . . User A1 (Talk | contribs) . . 1270×907 (86,796 bytes) (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 paralell with the duct walls. The observed spike is mainly due to numerical limitatio)