Computer simulation for a solar heating system

Computer simulation for a solar heating system

Mathematics and Computers in Simulation XX (1978) 114-127 0 North-Holland Publishing Company COMPUTER SIMULATION FOR A SOLAR HEATING SYSTEM Adel H...

1016KB Sizes 1 Downloads 75 Views

Mathematics and Computers in Simulation XX (1978) 114-127 0 North-Holland Publishing Company



Adel H. ELTIMSAHY College of Engineering,

The University of Toledo, Toledo, Ohio 43606,


and Clifford H. COPASS, Electrical Engineering Department,

Carnegie-Mellon University, Pittsburgh, Pennsylvania 15213,


In this paper a simulation model for a basic solar heating system which is suitable for control and optimization studies B developed. A model for flat plate solar collectors using fluids for the heat transfer which allows the analysis of operation in great detail is presented. The simulation model for the solar collectors has also the characteristic that its inputs and outputs are quantities which are easily measurable in the real world. The effect of temperature stratification in the heat storage device model is considered by simulating two storage chambers of variable volume. The heat storage model requires relatively small time increments, a requirement which is not as stringent as that for the solar collectors. The tank’s model primary variables are the easily measurable quantities of temperature and flow rates with the addition of volume as an internal variable. Weather, load and control modeling is also presented. The gradient search method is then applied to study optimal sizing of solar heating components for The University of Toledo solar house.

1. Introduction Several off-line computer simulation studies on solar heating systems have been performed by many research groups. Examples of these investigations have

been published by groups at the University of Wisconsin and Colorado State University [l-4]. Most of these simulations perform their calculations over roughly hour long intervals of system operation. Excessive computer time becomes a disadvantage if a

Fig. 1. Solar collector/storage

tank system.


A.H. Eltimsahy, C.H. &pass /Solar heating system

small time increment is used. To study control effects in detail requires small time increments. In this paper a computer simulation model for a basic solar heating system is developed which has several advantages. First, the heat storage device is made of two chambers of variable volume. Only a moderate size minicomputer is needed to implement this simulation. This computer model with its small time increments is practical when used for “on-line” control and optimization studies using an existing system with similar characteristics as a “reference” model

[51* The heart of any solar heating system is the solar collector which collects and converts solar radiation. The simplest, least expensive and most common collectors are of the flat plate variety and are used in this study as shown in fig. 1. The heat collected by this device is carried off by a fluid and either used directly or stored for future use. Such heat can be put to a great many uses, the most obvious being simple space heating. The next major device is the heat storage tank. Such a tank can be designed so that cooler fluid from the bottom can be drawn out and sent through the collectors to be heated and returned to the upper part of the tank. Over a period of time the entire mass of fluid can be heated and made available to a load.

FR = a collector efficiency factor, G = radiant energy (BTU/ft* hr), Q, = heat output (BTU), t, = outside ambient temperature (OF), li = collector inlet temperature (OF), tin = mean collector plate temperature (OF), U, = a collector loss coefficient (BTU/ft* hr “F), a! = the absorption coefficient of the plate, 0 = time (hours), 7 = the transmission coefficient of the glass. To begin converting equation (1) to the desired form, it is necessary to convert energy terms into relations using temperature. The heat output rate may be written as:

dQu _ But Cp H’(to - ti) = FRA~GcYT - FRACUL(ti - ta)

Therefore dto z




[Cp W(ti - to) + FRAC

X UL(ta - ti) + FRAcGcu][~AcC’~] -’ ,


2. Flat plate collector or Flat plate collector models are fairly common and easy to find in published literature. Many of these models are based on heat flow rates with constant flow of heat exchange fluid when operating [6,7]. In this section a modified simulation model is derived from existing models which will allow variable fluid flow rate. The model also is to have temperatures as the primary variables instead of heat which will make it suitable for control studies. The basic relation of conservation of energy for a flat plate solar collector may be stated as follows [6] :

dQu _ -&-FRA(JGO!T-

UL(ti - fJ] -A(JCi d$

where A= = collector area (ft’), 6; = heat capacity of collector (BTU/ft*“F),



dt, _ z - 2 [Cp W(ti - to) f FRA~GcIT + UL(ta - ti)][A&]

-’ -2,


where t, = Cto+ ti)/2, C, = fluid specific heat (BTU/lb OF), t0 = collector outlet temperature (OF), IV = fluid flow rate (lb/hr). All of the parameters and values in equation (4) can be measured or quickly calculated except FR. This efficiency factor can be found from the physical design of the heat absorbing plate and water tubes. It also depends on the fluid characteristics and fluid flow in the water tubes [8]. The factor FR can be found as follows. Consider

A.H. Eltimsahy, C.H. Copass /Solar heating system


1 -+I t+-d Fig. 2. Collector





plate cross section.

2. Let:

W, - d



[email protected]) bL


’ @[email protected]/[email protected]

hi = 0.023


p=(-+ FR=s

~RH)o.2$.47 y

wt 4_



(d/I+‘,) + El - (d/I+‘,>] F ) [I -exp(-q)],

