Transfer functions of solar heating systems for dynamic analysis and control design

Renewable Energy 77 (2015) 64e78

Contents lists available at ScienceDirect

Renewable Energy journal homepage: www.elsevier.com/locate/renene

Transfer functions of solar heating systems for dynamic analysis and control design
rd Kicsiny* Richa
n University, Pa
ter K. u. 1., 2100 Go € do €llo }, Hungary Department of Mathematics, Institute for Mathematics and Informatics, Szent Istva

a r t i c l e i n f o

a b s t r a c t

Article history: Received 6 August 2014 Accepted 1 December 2014 Available online

Mathematical modelling is the theoretically established tool for developing solar heating systems, e.g. with using transfer functions. If we know the transfer functions of the system, the outlet temperature can be predicted as a function of the input variables (solar irradiance, inlet temperature, environment temperatures), dynamic analysis can be carried out, and stable system control can be effectively designed based on the well-tried methods of control engineering. For these purposes, new, validated transfer functions for solar heating systems are worked out in this study based on a mathematical model, which can be found in the literature and has been applied successfully in the field. The transfer functions are used for dynamic analysis and control design of solar heating systems. The dynamic analysis is presented and the efficiency of the proposed stable control is demonstrated with respect to a real solar heating system. © 2014 Elsevier Ltd. All rights reserved.

Keywords: Solar heating systems Transfer functions Dynamic analysis Control design PI control

1. Introduction Mathematical modelling is the theoretically established tool for developing solar heating systems, e.g. with using transfer functions. Ordinary differential equation (ODE) models are widely used as they are relatively simple and easy to handle. Among the current collector models, the nonlinear ODE model proposed by Perers and Bales [1] may be the most widely used one, which is the improved version of the quasi-dynamic model from the standard EN 12975 [2]. If there is an external heat exchanger in the system, it can be modelled with the well-known effectiveness-NTU method [3], or the separate sides of the heat exchanger can be assumed to have homogeneous temperatures and can be modelled with ODEs [4]. The pipes of the system may be modelled with ODEs (assuming homogeneous pipe temperatures), or partial differential equations (PDEs) (with the one-dimensional linear heat transfer PDE) [5,6]. Solar storages can be also modelled with ODEs. See Ref. [7] for ODEs of mixed storages and stratified storages. Georgiev [8] connected a distributed (PDE) collector model and a mixed (ODE) storage model to describe a collector-storage

* Tel.: þ36 28522000/1414; fax: þ36 28410804. E-mail address: [email protected]. http://dx.doi.org/10.1016/j.renene.2014.12.001 0960-1481/© 2014 Elsevier Ltd. All rights reserved.

system. After connecting the models of the working components, solar heating systems are generally modelled with ODEs [9e11]. In Refs. [4,12], collector-heat exchanger-storage systems are modelled with a linear (multidimensional) ODE, which will be used in the present paper. This model is validated [4,12,13] or partly validated [14] and accurate enough for different successful applications [10,12,13,15] as well as its improved version in Ref. [13], where the pipes of the system are also modelled with ODEs. The advantage of the basic linear model of Ref. [4] is that it is simpler and easier to use than its extended linear version [13], its nonlinear version [16] or the delay differential equation model of Ref. [15]. The extended linear model of Ref. [13] and the model of Ref. [15] are roughly the same precise and they are more precise than the basic model of Ref. [4]. On the other hand, the extended linear model [13] is simpler and easier to use than its nonlinear version [13] or the model of [15]. Furthermore, the mentioned nonlinear models are approximately the same precise as their linear versions (see Ref. [13], cf [12,16] or see Ref. [17]). Thus the models of [13,15] are the most advantageous in view of accuracy while the linear model of [4] is the most advantageous in view of simplicity. In addition, the latter model is the basis for all other models in Refs. [13,15,16].
s et al. [9], collector transfer functions From the model of Buza have been determined and applied for the dynamic analysis of real collectors [18,19]. The present work extends these results by determining transfer functions for whole solar heating systems and applying them for the dynamic analysis of a real system.

R. Kicsiny / Renewable Energy 77 (2015) 64e78

Nomenclature Ac Ah Ahe cc ci ch Ic khe mh t

collector surface area (m2) surface area of the heat transfer inside the heat exchanger (m2) surface area of the heat exchanger to the environment (m2) specific heat capacity of the collector fluid (J/(kg K)) specific heat capacity of the fluid in the inlet loop (J/ (kg K)) specific heat capacity of the heat exchanger material (J/ (kg K)) (global) solar irradiance on the plane of the collector (W/m2) heat loss coefficient of the heat exchanger to the environment (W/(m2 K)) mass of the empty heat exchanger (kg) time (s)

Generally speaking, transfer function based modelling is relatively new and rare in the analysis of solar heating systems, especially, in domestic case. Besides the latter two references, some examples are the following: Amer et al. [20] solved a collector model with time and one space dimension for the fluid temperature using Laplace transformation. Huang and Wang [21] wrote a nonlinear two-node collector model into Laplace transformed form to gain transfer functions. Bettayeb et al. [22] used a two-node model to determine collector transfer functions for the fluid temperature and the absorbed solar energy. The most prevalent and simple control strategy is the on/off control for solar heating systems in domestic hot water (DHW) production working with constant flow rate, see e.g. Refs. [7,23]. Several controls using pump flow modulation are used in solar heating systems: Winn [24] compared on/off, I (integral) and PID (proportional integral differential) controls, Hirsch [25] compared €f [26] discussed on/ on/off, P (proportional) and hybrid controls. Lo off, P, I, PID, adaptive and certain optimal controls. Optimal controls often maximize the overall energy gain by flow rate modulation. See Refs. [27e29] for the case of no heat exchanger and [30] for the case of a counter flow heat exchanger. P and PI (proportional integral) controls for collectors are given and a PI control is worked out in details for a real collector field in Ref. [19] based on proposed collector transfer functions. Based on studies in the literature, not many improvements on control for solar heating systems used in domestic applications have been established in the recent few decades. In particular, transfer function based control is rather rare. Pasamontes et al. [31] serve with a further example on the control of the collector field of a solar cooling system based on the transfer function for a mathematical model with time delay. Transfer function based control is more prevalent for industrial processes e.g. for solar power or desalination plants [32e35]. In the present study, new, validated transfer functions for solar heating systems used primarily for domestic purposes are proposed and used for dynamic analysis and control design. According to a there appointed future research, the present study extends the results of [19], where transfer functions, dynamic analysis and controls have been proposed for solar collectors. The below worked out transfer functions are unique concerning the linear ODE model for solar heating systems in Ref. [4]. If the method for working out transfer functions for this basic model is presented, corresponding

Tc Tce Thc The Thh Ti UL vc vi Vc Vh εkh

h0 F rc ri

65

collector temperature ( C) collector environment temperature ( C) cold side temperature in the heat exchanger ( C) environment temperature of the heat exchanger ( C) hot side temperature in the heat exchanger ( C) inlet (fluid) temperature of the system ( C) overall heat loss coefficient of the collector (W/(m2 K)) (pump) flow rate in the collector loop (m3/s) (pump) flow rate in the inlet loop (m3/s) collector volume (m3) volume of the heat exchanger (m3) heat transfer coefficient inside the heat exchanger (W/ (m2 K)) collector optical efficiency (-) heat exchanger effectiveness (-) mass density of the collector fluid (kg/m3) mass density of the fluid in the inlet loop (kg/m3)

transfer functions can be derived rather straightforwardly for the more precise extended linear model of [13] based on the same method. That is why the linear model of [4] is essential and applied in this study to work out transfer functions. (Because of limits in volume, the transfer functions for the extended model of [13] cannot be derived and detailed here.) More precisely, the truncated version of the model of [4] (without modelling the solar storage) is used in the present study, the validation of which is also presented below based on measured data. This means that the here proposed transfer functions are also validated, since they form an alternative representation of the same mathematical model. The advantages are considerable: Knowing the transfer functions, dynamic analysis of solar heating systems can be carried out and feedback control can be effectively designed based on the welltried methods of control engineering. A control determined in such a way is generally much simpler than nonlinear and optimal controls and can follow the reference signal much more precisely than the most frequently used on/off controls. Perhaps, the simple applicability is the main advantage of the linear approach concerning transfer functions. The contributions of the present work in details are the following: 1. Based on a validated and successfully applied ODE model for solar heating systems [4], new transfer functions are mathematically derived (at the end of Section 3.1) and validated. The applicability of the transfer functions are interpreted with the dynamic analysis of a real system. 2. As a further important application of the transfer functions, closed-loop (PI) control design is given with stability criteria for solar heating systems using the methods of control engineering. The efficiency of the proposed control design is demonstrated based on simulation results. The present paper proposes all the concepts of work [19] (transfer functions, dynamic analysis, control design) for whole solar heating systems and not for collectors alone, so the present contributions are even more important and more general than the results of [19]. The paper is organized as follows: Section 2 describes the model for solar heating systems, for which the transfer functions are

