basic tmy2 simulations SOLAR.KnowledgePublications.com

basic tmy2 simulations
12 dec 2003
hordes of 3rd graders used to descend from school buses after school and
take over the terminals at the kiewit computation center at dartmouth and
write their own basic programs. solar house designers can do this too :-)

this program generates a "winter" weather file from a philadelphia weather
data file on an nrel tmy2 cd-rom and displays a graph of outdoor temps for
december and january. tmy2 files are also available from nrel's web site.
the first line of the winter file is the name of the city and its latitude
and longitude. the next lines have a month number (october through january),
day, hour, dry bulb temp, global horizontal sun during the previous hour
in btu/ft^2, and sun that falls on south, west, north and east walls, with
an isotropic sky model.

winprep.bas

10 pi=4*atn(1):screen 9:key off:cls
40 line (0,0)-(639,349),,b:df=0.434
50 for tr=20 to 70 step 10'temp ref lines
60 line (0,349-5*(tr-10))-(639,349-5*(tr-10)):next
80 open "13739.tm2" for input as #1
90 open "winter" for output as #2
100 line input#1,s$'read header
110 city$=mid$(s$,8,25)
120 lat=val(mid$(s$,40,2))+val(mid$(s$,43,2))/60
130 lon=val(mid$(s$,48,3))+val(mid$(s$,52,2))/60
140 print#2,city$,lat,lon
150 gosub 200:pass=1:close #1
160 open "13739.tm2" for input as #1
170 line input#1,s$'ignore header
180 gosub 200
200 for h=1 to 8760'hour of year
210 line input#1,s$
220 month=val(mid$(s$,4,2))'month of year (1-12)
230 if pass=0 and month<10 goto 580
235 if pass=1 and month=2 then end
236 wh=h+8760*pass
240 day=val(mid$(s$,6,2))'day of month
250 hour=val(mid$(s$,8,2))-.5'hour of day
260 n=h/24'day of year (1 to 365)
270 tdb=val(mid$(s$,68,4))*.18+32'dry bulb temp (f)
275 pset(df*(wh-8030),349-5*(tdb-10))
280 igloh=val(mid$(s$,18,4))*.317'global horizontal radiation (btu/ft^2)
290 print#2,month;day;hour;tdb;igloh;
300 idif=val(mid$(s$,30,4))*.317'diffuse horizontal radiation (btu/ft^2)
310 idir=val(mid$(s$,24,4))*.317'direct normal radiation (btu/ft^2)
320 l=pi*lat/180'phila latitude (radians)
330 t=hour'solar time (est)
340 x=-sin(pi*23.45/180)*cos(2*pi*(n+10)/365.25)
350 d=atn(x/sqr(-x*x+1))'sin^-1(x) declination (radians)
360 w=2*pi*(t-12)/24'hour angle (radians)
370 x=cos(l)*cos(d)*cos(w)+sin(l)*sin(d)
380 thetas=-atn(x/sqr(-x*x+1))+pi/2'cos^-1(x) sun zenith angle (radians)
390 x=cos(d)*sin(w)/sin(thetas)
400 if abs(x+1)<.000001 then phis=-1.570796327#:goto 420
410 phis=atn(x/sqr(-x*x+1))'sin^-1(x) sun azimuth angle (radians)
420 for phipd=0 to 180 step 90'azimuth angle of plane (degrees)
430 phip=pi*phipd/180
440 x=sin(thetas)*cos(phis-phip)
450 thetai=-atn(x/sqr(-x*x+1))+pi/2'incidence angle to surface (radians)
460 if thetai>=pi/2 then thetai=pi/2
470 rhog=0.2'ground reflectance
480 iglop=idir*cos(thetai)+idif/2+igloh*rhog/2'radiation on surface (btu/ft^2)
490 print#2,iglop;
500 'if phipd=0 then pset(h-8030,349-iglop)
510 next phipd
520 phip=pi*270/180
530 x=sin(thetas)*cos(phis-phip)
540 thetai=-atn(x/sqr(-x*x+1))+pi/2'incidence angle to surface (radians)
550 if thetai>=pi/2 then thetai=pi/2
560 iglop=idir*cos(thetai)+idif/2+igloh*rhog/2'radiation on surface (btu/ft^2)
570 print#2,iglop
575 if hour=0.5 then line (df*(wh-8030),349)-(df*(wh-8030),345)
580 next h
590 return

here's a simple direct gain simulation using the "winter" file for an 8'
r20 cube with an 8'x8' r2 south window with 80% solar transmission. with
an upper 70 f living space temperature limit, this cube needed 13,300 btu/f
(41% water fill :-) to keep it 60 f min at 8 am on january 25th, after
a 10-day slump. this program graphs indoor and outdoor temps for december
and january and displays the min indoor temp and the time it occurred. 

winsimdg.bas

10 cls:screen 9:line (0,0)-(639,349),,b:df=.43
20 for tr=20 to 70 step 10'plot temp ref lines
30 line (0,349-5*(tr-10))-(639,349-5*(tr-10)):next
40 rv=20'wall r-value
50 c=13300!'house capacitance (btu/f)
60 th=70:thmin=1000'initial house temps (f)
70 open "winter" for input as #1
80 line input#1,h$
90 input#1,month,day,hour,ta,sh,ss,sw,sn,se
100 solgain=0.8*64*ss-(th-ta)*32'net south window solar gain (btu)
110 ih=solgain-(th-ta)*5*64/rv
120 th=th+ih/c'find new house temperature (f)
130 if th>70 then th=70'limit house temp
140 if th-(th-ta)*32/rv then ih=solgain else ih=-(th-ta)*32/rv
120 ih=ih-(th-ta)*5*64/rv
130 th=th+ih/c'find new house temperature (f)
140 if th>70 then th=70'limit house temp
150 if th-(ti-ta)*32/rv then ih=solgain else ih=-(ti-ta)*32/rv
120 ih=ih-(tc-ta)*64/rv-(ti-ta)*4*64/rv
130 tc=tc+ih/c'find new ceiling temperature (f)

if the cube above is only occupied for 8 hours per day, the simulation
fragment below indicates that it only needs 274 btu/f of ceiling mass.

winsimcs.bas

...
50 c=274'ceiling capacitance (btu/f)
90 if tc<65 then ti=tc else ti=65'room temp (f)
100 solgain=0.8*64*ss-(tc-ta)*32'potential net south window solar gain (btu)
110 if solgain>-(ti-ta)*32/rv then ih=solgain else ih=-(ti-ta)*32/rv
120 ih=ih-(tc-ta)*64/rv'heatflow from ceiling (btu)
130 if hour>8 and hour<16 then ih=ih-(ti-ta)*4*64/rv'occupied hours
140 tc=tc+ih/c'find new ceiling temperature (f)

each program will fit on a single page. the assumptions are clear, and the
programs are simple to understand, except for the first one, which might be
rerun with a different ground reflectance (line 470.) they might be expanded
in many ways, or run with 30-year hourly vs tmy2 files, or become part of
a spreadsheet in visual basic...

nick