After seeing many of your posts, I expected the proper BASIC to be included
in your response.
The program below changes units, ignores lots of irrelevant data,
calculates sun on 4 walls, and makes a "winter" file that starts
in October. TMY2 files have dew point temps instead of humidity
ratios (w). This program doesn't happen to read those.
Article 16881 of alt.solar.thermal:
From:
[email protected] (Nick Pine)
Subject: BASIC TMY2 simulations
Date: 12 Dec 2003 07:35:25 -0500
Organization: Villanova University
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
F=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
ASS=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
F=.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<THMIN THEN THMIN=TH:MMIN=MONTH
MIN=DAY:HMIN=HOUR
150 IF MONTH=10 OR MONTH=11 GOTO 210'skip Oct and Nov display
160 SST=SST+SS'accumulate south wall sun
170 PSET(DF*H,349-5*(TA-10))'plot ambient temp
180 PSET(DF*H,349-5*(TH-10))'plot house temp
190 IF HOUR=0.5 THEN LINE (DF*H,349)-(DF*H,345)
200 H=H+1
210 IF MONTH=1 AND DAY=31 AND HOUR=23.5 THEN GOTO 220 ELSE GOTO 90
220 PRINT MMIN;DMIN;HMIN;THMIN;C;SST/62
The indirect gain version below puts an R20 wall between the cube and
the window and allows airflow between the living space and the "sunspace"
during the day and stops the airflow at night. With an upper 70 F living
space temperature limit, this cube needed 3050 Btu/F of capacitance to
avoid dipping below 60 F at 10 AM on 12/28.
winsimss.bas
....
40 RV=20'wall R-value
50 C=3050'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'potential net south window solar gain (Btu)
110 IF SOLGAIN>-(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<THMIN THEN THMIN=TH:MMIN=MONTH
MIN=DAY:HMIN=HOUR
160 IF MONTH=10 OR MONTH=11 GOTO 220'skip Oct and Nov display
170 SST=SST+SS'accumulate south wall sun
180 PSET(DF*H,349-5*(TA-10))'plot ambient temp
190 PSET(DF*H,349-5*(TH-10))'plot house temp
200 IF HOUR=0.5 THEN LINE (DF*H,349)-(DF*H,345)
210 H=H+1
220 IF MONTH=1 AND DAY=31 AND HOUR=23.5 THEN GOTO 230 ELSE GOTO 90
230 PRINT MMIN;DMIN;HMIN;THMIN;C;SST/62
This modification below puts mass in the ceiling where it can be warmer than
the living space and maintain the living space at a (mostly) constant 65 F
temperature. The "sunspace" has the ceiling mass temperature when air flows,
which may be a conservative assumption. The simulation indicates that this
cube only needed 616 Btu/F to avoid dipping below 60 F on 12/28. According
to the simulation, the ceiling mass was 242 F on 1/15...
winsimcm.bas
....
50 C=616'ceiling capacitance (Btu/F)
90 IF TC<65 THEN TI=TC ELSE TI=65
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-(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