66

R. Kicsiny / Renewable Energy 77 (2015) 64e78

determined in Section 3 and applied for the dynamic analysis of a real solar heating system. The transfer functions are used for designing a stable feedback control in Section 4, where the proposed control is applied and evaluated for the same real system. Section 5 contains final conclusions and future research proposals. See Ref. [36] for the used concepts of control engineering (Laplace transformation, transfer function, step response, P, PI controls, static error, stability, etc.). Maple [37], with its symbolic solver, is used below for mathematical derivations (e.g. by means of Laplace and inverse Laplace transformations) mainly to gain the transfer functions. Matlab [38] is used for validation and to implement and simulate the PI control for a real solar heating system.

This section presents the mathematical model for solar heating systems from Ref. [4] without the equation corresponding to the solar storage. The proposed transfer functions will be determined based on this truncated model, the validation of which is also carried out. 2.1. Mathematical model Consider the solar heating system in Fig. 1. The mathematical model of the system in Fig. 1 is the following (see also [4]):

dTc ðtÞ Ac h0 U Ac vc ¼ Ic ðtÞ þ L ðTce ðtÞ  Tc ðtÞÞ þ ðThh ðtÞ  Tc ðtÞÞ; dt rc cc Vc rc cc Vc Vc (1a) dThh ðtÞ r cc vc εk A ¼cm c ðTc ðtÞ  Thh ðtÞÞ þ c m h h V ðThc ðtÞ Vh h h h h dt þ rc cc þ rc cc h 2

2

 Thh ðtÞÞ þ c

Ahe khe =2 h mh

2

þ rc cc V2h

2

ðThe ðtÞ  Thh ðtÞÞ; (1b)

dThc ðtÞ rcv εk A ¼ c m i i i V ðTi ðtÞ  Thc ðtÞÞ þ c m h h V ðThh ðtÞ h h h h dt þ ri ci h þ ri ci h 2

2

 Thc ðtÞÞ þ c

2

Ahe khe =2 h mh

2

þ ri ci V2h

1. Here, vc is always maximal (constant) to keep the collector temperature at an all-time minimal level. This maximizes the efficiency of the collector and the solar heating system in case of any fixed value of vi. (Otherwise, vc could be variable in the models.) vi will be variable in Section 4, since it will be used as manipulated variable in the proposed control. 2. Basically, the effect of wind speed is not included in the model (1aec), nevertheless, one may consider it within the coefficient UL (see e.g. Ref. [39]). 2.2. Validation

2. Mathematical model and validation

2

Remark 2.1

2

ðThe ðtÞ  Thc ðtÞÞ: (1c)

Model (1aec) is applied for a real solar heating system used for
n UniDHW production. The system is installed at the Szent Istva € do €llo }, Hungary [40]. The system (SZIU system) has versity (SZIU) Go the following parameter values [12]: Ac ¼ 33.3 m2, h0 ¼ 0.74, rc ¼ 1034 kg/m3, cc ¼ 3623 J/(kg K), Vc ¼ 0.027 m3, UL ¼ 7 W/(m2 K), ch ¼ 464.76 J/(kg K), mh ¼ 37 kg, Vh ¼ 0.005 m3, εkh ¼ 2461.5 W/ (m2 K) Ah ¼ 2 m2, khe ¼ 5 W/(m2 K), Ahe ¼ 0.24 m2, ri ¼ 1000 kg/m3, ci ¼ 4200 J/(kg K), vc ¼ 0 or 16.3 l/min (¼0.000272 m3/s), vi ¼ 0 or 10.5 l/min (¼0.000175 m3/s). The calculations have been done in Matlab (and Matlab Simulink). Fig. 2 shows the Simulink diagram of model (1aec). For validation, the measured values of Ti, Ic, Tce, The, vc and vi are fed into the computer model of (1aec) as inputs as well as the initial values of Tc, Thh and Thc then the modelled and measured values of the outlet temperature Thc are compared. (Ti and Thc are measured on the inlet and outlet pipe just before and after the storage side of the heat exchanger, respectively.) Fig. 3 shows the comparison of the modelled and measured temperatures for two days: 22nd September 2012 and 2nd November 2012. The simultaneous operating states of the pumps (on/off) are also shown in the figure. The average of the absolute difference between the modelled and measured outlet temperatures is 2.6  C on 22nd September 2012 and 1.8  C on 2nd November 2012. In proportion to the difference between the maximal and minimal measured temperature values, the time average of the absolute difference (or absolute error) is 7.2% on 22nd September 2012 and 6.5% on 2nd November 2012 (c.f. [12]). It can be concluded that the model tracks the physical processes characteristically right with acceptable accuracy in view of several engineering aims (studying and developing solar heating systems). See e.g. Refs. [41], where the reasonability of such precision is reinforced for similar systems. Thus we can accept and apply the mathematical model (1aec). Remark 2.2 In fact, the outlet temperature Thc is difficult to model because of the small volume of the heat exchanger, which causes rather high and fast changes in the modelled temperature mainly when the pump flow rates are changing frequently (see Fig. 3). 3. Transfer functions In this section, the transfer functions are derived and applied for dynamic analysis. Here, vi is assumed to be constant. 3.1. Transfer function derivation

Fig. 1. Scheme of the solar heating system.

First, Eqs. (1aec) should be rewritten from time domain to Laplace domain with Laplace transformation:

R. Kicsiny / Renewable Energy 77 (2015) 64e78

67

Fig. 2. Computer model of the solar heating system.

 Ac h0 U Ac
sT c ðsÞ  Tc ð0Þ ¼ I c ðsÞ þ L T ce ðsÞ  T c ðsÞ rc cc Vc rc cc Vc  vc
þ T ðsÞ  T c ðsÞ ; Vc hh

sT hh ðsÞ  Thh ð0Þ ¼ c

rc cc vc h mh

2

þc

þ rc cc V2h


T c ðsÞ  T hh ðsÞ

εkh Ah h mh

2

þ

rc cc V2h

sT hc ðsÞ  Thc ð0Þ ¼ c

h mh

2

þ

þc

h mh

þ Wc3 ðsÞT ce ðsÞ;

þ ri ci V2h

T hh ðsÞ ¼ Whh0 ðsÞThh ð0Þ þ Whh1 ðsÞT c ðsÞ þ Whh2 ðsÞT hc ðsÞ þ Whh3 ðsÞT he ðsÞ;


T hh ðsÞ  T hc ðsÞ

 Ahe khe =2
T he ðsÞ  T hc ðsÞ ; Vh h mh þ r c i i2 2

(3a)

(2b)


T i ðsÞ  T hc ðsÞ

εkh Ah 2

þc

ri ci V2h

T c ðsÞ ¼ Wc0 ðsÞTc ð0Þ þ Wc1 ðsÞI c ðsÞ þ Wc2 ðsÞT hh ðsÞ


T hc ðsÞ  T hh ðsÞ

 A k =2
þ c m he he V T he ðsÞ  T hh ðsÞ ; h h h 2 þ rc cc 2 ri ci vi

(2a)

where the variables in Laplace domain are marked with overbars, s is the independent complex variable in Laplace domain, Tc(0), Thh(0), Thc(0) are the initial values of the state variables Tc, Thh, Thc. It is a great advantage that the system of linear ODEs (1aec) has been simplified to the system of linear algebraic equations (2aec). Rearranging equations (2aec), we get:

(3b)

T hc ðsÞ ¼ Whc0 ðsÞThc ð0Þ þ Whc1 ðsÞT i ðsÞ þ Whc2 ðsÞT hh ðsÞ þ Whc3 ðsÞT he ðsÞ; (2c)

Fig. 3. Validation of the used mathematical model.

(3c)

68

R. Kicsiny / Renewable Energy 77 (2015) 64e78

where

Wc0 ðsÞ ¼

tc ; tc s þ 1

Wc1 ðsÞ ¼

Wc2 ðsÞ ¼

tc vc $ ; t c s þ 1 Vc

Wc3 ðsÞ ¼

tc U Ac $ L ; tc s þ 1 rc cc Vc

