myprojects-hg-reborn

Fork

#! /usr/bin/env python
from scipy import *
import pylab


def v_th(T):
	k_B=1.38e-23
	m_air=28.8*1.66e-27
	v=sqrt((2.*k_B*T)/m_air)
	return v

def rho_T(T):
	rho=p*0.0289644/(8.314472*T)
	return rho



def nu_T(T):
	nu=-1.1555e-14*T**3. + 9.5728e-11*T**2. + 3.7604e-08*T - 3.448e-06
	return nu

def mu_T(T):
	mu=nu_T(T)*rho_T(T)
	return mu
	

def l_air(T):
	l=2.*mu_T(T)/(rho_T(T)*v_th(T))
	return l

def knu_T(T):
	lam_air=l_air(T)
	knu=2.*lam_air/diam_seq
	return knu

def slip(T):
	r=(1.+knu*(1.17+0.53*exp(-0.78/knu)))
	return r

def cond_air(T):
	conductivity=1.5207e-11*T**3.-4.8574e-8*T**2.+1.0184e-4*T-3.9333e-4
	return conductivity

	
	
def K_talbot_test(T,knu,Ra):
	lam_air=l_air(T)
	#knu=2.*lam_air/diam_seq
	cond_carbon=1.7
	#Ra=cond_air(T)/cond_carbon
	slip=1.+knu*(1.17+0.53*exp(-0.78/knu))
	C_s=1.17
	C_t=2.18
	C_m=1.14
	K_tb=2.*C_s*slip*(Ra+C_t*knu)/((1.+3.*C_m*knu)*(1.+2.*Ra+2.*C_t*knu))
	
	return K_tb
	
def K_talbot(T,diam_seq):
	lam_air=l_air(T)
	#print 'lam_air is',lam_air
	#print 'diam_seq is', diam_seq
	knu=2.*lam_air/diam_seq
	#print 'knu is', knu
	cond_carbon=1.7
	Ra=cond_air(T)/cond_carbon
	slip=1.+knu*(1.17+0.53*exp(-0.78/knu))
	C_s=1.17
	C_t=2.18
	C_m=1.14
	K_tb=2.*C_s*slip*(Ra+C_t*knu)/((1.+3.*C_m*knu)*(1.+2.*Ra+2.*C_t*knu))
	#print 'K_tb is', K_tb
	return K_tb

def V_thermo(T):
	K=K_talbot(T,diam_seq)
	nu=nu_T(T)
	Re=U*2.*R/nu
	Nusselt=0.023*Re**0.8*Prandtl**0.4
	vel= -K*nu/T*(T_w-T)/(2.*R/Nusselt)
	return vel


def f_th(T,t):
	theta_star=T_w/(T-T_w)
	K=K_talbot(T,diam_seq)
	nu=nu_T(T) #viscosity, NOT the Nusselt number!
	Re=U*2.*R/nu
	Nusselt=0.023*Re**0.8*Prandtl**0.4
	z=((theta_star+exp(-2.*Nusselt*U*t/(R*Re*Prandtl)))/(1+theta_star))**(Prandtl*K)
	return z

def f_th_asym(T,t):
	theta_star=T_w/(T-T_w)
	K=K_talbot(T,diam_seq)
	nu=nu_T(T)
	Re=U*2.*R/nu
	Nusselt=0.023*Re**0.8*Prandtl**0.4
	z=(theta_star/(1+theta_star))**(Prandtl*K)
	return z

	




#l_air_vec=vectorize(l_air)
T_seq=linspace(293.,450.,100)
#T_new=T
#p=zeros(len(T))
p=101325. #1atm in bars

t=linspace(0.,20.,400)

#################################################################################
#In the following calculations I vary the temperature and fix the particle size#######
#################################################################################

l_mean=map(l_air,T_seq)

diam_seq=70.e-9 #I fix the particle size

pylab.plot(T_seq,l_mean,linewidth=2.)
pylab.xlabel('Temperature[K]',fontsize=20.)
pylab.ylabel('air mean free path [m]',fontsize=20.)
#pylab.legend(('prey population','predator population'))
pylab.title('Air Mean Free Path',fontsize=20.)
pylab.grid(True)
pylab.savefig('./plots-air-mean-free-paths/air_free_path')

pylab.hold(False)

k_air=cond_air(T_seq)

pylab.plot(T_seq,k_air,linewidth=2.)
pylab.xlabel('Temperature[K]',fontsize=20.)
pylab.ylabel('air conductivity [W/mK]',fontsize=20.)
#pylab.legend(('prey population','predator population'))
pylab.title('Gas Conductivity',fontsize=20.)
pylab.grid(True)
pylab.savefig('./plots-air-mean-free-paths/air_conductivity')

pylab.hold(False)



mu=mu_T(T_seq)

pylab.plot(T_seq,mu,linewidth=2.)
pylab.xlabel('Temperature[K]',fontsize=20.)
pylab.ylabel('air dynamic viscosity [Kg/(m*s)]',fontsize=20.)
#pylab.legend(('prey population','predator population'))
pylab.title('Air dynamic viscosity',fontsize=20.)
pylab.grid(True)
pylab.savefig('./plots-air-mean-free-paths/air_dynamic_viscosity')

