- Email: [email protected]
An improved computer code for the simulation of solar heating systems
Solar Energy. Vol. 30, No. 2, pp. 9%108, 1983 Printed in Great Britain.
0038~92X/831020099-10503.00t0 © 1983 Pergamon Press Ltd.
AN IMPROVED COMPUTER CODE FOR THE SIMULATION OF SOLAR HEATING SYSTEMS P. B. HOWELLSand R. H. MARSHALL Solar Energy Unit, University College, Cardiff CF2 1TA, Wales
(Received 12 March 1982; accepted 16 April 1982) Abstraet--A generalised solar energy simulation code, which represents a significantimprovement over fixed step iterative simulation codes, is presented. The key features of the code are the simultaneous solution of all differential equations arising from closed loops in the system, and the technique for advancing the simulation from "control function" to "control function". Computer time is minimised by the use of sparse matrix techniques together with automatic time step control generated from an error bound estimate. The application of the improved code to the simulation of solar energy systems is demonstrated on two example problems. In the first problem, emphasis is placed on predicting the temperatures of the fluids within the system and the controller actions, while in the second problem the emphasis is on predicting long-term system performance. The simulation results of both these problems demonstrate that the improved code is suitable for both detailed system control and long-term system performance studies.
The governing system differential equations are formulated as sets of coupled differential equations to form a state matrix which is solved directly (i.e. without iteration) using a highly stable integration scheme codenamed RIKM[2, 3], with provision for alternative explicit and implicit integration schemes. Soft non-linearities, such as temperature dependent material properties, are treated by computing the state matrix equation at each time step. Automatic time step control generated from an error bound estimate is used to maintain the accuracy of the simulation to within a user-specified error tolerance. The choice of the error tolerance, as will be demonstrated, is usually a compromise between computer processor time and the closure of the final energy balances. The energies are determined within the program by assuming a constant specific heat capacity and then numerically integrating the computed temperature time trajectories using the trapezoidal rule. This energy accounting approach differs from that of fixed step iterative simulation codes, which use only a mean fluid temperature based on an assumed constant input to each component, because of the transient nature of the temperatures over the simulation time interval. The meteorological data are supplied in the form of constant time-averaged data, which may be averaged over time intervals of up to one hour maximum. The program assumes a constant value for each source term, e.g. H(t), for each time interval and does not attempt to link these values to form, for example, a ramp function as this process is unjustified when the raw data have been time-averaged. Two simulation example problems are presented which demonstrate the application of the improved code to the simulation of solar energy systems. In the first problem a domestic solar water heating system is considered. Here, the controller actions and the system temperature trajectories are demonstrated for a single day's operation. In the second problem, a more complex combined solar hot water and space heating system is considered where
I. INTRODUCTION
A number of simulation programs for predicting the long-term performance of different solar energy system designs are available. Of these, TRNSYS[1] is probably the most widely known program. This code is based on a modular approach with the system simulation performed by independently modelling the components whose governing equations are then solved in a serial fashion. The primary advantage of the modular approach is that it is flexible, and thus permits the rapid implementation of alternative system configurations. Since the governing equations are solved serially, simulation codes like TRNSYS require an iterative scheme in order to achieve an accuracy in the solution of the closed loop equations which approaches that of a direct simultaneous solution approach. The disadvantages in the use of an iterative scheme are two-fold: firstly, the governing equations have to be solved several times to ensure convergence of the result. Secondly, an indeterminate control status can occur when simulating components which have a hard nonlinear element such as an on/off type controller within the component, or whenever the controller is itself a component module which is interrogated during the iterative procedure. A simulation code is presented here which retains the desirable features of the modular approach, but with significant improvements over fixed step iterative simulation codes in the manner in which the system equations are solved and the controller actions are treated. A key feature of the code is that the simulation is advanced from "control function" to "control function" or to the end of the simulation time step, whichever is shorter. The program, although not yet commercially available, had been written in FORTRAN IV and has run successfully on the ICL 4-70, the CDC 7600 main frame computers as well as on the GEC 4070 minicomputer. The program requires approx. 50-60 k bytes of store on the CDC 7600 computer, although only 20 k bytes are actually used in a typical simulation run. 99
100
P. B. HOWELLSand R. H. MARSHALL
emphasis is placed on predicting the long-term system performance for different values of the user-specified error tolerance per time step. 2. SIMULATION PROGRAM STRUCTURE
All system components and subsystems, called modules, exist as FORTRAN subroutines. A list of the currently available modules is given in Appendix 1. The input data, aside from the meteorological data which are read from a magnetic tape or data file, are entered by means of assignment statements within the main program. The output appears on the line printer. The system parameters declared in the main program are transferred automatically from the main program to the component modules by means of two data arrays, each beginning with XIN and PAR. These arrays are unique to each module. The XIN array is used for the transfer of time dependent data, e.g. the mass flowrate through the collector, while the PAR array is used for the transfer of time independent data, e.g. the collector heat loss coefficient. A third array, beginning with OUT, is used to transfer the result of any internal calculations that may be made within the module to the main program, for example, the auxiliary energy rate required by an instantaneous water heater. Within the component modules, the coefficients of the governing state algebraic and differential equations are formulated automatically from the user's supplied input data. These coefficients are then transferred to a single operational module called a Central Processor, in which the new system temperatures which characterise the system state are evaluated. Thus, in the present approach the system temperatures do not appear at the output of the respective component modules as in many other simulation codes, but at the output of the Central Processor. The Central Processor incorporates the integration routine with access to the system control module with information feedback between the two. The integration routine advances the simulation through the solution of
® T 2 , 1 ~ T3,1 T/.1- -
System
the governing algebraic and differential equations, after which the control module checks for a change in control status. A fast matrix inversion routine as well as the automatic time stepping procedure are also included in the processing module to minimise computer time. 3. FORMULATION AND SOLUTION OF THE SYSTEM EQUATIONS
To demonstrate how the modular scheme is coupled with the direct simultaneous solution approach, consider the combined solar hot water and space heating system illustrated in Fig. 1. The constituent components and subsystems are shown as blocks in dotted lines with all capacitive elements shown as C.E. The advantage of coding subsystems as modules, i.e. those which represent a plurality of components, is that is simplifies the user's task in specifying the system configuration. The temperatures of the fluids within the system are labelled as T~.j, where i refers to both the number of the module as well as the order in which the module is called by the main program and j is the associated internal fluid temperature. To achieve a fast solution to the state matrix equation, all thermal storage devices are to appear last in the main program. As each module is called the, state matrix equation, {A}(J') + {B}(T) = (C)
is generated with all algebraic equations expressed in the form (T), = {D}(T)
~T ................
.......................
(2)
where (T) and (T)z are arrays containing the primary and secondary fluid temperatures respectively. The primary fluid temperatures are associated with the fluids which transfer energy across the component module boundaries, while the secondary fluid temperatures are associated with the fluids which distribute the energy within the module. Referring to Fig. 1, the primary fluid temperatures are T,,,, T2.1, T3,,, and T4., and the secondary fluid temperatures are T1.2T,,3T:.2 and T2,3. The {D}
[0 / t...............
(l)
o,01
I~L--I:-Y
I 7o~ver~ _-~_
I, i
~_. . . . . . .
Key (~)-Co&lecto~ Loop ~ - ~ o r Hot Water Heoting Loop
House (Fonco~l ,~ace Heoting System)
® Fig. 1. Combinedsolar hot water and space heating system.
Loop (~-C,~d Woter Tonk (Open Loop) ~)-Auxiioty Wott¢ Heoter IOpen Loopl (~-System Co~roCer
An improved computer code
for the simulationof solar heatingsystems
matrix, eqn (2), contains the functional relationships which exist between the primary and secondary fluid temperatures, or fluid temperatures internal and external to the module. As an example of casting the component equations in the form given by eqns (1) and (2), consider the combined collector/heat exchanger module illustrated in Fig. 2. This combined module has been written assuming that the solar collector operates in a quasi-transient mode, although the steady state version would have been sufficient. The first order linear differential equation describing the collector is given by[4]; = (tiC,,)r(T(t),.z - T(t),,,)
C, %
f A,F,z{(~a),rrH(t)-
(54
(54 Equations (3a) and (5a) have been combined to yield the following generalised differential equation; a,~ta,~+b,T,,tb,T,,,=c~
U,(T(t),,z-
F = WPh X(l-exp(w)).
T,(t))1 (3a)
(3b)
R
The heat exchanger is assumed to operate as a constant effectiveness device[S] with the heat transfer across the heat exchanger given by: Wh)
= (~C,MT(th,,
- T(thv)
= 5(tiC,U,,,(t)
- T,N)
(4)
where (tic,) is the minimum capacity fluid flowrate (i.e. the minimum of (ri~c,), and (tic,),). Equation (4) relates the fluid temperatures on either side of the heat exchanger which in {D} matrix form can be written as,
where D = ,_(~~E+WP)E II
(SIR
(G,)r
Ub)
D _ FG)E _ W,)E
I2 (+G)5 1-
(tic,),
(ha)
where
where
Qwrx = (IC,h(T(th.,-
101
(5c)
a, = CJL
(6b)
a2 = Cc&2
(6~)
b, = W,MDZI
- D,J + A~FRULD,,
bz = (tiC,M&-
Dn)+ AJGLJLDIz
cl = A,FR ((m)&(f)
+ Ur_Ta(t)).
(64 (be) (60
The above coefficient (i.e. a’s, b’s and cl) form elements within the state matrix equation, eqn (l), and are calculated automatically by the program. In the above example setting the effectiveness, 5, of the heat exchanger to unity reduces the combined collector/heat exchanger module to a simple flat plate collector module operating in the direct mode. Whenever a new module is required, the user needs to follow closely the above procedure to ensure that the governing equations are cast in a form compatible with the simulation code. Although the integration scheme can cope with modules in which some components are represented by algebraic equations, it is desirable to model a component subsystem with at least one capacitive element within each module, even if it is small in relation to the other capacitive elements in the system. The capacitance ensures that the state matrix equation associated with the integration scheme (discussed below) is diagonally dominant thus ensuring a stable solution.
-
-I
Fig. 2. Solar collector/heat
exchanger module.
102
P. B. HOWELLSand R. H. MARSHALL
After each module is called and the coefficients of all the algebraic and differential equations have been calculated, i.e. equations of the form (5b)-(5e) and (6b)-(6f), control is transferred automatically to the Central Processing module. At this point the state matrix equation (1), is converted into a form compatible with the RIKM integration code. Details of the procedure are reported by Howells and Marshall[2, 3]. Basically the procedure consists of (a) replacing the derivative terms (T) by the following finite difference approximation, T n+l- T"
7~ = ~
-~ ~o
The time-wise error bound estimate is given by: E = h ( t o l - ~9 oJ3+4to4-~1 ~o~)
(11)
and is used to generate automatic time step control. Each successive time interval is generated from the ratio of the estimated error per time step, eqn (11), and the user specified maximum error, E . . . . as follows: hn+l =
1t4hn"
Emax
(12)
(7a) 1
where n and n + 1 refer to the beginning and end of the time step of length h. (b) replacing each "T" on the r.h.s, of eqn (1) with weighted values of the variable at the n and n + 1 time levels as follows;
T=¢T"+I+(1-rr)T';
0~<~r~< 1.0.
(7b)
The transformed state matrix equation, eqn (1) then becomes {A'}(o~) = -{BI(T)" + (C)
(8)
which on performing the matrix multiplication on the r.h.s., reduces to, {A'}(~o) = (F).
(9)
Using the RIKM code, eqn (9) is solved in succession for the derivatives ~o~,o~2,..., ~o5 according to the equations: ~o = ~ (~o, + 4o~4 + ~o,) + 0(h) 5
0a,
1
The power of a has been selected rather than the (although the Merson scheme [2] claims an error per time step proportional to the time step raised to the fifth power, see eqn (10a) and not the fourth) to ensure that a conservative estimate of the new time step size is made. Using the derivative vector (~o), the inlet and outlet fluid temperatures of each module at the new time level are determined from the expression: (T) "+' = (T)" + (¢o) x h.
(13)
Finally, all algebraic equations are solved using eqn (2). Varying the weighting function (~r) and the number of function evaluations, i.e. ~o,,~o2. . . . . oJ,, for each time step permits access to the integration schemes listed in Table 1, where MERS is used to select the integration scheme for the time step. The modified state matrix {A'}, eqn (9), which is generated within the Central Processor from the state eqns (1) and (2) is, for the system illustrated in Fig. 1, of the form;
{A'}=
with
I
0°"41
a22 0 02.4 0
a3,3 a3.4
(14)
[. O,*.1 an.z 04.3 a4,4 ]
~o, = f(t, T)
(10b) (10c)
h T + h o,)
~3 = f ( t + ~h, T + ~ h (~o, + ¢o2))
and is sparse because all modules, other than the thermal module, are single input/output devices. To minimise computer time, use is made of the sparseness of the {A'} matrix when performing the matrix inversion. The governing equations are solved as follows:
(lOd)
[,, _ ~ , N-, f, aN,, i= 1
~oN = 2
h T+~(~o, h ~o`*= f(t +~,
+ 3w3))
aN.m -
(10e)
i=1
~o, = 1 ( f , - a, N~ON)' h o~5= [ (t + h, T + ~(~o~ - 3~3 + 4o~`*))
(lOf)
where f(t, T) is the r.h.s, of eqn (8). It can be seen from eqns (10a-10f) that, although the source terms such as irradiance and ambient temperature have been taken to be constant over the time interval, the integration scheme can accommodate time varying forcing functions, e.g. ramp inputs.
aid
"
di,i
(15a)
N - 1 aN.i a i . N
all
i=2...(N-
1).
(15b)
'
4. SIMULATIONOF CONTROLACTIONS At the end of the simulation time step control is transferred from the integration routine to the system control module, which in its present form simulates only simple on/off control with hysterisis. However, since the state variable derivative and integrated values are available, alternative control modules which simulate more
An improved computer code for the simulation of solar heating systems
103
Table 1. Range of methods treated with RIKM MERS
i.
2.
3.
4.
Integration Scheme
Implicit Runge Kutta M e r s o n (RIKM) (Fixed and v a r i a b l e step)
TL,3< T4,,;
0.5
N u m b e r of Functional Evaluations
to ~ w 1 5
Explicit Runge K u t t a M e r s o n (Fixed and v a r i a b l e step)
l
Crank Nicholson (fixed s t e p only)
0.5
First order E u l e r (fixed step only)
O
sophisticated forms of control, such as PDI control (proportional, derivative, integral), can be created. The on/off type control is simulated by means of two binary mode arrays, KP and KV, i.e. KP (i)= 1 or 0 and KV (i) = 1 or 0. These arrays are used to control the various pumps and valves throughout the system, respectively, by merely multiplying the pump capacity flowrate, or input mass flowrate in the case of a valve, by the appropriate binary operator. The reference temperature and dead bands associated with each control element are declared by the user in the main program. This information is then transferred to the control module via the corresponding XIN array. As an example of simulating an on/off control element, consider the control of the pumps P, and P2 within the collector/heat exchanger module shown in Fig. 1. It will be assumed that the control functions for these pumps are KP(1) and KP(2) respectively, and that control logic is
(T,,3+5°C)>-T4,,;
Weighting Function (o)
KP(1) = KP(2) = 1
5
(one step)
0o
1
(one step)
simulated, while each set of four of the remaining locations is used for the transfer of the associated reference temperatures and controller dead bands. Within the control module the program proceeds to execute the following logic statements for each control element. IF ((T1-T2). GE. DTON. and KP (I) = 0) KP (I)= 1 IF ((TI -T2). LT. DTOFF. and KP (I) = 1) KP (I) = 0 where I refers to the control element and T1 and T2 are the reference temperatures, which for the problem under consideration are T,,3 and T4: respectively. A similar procedure to the above is adopted for the control of the valves within the system. A key feature of the code is the detection of when a change in control status occurs through interpolating between the system states at the n and n + 1 time levels. For example, if in the above problem the pump switchon condition is satisfied, then the point in time when the switch occurred is given by
KP(1) = KP(2) = 0. n+h
8h = T~'Th -(T4., For this particular problem, the user declares in the main program the following assignment statements. Source code XINCMP (1) = 2.0 XINCMP (2)= 13 XINCMP (3)-- 41 XINCMP (4) = 5.0 XINCMP (5) ---0.0 XINCMP (6)= 13 XINCMP (7) = 41 XINCMP (8) = 5.0 XINCMP (9) = 0.0
Interpretation 2 control elements to be simulated Tl,3 T4,1 1st Controller DTON DTOFF Tt.3 2nd controller
T4: DTON DTOFF
As shown above, the first array location, XINCMP (1), is used to flag the number of on/off controllers to be
+DTON)
to3,1-- ~4,1
(16)
where 8h is the time interval measured back from the n + h time level. The corresponding system state at this point is approximated by: (T) "+"-8" = (T) "÷" - (to) x ~h.
(17)
Computation will now proceed from the point in time when the change occurred to the end of the time interval h, using as the new time step ~h. When several such changes in control status occur, arising from multiple control elements, the program selects the larger of the ~h's, that is the solution advances to the first occurrence of a change in control.
104
P.B. HOWELLSand R. H. MARSHALL 5. SIMULATION EXAMPLE PROBLEMS
GO.Om
(i) Example 1 A particularly severe test for any simulation program is an accurate representation of on/off controller actions. These devices represent hard nonlinear elements which, for simulation codes that use a fixed time step size, are allowed to change status only at an interger number of the simulation time step. The problem becomes even more severe when the time interval between switches is very short. Therefore, in this particular example attention is focussed on the technique of advancing the simulation from control function to control function in order to simulate such devices as on/off controllers independently of the time interval between switches. In this first example, for comparison between techniques, the temperature trajectories for one day's system operation are predicted using the simultaneous closed loop solution approach as well as a serial solution approach. In the latter, the temperature of the fluid at the output of the collector was determined under conditions of continuous (i.e. "imaginary") flow. Figure 3 illustrates a typical solar water heating system incorporating two on/off control elements: the on/off controller which controls the fluid flow in the collector loop and the on/off controller which controls the fluid temperature at the hot water tank. The control logic for both these controllers its assumed to be;
50.0
z.O.O
30.0 E 20.0
10.0" I
o ~ L
8
10
12
1/,.
16
18 20
22
2l-
Time (hr)
61.0 60-0
Off
Off
On
On
_
Off Off _ _
Upoer Temp
On On On
Lower Temp Limit
56.0 57.04 ~, 56.01 -~ 55.0 Io 54 O[ 2
z,
6
10
13
12
14
16
t8
20
22
24
Time (hr)
Fig. 4. Hot water demand pattern and hot water tank fluid temperature trajectory.
T~,1 >- T2,1 o f f - - ) o n ] ~ ,, T,,~ < T2.j on~ off.~ t~ouector loop 1
periods as short as rahr between switches. To achieve the same order of accuracy using a fixed step the simu1 lation would have been constrained to steps of ~ hr or less in order to detect the immersion heater controller actions, and thus avoid errors in the energy balances. Hence, by advancing from one control function to the next, the improved computer code is able to adjust the time step size to cope with rapid changes in control status. For the domestic solar water heating system illustrated in Fig. 3, the preheat tank and collector fluid temperature
T . ~<55°C off ~ on] TH >/60°C o n ~ offJ Hot water tank and is incorporated within the control module. Further simulation data for the system illustrated in Fig. 3 are given in Table A2.1, Appendix 2. The temperature trajectory at the hot water tank for one day's operation is presented in Fig. 4. The time period of interest lies between 18:00 and 24:00hr when as a result of rapid changes in demand (also shown in Fig. 4) there exist
Collector/Preheat Tank Pump Controller
iTc~d
I ..... 7
Demond
TDemand
~
I
J I
Co
J
i
,I,TZI
T2.1'
Preheat Tank
IH .....
I I ....
J
Hot Water Tank
I TCotd
Fig. 3. Domesticsolar water heating system.
J "1
J Auxiliary I H~ter
J J Controller
105
An improved computer code for the simulationof solar heating systems
hour meteorological data for different values of the userspecified maximum error tolerance per time step, Emax (see eqn (12)). At each time step, the energy transfers within the system are accounted for and used to determine the solar contribution to the domestic hot water and space heating loads. The accuracy of the simulation is judged by the degree to within which the energy balances agree, called "closure" by Klein et a/.[1]. For the system presented in Fig. 1, the largest energy imbalance was found to occur at the thermal store. The degree of imbalance, as a function of the user maximum error tolerance per time step (Emax), is presented in Fig. 6, together with the required computer processor times for the CDC 7600 computer. The results indicate that, for the system under consideration, the optimum error tolerance is approximately 0.01°C, and yields a yearly energy balance correct to within 1% while requiring only 9sec of computer time. The predicted system performance for the 0.01°C error tolerance is presented in Table 2. This table shows that 58 and 27 percent of the domestic hot water and space heating loads respectively are supplied by solar energy. Thus, as demonstrated in this problem, the improved computer code can be used for long-term system performance studies as well as detailed system control studies.
trajectories are seen in Fig. 5 for a single day of a yearly simulation using hourly meteorological data. These trajectories have been computed using the improved code as well as using a code based on a serial solution approach (i.e. one in which each system component is assumed to be decoupled from the rest of the system subject to a constant input). It is seen that both programs yield the same trajectories, but that the improved code can use a time step four times that permitted by the program based on the serial solution approach. The result is a 43 percent saving in computer time for an entire year's simulation on the CDC 7600 computer. Thus, the improved code may be used in situations where other fixed step iterative programs may prove to be prohibitively expensive. In the above example problem the use of a 1 hr time step using the program based on the serial solution approach was attemped. However, gross inaccuracies in both the computed temperature trajectories and the final energy balances were observed as a result of the large number of immersion heater switch-on's and errors propagating around the system as a result of assuming a constant input to each module. An alternative approach to reducing the time step from 1 to 41hr would have been to employ an iterative procedure at the end of each time step, as used in most simulation codes where the input to the module is assumed constant throughout the ensuing time step. However, in either case the resulting computer operating time would have been greater than that required by the improved simulation code.
6. CONCLUSION
An improved computer code which represents a significant improvement over existing fixed step iterative simulation programs has been presented. The key features of the program are the simultaneous solution of all differential equations appearing in closed loop form, and the technique for advancing the simulation form "control function" to "control function". Two simulation example problems demonstrating the application of the improved code to studies in control and long-term system performance were also presented.
(ii) Example 2 In the next problem the long-term performance of the combined solar hot water and space heating system illustrated in Fig. 1. is determined. The simulation data are given in Table A2.2, Appendix 2. The system is simulated for a single year of operation using hour by Key
60.0-
o-Improved Coml:xJter Code (h=lhr) =-Fixed Step Computer Code (h=l/4hr) (Imoginory Flow Conditon}
Weother Doto KEW (1963-1964) lZ.th Moy Pump
Off
50.0-
40.0-
Pump
On Z
~ 30-0~ I
,• /•"
"
Preheot Tonk Temp Trajectory
o - - -o
-- - o -
-- o - - - o - "
\ ,, \
• X \ \
,)
k
? ..,.."-..---Collector FluidOutlet
20.0-
"\
Temp.
o,
% ~o ~'o
Trajectory
o
o
10.0-
Time (hrs)
1i 1'~ 15 16 17 lb 19 20 2i
Fig. 5. Collector and preheat tank fluid temperature trajectories.
2'2 2'3 2t.
106
P.B. HOWELLSand R. H. MARSHALL 30.0-
6"L°
28.0' 26.0-
/
5.0
2A.O.
Key
,/
Computer "Dme {CDC 7600 Computer)
22.020.018.016.0-
,i,1'1!
,, <.o-!
% Imbolonce
- - - -
14.012.0IO.OE o
8.0.
~ 2.0-
6.0-
?
4.0-
°c w
/ 1.0
/
2.01.0
I
0,0001
0,001
001
0.1
Maximum Error Tolerance Per Time Step (£Mox)
Fig. 6. Accuracy and computer time as a function of the user-specified error tolerance (Emax).
Table 2. Combined solar hot water and space heating system performance data (Emax= 0.01°C)
Collector input
177.3
Collector output
GJ
32.87 GJ
Thermal store input
32.87 GJ
Thermal store output
28.6
GJ
Thermal store loss
4.3
GJ
Preheat tank input
8.4
GJ
Preheat tank output
8.0
GJ
Preheat tank loss
0.4
GJ
Auxiliary energy required by hot water system
5.4
GJ
Hot water load
13.7
GJ
Space heating heat exchanger input
20.2
GJ
Auxiliary energy required by space heating system.
53.6
GJ
Space heating load
73.9
GJ
% Hot water load supplied by solar energy
58%
% Space heating load supplied by solar energy
27%
In the first problem, a domestic solar water heating system was considered in which the temperature of the fluid at the hot water tank is controlled by an immersion heater and an on/off controller. The rapid switching
action of the immersion heater controller posed a very severe test. However, by advancing from each control action to the next (or to the end of the hour) the improved code resulted in an accurate simulation of the
An improved computer code for the simulation of solar heating systems controller actions without constraining the nominal time step to be less than the shortest time interval between 1 switches, calculated to be io hr. For the same problem, the system temperature trajectories were also computed and compared with those of a fixed step simulation program in which the input to each module was assumed to be constant during the ensuing time step. This comparison showed that the trajectories agree well even at a nominal 1 hr time step for the improved code, compared with the "hr time step used by the fixed step program. As a result of the larger nominal time step size a 43 percent saving in computer time is realised for the improved code. This saving is a direct result of the simultaneous solution approach which allows for time varying input to each module, thus avoiding any reduction in the simulation time step or the need for any form of iterative procedure to ensure the accuracy of the result. In the second problem a more complex combined solar hot water and space heating system was considered. The system was simulated for a single year's operation using h o u / b y hour meteorological data for different values of the user-specified error tolerance, E . . . . Based on a comparison between computer time and the largest energy transfer, which occurred at the thermal store, the optimum error tolerance was found to be 0.01°C. This value yielded a yearly energy balance correct to within 1 percent while requiring only 9 sec of computer operating time for the year's simulation. These example problems demonstrated that the improved code can be used for detailed system control as well as long-term system performance studies. Furthermore, because the system equations are solved simultaneously marching from "control function" to "control function" a high degree of accuracy for a reduced computer time is realised. NOMENCLATURE Ac collector area, m2
{A} matrix elements C~, specific heat capacity, kJ kg ~K {D} matrix elements Oi.i E error bound estimate Em~x user-specified maximum error per time step f~ {F} matrix elements F' plate efficiency factor h time step size, sec HO) irradiance on the tilted surface, Wm 2 MERS integration method selection parameter mass flowrate, kg sec mCp minimum capacity fluid flowrate, WK N matrix dimension Cc collector heat capacity, JK (T) fluid temperatures external to the component modules °C (T)r fluid temperatures internal to the component modules °C T,(t) ambient air temperature, °C Tm inlet fluid temperature, °C Tout outlet fluid temperature, °C rt fluid temperature at the preheat tank (domestic solar water heating system), °C T2 fluid temperatures at the output of the collector (domestic solar water heating system), °C r . fluid temperature at the hot water tank, °C rsH space heating demand temperature, °C To domestic hot water demand temperature, °C ai.i
107
TRM room temperature, °C
UL overall collector heat loss coefficient, WK ~ m 2 U heat loss coefficient UABLDC; building heat loss coefficient/area product, WK-t V(n volume, me Greek ol solar absorptance 8h time when a change in control function occurs measured back from the end of the time step "h", sec heat exchanger effectiveness o angle between the beam radiation and the normal to the collector t7 weighting function "7d diffuse radiation transmittance coefficient beam radiation transmittance coefficient "rb (Ol --->oj ~ RIKM integration coefficients
Subscripts I internal loop parameters E external loop parameters CL collector loop DL domestic hot water loop SL space heating loop PT preheat tank TS thermal store HT hot water tank Superscripts n , n + l time level Other { } square matrix ( ) column matrix rate of change with time
REFERENCES
1. TRNSYS, TRNSYS Simulation Manual. Solar Energy Laboratory, University Wisconsin-Madison (1979). 2. R. H. Marshall, A. E. Churches and J. A. Reizes, A hybrid integration scheme for the solution of viscous compressible flows. Appl. Math. Model. 3, 459--465. (1979). 3. P. B. Howells and R. H. Marshall, A fast and accurate integration scheme for coupled stiff differential equations, pp. 45-59. Proc. Conf. Numerical Methods for Coupled Problems. Pineridge Press, London (1981). 4. J. A. Duffle and W. A. Beckman, Solar Energy Thermal Processes. Wiley, New York (1974). 5. W. M. Kays and A. L. London, Compact Heat Exchangers. McGraw-Hill, New York (1958).
APPENDIX 1
Table AI.1. Component modules currently available (a) (b) (c) (d) (e) (f) (g) (h) (i) (j) (k) (I) (m) (n)
Solar radiation processing module Pumpmodule Cold water tank module Combinedflat plate collector/heat exchanger module Combinedpreheat tank/heat exchanger module Hot water tank module Fan coil space heating module Water radiator space heating system module Fully mixed thermal store module Stratifiedstore module Phasechange thermal store module Clock module Energy meter module On/Offcontrol module
P.B. HOWELLSand R. H. MARSHALL
108 APPENDIX 2
Table A2.2. Combined solar hot water and space heating system simulation data
Table A2.1. Domestic solar water heating system simulation data Parameter Ac ~ICL
F'
C 7d 7b Ot
G VOLer
Uer VOLHr
Um Integration scheme
Emax Auxiliary water heater rating Weather data
Parameter
Value
Ac
50 m 2 6.0WK ~m 2 0.017 kg sec i m-2 0.9 10.0 kJ K -t m -2 0.76 0.84 (-0.0119 + 3.42 cos 0 -3.97 cos: 0 + 1.56 cos 3 0) 0.9 4.2kJkg 1K i 1.0 2001 0.4Wm 2K h 0.068 kg sec 60°C 751m 2 0.4Wm 2K 20°C 20°C 400 WK t 0.8 0.017 kg sec -1 m -2 1.007 kJ kg K t 0.008 kg sec -t m -2 RIKM Kew (England) 1963-64
Value 4.0m 2 6.0 WK-t m 2 0.017 kg sec ~m 2 0.9 10kJ K-I m 2 0.76 0.84 (-0.0119+ 3.142 COS 0 -3.97 COS2 0 + 1.56 COS3 0) 0.9 4.2 kJ kg I K I 2001 0.4Wm 2K 1301 0.4 Wm -2 K-i RIKM 0.01 K 3 KW Kew (England) 1963-64.
Uc ~hcc F' C~ Td "rb Ol
G ~.cc
VOLpT UpT
]]//DL TD VOLts
Urs
TsH TRy UABLtm ~.SL
tilsg
Cpair r/lair Integration scheme Weather data
0038~92X/831020099-10503.00t0 © 1983 Pergamon Press Ltd.
AN IMPROVED COMPUTER CODE FOR THE SIMULATION OF SOLAR HEATING SYSTEMS P. B. HOWELLSand R. H. MARSHALL Solar Energy Unit, University College, Cardiff CF2 1TA, Wales
(Received 12 March 1982; accepted 16 April 1982) Abstraet--A generalised solar energy simulation code, which represents a significantimprovement over fixed step iterative simulation codes, is presented. The key features of the code are the simultaneous solution of all differential equations arising from closed loops in the system, and the technique for advancing the simulation from "control function" to "control function". Computer time is minimised by the use of sparse matrix techniques together with automatic time step control generated from an error bound estimate. The application of the improved code to the simulation of solar energy systems is demonstrated on two example problems. In the first problem, emphasis is placed on predicting the temperatures of the fluids within the system and the controller actions, while in the second problem the emphasis is on predicting long-term system performance. The simulation results of both these problems demonstrate that the improved code is suitable for both detailed system control and long-term system performance studies.
The governing system differential equations are formulated as sets of coupled differential equations to form a state matrix which is solved directly (i.e. without iteration) using a highly stable integration scheme codenamed RIKM[2, 3], with provision for alternative explicit and implicit integration schemes. Soft non-linearities, such as temperature dependent material properties, are treated by computing the state matrix equation at each time step. Automatic time step control generated from an error bound estimate is used to maintain the accuracy of the simulation to within a user-specified error tolerance. The choice of the error tolerance, as will be demonstrated, is usually a compromise between computer processor time and the closure of the final energy balances. The energies are determined within the program by assuming a constant specific heat capacity and then numerically integrating the computed temperature time trajectories using the trapezoidal rule. This energy accounting approach differs from that of fixed step iterative simulation codes, which use only a mean fluid temperature based on an assumed constant input to each component, because of the transient nature of the temperatures over the simulation time interval. The meteorological data are supplied in the form of constant time-averaged data, which may be averaged over time intervals of up to one hour maximum. The program assumes a constant value for each source term, e.g. H(t), for each time interval and does not attempt to link these values to form, for example, a ramp function as this process is unjustified when the raw data have been time-averaged. Two simulation example problems are presented which demonstrate the application of the improved code to the simulation of solar energy systems. In the first problem a domestic solar water heating system is considered. Here, the controller actions and the system temperature trajectories are demonstrated for a single day's operation. In the second problem, a more complex combined solar hot water and space heating system is considered where
I. INTRODUCTION
A number of simulation programs for predicting the long-term performance of different solar energy system designs are available. Of these, TRNSYS[1] is probably the most widely known program. This code is based on a modular approach with the system simulation performed by independently modelling the components whose governing equations are then solved in a serial fashion. The primary advantage of the modular approach is that it is flexible, and thus permits the rapid implementation of alternative system configurations. Since the governing equations are solved serially, simulation codes like TRNSYS require an iterative scheme in order to achieve an accuracy in the solution of the closed loop equations which approaches that of a direct simultaneous solution approach. The disadvantages in the use of an iterative scheme are two-fold: firstly, the governing equations have to be solved several times to ensure convergence of the result. Secondly, an indeterminate control status can occur when simulating components which have a hard nonlinear element such as an on/off type controller within the component, or whenever the controller is itself a component module which is interrogated during the iterative procedure. A simulation code is presented here which retains the desirable features of the modular approach, but with significant improvements over fixed step iterative simulation codes in the manner in which the system equations are solved and the controller actions are treated. A key feature of the code is that the simulation is advanced from "control function" to "control function" or to the end of the simulation time step, whichever is shorter. The program, although not yet commercially available, had been written in FORTRAN IV and has run successfully on the ICL 4-70, the CDC 7600 main frame computers as well as on the GEC 4070 minicomputer. The program requires approx. 50-60 k bytes of store on the CDC 7600 computer, although only 20 k bytes are actually used in a typical simulation run. 99
100
P. B. HOWELLSand R. H. MARSHALL
emphasis is placed on predicting the long-term system performance for different values of the user-specified error tolerance per time step. 2. SIMULATION PROGRAM STRUCTURE
All system components and subsystems, called modules, exist as FORTRAN subroutines. A list of the currently available modules is given in Appendix 1. The input data, aside from the meteorological data which are read from a magnetic tape or data file, are entered by means of assignment statements within the main program. The output appears on the line printer. The system parameters declared in the main program are transferred automatically from the main program to the component modules by means of two data arrays, each beginning with XIN and PAR. These arrays are unique to each module. The XIN array is used for the transfer of time dependent data, e.g. the mass flowrate through the collector, while the PAR array is used for the transfer of time independent data, e.g. the collector heat loss coefficient. A third array, beginning with OUT, is used to transfer the result of any internal calculations that may be made within the module to the main program, for example, the auxiliary energy rate required by an instantaneous water heater. Within the component modules, the coefficients of the governing state algebraic and differential equations are formulated automatically from the user's supplied input data. These coefficients are then transferred to a single operational module called a Central Processor, in which the new system temperatures which characterise the system state are evaluated. Thus, in the present approach the system temperatures do not appear at the output of the respective component modules as in many other simulation codes, but at the output of the Central Processor. The Central Processor incorporates the integration routine with access to the system control module with information feedback between the two. The integration routine advances the simulation through the solution of
® T 2 , 1 ~ T3,1 T/.1- -
System
the governing algebraic and differential equations, after which the control module checks for a change in control status. A fast matrix inversion routine as well as the automatic time stepping procedure are also included in the processing module to minimise computer time. 3. FORMULATION AND SOLUTION OF THE SYSTEM EQUATIONS
To demonstrate how the modular scheme is coupled with the direct simultaneous solution approach, consider the combined solar hot water and space heating system illustrated in Fig. 1. The constituent components and subsystems are shown as blocks in dotted lines with all capacitive elements shown as C.E. The advantage of coding subsystems as modules, i.e. those which represent a plurality of components, is that is simplifies the user's task in specifying the system configuration. The temperatures of the fluids within the system are labelled as T~.j, where i refers to both the number of the module as well as the order in which the module is called by the main program and j is the associated internal fluid temperature. To achieve a fast solution to the state matrix equation, all thermal storage devices are to appear last in the main program. As each module is called the, state matrix equation, {A}(J') + {B}(T) = (C)
is generated with all algebraic equations expressed in the form (T), = {D}(T)
~T ................
.......................
(2)
where (T) and (T)z are arrays containing the primary and secondary fluid temperatures respectively. The primary fluid temperatures are associated with the fluids which transfer energy across the component module boundaries, while the secondary fluid temperatures are associated with the fluids which distribute the energy within the module. Referring to Fig. 1, the primary fluid temperatures are T,,,, T2.1, T3,,, and T4., and the secondary fluid temperatures are T1.2T,,3T:.2 and T2,3. The {D}
[0 / t...............
(l)
o,01
I~L--I:-Y
I 7o~ver~ _-~_
I, i
~_. . . . . . .
Key (~)-Co&lecto~ Loop ~ - ~ o r Hot Water Heoting Loop
House (Fonco~l ,~ace Heoting System)
® Fig. 1. Combinedsolar hot water and space heating system.
Loop (~-C,~d Woter Tonk (Open Loop) ~)-Auxiioty Wott¢ Heoter IOpen Loopl (~-System Co~roCer
An improved computer code
for the simulationof solar heatingsystems
matrix, eqn (2), contains the functional relationships which exist between the primary and secondary fluid temperatures, or fluid temperatures internal and external to the module. As an example of casting the component equations in the form given by eqns (1) and (2), consider the combined collector/heat exchanger module illustrated in Fig. 2. This combined module has been written assuming that the solar collector operates in a quasi-transient mode, although the steady state version would have been sufficient. The first order linear differential equation describing the collector is given by[4]; = (tiC,,)r(T(t),.z - T(t),,,)
C, %
f A,F,z{(~a),rrH(t)-
(54
(54 Equations (3a) and (5a) have been combined to yield the following generalised differential equation; a,~ta,~+b,T,,tb,T,,,=c~
U,(T(t),,z-
F = WPh X(l-exp(w)).
T,(t))1 (3a)
(3b)
R
The heat exchanger is assumed to operate as a constant effectiveness device[S] with the heat transfer across the heat exchanger given by: Wh)
= (~C,MT(th,,
- T(thv)
= 5(tiC,U,,,(t)
- T,N)
(4)
where (tic,) is the minimum capacity fluid flowrate (i.e. the minimum of (ri~c,), and (tic,),). Equation (4) relates the fluid temperatures on either side of the heat exchanger which in {D} matrix form can be written as,
where D = ,_(~~E+WP)E II
(SIR
(G,)r
Ub)
D _ FG)E _ W,)E
I2 (+G)5 1-
(tic,),
(ha)
where
where
Qwrx = (IC,h(T(th.,-
101
(5c)
a, = CJL
(6b)
a2 = Cc&2
(6~)
b, = W,MDZI
- D,J + A~FRULD,,
bz = (tiC,M&-
Dn)+ AJGLJLDIz
cl = A,FR ((m)&(f)
+ Ur_Ta(t)).
(64 (be) (60
The above coefficient (i.e. a’s, b’s and cl) form elements within the state matrix equation, eqn (l), and are calculated automatically by the program. In the above example setting the effectiveness, 5, of the heat exchanger to unity reduces the combined collector/heat exchanger module to a simple flat plate collector module operating in the direct mode. Whenever a new module is required, the user needs to follow closely the above procedure to ensure that the governing equations are cast in a form compatible with the simulation code. Although the integration scheme can cope with modules in which some components are represented by algebraic equations, it is desirable to model a component subsystem with at least one capacitive element within each module, even if it is small in relation to the other capacitive elements in the system. The capacitance ensures that the state matrix equation associated with the integration scheme (discussed below) is diagonally dominant thus ensuring a stable solution.
-
-I
Fig. 2. Solar collector/heat
exchanger module.
102
P. B. HOWELLSand R. H. MARSHALL
After each module is called and the coefficients of all the algebraic and differential equations have been calculated, i.e. equations of the form (5b)-(5e) and (6b)-(6f), control is transferred automatically to the Central Processing module. At this point the state matrix equation (1), is converted into a form compatible with the RIKM integration code. Details of the procedure are reported by Howells and Marshall[2, 3]. Basically the procedure consists of (a) replacing the derivative terms (T) by the following finite difference approximation, T n+l- T"
7~ = ~
-~ ~o
The time-wise error bound estimate is given by: E = h ( t o l - ~9 oJ3+4to4-~1 ~o~)
(11)
and is used to generate automatic time step control. Each successive time interval is generated from the ratio of the estimated error per time step, eqn (11), and the user specified maximum error, E . . . . as follows: hn+l =
1t4hn"
Emax
(12)
(7a) 1
where n and n + 1 refer to the beginning and end of the time step of length h. (b) replacing each "T" on the r.h.s, of eqn (1) with weighted values of the variable at the n and n + 1 time levels as follows;
T=¢T"+I+(1-rr)T';
0~<~r~< 1.0.
(7b)
The transformed state matrix equation, eqn (1) then becomes {A'}(o~) = -{BI(T)" + (C)
(8)
which on performing the matrix multiplication on the r.h.s., reduces to, {A'}(~o) = (F).
(9)
Using the RIKM code, eqn (9) is solved in succession for the derivatives ~o~,o~2,..., ~o5 according to the equations: ~o = ~ (~o, + 4o~4 + ~o,) + 0(h) 5
0a,
1
The power of a has been selected rather than the (although the Merson scheme [2] claims an error per time step proportional to the time step raised to the fifth power, see eqn (10a) and not the fourth) to ensure that a conservative estimate of the new time step size is made. Using the derivative vector (~o), the inlet and outlet fluid temperatures of each module at the new time level are determined from the expression: (T) "+' = (T)" + (¢o) x h.
(13)
Finally, all algebraic equations are solved using eqn (2). Varying the weighting function (~r) and the number of function evaluations, i.e. ~o,,~o2. . . . . oJ,, for each time step permits access to the integration schemes listed in Table 1, where MERS is used to select the integration scheme for the time step. The modified state matrix {A'}, eqn (9), which is generated within the Central Processor from the state eqns (1) and (2) is, for the system illustrated in Fig. 1, of the form;
{A'}=
with
I
0°"41
a22 0 02.4 0
a3,3 a3.4
(14)
[. O,*.1 an.z 04.3 a4,4 ]
~o, = f(t, T)
(10b) (10c)
h T + h o,)
~3 = f ( t + ~h, T + ~ h (~o, + ¢o2))
and is sparse because all modules, other than the thermal module, are single input/output devices. To minimise computer time, use is made of the sparseness of the {A'} matrix when performing the matrix inversion. The governing equations are solved as follows:
(lOd)
[,, _ ~ , N-, f, aN,, i= 1
~oN = 2
h T+~(~o, h ~o`*= f(t +~,
+ 3w3))
aN.m -
(10e)
i=1
~o, = 1 ( f , - a, N~ON)' h o~5= [ (t + h, T + ~(~o~ - 3~3 + 4o~`*))
(lOf)
where f(t, T) is the r.h.s, of eqn (8). It can be seen from eqns (10a-10f) that, although the source terms such as irradiance and ambient temperature have been taken to be constant over the time interval, the integration scheme can accommodate time varying forcing functions, e.g. ramp inputs.
aid
"
di,i
(15a)
N - 1 aN.i a i . N
all
i=2...(N-
1).
(15b)
'
4. SIMULATIONOF CONTROLACTIONS At the end of the simulation time step control is transferred from the integration routine to the system control module, which in its present form simulates only simple on/off control with hysterisis. However, since the state variable derivative and integrated values are available, alternative control modules which simulate more
An improved computer code for the simulation of solar heating systems
103
Table 1. Range of methods treated with RIKM MERS
i.
2.
3.
4.
Integration Scheme
Implicit Runge Kutta M e r s o n (RIKM) (Fixed and v a r i a b l e step)
TL,3< T4,,;
0.5
N u m b e r of Functional Evaluations
to ~ w 1 5
Explicit Runge K u t t a M e r s o n (Fixed and v a r i a b l e step)
l
Crank Nicholson (fixed s t e p only)
0.5
First order E u l e r (fixed step only)
O
sophisticated forms of control, such as PDI control (proportional, derivative, integral), can be created. The on/off type control is simulated by means of two binary mode arrays, KP and KV, i.e. KP (i)= 1 or 0 and KV (i) = 1 or 0. These arrays are used to control the various pumps and valves throughout the system, respectively, by merely multiplying the pump capacity flowrate, or input mass flowrate in the case of a valve, by the appropriate binary operator. The reference temperature and dead bands associated with each control element are declared by the user in the main program. This information is then transferred to the control module via the corresponding XIN array. As an example of simulating an on/off control element, consider the control of the pumps P, and P2 within the collector/heat exchanger module shown in Fig. 1. It will be assumed that the control functions for these pumps are KP(1) and KP(2) respectively, and that control logic is
(T,,3+5°C)>-T4,,;
Weighting Function (o)
KP(1) = KP(2) = 1
5
(one step)
0o
1
(one step)
simulated, while each set of four of the remaining locations is used for the transfer of the associated reference temperatures and controller dead bands. Within the control module the program proceeds to execute the following logic statements for each control element. IF ((T1-T2). GE. DTON. and KP (I) = 0) KP (I)= 1 IF ((TI -T2). LT. DTOFF. and KP (I) = 1) KP (I) = 0 where I refers to the control element and T1 and T2 are the reference temperatures, which for the problem under consideration are T,,3 and T4: respectively. A similar procedure to the above is adopted for the control of the valves within the system. A key feature of the code is the detection of when a change in control status occurs through interpolating between the system states at the n and n + 1 time levels. For example, if in the above problem the pump switchon condition is satisfied, then the point in time when the switch occurred is given by
KP(1) = KP(2) = 0. n+h
8h = T~'Th -(T4., For this particular problem, the user declares in the main program the following assignment statements. Source code XINCMP (1) = 2.0 XINCMP (2)= 13 XINCMP (3)-- 41 XINCMP (4) = 5.0 XINCMP (5) ---0.0 XINCMP (6)= 13 XINCMP (7) = 41 XINCMP (8) = 5.0 XINCMP (9) = 0.0
Interpretation 2 control elements to be simulated Tl,3 T4,1 1st Controller DTON DTOFF Tt.3 2nd controller
T4: DTON DTOFF
As shown above, the first array location, XINCMP (1), is used to flag the number of on/off controllers to be
+DTON)
to3,1-- ~4,1
(16)
where 8h is the time interval measured back from the n + h time level. The corresponding system state at this point is approximated by: (T) "+"-8" = (T) "÷" - (to) x ~h.
(17)
Computation will now proceed from the point in time when the change occurred to the end of the time interval h, using as the new time step ~h. When several such changes in control status occur, arising from multiple control elements, the program selects the larger of the ~h's, that is the solution advances to the first occurrence of a change in control.
104
P.B. HOWELLSand R. H. MARSHALL 5. SIMULATION EXAMPLE PROBLEMS
GO.Om
(i) Example 1 A particularly severe test for any simulation program is an accurate representation of on/off controller actions. These devices represent hard nonlinear elements which, for simulation codes that use a fixed time step size, are allowed to change status only at an interger number of the simulation time step. The problem becomes even more severe when the time interval between switches is very short. Therefore, in this particular example attention is focussed on the technique of advancing the simulation from control function to control function in order to simulate such devices as on/off controllers independently of the time interval between switches. In this first example, for comparison between techniques, the temperature trajectories for one day's system operation are predicted using the simultaneous closed loop solution approach as well as a serial solution approach. In the latter, the temperature of the fluid at the output of the collector was determined under conditions of continuous (i.e. "imaginary") flow. Figure 3 illustrates a typical solar water heating system incorporating two on/off control elements: the on/off controller which controls the fluid flow in the collector loop and the on/off controller which controls the fluid temperature at the hot water tank. The control logic for both these controllers its assumed to be;
50.0
z.O.O
30.0 E 20.0
10.0" I
o ~ L
8
10
12
1/,.
16
18 20
22
2l-
Time (hr)
61.0 60-0
Off
Off
On
On
_
Off Off _ _
Upoer Temp
On On On
Lower Temp Limit
56.0 57.04 ~, 56.01 -~ 55.0 Io 54 O[ 2
z,
6
10
13
12
14
16
t8
20
22
24
Time (hr)
Fig. 4. Hot water demand pattern and hot water tank fluid temperature trajectory.
T~,1 >- T2,1 o f f - - ) o n ] ~ ,, T,,~ < T2.j on~ off.~ t~ouector loop 1
periods as short as rahr between switches. To achieve the same order of accuracy using a fixed step the simu1 lation would have been constrained to steps of ~ hr or less in order to detect the immersion heater controller actions, and thus avoid errors in the energy balances. Hence, by advancing from one control function to the next, the improved computer code is able to adjust the time step size to cope with rapid changes in control status. For the domestic solar water heating system illustrated in Fig. 3, the preheat tank and collector fluid temperature
T . ~<55°C off ~ on] TH >/60°C o n ~ offJ Hot water tank and is incorporated within the control module. Further simulation data for the system illustrated in Fig. 3 are given in Table A2.1, Appendix 2. The temperature trajectory at the hot water tank for one day's operation is presented in Fig. 4. The time period of interest lies between 18:00 and 24:00hr when as a result of rapid changes in demand (also shown in Fig. 4) there exist
Collector/Preheat Tank Pump Controller
iTc~d
I ..... 7
Demond
TDemand
~
I
J I
Co
J
i
,I,TZI
T2.1'
Preheat Tank
IH .....
I I ....
J
Hot Water Tank
I TCotd
Fig. 3. Domesticsolar water heating system.
J "1
J Auxiliary I H~ter
J J Controller
105
An improved computer code for the simulationof solar heating systems
hour meteorological data for different values of the userspecified maximum error tolerance per time step, Emax (see eqn (12)). At each time step, the energy transfers within the system are accounted for and used to determine the solar contribution to the domestic hot water and space heating loads. The accuracy of the simulation is judged by the degree to within which the energy balances agree, called "closure" by Klein et a/.[1]. For the system presented in Fig. 1, the largest energy imbalance was found to occur at the thermal store. The degree of imbalance, as a function of the user maximum error tolerance per time step (Emax), is presented in Fig. 6, together with the required computer processor times for the CDC 7600 computer. The results indicate that, for the system under consideration, the optimum error tolerance is approximately 0.01°C, and yields a yearly energy balance correct to within 1% while requiring only 9sec of computer time. The predicted system performance for the 0.01°C error tolerance is presented in Table 2. This table shows that 58 and 27 percent of the domestic hot water and space heating loads respectively are supplied by solar energy. Thus, as demonstrated in this problem, the improved computer code can be used for long-term system performance studies as well as detailed system control studies.
trajectories are seen in Fig. 5 for a single day of a yearly simulation using hourly meteorological data. These trajectories have been computed using the improved code as well as using a code based on a serial solution approach (i.e. one in which each system component is assumed to be decoupled from the rest of the system subject to a constant input). It is seen that both programs yield the same trajectories, but that the improved code can use a time step four times that permitted by the program based on the serial solution approach. The result is a 43 percent saving in computer time for an entire year's simulation on the CDC 7600 computer. Thus, the improved code may be used in situations where other fixed step iterative programs may prove to be prohibitively expensive. In the above example problem the use of a 1 hr time step using the program based on the serial solution approach was attemped. However, gross inaccuracies in both the computed temperature trajectories and the final energy balances were observed as a result of the large number of immersion heater switch-on's and errors propagating around the system as a result of assuming a constant input to each module. An alternative approach to reducing the time step from 1 to 41hr would have been to employ an iterative procedure at the end of each time step, as used in most simulation codes where the input to the module is assumed constant throughout the ensuing time step. However, in either case the resulting computer operating time would have been greater than that required by the improved simulation code.
6. CONCLUSION
An improved computer code which represents a significant improvement over existing fixed step iterative simulation programs has been presented. The key features of the program are the simultaneous solution of all differential equations appearing in closed loop form, and the technique for advancing the simulation form "control function" to "control function". Two simulation example problems demonstrating the application of the improved code to studies in control and long-term system performance were also presented.
(ii) Example 2 In the next problem the long-term performance of the combined solar hot water and space heating system illustrated in Fig. 1. is determined. The simulation data are given in Table A2.2, Appendix 2. The system is simulated for a single year of operation using hour by Key
60.0-
o-Improved Coml:xJter Code (h=lhr) =-Fixed Step Computer Code (h=l/4hr) (Imoginory Flow Conditon}
Weother Doto KEW (1963-1964) lZ.th Moy Pump
Off
50.0-
40.0-
Pump
On Z
~ 30-0~ I
,• /•"
"
Preheot Tonk Temp Trajectory
o - - -o
-- - o -
-- o - - - o - "
\ ,, \
• X \ \
,)
k
? ..,.."-..---Collector FluidOutlet
20.0-
"\
Temp.
o,
% ~o ~'o
Trajectory
o
o
10.0-
Time (hrs)
1i 1'~ 15 16 17 lb 19 20 2i
Fig. 5. Collector and preheat tank fluid temperature trajectories.
2'2 2'3 2t.
106
P.B. HOWELLSand R. H. MARSHALL 30.0-
6"L°
28.0' 26.0-
/
5.0
2A.O.
Key
,/
Computer "Dme {CDC 7600 Computer)
22.020.018.016.0-
,i,1'1!
,, <.o-!
% Imbolonce
- - - -
14.012.0IO.OE o
8.0.
~ 2.0-
6.0-
?
4.0-
°c w
/ 1.0
/
2.01.0
I
0,0001
0,001
001
0.1
Maximum Error Tolerance Per Time Step (£Mox)
Fig. 6. Accuracy and computer time as a function of the user-specified error tolerance (Emax).
Table 2. Combined solar hot water and space heating system performance data (Emax= 0.01°C)
Collector input
177.3
Collector output
GJ
32.87 GJ
Thermal store input
32.87 GJ
Thermal store output
28.6
GJ
Thermal store loss
4.3
GJ
Preheat tank input
8.4
GJ
Preheat tank output
8.0
GJ
Preheat tank loss
0.4
GJ
Auxiliary energy required by hot water system
5.4
GJ
Hot water load
13.7
GJ
Space heating heat exchanger input
20.2
GJ
Auxiliary energy required by space heating system.
53.6
GJ
Space heating load
73.9
GJ
% Hot water load supplied by solar energy
58%
% Space heating load supplied by solar energy
27%
In the first problem, a domestic solar water heating system was considered in which the temperature of the fluid at the hot water tank is controlled by an immersion heater and an on/off controller. The rapid switching
action of the immersion heater controller posed a very severe test. However, by advancing from each control action to the next (or to the end of the hour) the improved code resulted in an accurate simulation of the
An improved computer code for the simulation of solar heating systems controller actions without constraining the nominal time step to be less than the shortest time interval between 1 switches, calculated to be io hr. For the same problem, the system temperature trajectories were also computed and compared with those of a fixed step simulation program in which the input to each module was assumed to be constant during the ensuing time step. This comparison showed that the trajectories agree well even at a nominal 1 hr time step for the improved code, compared with the "hr time step used by the fixed step program. As a result of the larger nominal time step size a 43 percent saving in computer time is realised for the improved code. This saving is a direct result of the simultaneous solution approach which allows for time varying input to each module, thus avoiding any reduction in the simulation time step or the need for any form of iterative procedure to ensure the accuracy of the result. In the second problem a more complex combined solar hot water and space heating system was considered. The system was simulated for a single year's operation using h o u / b y hour meteorological data for different values of the user-specified error tolerance, E . . . . Based on a comparison between computer time and the largest energy transfer, which occurred at the thermal store, the optimum error tolerance was found to be 0.01°C. This value yielded a yearly energy balance correct to within 1 percent while requiring only 9 sec of computer operating time for the year's simulation. These example problems demonstrated that the improved code can be used for detailed system control as well as long-term system performance studies. Furthermore, because the system equations are solved simultaneously marching from "control function" to "control function" a high degree of accuracy for a reduced computer time is realised. NOMENCLATURE Ac collector area, m2
{A} matrix elements C~, specific heat capacity, kJ kg ~K {D} matrix elements Oi.i E error bound estimate Em~x user-specified maximum error per time step f~ {F} matrix elements F' plate efficiency factor h time step size, sec HO) irradiance on the tilted surface, Wm 2 MERS integration method selection parameter mass flowrate, kg sec mCp minimum capacity fluid flowrate, WK N matrix dimension Cc collector heat capacity, JK (T) fluid temperatures external to the component modules °C (T)r fluid temperatures internal to the component modules °C T,(t) ambient air temperature, °C Tm inlet fluid temperature, °C Tout outlet fluid temperature, °C rt fluid temperature at the preheat tank (domestic solar water heating system), °C T2 fluid temperatures at the output of the collector (domestic solar water heating system), °C r . fluid temperature at the hot water tank, °C rsH space heating demand temperature, °C To domestic hot water demand temperature, °C ai.i
107
TRM room temperature, °C
UL overall collector heat loss coefficient, WK ~ m 2 U heat loss coefficient UABLDC; building heat loss coefficient/area product, WK-t V(n volume, me Greek ol solar absorptance 8h time when a change in control function occurs measured back from the end of the time step "h", sec heat exchanger effectiveness o angle between the beam radiation and the normal to the collector t7 weighting function "7d diffuse radiation transmittance coefficient beam radiation transmittance coefficient "rb (Ol --->oj ~ RIKM integration coefficients
Subscripts I internal loop parameters E external loop parameters CL collector loop DL domestic hot water loop SL space heating loop PT preheat tank TS thermal store HT hot water tank Superscripts n , n + l time level Other { } square matrix ( ) column matrix rate of change with time
REFERENCES
1. TRNSYS, TRNSYS Simulation Manual. Solar Energy Laboratory, University Wisconsin-Madison (1979). 2. R. H. Marshall, A. E. Churches and J. A. Reizes, A hybrid integration scheme for the solution of viscous compressible flows. Appl. Math. Model. 3, 459--465. (1979). 3. P. B. Howells and R. H. Marshall, A fast and accurate integration scheme for coupled stiff differential equations, pp. 45-59. Proc. Conf. Numerical Methods for Coupled Problems. Pineridge Press, London (1981). 4. J. A. Duffle and W. A. Beckman, Solar Energy Thermal Processes. Wiley, New York (1974). 5. W. M. Kays and A. L. London, Compact Heat Exchangers. McGraw-Hill, New York (1958).
APPENDIX 1
Table AI.1. Component modules currently available (a) (b) (c) (d) (e) (f) (g) (h) (i) (j) (k) (I) (m) (n)
Solar radiation processing module Pumpmodule Cold water tank module Combinedflat plate collector/heat exchanger module Combinedpreheat tank/heat exchanger module Hot water tank module Fan coil space heating module Water radiator space heating system module Fully mixed thermal store module Stratifiedstore module Phasechange thermal store module Clock module Energy meter module On/Offcontrol module
P.B. HOWELLSand R. H. MARSHALL
108 APPENDIX 2
Table A2.2. Combined solar hot water and space heating system simulation data
Table A2.1. Domestic solar water heating system simulation data Parameter Ac ~ICL
F'
C 7d 7b Ot
G VOLer
Uer VOLHr
Um Integration scheme
Emax Auxiliary water heater rating Weather data
Parameter
Value
Ac
50 m 2 6.0WK ~m 2 0.017 kg sec i m-2 0.9 10.0 kJ K -t m -2 0.76 0.84 (-0.0119 + 3.42 cos 0 -3.97 cos: 0 + 1.56 cos 3 0) 0.9 4.2kJkg 1K i 1.0 2001 0.4Wm 2K h 0.068 kg sec 60°C 751m 2 0.4Wm 2K 20°C 20°C 400 WK t 0.8 0.017 kg sec -1 m -2 1.007 kJ kg K t 0.008 kg sec -t m -2 RIKM Kew (England) 1963-64
Value 4.0m 2 6.0 WK-t m 2 0.017 kg sec ~m 2 0.9 10kJ K-I m 2 0.76 0.84 (-0.0119+ 3.142 COS 0 -3.97 COS2 0 + 1.56 COS3 0) 0.9 4.2 kJ kg I K I 2001 0.4Wm 2K 1301 0.4 Wm -2 K-i RIKM 0.01 K 3 KW Kew (England) 1963-64.
Uc ~hcc F' C~ Td "rb Ol
G ~.cc
VOLpT UpT
]]//DL TD VOLts
Urs
TsH TRy UABLtm ~.SL
tilsg
Cpair r/lair Integration scheme Weather data