re: thermal performance model SOLAR.KnowledgePublications.com

re: thermal performance model
8 feb 2005
here's a simulation for a philadelphia row house using nrel's hourly
typical meteorological year (tmy2) weather data... the first program
makes a simplified output file called "ecayear." 

10 screen 9:key off:cls:pi=4*atn(1)
20 line (0,0)-(639,349),,b:xdf=.073:ydf=3.88
30 for tr= 60 to 80 step 10'temp ref lines
40 line (0,349-ydf*(tr-10))-(639,349-ydf*(tr-10)):next
50 open "13739.tm2" for input as #1'nrel tmy2 file name (phila)
60 open "ecayear" for output as #2
70 line input#1,s$'read header
80 city$=mid$(s$,8,25)
90 lat=val(mid$(s$,40,2))+val(mid$(s$,43,2))/60
100 lon=val(mid$(s$,48,3))+val(mid$(s$,52,2))/60
110 print#2,city$,lat,lon
120 for h=1 to 8760'hour of year
130 line input#1,s$
140 month=val(mid$(s$,4,2))'month of year (1-12)
150 day=val(mid$(s$,6,2))'day of month
160 hour=val(mid$(s$,8,2))-.5'hour of day
170 n=1+h/24'day of year (1 to 365)
180 tdb=val(mid$(s$,68,4))*.18+32'dry bulb temp (f)
190 tdp=val(mid$(s$,74,4))*.18+32'dew point temp (f)
200 pset(xdf*h,349-ydf*(tdb-10))
210 if day=1 and hour=.5 then line (xdf*h,349)-(xdf*h,345)'tick months
220 wind=val(mid$(s$,96,3))*.22371'wind velocity (mph)
230 igloh=val(mid$(s$,18,4))*.317'global horizontal radiation (btu/ft^2)
240 print#2,month;day;hour;tdb;wind;tdp;igloh;
250 idif=val(mid$(s$,30,4))*.317'diffuse horizontal radiation (btu/ft^2)
260 idir=val(mid$(s$,24,4))*.317'direct normal radiation (btu/ft^2)
270 l=pi*lat/180'phila latitude (radians)
280 t=hour'solar time (est)
290 x=-sin(pi*23.45/180)*cos(2*pi*(n+10)/365.25)
300 d=atn(x/sqr(-x*x+1))'sin^-1(x) = declination (radians)
310 w=2*pi*(t-12)/24'hour angle (radians)
320 x=cos(l)*cos(d)*cos(w)+sin(l)*sin(d)
330 thetas=-atn(x/sqr(-x*x+1))+pi/2'cos^-1(x) = sun zenith angle (radians)
340 x=cos(d)*sin(w)/sin(thetas)
350 if x^2 >=1 then phis=-1.570796327#:goto 370
360 phis=atn(x/sqr(-x*x+1))'sin^-1(x) = sun azimuth angle (radians)
370 for phipd=0 to 180 step 90'azimuth angle of plane (degrees)
380 phip=pi*phipd/180
390 x=sin(thetas)*cos(phis-phip)
400 thetai=-atn(x/sqr(-x*x+1))+pi/2'incidence angle to surface (radians)
410 if thetai>=pi/2 then thetai=pi/2
420 rhog=.2'ground reflectance
430 iglop=idir*cos(thetai)+idif/2+igloh*rhog/2'radiation on surface (btu/ft^2)
440 print#2,iglop;
450 'if phipd = 0 then pset(xdf*h,349-iglop)
460 next phipd
470 phip=pi*270/180
480 x=sin(thetas)*cos(phis-phip)
490 thetai=-atn(x/sqr(-x*x+1))+pi/2'incidence angle to surface (radians)
500 if thetai>=pi/2 then thetai=pi/2
510 iglop=idir*cos(thetai)+idif/2+igloh*rhog/2'radiation on surface (btu/ft^2)
520 print#2,iglop
530 next h
430 iglop=idir*cos(thetai)+idif/2+igloh*rhog/2'radiation on surface (btu/ft^2)

here's a sample of its ecayear output file, with the month, day, hour, dry bulb
temperature (f), wind velocity (mph), dew point temp (f) and btu/ft^2 of sun
that falls on the ground and south, west, north, and east walls for every hour
of the year. it starts just past midnight on new year's day...

