Optimal Design of a Solar Desalination Unit with Heliostats

1Laboratoire des Systèmes Avancés (LSA), Ecole Polytechnique de Tunisie, BP 740, 2078 La Marsa, Tunisia 2Faculté des Sciences de Bizerte, 7021 Zarzouna, Tunisia 3Laboratoire d’Energétique et des Transferts Thermique et Massique de Tunis, Faculté des Sciences de Tunis, Campus Universitaire, 1060 Tunis, Tunisia 4Laboratoire Réactions et Génie des Procédés (LRGP), CNRS-ENSIC, Université de Lorraine, 1 rue Grandville, BP 20451, 54001 Nancy Cedex, France


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  .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][10][11][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][15][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 Section 1, 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.

First-Principles Dynamic Model of the Distillation Cell
The distillation cell, later denoted as the plant, illustrated in Figure 1, 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 , spatial coordinate in the horizontal direction, and , spatial coordinate in the vertical direction, are from left to right in the solid domain, heated plate, from left to right in the liquid film domain,  ≤   ≤ , from left to right in the gas domain. (1) 2.1.Heated Metal Plate.The conductive heat transfer through the metal plate (subscript "") heated at the outside by solar radiation follows Fourier law as where   is the spatial coordinate in the transverse direction of the solid plate and  the vertical direction.

Liquid Water Film.
The liquid water film (subscript "") falls on the internal surface of the heated plate.Its behaviour is described by momentum and energy balances as where   is the spatial coordinate in the transverse direction of the liquid water film.

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 "" of gas is omitted for convenience in all variables , , V, , and .The movement of the gas mixture in the cavity is described by the following conservation equations: Momentum equations Energy equation Transport of the diffusing species It is assumed that the two horizontal faces of the cell are adiabatic.
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 Similarly, the heat flux must verify the continuity equation , ∀, ∀.
The adherence condition of the water film is 2.4.3.In the Water Film.The temperature of the water film at the inlet of the cell is given by 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, The continuity equation for the shear stresses can be written as The heat flux also verifies the continuity equation where the second term in the right-hand stands for the liquidgas 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].
where   is the mass fraction of water vapor at the interface and D is the diffusion coefficient of water in air.The mass flux density ṁV exchanged at the liquid-gas interface results as [17] ṁV = (    ) 2.4.5.In the Gas Phase.It is assumed that the horizontal walls are perfectly adiabatic.
On the horizontal walls, the adherence condition is verified as follows: On the condensation plate, the adherence condition is verified as follows: On the condensation plate, the transverse velocity is deduced from Fick law similarly to (17) as follows: On the condensation plate, a uniform temperature is imposed and consequently the corresponding equilibrium concentration results.
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 density ṁΦ is obtained by application of Fick law at the contact surface between the vapor and the condensation plate as follows: where   is the mass fraction of water vapor at the condensation plate, D the diffusion coefficient of water vapor in air, and  the concentration of water vapor.The distillation flow rate   of the cell proportional to the concentration gradient of the diffusing species at the interface between the gas phase and the condensation plate is where  is the length of the film.

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.

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 density q   , the water feed flow rate   , the water feed temperature   , and the condensation temperature   .where the three first inputs are manipulated variables and the fourth input is a disturbance.The scalar controlled output is the distillation flow rate   .

