myprojects-hg-reborn

Fork

rm(list=ls())

# the formulas are taken from the paper by Housiadas & Yannis

# Calculation of the dimensionless temperature profile for a turbulent flow.

#Eq. 14 of the paper

f<-function(Re)
{
factor<-(1/(((8/Re)^10+(Re/36500)^20))^0.5+(2.21*log(Re/7))^10)^(1/5)
myf<-2/factor
return(myf)
}

#eq. 13 in the paper

Nusselt<-function(Re,Prandtl)
{
(f(Re)/2)*(Re-1000)*Prandtl/(1+12.7*(f(Re)/2)^0.5*(Prandtl^(2/3)-1))
}

Nu_Hou<-function()
{
Nusselt(Re,Prandtl)
}


# now I calculate the Nusselt number using another correlation

Nu_dittus_boelter<-function(Re,Prandtl)
{
0.023*Re^0.8*Prandtl^0.4
}

Nu_db<-function()
{
Nu_dittus_boelter(Re,Prandtl)
}


# turbulent dimensionless temperature profile; see page 307, top right

theta_m<-function(Re,Prandtl,x,R)
{
exp(-2*Nusselt(Re,Prandtl)*x/(R*Re*Prandtl))

}

theta_wrap<-function()
{
theta_m(Re,Prandtl,x,R)
}


theta_m_dittus_boelter<-function(Re,Prandtl,x,R)
{
exp(-2*Nu_dittus_boelter(Re,Prandtl)*x/(R*Re*Prandtl))

}


theta_db_wrap<-function()
{
theta_m_dittus_boelter(Re,Prandtl,x,R)
}



Prandtl<-0.679        #0.69

R<-0.05
x_fin<-40
x<-seq(0,x_fin,length=1001)

K<-1

vel_cal<-0  # only if !=0 the code performs also the velocity calculations

# Now I calculate the Reynolds #
u_av<-2.12
D_pipe<-2*R
rho<-0.93
mu<-2.41e-5
Re<-rho*u_av*D_pipe/mu



Nu_housiadas<-Nu_Hou()
print("The Nusselt estimated as in the paper by Housiadas is")
print(Nu_housiadas)


Nu_dittus<-Nu_db()
print("The Nusselt estimated as in the paper by Dittus-Boelter correlation is")
print(Nu_dittus)

print("Re is")
print(Re)


T_prof<-theta_wrap()

T_prof_db<-theta_db_wrap()

pdf("theta_m_comparison.pdf")
plot(x/R,T_prof,"l",lwd=1,col="blue",xlab=expression("distance from the inlet/tube radius"),
ylab=expression("Dimensionless T"),ylim=range(c(min(T_prof,T_prof_db),max(T_prof,T_prof_db))))
lines(x/R,T_prof_db,col="red",lwd=1)
legend(50, 0.85, c("Correlation in paper by Christos", "Dittus-Boelter"), 
lwd= c(1,1),col=c("blue","red" ))
dev.off()


#Now I revert to dimensional units to check how the measured T compares to the theoretical
# prediction

T_0=380.
T_w=303.

T_dim=T_w+(T_0-T_w)*T_prof
T_dim_db=T_w+(T_0-T_w)*T_prof_db

pdf("dimensional_T_comparison.pdf")
plot(x/R,T_dim,"l",lwd=1,col="blue",xlab=expression("distance from the inlet/tube radius"),
ylab=expression("T"),ylim=range(c(min(T_dim,T_dim_db),max(T_dim,T_dim_db))))
lines(x/R,T_dim_db,col="red",lwd=1)
legend(50, 330, c("Correlation in paper by Christos", "Dittus-Boelter"), 
lwd= c(1,1),col=c("blue","red" ))
dev.off()


pdf("dimensional_T_comparison_celsius.pdf")
plot(x/R,T_dim-273.,"l",lwd=1,col="blue",xlab=expression("distance from the inlet/tube radius"),
ylab=expression("T"),ylim=range(c(min(T_dim,T_dim_db)-273.,max(T_dim,T_dim_db)-273.)))
lines(x/R,T_dim_db-273.,col="red",lwd=1)
legend(50, 60, c("Correlation in paper by Christos", "Dittus-Boelter"), 
lwd= c(1,1),col=c("blue","red" ))
dev.off()


theta_star=T_w/(T_0-T_w) #I define the crucial dimensionless constant to be used in the 
# following


# Now I try estimating the deposition coefficients
f_thermo<-function(Re,Prandtl,x,R,K)
{
((theta_star+exp(-2*Nusselt(Re,Prandtl)*x/(R*Re*Prandtl)))/
(1+theta_star))^(Prandtl*K)
}

f_th<-function()
{
f_thermo(Re,Prandtl,x,R,K)
}

thermo<-f_th()


