# Optimal Design of a Solar Desalination Unit with Heliostats

Authors: M. Abidi; S. Ben Jabrallah; J. P. Corriou
Journal: International Journal of Chemical Engineering
 (2017)

---

## Abstract

The objective is to improve the yield of a solar desalination cell using concentration of solar rays by means of automatically controlled heliostats. The vertical cell is orientated towards the north. It is mainly composed of two plates; the one being heated by the solar rays reflected by the mirrors is used for evaporation of a falling water film; the other one is used for water vapor condensation. Each heliostat consists of an altitude-azimuth mount having two degrees of freedom and supporting a plane mirror. The heliostat permanently follows the sun trajectory and reflects the solar rays on the cell by means of automatic control implemented in a control card based on a microcontroller. Model predictive control allows us to maximize the distilled water production.


---

## Body

## 1. Introduction

In many countries, the shortage of potable water becomes a huge problem for all individuals. On another side, given the society development, water consumption increases tremendously. Water supply in arid and semiarid regions becomes such an important issue that it motivates researchers to design desalination systems that exploit renewable energy sources. Most autonomous desalination systems use solar radiation. The studies by [1, 2] show that the production flow rate of the cell is significantly improved with the increase of the heat flux density q˙f. The solution requires the increase of the number of solar sources by means of a simple multiplication of the number of reflecting mirrors [3, 4] which are used for various purposes such as energy production [5], material study at high temperatures, hydrogen production [6], and heating [7]. Many types of heliostats exist [8]. Their control is the subject of much industrial activity, with both open loop or closed loop choices [6, 9–12]. Patents are often taken [13]. In some cases [12], four degrees of freedom, azimuth, elevation, rotation, and pitch, are even used. In some applications, heliostats with tower plants are even used [14–16].The objective of the present paper is to improve the production of a solar distillation cell that is strongly dependent on the solar flux. The use of several reflecting mirrors and the control of their individual position to perform an adequate tracking of the sun will allow us to obtain a maximum reflection of the solar rays on the cell and thus an increase of the incident heat flux density to the cell. Thus, the production of distilled water will be maximized. In Section1, the objectives were presented and the literature reviewed. Section 2 deals with the description of the complex model of the cell based on partial differential equations describing mass and heat balances are described. In Section 3, a reduced dynamic model is built in order to catch the main dynamic characteristics of the cell for future use in Model Predictive Control (MPC). In Section 4, the experimental setup is described briefly. In Section 5, MPC is briefly described in general and under the form which was used for this work. In Section 6, it is explained how the production of distilled water can be maximized by MPC and experimental results are given and commented on. Finally, a conclusion ends the paper. The positioning of the heliostats is described with more details in Appendix A.


## 2. First-Principles Dynamic Model of the Distillation Cell

The distillation cell, later denoted as the plant, illustrated in Figure1, is a parallelepiped cavity, of form factor height/width equal to 10, mainly consisting of two vertical and parallel metal plates (of dimensions 0.5 m × 0.5 m). One plate uniformly heated by solar energy is used for evaporating the falling film; the other one is at a lower temperature for condensing the water vapor. The steady state process was modelled by Ben Jabrallah et al. [2] using two-dimensional mass, energy, and momentum balances. The influence of main operating parameters could be partly explained by the fine description of the fluid circulation in the cell. The steady state model was later experimentally validated by the same authors [17] and the influence of the different operating parameters on the yield was compared to theoretical predictions with a good agreement. However, in order to perform the present study along days, the previous model has been extended as a dynamic model by Abidi et al. [1] composed of partial differential equations to represent the transient behaviour. In the same way as the steady state model, the model describing the dynamic behaviour of the distillation cell is divided into three parts corresponding, respectively, to the heated wall, the liquid phase, and the gas phase (Figure 1). The conservation equations are completed by initial and boundary conditions. The respective domains of variation of the coordinates x, spatial coordinate in the horizontal direction, and y, spatial coordinate in the vertical direction, are(1)0≤y≤h,0=bottom,h=top,0≤xs≤Lp,from left to right in the solid domain, heated plate,0≤xw≤δ,from left to right in the liquid film domain,δ≤xg≤b,from left to right in the gas domain.Figure 1
Scheme of the distillation cell for modelling.


### 2.1. Heated Metal Plate

The conductive heat transfer through the metal plate (subscript “s”) heated at the outside by solar radiation follows Fourier law as(2)∂Ts∂t=λsρsCps∂2Ts∂xs2+∂2Ts∂y2,where xs is the spatial coordinate in the transverse direction of the solid plate and y the vertical direction.


### 2.2. Liquid Water Film

The liquid water film (subscript “w”) falls on the internal surface of the heated plate. Its behaviour is described by momentum and energy balances as(3)∂∂xwμw∂vw∂xw+ρwg=ρw∂vw∂t,∂Tw∂t+vw∂Tw∂y=λwρCpw∂2Tw∂xw2,where xw is the spatial coordinate in the transverse direction of the liquid water film.


### 2.3. Gas Phase

Some simplifying hypotheses are adopted for the gas phase, that is, the flow is assumed to be only two-dimensional; the densities of the water vapor and air gas mixture are constant except in the buoyancy force term (Boussinesq approximation). In the following equations, the subscript “g” of gas is omitted for convenience in all variables x, u, v, T, and C. The movement of the gas mixture in the cavity is described by the following conservation equations:Continuity equation(4)∂u∂x+∂v∂y=0.Momentum equations(5)ρ∂u∂t+u∂u∂x+v∂u∂y=-∂p∂x+μ∂2u∂x2+∂2u∂y2,ρ∂v∂t+u∂v∂x+v∂v∂y=-∂p∂y+μ∂2v∂x2+∂2v∂y2-ρg.Energy equation(6)∂T∂t+u∂T∂x+v∂T∂y=λρCp∂2T∂x2+∂2T∂y2.Transport of the diffusing species(7)∂C∂t+u∂C∂x+v∂C∂y=D∂2C∂x2+∂2C∂y2.


