Abstract

Heat transfer in counterflow heat exchangers is modeled by using transport and balance equations with the temperatures of cold fluid, hot fluid, and metal pipe as state variables distributed along the entire pipe length. Using such models, boundary value problems can be solved to estimate the temperatures over all the length by means of measurements taken only at the boundaries. Conditions for the stability of the estimation error given by the difference between the temperatures and their estimates are established by using a Lyapunov approach. Toward this end, a method to construct nonlinear Lyapunov functionals is addressed by relying on a polynomial diagonal structure. This stability analysis is extended in case of the presence of bounded modeling uncertainty. The theoretical findings are illustrated with numerical results, which show the effectiveness of the proposed approach.

1. Introduction

The mathematical modeling of heat exchangers is really important since these devices are employed in chemical processes, petroleum refineries, power systems, and heating/refrigeration of air-conditioning plants [1, 2]. In this paper, we address the modeling of a counterflow heat exchanger based on transport and balance PDEs and use such models for the purpose of monitoring the thermodynamic process by solving boundary value problems that exploit only few temperature measurements at the boundaries. A rigorous stability analysis is addressed to prove the effectiveness of the resulting temperature estimates from the theoretical and numerical points of view in line with the widespread methods based on the Lyapunov approach in a number of different applications [36].

Heat exchangers for industrial processes allow recovering lost energy from hot fluid streams or to heat a cold fluid for the purpose of air conditioning. Finite-dimensional models of such processes are used for fouling detection by using Kalman filtering methods [1, 2]. As compared with these models, the proposed modeling framework is in an infinite dimension since it relies on hyperbolic PDEs. Based on such equations, first we will focus on a model with the dynamics of the temperatures of cold fluid, hot fluid, and metal pipe. In addition, a reduced model is considered with only the temperatures of cold and hot fluid as state variables in line with [7]. Since such temperatures cannot be determined by direct, distributed measurements, it is fundamental to estimate them by using a few numbers of measures. Toward this end, using such models, we address the solution of boundary value problems to estimate the distributed state, which provides such estimate by having at our disposal measurements of temperature only at the boundaries. We will refer to them also as boundary-output Luenberger observers or simply state observers (see [8] and the references therein). Stability conditions on the estimation error (i.e., the difference between the temperatures and their estimates) will be presented to construct such estimators in case of perfect modeling and assuming bounded uncertainty. In such a case, we will rely on the notion of quadratic boundedness (QB, for short) [9]. Finally, numerical results will be given to showcase the theoretical findings.

The stability of the estimation error is analyzed in two different settings. First, we will consider perfect modeling and hence asymptotic stability results will be provided. Next, we will assume the presence of additive noises, which account for bounded modeling uncertainties and stability is studied by using QB. Generally speaking, QB allows dealing with positively invariant sets and derive upper bounds on the trajectories of the state of a dynamic system subject to bounded disturbances [9]. Thus, QB has been successfully used for both output feedback control [10] and state estimation [11]. Here, we will extend the use of QB for the purpose of estimation in the considered infinite dimensional context. It is worth noting that, in lieu of searching for exact solutions of the flow equations (see, e.g., [12, 13]), here we deal with the construction of boundary-output observers with guaranteed stability properties on the estimation error.

All the stability conditions that will be presented demand the selection of nonlinear Lyapunov functionals. In this paper, we will address this task by using the SOS (sum-of-squares) approach [14, 15], which was originally developed for the analysis and control of polynomial systems described by ODEs (see [16] and the references therein) and recently applied also to systems described by PDEs [1719]. SOS methods are quite computationally efficient since they are based on semidefinite programming (SDP) [20, 21]. In this respect, the SOS approach allows finding Lyapunov functionals just like the methods based on LMIs (linear matrix inequalities) enable computing Lyapunov functions for linear systems based on ODE models [22]. As compared with the search of classical Lyapunov functions to investigate local stability [23], we do not rely on linearized models but provide conditions of global stability.