thh rc c c v c thh εkh Ah $ $ ; Whh2 ðsÞ ¼ ; Vh V thh s þ 1 ch mh thh s þ 1 ch mh þ rc cc þ rc cc h 2 2 2 2 thh Ahe khe =2 thc thc ri ci vi $ ; Whc1 ðsÞ ¼ $ ; Whc0 ðsÞ ¼ ; Whh3 ðsÞ ¼ V V thh s þ 1 ch mh thc s þ 1 thc s þ 1 ch mh þ rc cc h þ ri ci h 2 2 2 2 thc εkh Ah thc Ahe khe =2 $ $ ; Whc3 ðsÞ ¼ ; Whc2 ðsÞ ¼ V V thc s þ 1 ch mh thc s þ 1 ch mh þ ri ci h þ ri ci h 2 2 2 2 Whh0 ðsÞ ¼

thh ; thh s þ 1

tc Ac h0 $ ; tc s þ 1 rc cc Vc

Whh1 ðsÞ ¼

where tc, thh, thc are the so-called time constants of the collector and the hot and cold sides of the heat exchanger, respectively:

tc ¼

1 UL Ac rc cc Vc

thc ¼

þ Vvcc

; thh ¼

ch mh þ rc cc Vh ; 2ðrc cc vc þ εkh Ah Þ þ Ahe khe

Thh(0), Thc(0) are assumed to be zero. Based on Eq. (4), the transfer functions for each input are the following:

T hc ðsÞ T ðsÞ T ðsÞ T ðsÞ ¼ W1 ðsÞ; hc ¼ W2 ðsÞ; hc ¼ W3 ðsÞ; hc ¼ W4 ðsÞ: T i ðsÞ I c ðsÞ T ce ðsÞ T he ðsÞ

ch mh þ ri ci Vh : 2ðri ci vi þ εkh Ah Þ þ Ahe khe

The response of the output temperature to the initial conditions can be similarly characterized with the following functions:

Solving the algebraic equations (3aec) for T hc ðsÞ, we get:

T hc ðsÞ ¼ Wi1 ðsÞTc ð0Þ þ Wi2 ðsÞThh ð0Þ þ Wi3 ðsÞThc ð0Þ þ W1 ðsÞT i ðsÞ þ W2 ðsÞI c ðsÞ þ W3 ðsÞT ce ðsÞ þ W4 ðsÞT he ðsÞ; (4) where

T hc ðsÞ T hc ðsÞ T hc ðsÞ ¼ Wi1 ðsÞ; ¼ Wi2 ðsÞ; ¼ Wi3 ðsÞ: Tc ð0Þ Thh ð0Þ Thc ð0Þ According to the superposition principle of linear systems, the resultant effect of the inputs and the initial temperatures on the output is the simple sum of the single effects of each input and initial temperatures. Eq. (4) expresses this principle in Laplace domain.

Wi1 ðsÞ ¼

Whc2 Whh1 Wc0 Whc2 Whh0 Whc0 ð  1 þ Whh1 Wc2 Þ ; Wi2 ðsÞ ¼ ; Wi3 ðsÞ ¼ ; 1 þ Whh1 Wc2 þ Whc2 Whh2 1 þ Whh1 Wc2 þ Whc2 Whh2 1 þ Whh1 Wc2 þ Whc2 Whh2

W1 ðsÞ ¼

Whc1 ð  1 þ Whh1 Wc2 Þ Whc2 Whh1 Wc1 Whc2 Whh1 Wc3 ; W2 ðsÞ ¼ ; W3 ðsÞ ¼ ; 1 þ Whh1 Wc2 þ Whc2 Whh2 1 þ Whh1 Wc2 þ Whc2 Whh2 1 þ Whh1 Wc2 þ Whc2 Whh2

W4 ðsÞ ¼

Whh3 Whc2  Whc3 þ Whc3 Whh1 Wc2 ; 1 þ Whh1 Wc2 þ Whc2 Whh2

where, for simplicity, the independent complex variable s was not indicated everywhere. In systems engineering approach, the solar heating system has an output variable Thc and input variables according to Fig. 4. The transfer functions are the Laplace transformed form of the output T hc ðsÞ divided with the Laplace transformed form of the inputs T i ðsÞ, I c ðsÞ, T ce ðsÞ, T he ðsÞ. When the transfer function corresponding to a given input is determined, the other inputs and Tc(0),

Fig. 4. Block diagram of the solar heating system.

Remark 3.1 The above proposed transfer functions are validated, since they form an alternative representation of the mathematical model (1aec) validated in Section 2.2.

3.2. Dynamic analysis The transfer functions can be used for the dynamic analysis of solar heating systems. The dynamic features of a system can be well characterized with its unit step responses. The unit step response corresponding to a given input is the response (output) of the system to the input in time domain, supposing that the input is the unit step input, while the other inputs and the initial values of the state variables are zero. The unit step input is the following:

R. Kicsiny / Renewable Energy 77 (2015) 64e78

 InputðtÞ ¼

0; 1;

t < 0; ; t  0;

69

(5)

the Laplace transformed form of which is the following:

1 InputðsÞ ¼ : s

(6)

The unit step response as output can be given with the following formula in Laplace domain using the transfer function for the corresponding input (W(s)): Fig. 5. Output response of the SZIU system to the unit step of Ti.

1 OutputðsÞ ¼ WðsÞ : s

(7)

The unit step response in time domain can be determined from OutputðsÞ:

  1 OutputðtÞ ¼ L1 WðsÞ ; s

(8)

where L1 denotes the inverse Laplace transformation. As a part of the dynamic analysis, the effect of the initial values of the state variables on the system can be also investigated in time domain based on the inverse Laplace transformed form of the product of the given initial condition and the corresponding transfer function (Wi(s)) (here, the inputs and the other initial values are zero again):

¼ Initial condition$L

½Wi ðsÞ:

(11) The unit step response corresponding to Tce is the following:

Thc ðtÞ ¼ 0:2  0:00019e0:58t þ 0:025e0:054t  0:223e0:0055t : (12) The unit step response corresponding to The is the following:

Thc ðtÞ ¼ 0:0013 þ 0:0000026e0:58t  0:00054e0:054t  0:00075e0:0055t :

OutputðtÞ ¼ L1 ½Wi ðsÞ$Initial condition 1

Thc ðtÞ ¼ 0:0210:00002e0:58t þ0:0026e0:054t 0:024e0:0055t :

(9)

Apply the described dynamic analysis for the solar heating system. The unit step responses of the system corresponding to the inputs Ti, Ic, Tce, The are the following:

The response corresponding (Tc(0) ¼ 1  C) (see also Fig. 7):

(13) to Tc(0)

is

the

following

Thc ðtÞ ¼ 0:05e0:58t  0:577e0:054t þ 0:528e0:0055t :

(14)

    1 1 Thc ðtÞ ¼ L1 W1 ðsÞ ; Thc ðtÞ ¼ L1 W2 ðsÞ ; s s     1 1 and Thc ðtÞ ¼ L1 W4 ðsÞ : Thc ðtÞ ¼ L1 W3 ðsÞ s s The responses corresponding to the initial conditions Tc(0), Thh(0) and Thc(0) are the following:

Thc ðtÞ ¼ Tc ð0ÞL1 ½Wi1 ðsÞ; Thc ðtÞ ¼ Thh ð0ÞL1 ½Wi2 ðsÞ and Thc ðtÞ ¼ Thc ð0ÞL1 ½Wi3 ðsÞ:

3.2.1. Application for the SZIU system The above dynamic analysis is applied for the SZIU system. The tap water is led into the storage side of the heat exchanger with temperature Ti. The output fluid is the required DHW with temperature Thc (see Fig. 1). Here, the pumps are considered switched on. The response functions to each input are provided below. Some of them are also depicted (but not all because of limits in volume). The unit step response corresponding to Ti is the following (see also Fig. 5):

Fig. 6. Output response of the SZIU system to the unit step of Ic.

Thc ðtÞ ¼ 0:8  0:031e0:58t  0:34e0:054t  0:43e0:0055t : (10) The unit step response corresponding to Ic is the following (see also Fig. 6):

Fig. 7. Output response of the SZIU system to Tc(0) ¼ 1  C.

70

R. Kicsiny / Renewable Energy 77 (2015) 64e78

2. The above dynamic analysis is essentially theoretical, since it requires a perfect unit step as one of the inputs while the others are constant, furthermore, the system must be in a perfect equilibrium at the beginning. Of course, such conditions do not occur in practice. Nevertheless, this dynamic analysis is a standard and widely used way for characterizing dynamic systems.