(8) -1 ’



where d = external tube diameter (ft), do = internal tube diamter (ft), K = plate thermal conductivity (BTU/ft hr “F), k = fluid thermal conductivity (BTU/ft hr “F), M = fin thickness (ft), RH = area of flow divided by perimeter of flow (ft), IV, = fluid flow rate per tube cross sectional area (lb/ft2 hr), W, = collector tube spacing (ft), /J = fluid viscosity (lb/ft hr). Relation (10) for FR was published by R.W. Bliss [8] as previously mentioned. Actual simulating a collector’s behavior can be easily accomplished using a computer to perform calculations over small discrete intervals of simulated time. Most simulation routines presently available perform their calculations over roughly hour long intervals of system operation. This is good for many uses but has great limitations when the effects of various control methods are to be compared. Generally, any digital simulation requires a constant control action over the full duration of the calculation interval. Thus, to study control effects in detail requires small time

increments allowing the control to change often. In real life it is not uncommon to find partly cloudy or other marginal conditions requiring even a simple on/ off controller to cycle every few minutes. It was desired when developing this set of simulation routines to study the effects of control, thus it was desirable to use small time increments. Over small intervals, discrete representation approaches the accuracy of a differential equation. Unfortunately over longer intervals, a discrete equation may become unstable, developing wide swings in the amplitudes of various values. The relations previously given have exhibited this’property. Over small time increments they give very good results, but over larger time intervals increasing instability is observed. One immediate disadvantage of this kind of behavior is the excessive computer time needed to perform extended simulations. Fortunately, only a moderate size minicomputer is needed to implement this simulation. Such a minicomputer is often available for the extended periods of time needed. Programming a computer to simulate a solar collector using these equations can easily be accomplished as follows. For purposes of efficiency, it is desirable to place the value of those parts of equations constant from iteration to iteration into a storage location for immediate use without recalculation. Other factors dependent upon temperature or flow rates must be recalculated for each iteration. Depending upon the fluid used, it may be necessary to calculate its viscosity as a function of temperature. Ethylene glycol in a 50 percent mix with water exhibits significant viscosity variation. This variation can be fitted with an equation by several standard means to the accuracy desired. One such equation is given below. /J = 25.9652 - 0.363198 - 0.380783 p=l

X lo-’


t, + 0.196965

X 1O-2 t”,

t: for t, < 215°F ,

(114 (1 lb)

A constant can be calculated for the physical contiguration of the collectors to be modeled. B, =

0.023 k213C;13 A$8(4RH)o.2


where AT represents the total tube cross sectional area

A.H. Mtimsahy, C.H. Copass /Solar heating system

(ft’), Two calculation point.

steps can be made from this

B,I = B,//.f’.47 ,

hi = w0.8B, 1

Two more constants B2



Two more constants B4

(13, 14)

are then needed:


= CpIh.&


B5 = - ULAcICP = 1/B4.





Equations step:

are needed for the next step.

Equations (17) through (19) are used to calculate FR:




[l -(d/W,)]




Three more constants step:

(15) and (16) allow another calculation



Bz/hi + 83 *

are needed for the next major

Be=%, AcCa



fluid Cal al ate .P cvi scosi t )

tm = tm+ At*


/ Return




to+ &to


Fig. 3. Collector simulation routine flow chart.


A.H. Eltimsahy, C.H. Copass /Solar heating system



Equations (2) through (23) are used to calculate At,. At, = {BbW(ti - to) + FRB7 + FRB7 [CBS + U~(ta - ti)]} A0 - Att .


Now some record keeping is needed prior to the next iteration. t, = (to + ti)/2 .