### 2.4. Boundary Conditions

#### 2.4.1. Heated Metal Plate

The heat flux densityq˙f corresponding to solar radiation is imposed on the external face of the heated metal plate of thickness Lp.(8)-λs∂Ts∂xsxs=0=q˙f,∀y,∀t.It is assumed that the two horizontal faces of the cell are adiabatic.(9)∂Ts∂y=0,aty=0,y=h,0≤xs≤Lp,∀t.


#### 2.4.2. Heated Plate-Water Film Interface

At the interface between the solid plate and the water film, the temperature should verify the continuity equation(10)Tsxs=Lp=Twxw=0,∀y,∀t.Similarly, the heat flux must verify the continuity equation(11)λs∂Ts∂xsxs=Lp=λw∂Tw∂xwxw=0,∀y,∀t. The adherence condition of the water film is(12)vwxw=0=0,∀y,∀t.


#### 2.4.3. In the Water Film

The temperature of the water film at the inlet of the cell is given by(13)Twy=h=Tf,0≤xw≤δ,∀t.


#### 2.4.4. Liquid-Gas Interface

At the liquid film-gas phase interface, the continuity condition for temperature and velocity and equilibrium condition for concentration are, respectively,(14)Twxw=δ=Tgxg=δ,∀y,∀t,vwxw=δ=vgxg=δ,∀y,∀t,Csatxw=δ=Cxg=δ,∀y,∀t.The continuity equation for the shear stresses can be written as(15)μw∂vw∂xxw=δ=μg∂vg∂xxg=δ,∀y,∀t.The heat flux also verifies the continuity equation(16)-λg∂Tg∂xgxg=δ=-λw∂Tw∂xwxw=δ-ρgugxg=δLv,∀y,∀t,where the second term in the right-hand stands for the liquid-gas change of state.The transverse gas velocity at the interface results from Fick diffusion law and it is assumed that the air-water interface is impermeable to dry air [18].(17)ugxg=δ=-D1-wxg=δ∂w∂xxg=δ,where wδ is the mass fraction of water vapor at the interface and D is the diffusion coefficient of water in air. The mass flux density m˙v exchanged at the liquid-gas interface results as [17](18)m˙v=ρgugxg=δ.


#### 2.4.5. In the Gas Phase

It is assumed that the horizontal walls are perfectly adiabatic.(19)∂Tg∂y=∂Cg∂y=0,aty=0,y=h,∀δ≤xg≤b,∀t.On the horizontal walls, the adherence condition is verified as follows:(20)ug=vg=0,aty=0,y=h,∀δ≤xg≤b,∀t.On the condensation plate, the adherence condition is verified as follows:(21)vg=0,atxg=b,∀y,∀t.On the condensation plate, the transverse velocity is deduced from Fick law similarly to (17) as follows:(22)ug=-D1-wgxg=b∂wg∂xxg=b.On the condensation plate, a uniform temperature is imposed and consequently the corresponding equilibrium concentration results.(23)Tg=Tc,Cg=Csat,atxg=b,∀y,∀t.


#### 2.4.6. Relation with Distillation Flow Rate

The flow of water vapor inside the cavity results from the thermal and mass gradients. The water vapor in the cavity condenses when it comes into contact with the condensation plate, maintained at a relatively low temperature. The exchanged mass flux densitym˙Φ is obtained by application of Fick law at the contact surface between the vapor and the condensation plate as follows:(24)m˙Φ=ρgugx=b=-ρgD1-wgx=b∂wg∂xx=b,where wg is the mass fraction of water vapor at the condensation plate, D the diffusion coefficient of water vapor in air, and C the concentration of water vapor.The distillation flow rateFd of the cell proportional to the concentration gradient of the diffusing species at the interface between the gas phase and the condensation plate is(25)Fd=∫0hm˙ϕldy,where l is the length of the film.


### 2.5. Numerical Solution and Validation

The desalination cell is an experimental prototype of size 0.5 m × 0.5 m that was described in detail in [1, 2].The steady state model of the plant, consisting of partial differential equations (PDEs) in two dimensions and boundary and initial conditions, was solved by means of the method of finite volumes [2, 19]. The latter method respects the final conservation of the balance equations opposite to the simpler finite differences method which does not guarantee it. The steady state model was validated experimentally by Ben Jabrallah et al. [17]. The dynamic model adds the time dimension to the steady state model and is more difficult to solve because the convergence should be verified at each time step. Thus, as the accuracy of the open loop results is not of prime importance and some limited error is acceptable, it was chosen to perform a spatial discretization according to the method of lines so that the partial differential equations with respect to time and space result in a large system of first-order ordinary differential equations (ODEs). This system of ODEs was written as a tridiagonal system according to the horizontal and vertical directions. To solve it, a double sweeping method was adopted which presents the advantage to rapidly propagate the influence of the boundary conditions towards the inside of the cell and to accelerate the convergence. With the solution being searched at each time step, the boundary conditions were updated at each iteration [1]. As the dynamic model is deduced from the steady state model by use of first principles, it is also assumed that the steady state validation is sufficient for its dynamic purpose. A Fortran program was thus developed which simulates the dynamic behaviour of the cell with respect to the variation of its different manipulated inputs. This model of the plant will serve as the reference along the study.