4. System control

Fig. 8. Output response of the SZIU system to Thc(0) ¼ 1  C.

The response corresponding to Thh(0) is the following (Thh(0) ¼ 1  C):

Thc ðtÞ ¼ 0:48e0:58t þ 0:419e0:054t þ 0:065e0:0055t :

(15)

The response corresponding to Thc(0) is the following (Thc(0) ¼ 1  C) (see also Fig. 8):

Thc ðtÞ ¼ 0:466e0:58t þ 0:473e0:054t þ 0:061e0:0055t :

Stable control can be designed for solar heating systems using the transfer functions and the methods of control engineering. The aim is to change the outlet temperature as controlled variable according to a given reference function in time by proper modulation of the inlet pump flow rate vi as manipulated variable. Here, vc is constant, vi is variable (see Note 1 in Remark 2.1). In this approach, functions Ti, Ic, Tce, The are disturbances. This control is summarized in Fig. 9. Now, vi(t) is not constant, so not all coefficients are constant in system (1aec) even (1aec) is not linear in the variables Thc(t), Ti(t) and vi(t), because of the products vi(t)Thc(t) and vi(t)Ti(t) in (1c). Therefore, the classical linear methods of control engineering cannot be directly applied. First of all, the model (1aec) should be linearized in a convenient operating point.

(16)

If the different inputs and initial conditions effect at the same time, the resultant output can be determined with a simple sum of functions (10)e(16) according to the superposition principle:

Thc ðtÞ¼1:0210:00002e0:58t þ0:0026e0:054t 0:024e0:0055t : (17) For visibility, the responses in Figs. 5e8 are presented for different time periods. Remark 3.2 1. According to functions (10)e(16), the largest effect of the inputs on the outlet temperature Thc is evoked by the unit change of Ti, and the lowest effect is evoked by the unit change of The (since khe ¼ 5 W/(m2 K) is very low because of the insulation of the heat exchanger). Among the initial conditions, Thc(0) has the largest and Thh(0) has the lowest effect. It seems that each input has rather low effect. For example, if Ic increases by 1 W/m2, Thc increases only by approximately 0.02  C (see Fig. 6). It is reasonable considering the rather high flow rates vc ¼ 16.3 l/min, vi ¼ 10.5 l/min. Nevertheless, the output is proportional to the inputs according to the linearity principle, so higher inputs involve (proportionally) higher output: e.g. if Ic increases by 100 W/m2, Thc increases by approximately 2  C.

4.1. Model linearization Practically, such an equilibrium of model (1aec) is chosen as operating point, which represents a sort of “average” operating condition, that is to say, in case of which functions Tc(t), Thh(t), Thc(t), Ti(t), Ic(t), Tce(t), The(t) are constants, where each constant is the approximate average between the upper and lower limits of its real 0 , T 0 , T 0 , I 0 , T 0 and T 0 denote the corpossible values. Let Tc0 , Thh c ce i hc he responding constants at such an operating point. The right hand side of (1aec) are zero at the operating point, since it is an equilibrium:

0¼

 v   Ac h0 0 U Ac  0 c 0 Tce  Tc0 þ Thh Ic þ L  Tc0 ; rc cc Vc rc cc Vc Vc

0 ¼c

rc cc V_ c h mh 2

þ

h mh

2

þ rc cc V2h

ri ci v0i h mh

2

þc

 0 Tc0  Thh þc

Ahe khe =2 

þc

0¼c



rc cc V2h

þ ri ci V2h



2

þ ri ci V2h

þ

rc cc V2h

0 0 Thc  Thh

 (18b)

0 0 The ;  Thh

 0 Ti0  Thc þc

Ahe khe =2 

h mh





εkh Ah h mh 2

(18a)



εkh Ah h mh

2



þ ri ci V2h

0 0 Thh  Thc

0 0 The ;  Thc



(18c)

The inlet flow rate at the operating point from Eq. (18c) is

v0i

   . 0  T0 0  T0 þ Ahe khe Thc 2 εkh Ah Thc hh he   ¼ : 0 rc Ti0  Thc i i

(19)

Eq. (1c) has the following form:

dThc ðtÞ ¼ f ðThh ðtÞ; Thc ðtÞ; Ti ðtÞ; The ðtÞ; vi ðtÞÞ; dt

Fig. 9. Block diagram of the solar heating system in view of control.

(20)

based on which the linearized version of (1c) at the operating point is the following:

R. Kicsiny / Renewable Energy 77 (2015) 64e78

71

   dThc ðtÞ vf  0 0 0 0 0   vf  0 0 0 0 0  0 0 0 0 ¼ f Thh Thh ; Thc ; Ti ; The ; vi $ Thh ðtÞ  Thh þ T ;T ;T ;T ;v $ ; Thc ; Ti0 ; The ; v0i þ dt vThh vThc hh hc i he i  vf      vf vf  0 0 0 0 0   0 0 0 0 0 ðThc ðtÞ  Thc ; Thc ; Ti0 ; The ; v0i $ Ti ðtÞ  Ti0 þ þ Thh Thh ; Thc ; Ti ; The ; vi $ The ðtÞ  The þ $ vTi vThe vvi . 
0      ri ci vi þ εkh Ah þ khe Ahe 2   εkh Ah 0 0 0 0 þ Thh ; Thc ; Ti0 ; The ; v0i $ vi ðtÞ  v0i ¼ 0  Thh ðtÞ  Thh $ ch mh V ch mh V þ ri ci h þ ri ci h 2 2 2 2   0 0    ri ci Ti  Thc     ri ci v0i khe Ahe =2  0 0 þ þ Ti ðtÞ  Ti0 þ The ðtÞ  The vi ðtÞ  v0i Thc ðtÞ  Thc ch mh V ch mh V ch mh V þ ri ci h þ ri ci h þ ri c i h 2 2 2 2 2 2

sT~ c ðsÞ ¼ Eq. (1a,b) are linear with respect to all time-dependent functions, so the coefficients in these equations remain the same in the linearized model below (22aec). 0 , T ~ ðtÞ ¼ T ðtÞ  T 0 , Let T~ c ðtÞ ¼ Tc ðtÞ  Tc0 , T~ hh ðtÞ ¼ Thh ðtÞ  Thh hc hc hc 0 0 0, ~ ~ I c ðtÞ ¼ Ic ðtÞ  Ic , T~ ce ðtÞ ¼ Tce ðtÞ  Tce T i ðtÞ ¼ Ti ðtÞ  Ti , 0 , ~ vi ðtÞ ¼ vi ðtÞ  v0i , with which the linearized T~ he ðtÞ ¼ The ðtÞ  The model is

 v   Ac h0 ~ U Ac  ~ c ~ Ic ðsÞþ L T ce ðsÞ T~ c ðsÞ þ T hh ðsÞ T~ c ðsÞ ; rc cc Vc r c c c Vc Vc

sT~ hh ðsÞ ¼c

(22a)

h mh

2

þc þc

~ ðsÞ ¼ W 1

~ ~ ~ ~ 1 þ W hh1 ðsÞW c2 ðsÞ þ W hc2 ðsÞW hh2 ðsÞ

þc

(22b)

þ rc cc V2h þ

(23a)

 T~ c ðsÞ  T~ hh ðsÞ 

εkh Ah h mh 2

rc cc V2h

T~ hc ðsÞ  T~ hh ðsÞ



 Ahe khe =2  ~ T he ðsÞ  T~ hh ðsÞ ; Vh 2 þ rc cc 2



þ

ri ci V2h

T~ i ðsÞ  T~ hc ðsÞ 

εkh Ah h mh 2

þ

ri ci V2h



T~ hh ðsÞ  T~ hc ðsÞ



 Ahe khe =2  ~ T he ðsÞ  T~ hc ðsÞ Vh h mh 2 þ ri ci 2   rc 0 ~vðtÞ; þ c m i i V Ti0  Thc h h h 2 þ ri ci 2

(23c)

þc

(22c)

from which the resultant effect of the inputs is the sum of the single effects according to the linear superposition principle:

~ ðsÞT~ ðsÞ þ W ~ ðsÞ~I c ðsÞ þ W ~ ðsÞT~ ce ðsÞ þ W ~ ðsÞT~ ðsÞ T~ hc ðsÞ ¼ W 1 2 3 4 i he ~ 5 ðsÞ~v ðsÞ; þW i

~ ðsÞ ¼ ; W 2

(24)

where

~ ~ ~ W hc2 ðsÞW hh1 ðsÞW c1 ðsÞ ~ ~ ~ ~ 1 þ W ðsÞW ðsÞ þ W W hh1

