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
end
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
end
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’)