t, = t, + At, ,


Considerable simplification results for the special case of zero fluid flow. For this condition, collector temperature becomes more uniform and can be represented by t,. energy outflow due to fluid transport becomes zero. At, = [GBs - UL(tm - ta)] $. a


With necessary record keeping: I, = t,,, + At,


1, = t Ill.


This sequence of calculation can be written into a computer program which repeatedly for each time interval, calculates weather conditions and building load, performs the calculations for modelling a collector and then simulates the action of a heat storage device. The modified collector simulation model would require as inputs during each iteration: t,, W, ti, t,, t,, Ati and G. Of these, At, and t, can be maintained internally, At, and G must be supplied by weather data or simulation, ti and Ati can be supplied by a storage simulation and W must be supplied by the control method used. Other values for A0 and physical characteristics can be supplied once per simulation run. Outputs are At, or At,,, which can be used to maintain the desired variable t,. A flow chart of the subroutine used to perform this simulation is outlined in fig. 3.

3. Heat storage tank Most of the usual storage tank models use one or more well mixed chambers. It is quite difficult to model the temperature differences from place to place in a real storage tank; however, it is possible to rough-

ly simulate the effects of temperature differences by assuming several chambers, each at a uniform temperature. Each chamber is considered well mixed to produce the uniform temperature desired. The simplest model is of course a single well mixed chamber using conservation of energy to determine the temperature. This kind of model is excellent for a heavily agitated storage tank or for an approximate tank model needing a minimum of computational effort. Generally however, any real heat storage tank will be significantly hotter at the top than at the bottom. This, of course, is due to the changes in density that water undergoes when heated. This effect can be put to good use by removing collector inlet water from the cooler bottom of a tank and removing warmer water from the top to drive a heating load. Using the coolest collector inlet water available lowers the average collector plate temperature allowing a minimum of heat loss. Drawing the hottest water for the load keeps the load supply hotter for a longer period if cool water is returned to the bottom and not allowed to mix with the hotter water above. Certainly the average temperature goes down by this method as heat is removed; however, not allowing the tank to cool uniformly makes best use of the heat available. The usual method of modelling this change in temperature vertically in the tank is to assume two or more cross sectional chambers of fHed volume. Each chamber is generally isolated from the others except for fluid flow one to another caused by externally pumping fluid out of one chamber and into another. Conservation of energy is used with each chamber to determine its temperature. This method generally is limited to three chambers for reasonable results. When more chambers are assumed, it becomes difficult to properly evaluate the effects of the heat loss out of the top chamber. This chamber generally has a disproportionately large heat loss area due to the top cover on the tank and thus the isolated chamber at the top may show excessive heat loss and temperature decline. Simulations using this method with three chambers often show two upper chambers at almost the same temperature with the lower chamber substantially cooler. For this reason and also because the real tank to be used would be well mixed when the collector pump was running, it was decided that the model to be developed should use no more than two chambers.

A.H. Eltimsahy, C.H. Copass /Solar heating system

A standard hot water heater tends to develop uniformly hot water when its heating unit is in operation and then as water is used, it developes an upper area of relatively uniform hot water and a lower area of uniformly cold wate-:. As more water is drawn off, the hot section becomes smaller and the cold section larger as the dividing line between moves upward in the tank. This kind of behavior seems completely reasonable for a solar heat storage device also. During a sunny day the collector pump will operate allowing the tank to receive hot water and also agitating the tank, all of which tends to produce a nearly uniform heated condition. Later, as hot water is drawn off, a cold area would develop near the bottom which would enlarge and rise up through the tank as hot water is continually drawn off. It was decided to approximate this behavior by assuming a two chamber tank with a variable partition between the chambers. The total volume of the two chambers would be required to equal the total tank volume at all times but one chamber would shrink and the other expand as water is pumped out of one and into the other. In addition, it was assumed that incoming water if introduced into hotter water would tend to sink through the less dense water with little mixing until it reached a layer of water roughly at its own temperature. This was simulated with a pair of valves on the collector return line which would send the water to whichever chamber was cooler than the incoming water with a preference to the upper chamber. Chamber size changes would be controlled by the overall fluid transfer between chambers during the selected time interval. As before with the collectors, it was desired to have the primary variables as flow rate and temperature rather than energy because they are quantities easily measurable in the real world. The tank model would also be used with the short time intervals that the collector equations need; thus, the tank model could also use discrete rather than differential quantities. The tank representation used is illustrated by figs. 4 through 6. If fluid transfer in one direction outweighs the other for an extended period of time, one chamber may increase in size until it includes nearly all of the tank volume while the other chamber shrinks toward zero. To prevent negative volumes from developing, it is necessary to stop the changing boundary