c2

hc2

hh2 ðsÞ

;

~ ~ ~ ~ ~ ~ ~ ~ ~ W hc2 ðsÞW hh1 ðsÞW c3 ðsÞ ~ ðsÞ ¼ W hh3 ðsÞW hc2 ðsÞ  W hc3 ðsÞ þ W hc3 ðsÞW hh1 ðsÞW c2 ðsÞ; ; W 4 ~ ~ ~ ~ ~ ~ ~ ~ 1 þ W hh1 ðsÞW c2 ðsÞ þ W hc2 ðsÞW hh2 ðsÞ 1 þ W hh1 ðsÞW c2 ðsÞ þ W hc2 ðsÞW hh2 ðsÞ   ~ ~ ~ W hc4 ðsÞ  1 þ W hh1 ðsÞW c2 ðsÞ ~ ðsÞ ¼ ; W 5 ~ ~ ~ ~ 1 þ W hh1 ðsÞW c2 ðsÞ þ W hc2 ðsÞW hh2 ðsÞ

~ ðsÞ ¼ W 3

(23b)

h mh

ri ci v0i h mh

2

Rewrite (22aec) (with T~ c ð0Þ ¼ T~ hh ð0Þ ¼ T~ hc ð0Þ ¼ 0  C) into Laplace domain:

  ~ ~ ~ W hc1 ðsÞ  1 þ W hh1 ðsÞW c2 ðsÞ

sT~ hc ðsÞ ¼c



rc cc vc

 v   dT~ c ðtÞ Ac h0 ~ U Ac  ~ c ~ ¼ I c ðtÞþ L T ce ðtÞ T~ c ðtÞ þ T ðtÞ T~ c ðtÞ ; dt rc cc Vc rc cc Vc Vc hh   dT~ hh ðtÞ r cc vc ¼c m c T~ c ðtÞ  T~ hh ðtÞ V h h h dt 2 þ rc cc 2   εk A þ c m h h V T~ hc ðtÞ  T~ hh ðtÞ h h h 2 þ rc cc 2  A k =2  þ c m he he V T~ he ðtÞ  T~ hh ðtÞ ; h h h 2 þ rc cc 2 0   ~ ri ci vi dT hc ðtÞ ¼c m T~ i ðtÞ  T~ hc ðtÞ V h h h dt 2 þ ri ci 2   εk A þ c m h h V T~ hh ðtÞ  T~ hc ðtÞ h h h 2 þ ri ci 2  A k =2  þ c mhe he V T~ he ðtÞ  T~ hc ðtÞ h h h 2 þ ri ci 2   rc 0 ~vi ðtÞ: þ c m i i V Ti0  Thc h h h 2 þ ri ci 2

(21)

72

R. Kicsiny / Renewable Energy 77 (2015) 64e78

where

~ ðsÞ ¼ Wc1 ðsÞ; W c1

~ ðsÞ ¼ W ðsÞ; W c2 c2

~ ðsÞ ¼ W ðsÞ; W c3 c3

~ W hh1 ðsÞ ¼ Whh1 ðsÞ;

according to the notation in Section 3.1, and

ri ci v0i ~ thc ~ $c m ; W hc1 ðsÞ ¼ ~ thc s þ 1 h h þ ri ci Vh 2

2

~ thc εkh Ah $ ; Whc2 ðsÞ ¼ ~ thc s þ 1 ch mh þ ri ci Vh 2

2

where

(25)

(26)

ðsÞ ¼

~ ðsÞ W 1 ; ~ ~ ðsÞ 1 þ W c ðsÞW 5

(27)

~

ðsÞ ¼

~ ðsÞ W 2 ; ~ 5 ðsÞ ~ 1 þ W c ðsÞW

(28)

~~ W T

hc ;I c

~

ðsÞ ¼

~ ðsÞ W 3

;

(29)

~ ðsÞ W 4 : ~ ðsÞ ~ 1 þ W c ðsÞW 5

(30)

~ 5 ðsÞ ~ c ðsÞW 1þW

~ c ðsÞ ¼ AP ; W

~ c ðsÞ ¼ AP 1 þ 1 ¼ AP ð1 þ sTI Þ; W sTI sTI

(31)

(32)

(33)

where AP and TI are constant. Based on Section 4.1, it can be derived ~ c ðsÞW ~ 5 ðsÞ is in accordance with the general form that the product W of Eq. (31) in case of each control type:

P:

PI :

~ c ðsÞW ~ 5 ðsÞ ¼ cc;P W ~ ðsÞ; W 0 s0 ~ 5 ðsÞ ¼ cc;PI W ~ c ðsÞW ~ ðsÞ: W 0 s1

(34)

(35)

Let us take reference inputs of the following sort (when the disturbances are zero):

T~ hc;r ðtÞ ¼ cr t j ;

~

hc ;T i

ðsÞ ¼

~ ð0Þ ¼ 1. where i and cc are constant and W 0 Consider the possibility of P and PI controls:

PI :

The aim is to realize a stable closed-loop control of the solar heating system (1aec) according to Figs. 9 and 10 such that the outlet temperature Thc(t) follows a predetermined reference input Thc,r(t) in time precisely enough. This requirement corresponds to 0 . that T~ hc ðtÞ follows T~ hc;r ðtÞ, where T~ hc;r ðtÞ ¼ Thc;r ðtÞ  Thc ~ c such that the In effect, the problem is the determination of W control is stable with conveniently small static errors corresponding to the inputs T~ hc;r , T~ i , ~I c , T~ ce , T~ he . The transfer functions of the control in Fig. 10 with respect to the reference input and the disturbances T~ i , ~I c , T~ ce , T~ he are the following:

~~ W T

hc ;T he

P:

4.2. Control design

~ ~ ~ ~ ~ ðsÞ ¼ W c ðsÞW 5 ðsÞ ; W T hc ;T hc;r ~ c ðsÞW ~ 5 ðsÞ 1þW

~~ W T

~

~ 5 ðsÞ ¼ cc W ~ c ðsÞW ~ ðsÞ; W 0 si

2

ch mh þ ri ci Vh  ~ thc ¼  : 2 ri ci v0i þ εkh Ah þ Ahe khe

hc ;T ce

The so-called loop gain of the controlled system (the gain in ~ c ðsÞW ~ 5 ðsÞ. Let us write it going around the feedback loop) is W below (in (34) and (35)) in the general transfer function form cc ~ W 0 ðsÞ: si

~ thc A k =2 $ he he ; Whc3 ðsÞ ¼ ~ thc s þ 1 ch mh þ ri ci Vh 2 2   0  T0 T r c i i ~ thc i hc $ ; Whc4 ðsÞ ¼ ~ thc s þ 1 ch mh þ ri ci Vh 2

~~ W T

~ ~ W hh2 ðsÞ ¼ Whh2 ðsÞ; W hh3 ðsÞ ¼ Whh3 ðsÞ;

(36)

where j and cr are constant. If j ¼ 0, T~ hc;r ðtÞ is step function, if j ¼ 1, it is ramp function (only t  0 is of importance in our case). If i > j holds for i and j in (31) and (36), the static error of the control corresponding to T~ hc;r is zero. The control is required to follow the step inputs (j ¼ 0) precisely, that is with zero static error, so the PI control is chosen, since i ¼ 1 > j ¼ 0 (see (35)) is fulfilled in this case. Thus cc ¼ cc,PI, where cc,PI can be found in the Appendix as Eq. (A1) (which can be derived based on Section 4.1). The static error of the PI control in case of ramp function as reference input (j ¼ 1) is the following:

  c r er;s ¼ lim T~ hc;r ðtÞ  T~ hc ðtÞ ¼ : t/∞ cc

(37)

Consider the disturbance T~ i ðtÞ in the following form (when ~ T hc;r ðtÞ and the other disturbances are zero):

T~ i ðtÞ ¼ c1 t k ; Fig. 10. Closed-loop control of the solar heating system.

(38)

R. Kicsiny / Renewable Energy 77 (2015) 64e78

where k and c1 are constant. If k ¼ 0, T~ i ðtÞ is step function, if k ¼ 1 it is ramp function. ~ ðsÞ corresponding to T~ should be The transfer function W 1 i considered in the following form:

~ ðsÞ ¼ W 1

cT~ sl

i

~

~ T i ðsÞ; W 0

(39) ~

