Simplified dynamical input–output modeling of plate heat exchangers – case study

Applied Thermal Engineering 98 (2016) 880–893

Contents lists available at ScienceDirect

Applied Thermal Engineering j o u r n a l h o m e p a g e : w w w. e l s e v i e r. c o m / l o c a t e / a p t h e r m e n g

Research Paper

Simplified dynamical input–output modeling of plate heat exchangers – case study Michal Fratczak, Pawel Nowak, Jacek Czeczot *, Mieczyslaw Metzger Faculty of Automatic Control, Electronics and Computer Science, Institute of Automatic Control, Silesian University of Technology, ul. Akademicka 16, 44100 Gliwice, Poland

H I G H L I G H T S

• • • •

Simplified dynamical modeling for plate heat exchangers was suggested. The optimization-based method for its tuning from the realistic measurement data was suggested. The suggested model was validated based on the measurement data. Stability and sensitivity analysis of the suggested model was provided.

A R T I C L E

I N F O

Article history: Received 3 September 2015 Accepted 2 January 2016 Available online 11 January 2016 Keywords: Plate heat exchangers Simplified dynamical modeling Orthogonal collocation method Practical validation

A B S T R A C T

In this paper, the simplified dynamical modeling of PHEs (plate heat exchangers) is suggested and validated based on realistic measurement data. The concept combines general distributed parameter model of double-pipe heat exchanger with its spatial discretization based on orthogonal collocation method. This approach preserves the distributed parameter nature and nonlinearities of the unit with relatively low numerical complexity and high modeling accuracy. For the suggested model, the tuning procedure is suggested and modeling accuracy is validated for realistic experimental measurement data. The stability of the model (both structural and numerical) is addressed and sensitivity analysis and comparison with finite difference approximation are presented. The results show that the suggested modeling and tuning method is useful for practical applications. © 2016 Elsevier Ltd. All rights reserved.

1. Introduction Nowadays, improving efficiency of industrial processes is considered as one of the most significant challenges for both control and process engineering. Many possibilities of reaching this goal were presented suggesting the concepts of such improvements. Some of them deal with economic impact of improvement in control. Bauer and Craig [1] presented the survey of methods used for assessing economic justification for applying new advanced process control systems to a process. Some of these methods are based on performance functions that incorporate economic aspects of variations of process variables around optimum operating point. Rhinehart et al. [2] presented the survey of advanced control methods that can be successfully applied in the industrial practice to provide improvement in control performance that can be linked to economic benefits. Yuan et al. [3] discussed the potential benefits resulting from simultaneous design of chemical processes and their control systems. This approach is based on optimizing performance function that

* Corresponding author. Tel.: (+48) 32 2371473; fax: (+48) 32 2372127. E-mail address: [email protected] (J. Czeczot). http://dx.doi.org/10.1016/j.applthermaleng.2016.01.004 1359-4311/© 2016 Elsevier Ltd. All rights reserved.

combines decision variables defined for both process design and its control. Another important possibility for process efficiency improvement results from reducing energy consumption combined with high safety and environmental goals. For this purpose, heat exchange and distribution systems that constitute an important part of many industrial process must be properly designed and controlled. In his textbook, Kemp [4] provides the framework for pinch analysis methods that allow for minimizing energy consumption for chemical processes. Bonhivers et al. [5] proposed the method for heat exchange network (HEN) retrofit that provides significant energy savings. Semkov et al. [6] proposed the method for assessing potential heat energy savings by reducing waste heat and proved its usefulness in food industry. Jie et al. [7] noticed that operating cost of district heating systems can be minimized by reducing pumping and heat loss cost. The improvement in heat exchanger design that results in lowering operating cost can be also provided by proper adjusting design margins compensating for process uncertainties [8]. Morrison et al. [9] suggested the method for minimizing total annualized cost of HEN for non-continuous processes by its proper design (reducing total exchange surface area). Bakošová and Oravec [10] proposed the predictive control strategy for HEN that is based on simplified model-based approach. Isafiade et al. [11]

M. Fratczak et al./Applied Thermal Engineering 98 (2016) 880–893

also found the cost-effective synthesis of HENs as an important issue and they accounted for optimal design in the presence of multiperiod operating conditions including variations in process requirements and of stages of operation. In industrial heat exchange systems, heat exchangers play an important role and their modeling is crucial for both optimal design and operating. Different types of heat exchangers are applied in (HENs) but plate heat exchangers (PHEs) are currently one of the most popular devices due to their high heat exchange efficiency combined with compact design, in comparison to conventional shell and tube heat exchangers [12,13]. This efficiency results from specific combination of gasketed corrugated metal plates that are pressed together to form the whole unit. The pattern of plates perforation and the internal flow configuration determine the performance of a particular PHE but at the same time, these construction details require complex and nonlinear modeling with a relatively large number of physical parameters that must be determined based on PHE assemblage. Consequently, every model should be casedependent, which limits its generality. Various works dealt with the subject of PHEs modeling based on their construction details. Starting from simple thermal modeling [14,15], interesting new works presenting different trends were also published. Thermal steady-state modeling of PHEs accounting for different plate and flow configurations was presented by Gut and Pinto [16]. The dynamical modeling where thermal dynamical equations describing PHE were coupled with material balance for milk fouling description can be found in [17]. 3D modeling and simulation of PHE based on computational fluid dynamics (CFD) was presented by Sammeta et al. [18]. The complexity of detailed modeling of heat transfer in PHE channels can be found in [19] where authors provided semi-empirical relationship describing heat exchange phenomenon accounting for channels geometry, flow velocity and physical fluid properties. These modeling approaches are very useful for optimal design and operating of new PHEs but their complexity is rather high, even if they assumed some simplifications. Usually, detailed case-dependent considerations were limited to steady-state conditions. At the same time, dynamical models were derived based on simplified assumptions (constant heat transfer coefficient, no heat loss, etc.) but they still accounted for PHE construction details. Consequently, they are useless for optimal design of high-layer control of existing HENs and for deriving regulatory-layer model-based control of single PHE operating in an existing process. For the latter, sometimes, the construction details for modeled PHE can be unknown or even if they are known, their inclusion in the model requires very hard modeling effort and finally leads to complex form. Thus, the motivation for this work was to suggest the general methodology, which provides the possibility of simplified input–output PHE modeling based only on input– output measurement data, without accounting for construction details. Additionally, it was assumed that the suggested method should be able to provide relatively accurate PHE modeling based on limited amount of measurement data. It is extremely important when modeling is carried out in industrial conditions where the measurements are collected from operating process and it is impossible to apply any extensive process excitation, especially leading to significant changes in operating conditions. In such cases, the simplified modeling is required and one possibility is to apply simplified dynamical description in the form of the FOPDT (First Order Plus Dead Time) Laplace transfer functions, whose parameters can be easily identified from the process step response. This approach is very popular for approximating the significant dynamical properties of the industrial SISO (Single Input Single Output) processes [20]. However, this description is limited only for a single (control) forward path and its parameters must be time varying to represent the process nonlinearities resulting from changes in operating point. Moreover, it does not include the in-

881