## 2.1. Heated Metal Plate

The conductive heat transfer through the metal plate (subscript “s”) heated at the outside by solar radiation follows Fourier law as(2)∂Ts∂t=λsρsCps∂2Ts∂xs2+∂2Ts∂y2,where xs is the spatial coordinate in the transverse direction of the solid plate and y the vertical direction.


## 2.2. Liquid Water Film

The liquid water film (subscript “w”) falls on the internal surface of the heated plate. Its behaviour is described by momentum and energy balances as(3)∂∂xwμw∂vw∂xw+ρwg=ρw∂vw∂t,∂Tw∂t+vw∂Tw∂y=λwρCpw∂2Tw∂xw2,where xw is the spatial coordinate in the transverse direction of the liquid water film.


## 2.3. Gas Phase

Some simplifying hypotheses are adopted for the gas phase, that is, the flow is assumed to be only two-dimensional; the densities of the water vapor and air gas mixture are constant except in the buoyancy force term (Boussinesq approximation). In the following equations, the subscript “g” of gas is omitted for convenience in all variables x, u, v, T, and C. The movement of the gas mixture in the cavity is described by the following conservation equations:Continuity equation(4)∂u∂x+∂v∂y=0.Momentum equations(5)ρ∂u∂t+u∂u∂x+v∂u∂y=-∂p∂x+μ∂2u∂x2+∂2u∂y2,ρ∂v∂t+u∂v∂x+v∂v∂y=-∂p∂y+μ∂2v∂x2+∂2v∂y2-ρg.Energy equation(6)∂T∂t+u∂T∂x+v∂T∂y=λρCp∂2T∂x2+∂2T∂y2.Transport of the diffusing species(7)∂C∂t+u∂C∂x+v∂C∂y=D∂2C∂x2+∂2C∂y2.


## 2.4. Boundary Conditions

### 2.4.1. Heated Metal Plate

The heat flux densityq˙f corresponding to solar radiation is imposed on the external face of the heated metal plate of thickness Lp.(8)-λs∂Ts∂xsxs=0=q˙f,∀y,∀t.It is assumed that the two horizontal faces of the cell are adiabatic.(9)∂Ts∂y=0,aty=0,y=h,0≤xs≤Lp,∀t.


### 2.4.2. Heated Plate-Water Film Interface

At the interface between the solid plate and the water film, the temperature should verify the continuity equation(10)Tsxs=Lp=Twxw=0,∀y,∀t.Similarly, the heat flux must verify the continuity equation(11)λs∂Ts∂xsxs=Lp=λw∂Tw∂xwxw=0,∀y,∀t. The adherence condition of the water film is(12)vwxw=0=0,∀y,∀t.


### 2.4.3. In the Water Film

The temperature of the water film at the inlet of the cell is given by(13)Twy=h=Tf,0≤xw≤δ,∀t.


### 2.4.4. Liquid-Gas Interface

At the liquid film-gas phase interface, the continuity condition for temperature and velocity and equilibrium condition for concentration are, respectively,(14)Twxw=δ=Tgxg=δ,∀y,∀t,vwxw=δ=vgxg=δ,∀y,∀t,Csatxw=δ=Cxg=δ,∀y,∀t.The continuity equation for the shear stresses can be written as(15)μw∂vw∂xxw=δ=μg∂vg∂xxg=δ,∀y,∀t.The heat flux also verifies the continuity equation(16)-λg∂Tg∂xgxg=δ=-λw∂Tw∂xwxw=δ-ρgugxg=δLv,∀y,∀t,where the second term in the right-hand stands for the liquid-gas change of state.The transverse gas velocity at the interface results from Fick diffusion law and it is assumed that the air-water interface is impermeable to dry air [18].(17)ugxg=δ=-D1-wxg=δ∂w∂xxg=δ,where wδ is the mass fraction of water vapor at the interface and D is the diffusion coefficient of water in air. The mass flux density m˙v exchanged at the liquid-gas interface results as [17](18)m˙v=ρgugxg=δ.


### 2.4.5. In the Gas Phase

It is assumed that the horizontal walls are perfectly adiabatic.(19)∂Tg∂y=∂Cg∂y=0,aty=0,y=h,∀δ≤xg≤b,∀t.On the horizontal walls, the adherence condition is verified as follows:(20)ug=vg=0,aty=0,y=h,∀δ≤xg≤b,∀t.On the condensation plate, the adherence condition is verified as follows:(21)vg=0,atxg=b,∀y,∀t.On the condensation plate, the transverse velocity is deduced from Fick law similarly to (17) as follows:(22)ug=-D1-wgxg=b∂wg∂xxg=b.On the condensation plate, a uniform temperature is imposed and consequently the corresponding equilibrium concentration results.(23)Tg=Tc,Cg=Csat,atxg=b,∀y,∀t.


### 2.4.6. Relation with Distillation Flow Rate

The flow of water vapor inside the cavity results from the thermal and mass gradients. The water vapor in the cavity condenses when it comes into contact with the condensation plate, maintained at a relatively low temperature. The exchanged mass flux densitym˙Φ is obtained by application of Fick law at the contact surface between the vapor and the condensation plate as follows:(24)m˙Φ=ρgugx=b=-ρgD1-wgx=b∂wg∂xx=b,where wg is the mass fraction of water vapor at the condensation plate, D the diffusion coefficient of water vapor in air, and C the concentration of water vapor.The distillation flow rateFd of the cell proportional to the concentration gradient of the diffusing species at the interface between the gas phase and the condensation plate is(25)Fd=∫0hm˙ϕldy,where l is the length of the film.