~ T i ðsÞ ¼ 1. where cT~ i is constant and W 0 If i > k þ l is fulfilled for i, k and l in (31), (38) and (39), the static error of the control corresponding to T~ i is zero. ~ ðsÞ is in accordance with (39), It can be derived similarly that W 1 where l ¼ 0 and cT~ i is according to (A2) and (A3) in the Appendix. Thus the static error of the control is zero if T~ i ðtÞ is step function (k ¼ 0). If T~ i ðtÞ is ramp function (k ¼ 1), the static error corresponding to T~ i ðtÞ is the following:

  c~ T e1;s ¼ lim T~ hc;r ðtÞ  T~ hc ðtÞ ¼ i c1 : t/∞ cc

(40)

Consider ~I c ðtÞ, T~ ce ðtÞ and T~ he ðtÞ similarly as in (38):

~I c ðtÞ ¼ c t m ; T~ ce ðtÞ ¼ c t n ; T~ ðtÞ ¼ c t q ; 2 3 4 he where m, n, q, c2, c3 and c4 are constant. It can be shown like above that the static error corresponding to ~I c (e2,s), T~ ce (e3,s) or T~ (e4,s) are zero if m ¼ 0, n ¼ 0, q ¼ 0, he respectively. If ~I c ðtÞ, T~ ce ðtÞ or T~ he ðtÞ is ramp function (m ¼ 1, n ¼ 1 or q ¼ 1, respectively), the corresponding static errors are the following, respectively:

e2;s ¼ e3;s ¼

e4;s ¼

c~Ic

c ; ccT~cce 2 c ; cc 3 cT~

he

cc

(41) (42)

c4 ;

(43)

where c~Ic , cT~ ce and cT~ are according to (A3)e(A6) in the Appendix. he The control parameters AP and TI should be such that the absolute values of the above static errors are not higher than a prefixed positive limit E, based on the constants cr, cT~ i , c~Ic , cT~ ce , cT~ , c1, he c2, c3, c4 and that cc ¼ cc,PI. Furthermore, AP and TI should be such that the control is stable. ~~ ~ Consider the transfer function with respect to T~ hc;r that is W T hc ;T hc;r (see (26)). The controlled system is stable corresponding to T~ hc;r if ~~ ~ the real parts of the zeros of the denominator of W are T hc ;T hc;r ~ c ðsÞW ~ 5 ðsÞ negative. It can be derived that the denominator 1 þ W has the following form:

dr4 s4 þ dr3 s3 þ dr2 s2 þ dr1 s1 þ dr0 :

(44)

The constants dr4, dr3, dr2, dr1, dr0 can be found in the Appendix ((A7)e(A11)). The following system of conditions is sufficient for the stability according to the Routh-Hurwitz criterion:

dr4 > 0; dr3 > 0; dr2 > 0;   dr3 dr1     dr3 dr1      dr4 dr2  > 0;  dr4 dr2  0 dr3

dr1 > 0;  0  dr0  > 0; dr1 

73

In sum, the mathematical task of designing a PI control is the following:   determine   the free control  parameters    AP and TI such that er;s   E, e1;s   E, e2;s   E, e3;s   E, e4;s   E hold, and conditions (SC) hold. Remark 4.1 The Routh-Hurwitz criterion assures the stability of not only the linearized but the corresponding nonlinear controlled system (where vi(t) is not constant), since the latter is also stable if the real ~ ~ ~ are negative acparts of the zeros of the denominator of W T hc ;T hc;r cording to Lyapunov. 4.2.1. Application for the SZIU system Let us give an appropriate PI control for the SZIU system. 0 ¼ 55  C, which is Consider an end-of-spring period in May. Let Thc 0 ¼ 56  C (a bit high enough for general domestic purposes. Let Thh 0 ), T 0 ¼ 15  C (approximate average tap water temhigher than Thc i 0 ¼ 20  C (approximate average daytime environment perature), Tce 0 ¼ 25  C (approximate average temperature in the temperature), The maintenance house). From these data, the remaining equilibrium values Tc0 , Ic0 and v0i can be determined according to (18aec): Tc0 ¼ 60.85  C, Ic0 ¼ 586.9 W/m2 (which is in accordance with the approximate average daytime solar irradiance of a clear day in May in Hungary [42]), v0i ¼ 0.000032 m3/s (¼1.92 l/min), the maximal value of vi is 10.5 l/min (see Section 2.2). As further limitation, it is assumed that vi(t) can vary between zero and the maximal value in 3 s that is the absolute value of the speed of flow rate changing (the absolute value of the acceleration of the flow) is 0.000058 m3/s2. Let the absolute values of the static errors (37), (40)e(43) be less or equal than 0.1  C, which is convenient for a DHW heating application. Let the controlled system be stable that is conditions (SC) are satisfied. Suppose that such fast changing disturbances (with high amplitudes) act on the controlled system simultaneously, which are rare even separately in real circumstances, thus they are even more improbable in the same time. Let us check in such a case whether the controlled system is still able to follow precisely enough a relatively fast changing reference input (with high amplitude). If the controlled system can follow well such an extreme reference input under such extreme disturbances (which have small probabilities), it will work even more satisfactorily under real circumstances. Thus the practical applicability of the control will be thoroughly supported. Sinusoidal inputs are given as the mentioned reference input and disturbances (see Figs. 11 and 12) with initial values corresponding to the above operating point of the SZIU system: 0 ¼ 55 Thc;r ð0Þ ¼ T~ hc;r ð0Þ þ Thc gradient ¼ 0.5  C/min;

 C,

amplitude ¼ 5

 C,

d  r0 > 0;   dr3 dr1 0 0     dr4 dr2 dr0 0     0 dr3 dr1 0  > 0:    0 d  r4 dr2 dr0 (SC)

~ c ðsÞW ~ 5 ðsÞ, so (SC) The denominators in (27)e(30) are also 1 þ W is sufficient for the stability of the controlled system corresponding to T~ i , ~I c , T~ ce and T~ he as well.

Fig. 11. Reference input to the controlled solar heating system.

initial

74

R. Kicsiny / Renewable Energy 77 (2015) 64e78

Fig. 12. Disturbances to the controlled solar heating system.

Ti ð0Þ ¼ T~ i ð0Þ þ Ti0 ¼ 15  C, amplitude ¼ 5  C, gradient ¼ 3  C/min; Ic ð0Þ ¼ ~I c ð0Þ þ Ic0 ¼ 586.9 W/m2, amplitude ¼ 150 W/m2, gradient ¼ 20 W/(m2s); 0 ¼ 20  C, amplitude ¼ 8  C, Tce ð0Þ ¼ T~ ce ð0Þ þ Tce gradient ¼ 1  C/min; 0 The ð0Þ ¼ T~ he ð0Þ þ The ¼ 25  C, amplitude ¼ 3  C, gradient ¼ 0.5  C/min.

initial initial initial initial

The control parameters are fixed AP ¼ 0.0000086, TI ¼ 0.4,  er;s   with which the stability is assured, and the requirements         0.1  C, e1;s   0.1  C, e2;s   0.1  C, e3;s   0.1  C, e4;s   0.1  C hold

even in case of maximal speed of changing of the inputs (when cr ¼ 0.5  C/min, c1 ¼ 3  C/min, c2 ¼ 20 W/(m2 s), c3 ¼ 1  C/min, c4 ¼ 0.5  C/min). Let us apply and test the above PI control for the original e not linearized e model (1aec) in case of the SZIU system. The model and the control have been fed into computer in Matlab (and Matlab Simulink) [38] environment according to Fig. 13. The initial state variables are set Tc(0) ¼ Thc(0) ¼ Thh(0) ¼ 20  C, hence the control has to reconcile relatively high initial error: Thc,r(0)  Thc(0) ¼ 35  C. Figs. 14 and 15 show the simulation results. Fig. 14 shows the reference input temperature Thc,r(t), the controlled outlet temperature Thc(t) and the manipulated pump flow rate vi(t).

Fig. 13. Computer model of the controlled solar heating system.

R. Kicsiny / Renewable Energy 77 (2015) 64e78

75

Fig. 14. Thc,r(t), Thc(t) (controlled variable) and vi(t) (manipulated variable).

Fig. 15 shows the control error Thc,r(t)  Thc(t). The results show that the absolute value of the control error decreases permanently below 5% of the initial error (below 1.75  C) within 25.0 min, the absolute value of the error decreases permanently below the prescribed limit 0.1  C within 26.4 min (settling time). Based on these values, the worked out PI control is considered satisfactorily fast and accurate in view of the control purpose (even in case of extreme reference input and disturbances). Remark 4.2 1. On an average clear day in May, at least 440 W/m2 solar irradiance is still expected for 5e5.5 h after the above settling time 26.4 min [41]. 436.9 W/m2 is the minimal irradiance in the simulation (see Fig. 12), which proved to be still enough to maintain the outlet temperature at the desired level (55  C) according to simulation experiences. Thus the controlled outlet