f_th_asym<-(theta_star/(1+theta_star))^(Prandtl*K)





th_penetr<-seq(1:length(x))
th_penetr[]<-f_th_asym

pdf("Thermal_penetration.pdf")
plot(x/R,thermo,"l",lwd=1,col="blue",xlab=expression("distance from the inlet/tube radius"),
ylab=expression("Thermal penetration"),ylim=range(c(min(thermo,th_penetr),max(thermo,th_penetr))))
lines(x/R,th_penetr,lwd=1,col="red")
dev.off()

# Now I read the files with the experimental temperature and I carry out the integration


def_integral<-function(integrand,h)
{
n<-length(integrand)
myseq<-seq(1,n-4,4)
integral<-0
for (i in myseq)
{
integral<- integral + h/45.0*(14.0*integrand[i]+64.0*integrand[i+1] +24.0*integrand[i+2]+64.0*integrand[i+3]+14.0*integrand[i+4])

}
return(integral)
}


temp_cal<-0

{
if (temp_cal !=0)
{

lendat<-10
for(i in 1:lendat)
{
 fn <- paste("temp",i-1,".txt",sep="")
 dat <- read.table(fn,header=TRUE,nrow=20000)
 # ... do your stuff on "dat" here ...
 dat_temp<-as.matrix(dat)
{
if (i ==1) {mymat<-dat_temp}

else
mymat<-cbind(mymat,dat_temp)
}
}

print("OK reading the first files")

mymat_save<-mymat # It is handy to save these data into a separate object and work with mymat only in the following


h<-abs(mymat[200,1]-mymat[201,1])


# Now I need to multiply each file times 2*pi*r

r_seq<-seq(0,0.05,length=20000)

for (i in 1:lendat)
{
mymat[ ,2*i]<-mymat[ ,2*i]*r_seq*2*pi
}

x_seq<-seq(0,(lendat-1))

int_value<-seq(1:lendat)


for(i in 1:lendat)
{
myseq2<-mymat[ ,2*i]
int_value[i]<-def_integral(myseq2,h)
}

int_value<-int_value/(R^2*pi) # I finally have the values after the radial integration

# now I build the new dimensionless sequence

x_seq<-x_seq/R

T_w<-343  # wall temperature

num_theta_m<-(int_value-T_w)/(int_value[1]-T_w)


pdf("thermal_penetration_comparison.pdf")
plot(x/R,T_prof,"l",lwd=1,col="blue",xlab=expression("distance from the inlet/tube radius"),
ylab=expression("Dimensionless T"),ylim=range(c(min(T_prof,num_theta_m),max(T_prof,num_theta_m))))
lines(x_seq,num_theta_m,"b",col="red")
dev.off()

}
}

{ if (vel_cal!=0)
{
# now a new part starts in which I calculate the correlations for the velocity
rm(mymat)

# Now I read the velocity profiles

lendat<-10
for(i in 1:lendat)
{
 fn <- paste("vel",i-1,".txt",sep="")
 dat <- read.table(fn,header=TRUE,nrow=20000)
 # ... do your stuff on "dat" here ...
 dat_temp<-as.matrix(dat)
{
if (i ==1) {mymat<-dat_temp}

else
mymat<-cbind(mymat,dat_temp)
}
}


h<-abs(mymat[200,1]-mymat[201,1])

for(i in 1:lendat)
{
myseq2<-mymat[ ,2*i]*r_seq*2*pi
int_value[i]<-def_integral(myseq2,h)
}

int_value<-int_value/(R^2*pi) # I finally have the values after the radial integration

# At this point I have worked out the integral of the velocity on the tube cross-section

turbul_profile_1_n<-function(r_seq,R,Re)
{

minvector<-c(abs(Re-4e3),abs(Re-2.3e4),abs(Re-1.1e5))
mymin<-which(minvector==min(minvector))
{
if (mymin==1)
n<-6
if (mymin==2)
n<-6.6
if (mymin==3)
n<-7
}
print("n is")
print(n)
z<-(1-r_seq/R)^(1/n)

return(z)
}

t_prof1<-function()
{
turbul_profile_1_n(r_seq,R,Re)
}

r_seq<-seq(0,R,length=20000)

one_seven<-t_prof1()


for(i in 1:lendat)
{
 fn <- paste("u_and_1_n_law",(i-1),".pdf",sep="")
 pdf(fn)
plot(mymat[ ,(2*i-1)]/R,mymat[ ,2*i]/max(mymat[ ,2*i]),"l",col="blue",lwd=1,
ylim=range(c(min(mymat[ ,2*i]/mymat[ ,2*i],min(one_seven)),1)),ylab=expression("u/u_max"),xlab=expression(r/R))
lines(r_seq/R,one_seven,col="red",lwd=1)
dev.off()
}
 }

}



print("So far so good")