inter3.f
c <html>
c <head><title></title></head>
c <body>
c <pre>
c
program interpolate
implicit none
c
c Demonstation of 2 simple interpolation methods for
c smoothly joining results from 2 disjoint regions
c
c John Mahaffy 2/3/96
c
c Open up two files in your current directory to receive results
c
open(10,file='nuwlin')
open(11,file='nuwcube')
c
c First a linear weighted transition
c
call nuwlin
c
c Next a cubic weighted transition
c
call nuwcube
c
c If you are in Hammond Lab, change this subroutine to plot results
c
call plotit
c
stop
end
c
subroutine nuwlin
implicit none
real Nulam,Nuturb,Pr,Re,w,Nu
real Re1,Re2,rwfac
c
c Demonstration of linearly weighted transition region
c Calculation of Nusselt Number with weighted transition
c in Re1<Re<Re2
c
c
c John Mahaffy 2/3/96
c
c Nu - Nusselt Number
c Nulam - laminar Nusselt Number
c Nuturb - Turbulent Nusselt Number
c Re - Reynolds number
c w - weighting function in transition zone
c <a name="para"><font color="FF0000">
parameter (Pr=1.0,Nulam=4.0,Re1=640,Re2=2000.,rwfac=1./(Re2-Re1))
c </font></a>
c
Re=0.0
c
c Evaluate the Nusselt number over the range 0.<=Re<=3000.
c
do 100 while (Re.le.3000.)
if(Re.lt.Re1) then
Nu = Nulam
else if (Re.gt.Re2) then
Nu=Nuturb(Re,Pr)
else
w = (Re - Re1)*rwfac
Nu = (1.-w)*Nulam + w*Nuturb(Re,Pr)
endif
write(10,*)Re,Nu
Re=Re+10.
100 continue
return
end
subroutine nuwcube
implicit none
real Nuturb,Nulam,Pr,Re,w,SRe,Re1,Re2, Nu
c
c Demonstration of cubic weighting in transition zone
c Calculation of Nusselt Number with weighted transition
c in Re1<Re<Re2
c
c
c
c John Mahaffy 2/3/96
c
c Nu - Nusselt Number
c Nulam - laminar Nusselt Number
c Nuturb - Turbulent Nusselt Number
c Re - Reynolds number
c SRe - Scaled Reynolds number, ranges from 0 to 1 over interpolation range
c w - weighting function in transition zone
c
parameter (Pr=1.0,Nulam=4.0,Re1=640,Re2=2000.)
c
c Evaluate the Nusselt number over the range 0.<=Re<=3000.
c
Re=0.
do 200 while (Re.le.3000.)
if(Re.lt.Re1) then
Nu=Nulam
else if (Re.gt.Re2) then
Nu=Nuturb(Re,Pr)
else
SRe=(Re-Re1)/(Re2-Re1)
w= (3. - 2.*SRe)*SRe**2
Nu= (1.-w)*Nulam + w*Nuturb(Re,Pr)
endif
write(11,*)Re,Nu
Re=Re+10.
200 continue
return
end
function Nuturb(Re,Pr)
c
c Calculate a Turbulent heat transfer coefficient
c obtained from a Dittus-Boelter correlation
c
c John Mahaffy 2/2/96
c
implicit none
real Re,Pr,Nuturb
c
c Nuturb - Turbulent Nusselt number (Dittus-Boelter correlation)
c Re - Reynolds number
c Pr - Prandl number
c
c Nuturb=0.023*Re**.8*Pr**0.4
c
Nuturb=0.023*exp(log(Re)*0.8+log(Pr)*0.4)
return
end
subroutine plotit
c
c This subroutine is for those of you who are getting bored. At this
c stage of the class you will not be held responsible new material here.
logical fexist
c fexist - flag indicating whether or not the gnuin file exists
c First check to see if a file containing commands to the graphic program
c exists in your local directory
c
inquire (file='gp-int3',exist=fexist)
c
c If it doesn't then issue a unix command to copy it from my directory
c the subroutine "system" is not a standard part of Fortran, but is
c provided by some Unix systems to give a crude form of what is called
c a "system call". Any computer will supply an equivalent command, and
c other functions or subroutines that give more direct access to system
c support.
c
if (.not.fexist) then
call system('cp ~jhm/201/gp-int3 .')
print *, ' Graphics input file copied'
endif
c
c Close the data files so they are ready for reading by gnuplot
c
close ( 11 )
close ( 10 )
c
c If you are sitting at a Hammond lab work station remove the "c" in
c column 1 of the following "call system" to plot your results. If you
c are running from NCSA/BYU Telnet on a Mac, you can see plots by
c making the same change, and by following instructions in "gp-int3"
c to change that file for Tektronics emulation in gnuplot.
c
c call system ('gnuplot gp-int3')
c
return
end
c </pre>
c </body>
c </html>
c