## 2.4.1. Heated Metal Plate

The heat flux densityq˙f corresponding to solar radiation is imposed on the external face of the heated metal plate of thickness Lp.(8)-λs∂Ts∂xsxs=0=q˙f,∀y,∀t.It is assumed that the two horizontal faces of the cell are adiabatic.(9)∂Ts∂y=0,aty=0,y=h,0≤xs≤Lp,∀t.


## 2.4.2. Heated Plate-Water Film Interface

At the interface between the solid plate and the water film, the temperature should verify the continuity equation(10)Tsxs=Lp=Twxw=0,∀y,∀t.Similarly, the heat flux must verify the continuity equation(11)λs∂Ts∂xsxs=Lp=λw∂Tw∂xwxw=0,∀y,∀t. The adherence condition of the water film is(12)vwxw=0=0,∀y,∀t.


## 2.4.3. In the Water Film

The temperature of the water film at the inlet of the cell is given by(13)Twy=h=Tf,0≤xw≤δ,∀t.


## 2.4.4. Liquid-Gas Interface

At the liquid film-gas phase interface, the continuity condition for temperature and velocity and equilibrium condition for concentration are, respectively,(14)Twxw=δ=Tgxg=δ,∀y,∀t,vwxw=δ=vgxg=δ,∀y,∀t,Csatxw=δ=Cxg=δ,∀y,∀t.The continuity equation for the shear stresses can be written as(15)μw∂vw∂xxw=δ=μg∂vg∂xxg=δ,∀y,∀t.The heat flux also verifies the continuity equation(16)-λg∂Tg∂xgxg=δ=-λw∂Tw∂xwxw=δ-ρgugxg=δLv,∀y,∀t,where the second term in the right-hand stands for the liquid-gas change of state.The transverse gas velocity at the interface results from Fick diffusion law and it is assumed that the air-water interface is impermeable to dry air [18].(17)ugxg=δ=-D1-wxg=δ∂w∂xxg=δ,where wδ is the mass fraction of water vapor at the interface and D is the diffusion coefficient of water in air. The mass flux density m˙v exchanged at the liquid-gas interface results as [17](18)m˙v=ρgugxg=δ.


## 2.4.5. In the Gas Phase

It is assumed that the horizontal walls are perfectly adiabatic.(19)∂Tg∂y=∂Cg∂y=0,aty=0,y=h,∀δ≤xg≤b,∀t.On the horizontal walls, the adherence condition is verified as follows:(20)ug=vg=0,aty=0,y=h,∀δ≤xg≤b,∀t.On the condensation plate, the adherence condition is verified as follows:(21)vg=0,atxg=b,∀y,∀t.On the condensation plate, the transverse velocity is deduced from Fick law similarly to (17) as follows:(22)ug=-D1-wgxg=b∂wg∂xxg=b.On the condensation plate, a uniform temperature is imposed and consequently the corresponding equilibrium concentration results.(23)Tg=Tc,Cg=Csat,atxg=b,∀y,∀t.


## 2.4.6. Relation with Distillation Flow Rate

The flow of water vapor inside the cavity results from the thermal and mass gradients. The water vapor in the cavity condenses when it comes into contact with the condensation plate, maintained at a relatively low temperature. The exchanged mass flux densitym˙Φ is obtained by application of Fick law at the contact surface between the vapor and the condensation plate as follows:(24)m˙Φ=ρgugx=b=-ρgD1-wgx=b∂wg∂xx=b,where wg is the mass fraction of water vapor at the condensation plate, D the diffusion coefficient of water vapor in air, and C the concentration of water vapor.The distillation flow rateFd of the cell proportional to the concentration gradient of the diffusing species at the interface between the gas phase and the condensation plate is(25)Fd=∫0hm˙ϕldy,where l is the length of the film.


## 2.5. Numerical Solution and Validation

The desalination cell is an experimental prototype of size 0.5 m × 0.5 m that was described in detail in [1, 2].The steady state model of the plant, consisting of partial differential equations (PDEs) in two dimensions and boundary and initial conditions, was solved by means of the method of finite volumes [2, 19]. The latter method respects the final conservation of the balance equations opposite to the simpler finite differences method which does not guarantee it. The steady state model was validated experimentally by Ben Jabrallah et al. [17]. The dynamic model adds the time dimension to the steady state model and is more difficult to solve because the convergence should be verified at each time step. Thus, as the accuracy of the open loop results is not of prime importance and some limited error is acceptable, it was chosen to perform a spatial discretization according to the method of lines so that the partial differential equations with respect to time and space result in a large system of first-order ordinary differential equations (ODEs). This system of ODEs was written as a tridiagonal system according to the horizontal and vertical directions. To solve it, a double sweeping method was adopted which presents the advantage to rapidly propagate the influence of the boundary conditions towards the inside of the cell and to accelerate the convergence. With the solution being searched at each time step, the boundary conditions were updated at each iteration [1]. As the dynamic model is deduced from the steady state model by use of first principles, it is also assumed that the steady state validation is sufficient for its dynamic purpose. A Fortran program was thus developed which simulates the dynamic behaviour of the cell with respect to the variation of its different manipulated inputs. This model of the plant will serve as the reference along the study.


## 3. Reduced Dynamic Model of the Cell