from moving further. This is illustrated by figs. 5 and 6. The changes in flow configuration shown prevent the smaller chamber from shrinking any further but encourages it to reenlarge when conditions permit. When simulating a sunny day, the upper chamber will enlarge as heated fluid from the collectors fills it. Eventually the upper chamber fills the entire tank and a single well mixed chamber forms the effective simulation (see fig. 6). Later, as heat is drawn out, the upper chamber shrinks while returned cool water enlarges the lower chamber. Eventually the lower chamber fills the entire tank meaning that little or no heat remains (see fig. 5). Computing the volume of each chamber after every time increment is done very simply by adding the volumes of incoming flows and subtracting outgoing flows. V 1new = Vlo,d+AV,=


V1,,ld+------, 62.44

where VI = volume of the upper tank chamber (ft3), V, = would represent volume of the lower tank chamber (ft3), WL = fluid flow rate to the heat load (lb/hr), e = time (hours). The only requirements on equation (30) is that the flow rate WL must be constant over the entire interval and that AV, must not be too small compared to VI. This last requirement is to prevent excessive loss of accuracy in the calculations. Temperature changes due to flows can be computed as follows:

v, &jtl t1new =

.,d62.44 + WAet,

Vr62.44 + WA0


where tl = temperature in the upper tank chamber (OF), t2 = would represent temperature in the lower tank chamber (“F). The approximate heat loss from each chamber may be calculated as follows: AQi = RT AOA,i(tj - ti) , where Qi = heat from a chamber (BTU), RT = heat loss coefficient (BTU/ft* hr “F),



A.H. Eltimsahy, C.H. Copass /Solar heating system

Chanber t1 /


Chanber 12’

2 v2



J -From



Fig. 4. Normal storage tank model.

TO Load





Fig. 5. Storage tank model chamber 2 large.

Fig. 6. Storage tank model chamber 1 large.


A.H. Eltimsahy, C.H. Copass /Solar heating system

rj tI,

= temperature in a chamber (OF), = inside ambient temperature (OF), Ari = the total outside area of one chamber of the tank (ft’). The temperature change caused by the heat loss may now be calculated as follows:


= tjold


AQi cp



The next step is to compute the effects of the collector circulation upon the volumes and temperatures. 1;


set Vi = Vr -J2




set Vi= V2+J2PC1


set ti = tr


set t2 = tl


If ISW2 = 9 or -1 ; set Vi = Vi +Jr PC1 The temperature loss term should ideally be quite small and may become smaller than the computer can accurately subtract. On the other hand, heat loss may altogether be of minimal significance in many cases and good accuracy may not be needed. In some cases it may be desired to calculate the effects of heat conduction through the tank and internal mixing between chambers. Heat conduction through the walls and fluid can be roughly estimated but generally has such small effects that it can be ignored. Internal mixing on the other hand can have a more pronounced effect in many situations, however it is extremely difficult to make more than an educated guess as to the degree of its action and how to best represent it. A solar heating system should be designed to minimize hot water storage mixing for best efficiency, thus mixing should not be an important terms in a well designed system. Many variations exist for implementing the previous equations in a computer program. Calculations may be performed in several possible order of sequence and newly calculated values may be immediately substituted or saved and substituted later. The sequence given below is but one of many possibilities. The first operation chosen was to set the valves (Jr and J,), computational switch (ISW2) and precompute frequently used quantities. If t, > tl , If t, < t 1 ,

set Jr = 1, set J2 = 1 ,

J, = @J J1 = @


PC1 = IV Al3/62.44,


PC2 = WL A0/62.44,

(34) (35)

If V2 large, set ISW2 = -1 , If normal, set ISW2 = $J , If Vr large, set ISW2 = 1 .


set Vi= V2-Jr PCl


set ti = t* ,


where ti


collector inlet temperature


f,V, tJrPClt, t11 =


Vr tJrPC1 t* V, + J*PCl t,




VI = v; )



v, = v;

Step three computes

(47,4S) the effect of heat losses.

