Openopt thermal pv

sq_thermal.py

from __future__ import division
from openopt import *
from FuncDesigner import *

# Constants 
kb = 1.3806503e-23      # Boltzmann's Constant (m2 kg s-2 K-1)
elq = 1.60217646e-19    # Elemental charge (C)
h = 6.626068e-34        # Planck's constant (m2 kg s-1)
c = 299792458.0         # Speed of light (m/s)
fs = 6.8e-5             # Solar etendue (fraction of pi)
Ts = 5762.0             # Temperature of the sun (K)
n = 3.4                 # Refractive index (--)
X = 1.0                 # Solar concentration factor (--)
pi = 3.141592653589793  # duh

# Change units of k and h
kbev = kb/elq           # (J --> electonvolts)
hev = h/elq             # (J --> electonvolts)


# Openopt vars
Teh = oovar()         # Temperature of absorber (want to find this)
Ea, Eb = oovars(2)      # Some range of the solar spectrum which we can absorb


# Make these function of commonly used terms in the equations
p = lambda Ex, betax: (Ex**2. * betax**2. + 2.*Ex * betax + 2.) / betax**3                      # particle prefactor function
q = lambda Ex, betax: (Ex**3. * betax**3. + 3.*Ex**2*betax**2. + 6.*Ex*betax + 6.) / betax**4.  # heat prefactor function


# Equation Q1, which are the absorbed solar photon irradiance ( W / m^2 )
betas = (1./(kbev*Ts))
c1 = (2. * X * fs / (hev**3. * c**2.))
Q1 = elq * (c1 * q(Ea, betas) * exp(-Ea*betas) - c1 * q(Eb, betas) * exp(-Eb*betas))

# Equation Q2 which are the returned photon irradiance ( W / m^2 )
beta = 1./(kbev * Teh)
c2 = (2. * pi / (hev**3. * c**2.))
N2 = c2 * p(Ea, beta) * exp(-Ea*beta) - c2 * p(Eb,beta) * exp(-Eb*beta)
Q2 = elq * ( c2 * q(Ea,beta) * exp(-Ea*beta) - c2 * q(Eb,beta) * exp(-Eb*beta) )

# Equation Q3, which are the collected photon irradiance ( W / m^2 )
c3 = (2. * pi / (hev**3 * c**2))
N3 = elq * c3 * p(Ea,beta) * exp(-Ea*beta)
Q3 = elq * c3 * q(Ea,beta) * exp(-Ea*beta)


# # Check the we get sensible numbers
Teh_guess = 500;
Ea_fix = 1.
Eb_fix = 10.
# print Q1({Ea:Ea_fix, Eb:Eb_fix})
# print Q2({Ea:Ea_fix, Eb:Eb_fix, Teh:Teh_guess})
# print Q3({Ea:Ea_fix, Eb:Eb_fix, Teh:Teh_guess})

# Run the solver
startPoint = { Teh:Teh_guess, Ea:1, Eb:1.4 }
fixedVars=(Ea, Eb )
F = [ (Q1 == Q2 + Q3) ]
solver='interalg'
p = SNLE(F, startPoint)
p.constraints = (Teh>300., Teh<Ts )
r = p.solve(solver, maxSolutions = 1, fixedVars=fixedVars, iprint=True, plot=True)