The paper is organized as follows: Section 2 reports the basic definitions that will be used in the following discussions. Modeling of thermal flows in one-dimensional heat exchangers is given in Section 3. The proposed estimation approach is presented in Section 4, where stability is analyzed in a noise-free modeling framework. This analysis is extended to models affected by bounded disturbances in Section 5. A brief account on the use of the SOS method to find the Lyapunov functionals and thus ensure stability is given in Section 6. Section 7 illustrates the numerical results, while conclusions are drawn in Section 8.

2. Notations and Definitions

The set of real numbers is denoted by and hence will be used for the set of strictly positive real numbers. The set of the integer numbers equal to or greater than zero is denoted by .

The minimum and maximum eigenvalues of a symmetric matrix are denoted by and , respectively. Moreover, () means that it is also positive (negative) definite. Given a generic matrix , , and hence, for a vector , is its Euclidean norm. denotes the ring of real polynomials in ; is the set of real polynomials of a degree equal or less than .

Given the compact set , denotes the Hilbert space of square integrable functions with the norm for all . The solution of a PDE with the initial condition is said to be (i) is stable (to zero) if, for all , there exists such that for all (ii) is asymptotically stable (to zero) if it is stable and (iii) is exponentially stable (to zero) if there exists such that for some and all

3. Modeling Heat Transfer in Heat Exchangers

Consider the one-dimensional model of counterflow heat exchangers with three distributed state variables, i.e., , , and for and , as depicted in Figure 1; such variables represent the temperatures of cold fluid, hot fluid, and pipe, respectively (the subscript indicating that it is made of some metal in such a way to maximize the heat transfer). Based on the balance of enthalpy for the hot fluid, we get where is the density (in ), is the specific heat (in J/(kg K)), is the mass flow rate (in kg/s), is the transfer coefficient (in W/(m2K)), (in m2) is the internal surface of the pipe, and is its perimeter (in m). Concerning the cold fluid, the balance of enthalpy is given by where is density (in kg/m3), is the specific heat (in J/(kg K)), and is the transfer coefficient (in W/(m2)). Moreover, (in m2) is the external surface of the pipe and is its perimeter (in m). Finally, for the metal pipe it follows that where is the density (in kg/m3) of the metal, is the specific heat (in J/(kg K)) of the metal, and (in m2) is the surface of the pipe section.

After dividing each of the three previous equations for the first coefficient in each respective l.h.s., we get where or, more concisely, as follows: where and

The system (11) is solved with initial conditions and Dirichlet boundary conditions for all . The boundary conditions on are not necessary since equation (9) does not contain spatial derivatives of .

Likewise, in [7], a reduced model can be easily derived from (11) by neglecting the dynamics of the temperature in the pipe, i.e., assuming that such a temperature is fixed by the values of the temperatures of cold and hot fluids, i.e., imposing that the r.h.s. of (6) is null. Therefore, we obtain or, more compactly, where and with initial conditions and Dirichlet boundary conditions for all .

Model (11) is of a hyperbolic type but not strictly hyperbolic since matrix has a zero eigenvalue. By contrast, (16) is strictly hyperbolic since the eigenvalues of have a non-null-real part. In Section 4, we will consider the problem to reconstruct all the temperatures by using a generic hyperbolic model.

4. Estimation for Heat Exchangers with Stable Estimation Error

Consider the hyperbolic equation where , diagonal, and . Consider also the observer for (20) given by where is the estimate of . Such PDEs need some suitable boundary conditions, which will be given later. The stability of the estimation error can be guaranteed under suitable conditions as follows.

Theorem 1. The estimation error is asymptotically stable if there exists a diagonal matrix with for all such that and for all .

Proof. Consider the Lyapunov functional and note that the asymptotic stability of the estimation error is ensured if for all since for all . Since does not depend on time, we get and, owing to the diagonal structure of and , Thus, it follows and hence (25) holds owing to (22) and (23).

Note that if, instead of (23), we consider with for all , the estimation error is exponentially stable with a rate of decrease equal to . Thus, a design goal may consist in maximizing , which can be regarded as a sort of generalized eigenvalue problem for Lyapunov functionals instead of Lyapunov functions [22].

In case of the reduced model (16) (i.e., with and ), the stability conditions of Theorem 1 can be greatly simplified according to [7] by using and boundary conditions for all with to be suitably chosen. More specifically, in [7] it is proposed to satisfy (23) by using an LMI-based discretized approach to get and then select such that by choosing, for example,