philadelphia           pa    39.88334      75.25
 1  1  .5  33.98  4.69791   30.92  0  0  0  0  0
 1  1  1.5 33.98  9.172109  30.92  0  0  0  0  0
 1  1  2.5 33.98  10.29066  28.94  0  0  0  0  0
 1  1  3.5 32     23.04213  19.94  0  0  0  0  0
 1  1  4.5 28.04  19.68648  15.98  0  0  0  0  0
 1  1  5.5 24.98  20.80503   8.06  0  0  0  0  0
 1  1  6.5 21.92  19.68648  10.04  0  0  0  0  0
 1  1  7.5 15.98  24.16068 -.940      7.92  26.5    2.6   2.6   41.2
 1  1  8.5 14     25.27923 -5.08     41.5  131.8    9.0   9.0  145.2
 1  1  9.5 15.08  21.92358 -5.08     88.4  205.2   15.6  15.6  154.0
 1  1  10.5 15.08  26.39778 -5.08    83.6  102.9   33.8  33.8   62.7
 1  1  11.5  15.08  25.27923 -4.00   66.2   54.4   34.8  34.8   37.5
 1  1  12.5  17.06  25.27923 -5.08  109.3  150.8   51.9  36.4   36.4
 1  1  13.5  17.06  23.04213 -2.92   99.8  146.3   79.9  32.4   32.4
 1  1  14.5  15.98  25.27923 -2.02   64.3  103.3   81.3  21.8   21.8
 1  1  15.5  15.08  25.27923 -4.00   44.6  118.3  129.9  11.6   11.6
 1  1  16.5  12.02  21.92358 -7.06   10.1   32.1   49.8   3.7    3.7
 1  1  17.5  10.94  23.04213 -7.96  0  0  0  0  0
 1  1  18.5  10.94  23.04213 -7.96  0  0  0  0  0
 1  1  19.5  10.94  23.04213 -5.98  0  0  0  0  0
 1  1  20.5  10.94  18.34422 -5.08  0  0  0  0  0
 1  1  21.5  10.94  18.34422 -4.00  0  0  0  0  0
 1  1  22.5  12.02  20.80503 -2.92  0  0  0  0  0
 1  1  23.5  12.92  20.80503 -.940  0  0  0  0  0
 ...

skipping to the 4th of july...

 7  4  .5  64.94   3.35565  55.94  0  0  0  0  0
 7  4  1.5  64.04  3.35565  57.02  0  0  0  0  0
 7  4  2.5  62.96  4.69791  57.92  0  0  0  0  0
 7  4  3.5  62.06  3.35565  59     0  0  0  0  0
 7  4  4.5  60.98  2.2371   59        1.9    1.3    1.1   1.1    1.4
 7  4  5.5  62.96  2.2371   60.0     21.8   26.7    9.9   9.9   50.7
 7  4  6.5  66.92  5.81646  59       72.2   50.7   22.1  22.1  140.5
 7  4  7.5  71.06  2.2371   60.0    134.4   47.7   34.6  34.6  189.0
 7  4  8.5  73.94  5.81646  53.06   195.9   59.3   46.2  46.2  202.4
 7  4  9.5  75.92  9.172109 53.06   246.6   94.5   54.9  54.9  185.1
 7  4  10.5  78.98  6.93501  53.06  283.7  121.4   61.0  61.0  146.9
 7  4  11.5  80.06  8.05356  51.98  304.3  136.0   64.1  64.1   94.2
 7  4  12.5  82.04  8.05356  48.92  306.5  136.7   94.8  64.7   64.7
 7  4  13.5  82.04  8.05356  50     289.1  122.7  148.4  61.8   61.8
 7  4  14.5  82.94  10.29    51.08  254.8   96.5  188.3  56.2   56.2
 7  4  15.5  82.94  5.81646  50     206.0   61.4  208.0  47.8   47.8
 7  4  16.5  82.04  10.29    51.98  146.1   50.5  198.3  36.9   36.9
 7  4  17.5  80.06  11.63292 53.06   82.7   55.5  153.7  24.2   24.2
 7  4  18.5  75.92  10.29066 51.08   28.8   35.2   68.6  11.9   11.9
 7  4  19.5  71.06   8.05356 53.06    3.8    2.7    3.2   2.1    2.1
 7  4  20.5  68  6.93501  53.06        0  0  0  0  0
 7  4  21.5  66.02001  5.81646  55.04  0  0  0  0  0
 7  4  22.5  64.94  5.81646  55.04     0  0  0  0  0
 7  4  23.5  62.96001  5.81646  57.02  0  0  0  0  0
 ...

here's a program that simulates the house, using the ecayear weather file
as an input. we can try out various improvements to make it more comfortable
and save energy, eg airsealing (line 80), more roof insulation (line 90),
roof reflectance (line 100), a whole-house fan (line 110) and making some
of the south wall into a simple trombe wall (line 120.) trombe walls are
inefficient solar heaters, but they might be cheap retrofits. 

10 screen 9:key off:cls
20 daystart=0'display start time (days)
30 ds=daystart*24'display start time (hours)
40 range=8760'display range (hours)

