提交: 13f9edc99700dc12b2e211e1cc2cb7d0e955f2be - myprojects-hg-reborn (hg)

Commit MetaInfo

• 修订版: 13f9edc99700dc12b2e211e1cc2cb7d0e955f2be (tree)
• 时间: 2008-03-07 02:50:37
• 作者: iselllo
• Commiter: iselllo

Log Message

I added lambda_brownian.py which is a better version of test_kernel.py.
A few things have to be said: I am using a power-law effective density
which does not become a constant for small particles. For fractal
dimension different from 3, I am taking the constant a_0 in the
free-molecular (fractal) kernel as d_0/2, where d_0 is the monomer
diameter. However, it is not clear at all how to extend that to the case
of a constant effective density below a threshold which is NOT the
monomer diameter.

更改概述

• add: Python-codes/lambda_brownian.py

差异

diff -r 77a27c4a75fb -r 13f9edc99700 Python-codes/lambda_brownian.py
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Python-codes/lambda_brownian.py	Thu Mar 06 17:50:37 2008 +0000
@@ -0,0 +1,326 @@
+#! /usr/bin/env python
+# from scipy import *
+# import pylab  # used to read the .csv file
+
+
+import scipy as s
+import numpy as n
+import pylab as p
+
+
+
+
+
+
+def Brow_ker_fuchs_class_optim_slip(Vlist,diam_seq):
+	knu=2.*lam_air/diam_seq # vector with the Knudsen numbers
+	k_B=1.38e-23
+	#print 'knu is', knu
+#         if (save_knu==1):
+#             pylab.save("knudsen",knu)
+	if (slip_yannis==0):
+		#k_B=1.38e-23
+		Diff=k_B*T_0/(3.*s.pi*mu*diam_seq)
+	elif (slip_yannis==1):
+		#k_B=1.38e-23
+		#print 'k_B is ', k_B
+		#print 'T is, ', T_0
+		Diff=k_B*T_0/(3.*s.pi*mu*diam_seq)*(1.+knu*(1.17+0.53*s.exp(-0.78/knu)))
+	#print 'Diff is', Diff
+	m=rho_p*Vlist # this holds in general (i.e. for Vlist !=3.)
+	## as long as Vlist is the volume of solid
+	##and not the space occupied by the agglomerate structure
+	c=(8.*k_B*T_0/(s.pi*m))**0.5
+	#print 'c is', c
+	l=8.*Diff/(s.pi*c)
+        # if (save_knu==1):
+#             pylab.save("lambda_part",l)
+	#print 'l is', l
+	g=1./(3.*diam_seq*l)*((diam_seq+l)**3.-(diam_seq**2.+l**2.)**1.5)-diam_seq
+
+	beta=(\
+	(diam_seq[:,s.newaxis]+diam_seq[s.newaxis,:])   \
+	/(diam_seq[:,s.newaxis]+diam_seq[s.newaxis,:]\
+	+2.*(g[:,s.newaxis]**2.+g[s.newaxis,:]**2.)**0.5)\
+	+ 8.*(Diff[:,s.newaxis]+Diff[s.newaxis,:])/\
+	((c[:,s.newaxis]**2.+c[s.newaxis,:]**2.)**0.5*\
+	(diam_seq[:,s.newaxis]+diam_seq[s.newaxis,:]))\
+	)**(-1.)
+	
+	#print 'the mean of beta is', mean(beta)
+	
+	## now I have all the bits for the kernel matrix
+	kern_mat=Brow_ker_cont_optim_slip(Vlist,diam_seq)*beta
+	if (D_f==3. and slip_yannis==1):
+		p.save("beta_coalescing.dat", s.diag(beta))
+	elif (D_f !=3. and slip_yannis == 1):
+		p.save("beta_fractal.dat", s.diag(beta))
+	#print 'beta is', beta
+
+
+	
+	return kern_mat
+
+
+def Brow_ker_cont_optim_slip(Vlist,diam_seq):
+	k_B=1.38e-23
+	# same expression for the kernel in the continuous limit
+	# as the one used by Maricq.
+	knu=2.*lam_air/diam_seq # vector with the Knudsen numbers
+	if (slip_yannis ==0):
+		Slip=s.zeros(len(knu))
+		Slip[:]=1.
+	elif (slip_yannis==1):
+		Slip=1.+knu*(1.17+0.53*s.exp(-0.78/knu))
+	kern_mat=2.*k_B*T_0/(3.*mu)*(Vlist[:,s.newaxis]**(1./D_f)+\
+	Vlist[s.newaxis,:]**(1./D_f))* \
+	(Slip[:,s.newaxis]/Vlist[:,s.newaxis]**(1./D_f)+\
+	Slip[s.newaxis,:]/Vlist[s.newaxis,:]**(1./D_f))
+	return kern_mat
+
+def Brow_ker_free_optim(Vlist):
+	k_B=1.38e-23
+	a_0=d_0/2.
+	#a_0=1e-8
+	lam=2./D_f-0.5
+	kern_mat=(3./(4.*s.pi))**lam*(6.*k_B*T_0/rho_p)**0.5*a_0**(2.-6./D_f)*\
+	(1./Vlist[:,s.newaxis]+1./Vlist[s.newaxis,:])**0.5*\
+	(Vlist[:,s.newaxis]**(1./D_f)+Vlist[s.newaxis,:]**(1./D_f))**2.
+	return kern_mat
+
+
+
+
+def beta_calc(Vlist,diam_seq):
+	knu=2.*lam_air/diam_seq # vector with the Knudsen numbers
+	#print 'knu is', knu
+	if (slip_yannis==0):
+		Diff=k_B*T_0/(3.*s.pi*mu*diam_seq)*((5.+4.*knu+6.*knu**2.+18.*knu**3.)/\
+		(5.-knu+(8.+s.pi)*knu**2.))
+	elif (slip_yannis==1):
+		#print 'k_B is ', k_B
+		#print 'T is, ', T_0
+		Diff=k_B*T_0/(3.*s.pi*mu*diam_seq)*(1.+knu*(1.17+0.53*exp(-0.78/knu)))
+	#print 'Diff is', Diff
+	m=rho_p*Vlist # this holds in general (i.e. for Vlist !=3.)
+	## as long as Vlist is the volume of solid
+	##and not the space occupied by the agglomerate structure
+	c=(8.*k_B*T_0/(s.pi*m))**0.5
+	#print 'c is', c
+	l=8.*Diff/(s.pi*c)
+	#print 'l is', l
+	g=1./(3.*diam_seq*l)*((diam_seq+l)**3.-(diam_seq**2.+l**2.)**1.5)-diam_seq
+
+	beta=(\
+	(diam_seq[:,s.newaxis]+diam_seq[s.newaxis,:])   \
+	/(diam_seq[:,s.newaxis]+diam_seq[s.newaxis,:]\
+	+2.*(g[:,s.newaxis]**2.+g[s.newaxis,:]**2.)**0.5)\
+	+ 8.*(Diff[:,s.newaxis]+Diff[s.newaxis,:])/\
+	((c[:,s.newaxis]**2.+c[s.newaxis,:]**2.)**0.5*\
+	(diam_seq[:,s.newaxis]+diam_seq[s.newaxis,:]))\
+	)**(-1.)
+	
+	#print 'the mean of beta is', mean(beta)
+	
+	## now I have all the bits for the kernel matrix
+	
+
+	#print 'beta is', beta
+	
+
+	
+	return beta
+
+
+
+
+
+def lambda_Brownian(T_0,diam_seq,V_list,rho_p,lam_air):
+	knu=2.*lam_air/diam_seq
+	k_B=1.38e-23
+	Diff=k_B*T_0/(3.*s.pi*mu*diam_seq)*(1.+knu*(1.17+0.53*s.exp(-0.78/knu)))
+	m=rho_p*V_list
+	c=(8.*k_B*T_0/(s.pi*m))**0.5
+	l=8.*Diff/(s.pi*c)
+
+	return l
+
+
+def v_th(T):
+	k_B=1.38e-23
+	m_air=28.8*1.66e-27
+	v=s.sqrt((2.*k_B*T)/m_air)
+	return v
+
+
+
+T_0=300.
+
+d_0=50.e-9  # primary size in m
+
+
+D_f=3.
+diam_seq=s.logspace(s.log10(1.),s.log10(1000.),1000)
+diam_seq=diam_seq*1e-9 #to go from nm to m
+V_list=s.pi/6.*diam_seq**D_f
+rho_p=1300. # kg/m^3
+mu=184.6e-7
+v_thermal=v_th(T_0)
+rho_air=1.16 #kg/m^3
+lam_air=2.*mu/(rho_air*v_thermal)
+
+print "V_thermal is, ", v_thermal
+
+
+lam_br=lambda_Brownian(T_0,diam_seq,V_list,rho_p,lam_air)
+	
+	
+p.loglog(diam_seq*1e9,lam_br*1e9,linewidth=2.)
+p.xlabel('Particle Diameter [nm]',fontsize=20.)
+p.ylabel('Lambda Brownian [nm]',fontsize=20.)
+#pylab.legend(('prey population','predator population'))
+p.title('Lambda Brownian particle',fontsize=20.)
+p.grid(True)
+p.savefig('lambda_brownian.pdf')
+p.hold(False)
+p.clf()
+
+
+
+slip_yannis=1
+
+
+
+my_kernel=Brow_ker_fuchs_class_optim_slip(V_list,diam_seq)
+
+
+
+
+
+free_kernel=Brow_ker_free_optim(V_list)
+
+cont_kernel=Brow_ker_cont_optim_slip(V_list,diam_seq)
+
+slip_yannis=0
+
+cont_kernel_no_slip=Brow_ker_cont_optim_slip(V_list,diam_seq)
+
+
+
+
+
+
+
+
+p.loglog(diam_seq*1e9,s.diag(my_kernel),diam_seq*1e9,s.diag(free_kernel)\
+	     ,diam_seq*1e9,s.diag(cont_kernel),diam_seq*1e9,s.diag(cont_kernel_no_slip))
+#pylab.axvline(x=select_diam*1e9)
+p.xlabel('diameter [nm]')
+p.ylabel('kernel value along diagonal')
+p.legend(('Fuchs','Free Molecular','Continuum','continuum no slip'))
+p.title('Kernel number vs diameter')
+p.grid(True)
+p.savefig('test_kernel_diagonal_coalescing_particles.pdf')
+p.hold(False)
+p.clf()
+
+
+
+
+
+
+
+
+
+
+
+D_f=2.1 #I now use a non-trivial fractal dimension
+
+
+#V_list=s.pi/6.*diam_seq**3.*(diam_seq/d_0)**(D_f-3.)*(diam_seq>=d_0)+ \
+	#s.pi/6.*diam_seq**3.*(diam_seq<d_0)
+
+V_list=s.pi/6.*diam_seq**3.*(diam_seq/d_0)**(D_f-3.)
+
+
+lam_br=lambda_Brownian(T_0,diam_seq,V_list,rho_p,lam_air)
+	
+	
+p.loglog(diam_seq*1e9,lam_br*1e9,linewidth=2.)
+#p.ylim(4.,120.)
+p.xlabel('Particle Diameter [nm]',fontsize=20.)
+p.ylabel('Lambda Brownian [nm]',fontsize=20.)
+#pylab.legend(('prey population','predator population'))
+p.title('Lambda Fractal Brownian particle',fontsize=20.)
+p.grid(True)
+p.savefig('lambda_brownian_fractal.pdf')
+p.hold(False)
+p.clf()
+
+
+slip_yannis=1
+
+
+
+my_kernel=Brow_ker_fuchs_class_optim_slip(V_list,diam_seq)
+
+
+
+
+
+free_kernel=Brow_ker_free_optim(V_list)
+
+cont_kernel=Brow_ker_cont_optim_slip(V_list,diam_seq)
+
+slip_yannis=0
+
+cont_kernel_no_slip=Brow_ker_cont_optim_slip(V_list,diam_seq)
+
+
+
+
+
+p.loglog(diam_seq*1e9,s.diag(my_kernel),diam_seq*1e9,s.diag(free_kernel)\
+	     ,diam_seq*1e9,s.diag(cont_kernel),diam_seq*1e9,s.diag(cont_kernel_no_slip))
+#pylab.axvline(x=select_diam*1e9)
+p.xlabel('diameter [nm]')
+p.ylabel('kernel value along diagonal')
+p.legend(('Fuchs','Free Molecular','Continuum','continuum no slip'))
+p.title('Kernel number vs diameter')
+p.grid(True)
+p.savefig('test_kernel_diagonal_fractal_particles.pdf')
+p.hold(False)
+p.clf()
+
+beta_coal=p.load("beta_coalescing.dat")
+
+
+p.loglog(diam_seq*1e9,beta_coal, linewidth=2.)
+#pylab.axvline(x=select_diam*1e9)
+p.xlabel('diameter [nm]')
+p.ylabel('beta')
+#p.legend(('Fuchs','Free Molecular','Continuum','continuum no slip'))
+p.title('Beta vs diameter for coalescing Particles')
+p.grid(True)
+p.savefig('beta_coalescing_particles.pdf')
+p.hold(False)
+p.clf()
+
+
+
+beta_frac=p.load("beta_fractal.dat")
+
+
+p.loglog(diam_seq*1e9,beta_frac, linewidth=2.)
+#pylab.axvline(x=select_diam*1e9)
+p.xlabel('diameter [nm]')
+p.ylabel('beta')
+#p.legend(('Fuchs','Free Molecular','Continuum','continuum no slip'))
+p.title('Beta vs diameter for fractal Particles')
+p.grid(True)
+p.savefig('beta_fractal_particles.pdf')
+p.hold(False)
+p.clf()
+
+
+print "So far so good"