The construction of Lyapunov functionals is not an easy task to solve in general and will be addressed later in Section 6.

5. Estimation under Modeling Uncertainties

Let us focus on estimation in the presence of bounded disturbances. More specifically, instead of (20) we consider where and such that , for all and without loss of generality (different upper bounds can be taken into account by suitably scaling the coefficients of ), where is an integer number no larger than . Estimation will be performed by using the same observer adopted in the noise-free case, i.e., (21).

The presence of the system noises prevents from ensuring asymptotic stability, but since such disturbances are bounded, an invariant set for the estimation error exists and can be studied by using QB [9].

Toward this end, first of all we need to define QB in the sense. More specifically, the estimation error is said to be quadratically bounded with the Lyapunov functional as defined in (24) if

As a matter of fact, any positive constant can be chosen instead of , but with such a choice we reduce the notational burden. This does not entail loss of generality since the conditions ensuring QB are homogeneous in the design parameters, which scales with and thus allows for this simplification.

Owing to (36), the set turns out to be positively invariant and it is attractive (i.e., if the error is out of , it approaches asymptotically).

Theorem 2. The estimation error is quadratically bounded if there exists a diagonal matrix with for all and , such that (22) and hold for all .

Proof. Likewise in the proof of Theorem 1, we easily compute the time derivative of the Lyapunov functional (24) and get that subject to (22) is equivalent to Since implies and holds, using the well-known S-procedure (see [22], p. 23) we obtain that QB holds if for some and or equivalently if (38) and (39) hold.

Clearly, (38) implies for all , which ensures that the estimation error in the absence of noises is exponentially stable with a rate of decrease equal to . Thus, a convenient design goal may consist in maximizing , as pointed out in [7]. Another, quite popular objective of design is related to the gain between disturbance and estimation error.

Theorem 2 provides conditions ensuring that the set is attractive. Thus, one can design the estimator in such a way to keep as small as possible. Toward this end, note that for all and hence, as a design objective, we may maximize the minimum eigenvalue of since for all , . The inequality above allows proving the boundedness in the sense as for all .

To compare the effectiveness of a given observer in terms of rejection of the noise, it is quite popular to rely on the notion of the gain by assuming finite-energy noises, i.e., . More specifically, we say that the proposed observer admits and gain between disturbance and estimation error is equal to if for . It is straightforward to show that the bound of the gain holds if

Note also that this condition is quite similar to (38). Moreover, minimizing is equivalent to maximizing in accordance with the goal of maximizing .

6. Search of Lyapunov Functionals Using the SOS Approach

Stability conditions such as (23), (29), (38), or (47) (each together with ) demand the satisfaction of LMIs in for all , i.e., the solution of semi-infinite programming problems. Such problems arise in a variety of applications and are usually approximately solved in discretized form on a sufficiently fine mesh of points (see, e.g., [24]). The main difficulty to address semi-infinite programming problems concerns both the choice of a sufficiently large number of points and especially the local minima trapping, by which nonlinear programming solvers may be affected in trying to find the solution. A more appealing way to solve such problems consists in resorting to the SOS paradigm, which enables turning the construction of a Lyapunov functional into a convex problem that can be efficiently solved by using SDP without encountering such issues due to local minima.

The idea behind such an approach is the SOS decomposition of a candidate Lyapunov functional as well as of the opposite of its time derivative by using a positivity certification, which does not depend on the characteristics of the chosen polynomial [21, 25]. More specifically, the following result holds.

Theorem 3. A polynomial in has sum-of-squares decomposition (or is said to be SOS) if and only if there exists a real symmetric and positive semidefinite matrix such that , where is the vector of all the monomials in the components of of degree equal to or less than , i.e., of dimension

Proof. See Proposition 2.1 in [26] (p. 17).

Consider, for example, the estimation of the state variables of (7), (8), and (9) by using the results of Theorem 1. Based on Theorem 3, a procedure can be adopted to find suitable SOS polynomials as diagonal elements of , all to be taken in with for increasing values of . Toward this end, we say that is -SOS polynomial if with , is SOS for some “small” tolerance [27]. In addition, a polynomial, square matrix with for is said to be an -SOS matrix if is -SOS in with and .