we can expand the display range above, to look at a few days or hours...

50 line (0,0)-(639,349),,b:xdf=640/range:ydf=3.88
60 for tr=60 to 80 step 10'temp ref lines
70 line (0,349-ydf*(tr-10))-(639,349-ydf*(tr-10)):next
80 ach=.3'air infiltration (house volumes per hour)
90 rroof=40'roof r-value (f-h/btu)
100 refroof=.1'roof reflectance (1=mirror)
110 cfm=2470'whole-house fan cfm (lasko/air king 9600)
120 atrombe=64'area of trombe wall (ft^2, with no int. insul.)
130 aroof=30*40'roof area (ft^2)
140 groof=aroof/rroof'roof conductance (btu/h-f)
150 awall=2*30*16'wall area (ft^2, excluding common and basement)
160 rwall=2'wall r-value (f-h/btu, including windows)
170 aintwall=0'area of interior wall insulation (ft^2, with no ext. insul.)
180 rintwall=10'r-value of interior wall insulation (f-h/btu)

adding insulation inside the brick walls is an option...

190 gint=aintwall/(gwall+rintwall)'int. wall section conductance (btu/h-f)
200 gwall=gint+(awall-atrombe-aintwall)/rwall'wall conductance (btu/h-f)
210 gair=ach*30*40*24/60'air infiltration conductance (btu/h-f)
220 gh=gwall+gair'house conductance (btu/h-f, excluding roof)
230 ch=15000+(2*30*24-aintwall)*5!'house capacitance (btu/f)
240 gascap=50000!'furnace capacity (btu/h)
250 ta=30'initial average outdoor temp (f)
260 th=65'initial house temp (f)
270 open "ecayear" for input as #1:line input#1,h$
280 open "ecaout" for output as #2

this also makes an output file for a spreadsheet...

290 for h=1 to 8760'hours of typical (tmy2) year
300 input#1,month,day,hour,tdb,wind,tdp,igloh,ss,ws,ns,es
310 ta=.99*ta+.01*tdb'3-day moving average outdoor temp (f)
320 tt=86-.32*ta'target temp (f)

the fan target temp varies from 75 in winter to 65 in summer... 

330 gwind=2+wind/2'roof airfilm conductance (btu/h-f-ft^2)
340 tsar=tdb+igloh*(1-refroof)/gwind'sol-air roof temp (f)
350 rgain=(tsar-th)/(1/(aroof*gwind)+1/groof)'roof solar gain (btu)
360 if month>4 and month<9 then ss=0'shade trombe wall in summer
370 tgain=atrombe*(tdb+.9*ss-th)/(1+rwall)'trombe wall gain (btu)
380 q=rgain+tgain-(th-tdb)*gh'passive heat gain (btu)
390 th=th+q/ch'house temp before heating or cooling (f)
400 if tt<65 then tt=65
410 if tt>75 then tt=75
420 qfan=0:qgas=0:fantime=0'initializations
430 if th>tt then p=5 else p=1
440 if th>65 then p=p+2
450 if th>tdb then p=p+1
460 on p gosub 690,770,790,840,840,840,840,850
470 gasheat=gasheat+qgas
480 if qgas>0 then lastheat=0 else lastheat=lastheat+1
490 if qfan<0 then fancool=fancool-qfan:goto 510
500 fanheat=fanheat+qfan
510 fanhours=fanhours+fantime
520 th=th+(qfan+qgas)/ch'new house temp (f)
530 if th>80 then hothours=hothours+1

count the number of hours per year over 80 f, assuming no ac...

540 if hds+range then goto 610
550 'pset(xdf*(h-ds),349-ydf*(tdb-10))
560 'pset(xdf*(h-ds),349-ydf*(tt-10))
570 pset(xdf*(h-ds),349-ydf*(th-10))
580 if qgas>0 then line (xdf*(h-ds),300-qgas/1000!)-(xdf*(h-ds),300)'gas energy
590 line (xdf*(h-ds),320-qfan/2000)-(xdf*(h-ds),320)'mark fan energy

display furnace and fan usage...