AR = 2/r ,


A, =ART/‘, +AE,


At, = t,, =


- tl)


A,=ARVz+AE, At, = t2* =

(51) (52)

RTAOAT(t; - t2) Cp V,62.44


Step four is to find the effect of heat load. t 23

V&62.44 =

t wLAotLR

V262.44 + WLAO

If ISW2 = # or 1 ; (36)



set V, = VI - PC2


set V2 = V2 t PC2


where tLn represents heat load return fluid temperature (“F). The fifth step is included to reduce accumulated computational error in the chamber volumes.

A.H. Eltimsahy, C.H. Copass /Solar heating system




Vtot -

Vl v2 = V*ot

v, + v, ’



where Vtot represents total storage volume (f?).




Step six uses the previously calculated and stored temperature changes to calculate the new rate of temperature change going into the collectors for the collector simulation.



Calculate Calculate









Chanber 1 Calculate




Changes Calculate Chage Lke

Temoerature in







Fluid Calculate





Changes 1 \1 Calculate Outside




tank simulation


Upper Area

flow chart.


A.H. Eltimsahy, C.H. Copass /Solar heating system

[email protected]

set At, = tzl + tzz + t23 - 2tz .


(59) IfISW2=

set Atr = tll + tlz - II .



Step seven is similar to step 6 but is used to compute the new chamber temperatures. t, = t11 + t12 >


t2 = t21 + t22 + t23 - 12 .


Step eight sets the value for the heat load supply temperature tLO. If ISW2 = -1 ,

set tLo = t2 ,

If ISW2 = $Jor 1 ,

set tL0 = t1 .


Figure 7 shows the storage tank simulation routine flow chart. The primary variables needed at the beginning of each iteration are: t,, tLR, W, WL, tl, t2, V,, V2 and t6. Of these, tl, t2, V, and V2 can be maintained internally; t, and W must be supplied by the collector simulation and its control; tLR and WL must be supplied by the heat load simulation and t: must be supplied by some of weather and ambient condition generator. Primary outputs are tLO for the heat load simulation and ti with Atr for the collector simulation.

4. Weather, load and control modeling Any solar heating system simulation requires weather and heat load driving functions. Sun position, cloud cover, outside temperature, wind and even inside temperature affect the operation of a solar heating system. Sun position and interior temperature follow well known predictable patterns; however, cloud cover, wind and outside temperature have no regular predictable pattern but vary greatly and often irrefularly. Cloud cover and exterior temperature are particularly important because cloud cover determines how much solar energy is available to ground based units and outside temperature is the primary determining factor for how much heat energy is needed. Generally, solar energy intensity is used as a driving function when modeling solar heating systems thus combining the effects of sun position and cloud cover. Two major methods are used to supply solar inten-


sity and outside air temperature to a simulation program. These methods are tabular data or function generation. Combinations of these also exist. Tabular data means that for each iteration or group of iterations, the weather conditions to be used are recorded and only need to be read off. This method requires extensive lists of detailed data which may consume considerable storage space. On the other hand, an equation to generate weather data may be quite compact but is usually of questionable accuracy. The accuracy problem for both methods is compounded by the fact that the United States Weather Bureau keeps essentially no useful solar intensity data. Thus, the normal source of detailed weather data can supply only half of the needed information. Occasionally a private group will have records of solar intensity for an area; however, this is seldom the case. If no external sources for this data can be located, it is by far the best to simply begin making accurate solar intensity measurements in the area under study. Visual estimates of solar intensity are very difficult to make with any accuracy but are better than nothing at all. The appearance of graphically recorded solar intensity in Toledo is that of a positive going half rectified sine wave of random amplitude for the gross daily pattern added to a strong negative going noise signal representing partly cloudy periods. Such a complex pattern could cause great difficulty when modeling a solar heating system. It has been shown, however, by D.J. Close [6] that these rapid intensity fluctuations caused by broken cloud cover tend to average out and produce nearly the same effect as a uniformly hazy day with the same total energy availability. Thus average or smoothed data is adequate for almost all simulations. One method of simulating solar intensity that has several advantages is to use a function such as a half rectified sine wave with amplitude and duration parameters recorded in a data tile on a daily basis. This method can simulate the random day to day variations with reasonable accuracy without the need for extended tabular storage of point by point records. Exterior temperature simulation has the simplifying advantages of less abrupt changes and detailed records are available at the Weather Bureau. This makes it possible to keep stored tables of the temper-


