Search All of the Math Forum:

Views expressed in these public forums are not endorsed by NCTM or The Math Forum.

Topic: analytical solution laminar flow in rectangular microchannel
Replies: 29   Last Post: Mar 21, 2013 4:58 AM

 Messages: [ Previous | Next ]
 Torsten Hennig Posts: 2,419 Registered: 12/6/04
Re: analytical solution laminar flow in rectangular microchannel
Posted: Nov 5, 2012 2:47 AM

> Yes! I am interested. Thank you.
>
> Best regards,
> Anne

Here is the FORTRAN program.
By the comments included, it should be self-explaining.
If you have difficulties to understand and/or
translate part of the code, do not hesitate
to contact me.

Best wishes
Torsten.

program microflow
implicit none
C
integer nxmax, nymax
parameter (nxmax = 100, nymax = 100)
double precision Lx, Ly, Vpunkt
double precision x, deltax, y, deltay
integer n, nmin
double precision m
C
integer i,j,npoints_x,npoints_y
double precision pi,reihenglied,sum,konstante,
& Vpunkt_check
double precision velocity(nxmax,nymax)
C
C Length(Lx) and height(Ly) of the microchannel [m]
C The channel is considered to be centered at (0,0) so that
C velocity output is for -Lx/2 <= x <= Lx/2 and -Ly/2 <= y <= Ly/2
Lx = 5.0D2*1.0D-6
Ly = 5.0D2*1.0D-6

C Volume flow through the microchannel [m^3/s]
Vpunkt = 3.0009375D-8
C
C Number of points in x- and y-direction for which velocity is
C calculated (npoints_x, npoints_y <= 100)
npoints_x = 40
npoints_y = 40
C
C Minimum number of series terms to be considered
nmin = 10
C
C
C USUALLY NO CHANGES ARE REQUIRED BELOW THIS LINE
C OUTPUT (x-coordinate, y-coordinate, velocity in (x,y))
C is written to a file named 'Velocity')
C
C
pi = dacos(-1.0D0)
C
open(12,file = 'Velocity')
C
sum = 0.0D0
n = 0
10 continue
n = n + 1
m = pi/Ly*dble(2*n-1)
reihenglied = dtanh(5.0D-1*Lx*m)/m**5
sum = sum + reihenglied
if (dabs(reihenglied).le.1.0D-18.and.n.gt.nmin) goto 20
goto 10
20 continue
konstante = 4.0D0*Vpunkt/(Lx*Ly**3)/
& (1.0D0/3.0D0-6.4D1/(Lx*Ly**4)*sum)

deltax = Lx/dble(npoints_x)
deltay = Ly/dble(npoints_y)
x = -Lx/2.0D0 - deltax/2.0D0
i = 0
40 continue
i = i + 1
x = x + deltax
y = -Ly/2.0D0 - deltay/2.0D0
j = 0
50 continue
j = j + 1
y = y + deltay
sum = 0.0D0
n = 0
70 continue
n = n + 1
m = pi/Ly*dble(2*n-1)
reihenglied = dcos(m*y)/m**3*(dcosh(m*x)/dcosh(m*5.0D-1*Lx)-1.0D0)
if (mod(n,2).eq.1) reihenglied = -reihenglied
sum = sum + reihenglied
if (dabs(reihenglied).le.1.0D-18.and.n.gt.nmin) goto 80
goto 70
80 continue
velocity(i,j) = konstante*4.0D0/Ly*sum
write(12,'(3e18.9E3)') x,y,velocity(i,j)
if (j.lt.npoints_y) goto 50
if (i.lt.npoints_x) goto 40
C
C Check whether velocity profile satisfies volume flow constraint
C
Vpunkt_check = 0.0D0
Do i = 1, npoints_x
Do j = 1, npoints_y
Vpunkt_check = Vpunkt_check + velocity(i,j)
end do
end do
Vpunkt_check = deltax*deltay*Vpunkt_check
write(6,*) Vpunkt_check, dabs(Vpunkt-Vpunkt_check)
C
close(12)
C
stop
end

Date Subject Author
10/28/12 Anne
10/29/12 Torsten Hennig
10/29/12 Anne
10/29/12 Torsten Hennig
10/29/12 Anne
10/30/12 Torsten Hennig
10/31/12 Torsten Hennig
11/3/12 Anne
11/3/12 Torsten Hennig
11/4/12 Anne
11/5/12 Torsten Hennig
11/16/12 Anne
11/16/12 Torsten Hennig
11/16/12 Anne
11/19/12 Torsten Hennig
11/19/12 Torsten Hennig
12/3/12 Anne
12/4/12 Torsten Hennig
12/4/12 Anne
12/4/12 Torsten Hennig
12/4/12 Anne
12/4/12 Torsten Hennig
12/4/12 Anne
12/4/12 Torsten Hennig
12/4/12 Anne
12/5/12 Torsten Hennig
12/5/12 Anne
12/5/12 Anne
12/6/12 Torsten Hennig
3/21/13 Anne