Fluid Viscosity

Solution

ibdullahnew.m

clear all;

%% Part 1 and 2

T=[300 400 500 400]; % all temperature values are kept in an array as given in part 2

miu=(2.424e-5)*(10.^(247.8./(T-140))); % given expression for miu

dpdx=[1.5 1.5 1.5 0.5]; % given values of dP/dx

index=length(T); % length of temperature array

Ri=-dpdx./miu; % coefficient

% boundary conditions

y1 = 0; alpha = 0;

y2 = 0.1*1e-3; beta = 0; % converted into meter

npoints = 10; % number of grid points

h = (y2-y1)/npoints; % compute interval width

% We will solve the system of miu*(d2u/dy2)=-dP/dx in central differencing

% method and will give it a form of linear matrix equation like AY=B, where

% A is a matrix of size (npoints-1)*(npoints-1) and B and Y are

% (npoints-1)*1. We will solve it for Y by matrix division and and then add

% the boundary conditions to it.

% Initialization and allocation of vectors

b = zeros(npoints-1,1);

A = zeros(npoints-1, npoints-1);

Y = zeros(npoints-1,1);

% first for loop is for different temperature and the 2nd one for each

% points along y direction

forkk=1:index

for i=1:(npoints-1)

Y(i) = y1 + i*h; % value of Y vector at each point

ri=Ri(kk); % value of coefficient

b(i) = h^2*ri; % value of b at each point

% here we calculate A matrix

for j = 1:(npoints-1)

if j == i % the diagonal

A(i,j) = -2;

elseif j == i-1 % left of the diagonal

A(i,j) = 1;

elseif j == i+1 % right of the diagonal

A(i,j) = 1;

else

A(i,j) = 0; % off the tri-diagonal

end

U=A\b; % matrix division

y(:,kk) = [y1; Y; y2]; % appending boundary values

u(:,kk) = [alpha; U; beta]; % appending boundary values

end

figure

plot(y,u,’Linewidth’,2)

xlabel(‘Distance between plates (m)’)

ylabel(‘fluid velocity (m/s)’)

xlim([0 1e-4])

legendCell = cellstr(num2str(T’, ‘T=%-dK’));

legend(legendCell);

%% part 3: comparison of part 2a result with exact solution

mu=miu(1);

d=0.1*1e-3;

x=y(:,1);

Pdrop=dpdx(1);

uexact = (Pdrop)*d^2/2/mu*(x/d-(x/d).^2); % here we find the exact solution. One can easily do it

figure

plot(x,uexact,’b-‘,x,u(:,1),’r–‘,’Linewidth’,2);

xlabel(‘Distance between plates (m)’)

ylabel(‘fluid velocity (m/s)’)

legend(‘Exact Solution’,’Numerical Solution’)

ibdullahpart4.m

clear all;

%% Part 4 comments are same as of the previous code. This is done with part 2a values

T=300;

miu=(2.424e-5)*(10.^(247.8./(T-140)));

dpdx=1.5;

ri=-dpdx./miu; % coefficient

y1 = 0; alpha = 0;

y2 = 0.1*1e-3; beta = 0; % converted into meter

% convergence test done for 4 sample grids

points =[20 30 40 50];

index=length(points);

forkk=1:index

npoints=points(kk);

% compute interval width

h = (y2-y1)/npoints;

H(kk)=h; % step size of grid

% preallocate and shape the b vector and A-matrix

b = zeros(npoints-1,1);

A = zeros(npoints-1, npoints-1);

Y = zeros(npoints-1,1);

for i=1:(npoints-1)

Y(i) = y1 + i*h; % value of Y vector at each point

b(i) = h^2*ri; % value of b at each point

% here we calculate A matrix

for j = 1:(npoints-1)

if j == i % the diagonal

A(i,j) = -2;

elseif j == i-1 % left of the diagonal

A(i,j) = 1;

elseif j == i+1 % right of the diagonal

A(i,j) = 1;

else

A(i,j) = 0; % off the tri-diagonal

end

U=A\b;

y = [y1; Y; y2];

u = [alpha; U; beta];

d=y2;

x=y;

Pdrop=dpdx;

uexact = (Pdrop)*d^2/2/miu*(x/d-(x/d).^2);

%%% Error is calculated as the sum of the difference between exact and

%%% obtained solution

error(kk)=abs(sum(uexact-u));

end

% order of grid convergence index is defined as the slope of the equation

% log(error)=log(Const)+ ORDER*log(H) where H is the step size of grid

ORDER=log(error)./log(H);

plot(points,ORDER,’r*-‘,’Linewidth’,2);

xlabel(‘Number of grid points’);

ylabel(‘Order of Grid Convergence’)