|
|
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
|
|