fluence of disturbances, even if they are measurable. Thus, the general (case-independent) PHE modeling of a relatively simple form preserving its most important dynamical properties (such as distributed parameter nature, nonlinearities, etc.) was suggested. It is based on first principle modeling approach that is usually more accurate for wide variations of operating conditions [21,22]. At the same time, the number of model parameters should be minimized and the reliable method for identification of their values should be ensured based only on (limited) measurement data collected from a particular system. Such a model should also incorporate the influence of measured disturbances and describe the internal couplings between process inputs and outputs. In this paper, such an approach to PHEs modeling is suggested based on conventional and widely accepted first principle modeling in the form of the partial differential equation derived from heat and mass conservation considerations. This approach is known for decades, not only for heat exchangers [23,24] but also for other dynamical systems, e.g. for distributed parameter pH processes [25], distributed parameter bioreactors [26], etc. It was widely applied in more or less complex and detailed forms for modeling the dynamics of double pipe heat exchangers [27], of helical baffle heat exchangers [28], etc. The form of this model is rather general and it describes the distributed parameter nature of heat exchange process and its nonlinearities. The concept suggested in this paper is based on this framework and it stands as the significant extension of the work of Fratczak et al. [29]. Namely, the very general and simplified form of the model was applied with only four model parameters lumping all unknown and case-dependent physical process parameters. Their values were identified from measurement data to ensure possibly highest input–output PHE modeling accuracy for wide range of operating conditions. Comparing to previous works of the authors, this paper addresses detailed model validation, sensitivity analysis and both structural and numerical stability considerations. 2. Simplified modeling of PHEs The concept suggested in this paper for simplified modeling of PHE is based on the assumption that the unit is perfectly insulated and operated in conventional setup presented schematically in Fig. 1. The unit operates in the counter-flow mode with two circuits: one for hot water and the second for cold water. In both circuits, it is assumed that measurement data of flow rates F1, F2 [m3/ s], inlet temperatures Tin1, Tin2 [°C] and outlet temperatures Tout1, Tout2 [°C] are accessible on-line at certain sampling time τs [s]. Fig. 2 shows simplified schematic flow diagram for example symmetric PHE with parallel configuration with five plates p1-p5 of length L [m]. Liquids flows in spaces (channels) between the plates formed by their perforation and gaskets between them. Heat ex-

Fig. 1. Conventional application diagram of heat exchange system.

882

M. Fratczak et al./Applied Thermal Engineering 98 (2016) 880–893

Fig. 2. Schematic flow diagrams for PHE in parallel configuration.

change takes place through the plates p2-p4 that form three contact surfaces for hot and cold water. For considered parallel configuration, counter current flow takes place through each plate. Configuration presented in Fig. 2 could be considered as symmetric two-channel unit. Its effective length and substitute cross sections of the channels depend on construction details (namely, on inner dimensions of the plates and their perforation), especially that fluid streams are separated between certain plates. Both liquids exchange the heat through the plates p2-p4 and in simplification (neglecting ports between the plates), these plates can be considered as walls separating the channels and substitute dimensions of these walls depend on construction details. Suggested modeling approach is inspired by the fact that the simplified substitute flow diagram of example PHE shown in Fig. 2 can represent the symmetric double-channel counter-flow heat exchanger and it can be modeled based on conventional doublepipe heat exchanger description, under the following assumptions and simplifications:

• • • •

• •

unit is perfectly insulated, parallel flow configuration, plug flow conditions are assumed in each channel (namely, perfect radial mixing is assumed, so the suggested model can assume one dimension flow), thermal capacity of the wall separating the channels is neglected (this assumption can be questionable for PHEs where thermal capacity of plates can be significant compared to thermal capacity of liquid residing within the unit; however, this assumption allows to simplify the model and it is justified by the fact that measurement data for temperature of plates is not accessible), no heat exchange in flow direction (again, one dimension flow is assumed), fluids are incompressible and single phase, their rheology does not vary significantly.

Boundary conditions are defined by inlet temperatures T1(0,t) = Tin1(t), T2(1,t) = Tin2(t) while initial profiles are described as T1(z,0) = T1,0(z), T2(z,0) = T2,0(z). The outputs are defined by outlet temperatures T1(1,t) = Tout1(t), T2(0,t) = Tout2(t). If the model given by Eqs. (1) was to be based on strict heat conservation considerations, it should include a number of physical parameters (such as overall heat transfer coefficient, liquid densities, specific heats, etc.) and of geometrical parameters (such as substitute length L, substitute cross section for the channels, etc.). Both groups of parameters represent case-dependent construction details, and direct correlation between their values, geometrical dimensions of PHE and its material specification is difficult to determine in general. For instance, substitute length L and channels cross section are difficult to determine even for simplest considered flow configuration, not even saying about more complex cases with higher number of plates or different flow configuration. For the considered example, substitute length surely depends on plates perforation so it cannot be determined as triple length of the plates. At the same time, overall heat transfer coefficient depends strictly on plate material specification and its construction (perforation). Consequently, the value of this coefficient should be the lumped approximation of convective heat transfer coefficient, plate thermal conductivity and its thickness. This value can also vary due to potential liquid fouling. Liquid densities and specific heats can also vary slightly, subject to liquid composition. Thus, in the suggested model given by Eq. (1), all these unknown or uncertain physical and construction parameters are lumped and represented by substitute positive model parameters p1 [1/m3], p2 [1/m3], h1 [1/s], h2 [1/ s]. The highest possible modeling accuracy is expected if their values are identified based on the input–output measurement data collected operating process, which stands as necessary stage of tuning the suggested simplified model. 3. Space discretization for numerical integration In this paper, Orthogonal Collocation Method (OCM) [30] is applied for spatial discretization of PHE simplified model given by Eq. (1). This discretization which transforms PDE description into the corresponding approximating Ordinary Differential Equation (ODE) system. This method provides high approximation accuracy for a relatively small number of (N+1) discretization points denoted as zi (for i = 0..N) and unevenly located along the normalized substitute length z ∈[0, 1]. Consequently, small number of discretization points results in relatively low dimension of approximating state vector x(t). The discretization points for OCM are chosen as two boundary points (z0 = 0 and zN = 1) and M = N-1 internal points determined as the roots of general orthogonal Jacobi polynomial:

PM (z ) = γ M z M − γ M −1z M −1 + γ M −2z M −2 − … + ( −1)

M

(2)

with coefficients defined by the following recursive formula:

M − p M +α + β + p +1 γp β + p +1 p +1

After introducing normalized space variable z = x/L [–], the simplified PHE model is suggested based on set of the Partial Differential Equations (PDE) describing dynamical evolution of temperatures T1(z,t), T2(z,t) that respectively denote the distribution of temperatures of hot and of cold water along effective substitute length of both channels. The form of this model is inspired by conventional framework based on the heat conservation considerations suggested for double-pipe counter current heat exchangers [23,24]:

where p = 0 .. M-1, γ0 = 1 and α > -1, β > -1. The choice of the parameters α, β determines the location of internal discretization points. In OCM, the Lagrange interpolation polynomial is applied for approximating the space derivatives of both temperatures T1(z,t), T2(z,t) at discretization point zi in the following way:

∂ T 1 (z , t ) ∂ T (z , t ) = − p 1F1 (t ) 1 − h1 (T1 (z , t ) − T 2 (z , t )) , ∂t ∂z

(1a)

∂ T 1 (z , t ) ≈ LT (z i )Tin 1 (t ) , ∂z zi

(4a)

∂ T 2 (z , t ) ∂ T (z , t ) = p 2F2 (t ) 2 + h2 (T1 (z , t ) − T 2 (z , t )) , ∂t ∂z

(1b)

∂ T 2 (z , t ) ≈ LT (z i )Tin 2 (t ) , ∂z zi

(4b)

γ p +1 =

(3)

M. Fratczak et al./Applied Thermal Engineering 98 (2016) 880–893 T

883

where Lˆ (z i ) = ⎡⎣Lˆ0 (z i ) , Lˆ1 (z i ) , … , LˆN (z i )⎤⎦ is (N+1)-element vector of the Lagrange polynomial components defined at the discretization point zi as:

Once OCM has been applied for spatial discretization of PHE model given by Eq. (1), the resulting approximating ODE system can be written in the following form:

⎛ N (z − z p ) ⎞⎟ ⎜∏ p =0 d ⎜ p≠ j ⎟ Lˆ j (z i ) = ⎟ , dz ⎜⎜ N (z j − z p ) ⎟ ∏ ⎜ p =0 ⎟ ⎝ p≠ j ⎠ zi

