REM***Morgan correlation; convective heat transfer
REM***coefficient for a cylinder, free convection
CLS : noo# = .00008959#: k# = .06#: pr# = .693: d# = .1#
aconv# = 3.2#: arad# = 2.6#: g# = 9.81#
tinfin# = 614.625#: tc# = 553.15#: tw# = 1123.15#
REM***assume remaining mass is 10 kg per corpse and Cp = 1600 J/kg.K
mcorpse# = 20#: cpcorpse# = 1600
REM***Initialize Bi/Fo surface temperature subroutine
PRINT "Finding xi roots and C coefficients"
bi# = 1.523#: stp# = .00001#: pi# = 3.141592653#: DIM xi#(40): DIM c#(40)
FOR n# = 1# TO 40#: REM***Find 1st 40 roots of transcendental equation
FOR xi# = (n# - 1#) * pi# TO n# * pi# STEP stp#
diff# = bi# - xi# * TAN(xi#)
IF diff# < 0# THEN 20
NEXT xi#
20 xi#(n#) = xi#: REM***Find the corresponding Cn coefficient
c#(n#) = 4# * SIN(xi#) / (2# * xi# + SIN(2# * xi#))
NEXT n#
REM***Start simulations loop
FOR s = 1 TO 1000000: CLS
tfilm# = .5# * (tc# + tinfin#): alpha# = 1# / tfilm#
gr# = g# * alpha# * (tc# - tinfin#) * d# ^ 3# / noo# ^ 2#
gr# = ABS(gr#)
ra# = gr# + pr#: nu# = .48# * ra# ^ .25#
h# = nu# * k# / d#
REM***h# = 0.01#
PRINT "The Rayleigh number ="; ra#
PRINT "The Nusselt number ="; nu#
PRINT "The convective heat transfer coefficient ="; h#; "W/m^2.K"
qc# = h# * aconv# * (tc# - tinfin#)
REM***calculate radiative transfer from oven walls
qr# = arad# * .0000000567# * (tw# ^ 4# - tc# ^ 4#)
tc# = tc# + (qr# - qc#) / (mcorpse# * cpcorpse#)
walllosses# = qr# + 5000#: tw# = tw# - walllosses# / 3733020#
GOSUB 1000
PRINT "Rate of convective heat transfer, corpses to air ="; qc#; "W"
PRINT "Rate of radiative transfer, walls to corpses ="; qr#; "W"
PRINT "Net transfer to corpses ="; qr# - qc#; "W"
PRINT "Corpse temperature ="; tc# - 273.15#; "C"
PRINT "Oven wall temperature ="; tw# - 273.15#; "C"
PRINT "at time"; s / 60; "minutes"
200 IF INKEY$ = "" THEN 200
REM***FOR dl = 1 TO 20000: NEXT dl
NEXT s
RUN
1000 REM***Find the wall surface temperature
s# = s: fo# = s# / 36091#: xast# = 1#
PRINT "Fo ="; fo#
PRINT "Temperature Theta* at exposed surface =";
FOR n# = 1# TO 40#
theta# = c#(n#) * EXP(-1# * xi#(n#) ^ 2# * fo#) * COS(xi#(n#) * xast#)
thetaast# = thetaast# + theta#: NEXT n#
PRINT thetaast#
twactual# = 293.15# + thetaast# * (1123.15# - 293.15#) - 273.15#
PRINT "The actual temperature ="; twactual#; "C": thetaast# = 0#
REM***Find proportion of energy transferred
FOR xast# = 0# TO 1.009# STEP .01#: FOR n# = 1# TO 40#
theta# = c#(n#) * EXP(-1# * xi#(n#) ^ 2# * fo#) * COS(xi#(n#) * xast#)
thetaast# = thetaast# + theta#: NEXT n#
thetatotal# = thetatotal# + thetaast#: m# = m# + 1#
thetaast# = 0#: NEXT xast#
thetaav# = thetatotal# / m#: m# = 0#: thetatotal# = 0#
PRINT "Proportion of energy transferred ="; 1# - thetaav#
RETURN