A.H. Eltimsahy, C.H. Copass /Solar heating system

ature with a few hours between points. Because temperature also has a more or less regular fluctuation through the day, it is possible to fit a function here also with parameters stored on a daily basis. Either method can reasonably simulate the semirandom nature of temperature changes. House heat load is another very important unpredictable quantity. It depends on house construction, wind, outside and inside temperatures, sunlight available through windows and house usage patterns. All of these items are important; however, the temperature difference inside to outside is the overriding factor when considering a well constructed dwelling. Such things as heat produced by cooking and heat lost by opening doors can be very important but are very difficult to simulate and often tend to cancel out. One simple method of determining house heat load is to relate it proportionately to the difference between exterior and interior temperatures. This is reasonable for most simulation pruposes but ignores surge load conditions that can be caused by house usage patterns. The American Society of Heating, Refrigeration and Air Conditioning publishes considerable information on this subject. Interior temperatures may be programmed to follow a day-night pattern to conserve heat. When the house occupants are asleep, the temperature can be lowered several degrees to conserve heat. This technique is sometimes used with standard heating systems. In many cases it is desired to simulate the behavior of a solar heating system during a reasonable typical period of time many times over to compare the performance when various component parameters are changed. It may be desired to find the system configuration best for the local house heat load and weather conditions. This brings up the problem of choosing a typical period of weather for use with the simulation. One period of weather may favor one system configuration over another and to find a best system it is necessary to find a typical weather period. This problem can be approached using statistical methods but it is important not to average out all of the variations from day to day. Once the range of potential system parameters has been reduced, it is useful to compare different typical periods of weather to find the truly best system. Combining weather and load simulation with a cd-

lector and storage simulation nearly completes a simple simulation package. Time periods such as one minute can be selected as the incremental time Ae. Simulation can be performed by selecting some set of initial conditions and then repeatedly in order calling the weather simulator, collector simulator, storage simulator and load simulator. One item of extreme importance not mentioned so far are the controllers. It is necessary to have a flow rate controller between the storage device and collector. Another controller is needed between the house and the storage device to keep the house temperature as desired. This is where having the primary variables to be flows and temperatures helps the most. Whatever control proposed for an actual system can be simulated directly on the computer. These controls can be simple on-off with predetermined flow rates, proportional flow or even a complex form such as derivative or optimal control. All of these controls essentially specify the flow rate of a hot water system as a function of temperatures and other possible controlling factors such as time of day. In the simulation, the heat laod controller program between house and storage has the additional task of specifying the return water temperature. This is done by simple conservation of energy. Once the load flow rate, supply temperature and heat requirements are known, it is a simple matter to calculate return water temperature. Auxiliary heat control may be independent or tied into the heat load control also. Several options are available here also. Auxiliary heat can totally take over when solar stored heat is exhausted, it can continuously supplement the solar heating system by various amounts or even a timed control may be used for various purposes. One simple collector control law is to turn on a pump with predetermined flow rate when the collector is a few degrees warmer than the storage tank and to turn it off when the collector is equal or cooler in temperature. The heat load control must use heat availability and house heat requirements as a guide while the auxiliary heat control supplies the difference between required heat and solar supplied heat. Adding the controllers completes the simulation presented. A computer can simulate several days of solar heating operating in only a few minutes. Making

A.H. Eltimsahy,

Cl ear

C.H. Copass /Solar




Parmeters end

and ’


Fi 1 es


I terstion


heating system









I -1

Add Energy to



For House Load, Supplied






Locations Solar and






1 Simulate 1













Simulation and

Data, Correct

Clear Possible

Summing Clock



Fig. 8. Off line simulation

5. Preliminary


Hour? Print

the load simulation more for example a heat pump or preheating of domestic contains a diagram of the cess.


complex allows simulating as an energy transfer agent hot water supplies. Figure 8 complete simulation pro-


A number of off line simulation runs were performed for The University of Toledo solar house.

flow chart.