dx (t ) = ( A (t ) + H1 ) x (t ) + (B (t ) + H 2 )Tin (t ) dt (5)

where x (t ) = [T1 (z 1, t ) , … , T1 (z N , t ) , T 2 (z 0, t ) , … , T 2 (z N −1, t )]T denotes (2*N) state vector of approximating ODE system and T Tin (t ) = [Tin1 (t )Tin 2 (t )] includes inlet temperatures (boundary conditions). Time-varying matrices A(t) and B(t) can be 0 ⎡ −F1 (t ) p 1A + ⎤ , written in the block form as A (t ) = ⎢ 0 F2 (t ) p 2 A − ⎥⎦ ⎣ + 0 ⎡ −F1 (t ) p1B ⎤ B (t ) = ⎢ with time-invariant matrices A+, A-, 0 F2 (t ) p 2B − ⎥⎦ ⎣

T Tin1 (t ) = [Tin1 (t ) , T1 (z 1, t ) , … T1 (z N , t )] is (N+1)-element vector containing the discretized spatial profile of temperature of hot water T 1 (z,t) completed with the inlet temperature T in1 , and T Tin 2 (t ) = [T 2 (z 0, t ) , T 2 (z 1, t ) , … Tin 2 (t )] is (N+1)-element vector containing the discretized spatial profile of temperature of cold water T2(z,t) completed with the inlet temperature Tin2.

⎡ Lˆ1 (z 1 ) Lˆ2 (z 1 ) Lˆ3 (z 1 ) ⎢ Lˆ2 (z 2 ) Lˆ3 (z 2 ) ⎢ Lˆ1 (z 2 ) ⎢ ˆ Lˆ2 (z 3 ) Lˆ3 (z 3 ) ⎢ L 1 (z 3 ) ⎢    ⎢ A + = ⎢ Lˆ1 (z i −1 ) Lˆ2 (z i −1 ) Lˆ3 (z i −1 ) ⎢ ⎢ Lˆ1 (z i ) Lˆ2 (z i ) Lˆ3 (z i ) ⎢   ⎢  ⎢ˆ ˆ ˆ ⎢L1 (z N −1 ) L 2 (z N −1 ) L 3 (z N −1 ) ⎢ Lˆ (z ) Lˆ (z ) Lˆ (z ) 2 N 3 N ⎣ 1 N

 

Lˆi −1 (z 1 ) Lˆ (z ) i −1

2

Lˆi −1 (z 3 )   Lˆi −1 (z i −1 )  Lˆ (z )

 

i −1



i

Lˆi (z 1 ) Lˆ (z )



B+, B- defined as:

LˆN −1 (z 1 ) Lˆ (z )

 i 2 N −1 2 Lˆi (z 3 )  LˆN −1 (z 3 )    Lˆi (z i −1 )  LˆN −1 (z i −1 ) Lˆ (z )  Lˆ (z ) i



N −1

i





i



 Lˆi −1 (z N −1 ) Lˆi (z N −1 )  LˆN −1 (z N −1 )  Lˆi −1 (z N ) Lˆi (z N )  LˆN −1 (z N )

⎡ Lˆ0 (z 0 ) Lˆ1 (z 0 ) Lˆ2 (z 0 ) ⎢ Lˆ1 (z 1 ) Lˆ2 (z 1 ) ⎢ Lˆ0 (z 1 ) ⎢ ˆ Lˆ1 (z 2 ) Lˆ2 (z 2 ) ⎢ L 0 (z 2 ) ⎢    ⎢ − ˆ ˆ ˆ A = ⎢ L 0 (z i ) L 1 (z i ) L 2 (z i ) ⎢ ˆ ˆ ˆ ⎢ L 0 (z i +1 ) L1 (z i +1 ) L 2 (z i +1 ) ⎢    ⎢ ⎢ˆ ˆ ˆ ⎢L 0 (z N −2 ) L1 (z N −2 ) L 2 (z N −2 ) ⎢ Lˆ (z ) Lˆ (z ) Lˆ (z ) 1 N −1 2 N −1 ⎣ 0 N −1

    

Lˆi (z 0 ) Lˆ (z )

Lˆi +1 (z 0 ) Lˆ (z )



Lˆi (z 2 )  ˆ L (z )

Lˆi +1 (z 2 )  ˆ L (z )

 

i

i

1

i

i +1



1

i +1



i

 Lˆi (z i +1 ) Lˆi +1 (z i +1 )    ˆ ˆ  L i (z N −2 ) L i +1 (z N −2 )  Lˆ (z ) Lˆ (z ) i

N −1

i +1

N −1

LˆN −2 (z 0 ) Lˆ (z ) N −2

1

LˆN −2 (z 2 )  ˆ L (z ) N −2

i

LˆN −2 (z i +1 )  ˆ  L N − 2 (z N − 2 )  Lˆ (z )  

N −2

N −1

LˆN (z 1 ) ⎤ ⎥ LˆN (z 2 ) ⎥ ⎥ LˆN (z 3 ) ⎥ ⎥  ⎥ LˆN (z i −1 ) ⎥ , ⎥ LˆN (z i ) ⎥ ⎥  ⎥ ⎥ LˆN (z N −1 )⎥ LˆN (z N ) ⎦⎥ ⎡ Lˆ0 (z 1 ) ⎤ ⎡ LˆN (z 0 ) ⎤ LˆN −1 (z 0 ) ⎤ ⎢ ⎢ ⎥ ⎥ ⎥ ˆ ⎢ L 0 (z 2 ) ⎥ ⎢ LˆN (z 1 ) ⎥ LˆN −1 (z 1 ) ⎥ ⎢ ⎢ ⎥ ⎥ ⎥ ⎢ Lˆ0 (z 3 ) ⎥ ⎢ LˆN (z 2 ) ⎥ LˆN −1 (z 2 ) ⎥ ⎢ ⎢ ⎥ ⎥ ⎥    ⎢ ⎥ ⎢ ⎥ ⎥ ⎢ ⎥ ⎢ ⎥ + − LˆN −1 (z i ) ⎥ , B = ⎢ Lˆ0 (z i −1 ) ⎥ , B = ⎢ LˆN (z i −2 ) ⎥ . ⎥ ⎢ Lˆ0 (z i ) ⎥ ⎢ LˆN (z i −1 ) ⎥ LˆN −1 (z i +1 ) ⎥ ⎢ ⎥ ⎢ ⎥ ⎥   ⎢ ⎥ ⎢ ⎥  ⎥ ⎢ˆ ⎥ ⎢ˆ ⎥ ⎥ ˆ L N −1 (z N −2 )⎥ ⎢L 0 (z N −1 )⎥ ⎢L N (z N −2 )⎥ ⎢ ˆ ⎥ ⎢ˆ ⎥ LˆN −1 (z N −1 ) ⎥⎦ ⎢⎣ L N (z N −1 ) ⎥⎦ ⎣⎢ L 0 (z N ) ⎥⎦

Time-invariant matrices H1 and H2 consist of model parameters h1 and h2 and are defined in the following way:

0 ⎡ −h1 0 ⎢ 0 −h 0 1 ⎢ ⎢ 0 0 −h1 ⎢   ⎢  ⎢ 0 0 0 ⎢ 0 0 ⎢ 0 ⎢    ⎢ ⎢ 0 0 0 ⎢ 0 0 0 H1 = ⎢ ⎢ 0 0 0 ⎢ 0 0 ⎢ h2 ⎢ 0 0 h2 ⎢ ⎢    ⎢ 0 0 ⎢ 0 ⎢ 0 0 0 ⎢   ⎢  ⎢ 0 0 0 ⎢ ⎢⎣ 0 0 0

0  0 0  0 0  0     −h1 0  0 −h1     0 0  0 0 0  0 0  0  0 0     0 h2  0 0     0 0  0 0