temperature Thc(t) should be able to follow the reference input (precisely enough) for at least 5e5.5 h. 2. Based on simulation data, 591.5 L DHW, at an average temperature of 55  C, has been produced in 5.5 h after the above settling time 26.4 min (on an average clear day in May), which covers 30% of the daily demand (1990 l [6]) of the consumer of the SZIU system (a kindergarten). The produced DHW can be stored in a water storage, which can be discharged during the day according to the all-time hot water demand. 3. The gained simulation results, where the control error converges to zero, underlie the statement of Remark 4.1 corresponding to the stability of the nonlinear controlled system (where vi(t) is not constant), since the nonlinear system model was controlled above (see Fig. 13). 4. The control error converges to zero with oscillation (see Fig. 15), which is normal in case of PI controllers, since the manipulated variable vi(t) is partially proportional to the time integral of the error (not only directly to the error). This causes recurring overruns/oscillation. In fact, vi(t) cannot be proportional to (the integral of) the error every time because of its lower (0 l/min) and upper (10.5 l/min) limits (see Fig. 14). That is why the graph of the control error is not smooth but has breaking points (see Fig. 15). 5. Conclusion

Fig. 15. Control error Thc,r(t)  Thc(t).

Generally speaking, transfer function based modelling is relatively new and rare in the analysis of solar heating systems, especially, in domestic case. In addition, not many improvements on control for solar heating systems used in domestic applications have been established in the recent few decades. In particular, transfer function based control design is rather rare despite of the simple applicability, which is the main advantage of the linear approach concerning transfer functions. Controls based on transfer functions are generally much simpler than nonlinear and optimal controls and can follow the reference signal much more precisely than the most frequently used on/off controls. This paper intended to contribute to fulfil these research gaps by proposing new,

76

R. Kicsiny / Renewable Energy 77 (2015) 64e78

validated transfer functions based on a validated system model and designing stable controls (in particular, a closed-loop PI control) by means of the transfer functions. The transfer functions have been also applied for the dynamic analysis of a real solar heating system, for which the proposed stable PI control has been applied to make the outlet temperature follow a reference input. According to a there appointed future research, this study has extended the results of [19], where transfer functions, dynamic analysis and controls have been proposed for solar collectors. Since the present paper proposes all of these concepts for whole solar heating systems, its contributions are even more important and more general than the results of [19]. As the main findings of the present paper, it can be stated that the proposed transfer functions are successful and relatively easy to apply (see the applications in the paper) for dynamic analysis and stable control design. Furthermore, the proposed PI control is appropriate for the corresponding control purpose because of its precision and rapidity even in case of fast changing reference input and disturbances with high amplitudes.

points, and can be also used for dynamic analysis and control design. Acknowledgement The author thanks the Editor for the encouraging help in the submission process, the anonymous Referees for their valuable
n Farkas and the Department of Physics comments and Prof. Istva and Process Control (SZIU) for the measured data. The author also thanks his colleagues in the Department of Mathematics in the Faculty of Mechanical Engineering (SZIU) for their contribution. Appendix The equations below have been derived based on the relations in the main text (with the symbolic solver of Maple) and are referred to in the text, where needed. For example, cc,PI has been ~ 5 ðsÞ and W ~ c ðsÞ in Sections 4.1 derived based on the expressions of W and 4.2.

  .  0 Ahe khe cc rc vc þ 2εkh Ah UL Ac þ 2εkh Ah rc cc vc þ 2cc rc vc UL Ac þ Ahe khe UL Ac cc;PI ¼  2AP ri ci  Ti0 þ Thc   TI 4rc cc vc UL Ac εkh Ah þ 2rc cc vc UL Ac khe Ahe þ 4rc cc vc UL Ac ci ri v0i þ 4εkh Ah rc cc vc khe Ahe þ 2ri ci $

(A1)

UL Ac khe Ahe v0i þ A2he k2he UL Ac þ 2khe Ahe rc cc vc ri ci v0i þ A2he k2he rc cc vc þ 4εkh Ah rc cc vc ri ci v0i þ 4εkh $  : UL Ac khe Ahe þ 4εkh Ah UL Ac ri ci v0i

In essence, the introduced dynamic analysis can be easily adapted for any solar heating system with an external heat exchanger. The proposed control design can be also comprehensively applied for solar heating applications, where the output temperature should follow some reference signal in time (e.g. solar

den ¼

cT~ i ¼ð2ri ci v0i ððAhe khe þ 2εkh Ah þ 2UL Ac Þcc rc vc þ ðAhe khe þ 2εkh Ah ÞUL Ac ÞÞ=den;

(A2)

where

     2Ahe khe v0i þ 4UL Ac v0i þ 4εkh Ah v0i cc rc vc þ 2Ahe khe v0i þ 4εkh Ah v0i UL Ac ci ri þ 4εkh Ah     þ2Ahe khe Ac UL þ A2he k2he þ 4εkh Ah Ahe khe cc rc vc þ A2he k2he þ 4εkh Ah Ahe khe Ac UL :



power and desalination plants). The manipulated variable was the flow rate in the inlet loop. The inlet temperature may be also used as manipulated variable, although, this solution is very rare. On the contrary, variable flow rate pumps are comprehensively applied in the practice of solar heating, and more generally, in building systems. Further research work may deal with the derivation of transfer functions corresponding to the more precise extended linear model of [13]. Such a work can follow in essence the steps of the derivation for the essential basic model of [4] presented above. Furthermore, so-called describing functions may be determined corresponding to nonlinear mathematical models of solar heating systems. Describing functions can be gained by means of harmonic linearization, which is a linearization method (differing from the one applied in this paper) corresponding to convenient operating

(A3)

c~Ic ¼ 4vc εkh Ah h0 Ac rc cc =den;

(A4)

cT~ ce ¼ 4vc εkh Ah UL Ac rc cc =den;

(A5)

cT~ ¼ðððAhe khe þ 4εkh Ah þ 2UL Ac Þcc rc vc he

þ ðAhe khe þ 4εkh Ah ÞUL Ac ÞAhe khe Þ=den: dr4 ¼ TI Vc cc ðch mh þ ri ci Vh Þrc ðch mh þ rc cc Vh Þ;

(A6) (A7)

R. Kicsiny / Renewable Energy 77 (2015) 64e78

dr3 ¼

dr2 ¼

dr1 ¼

    0 0 2 2 2ch mh rc cc Vc þ 2r2c c2c Vh Vc Ti0  2Thc ch mh rc cc Vc  2Thc rc cc Vh Vc ci ri AP þ 2Vc v0i þ vc Vh $        rc cc þ UL Ac Vh ch mh þ 2vc þ 2v0i Vh Vc þ Vh2 vc c2c r2c þ Ahe khe þ 2εkh Ah ÞVh Vc þ Vh2 UL Ac rc cc ci $       ri þ rc cc vc þ UL Ac c2h m2h þ vc Vh þ 2vc Vc Þc2c r2c þ 2Ahe khe þ 4εkh Ah ÞVc þ UL Ac Vh rc cc ch mh    þ Ahe khe þ 2εkh Ah Vh Vc c2c r2c TI ;

77



(A8)

       2UL Ac þ 2rc cc vc ch mh þ 2vc Vh þ 4vc Vc c2c r2c þ 2Ahe khe þ 4εkh Ah Vc þ 2UL Ac Vh rc cc $       0 0 0 0 0 mh ch þ  4Thc Ti0 þ  2rc cc vc Thc  2UL Ac Thc vc Vc  2Thc vc Vh c2c r2c þ  4εkh Ah  2Ahe khe Thc Vc       0 2Thc Vh UL Ac rc cc Þri ci AP þ 2rc cc vc v0i þ 2UL Ac v0i ch mh þ 2Vh vc v0i þ 4vc Vc v0i c2c r2c þ 4εkh Ah v0i         þ2Ahe khe v0i Vc þ 2vc þ 2v0i Ac UL þ 2εkh Ah vc þ Ahe khe vc Vh cc rc þ Ahe khe þ 2εkh Ah Ac UL Vh ci ri      þ4εkh Ah vc þ 2Ahe khe vc þ 2vc UL Ac cc rc þ 4εkh Ah þ 2Ahe khe Ac UL mh ch þ 4εkh Ah vc þ 2Ahe khe vc Vc          þ 2εkh Ah vc þ 2Ahe khe vc Vh c2c r2c þ A2he k2he þ 4εkh Ah Ahe khe Vc þ Ahe khe þ 2εkh Ah Ac UL Vh cc rc TI    0 0 2 2 þ 2ch mh rc cc Vc þ 2c2c r2c Vh Vc Ti0  2Thc ch mh rc cc Vc  2Thc cc rc Vh Vc ci ri AP ;