The model of PDEs with initial and boundary conditions is already difficult to solve. Thus, it is inadequate for control purposes and, for that reason, a reduced dynamic model is necessary. However, the main PDE model will constantly be used to represent the plant and provide the real outputs. Several inputs influence the behaviour of the cell: the heat flux densityq˙fc, the water feed flow rate Ff, the water feed temperature Tf, and the condensation temperature Tc.(26)u=TfFfTcq˙fcT,where the three first inputs are manipulated variables and the fourth input is a disturbance. The scalar controlled output is the distillation flow rate Fd.(27)y=Fd,so that the cell can be outlined as in Figure 2.Figure 2
Black-box model of the cell.Among advanced control possibilities, Model Predictive Control (MPC) has been used because of its multivariable character and the possibility of handling constraints and its intrinsic robustness. As linear MPC will be used to control the plant and optimize its distilled water production, it is desirable to obtain a linear model of the plant around some operating point. The chosen steady state operating regime corresponds to the following inputs:q˙fc=600 W·m−2, Tf = 25°C, Ff=1.46 kg·h−1, and Tc = 20°C. Furthermore, MPC is used under the Quadratic Dynamic Matrix Form (QDMC) [20] based on step response coefficients. Thus, around the steady state operating regime, step moves of 10% of the manipulated inputs, each at its turn, are imposed and the corresponding dynamic variations of the controlled output are recorded. These nonlinear responses are then fitted by means of continuous first-order transfer functions (Table 1). Finally, the step response coefficients are deduced from the response of these first-order transfer functions to unit step inputs. The comparison (Figure 3) between the step response of the rigorous PDE model of the plant and the reduced model given by the transfer functions shows a good agreement. A remarkable feature of this study is that the full PDE model can be well represented by a system of first-order transfer function.Table 1
Identified first-order transfer functions (units are SI).
Input T f F f T c q ˙ f Transfer function 79.5 × 1 0 - 4 4500 s + 1 - 0.0624 4000 s + 1 - 9 × 1 0 - 5 2800 s + 1 9.26 × 1 0 - 4 4000 s + 1Figure 3
(a) Step response ofFd versus Tf of the PDE model of distillation cell (circles) and of continuous transfer function (continuous line) for ΔTf=4°C. (b) Step response of Fd versus Ff of the PDE model of distillation cell (circles) and of continuous transfer function (continuous line) for ΔFf=0.4 kg·h−1. (c) Step response of Fd versus Tc of the PDE model of distillation cell (circles) and of continuous transfer function (continuous line) for ΔTc=2°C. (d) Step response of Fd versus q˙f of the PDE model of distillation cell (circles) and of continuous transfer function (continuous line) for Δq˙f=50 W/m2.
(a) (b) (c) (d)


## 4. Solar Concentration by Heliostats and Experimental Setup

In the previous discussion about distillation cell modelling, the heating flux could be of any type, even electrical at the laboratory scale [2]. Of course, in real practice, solar power must be used. Instead of a simple direct use, heliostats tracking the sun position can be used to concentrate the solar rays on the heating plate of the cell. A heliostat is a plane mirror, mechanically equipped with two degrees of freedom. Its position is defined by two incidence angles i and i′ that simultaneously depend on the azimuth and the sun elevation and also on its position with respect to the desalination cell defined by the two angles γ and ϕ (Figure 6). The solar heat flux density reflected by the heliostat depends on these four angles, on the solar heat flux density q˙f and on the state of the reflecting surface of the heliostat mirror. To increase the heat flux density reflected towards the cell, two independent operations are necessary, the multiplication of the number of heliostats and their optimal orientation in order to track the sun.In the proposed setup, three heliostats are used. In order to correctly control each heliostat, it is necessary to know its parameters and to determine its position with respect to the desalination cell, that is, its elevation angle (north-south) and its rotation angle (east-west). In a given place and at a precise instant, the sun position is known, generally by its azimuthθ and elevation or height h. From the classical relations of celestial mechanics [21, 22], the latitude, longitude, and altitude define the place; furthermore, with the day and hour, the declination of the earth δ and the clockwise angle ω can be calculated. By convention, it is admitted that the azimuth of the sun is positive during the morning and equal to 180∘ at solar noon. The addition of an additional aiming system would be beneficial for the exact tracking of the sun but was not installed. Two jacks are used per heliostat for the azimuthal move and the elevation. Each motor is equipped with a speed reducer, coupled to a jack with an endless screw. A sensor allows us to send impulses during the rotation of the motor. More details about heliostats and their control are given in Appendix A.


## 5. Model Predictive Control