𝑢 = [𝑇
so that the cell can be outlined as in Figure 2.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   = 600 W⋅m −2 ,   = 25 ∘ C,   = 1.46 kg⋅h −1 , and   = 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.

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  and   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  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 ℎ.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.

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   .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 horizon   .
(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 horizon   .
(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][26][27][28][29][30][31][32] and textbooks [33][34][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: where ℎ  's are the step response coefficients at instant ,  ss is the steady state output,  is the disturbance, and  is the model horizon.
The predicted disturbance is considered constant and takes into account the measured output as 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 where ŷ( +  | ) are the predicted outputs (  = 1) based on past and future input variations, ( + ) the reference trajectory, and Δ  () the variations of the manipulated inputs (  = 4).  and   are the control and prediction horizons, respectively. and   are tuning weights.This criterion is minimized subject to constraints which are summarised as hard constraints on manipulated variables hard constraints on the moves of the manipulated variables and soft constraints on the predicted outputs All these constraints can be gathered as a set of linear inequalities [20,24,33] 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].

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 (Figure 4).The first subsystem deals with the automatic individual control of the three 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  which is a controlled output for the heliostat subsystem but is a disturbance for the cell subsystem.

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 Figure 5.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.
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 , , incident on the plane mirror of the heliostat #1 depends on the azimuthal angle  1 and on the elevation angle ℎ 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: International Journal of Chemical Engineering These expressions of the incidence angles,  1 in the azimuthal plane and   1 in the elevation plane, are The maximum incident flux is reached at solar noon and is equal to The incident flux q   to the desalination cell is 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.

MPC Control of the Cell.
The manipulated variables for the cell control (Figure 2) are the brine feed temperature   , the condensation plate temperature   , and the brine flow rate   .On the opposite, the input variable q   is considered as a disturbance as it mainly depends on the heat flux density q  .The controlled output is the produced distilled water flow rate   .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 (  −   ) 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   greatly influences the yield of the cell, and it is fixed to twice the distilled flow rate   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: 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   is equal to 20 and the control horizon   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 Figure 8 that the difference   −   is maximum,   being at the maximum constraint and   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.
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 (Figure 10).It must be noted that whereas the solar heat flux density (Figure 5) is a bellshaped 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.
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  flux q , 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  was represented.The density of the incident solar flux q , 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  real (-) and q , incident to the cell (--), orientated towards the north, and with one heliostat.

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.

A. Positioning of the Heliostats
In the present case, by means of three heliostats equipped with plane mirrors (Figure 12), the solar rays are concentrated towards the sun as shown in Figure 13.
These heliostats are based on an altazimuth, that is, a two-axis frame equipped with two motors.One motor allows fixing the azimuth of the heliostat and the other one allows making its elevation vary.Both motors for each heliostat are controlled by a programmable electronic card based on the Microchip PIC18F4620 microcontroller.This control adjusts the direct reflection of the solar rays on the cell during a whole day for all the seasons of the year.The mirror used to reflect the solar rays is plane.The losses by reflection of the incident flux outside of the heating plate of the cell are reduced by improving the focusing of the heliostats and the accuracy of the aiming.The knowledge of the geographic location of the heliostat (latitude and altitude) as well as the date and hour allows us to calculate the azimuth and the elevation of the sun.Furthermore, the adequate position of the heliostat is calculated to allow an optimum reflection of the solar rays on the cell.The reflected flux density is measured by a flux sensor with a thin wire of HFS-4 type coupled with a processor meter OMEGA DP41-E equipped with a series output that delivers the measured value to the central computer.
At night, heliostat #1 is automatically brought back in the north-south direction and positioned vertically so that the reflecting face is orientated towards south.Each morning at a precise time, according to an already prepared off-line table, the time is read from the real time clock; then the azimuth  1 and the height ℎ of the sun are calculated.
A central computer allows us on one side to monitor the positions of the heliostats by the measurement of the reflected flux density on the desalination cell, on the other side to control by MPC an electrovalve which acts on the feed flow rate of the cell.It must be noted that the central computer does not control the heliostats, but it receives the information about the angles  and   from each electronic card equipping the heliostat.Then, after simultaneous acquisition of the real heat flux density and the heat flux density incident to the cell, it verifies if a heliostat is ill positioned.
Figure 14 shows the variation of the azimuth angle  1 and of elevation angle ℎ of the sun during all the day of July 7th that are calculated by our code according to the date and the time of the day.
The variations of  1 and of   1 that are calculated along the day are presented in Figure 15.They are used to control the corresponding heliostat, which allowed us to obtain the density of the incident solar flux q , to the cell presented in Figure 11.Heliostat Symbols  at : Attenuation factor ℎ: Eleva tio na n gle  1 ,  2 ,  3 : Incidence angles of the three heliostats in the azimuthal plane  1 ,  2 ,  3 : Reflection angles of the three heliostats in the azimuthal plane   1 ,   2 ,   3 : Incidence angles of the three heliostats in the elevation plane   1 ,   2 ,   3 : Reflection angles of the three heliostats in the elevation plane : Normal unitary vector to the surface of a plane mirror q , : Density of the incident heat flux to a mirror (W⋅m −2 ) q , : Density of the reflective heat flux to a mirror (W⋅m −2 ) q , : Density of the incident heat flux to the desalination cell (W⋅m −2 ) : Angle of incidence of the cell : Position angle of a heliostat with respect to north-south direction of the cell : Reflective power of a mirror surface : Azimuthal angle.

Figure 2 :
Figure 2: Black-box model of the cell.

Figure 3 :
Figure 3: (a) Step response of   versus   of the PDE model of distillation cell (circles) and of continuous transfer function (continuous line) for Δ  = 4 ∘ C. (b) Step response of   versus   of the PDE model of distillation cell (circles) and of continuous transfer function (continuous line) for Δ  = 0.4 kg⋅h −1 .(c) Step response of   versus   of the PDE model of distillation cell (circles) and of continuous transfer function (continuous line) for Δ  = 2 ∘ C. (d) Step response of   versus q  of the PDE model of distillation cell (circles) and of continuous transfer function (continuous line) for Δ q  = 50 W/m 2 .

Figure 4 :
Figure 4: Principle of the desalination cell with heliostats.

Figure 5 :Figure 6 :
Figure 5: Solar heat flux density: real (-) and received (--) by the cell in a vertical position and orientated towards the south, without heliostat.

− 1 )Figure 10 :
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 (-).

Figure 13 :
Figure 13: Disposition of the three heliostats in the azimuthal plane.
The heat flux density q  corresponding to solar radiation is imposed on the external face of the heated metal plate of thickness   .

Table 1 :
Identified first-order transfer functions (units are SI).