|
|
Re: analytical solution laminar flow in rectangular microchannel
Posted:
Nov 19, 2012 4:00 AM
|
|
> This code is written by a senior. The velocity / flow
> rate is set as 0.15u l/min^-1. The value is then
> convert to form of (m/s).
>
> Regards,
> Anne
I don't know whether the formula you used to calculate
the axial velocity in your code is correct since
I could not find it in the literature.
The formula I used can be found in
http://www.amazon.de/Str%C3%B6mungslehre-Einf%C3%BChrung-Theorie-Str%C3%B6mungen-Springer-Lehrbuch/dp/3642131425
Chapter 6.1.6: Strömung durch nichtkreisförmige Rohre
(Flow through non-circular tubes)
I enclose my attempt to translate the FORTRAN-code
to MATLAB.
Best wishes
Torsten.
%Length(Lx) and height(Ly) of the microchannel [m]
%The channel is considered to be centered at (0,0) so that
%velocity output is for -Lx/2 <= x <= Lx/2 and -Ly/2 <= y <= Ly/2
Lx = 500e-6;
Ly = 500e-6;
%Volume flow through the microchannel [m^3/s]
Vpunkt = 3.0009375e-8;
%Number of points in x- and y-direction for which velocity is
%calculated (npoints_x, npoints_y <= 100)
npoints_x = 40;
npoints_y = 40;
%Minimum number of series terms to be considered
nmin = 10;
%USUALLY NO CHANGES ARE REQUIRED BELOW THIS LINE
%OUTPUT (x-coordinate, y-coordinate, velocity in (x,y))
%is written to a file named 'Velocity')
x_vector = zeros(npoints_x);
y_vector = zeros(npoints_y);
velocity = zeros(npoints_y,npoints_x);
sum = 0;
n = 0;
reihenglied = 1.0;
while (abs(reihenglied) > 1e-18) || (n <= nmin)
n = n + 1;
m = pi/Ly*(2*n-1);
reihenglied = tanh(0.5*Lx*m)/m^5;
sum = sum + reihenglied;
end
konstante = 4.0D0*Vpunkt/(Lx*Ly^3)/(1.0/3.0-64.0/(Lx*Ly^4)*sum)
deltax = Lx/npoints_x;
deltay = Ly/npoints_y;
x = -Lx/2.0 - deltax/2.0;
i = 0;
While (i < npoints_x)
i = i + 1;
x = x + deltax;
x_vector(i) = x;
y = -Ly/2.0 - deltay/2.0;
j = 0;
While (j < npoints_y)
j = j + 1;
y = y + deltay;
y_vector(j) = y;
sum = 0.0;
n = 0 ;
reihenglied = 1.0;
While (abs(reihenglied) > 1e-18) || (n <= nmin)
n = n + 1;
m = pi/Ly*(2*n-1);
reihenglied = cos(m*y)/m^3*(cosh(m*x)/cosh(0.5*m*Lx)-1.0);
if (mod(n,2) == 1)
reihenglied = -reihenglied;
end
sum = sum + reihenglied;
end
velocity(j,i) = konstante*4.0/Ly*sum;
end
end
surf(x_vector,y_vector,velocity);
%
%Check whether velocity profile satisfies volume flow constraint
%
Vpunkt_check = 0.0;
For i = 1:npoints_x
For j = 1:npoints_y
Vpunkt_check = Vpunkt_check + velocity(i,j);
end
end
Vpunkt_check = deltax*deltay*Vpunkt_check;
Vpunkt;
Vpunkt_check;
|
|