Model Predictive Control (MPC) owes its good reputation, first in petrochemical industry [23] then in all industrial domains, to its capacity to handle multivariable systems, even of large size, while taking into account linear constraints on controlled variables and their moves and even possibly nonlinear constraints on controlled variables. MPC is stated under several forms using the same philosophy [24]. MPC lies on a model of the predicted outputs based on past and future inputs over a prediction horizon Hp. This model can be based on impulse response coefficients, step response coefficients, or state space approach. At each sampling time, it performs the following:(i)
The calculation of the predicted outputs based on past inputs over the prediction horizon.(ii)
The calculation of a reference trajectory to be followed over a prediction horizonHp.(iii)
The estimation of future disturbances which take into account the output measurements.(iv)
From the minimization of a quadratic criterion subject to hard and soft constraints, the vector of future manipulated input variations is calculated over a control horizonHc.(v)
Only the first element of the vector of future manipulated input variations is used to calculate the control law applied on the system at the sampling time.After implementation of the control law, the sequence of calculation is repeated at the next sampling instant according to the receding horizon methodology.Many synthesis articles [20, 25–32] and textbooks [33–35] have been published which show the development of MPC and its application in the industrial world [36]. Compared to most usual control possibilities, the main advantages of MPC are the multivariable character by means of the internal model based on step or pulse responses or state space model and the consideration of constraints [24]. In our Fortran code, the internal model based on step responses is used as follows:(28)yk+1=yss+∑i=1MhiΔuk-i+1+dk+1,where hi’s are the step response coefficients at instant i, yss is the steady state output, d is the disturbance, and M is the model horizon.The predicted output of the model is(29)y^k+l∣k=yss+∑i=l+1MhiΔuk+l-i︸effects of past inputs+∑i=1lhiΔuk+l-i︸effects of future inputs+d^k+l∣k︸disturbances. The predicted disturbance is considered constant and takes into account the measured output as(30)d^k+l∣k=d^k∣k=ymk-y∗k∣k.The quadratic criterion generally considered in MPC algorithms is the sum of two contributions, a term of performance related to the deviations of the controlled outputs and a term of energy related to the variations of the manipulated inputs(31)JΔu=∑j=1ny∑i=1Hpγj2y^jk+i∣k-rjk+i2+∑j=1nu∑i=1HcλjΔujk+i-12,where y^(k+i∣k) are the predicted outputs (ny=1) based on past and future input variations, r(k+i) the reference trajectory, and Δuj(k) the variations of the manipulated inputs (nu=4). Hc and Hp are the control and prediction horizons, respectively. γ and λj are tuning weights. This criterion is minimized subject to constraints which are summarised as hard constraints on manipulated variables(32)umin⁡≤u≤umax⁡,hard constraints on the moves of the manipulated variables(33)Δumin⁡≤Δu≤Δumax⁡,and soft constraints on the predicted outputs(34)ymin⁡≤y^≤ymax⁡.All these constraints can be gathered as a set of linear inequalities [20, 24, 33](35)BΔuk≤Ck+1.In the absence of constraints, the criterion (31) yields an analytical solution. In the presence of linear constraints (35), the criterion can be easily minimized numerically using a QP code. When nonlinear soft constraints on the outputs are present, the problem becomes nonquadratic and it is solved by nonlinear sequential quadratic programming with the code NLPQL [37]. A Fortran90 MPC code, based on step inputs, able to work with any number of manipulated inputs and controlled outputs, with all types of hard or soft constraints was used [24, 38].


## 6. Maximization of Water Production

The increase of water production requires working on two aspects, one dealing with the increase of the solar heat flux density reflected on the cell by automatic control of heliostat orientation, followed by application of MPC control of the water desalination cell.The global system formed by the desalination cell and the heliostats consists, from the point of view of control, of two distinct subsystems (Figure4). The first subsystem deals with the automatic individual control of the three heliostats.Figure 4
Principle of the desalination cell with heliostats.The second subsystem aims at providing a control able to maximize the production of distilled water. As the desalination cell is modelled by a distributed parameter system [1, 2], influenced by both manipulated inputs and disturbances, the choice of the control technique for maximizing the yield is decisive. The subsystems will be controlled separately but with an interaction represented by the total incident heat flux density q˙fc which is a controlled output for the heliostat subsystem but is a disturbance for the cell subsystem.


### 6.1. Automatic Positioning of Heliostats

In the case of the cell considered alone, without heliostat, placed in a vertical position and orientated towards the south, during the chosen day, July 21st, the solar flux density is represented in Figure5. It can be noticed that this configuration is not optimum as the solar rays present a variable angle with respect to the cell during the whole day. During a large part of the day, the incidence angle of the solar rays is large, which generates a smaller received flux. The optimum would be reached if the solar rays were perpendicular to the cell by means of a rotation in the azimuthal plane and another rotation in the elevation plane. With the cell being fixed, the latter solution is not optimal.Figure 5
Solar heat flux density: real (-) and received (- -) by the cell in a vertical position and orientated towards the south, without heliostat.Figure 6
Disposition of a heliostat in the azimuthal plane.In the present study, a new solution is proposed which consists in orientating the heating plate of the cell towards the north and in a vertical position. When heliostats are used, the cell is orientated towards the north in front of the mirrors that reflect the solar rays. In this way, it is possible to increase the flux incident to the cell by multiplying the number of mirrors. This is the technique used in solar towers.The solar flux,q˙f,i, incident on the plane mirror of the heliostat #1 depends on the azimuthal angle θ1 and on the elevation angle h of the sun. The heliostat is placed in the north-south direction in front of the cell with an angle γ between both horizontal planes supported by the two normal vectors to the cell and heliostat #1. From Figures 6 and 7, the expression of the incident flux on the plane mirror plan of heliostat #1 results in the following:(36)q˙f,i=q˙fcos⁡i1cos⁡i1′.These expressions of the incidence angles, i1 in the azimuthal plane and i1′ in the elevation plane, are(37)i1=180∘-θ12,i1′=h-γ2.The maximum incident flux is reached at solar noon and is equal to(38)q˙f,i=q˙fcos⁡i1′.The incident flux q˙fc to the desalination cell is(39)q˙fc=q˙f,rcos⁡γ.This value can be improved by increasing up to the possible maximum angle γ as it appears in Figure 7. An automatic control of the three heliostats is performed by electronic cards based on PIC 18F4620 microchip.Figure 7
Disposition of a heliostat in the elevation plane.


### 6.2. MPC Control of the Cell