These were performed using the gradient search method. A nine day period was selected to form the weather data base. These nine days included a period of unusually severe winter weather intermixed with more moderate but more cloudy winter weather which altogether would tax any heating system to the limit. Auxiliary heating rates were based on a recent electricity bill assuming resistance heating (0.0388 dollars/kilowatt hour). Solar heating equipment cost was priced at current rates (seven dollars per square foot of collector with representative tank and antifreeze


A.H. Eltimsahy, C.H. Copass /Solar heating system

combination of costs proved to be remarkably stable regardless of the amount of solar energy equipment assumed. It was found, however, that solar heating could slightly reduce the total cost. This cost reduction increased with collector size until approximately 550 square feet of collector with a 150 cubic foot tank was used. Such a large collector could not be placed on the current roof without significant modifications. Cheaper energy such as that currently supplied by natural gas is less expensive than solar equipment; however, it is in very short supply. In the near future as fuel costs rise and solar heating equipment becomes cheaper, a much more significant cost benefit will be realized.

Table 1. Weather data for simulations Date

l-lo-72 1-11-72 l-12-72 1-13-72 1-14-72 1-15-72 1-16-72 I-17-72 1-18-72

Temperature (“F) Mean


Max solar intensity (BTU/hr ft*)

39 31 35 28 8 -7 -3 22 39

45 39 44 44 21 -2 11 33 46

298.6 26.5 199.0 22.0 321.0 221.0 321.0 336.0 17.0

Day length (hr) 7.5 8.5 7.0 7.0 8.5 8.75 9.0 9.0 6.0

costs) assuming a twenty year life and an eight percent loan. Thus, 13% of the initial cost per year must be absorbed. The house load was calculated to be 2500 BTU/hr “F. This value, when multiplied by the time interval in hours and the difference between indoor and outdoor temperature, gives the amount of heat required during that time interval. Table 1 lists the weather data used by the simulation for each day of the run. Table 2 shows some of the gradient search results for various tank and collector sizes. Base size is the component size used to generate a reference cost. Increment is the amount that the component size was changed to test the effect on the total heating cost. Cost reduction ratio is the relative amount that the total heating cost was reduced when the increment value was added to the base component size. This

6. Conclusion This paper dealt with the mathematical modelling and simulation of a solar heating system. The model developed uses temperatures as variables instead of energies which would make it suitable for control and optimization studies. The effect of temperature stratification vertically in the heat storage tank was approximated by simulating two chambers of variable volume. Preliminary off line simulation runs were performed for The University of Toledo solar house, and the gradient search method was applied in order to find the optimal sizing for the solar heating components such as the solar collector and the heat storage device.

Table 2. Simulation results Collector


Base size (ft*)

Increment (ft*)

500 600 800 1000

50 50 100 50

Cost reduction ratio 0.011317 -0.02622 -0.0256 -0.0444

Total cost

Base size (ft3)

Increment (ft3)

150 150 200 300

20 20 20 10

Cost reduction ratio -0.0804 0.0188 -0.052 -0.2239

277.64 280.16 284.30 288.60

A.H. Eltimsahy, C.H. Copass /Solar heating system


[l] L.A. Butz, W.A. Beckman and J.A. Duffie, Simulation of a Solar Heating and Cooling System, Solar Energy, Vol. 16,pp. 129-136,1974. [2] S. Klein, P. Cooper, T. Freeman, D. Beckman, W. Beckman, and J. Duffie, A Method for Simulation of Solar Processes and its Application, Solar Energy, Vol. 17, pp. 29-37,1975. [ 31 S.A. Klein, W.A. Beckman, and J.A. Duffie, A Design Procedure for Solar Heating Systems, Solar Energy, Vol. 18,~~. 113-127,1976. [4] C.B. Winn and G.R. Johnson, Dynamic Simulation for Performance Analysis of Solar Heated and Cooled Build-






ings, presented at the American Society of Mechanical Engineers Annual Meeting, New York, November 17-22, 1974. A.H. Eltimsahy, and C.H. Copass, On-line Simulation of Solar Heating Systems, Simulation, Vol. 26, No. 5, May 1976. D.J. Close, A Design Approach for Solar Processes, presented at the Solar Energy Conference, Boston, Massachusetts, March 20-22,1967. N.R. Sheridan, K.J. Bullock, and J.A. Duffie, Study of Solar Processes by Analog Computers, Solar Energy, Vol. 11, No. 2,1967, pp. 69-77. R.W. Bliss, The Derivation of Several Plate-Efficiency Factors Useful in the Design of Flat-Plate Solar Heat Collectors, Solar Energy, December 1959, pp. 55-64.