Given with and , the problem to solve when dealing, for example, with the stability conditions of Theorem 1 is the following:

Problem 4. Find such that for all .

If Problem 4 admits a solution, we can follow the same reasoning that leads to (34) to ensure the satisfaction of (22) by choosing

Similar problems may be formulated by replacing the SOS description of (52) with those corresponding to (29), (38), and (47).

7. Numerical Results

For the sake of brevity, in the following discussion, we will describe only the numerical solution of the boundary value problem that results from the discretization of the model (11) since this includes that of the reduced model (16) as a special case. We have discretized the domain with equally spaced points for , and both (16) and (21) are treated by using a finite difference method. Equations (7) and (8) composed of an advection term ( and ) and a reaction one (the right side of the equalities) are semidiscretized in space by an upwind finite difference scheme [2830]. The choice of an upwind scheme is more appropriate for the convective term in which there is a propagation direction. The boundary conditions (19) for (16) (and (31) and (32) for (21)) are given for and only in the left () and right () parts of the domain, respectively. The term in (7) is therefore discretized by using backward finite difference since is positive, whereas the term in (8) is discretized by a forward one since is negative.

Setting by the solution at the discretization point at time (and analogously for and ), we obtain the following system of ODEs: with initial conditions

The boundary conditions in (19) are translated in the system by setting

The temperatures and at the initial time are chosen constant and equal to the values at the boundary points where their measures are available, i.e., equals to in and in such as with a constant value belonging to the interval taken as initial conditions for , i.e.,

The scheme for the observer (21) is the same by replacing with , with and so on, except for the boundary conditions, which are the following:

For the observer, the initial conditions in the simulation are chosen as:

The system of ODEs (52) is solved by using a fourth-order Runge-Kutta method (implemented in the ode45 function of MATLAB) with variable time step. In the first instants of the simulation, the fast dynamics of the convective part prevails, while later the reaction term predominates. To accurately simulate the model, it is therefore necessary to select a very small time step in the first time instants, while it may be chosen larger afterwards, but of course, always respecting the stability condition of the numerical scheme. Implicit methods (e.g., the scheme implemented in the routine ode15s of MATLAB) are tested, since they generally allow using a larger time step. This resulted in a longer running time, due to the higher complexity to complete each iteration. Then, explicit methods seem to be the most appropriate choice for this problem.

In the numerical case study, we have fixed  m,  m/s,  m/s, 1/s, 1/s, 1/s, and 1/s. Concerning the proposed estimation approach with (11) and (16), we have solved Problem 4 by using the SOS toolbox [27] and the numerical solver Yalmip [31] with . In the first case (i.e., and ), we have got a solution to such a problem with given by and select according to (52). The same has been done for the reduced model with and by obtaining and taking again according to (52). We have chosen , , , and .

Figure 2 shows the behavior of the temperatures , , and and their estimates given by , , and at different time instants, all based on (11). Figure 3 shows the behavior of the temperatures , and the corresponding estimates , , based on the reduced model (16).

8. Conclusions

Two models have been presented to account for the temperature dynamics in heat exchangers. Such models are based on hyperbolic PDEs with the most complete having the temperatures of cold fluid, hot fluid, and pipe as state variables, whereas the reduced one relies only on the temperatures of cold and hot fluids. Estimators resulting from the solution of boundary value problems based on such models have been proposed and provided with a rigorous stability analysis with the Lyapunov approach. To this end, the existence of nonlinear polynomial Lyapunov functionals has been addressed by using the SOS approach.

Future work will concern the investigation on the stability of systems described by nonlinear hyperbolic equations (see, e.g., [32]). Another direction of research is the study of stability of cascaded hyperbolic and other nonlinear PDEs such as the normal flow equation to deal with multiphase problems, which turn out to be increasingly difficult but with a wide range of potential applications [33].

Data Availability

The simulation code is available upon request to Angelo Alessandri, email: [email protected].

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This work was supported by Intermarine SpA under contract no. 2967/2017.