The manipulated variables for the cell control (Figure2) are the brine feed temperature Tf, the condensation plate temperature Tc, and the brine flow rate Ff. On the opposite, the input variable q˙fc is considered as a disturbance as it mainly depends on the heat flux density q˙f. The controlled output is the produced distilled water flow rate Fd. As the objective of the plant is the maximization of distilled water quantity, which can be expressed as the integral of the instantaneous distilled water flow rate that is measured, it theoretically results in a problem of dynamic optimization. However, this latter is most often performed off-line as an open loop problem by methods issued from variational calculus such as Hamilton-Jacobi method, Pontryagin’s maximum principle, or iterative dynamic programming [24]. Moreover, with the PDE model with initial and boundary conditions, dynamic optimization is a difficult problem that only calculates the set points for future on-line control. Thus, a different manner of handling the maximization problem was preferable in the present context. Indeed, by use of physical thoughts about the plant behaviour, it is possible to adopt a strategy that finally will maximize the distilled water quantity rather than classically follow a reference trajectory. The temperature of the condensation plate must be maintained as low as possible, that is, the difference (Tf-Tc) must be as large as possible. In our case, the reference was chosen equal to 1.5 times the maximum physical flow rate that could be delivered by the cell, so that the reference was impossible to reach. Similarly, the feed flow rate Ff greatly influences the yield of the cell, and it is fixed to twice the distilled flow rate Fd considered as correct [39] and previously determined by laboratory experiments. The feed temperature can be increased by making the water circulate in pipes exposed to the sun. All these constraints are expressed as follows:(40)Tc+10≤Tf≤333K,0.5≤Ff≤5kg·h-1,293≤Tc≤303K,Tc≤Tf,Fd≤Ff.The experiments on the distillation cell are presented for a particular day, July 21st. The sampling period is taken equal to 900s. After experiments performed on the cell, the weights of the quadratic criterion that give rise to the maximum production were chosen equal toΓ=1 and Λ=0.0001; that is, the priority is given to the performance term, thus to the reference tracking. The prediction horizon Hp is equal to 20 and the control horizon Hc is equal to 3, as in general a value between 1 and 6 is recommended. The reflected flux density is considered as a disturbance in MPC control as it depends on climate conditions.It must be noted in Figure8 that the difference Tf-Tc is maximum, Tf being at the maximum constraint and Tc at the minimum constraint, which allows us to improve the production of drinkable water. It can be noticed (Figure 9) that the heat flux density resulting from the use of the three heliostats is three times larger than the heat flux density obtained when the cell is fixed, orientated towards the south, and perpendicular to the ground.Figure 8
Temperatures: feed temperature (-) and condensation temperature (-·).Figure 9
Different densities of incident heat flux to the cell: sum of the densities of reflected heat flux by the three heliostats (∗), density of the heat flux of solar rays (-), density of reflected heat flux by heliostat #2 (-·), density of reflected heat flux by heliostat #1 (···), and density of reflected heat flux by heliostat #3 (- -).The produced flow rate resulting from the reflection of the solar rays by the three mirrors is promising as the reflection is twice more important at noon (Figure10). It must be noted that whereas the solar heat flux density (Figure 5) is a bell-shaped curve close to symmetrical, the flow rate is much distorted due to the delay of the distillation cell with respect to the solar heat flux density.Figure 10
Outlet flow rates: flow rate obtained with the three 3 mirrors (··) and flow rate obtained if the cell is perpendicular to the ground and orientated towards the south (-).A comparison is made between a conventional cell orientated towards the south in the absence of heliostat and the proposed setup with the cell orientated towards the north and equipped with heliostats. The density of solar fluxq˙f,c that strikes a cell conventionally orientated towards the south and set vertically was presented in Figure 5. In the same Figure, the heat flux density of the day q˙f was represented. The density of the incident solar flux q˙f,c to the cell orientated towards the north and equipped with one heliostat is presented in Figure 11. It appears that the received flux densities are close in both cases.Figure 11
Heat flux densities:q˙f real (-) and q˙f,c incident to the cell (- -), orientated towards the north, and with one heliostat.


## 6.1. Automatic Positioning of Heliostats

In the case of the cell considered alone, without heliostat, placed in a vertical position and orientated towards the south, during the chosen day, July 21st, the solar flux density is represented in Figure5. It can be noticed that this configuration is not optimum as the solar rays present a variable angle with respect to the cell during the whole day. During a large part of the day, the incidence angle of the solar rays is large, which generates a smaller received flux. The optimum would be reached if the solar rays were perpendicular to the cell by means of a rotation in the azimuthal plane and another rotation in the elevation plane. With the cell being fixed, the latter solution is not optimal.Figure 5
Solar heat flux density: real (-) and received (- -) by the cell in a vertical position and orientated towards the south, without heliostat.Figure 6
Disposition of a heliostat in the azimuthal plane.In the present study, a new solution is proposed which consists in orientating the heating plate of the cell towards the north and in a vertical position. When heliostats are used, the cell is orientated towards the north in front of the mirrors that reflect the solar rays. In this way, it is possible to increase the flux incident to the cell by multiplying the number of mirrors. This is the technique used in solar towers.The solar flux,q˙f,i, incident on the plane mirror of the heliostat #1 depends on the azimuthal angle θ1 and on the elevation angle h of the sun. The heliostat is placed in the north-south direction in front of the cell with an angle γ between both horizontal planes supported by the two normal vectors to the cell and heliostat #1. From Figures 6 and 7, the expression of the incident flux on the plane mirror plan of heliostat #1 results in the following:(36)q˙f,i=q˙fcos⁡i1cos⁡i1′.These expressions of the incidence angles, i1 in the azimuthal plane and i1′ in the elevation plane, are(37)i1=180∘-θ12,i1′=h-γ2.The maximum incident flux is reached at solar noon and is equal to(38)q˙f,i=q˙fcos⁡i1′.The incident flux q˙fc to the desalination cell is(39)q˙fc=q˙f,rcos⁡γ.This value can be improved by increasing up to the possible maximum angle γ as it appears in Figure 7. An automatic control of the three heliostats is performed by electronic cards based on PIC 18F4620 microchip.Figure 7
Disposition of a heliostat in the elevation plane.