0 0  0 0 0  0 0 0  0      0 0 0  0 0 0      −h1 0 0  0 −h1 0 0 −h2  0 0 0  0  0 0 0      0 0 0  0 0 0      0 0 0  h2 0 0

h1 0 0  0 0  0 0 0 −h2 0  0 0  0 0

0 h1 0  0 0  0 0 0 0 −h2  0 0  0 0

(6)

 0  0  0    0  0    0  0  0  0  0    −h2  0    0  0

0 0 0  h1 0  0 0 0 0 0  0 −h2  0 0

 0  0  0    0  0    0  0  0  0  0    0  0    −h2  0

0 ⎤ ⎡0 0⎤ ⎢0 0⎥ 0 ⎥ ⎢ ⎥ ⎥ ⎢0 0⎥ 0 ⎥ ⎢ ⎥ ⎥  ⎥ ⎥ ⎢ ⎢0 0⎥ 0 ⎥ ⎢ ⎥ ⎥ 0 ⎥ ⎢0 0⎥ ⎢ ⎥   ⎥ ⎢ ⎥ ⎥ ⎢0 0⎥ h1 ⎥ ⎢ ⎥ ⎥ 0 ⎥ 0 h1 ⎥ , H2 = ⎢ . ⎢h2 0 ⎥ 0 ⎥ ⎢ ⎥ ⎥ 0 ⎥ ⎢0 0⎥ ⎢0 0⎥ 0 ⎥ ⎢ ⎥ ⎥ ⎢  ⎥ ⎥ ⎢ ⎥ ⎥ 0 ⎥ ⎢0 0⎥ ⎢0 0⎥ 0 ⎥ ⎢ ⎥ ⎥  ⎥ ⎥ ⎢ ⎢0 0⎥ 0 ⎥ ⎢ ⎥ ⎥ ⎢⎣ 0 0 ⎥⎦ −h2 ⎥⎦

884

M. Fratczak et al./Applied Thermal Engineering 98 (2016) 880–893

19 [°C]). Figs. 4–6 show the variations of PHE inputs collected during the aforementioned experiments with sampling time τS = 1 [s] and in the presence of denoted step changes of the adjustable process variables. It should be also noted that due to configuration of experimental setup shown in Fig. 3, step changes of the set-point F1SP cause both the variations of flow rate F1 and also the variations of inlet temperature Tin1 (Ph = const during this experiment), which is presented in Fig. 4.Concept of the tuning procedure is presented in Fig. 7. After adjusting the initial values of model parameters, input measurements presented in Figs. 4–6 were separately applied to discretized model (6) to compute modeled outputs. These outputs were compared with corresponding measurement data of PHE outputs. For each excitation signal, the modeling accuracy of the suggested simplified model was quantified by mean square modeling error MSME [°C] defined as:

Then, this approximating ODE system given by Eq. (6) can be solved by any numerical integration method, e.g. the simplest Euler method or more complex Runge–Kutta method, depending on required accuracy and acceptable computational complexity. 4. Tuning of simplified PHE model As it was said, the suggested simplified PHE model given by Eq. (1) requires tuning that consists in adjusting the values of substitute model parameters p1, p2, h1, h2. This tuning was completed based on measurement data of the PHE inputs (Tin1, Tin2, F1, F2) and its corresponding outputs (Tout1, Tout2) collected during the experiments carried out for the laboratory setup presented schematically in Fig. 3. It represents simple heat exchange and distribution system with heat source (electric flow heater) and heat receiver (PHE) with varying heat consumption set by adjusted flow rate F2. In this system, PHE LM25-6 made by TAU Energy operates in application presented in Fig. 1 and it consists of six plates asymmetrically forming two channels: for hot water of the volume of 0.19 [L] and for cold water of the volume of 0.29 [L]. Temperature measurements (Tin1, Tin2, Tout1, Tout2) are provided by platinum-wire resistance temperature devices (RTD) directly-coupled to 0–10 [V] transmitters. Flow rates F1, F2 are measured by magnetic flow sensors SIEMENS© Sitrans fm mag 1100 and can be adjusted respectively by setting the setpoints F1SP, F2SP for PI-based flow controllers operating the control valves V1, V2. Hot water is supplied by electric flow heater with power supply Ph adjusted within the percentage range 0–100 [%] of the nominal power Pnom = 4 [kW] by thyristor-based unit controlled by PWM (Pulse Width Modulation) algorithm. All process signals are connected to analog I/O modules of National Instruments Compact Fieldpoint cFP2020 that is interfaced to PC computer with LabVIEW-based SCADA system through Ethernet protocol. The tuning procedure was based on possibly most realistic experiment scenarios with limited process excitation. It was decided to apply three separated step changes of the process variables F1SP, F2SP and Ph and use the corresponding measured variations of PHE inputs (Tin1, Tin2, F1, F2) and outputs (Tout1, Tout2) as three different sets of measurement data for model tuning. Readers should note that some of PHE inputs (F1, F2, Tin1) were adjusted indirectly by respective changes of F1SP, F2SP and Ph while one input Tin2 was only measured and could not be adjusted (municipal water supply of approximately constant temperature slightly varying within the range of 17–

PI

MSME =

1

τ max

τ max

∫ e (t )

T

Qe (t )dt ,

(7)

0

where e (t ) = y (t ) − y m (t ) , y (t ) = [Tout 1 (t )Tout 2 (t )] represents meaT m m sured output temperatures, y m (t ) = [Tout 1 (t )Tout 2 (t )] contains the respective temperatures computed by the model (6), Q = diag(0.4, 0.6) is weighting matrix and τmax denotes the time of experiment. In the considered case, matrix Q was intentionally defined as asymmetric because outlet temperature Tout2(t) is considered as the major system output (e.g. for control purposes) and thus the accuracy of its modeling is more important. The tuning procedure was based on minimization of MSME value defined by Eq. (7) as the objective (quality factor) depending on set of decision variables. The minimization procedure was numerically solved by Nelder and Mead Simplex algorithm [31,32] that provides very good numerical properties, even for functions with discontinuities or for which the gradient-based search fails. Additionally, this method is implemented in popular computing environments, both commercial and non-commercial, e.g. Matlab, Octave, etc., which is important practical advantage of suggested approach. Generally, for OCM discretization method, set of decision variables includes substitute model parameters p1, p2, h1, h2, a number of discretization points N and OCM parameters α > -1, β > -1 that determine the location of internal discretization points computed by formula (2). However, in such a case, numerical complexity of the minimization problem is relatively high. At the same time, the T

F1SP Tin1

V1 F

Tout2

T

F2 T

F

F1

Ph

electric flow heater

heat exchanger

municipal water supply

Tout1

PI

F2SP municipal water supply

T T

Tin2

V2

Fig. 3. Schematic diagram of heat exchange and distribution laboratory setup.

M. Fratczak et al./Applied Thermal Engineering 98 (2016) 880–893

885

Fig. 4. Variations of flow rate F1 and corresponding variations of inlet temperature Tin1 in the presence of denoted step changes of F1SP.

local minima are more likely, which is a difficulty when the optimization problem is solved numerically. Thus, it was decided to simplify computational complexity of optimization method. Namely, the general formula (2) was considered only on three widely accepted cases that represent Chebyshev polynomials of 1st and of 2nd kind (respectively denoted as Che1 with α = β = 0.5 and Che2 with α = β = -0.5) and Legendre polynomial (denoted as Leg with α = β = 0). This choice is justified by the fact that these polynomials represent the optimal sensors locations for numerical integration and for design of OCM discretization-based state observers for distributed parameter systems [26,33,34]. At the same time, for chosen locations, the internal discretization points are always distributed symmetrically around z = 0.5, which is a very desired property for spatial discretization of PHE model for counter current flow. Once the values of the parameters α, β are fixed, the set of decision variables is limited but still there is a difficulty with optimal choice of a number of discretization points N. The parameters p1, p2, h1, h2 are positive and real while N is natural. Thus, it was decided to solve simplified minimization problem by successive solving for min (MSME ) , for N chosen within the range N = 2..10 to deterp1 ,p 2 ,h1 ,h2

mine the compromise between the number of the discretization points and corresponding MSME value. Starting from N = 2, mini-

mization problem

min (MSME ) was solved for each location of

p1 ,p 2 ,h1 ,h2

discretization points determined by Che1, Che2, Leg forms of the formula (2), using three excitation signals shown in Figs. 4–6. The results of corresponding MSME values are presented in Fig. 8. First, it should be said that for each considered case, the corresponding MSME values were higher from the outputs variability resulting from measurement noise, which was quantified by noise variance σ = 0.03 [°C] computed for example operating steady states. Thus, the computed MSME values really represent the impact of measurement inaccuracy subject to different OCM parametrization. At the same time, computed noise variance allows for assessment if differences between considered MSME values are statistically significant. When tuning was based on variations of flow rate F1, the best modeling accuracy was obtained for N = 4, practically regardless of the choice between Che1, Che2, Leg. The similar case takes place for the excitation by flow rate F2. The choice of N for exciting by temperature Tin1 is not so clear but still N = 4 seems to be acceptable. Additionally, the results presented in Fig. 8 show that the choice between Che1, Che2 and Leg is also not clear because the differences between corresponding MSME values are rather insignificant. Namely, the level of these differences is comparable with com-

886

M. Fratczak et al./Applied Thermal Engineering 98 (2016) 880–893

Fig. 5. Variations of flow rate F2 in the presence of denoted step changes of F2SP.

puted noise variance σ = 0.03 [°C] and can vary, subject to impact of the measurement noise in data used for tuning procedure. Thus, it was finally decided to continue further tuning based on location of (N+1) = 5 discretization points with application of Che2. However, this choice still corresponds to different values of PHE model parameters, subject to excitation signal used for tuning, which is shown in Table 1. Thus, for each defined parameters set SF1, SF2, STin1, OCM-discretized model (6) was excited once again by the same input signals shown in Figs. 4–6. Then, for each case, the respective MSME values were computed and they are presented in Fig. 9 with additional bar representing the mean MSME value calculated separately for each parameters set. Once again, differences in modeling accuracy are not very significant comparing to impact of measurement noise. Thus, the final choice was made for SF1 (bolded in Table 1) because it provides highest process excitation (simultaneous variations of F1 and Tin1, see Fig. 4). The modelling accuracy for this choice is presented in Figs. 10–12 in comparison with the measured PHE output signals Tout1, Tout2. Additional comment should be made on values of model parameters shown in Table 1. For each set, p1 ≠ p2 and h1 ≠ h2 appear

while their (approximate) equality was expected assuming symmetric example PHE shown in Fig. 2. However, readers should note that (I) PHE used in experimental setup is not symmetric and (II) the model (6) was derived based on assumptions that are not very realistic. For example, the unit is not perfectly insulated and the heat capacity of plates (including p1 and p5, see Fig. 2) was neglected. Consequently, tuning procedure adjusted the model parameters to ensure a compromise between possibly highest modeling accuracy and compensation for modeling simplifications. Fig. 9 also shows the influence of process operating conditions on the tuning procedure. Fig. 4–6 show the excitation signals used for tuning and the lowest MSME values were obtained for each parameters set SF1, SF2, STin1, which results from the fact that exciting variations of Tin1 are very smooth. Consequently, the highest MSME values can be observed for each case for exciting by F2, whose variations are oscillatory due to faulty control valve V2 (see, Fig. 3). However, it does not question the tuning results – tuning procedure is still reliable but as in every experimental case, the results suffer from any faults and quality of measurement data.

Fig. 6. Variations of inlet temperature Tin1 in the presence of denoted step changes of Ph.

M. Fratczak et al./Applied Thermal Engineering 98 (2016) 880–893

crement of 1%. Then, for each case, respective MSME value was computed and the results are presented in Fig. 14. This figure shows the relationship between percentage variation of each model parameter and corresponding percentage increment of MSME value (calculated with respect to its nominal value computed for experiment shown in Fig. 13). These results lead to following important conclusions:

input measurements Tin1, Tin2, F1, F2

•

PHE model

experimental setup

•

modelled outputs ym = [Tmout1, Tmout2]

887

output measurements y = [Tout1, Tout2]

_ + modelling error e Fig. 7. Data flow diagram of experimental tuning procedure.

5. Final validation of simplified PHE model Tuning of OCM-based simplified PHE model (6) was carried out based on a relatively short series of measurement data of inputs shown in Figs. 4–6. However, the extensive validation was carried out during operating the laboratory setup shown in Fig. 3 to simulate realistic conditions and to provide possibly wide variations of operating points. During numerous experiments, the model (6) with suggested values of model parameters SF1 was computed online and excited with measured inputs F1, F2, Tin1, Tin2 and computed outlet temperatures Tout1, Tout2 were compared with measurement data. The results were satisfying, which can be seen for example experiment run presented in Fig. 13. Corresponding MSME value computed over the whole experiment is 0.36 and the results show that modeling accuracy is acceptable from practical viewpoint, both for PHE dynamics and its statics. Stability of OCM-based approximating model (6) and of its numerical integration procedure are discussed in details in Appendix A, jointly with analysis of dynamical properties. After tuning, apart from final validation for different operating conditions, each model should be analyzed in terms of sensitivity to variations of its parameters. Formally, such a sensitivity analysis should be carried out analytically by separate calculating partial derivatives of model output with respect to model parameters. However, in the considered case, approximating OCM model (6) has two outputs (outlet temperatures Tout1, Tout2) and analytical approach is very difficult due to complexity of OCM approximation. Thus, instead it was decided to apply numerical sensitivity analysis. This approach consists in successive validation of the suggested model (6) based on the same measurement data presented in Fig. 13 but with model parameters changed around their nominal values SF1 bolded in Table 1. The parameters changes Δp1, Δp2, Δh1, Δh2 were applied separately within normalized range of [−10%, 10%] with in-

•

even if the suggested set of model parameters SF1 was chosen by tuning based only on short series of measurement data shown in Figs. 4–6, this set is still optimal for more realistic example experiment presented in Fig. 13 because SF1 provides approximately the minimal MSME value, even if percentage increment of MSME value is slightly asymmetric for Δp2 and Δh2, approximated sensitivity functions have quadratic and modeling accuracy deteriorates significantly with variation of each model parameter, it can be seen that the model (6) is more sensitive to variations of parameters p2 and h2, which probably results from asymmetric choice of weighting matrix Q defined for MSME computed by Eq. (7).

Readers should also note that both tuning procedure and model validation were carried out based on real measurement data that is always noisy. Thus, the nominal values of model parameters can vary slightly subject to impact of this random noise. However, presented sensitivity analysis allows for determining accuracy of parameters adjustment to ensure desired modeling accuracy quantified by acceptable percentage MSME increment. Suggested OCMbased simplified PHE modeling was additionally compared with more conventional approach that consists in spatial discretization of dynamical model (1) the Finite Difference Method (FDM) [35]. For this method, the approximations of space derivatives of temperature profiles are calculated for (N+1) evenly located discretization point z i as:

∂T1 (z , t ) ∂z zi

≈ T1 (z i , t )−ΔTz1 (z i −1, t ) ,

∂T 2 (z , t ) ∂z zi

≈ T2 (z i +1, tΔ)z−T2 (z i , t ) , where

Δz = 1/N denotes constant space discretization instant. This approximation leads to the same general form of discrete-space model (6) but matrices A+, A-, B+, B- have different (simpler) forms. Two cases were considered: high number of discretization points (N+1) = 100 that is suggested for FDM spatial approximation and the same number (N+1) = 5 that was adjusted for OCM-based modeling. For both FDM-based cases, the same tuning procedure described in Section 4 was separately applied and corresponding modeling accuracy was compared with accuracy of suggested OCM-based PHE modeling. Equivalent modeling accuracy was obtained for (N+1) = 100 (MSME = 0.36 for example experiment shown in Fig. 13). For the case with (N+1) = 5, modeling accuracy of FDM-based approximation was significantly lower compared to suggested OCM-based approach. Even if both discretized models were separately tuned and computed on-line with optimal values of their model parameters p1, p2, h1, h2, MSME = 0.42 was obtained for FDM-based approximation for example experiment run shown in Fig. 13, which is approximately 16% increment comparing to OCM-based modeling. 6. Potential applicability and extension to other cases Due to its relative simplicity resulting from high modeling accuracy obtained for a low number of spatial discretization points, the suggested modeling approach potentially has a wide range of applications:

•

it allows for dynamical simulation of complex HENs and consequently it supports their design, optimization and control, including virtual commissioning of industrial control systems, which is a very promising approach that allows for significant

M. Fratczak et al./Applied Thermal Engineering 98 (2016) 880–893

0.35

Che1 Che2 Leg

888

MSME for excitation signal F1 Che1 Che2 Leg

Che1 Che2 Leg

Che1 Che2 Leg

Che1 Che2 Leg

N=4

Che1 Che2 Leg

N=3

Che1 Che2 Leg

0.2

Che1 Che2 Leg

0.25

Che1 Che2 Leg

0.3

N=5

N=6

N=7

N=8

N=9

N=10

0.15 0.1 0.05 0

Che1 Che2 Leg

Che1 Che2 Leg

Che1 Che2 Leg

Che1 Che2 Leg

Che1 Che2 Leg

0.25

Che1 Che2 Leg

MSME for excitation signal F2 Che1 Che2 Leg

0.3

Che1 Che2 Leg

0.35

Che1 Che2 Leg

N=2

N=4

N=5

N=6

N=7

N=8

N=9

N=10

0.2 0.15 0.1 0.05 0 N=2

N=3

0.35

Che1 Che2 Leg

Che1 Che2 Leg

Che1 Che2 Leg

Che1 Che2 Leg

Che1 Che2 Leg

Che1 Che2 Leg

0.2 0.15

Che1 Che2 Leg

Che1 Che2 Leg

0.25

MSME for excitation signal Tin1 Che1 Che2 Leg

0.3

N=4

N=5

N=6

N=7

N=8

N=9

N=10

0.1

0.05 0 N=2

N=3

Fig. 8. Optimal MSME values for chosen range of N, for denoted three excitation signals and three different locations of discretization points determined by Che1, Che2, Leg.

0.25

MSME

Model parameters Excitation signal

p1 [1/m ]

p2 [1/m ]

h1 [1/s]

h2 [1/s]

Parameters set SF1 Optimal for F1 excitation (Fig. 4) Parameters set SF2 Optimal for F2 excitation (Fig. 5) Parameters set STin1 Optimal for Tin1 excitation (Fig. 6)

902.039

826.651

0.0901

0.0815

869.565

833.194

0.0869

0.0821

895.496

823.655

0.0894

0.0812

0.1 0.05 0

SF1

SF2

mean MSME; 0,17

3

F1 F2 Tin1

3

F1 F2 Tin1

0.15 Table 1 Three optimal sets of model parameters for three exciting signals for OCM discretization with (N+1) = 5 and Che2.

mean MSME; 0,18

0.2

F1 F2 Tin1

•

Finally, it is very likely that suggested modeling approach can be extended for modeling of any liquid-to-liquid plate heat exchangers, without accounting for their construction details and operating mode. The optimization-based tuning procedure is

mean MSME; 0,17

•

reduction of real commissioning time of the PLC-based control systems and SCADA applications [36–38], it can be used as a basis for deriving nonlinear (adaptive) modelbased controllers according to general suggestions presented in [39–42] – this approach allows for very simple feedforwarding from measured disturbances due to multi-input multi-output form of the suggested model, such a modeling can be also applied as a basis for design of dynamical state observers that support the advanced control and SCADA systems [43–45].

STin1

Fig. 9. MSME values for three OCM-discretized model parameters sets with (N+1) = 5 and Che2.

M. Fratczak et al./Applied Thermal Engineering 98 (2016) 880–893

Fig. 10. Modeling accuracy for excitation signal F1 from Fig. 4 (N = 4, Che2, SF1).

Fig. 11. Modeling accuracy for excitation signal F2 from Fig. 5 (N = 4, Che2, SF1).

Fig. 12. Modeling accuracy for excitation signal Tin1 from Fig. 6 (N = 4, Che2, SF1).

889

890

M. Fratczak et al./Applied Thermal Engineering 98 (2016) 880–893

Fig. 13. Example modeling accuracy for final validation of OCM-based PHE model in the presence of denoted step changes of process variables (N = 4, Che2), MSME = 0.355.

reliable and it allows for relatively easy adjustment of values of the model parameters p1, p2, h1, h2 based only on limited measurement data collected from normally operating process without its extensive excitation. In the considered experimental case, time-invariant form of the suggested model given by Eq. (1) jointly with optimization-based tuning procedure ensured acceptable modeling accuracy. However, if case-dependent results are not satisfying, the suggested timeinvariant form can be easily rearranged into the time-varying form if by setting the model parameters p1, p2, h1, h2 as functions of operating conditions. Namely, influence of flow rates F1, F2 (defining process operating point), of physical parameters of liquids (if they vary significantly) or of ambient temperature Tamb (as process disturbance if it is measured on-line) can be included. This inclusion can be made by using appropriate formulas, e.g. p1 = fp1(liquids rheology), p2 = fp2(liquids rheology), h1 = fh1(F1, F2, Tamb, liquids rheology), h2 = fh2(F1, F2, Tamb, liquids rheology). Such formulas can be taken from literature but it always require case-dependent construction details. They can be also approximated from process measurement data but in such a case, much more extensive process excitation is required to collect a set of measurements representing steady-state of different operating points spread within whole operating region. Then,

p1

p2

suggested tuning procedure should be repeated separately for each operating point based on selected measurement data to obtain locally-optimal sets of model parameters. Finally, the aforementioned formulas should be approximated based on these optimal sets. However, readers should note that both approaches complicate PHE modeling and it should be decided if application of timeinvariant form of the suggested model can improve overall modeling accuracy. 7. Conclusions In this paper, simplified dynamical modeling of PHEs is suggested. The concept combines general form of first principle distributed parameter model of double-pipe heat exchanger with its OCM-based spatial discretization and optimization-based tuning procedure. Such an approximating PHE model was tuned based on limited measurement data and then successfully validated based on measurement data representing different operating points. Its stability and dynamical properties were also addressed. The results show that suggested modeling and tuning method is reliable and it can be useful for practical applications where simplified modeling is preferred.

h1

h2

percentage MSME increment

80 70 60 50 40 30 20 10 0

-10 -9

-8

-7

-6

-5

-4 -3 -2 -1 0 1 2 3 4 5 percentage variation of model parameter

6

Fig. 14. Results of sensitivity analysis of the suggested OCM-based model.

7

8

9

10

M. Fratczak et al./Applied Thermal Engineering 98 (2016) 880–893

891

Acknowledgements M. Fratczak and P. Nowak were supported under “DoktoRIS Scholarship for innovative Silesia”, co-financed by European Union under European Social Fund. J. Czeczot and M. Metzger were supported by Ministry of Science and Higher Education under grant BKUiUA. Calculations were done with the use of GeCONiI infrastructure (PO IG 02.03.01-24-099). Authors would like to thank anonymous reviewers for their help in preparing the final version of this manuscript. Appendix A: Dynamical properties of OCM-based approximation Potentially, analysis of dynamical properties of OCM-based spatially dicretized approximating model of PHE (6) can be carried out by two methods. First approach can be based on second Lyapunov method but in this case, only global stability can be proved and nothing can be said about any other dynamical properties of the suggested model. For deeper insight into these properties, more detailed analysis is required and it is based on Eq. (6) rewritten assuming constant flow rates F10 > 0, F20 > 0:

dx (t ) = ( A 0 + H1 ) x (t ) + (B0 + H 2 )Tin (t ) , dt

(A1)

0 ⎤ 0 ⎤ ⎡ −F10 p1A + ⎡ −F10 p1B + , B0 = ⎢ with the where A 0 = ⎢ −⎥ F20 p 2 A ⎦ F20 p 2B − ⎥⎦ 0 0 ⎣ ⎣ same time-invariant matrices A+, A-, B+, B- defined in Section 3. From Eq. (A1), the equilibrium point x0 can be calculated assuming T constant inlet temperatures (boundary conditions) Tin 0 = [Tin10Tin 20 ] as:

x 0 = − ( A 0 + H1 )

−1

(B0 + H 2 )Tin 0.

Fig. A2. Distribution of zeros and poles of the example transfer function T 2 (s ) . K 12 (s ) = Tout in 1 (s )

F20 varying within the range [8.33·10−6, 6.67·10−5] with the instant of 8.33·10−6. The analysis of Fig. A1 allows for drawing the following conclusions:

• •

(A2)

Eq. (A1) has the form of linear dynamic system with state matrix (A0 + H1) depending on operating point defined by flow rates F10, F20. Thus, these flow rates determine dynamical properties of the model (6). Assuming their range of variations (usually taken from technical specification of PHE installation setup), for desired flow rates F10, F20, the series of state matrices can be generated and the properties of these matrices can be investigated from their eigenvalues λ. The example is presented in Fig. A1, which shows the distribution of the eigenvalues for state matrix (A0 + H1) computed for F10,

•

OCM-based approximation (6) is stable for assumed variations of flow rates because real parts of all eigenvalues are strictly negative, approximating system (6) provides oscillatory behavior defined by tg(α) = b/a, which is a characteristic feature of OCM-based spatial discretization; however, the oscillations are not very significant, which can be also seen in Figs. 10-13 showing the validation of the model (6), the eigenvalue of largest negative real part (denoted in Fig. A1 as λm = −1.228) determines the smallest time constant of the dynamical approximating system (6), which is 0.814 [s]; consequently, this value determines the largest integration step that can be applied for numerical integration of the model (6) to ensure its numerical stability (e.g., according to the rule of thumb, for Euler first order integration method this step should not exceed 0.08 [s], which is ensured in this paper).

Surely, for higher flow rates, the value of λm becomes more negative so the acceptable integration step decreases but still OCMbased model (6) is stable. Additionally, based on Eq. (A1), Laplace analysis can be also provided to investigate chosen transfer functions. After Laplace transform, Eq. (A1) is rewritten as follows:

dx (s ) = ( A 0 + H1 ) x (s ) + (B0 + H 2 )Tin (s ) , dt

(A3)

and for output vector y (s ) = [Tout 1 (s )Tout 2 (s )] , the transfer function matrix K(s) can be calculated as: T

y (s ) = C (s I − ( A 0 + H1 )) (B0 + H 2 )Tin (s ) .   −1

(B4)

K(s )

Fig. A1. Example distribution of the eigenvalues for state matrix (A0 + H1).

Fig. A2 shows the distribution of zeros (circles) and poles (crosses) T 2 (s ) of the example transfer function K 12 (s ) = Tout for operating in 1 (s ) −5 point defined by flow rates F10 = 4.17·10 , F20 = 3.33·10−5. This distribution shows two important features of suggested OCM-based approximating model (6):

892

• •

M. Fratczak et al./Applied Thermal Engineering 98 (2016) 880–893

relative order of this transfer function is unitary, which is a big advantage when the model (6) is to be applied as a basis for model-based controller design, existence of positive zero determines the inverse response of the approximating model (6); however, this phenomenon is not very significant, which can be seen in Figs. 10-13 showing the validation of the suggested model.

This analysis of dynamical properties of the model (6) is simplified to the case when flow rates are constant. However, this approach can be extended for cases assuming variability of flow both rates but at the price of more complex analysis. Such an approach requires formal linearization of the model (6) but even then, state matrix (A0 + H1) is the same and defined by operating point determined by flow rates F10, F20 and consequently, the conclusions resulting from presented simplified analysis are still valid. Nomenclature A A+ AB B+ Be F FSP h H1 H2 i I K L Lˆ j (z i ) Lˆ (z i ) M MSME N+1 p Ph Pnom Q s t T Tin Tin Tout x y z

OCM-based state block matrix of the discretized model Upper-left component matrix of A Lower-right component matrix of A OCM-based input block matrix of the discretized model Upper-left component matrix of B Lower-right component matrix of B Vector of modeling error Flow rate [m3/s] Set-point for flow controller Substitute model parameter representing heat exchange [1/s] State matrix of the heat exchange parameters of OCMbased discretized model Input matrix of the heat exchange parameters of OCMbased discretized model Space discretization index Unitary matrix Transfer function matrix Substitute length [m] Lagrange polynomial component Vector of Lagrange polynomial components Number of internal space discretization points Mean square modeling error [°C] Number of space discretization points Substitute model parameter representing flow conditions [1/m3] Percentage of heater power supply range [%] Nominal power of heater [W] Diagonal weighting matrix Laplace transform parameter Time [s] Fluid temperature Inlet temperature [°C] Discretized spatial profile of temperature Outlet temperature [°C] State vector of discretized model Output vector of discretized model Normalized space variable [–]

Chosen abbreviations Che1 Chebyshev polynomial of 1st kind Che2 Chebyshev polynomial of 2nd kind FDM Finite difference method HEN Heat exchange network Leg Legendre polynomial OCM Orthogonal collocation method ODE Ordinary differential equation(s)

PDE PHE SF1 SF2 STin1

Partial differential equation(s) Plate heat exchanger Set of discretized model parameters optimal for excitation F1 Set of discretized model parameters optimal for excitation F2 Set of discretized model parameters optimal for excitation Tin1

Greek symbols α, β Parameters of orthogonal collocation method Δ Increment γ Coefficient of orthogonal Jacobi polynomial λ Eigenvalue Eigenvalue of the largest real part λm σ Measurement noise variance [°C] Experiment time [s] τmax Sampling time [s] τs Subscripts 0 Steady state 1 Hot fluid 2 Cold fluid Superscript m Measured output References [1] M. Bauer, I.K. Craig, Economic assessment of advanced process control – a survey and framework, J. Process Control 18 (2008) 2–18. [2] R.R. Rhinehart, M.L. Darby, H.L. Wade, Editorial – choosing advanced control, ISA Trans. 50 (2011) 2–10. [3] Z. Yuan, B. Chen, G. Sin, R. Gani, State-of-the-art of progress in the optimizationbased simultaneous design and control for chemical processes, AIChE J. 58 (6) (2012) 1640–1659. [4] I.C. Kemp, Pinch Analysis and Process Integration: A User’s Guide to Process Integration for the Efficient Use of Energy, Butterworth-Heinemann, Oxford, 2007. [5] J.C. Bonhivers, B. Srinivasan, P.R. Stuart, New analysis method to reduce the industrial energy requirements by heat-exchanger network retrofit: part 1 – concepts, Appl. Therm. Eng. (2014) http://dx.doi.org/10.1016/ j.applthermaleng.2014.04.078. [6] K. Semkov, E. Mooney, M. Connolly, C. Adley, Efficiency improvement through waste heat reduction, Appl. Therm. Eng. 70 (2014) 716–722. [7] P. Jie, N. Zhu, D. Li, Operation optimization of existing district heating systems, Appl. Therm. Eng. 78 (2015) 278–288. [8] C.A. Bennett, R.S. Kistler, T.G. Lestina, D.C. King, Improving heat exchanger designs, Chem. Eng. Prog. (2007) 40–45. [9] A.S. Morrison, M.J. Atkins, M.R.W. Walmsle, Ensuring cost-effective heat exchanger network design for non-continuous processes, Chem. Eng. Trans. 29 (2012) 295–300. [10] M. Bakošová, J. Oravec, Robust model predictive control for heat exchange network, Appl. Therm. Eng. 73 (2014) 922–928. [11] A. Isafiade, M. Bogataj, D. Fraser, Z. Kravanja, Optimal synthesis of heat exchanger networks for multi-period operations involving single and multiple utilities, Chem. Eng. Sci. 127 (2015) 175–188. [12] L. Wang, B. Sundén, R.M. Manglik, Plate Heat Exchangers. Design, Applications and Performance, WIT Press, Southampton, Boston, 2007. [13] S. Kakaç, H. Liu, A. Pramuanjaroenkij, Heat Exchangers: Selection, Rating, and Thermal Design, third ed., CRC Press Taylor & Francis Group, 2012. [14] A.A. McKillop, W.L. Dunkley, Plate heat exchangers: heat transfer, Ind. Eng. Chem. 52 (9) (1960) 740–744. [15] M. Masubuchi, A. Ito, Dynamic analysis of a plate heat exchanger system, Bull. JSME 20 (142) (1977) 434–441. [16] J.A.W. Gut, J.M. Pinto, Modeling of plate heat exchangers with generalized configurations, Int. J. Heat Mass Transf. 46 (2003) 2571–2585. [17] M.C. Geordiakis, S. Macchietto, Dynamic modelling and simulation of plate heat exchangers under milk fouling, Chem. Eng. Sci. 55 (2000) 1605–1619. [18] H. Sammeta, K. Ponnusamy, M.A. Majid, K. Dheenathayalan, Effectiveness charts for counter flow corrugated plate heat exchanger, Simul. Model. Pract. Theory 19 (2011) 777–784. [19] O.P. Arsenyeva, L.L. Tovashnyanskyy, P.O. Kapustenko, O.V. Demirskiy, Generalized semi-empirical correlation for heat transfer in channels of plate heat exchanger, Appl. Therm. Eng. 61 (2) (2014) 1208–1215. [20] J.D. Álvarez, L.J. Yebra, M. Berenguel, Repetitive control of tubular heat exchangers, J. Process Control 17 (2007) 689–701.

M. Fratczak et al./Applied Thermal Engineering 98 (2016) 880–893

[21] R.R. Rhinehart, Editorial viewpoint. The century’s greatest contributions to control practice, ISA Trans. 39 (2000) 3–13. [22] R. Murray-Smith, T.A. Johansen, Multiple Model Approaches to Modelling and Control, Taylor & Francis, 2010. [23] Y. Takahashi, in: A. Tustin (Ed.), Transfer Function Analysis of Heat Exchange Processes, Automation & Manual Control, Academic Press Inc., New York, 1952, p. 1952. [24] H.M. Paynther, Y. Takahashi, A new method of evaluating dynamic response of counterflow and parallel flow heat exchangers, Trans. ASME 78 (1956) 749–758. [25] K. Stebel, M. Metzger, Distributed parameter model for pH process including distributed continuous and discrete reactant feed, Comput. Chem. Eng. 38 (2012) 82–93. [26] N. Tali-Maamar, T. Damak, J.P. Babary, Application of the collocation method for simulation of distributed parameter bioreactors, Math. Comput. Simul. 35 (1993) 303–319. [27] M.A. Abgeldhani-Idrissi, F. Bagui, L. Estel, Analytical and experimental response time rate step along a counter flow double pipe heat exchanger, Int. J. Heat Mass Transf. 44 (2001) 3721–3730. [28] C. Dong, Y.-P. Chen, J.-F. Wu, Flow and heat transfer performances of helical baffle heat exchangers with different baffle configurations, Appl. Therm. Eng. 80 (2015) 328–338. [29] M. Fratczak, J. Czeczot, P. Nowak, Simplified modeling of plate heat exchangers, 2014. 19th International Conference on Methods and Models in Automation and Robotics (MMAR) 578-583, Mie˛dzyzdroje, Poland. [30] J.V. Villadsen, M.L. Michelsen, Solution of Differential Equation Models by Polynomial Approximation, Prentice Hall, Englewood Cliffs, NJ, 1978. [31] J.A. Nelder, R. Mead, A simplex method for function minimization, Comput. J. 7 (1965) 308–313. [32] J.C. Lagarias, J.A. Reeds, M.H. Wright, P.E. Wright, Convergence properties of the Nelder-Mead simplex method in low dimensions, SIAM J. Optim. 9 (1) (1998) 112–147. [33] W. Waldraff, D. Dochain, S. Bourrel, A. Magnus, On the use of the observability measures for sensor location in tubular reactor, J. Process Control 8 (5–6) (1998) 497–505.

893

[34] D. Dochain, State observers for tubular reactors with unknown kinetics, J. Process Control 10 (2000) 259–268. [35] M.B. Carver, H.W. Hinds, The method of lines and advective equation, Simulation 31 (1978) 59–69. [36] C.G. Lee, S.C. Park, Survey on the virtual commissioning of manufacturing systems, J. Comput. Des. Eng. 1 (3) (2014) 213–222. [37] P. Laszczyk, J. Czeczot, R. Czubasiewicz, K. Stebel, Practical Verification of the Control Strategies for the Counter-Current Heat Exchanger, 2012. 17th International Conference on Methods and Models in Automation and Robotics (MMAR) 96–101, Mie˛dzyzdroje, Poland. [38] M. Fratczak, P. Nowak, T. Klopot, J. Czeczot, S. Bysko, B. Opilski, Virtual commissioning for the control of the continuous industrial processes – case study, 2015. 20th International Conference on Methods and Models in Automation and Robotics (MMAR), Mie˛dzyzdroje, Poland. [39] J.P. Babary, S. Julien, M.T. Nihtila, J. Czeczot, M. Metzger, New boundary conditions and adaptive control of fixed-bed bioreactors, Chem. Eng. Process 38 (1999) 35–44. [40] A. Maidi, M. Diaf, J.P. Corriou, Boundary geometric control of a counter-current heat exchanger, J. Process Control 19 (2009) 297–313. [41] A. Maidi, M. Diaf, J.P. Corriou, Boundary control of a parallel-flow heat exchanger by input-output linearization, J. Process Control 20 (2010) 1161–1174. [42] M. Fratczak, J. Czeczot, P. Nowak, Simulation validation of the model-based control of the plate heat exchanger with on-line compensation for modelling inaccuracies, in: Proc. of the 4th International Conference on Simulation and Modeling Methodologies, Technologies and Applications (SIMULTECH-2014), SCITEPRESS (Science and Technology Publications, Lda.), 2014, pp. 657–665, doi:10.5220/0005015006570665. ISBN: 978-989-758-038-3. [43] P. Kadlec, B. Gabrys, S. Strandt, Data-driven soft sensors in the process industry, Comput. Chem. Eng. 33 (4) (2009) 795–814. [44] C. Kravaris, J. Hahn, Y. Chu, Advances and recent developments in state and parameter estimation, Comput. Chem. Eng. 51 (2013) 111–123. [45] P. Nowak, J. Czeczot, Observer-Based Cascade Control of the Heat Distribution System, 2013. 18th International Conference on Methods and Models in Automation and Robotics (MMAR) 633-638, Mie˛ dzyzdroje, Poland.