pylab.hold(False)

rho=rho_T(T_seq)


pylab.plot(T_seq,rho,linewidth=2.)
pylab.xlabel('Temperature[K]',fontsize=20.)
pylab.ylabel('air density [Kg/m**3]',fontsize=20.)
#pylab.legend(('prey population','predator population'))
pylab.title('Air density',fontsize=20.)
pylab.grid(True)
pylab.savefig('./plots-air-mean-free-paths/air_density')

pylab.hold(False)

pylab.plot(T_seq,mu/rho,linewidth=2.)
pylab.xlabel('Temperature[K]',fontsize=20.)
pylab.ylabel('air kinematic viscosity [Kg/(m*s)]',fontsize=20.)
#pylab.legend(('prey population','predator population'))
pylab.title('Air kinematic viscosity',fontsize=20.)
pylab.grid(True)
pylab.savefig('./plots-air-mean-free-paths/air_kinematic_viscosity')

pylab.hold(False)

#print 'rho and T are', rho, T




######################################################################
#Now I fix T (i.e. T_0) and consider later on a set of particle sizes#########
####################################################################


#Now I choose a particular value of T and save the properties of air
input_air=pylab.load('input_air')
input_air=input_air*1.0 # to get a floating point array
T_0=input_air[0] # careful about what I choose here to get a sensible estimate of air properties
T_w=273.+70.
print 'T_0 and T_w are', T_0,T_w
data_save=zeros(6)



data_save[0]=T_0
data_save[1]=l_air(T_0)
data_save[2]=nu_T(T_0) 
data_save[3]=mu_T(T_0)
data_save[4]=rho_T(T_0)
# although it is not calculated here at all, but it is simply an experimental datum,
# I add here the mean axial velocity, which is also included into the input file
U=input_air[1]
Prandtl=0.69 # treated a constant
R=0.05 # radius of the anaconda
data_save[5]=U

#Now some calculations of the depositio efficiency: I fixed T_0 and the particle size
#is still fixed

f_time=f_th(T_0,t)

f_no_time=zeros(len(f_time))
f_no_time[:]=f_th_asym(T_0,t)

#print 'f_no_time is', f_no_time

pylab.plot(t,f_time,t,f_no_time,linewidth=2.)
pylab.xlabel('Ut [m]',fontsize=20.)
pylab.ylabel('Thermophoretic Deposition Efficiency',fontsize=20.)
#pylab.legend(('prey population','predator population'))
pylab.title('Thermophoretic Deposition',fontsize=20.)
pylab.grid(True)
pylab.savefig('./plots-air-mean-free-paths/thermophoretic_deposition_efficiency')

#pylab.hold(False)

pylab.save("output_air", data_save)

print 'data_save is', data_save

nu_2=-1.1555e-14*T_0**3. + 9.5728e-11*T_0**2. + 3.7604e-08*T_0 - 3.448e-06

#print 'the two estimates of nu are', data_save[2], nu_2

#Now a calculation of the thermophoretic coefficient
knudsen=logspace(-2.,1.,10)
Ra=1e-3

K_thermo=K_talbot_test(T_0,knudsen,Ra)

#pylab.semilogx(knudsen,K_thermo,linewidth=2.)
#pylab.xlabel('knudsen number',fontsize=20.)
#pylab.ylabel('K_th',fontsize=20.)
#pylab.title('Thermophoretic coefficient',fontsize=20.)
#pylab.grid(True)
#pylab.savefig('./plots-air-mean-free-paths/thermophoretic_coefficient')

#pylab.hold(False)


#diam_seq=linspace(2.,600.,100)

diam_seq=pylab.load("D_mean_seq")  #now diam_seq stands for a list of diameters

diam_seq=diam_seq/1e9 

print 'the shape of diam_seq is', shape(diam_seq)

#print "diam_seq is", diam_seq

K_bis=K_talbot(T_0,diam_seq)
#print 'K_bis is', K_bis
print 'the shape of k_bis is', shape(K_bis)

#pylab.semilogx((diam_seq*1e9),K_bis,linewidth=2.)
pylab.plot((diam_seq*1e9),K_bis,linewidth=2.)
pylab.xlabel('diameter [nm]',fontsize=20.)
pylab.ylabel('K_th',fontsize=20.)
pylab.title('Thermophoretic coefficient at fixed T-T_w',fontsize=20.)
pylab.grid(True)
pylab.savefig('./plots-air-mean-free-paths/thermophoretic_coefficient_vs_particle_diameter')

pylab.hold(False)

V_th=V_thermo(T_0)



pylab.plot((diam_seq*1e9),V_th,linewidth=2.)
pylab.xlabel('diameter [nm]',fontsize=20.)
pylab.ylabel('K_th',fontsize=20.)
pylab.title('Thermophoretic velocity at fixed T-T_w',fontsize=20.)
pylab.grid(True)
pylab.savefig('./plots-air-mean-free-paths/thermophoretic_velocity_vs_particle_diameter')

pylab.hold(False)

pylab.save('v_thermophoresis',V_th)

#Now I calculate the deposition efficiency

print 'So far so good'