#-------------------------------------------------------------------------------
# Name: Emissivity Calculations
# Purpose: numerical integration of Planck Law
# Version: 1.2
# Python: 2.7
#-------------------------------------------------------------------------------
# Changelog:
# v1.1 -> v1.2
# add emissivity correction for reactor ridges.
# every real emissivity is corrected by e' = e/(1 - F(1-e))
# however note that the report emissivity is still used since it was entered
# into Optris without any correction for ridges and affects report results
# v1 -> v1.1
# Bolometer function added to calculation (instead of assuming constant
# in range 7.5u - 13u)
# Description added to power adjustment approximations
# Results printed for a range of possible inputs:
# band emissivity, bolometer characteristics.
# dummy run "ballpark figure" for temp explained and made weighted av of
# surface temps. NB this is still not done properly!
from math import exp,sin,pi
from functools import partial
import scipy.interpolate as sp
import numpy
import matplotlib.pyplot as plt
#Physical constants (SI units)
c = 2.997e8
k = 1.38e-23
h = 6.626e-34
NUM_PTS = 100 # number of points for numerical integration and interpolation
BANDL = 7.5e-6
BANDH = 14e-6
# Play with these constants to show sensitivity of results to emissivity
BAND_E_BIAS = 1.0 # 1.0 => no change, multiples (1-e) spectral values in band
# e.g. 2 moves values up from 0.9 to 0.95, 0.8 to 0.9 etc
# e.g. 0.5 moves values down 0.9 to 0.8, 0.8 to 0.6 etc
REAL_E_BIAS = 1.0 # 1.0 => no change, multiplies low e values
BOLOMETER_SKEW = '' # '' => normal, 'high' => high lambda bias,
# 'low' => low lambda bias
ROOM_TEMP = 21 # Centigrade
data = {
# datasets from report: Temp/C, Pin, Prad, Pconv, Prods, Pjoule
# each item has as key the number of the relevant active run file or 0 for dummy
12:[1401, 907, 2398, 428, 353, 42 ],
3: [1256, 791, 1725, 385, 308, 36],
0: [400, 486, 183, 133, 130, 7] #NB see note below
}
# Note on dummy figures.
# The data for the dummy run is given in different form from the active runs
# with details of temoeratures from different segments. Rather than complicate
# the analysis we take a surface weighted average temperature here.
# Because convection is much
# more important in the dummy this average is more correct for the calculation
# of the adjustment factor than the max recator body figure.
# It should be noted that the dummy run recalculation is therefore more "hand
# waving than the active runs recalculation. A better job could be done by
# considering each part separately for the dummy run - maybe that would deal
# with the remaining difference between dummy and active runs.
# points defining curve of alumina spectral emissivity in IR band
# NB there is a small change in this curve with temperature which is ignored
alumina_spectral_e_pts = sorted({ 7e-6:0.92, 8e-6:0.94, 10e-6:0.95, 11e-6:0.90,
12e-6:0.52, 13e-6:0.53, 15.5e-6:0.28}.items())
# points defining spectral sensitivity of IR camera within its band
bol_sens_pts = sorted({7e-6:0.1, 7.5e-6:0.2, 8.2e-6:0.9,
8.7e-6:1.0,9.6e-6:0.92, 10.1e-6:0.96, 11e-6:0.85,12e-6:0.5,
13e-6:0.41 ,13.5e-6:0.31,14.2e-6:0.34 }.items())
def bolometer_sensitivity_pts():
return [(l, 0.0 if (l < 10e-6 and BOLOMETER_SKEW=='high') or
(l > 10e-6 and BOLOMETER_SKEW=='low')
else x) for (l,x) in bol_sens_pts]
# points defining the total emissivity curve for alumina used in the report
# note this does not need to be correct, just what the report testers used
# NB temperatures are in C not K
rep_alumina_tot_e_pts = sorted({ 20:0.66, 200:0.7,350:0.79, 480:0.7, 900:0.48,
1250:0.41, 1510:0.4 }.items())
#
# depricated !
# view_factor = 1.0 - sin(pi/4.0) # correct view factor for ridges
#
# New view factor can be found below
view_factor = 0.42816
#------------------------------------------------------------------------------
# UTILITY FUNCTIONS
# -----------------------------------------------------------------------------
def e_adjust(e):
# adjust e by view factor
return e / (1.0 - view_factor*(1.0-e))
def temp_to_k(tempc):
return tempc+273
def frange(x, y, jump):
while x < y:
yield x
x += jump
def weighted_av(a, b, yfun, wfun):
# form weighted av of yfun weighted by wfun over range a to b
ysum = 0
wsum = 0
for y in frange(a, b, (b-a)/NUM_PTS):
w = wfun(y)
wsum += w
ysum += yfun(y)*w
return float(ysum)/wsum
def interp(t1, pts):
if t1 < pts[0][0] or t1 > pts[-1][0]:
print 'out of range interpolation input',t1, pts
raise
for i, (temp, e) in enumerate(pts):
if temp > t1:
temp0,e0 = pts[i-1]
frac = (t1-temp0)/float(temp-temp0)
return frac*e+(1-frac)*e0
# subfunction used to generate point data
def fun_to_pts(a, b, fun):
res = {}
for x in frange(a,b,(b-a)/float(NUM_PTS)):
res[x]=fun(x)
return sorted(res.items())
#------------------------------------------------------------------------------
# FUNCTIONS CALCULATED FROM DATA
# -----------------------------------------------------------------------------
def rep_alumina_tot_e_fun(tp):
# not adjusted for ridges since Report did not do this
return interp(tp,rep_alumina_tot_e_pts)
def real_alumina_tot_e_fun(tp):
# not adjusted for ridges since Report did not do this
e = rep_alumina_tot_e_fun(tp)
if e < 0.65:
e *= REAL_E_BIAS
if e > 0.65: e = 0.65
if e < 0: e = 0
return e
def alumina_spectral_e_fun(l):
# used only for alumina emissivity in bolometer bandwidth
# adjusted for ridges
return e_adjust(1-((1- interp(l, alumina_spectral_e_pts))*BAND_E_BIAS))
def band_sens_weighting_fun(tp, l):
return interp(l, bolometer_sensitivity_pts())*planck(l,tp)
# NB this is a function used to generate points on curve.
# These are then interpolated.
def alumina_band_e_fun_temp(tp):
return weighted_av(BANDL, BANDH, alumina_spectral_e_fun,
partial(band_sens_weighting_fun, tp))
# Interpolates points generated from above function - for efficiency
def alumina_band_e_fun(tp):
return interp(tp, alumina_band_e_pts)
# the ideal black body spectrum function from theory
# constant multiplier here is arbitrary
# tp is temperature in C
# l is wavelength in metres
# result is divided by 1e30 to keep numbers within double bounds.
def planck(l, tp):
return 1/((1e6*l)**5*(exp(c*h/(l*k*temp_to_k(tp)))-1.0))
def planck_band_ratio(tp1, tp2):
vals1 = [band_sens_weighting_fun(tp1,l)
for l in frange(BANDL,BANDH,(BANDH-BANDL)/float(NUM_PTS))]
vals2 = [band_sens_weighting_fun(tp2,l)
for l in frange(BANDL,BANDH,(BANDH-BANDL)/float(NUM_PTS))]
return sum(vals1)/sum(vals2)
def solve(tp):
# iterate to get temperature solving eqn
r_target = rep_alumina_tot_e_fun(tp)/alumina_band_e_fun(tp)
t1 = tp
while planck_band_ratio(t1,tp) > r_target:
r_target = rep_alumina_tot_e_fun(tp)/alumina_band_e_fun(t1)
t1 *= 0.95
while planck_band_ratio(t1,tp) < r_target:
t1 *= 1.01
while planck_band_ratio(t1,tp) > r_target:
t1 *= 0.999
return t1
#------------------------------------------------------------------------------
# DISPLAY FUNCTIONS (REQUIRES NUMPY, SCIPY, PYLAB)
# -----------------------------------------------------------------------------
def interp1(pts):
xv = numpy.asarray([x for (x,y) in pts])
yv = numpy.asarray([y for (x,y) in pts])
xnew = numpy.linspace(pts[0][0], pts[-1][0], 100)
tck = sp.splrep(xv, yv, s=0,k=1)
ynew = sp.splev(xnew, tck,der=0)
plt.plot(xnew, ynew)
plt.title('Band Emissivity')
plt.show()
#------------------------------------------------------------------------------
# MAIN CODE
# -----------------------------------------------------------------------------
def pconv(dats):
[t1,pin, prad,pconv,prods,pjoule] = dats
global BOLOMETER_SKEW, BAND_E_BIAS
global body_area, fin_area
body_area = 12.5e-3
fin_area = 26.3e-3
# print the results for normal, biassed low lambda,
# biassed high lambda, bolometer sensitivity (bskew)
# and for band emissivity (band_bias) shifted higher (0.5) AND LOWER (1.5)
# from nominal
for band_bias in [1.0, 1.5, 0.5]:
for bskew in ['', 'low', 'high']:
BOLOMETER_SKEW = bskew
BAND_E_BIAS = band_bias
if not bskew: bskew = 'none'
t2 = solve(t1)
# Note that solve(t1) finds the real temperature t2 for a given
# reported temperature t1. We then calculate power adjustment
# factors (arad,aconv) = power_adj(t1,t2) for the ratios of real
# rad or conv powers to reported powers
# We apply these factors to the report powers to find the real powers
#
#
#
# The problem is that not all the reactor is at the same temperature
# Available data for power is given as radiation, convection, rods.
# The joule heating is a very small correction to pin.
# The reactor body is at the headline temperature, the caps are at a
# lower temperature. The surface areas are:
#
# depricated
# body = 10*1.25*e-3 = 12.5e-3 m^2
# caps = 6*1.67e-3 = 10e-3 m^2
#
# now assigned as variables to be able to correct for fin area
# see above
# we adjust radiation and convection using the factor for the reactor
# body which makes 55% of the surface area, nearly all of the
# radiated power (because of the high temp dependence) and more than
# half of the convected power. The caps power thus has the wrong
# conversion factor applied: at a lower temperature the reduction
# factor will be less. However for active tests this is a small
# effect, because Prad dominates and the contribution here from caps
# is small.
#
# Note the way the dummy headline temperature is given as a ballpark
# figure between body and caps, weighted by area.
#
# For dummy test the caps contribution is larger than for Prad but
# here the difference between the power_adj factors for the two
# areas is small. Overall the error due to this approximation is
# similarly small. A more accurate calculation could be done for
# the dummy test where more data is available. However this is
# very close for the active test powers because radiation depends
# as T^4 and is nearly all from the reactor body.
# See elsewhere for a fuller error analysis.
#
# The rods temperature varies. We average the rad and conv factors
# for a rough estimate. Given the rods contribution is small
# this will suffice. We use the adjustment factor from this average
# temperature for all the rods power.
pout1 = prad+pconv+prods
ar,ac = power_adj(t1,t2)
ar2,ac2 = power_adj(t1/2, solve(t1/2))
pout2 = prad*ar+pconv*ac+prods*(ar2+ac2)/2
print(
"Bias=%4s,%3.1f: %dC, %dC(real) Pin=%.1f, Pout=%.0f, Pout(real)=%.0f,"\
"COP(real)=%.3f" % (bskew, band_bias, int(t1), int(t2), pin-pjoule, pout1,
round(pout2), pout2/float(pin-pjoule)))
# returns pair of adjustment factors for radiation and convection
# for given report (t1) and real (t2) temps
# NB e is adjusted for ridge effect
def power_adj(t1, t2):
# NB - this is not the right adjustment! it over-adjusts for ridges.
# in reality e_adjust should be applied to spectral emissivity
# of alumina and the results integrated weighted by the Plank Law
# to get adjusted total emissivity.
# still this is an OK approximation, better than no adjustment.
prad_amb = e_adjust(real_alumina_tot_e_fun(ROOM_TEMP))*temp_to_k(ROOM_TEMP)**4
a_rad = (e_adjust(real_alumina_tot_e_fun(t2))*temp_to_k(t2)**4 - prad_amb) / \
(rep_alumina_tot_e_fun(t1)*temp_to_k(t1)**4-prad_amb)
#
# add the fin area correction to a_rad
#
a_rad = ((1 - view_factor)*fin_area/body_area)*a_rad
a_conv = ((t2-ROOM_TEMP)/(t1-ROOM_TEMP))**1.25
return (a_rad,a_conv)
def main():
# This is a calculated set of points on the curve for the alumina effective
# band emissivity with temperature.
# This varies with temperature slightly due to different weighting by bolometer
# sensitivity and black body spectrum
global alumina_band_e_pts
alumina_band_e_pts = fun_to_pts(100, 1600, alumina_band_e_fun_temp)
# print the conversion data and compare with report data
for (n,dats) in data.items():
print n
pconv(dats)
print('e-adjust', 0.4, e_adjust(0.4), e_adjust(0.4)/0.4)
interp1(alumina_band_e_pts)
return
if __name__ == '__main__':
main()