## 6.2. MPC Control of the Cell

The manipulated variables for the cell control (Figure2) are the brine feed temperature Tf, the condensation plate temperature Tc, and the brine flow rate Ff. On the opposite, the input variable q˙fc is considered as a disturbance as it mainly depends on the heat flux density q˙f. The controlled output is the produced distilled water flow rate Fd. As the objective of the plant is the maximization of distilled water quantity, which can be expressed as the integral of the instantaneous distilled water flow rate that is measured, it theoretically results in a problem of dynamic optimization. However, this latter is most often performed off-line as an open loop problem by methods issued from variational calculus such as Hamilton-Jacobi method, Pontryagin’s maximum principle, or iterative dynamic programming [24]. Moreover, with the PDE model with initial and boundary conditions, dynamic optimization is a difficult problem that only calculates the set points for future on-line control. Thus, a different manner of handling the maximization problem was preferable in the present context. Indeed, by use of physical thoughts about the plant behaviour, it is possible to adopt a strategy that finally will maximize the distilled water quantity rather than classically follow a reference trajectory. The temperature of the condensation plate must be maintained as low as possible, that is, the difference (Tf-Tc) must be as large as possible. In our case, the reference was chosen equal to 1.5 times the maximum physical flow rate that could be delivered by the cell, so that the reference was impossible to reach. Similarly, the feed flow rate Ff greatly influences the yield of the cell, and it is fixed to twice the distilled flow rate Fd considered as correct [39] and previously determined by laboratory experiments. The feed temperature can be increased by making the water circulate in pipes exposed to the sun. All these constraints are expressed as follows:(40)Tc+10≤Tf≤333K,0.5≤Ff≤5kg·h-1,293≤Tc≤303K,Tc≤Tf,Fd≤Ff.The experiments on the distillation cell are presented for a particular day, July 21st. The sampling period is taken equal to 900s. After experiments performed on the cell, the weights of the quadratic criterion that give rise to the maximum production were chosen equal toΓ=1 and Λ=0.0001; that is, the priority is given to the performance term, thus to the reference tracking. The prediction horizon Hp is equal to 20 and the control horizon Hc is equal to 3, as in general a value between 1 and 6 is recommended. The reflected flux density is considered as a disturbance in MPC control as it depends on climate conditions.It must be noted in Figure8 that the difference Tf-Tc is maximum, Tf being at the maximum constraint and Tc at the minimum constraint, which allows us to improve the production of drinkable water. It can be noticed (Figure 9) that the heat flux density resulting from the use of the three heliostats is three times larger than the heat flux density obtained when the cell is fixed, orientated towards the south, and perpendicular to the ground.Figure 8
Temperatures: feed temperature (-) and condensation temperature (-·).Figure 9
Different densities of incident heat flux to the cell: sum of the densities of reflected heat flux by the three heliostats (∗), density of the heat flux of solar rays (-), density of reflected heat flux by heliostat #2 (-·), density of reflected heat flux by heliostat #1 (···), and density of reflected heat flux by heliostat #3 (- -).The produced flow rate resulting from the reflection of the solar rays by the three mirrors is promising as the reflection is twice more important at noon (Figure10). It must be noted that whereas the solar heat flux density (Figure 5) is a bell-shaped curve close to symmetrical, the flow rate is much distorted due to the delay of the distillation cell with respect to the solar heat flux density.Figure 10
Outlet flow rates: flow rate obtained with the three 3 mirrors (··) and flow rate obtained if the cell is perpendicular to the ground and orientated towards the south (-).A comparison is made between a conventional cell orientated towards the south in the absence of heliostat and the proposed setup with the cell orientated towards the north and equipped with heliostats. The density of solar fluxq˙f,c that strikes a cell conventionally orientated towards the south and set vertically was presented in Figure 5. In the same Figure, the heat flux density of the day q˙f was represented. The density of the incident solar flux q˙f,c to the cell orientated towards the north and equipped with one heliostat is presented in Figure 11. It appears that the received flux densities are close in both cases.Figure 11
Heat flux densities:q˙f real (-) and q˙f,c incident to the cell (- -), orientated towards the north, and with one heliostat.


## 7. Conclusion

The objective of the present work is to maximize the production of distilled water by a solar desalination cell. The solar rays are reflected by three heliostats. The complete system to be controlled thus consists of two interacting subsystems, the distillation cell and the heliostats. The density of the heat flux incident to the cell depends on azimuth and elevation of the sun and on the incidence angles that depend on the orientation of the heliostats. Based on the rigorous model of the system formed by partial differential equations which represent the plant, a step response model of the system has been identified for future MPC control, resulting in a model of transfer functions with four inputs, three of them being manipulated and the fourth one a disturbance. The maximization of the production of drinkable water is performed by MPC with automatic orientation of the heliostats. The density of reflected heat flux is more than three times larger than the solar flux density received by a vertical cell orientated towards the south and the production flow rate is two times larger than that of the previous vertical cell. This result is promising as it shows that a concentration of the solar rays by automatically controlled heliostats allows us to optimize water desalination with zero cost energy, and is moreover available a large part of the year.



---

*Source: 5924173-2017-01-16.xml*