get_ipython().run_cell_magic(u'capture', u'storage', u"# The above redirects all output of the below commands to the variable 'storage' instead of displaying it.\n# It can be viewed by typing: 'storage()'\n# Setting up worksheet and importing equations from other worksheets\nload('temp/Worksheet_setup.sage')\nload_session('temp/leaf_enbalance_eqs.sobj')\nfun_loadvars(dict_vars) # any variables defined using var2() are saved with their attributes in dict_vars. Here we re-define them based on dict_vars from the other worksheet.")
def fun_SS(vdict1):
'''
Steady-state T_l, R_ll, H_l and E_l under forced conditions.
Parameters are given in a dictionary (vdict) with the following entries:
a_s, a_sh, L_l, P_a, P_wa, R_s, Re_c, T_a, g_sw, v_w
'''
vdict = vdict1.copy()
vdict[nu_a] = eq_nua.rhs().subs(vdict)
vdict[Re] = eq_Re.rhs().subs(vdict)
vdict[Nu] = eq_Nu_forced_all.rhs().subs(vdict)
vdict[k_a] = eq_ka.rhs().subs(vdict)
vdict[h_c] = eq_hc.rhs().subs(vdict)
vdict[D_va] = eq_Dva.rhs().subs(vdict)
vdict[alpha_a] = eq_alphaa.rhs().subs(vdict)
vdict[rho_a] = eq_rhoa.rhs().subs(vdict)
vdict[Le] = eq_Le.rhs().subs(vdict)
vdict[g_bw] = eq_gbw_hc.rhs().subs(vdict)
vdict[R_ll] = eq_Rll.rhs().subs(vdict)
vdict[H_l] = eq_Hl.rhs().subs(vdict)
vdict[g_tw] = eq_gtw.rhs().subs(vdict)
vdict[C_wa] = eq_Cwl.rhs()(P_wl = P_wa, T_l = T_a).subs(vdict)
vdict[P_wl] = eq_Pwl.rhs().subs(vdict)
vdict[C_wl] = eq_Cwl.rhs().subs(vdict)
vdict[E_lmol] = eq_Elmol.rhs().subs(vdict)
vdict[E_l] = eq_El.rhs().subs(eq_Elmol).subs(vdict)
try:
vdict[T_l] = find_root((eq_Rs_enbal - R_s).rhs().subs(vdict), 273, 373)
except:
print 'too many unknowns for finding T_l: ' + str((eq_Rs_enbal - R_s).rhs().subs(vdict).args())
Tlss = vdict[T_l]
for name1 in [C_wl, P_wl, R_ll, H_l, E_l, E_lmol]:
vdict[name1] = vdict[name1].subs(T_l = Tlss)
if n((E_l + H_l + R_ll - R_s).subs(vdict))>1.:
return 'error in energy balance: El + Hl + Rll - R_s = ' + str(n((E_l + H_l + R_ll - R_s).subs(vdict)))
return vdict
vdict = cdict.copy()
vdict[a_s] = 1.0
vdict[g_sw] = 0.01
vdict[T_a] = 273 + 25.5
vdict[T_w] = vdict[T_a]
vdict[P_a] = 101325
rha = 0.5
vdict[P_wa] = rha*eq_Pwl.rhs()(T_l = T_a).subs(vdict)
vdict[L_l] = 0.03
vdict[Re_c] = 3000
vdict[R_s] = 0.
vdict[v_w] = 1
dict_print(fun_SS(vdict))
width = 0.05
height = 0.03
volume = 310/(100^3)
print 'Volume = ' + str((volume*1000).n()) + ' l'
print 'min flow rate for flushing = ' + str((volume*3/100^3/60).n()) + ' m3/s'
print 'min flow rate for flushing = ' + str((volume*3*1000).n()) + ' l/min'
print 'flow rate for 1 m/s direct wind = ' + str(width*height*1) + ' m3/s'
print 'flow rate for 1 m/s direct wind = ' + str(width*height*1*1000*60) + ' l/m'
print 'flow rate for 5 m/s direct wind = ' + str(width*height*5) + ' m3/s'
print 'flow rate for 5 m/s direct wind = ' + str(width*height*5*1000*60) + ' l/m'
width = 0.05
height = 0.03
volume = 400/(100^3)
print 'Volume = ' + str((volume*1000).n()) + ' l'
print 'min flow rate for flushing = ' + str((volume*3/100^3/60).n()) + ' m3/s'
print 'min flow rate for flushing = ' + str((volume*3*1000).n()) + ' l/min'
print 'flow rate for 1 m/s direct wind = ' + str(width*height*1) + ' m3/s'
print 'flow rate for 1 m/s direct wind = ' + str(width*height*1*1000*60) + ' l/m'
var2('c_pi', 'Heat capacity of insulation material', joule/kilogram/kelvin, latexname='c_{pi}')
var2('lambda_i', 'Heat conductivity of insulation material', joule/second/meter/kelvin)
var2('rho_i', 'Density of insulation material', kilogram/meter^3)
var2('L_i', 'Thickness of insulation material', meter)
var2('A_i', 'Conducting area of insulation material', meter^2)
var2('Q_i', 'Heat conduction through insulation material', joule/second)
var2('dT_i', 'Temperature increment of insulation material', kelvin)
assumptions(c_pi)
eq_Qi = Q_i == dT_i*lambda_i*A_i/L_i
units_check(eq_Qi)
eq_Li = solve(eq_Qi, L_i)[0]
units_check(eq_Li)
units_check(c_pi*rho_i*dT_i*L_i)
(c_pi*rho_i*dT_i*L_i).subs(eq_Li)
vdict = {}
vdict[lambda_i] = 0.022
vdict[c_pi] = 1400
vdict[rho_i] = 30
(c_pi*rho_i*dT_i*L_i).subs(eq_Li).subs(vdict)(A_i = 0.3, dT_i = 0.1, Q_i = 0.01)
vdict = {}
vdict[lambda_i] = 0.038
vdict[c_pi] = 1400
vdict[rho_i] = 15
(c_pi*rho_i*dT_i*L_i).subs(eq_Li).subs(vdict)(A_i = 0.3, dT_i = 0.1, Q_i = 0.01)
units_check(lambda_i*A_i*dT_i*L_i)
eq_Li(A_i = 0.3*0.1*2 + 0.3*0.05*2, dT_i = 0.1, Q_i = 0.01).subs(vdict)
vdict = {}
vdict[lambda_i] = 0.033
vdict[c_pi] = 1400
vdict[rho_i] = 30
(c_pi*rho_i*dT_i*L_i).subs(eq_Li).subs(vdict)(A_i = 0.3, dT_i = 0.1, Q_i = 0.01)
var2('alpha_l', 'Leaf albedo, fraction of shortwave radiation reflected by the leaf', watt/watt)
var2('R_d', 'Downwelling global radiation', joule/second/meter^2)
var2('R_la', 'Longwave radiation absorbed by leaf', joule/second/meter^2, latexname='R_{la}')
var2('R_ld', 'Downwards emitted/reflected global radiation from leaf', joule/second/meter^2, latexname='R_{ld}')
var2('R_lu', 'Upwards emitted/reflected global radiation from leaf', joule/second/meter^2, latexname='R_{lu}')
var2('R_u', 'Upwelling global radiation', joule/second/meter^2)
var2('S_a', 'Radiation sensor above leaf reading', joule/second/meter^2)
var2('S_b', 'Radiation sensor below leaf reading', joule/second/meter^2)
var2('S_s', 'Radiation sensor beside leaf reading', joule/second/meter^2)
eq_Rs_Rd = R_s == (1-alpha_l)*R_d
eq_Sa = S_a == R_d - R_lu
eq_Sb = S_b == R_ld - R_u
eq_Ss = S_s == R_d - R_u
eq_assRldRlu = R_ld == R_lu
solve([eq_Sa, eq_Sb.subs(eq_assRldRlu), eq_Ss], R_d, R_lu, R_u)
eq1 = solve(eq_Sa, R_d)[0]
eq2 = solve(eq_Sb, R_ld)[0].subs(eq_assRldRlu)
eq3 = solve(eq_Ss, R_u)[0]
solve(eq1.subs(eq2).subs(eq3), S_s)
eq_assRu0 = R_u == 0
soln = solve([eq_Sa, eq_Sb, eq_Ss, eq_assRu0], R_d, R_lu, R_ld, R_u)
print soln
eq_Rs_Rll = R_s - R_ll == R_d + R_u - R_lu - R_ld
eq_Rbalance = R_s - R_ll == S_a - S_b
solve([eq_Sa, eq_Sb, eq_Ss, R_d + R_u - R_lu - R_ld == S_a - S_b], R_d, R_lu, R_ld, R_u)
var2('B_l', 'Boundary layer thickness', meter)
vdict = cdict.copy()
Ta = 300
vdict[T_a] = Ta
vdict[T_l] = Ta
vdict[L_l] = 0.15
vdict[v_w] = 0.5
vdict[Re_c] = 3000
vdict[a_s] = 1
vdict[P_a] = 101325
vdict[P_wa] = 0
eq_Bl = B_l == 4.91*sqrt(nu_a*L_l/v_w)
print eq_Bl.subs(eq_nua).subs(vdict)
eq_gbw_hc.subs(eq_hc).subs(eq_Nu_forced_all).subs(eq_Re).subs(eq_rhoa).subs(eq_Le).subs(eq_Dva).subs(eq_alphaa, eq_nua, eq_ka).subs(vdict)
vdict[L_l] = 0.03
vdict[v_w] = 8
print eq_Bl.subs(eq_nua).subs(vdict)
eq_gbw_hc.subs(eq_hc).subs(eq_Nu_forced_all).subs(eq_Re).subs(eq_rhoa).subs(eq_Le).subs(eq_Dva).subs(eq_alphaa, eq_nua, eq_ka).subs(vdict)
vdict[v_w] = v_w
P = plot(eq_Bl.rhs().subs(eq_nua).subs(vdict), (v_w, 0.5,5))
P.axes_labels(['Wind speed (m s$^{-1}$)', 'Boundary layer thickness (m)'])
P
600*0.03^2
var2('W_c', 'Chamber width', meter)
var2('L_c', 'Chamber length', meter)
var2('H_c', 'Chamber height', meter)
var2('V_c', 'Chamber volume', meter^3)
var2('n_c', 'molar mass of gas in chamber', mole)
var2('F_in_v', 'Volumetric flow rate into chamber', meter^3/second, latexname='F_{in,v}')
var2('F_in_mola', 'Molar flow rate of dry air into chamber', mole/second, latexname='F_{in,mol,a}')
var2('F_in_molw', 'Molar flow rate of water vapour into chamber', mole/second, latexname='F_{in,mol,w}')
var2('F_out_mola', 'Molar flow rate of dry air out of chamber', mole/second, latexname='F_{out,mol,a}')
var2('F_out_molw', 'Molar flow rate of water vapour out of chamber', mole/second, latexname='F_{out,mol,w}')
var2('F_out_v', 'Volumetric flow rate out of chamber', meter^3/second, latexname='F_{out,v}')
var2('T_d', 'Dew point temperature of incoming air', kelvin)
var2('T_in', 'Temperature of incoming air', kelvin, latexname='T_{in}')
var2('T_out', 'Temperature of outgoing air (= chamber T_a)', kelvin, latexname='T_{out}')
var2('T_room', 'Lab air temperature', kelvin, latexname='T_{room}')
var2('P_w_in', 'Vapour pressure of incoming air', pascal, latexname='P_{w,in}')
var2('P_w_out', 'Vapour pressure of outgoing air', pascal, latexname='P_{w,out}')
var2('R_H_in', 'Relative humidity of incoming air', latexname='R_{H,in}')
var2('L_A', 'Leaf area', meter^2)
eq_V_c = fun_eq(V_c == W_c*L_c*H_c)
eq_F_in_mola = fun_eq(F_in_mola == (P_a - P_w_in)*F_in_v/(R_mol*T_in))
eq_F_in_molw = fun_eq(F_in_molw == (P_w_in)*F_in_v/(R_mol*T_in))
eq_F_out_mola = fun_eq(F_out_mola == (P_a - P_w_out)*F_out_v/(R_mol*T_out))
eq_F_out_molw = fun_eq(F_out_molw == (P_w_out)*F_out_v/(R_mol*T_out))
eq_F_out_v = fun_eq(F_out_v == (F_out_mola + F_out_molw)*R_mol*T_out/P_a)
eq_Foutv_Finv_Tout = eq_F_out_v.subs(F_out_mola = F_in_mola, F_out_molw = F_in_molw + E_lmol*L_A).subs(eq_F_in_mola, eq_F_in_molw).simplify_full()
units_check(eq_Foutv_Finv_Tout)
eq_Foutv_Finv_Tout
eq_Foutmolw_Finmolw_Elmol = F_out_molw == (F_in_molw + E_lmol*L_A)
units_check(eq_Foutmolw_Finmolw_Elmol)
units_check(eq_F_out_v)
var2('M_air', 'Molar mass of air (kg mol-1)', kilogram/mole, value = 0.02897, latexname='M_{air}')
var2('c_pv', 'Specific heat of water vapour at 300 K', joule/kilogram/kelvin, latexname = 'c_{pv}', value = 1864)
var2('Q_in', 'Internal heat sources, such as fan', joule/second, latexname = 'Q_{in}')
eq_chamber_energy_balance = 0 == H_l*L_A + Q_in + T_in*(c_pa*M_air*F_in_mola + c_pv*M_w*F_in_molw) - (T_out*(c_pa*M_air*F_out_mola + c_pv*M_w*F_out_molw)); show(eq_chamber_energy_balance)
eq_Hl_enbal = solve(eq_chamber_energy_balance.subs(F_out_mola == F_in_mola, F_out_molw == F_in_molw + L_A*E_lmol), H_l)[0].expand()
print units_check(eq_Hl_enbal)
soln = solve(eq_chamber_energy_balance.subs(eq_Foutmolw_Finmolw_Elmol, F_out_mola == F_in_mola), T_out)
print soln
eq_Tout_Finmol_Tin = soln[0]
units_check(eq_Tout_Finmol_Tin).simplify_full().convert().simplify_full()
eq_Tout_Finv_Tin = eq_Tout_Finmol_Tin.subs(eq_F_in_mola, eq_F_in_molw).simplify_full()
print eq_Tout_Finv_Tin
show(eq_Tout_Finv_Tin)
units_check(eq_Tout_Finv_Tin).simplify_full().convert()
eq1 = (eq_F_out_molw + eq_F_out_mola).simplify_full(); show(eq1)
eq2 = eq1.subs(F_out_mola == F_in_mola, eq_Foutmolw_Finmolw_Elmol).subs(eq_F_in_mola, eq_F_in_molw); show(eq2)
soln = solve(eq2,F_out_v); print soln
eq_Foutv_Finv = soln[0]; show(eq_Foutv_Finv)
units_check(eq_Foutv_Finv)
eq_Foutv_Finv
eq_Foutv_Finv_Tout
eq_Tout_Finmol_Tin
assume(F_in_v > 0)
assume(F_out_v > 0)
assume(E_lmol >=0)
show(eq_chamber_energy_balance)
soln = solve(eq_chamber_energy_balance.subs(F_out_mola == F_in_mola).subs(eq_Foutmolw_Finmolw_Elmol).subs(eq_F_in_mola, eq_F_in_molw), T_in)
print soln
eq_T_in_ss = soln[0]
show(eq_T_in_ss)
soln = solve(eq_T_in_ss,Q_in)
print soln
eq_Qin_Tin_Tout = soln[0]
soln = solve(eq_T_in_ss,H_l)
print soln
eq_Hl_Tin_Tout = soln[0]
eq_Hl_Tin_Tout.show()
soln = solve(eq_chamber_energy_balance.subs(F_out_mola == F_in_mola), H_l)
print soln
eq_Hl_Tin_Tout_Fmol = soln[0].simplify_full()
eq_Hl_Tin_Tout_Fmol.show()
var2('T0', 'Freezing point in kelvin',kelvin, value = 273.15)
eq_Pwin_Tdew = P_w_in == 611.21*exp((18.678 - (T_d - T0)/234.5)*((T_d - T0)/(257.14-T0+T_d)))
list_Tdew = np.array(srange(-40,50,6))+273.25
list_Pwin = [eq_Pwin_Tdew.rhs()(T_d = dummy).subs(cdict) for dummy in list_Tdew]
print list_Pwin
P = plot(eq_Pwin_Tdew.rhs().subs(cdict), (T_d, 253,303), frame = True, axes = False, legend_label = 'Arden Buck')
P += plot(eq_Pwl.rhs().subs(cdict), (T_l, 253,303), color = 'red', linestyle = '--', legend_label = 'Clausius-Clapeyron')
P += list_plot(zip(list_Tdew,list_Pwin), legend_label = 'Cellcraft')
P.axes_labels(['Dew point (K)', 'Saturation vapour pressure (Pa)'])
P
eq_F_in_molw.show()
eq_Finv_Finmol = F_in_v == (F_in_mola + F_in_molw)*R_mol*T_in/P_a
units_check(eq_Finv_Finmol)
var2('F_in_va_n', 'Volumetric inflow of dry air at 0oC and 101325 Pa', meter^3/second, latexname='F_{in,v,a,n}')
var2('P_r', 'Reference pressure',pascal, value = 101325)
var2('T_r', 'Reference temperature', kelvin)
eq_Finmola_Finva_ref = fun_eq(F_in_mola == F_in_va_n * P_r/(R_mol*T_r))
eq_Finmolw_Finv = F_in_molw == (P_w_in*F_in_v)/(R_mol*T_in)
print units_check(eq_Finmolw_Finv)
eq_Finv_Finmola = F_in_v == F_in_mola*R_mol*T_in/(P_a - P_w_in)
print units_check(eq_Finv_Finmola)
eq_Finmolw_Finmola_Pwa = fun_eq(eq_Finmolw_Finv.subs(eq_Finv_Finmola))
eq_Finv_Finva_ref = fun_eq(eq_Finv_Finmola.subs(eq_Finmola_Finva_ref))
vdict = cdict.copy()
vdict[F_in_va_n] = 10e-3/60
vdict[T_d] = 273.15 + 10
vdict[P_a] = 101325.
vdict[T_r] = 273.15
vdict[P_w_in] = eq_Pwin_Tdew.rhs().subs(vdict)
print vdict[P_w_in]
print vdict[F_in_va_n]
vdict[T_in] = 273.15+0
inflow0 = eq_Finv_Finva_ref.rhs().subs(vdict)
print 'Volumentric flow at 0 oC: ' + str(inflow0) + ' m3/s'
vdict[T_in] = 273.15+25
inflow25 = eq_Finv_Finva_ref.rhs().subs(vdict)
print 'Volumentric flow at 25 oC: ' + str(inflow25) + ' m3/s'
print '25oC/0oC: ' + str(inflow25/inflow0)
vdict[T_in] = 273.15+25
vdict[P_w_in] = 0.
inflow25 = eq_Finv_Finva_ref.rhs().subs(vdict)
print 'Volumentric flow at 25 oC without added vapour: ' + str(inflow25) + ' m3/s'
eq_Foutv_Finv_Tout.show()
eq_F_out_v.show()
eq_F_in_mola.subs(eq_Finv_Finva_ref).show()
eq_Finv_Finva_ref.show()
eq_F_in_molw.show()
eq_F_out_molw.show()
eq_Foutmolw_Finmolw_Elmol.show()
eq1 = eq_F_out_molw.rhs() == eq_Foutmolw_Finmolw_Elmol.rhs().subs(eq_F_in_molw)
eq1.show()
soln = solve(eq1, P_w_out)
eq_Pwout_Elmol = fun_eq(soln[0].subs(eq_Foutv_Finv_Tout.subs(eq_F_in_molw)).simplify_full())
units_check(eq_Pwout_Elmol)
soln[0].subs(eq_Foutv_Finv_Tout.subs(eq_F_in_molw)).simplify_full().show()
soln[0].subs(eq_F_out_v).subs(F_out_mola = F_in_mola, F_out_molw = F_in_molw + E_lmol*L_A).simplify_full().show()
eq_Pwout_Elmol.subs(eq_Finv_Finva_ref).simplify_full().show()
show(soln[0])
show(eq_Foutv_Finv_Tout)
show(eq_F_in_molw)
show(soln[0].subs(eq_Foutv_Finv_Tout.subs(eq_F_in_molw)).simplify())
show(soln[0].subs(eq_Foutv_Finv_Tout.subs(eq_F_in_molw)).expand())
eq_Elmol_El = E_lmol == E_l/(lambda_E*M_w)
print units_check(eq_Elmol_El)
eq_El_Elmol = E_l == E_lmol*lambda_E*M_w
soln = solve(eq_Pwout_Elmol, P_w_in)
eq_Pwin_Elmol = soln[0].simplify_full()
show(eq_Pwin_Elmol)
vdict = cdict.copy()
vdict[F_in_va_n] = 10e-3/60
vdict[T_d] = 273.15 + 10
vdict[P_a] = 101325.
vdict[T_r] = 273.15
vdict[P_w_in] = eq_Pwin_Tdew.rhs().subs(vdict)
print vdict[P_w_in]
print vdict[F_in_va_n]
vdict[T_in] = 273.15+0
inflow0 = eq_Finv_Finva_ref.rhs().subs(vdict)
print 'Volumentric flow at 0 oC: ' + str(inflow0) + ' m3/s'
vdict[T_in] = 273.15+25
inflow25 = eq_Finv_Finva_ref.rhs().subs(vdict)
vdict[F_in_v] = eq_Finv_Finva_ref.rhs().subs(vdict)
print 'Volumentric flow at 25 oC: ' + str(inflow25) + ' m3/s'
print '25oC/0oC: ' + str(inflow25/inflow0)
vdict[T_in] = 273.15+25
vdict[P_w_in] = 0.
inflow25 = eq_Finv_Finva_ref.rhs().subs(vdict)
print 'Volumentric flow at 25 oC without added vapour: ' + str(inflow25) + ' m3/s'
vdict[E_l] = 100.
vdict[L_A] = 0.03^2
vdict[E_lmol] = eq_Elmol_El.rhs().subs(vdict)
vdict[T_out] = 273+20.
eq_Pwout_Elmol.subs(vdict)
var2('L_s', 'Width of net radiation sensor', meter)
var2('L_ls', 'Distance between leaf and net radiation sensor', meter)
var2('F_s', 'Fraction of radiation emitted by leaf, absorbed by sensor', 1)
Wi = L_s/L_ls
Wj = L_l/L_ls
eq_Fs = F_s == (sqrt((Wi + Wj)^2 + 4) - sqrt((Wj - Wi)^2 + 4))/(2*Wi)
units_check(eq_Fs)
vdict = cdict.copy()
vdict[L_l] = 0.03
vdict[L_s] = 0.01
print eq_Fs.rhs().subs(vdict)(L_ls = 0.01)
P = plot(eq_Fs.rhs().subs(vdict), (L_ls, 0.001, 0.02))
P.axes_labels(['Distance to leaf (m)', 'Fraction of radiation captured'])
P
fun_export_ipynb('leaf_chamber_eqs', 'temp/')
save_session('temp/leaf_chamber_eqs.sobj')
var('m s J Pa K kg mol')
subsdict = {meter: m, second: s, joule: J, pascal: Pa, kelvin: K, kilogram: kg, mole: mol}
var('N_Re_L N_Re_c N_Le N_Nu_L N_Gr_L N_Sh_L')
dict_varnew = {Re: N_Re_L, Re_c: N_Re_c, Le: N_Le, Nu: N_Nu_L, Gr: N_Gr_L, Sh: N_Sh_L}
dict_varold = {v: k for k, v in dict_varnew.iteritems()}
variables = sorted([str(variable.subs(dict_varnew)) for variable in udict.keys()],key=str.lower)
tableheader = [('Variable', 'Description (value)', 'Units')]
tabledata = [('Variable', 'Description (value)', 'Units')]
for variable1 in variables:
variable2 = eval(variable1).subs(dict_varold)
variable = str(variable2)
tabledata.append((eval(variable),docdict[eval(variable)],fun_units_formatted(variable)))
table(tabledata, header_row=True)