600 if day=1 and hour=.5 then line (xdf*(h-ds),349)-(xdf*(h-ds),345)'months
610 print#2,month;day;hour;tdb;th;gasheat,fanheat,fancool
620 'if h>4310 then print igloh,tdb,tsar:end
630 next h
640 close #1
650 print ach,rroof,refroof,cfm,atrombe
660 print gasheat/100000!,fanheat/100000!,fancool/100000!,fanhours,hothours
670 list 80-120
680 end
690 'th<=tt and th<=65 and th<=tdb heat with fan towards tt, then gas to 65
700 qhmax=(tt-th)*ch'to heat all the way to tt
710 qfmax=(tdb-th)*cfm'max heating from fan in 1 hour
720 if qfmax>qhmax then fantime=qhmax/qfmax:qfan=qhmax:return
730 qhmax=(65-th)*cfm'heat to 65
740 if qfmax>qhmax then fantime=qhmax/qfmax:qfan=qhmax:return
750 fantime=(qhmax-gascap)/(qfmax-gascap):qgas=(1-fantime)*gascap
760 qfan=qfmax*fantime:return
770 'th<=tt and th<=65 and th> tdb heat with gas to 65
780 qgas=(65-th)*ch:return
790 'th<=tt and th> 65 and th<=tdb heat with fan towards tt
800 qhmax=(tt-th)*ch'to heat all the way to tt
810 qfmax=(tdb-th)*cfm'max heating from fan in 1 hour
820 if qfmax>qhmax then fantime=qhmax/qfmax:qfan=qhmax:return
830 fantime=1:qfan=qfmax
840 return
850 'th> tt and th> 65 and th> tdb cool with fan towards tt
860 if lastheat<80 then return'don't fan-cool if recently heated
870 qhmax=(th-tt)*ch'to cool all the way to tt
880 qfmax=(th-tdb)*cfm'max cooling from fan in 1 hour
890 if qfmax>qhmax then fantime=qhmax/qfmax:qfan=-qhmax:return
900 fantime=1:qfan=-qfmax:return

1. ach           rroof         refroof       fan cfm       trombe area
   .7            13            .1            0             0
   gas heat      fan heat      fan cool      fan hours     hot hours
   1099.73       0             0             0             298

case 1 above is the unimproved house, which uses about 1100 therms of gas
heat at 100% furnace efficiency and has 298 h/y of indoor temps over 80 f. 

2. ach           rroof         refroof       fan cfm       trombe area
   .3            13            .1            0             0
   gas heat      fan heat      fan cool      fan hours     hot hours
   859.8806      0             0             0             283

reducing air leakage without changing anything else lowers the gas bill
by 22% and reduces the number of hot hours slightly.

3. ach           rroof         refroof       fan cfm       trombe area
   .7            40            .1            0             0
   gas heat      fan heat      fan cool      fan hours     hot hours
   1040.083      0             0             0             207

adding roof insulation without changing anything else lowers the gas bill
slightly and reduces the number of hot hours significantly.

4. ach           rroof         refroof       fan cfm       trombe area
   .7            13            .9            0             0
   gas heat      fan heat      fan cool      fan hours     hot hours
   1122.53       0             0             0             192

making the roof white without changing anything else raises the gas bill
slightly and reduces the number of hot hours significantly.


5. ach           rroof         refroof       fan cfm       trombe area
   .7            13            .1            2470          0
   gas heat      fan heat      fan cool      fan hours     hot hours
   1101.208      7.445226      61.50649      895.0826      121

without changing anything else, a whole-house fan for heating and cooling
reduces the number of hot hours significantly.

6. ach           rroof         refroof       fan cfm       trombe area
   .7            13            .1            0             64
   gas heat      fan heat      fan cool      fan hours     hot hours
   1049.856      0             0             0             297

without changing anything else, a 64 ft^2 trombe wall reduces the gas bill
slightly.


7. ach           rroof         refroof       fan cfm       trombe area
   .3            40            .1            0             0
   gas heat      fan heat      fan cool      fan hours     hot hours
   800.3312      0             0             0             181

combined airsealing and roof insulation saves lots of gas.

8. ach           rroof         refroof       fan cfm       trombe area
   .3            40            .9            0             0
   gas heat      fan heat      fan cool      fan hours     hot hours
   807.8726      0             0             0             131

adding a white roof raises the gas bill slightly and keeps the house cooler.

9. ach           rroof         refroof       fan cfm       trombe area
   .3            40            .9            2470          0
   gas heat      fan heat      fan cool      fan hours     hot hours
   806.815       9.980721      48.20572      806.8318      34

adding a whole-house fan can keep the house a lot cooler.

10.ach           rroof         refroof       fan cfm       trombe area
   .3            40            .9            2470          64
   gas heat      fan heat      fan cool      fan hours     hot hours
   757.1798      8.498159      49.49792      797.0345      33

adding a trombe wall still saves some gas.

11.ach           rroof         refroof       fan cfm       trombe area
   .3            40            .1            2470          64
   gas heat      fan heat      fan cool      fan hours     hot hours
   750.3186      7.776429      52.50024      816.752       41

and leaving the roof black saves a bit more...

with a few more lines we can add internal heat gain from electrical use,
ventilation constraints based on dew point, and so on.

nick