SigmaPlot 4.01 Program for calculating the thermal conductivity of the bed of SLH 3-mm
alumina  spheroids, in static helium gas, at a pressure of 99 kPa over a range of temperatures
given in column 13 of the SigmaPlot spreadsheet.  Variables not defined below are as given in
paper submitted to J. Heat and Mass Transfer, 2001, March 2001, by A.J.Slavin, V. Arcas, et al.

By Alan Slavin, Trent University, 2001

T = col(13)
P=99

N=If(P=99,8,0)   'Puts the results in columns 80 to 89
N1=N*10
'Fixed constants
pi=3.14
k=1.38e-23
sig=5.67e-8
L=sqrt(8/3)*R       'Height of unit cell

Pf=0.58 '3mm SLH         'Packing fraction

alphavp=(0.74/Pf)-1
Acp=2*sqrt(3)*R^2        'Area of the close-packed fraction
Av = Acp*alphavp         'Area of the void fraction

Nc =1.5 '(everything else)    'No. of effective contact points per unit cell
thetamax = 60            'Maximum angle of integration 
alphar = 1
alphao=1
alphai=1
alphas = 2/1.5           'alphas/Nc

hr = 1.22e-6    'SLH 3mm 'Average surface roughness of the spheroids


'Solid specifications

Ks = 24.4-0.0341*T+1.38e-5*T^2  'Thermal conductivity of bulk alumina, from supplier

eps=0.75                 'Thermal emissivity

R=1.35/1000              'Average spheroid radius

a=1                 'Accommodation coefficient

'Gas specifications
Cv=1.5*k                                'He
                                                       
gamma=5/3                'Cp/Cv for He

n=P*1000/(k*T)           'Molecular number density
D=2.18e-10                                        'He molecular diameter


lambda=1/(sqrt(2)*n*pi*D*D)   'Molecule mean free path
v=sqrt(8*k*T/(pi*6.64e-27))              'He molecule average velocity

j=1.8*lambda                                       'jump distance in He 


'Different conductivities
Kg = 25*pi*Cv*v*lambda*n/64            'He gas at high pressure

col(N1+8) = if(lambda*2/3 < 2*hr, 0,lambda*2/3-2*hr) 
del=col(N1+8)
   thetalam=sqrt(del/R)  'theta (sub) lambda
   rlambda=sqrt(del*R)   'r (sub) lambda
   Gi=alphai*n*v*Cv*pi*rlambda^2*a*(gamma+1)/(16-(8*a))
col(N1+2)=Gi

Go=alphao*Kg*pi*(1-exp(-R/lambda))*((cos(thetamax)-cos(thetalam))*R 
          + (R+j + hr )*(ln(R+j+hr-R*cos(thetamax)) - ln(R+j+hr-R*cos(thetalam))))     
col(N1+1)=Go

Gs=alphas*Ks*pi*R^2/L
col(N1+3)=Gs

Grp=alphar*Acp*4*sig*T^3/((2/eps)-1)
col(N1+4)=Grp

Grv=alphar*Av*4*sig*T^3/((2/eps)-1)
col(N1+5)=Grv

Ggv=Kg*Av*(1-exp(-R/lambda))/L
col(N1+6) = Ggv

Gcontact=0
col(N1+7)=Gcontact
Gcp=Grp+Nc*(Gs*(Gi+Go+Gcontact))/(Gs+Gi+Go+Gcontact) ' Model 2-series model
Gv=Grv+Ggv
Geff = Gcp+Gv
Keff = (Geff)*(L/(Acp+Av))
col(N1)=Keff             'Overall bed conductivity in W/(m K)