(A9)

     2Ahe khe þ 4εkh Ah þ 4UL Ac vc cc rc þ 2Ahe khe þ 4εkh Ah UL Ac Ti0 þ 4UL Ac  2Ahe khe       0 0 ci ri AP þ 4Ac UL v0i þ 4εkh Ah v0i þ 2Ahe khe v0i $ 4εkh Ah Thc vc cc rc þ 4εkh Ah  2Ahe khe Ac UL Thc       vc cc rc þ 2Ahe khe v0i þ 4εkh Ah v0i Ac UL ci ri þ 2Ahe khe þ 4εkh Ah Ac UL þ A2he k2he þ 4εkh Ah Ahe khe vc $         2rc cc vc þ 2Ac UL mh ch þ 4Vc þ 2Vh vc c2c r2c þ 2Ahe $ cc rc þ A2he k2he þ 4εkh Ah Ahe khe Ac UL TI þ        0 0 0 0 UL Ac  2Thc rc cc vc mh ch þ  2Thc Vh  4Thc Vc vc c2c r2c $ khe þ 4εkh Ah Vc þ 2UL Ac Vh rc cc Ti0 þ  2Thc     0 0 Vc  2Thc Vh UL Ac rc cc ci ri AP ; þ  4εkh Ah  2Ahe khe Thc

(A10)





  0 Ahe khe cc rc vc þ 2εkh Ah UL Ac dr0 ¼ 2AP ri ci  Ti0 þ Thc  þ 2εkh Ah rc cc vc þ 2cc rc vc UL Ac þ Ahe khe UL Ac : (A11) References [1] Perers B, Bales C. A solar collector model for TRNSYS simulation and system testing. A report of IEA SHC Task 26 Solar Combisystems. 2002. [2] CEN. European Standard Draft. European Committee for Standardisation; 1997. prEN 12975e2. [3] Kays WM, London AL. Compact heat exchangers. New York: McGraw-Hill; 1984.
s J, Farkas I. Solar domestic hot water system simulation using block[4] Buza oriented software. In: The 3rd ISES-Europe solar World Congress (Eurosun 2000), Copenhagen, Denmark, CD-ROM Proceedings; 2000. p. 9. [5] Orbach A, Rorres C, Fischl R. Optimal control of a solar collector loop using a distributed-lumped model. Automatica 1981;17(3):535e9. [6] Kicsiny R, Farkas I. Improved differential control for solar heating systems. Sol Energy 2012;86:3489e98. [7] Duffie JA, Beckman WA. Solar engineering of thermal processes. New York: John Wiley and Sons; 2006.

[8] Georgiev A. Simulation and experimental results of a vacuum solar collector system with storage. Energy Convers Manag 2005;46:1423e42.
A, Ne
meth R. Modelling and simulation of a solar [9] Buz
as J, Farkas I, Biro thermal system. Math Comput Simul 1998;48:33e46. [10] Kumar R, Rosen MA. Thermal performance of integrated collector storage solar water heater with corrugated absorber surface. Appl Therm Eng 2010;30:1764e8. [11] Zeghib I, Chaker A. Simulation of a solar domestic water heating system. Energy Procedia 2011;6:292e301. [12] Kicsiny R, Varga Z. Real-time state observer design for solar thermal heating systems. Appl Math Comput 2012;218:11558e69.
ki Cs. Extended ordinary differential equation models [13] Kicsiny R, Nagy J, Szalo for solar heating systems with pipes. Appl Energy 2014;129:166e76. [14] Kicsiny R. Multiple linear regression based model for solar collectors. Sol Energy 2014;110:496e506. [15] Kicsiny R. New delay differential equation models for heating systems with pipes. Int J Heat Mass Transf 2014;79:807e15. [16] Kicsiny R, Varga Z. Real-time nonlinear global state observer design for solar heating systems. Nonlinear Anal Real World Appl 2013;14:1247e64. [17] Kicsiny R. Real-time state observers for solar heating systems. Mech Eng Lett 2013;10:60e6. [18] Buz
as J. Block-oriented modeling of solar thermal systems. Thesis of the €do €llo }: Szent Istva
n University; 2009. p. 33 doctoral (Ph.D.) dissertation. Go [04.11.2014], http://szie.hu//file/tti/archivum/Buzas_thezis.pdf. [19] Buz
as J, Kicsiny R. Transfer functions of solar collectors for dynamical analysis and control design. Renew Energy 2014;68:146e55.

78

R. Kicsiny / Renewable Energy 77 (2015) 64e78

[20] Amer EH, Nayak JK, Sharma GK. A new dynamic method for testing solar flat-plate collectors under variable weather. Energy Convers Manag 1999;40:803e23. [21] Huang BJ, Wang SB. Identification of solar collector dynamics using physical model-based approach. J Dyn Syst Meas Control 1994;116:755e63. [22] Bettayeb M, Nabag M, Al-Radhawi MA. Reduced order models for flat-plate solar collectors. In: GCC Conference and Exhibition (GCC), 2011 IEEE; 2011. p. 369e72. [23] Kalogirou SA. Solar energy engineering. Academic Press; 2009. [24] Winn CB. Controls in solar energy systems. Advances in solar energy. American Solar Energy Society Inc.; 1983. [25] Hirsch UT. Control strategies for solar water heating systems. Master of Science Thesis. University of Wisconsin-Madison; 1985. €f G. Active solar systems. MIT Press; 1993. [26] Lo [27] Kovarik M, Lesse PF. Optimal control of flow in low temperature solar heat collectors. Sol Energy 1976;18:431e5. [28] Winn CB, Hull DE. Optimal controllers of the second kind. Sol Energy 1979;23: 529e34. [29] Badescu V. Optimal control of flow in solar collector systems with fully mixed water storage tanks. Energy Convers Manag 2008;49:169e84. [30] Hollands KGT, Brunger AP. Optimum flow rates in solar water heating systems with a counterflow exchanger. Sol Energy 1992;48:15e9.
[31] Pasamontes M, Alvarez JD, Guzm
an JL, Lemos JM, Berenguel M. A switching control strategy applied to a solar collector field. Control Eng Pract 2011;19:135e45. [32] Ayala CO, Roca L, Guzman JL, Normey-Rico JE, Berenguel M, Yebra L. Local model predictive controller in a solar desalination plant collector field. Renew Energy 2011;36:3001e12.

[33] Roca L, Guzman JL, Normey-Rico JE, Berenguel M, Yebra L. Robust constrained predictive feedback linearization controller in a solar desalination plant collector field. Control Eng Pract 2009;17:1076e88. [34] Roca L, Guzman JL, Normey-Rico JE, Berenguel M, Yebra L. Filtered Smith predictor with feedback linearization and constraints handling applied to a solar collector field. Sol Energy 2011;85:1056e67. [35] Fontalvo A, Garcia J, Sanjuan M, Padilla RV. Automatic control strategies for hybrid solar-fossil fuel plants. Renew Energy 2014;62:424e31. [36] Bakshi UA, Bakshi VU. Linear control systems. India: Technical Publications Pune; 2007. [37] Maplesoft. Maple 9 learning guide. Waterloo Maple Inc.; 2003. [38] Etter DM, Kuncicky D, Moore H. Introduction to MATLAB 7. Springer; 2004. [39] Amrizal N, Chemisana D, Rosell JI, Barrau J. A dynamic model based on the piston flow concept for the thermal characterization of solar collectors. Appl Energy 2012;94:244e50.
s J, La
gym

r I, Kaboldy E, Nagy L. In: Frankovic B, [40] Farkas I, Buza anyosi A, Kalma editor. A combined solar hot water system for the use of swimming pool and kindergarten operation, energy and the environment, vol. I. Opatija: Croatian Solar Energy Association; 2000. p. 81e8. [41] Kalogirou SA. Applications of artificial neural-networks for energy systems. Appl Energy 2000;67:17e35. [42] http://energetika.13s.hu/pub/_epuletenergetika_szakirany_/megujulo% 20energiaforrasok/Naplopo/EPGEP_Napkollektorok-1%5B1%5D.ppt [in Hungarian, 04.11.2014].