Solar Energy Vol. 50, No. 6, pp. 469--476, 1993 Printed in the U.S.A.
0038--092X/93 $6.00+.00 Copyright © 1993 Pergamon Press Ltd.
NUMERICAL ANALYSIS OF GREENHOUSE-TYPE SOLAR STILLS WITH HIGH INCLINATION A. PALACIOand J. L. FERNANDEZ Institute of Engineering,UNAM, Ciudad Universitaria, Coyoacan 04510, Mexico Abstraet--A numerical analysishas been carried out to evaluate the technical feasibilityof singleeffect solar stills with a cover slope of up to 60 °. The aim of the research is to establishthe overall distillate productivity and the relative importance of the heat and mass transfer mechanisms when water diffusion is minimal and convection dominates. Experimental measurements of the distillate obtained from an already working unit are related to the numerical results by means of an "equivalent" Prandtl number that simplifiesthe analysis, allowing for the moist air in the inside of the cavity to be represented as a perfect gas.
1. INTRODUCTION Single effect solar stills have been well studiedand are currently being marketed in the dry regions of NW Mexico, where high-quality drinking water is a rare commodity and the solar option is economically viable, as proved by Sunwater T M solar stills. Most of these stills rely on diffusion to transport water from the hot basin to the colder glass cover. It is well known that a maximum distillate yield will result from a minimal air volume contained within the greenhouse envelope, attaining nearly ideal conditions for water diffusion. This design criterion results in generally narrow stills (about 1 m wide) which in turn makes the building of large solar stills a labor-intensive task. The largest height generally accepted for the greenhouse roof above the water surface is about 0.30 m[ 3 ], beyond which distillate yield diminishes abruptly. McCracken has successfully shown that a still with a maximum height of about 0.10 m will yield record production. The drawback is that stills must be small to perform efficiently. This article deals with stills that were developed to assess the feasibility of rising the height of the roof apex to enable the general cross section of the still to be enlarged. A pronounced slope is preferred to permit a minimal buildup of desert dust, a pervasive element in the region of interest. The experimental results can be coupled with a computer model that, in turn, sheds more light on the "ideal" cavity geometry that would enhance natural convection, thus partially compensating the loss in diffusion that stems from the largerthan-usual air volumes. 2. EXPERIMENTALINSTALLATION A solar still, 40 m long, 1.60 m wide, as measured at the shallow (2.5 cm) water basin bottom, and 0.46 m high at the middle (inside measurement) was built at the Centro de Investigaciones Biologicas near La Paz, Baja California Sur. The cover slope was approximately 24 °. A slight lengthwise slope, of about 1%, ensured that sea water would flow from the upper end to the lower one, while transversal ribs placed at the
basin bottom slowed down the water runoff and ensured the overall wetting of the basin bottom. Sea water was fed at a constant rate of about 25 l/h, as provided by a constant head tank, which in turn was filled with sea water by means of a wind pump. The installation is shown in Fig. I. Several ambient and performance data were gathered during a period of several weeks, during which the output of distillate was measured and related to other operation parameters. As expected [ 5 ], the distillate production didn't vary appreciably with respect to input water or ambient temperature, and occasional water feed irregularities didn't seem to have a measurable impact on distillate yield. Solar irradiance, on the other hand, seemed to be determining the rate of distilled water production. This same behavior has also been observed at the site in solar stills of other design characteristics [ 6 ]. Various other measurements were also taken, but their relevance to the present discussion is scarce and will thus be omitted. However, basin water temperature, glass temperature, and other basic performance indicators will be needed in the following sections. Temperature measurements were made with copper constantan 30 ga thermocouples with the aid of an automatic data recording system based on a PC computer. Temperatures are accurate within 0.3°C when compared with reference thermometer readings. A typical daily production cycle is shown in Fig. 2. A clear day operation during winter is registered by a calibrated pyranometer (global irradiance over a horizontal surface) and compared to the amounts of distillate accumulated during 15-min periods. The distillate was recorded manually by means of a calibrated container which was fed by two polyvinyl chloride (PVC) hoses attached to the lowermost end of each of the distillate collecting channels. The thermal inertia of the system is seen to displace the instant of maximum yield about 2 h after midday. Production stops for practical purposes 3 h after sunset and starts about l h after sunrise. The accumulated yield of the day shown is 223 l, which amounts to slightly less than 3.5 l / m 2 of water basin inside surface. A McCracken solar still operating
A. PALACIOand J. L. FERNANDEZ distribution inside the solar still. Additional simplifying assumptions and boundary conditions are presented in section 4.
Fig. 1. Photograph of the Centro de lnvestigacionesBiologicas (CIB) large channel solar still. nearby produced, under identical circumstances, nearly 6 l / m 2 per day. These stills have been in operation for over two years and have already undergone several external and internal cleaning and maintenance procedures. It is apparent, although it cannot be properly quantified, that maintenance for the larger still, where an operator can crawl inside, is quicker and cheaper, and probably more effective. A fact that is difficult to measure, but that also becomes apparent, is that dust buildup can reduce distillate production, and therefore the shallower still needs frequent external dusting, whereas the dust rolls by itself down the tall still.
3.1 Coordinates and variables The general problem to be considered is steady, twodimensional, axisymmetric, laminar flow of a mixture enclosed in a cavity. A schematic diagram of the system indicating the geometrical dimensions is shown in Fig. 3. In order to represent the still geometry in an adequate manner, a general curvilinear coordinate system is employed for which the independent variables in the vertical and horizontal directions are y and x, respectively. According to the symmetry of the problem, just half of the domain was considered for the solution. Figure 4 shows the grid generated to discretize the computational domain, with 30 cells employed in the x direction and 25 in the y direction. This grid was found satisfactory after a preliminary study of grid independence. The main dependent variables that characterize the present flow situation are: the velocity component in the vertical direction, v, the velocity component in the horizontal direction, u, the static temperature, T, and the turbulent kinetic energy and its dissipation rate, K and ~, respectively. These variables are complemented by two main auxiliary variables, namely the static pressure, p, and the density, p. 3.2 Conservation equations The set of partial differential conservation equations describing the flow examined herein can be compactly represented as follows: 3.2.1 Continuity equation. V.(p~') = o
where ~¢ represents the velocity vector and V is the vector differential operator. 3.2.2 Momentum equations. @
3. M A T H E M A T I C A L
V'(pufC)- V'(#gradu)=--~x+ Ffx
@ V . ( p v V ) - V . ( # grad v) = - ~y +
This section includes the basic equations that must be solved to know the velocity field and the temperature 8 •
xX X xX XXXox X X
Dote: 2 4 / J a n / 1 9 9 1 Total distillate produced =2231
• • x
X× 400 9
I 17 T i m e , hr
Fig. 2. Typical daily production/solar irradiance.
° °* 19
Greenhouse-type solar stills
3.2.4 Turbulence transport equations. For the calculations presented in this study, the turbulent viscosity #t is determined from the two-equation k - ~ model. In this model, the turbulent viscosity is computed from the local values of the turbulent kinetic energy, k, and its dissipation rate, ~, as follows:
#t = C , pk2/~
Fig. 3. Schematic diagram with dimensions. where grad is the gradient of the associated variable, and Ff is the body force due to buoyancy. Equations (2) and (3) represent the momentum balance in the coordinate directions x and y, respectively. The term referring to the buoyancy force may be expressed as follows:
Ff= ~(p - Po)
being the gravitational vector, and po a reference density value. Employing Boussinesq approximation, eqn (4) may be modified including the definition of the thermal expansion coefficient,/3, to yield:
Ff = -pop,[3( T -
where To is the reference static temperature. 3.2.3 Energy equation.
V = (-~pgrad T) = l~¢
where x is the thermal conductivity, Cp is the specific heat capacity at constant pressure, # is the dynamic viscosity, and cI, is the dissipation function. The latter term, which represents the mechanical dissipation effects, is neglected in the energy equation on the grounds that the Brinkman number defined as u V Z / ( x A T ) is very small (of the order of 10-6). This number provides a measure of the extent to which viscous heating is important relative to the heat flow resulting from the impressed temperature difference (see, e.g., Bird et al.).
where (7, is an empirical coefficient that is assigned the standard value of 0.09 [ 8 ]. The turbulent kinetic energy and its dissipation rate are obtained by solving the following modeled transport equations:
V - (pfCk) = V - ([# + #t/ak]grad k) + P k - p~
V - (pVe) = V - ([~ + ~,/z,]grad ~) + ~k(C~,Pk-- pC2,~)
and Pk is the volumetric stress-production rate of k. The turbulence-model coefficients aK, a,, C1,, and (72, are assigned the standard constant values of 1.0, 1.314, 1.44, and 1.92, respectively[ 8 ]. At solid walls, the wellknown wall-function approach outlined by Rodi[8 ] is adopted, which in effect means that the boundary conditions are not specified fight at the wall but at a point outside the viscous sublayer, where the logarithmic law of the wall prevails and the turbulence can usually be assumed to be in local equilibrium. 3.3 Solution method 3.3.1 The finite-volume equations. Each of the transport equations can be written in terms of a single general equation as follows, (p~b), + V . ( p V , ) = V.(F~V,) + S~
wherein ¢ is the dependent variable (u, v, T); 1~, is the exchange coefficient for ¢; S , the source of ~ per unit volume, and the subscript t denotes a partial derivative with respect to time. This form is well suited
Fig. 4. Domain discretization. Typical 35 × 25 cells grid.
A. PALACIOand J. L. FERNANDEZ
for numerical procedures employing the finite-volume method of discretization, such as that embodied in the PHOENICS computer program[9 ]. For the present calculation the first term, that is the transient one, vanishes for all the equations. The solution domain is subdivided into a number of control volumes, each associated with a grid point. The scalar variables and pressure are stored at the grid points, while the velocities are stored at staggered locations that lie between the pressure nodes. The control volumes for the velocities are staggered in relation to the control volumes for scalar variables. The finite-volume equations for each variable are derived by a control volume approach with the aid of assumptions about the distribution of the variables between nodes centered in each control volume. For the convection terms, all fluid properties are assumed uniform over cell faces and, except in respect of the velocities for which the face-center values are stored, the values prevailing at the cell face are determined by using upwind interpolation (see, e.g., Courant et al. [ 10 ] ). In diffusion terms, the property gradients and the transport properties that they multiply are uniform over cell faces. The gradients are based on the supposition that the properties vary linearly, and the transport properties are arithmetic averages of those on either side of the cell faces. In source terms, the nodal values are supposed to prevail over the whole of the cell volume. For any dependent variable, the partial differential equation is represented by a coupled set of algebraic equations of the form
Ap(~p = ~ Aicbi + S~
where the Ai's are the coefficients expressing the influence of convection and diffusion, S, is the integral source term, and the summation sign ~ stands for summation over the neighboring nodes. The effects of convection and diffusion are weighted using the hybrid formulation of Spalding, which cuts out the diffusive flux for cell Peclet numbers in excess of two. A complete statement of the PHOENICS finite-volume equations is given by Rosten and Spalding. 3.3.2 Solution procedure. The solutions are obtained with the computer code PHOENICS of Ludwig et al.. The solution algorithm is based upon the iterative "guess-and-correct" procedure of Patankar and Spalding[ 13 ], but is modified in accordance with the SIMPLEST algorithm of Spalding [ 14]. In general, the solution is obtained by sweeping repeatedly through the calculation domain, solving successively at each slab of cells for each dependent variable in turn. First, the equation of state is solved and then a modified version of Stone's [ 15 ] strongly implicit method is used to solve, in turn, the equation for T. Next, the velocities are obtained by solving the momentum equations using the old iterate pressures. The solution for the velocity variables proceeds by way of the SIMPLEST procedure. Then, continuity is enforced by solving a pressure-correction (p') equation
(see Patankar and Spalding [ 13 ]), which determines the required adjustments to the velocities and pressures. In the calculations, a compressible form of the p'equation is employed. The equation is developed by utilizing the equation of state, source terms from the finite-volume continuity equations, and coefficients from truncated forms of the momentum equations. The p' equation is solved in a "whole-field" manner by using a simultaneous procedure, which is also similar to Stone's method. The whole process is then repeated until the solution converges. It is worth mentioning here that the PHOENICS finite-volume equations are expressed in "correction" form prior to their solution, thus:
Ar4~p = ~ A,ch~ + Rp
where Rp is the residual defined by
Rv = Z A,¢~* - Avck~ + S,
and the 4~*'s are the in-store values of the ~b's,and the ~'s are the corrections that must be supplied in order to make all the equations balance. The advantages of this approach have been discussed by Rosten and Spalding [ 12 ]. Like all iterative procedures, the PHOENICS solution procedure is not unconditionally stable. Convergence of the procedure is obtained through careful linearization of source terms together with appropriate relaxation of the flow variables. Two types of relaxation are employed, namely inertial and linear. Inertial relaxation may be applied to any dependent variable 4~, by adding the following source term to the fight-hand side of the finite-volume equation
S,,, = r/(~bp,o,d-- ~bp)
where r/is the so-called inertia, defined by n = pVp/ty, where Vp is the cell volume and tf is a false time-step. A false time-step under relaxation is applied to the velocities, and the temperature. In addition to the dependent variables, other quantities, such as pressure, are relaxed by using linear relaxation, as follows: ~b = ot~bnew + ( 1 - ct ) ~bojd
( 15 )
where a is the linear relaxation factor.
4. N U M E R I C A L
4.1 Simplifying assumptions At the present stage of understanding the problem, it was decided that a first approximation to the "real" problem with a minimum of simplifying assumptions would be adequate. The problem to be solved mathematically was defined as steady-state, which would approach the real behavior around midday, and temperature boundary conditions were chosen accordingly.
Greenhouse-type solar stills Water basin temperature was chosen from measured values around noon, as well as air temperatures surrounding the still. The heat transfer between the cavity and the surroundings is calculated by means of a global heat transfer coefficient, U, that as derived from previous research, can be assigned values between 5 and 20 W / m 2 K for average operating conditions. These values are obtained considering that the resistance to heat flow, due to the internal film coefficient, is negligible when compared to the resistance due to conduction through the glass cover and the one due to the air surrounding the still. If a value of 0.76 W / m K is employed for the thermal conductivity of the glass (whose thickness is 6 mm), and the external convection coefficient is varied from 6 to 25 W / m 2 K (see, e.g., Roshenow et al.), the values for U lie between the range specified in the previous paragraph. The minimum value would be representative of a low velocity wind condition, whereas the maximum value would correspond to a very high wind velocity. The gas mixture inside the still was considered homogeneous, with a thermal conductivity k chosen in such a way that it would reflect the expected tendency of the water vapor to "slide" with respect to the air. Prandtl number for the mixture would then be smaller than for pure air or pure water vapor, and its value was found from iterations until an adequate calculation of the heat flow, as compared with the experimental measurements of the distillate (Fig. 2), was attained. In other words, assigning a low value to the Prandtl number was a simple practice to match the experimental results by artificially increasing the diffusion coefficient of the energy equation. Steady-state temperatures were chosen from experimental measurements as 95 °C for the water basin and 20°C for the surrounding air. It should be mentioned that this temperature for the water basin represents one of the highest values recorded in the apparatus. At peak production capacity, 2 to 3 h after noon, the solar still operating under the above mentioned conditions will have an average heat flux of about 0.3 k W / m 2 as measured from the condensate rate, referred to the water basin area. The heat flow from the con-
- i 0.52m/s
densate rate is determined by multiplying the mass flow rate of distillate by the corresponding heat of condensation. As it will be shown later, this heat flux is properly reproduced by the numerical model when a Prandtl number of O. 1 is used. For all cases presented, simulations considering separately both laminar and turbulent flow conditions were performed in order to assess the differences arising from the use of the turbulent equations in the mathematical model. These differences had their major impact in the calculation of the heat flux through the glass cover, as stated later in section 4.3. The figures presented in the remaining part of this article correspond to the simulations performed with the mathematical model. 4.2 Flow velocities and temperature distributions A typical flow pattern, corresponding to a global heat transfer coefficient U = 10 W / m 2 K, is shown in Fig. 5. The fluid flows upward along the central section of the cavity and descends along the cover (in this and the following figures only the fight half of the whole still cavity is shown ). Maximum flow velocities, of the order of 0.4 cm/s, are seen to occur at the right hand of the figure, where temperature gradients are largest, and along the vertical center line of the still, some 0.33 m from the water basin. Flow patterns at other U values are basically the same, with some variation in velocity magnitudes. A plot of the vertical velocity along the cavity center line for U values of 5, 10, and 20 W / m z K is shown in Fig. 6. A plot of the isotherms for U = 10 W / m 2 K is shown in Fig. 7. It can be noticed that the regions with highest temperatures are grouped toward the center of the cavity, and that temperature drops as the fluid flows upward. The lowest temperature is attained below the lowermost edge of the cover. Temperature contours are again similar for other U values within the stated range. Figure 8 shows the temperature distribution along the cover for the three values of U. The horizontal axis of the figure corresponds to the lowermost portion of the cover, and dis-
Min: 1.5874E-05 m/s
Max: 5.9935E-01 m/s
Fig. 5. Velocity flow field inside the solar still.
A. PALACIOand J. L. FERNANDEZ value of Pr, the numerical model cannot be used with confidence in still geometries and under conditions very different from the ones reported herein. However, some insight in the nature of the heat transfer mechanisms and their relation to the still geometry can be gained through extrapolation (see Section 4.3).
• U=5W/m2K + U=lOW/m2K U=20W/m2K
20 15 10 0.00
Fig. 6. Vertical centerline velocities dependence on heat transfer coefficient.
tance is measured horizontally along the water basin toward the midsection of the still. For the highest value of U, the temperature at the lower portion of the cover is almost equal to the ambient air temperature. As it was mentioned earlier, results in Figs. 5-8 were found for Pr = 0.1. A higher Prandtl number results in a diminished heat flux, and hence calculations depart from the experimental findings. It is interesting to show how the relative importance of the convective and diffusive heat transfer mechanisms evolves as the Prandtl number varies. Figure 9 illustrates this importance in terms of the thermal flux. For the calculation of the purely diffusive case, the convective terms that are present in the energy transport equation were turned off. Similarly, for the purely convective case, diffusive terms were deleted from the same equation. Notice then that the mere addition of the heat flux calculated with pure diffusion to the heat flux calculated with pure convection, does not yield the same result as the calculation with both mechanisms at work due to the interaction between mechanisms. The selection of a diminished value for Pr is justified solely, as stated before, on the basis of the experimentally observed heat flux and must be studied further. Perhaps a more elaborate justification of an equivalent conductivity for the single phase gas mixture approximation must be made, probably along geometrical considerations, as done by Kuehn and Goldstein [ 17 ]. Subject to the prevailing uncertainty regarding the
4.3 Heat flux dependence on cover slope In the case when the values of U = 10 W / m 2 K and Pr = 0.1 are retained, Fig. 10 shows the variation of the heat flux with cover inclination for the diffusiononly condition, for the case where convection exists without diffusion, and for the combined heat transfer mechanisms at work. It is seen that the relative importance of diffusion is small in all cases studied, and diminishes as the tilt angle grows from a value of 25 ° to an upper value of 60 ° . On the other hand, convection is seen to dominate throughout the slope range, and to increase with slope. It should be mentioned here, that the inclusion of the turbulence equations in the mathematical model had its major effects for the cases where the tilt angle took the largest values. Under these circumstances, the calculated heat flux was approximately 5% larger than that obtained for pure laminar flow conditions. The flow patterns were essentially the same as those obtained for laminar flow simulations, except for the appearance of additional weak vortices enclosed in very small secondary recirculation zones. These results are supported by the fact that the Rayleigh number based on the still height, is nearly 10 9 for the case where the tilt angle is assigned the lowest value of 24 ° , whereas for the angle of 60 ° it exceeds a value of 10 J0. Those results are not expected to be valid for lower inclinations, since an increased relative importance of diffusion will be coupled with a lower value of Pr. Indeed, if the condensate peak heat flux is measured in the McCracken still, a value of about 500 W / m 2 is found. This result is not adequately reproduced by the numerical model unless the value of Pr is further lowered. Figure 11 shows the results calculated for the McCracken still, retaining the value of U = 10 W / m 2 K and the boundary temperatures employed in the previous figures, but solving for the McCracken geometry, which has a cover tilt of about 4 °. The
Fig. 7. Isotherms distribution inside the cavity.
Greenhouse-type solar stills 55.0
= u= ~W/~K +
U= I O W / m Z K U = 2 0 W/m 2 K
200 150 25.0
,,+.--v" .,-..4--~--~ . IO0
.4.l,-a,- ~vr-.x,-.~ ~ 150 O0
45 Angle of inclination,
Fig. 8. Temperature variation along the still cover.
Fig. 10, Heat fluxdependenceon coverslope and heat transfer mechanisms. independent variable in the graph is the Prandtl number. The figure shows that in this case the diffusive effect is dominant. The calculated results match the measured heat flux at a Prandtl number between 0.01 and 0.03, but underestimate the heat flux by about 30% at P r = O . I . 5. CONCLUSIONS The results reveal that the vertex of a solar still can be higher than the 0.3 m or so (as suggested by previous authors) without incurring in the substantial loss of distillate yield that was expected. However, as the still gets taller, convection rapidly displaces diffusion as the major heat and mass transfer mechanism, and therefore the cavity geometry must be adapted to facilitate convection vortices to develop. In those tall stills where convection dominates, usually when still width exceeds 1.5 m and slope is steeper than 25 ° , flow velocities and temperature distributions can be studied numerically with the approximation of a two-dimensional cavity filled with a
homogeneous gas, if the assumption of an increased thermal conductivity is accepted. The numerical model can then be employed to design further physical experiments to validate temperature and velocity fields, and eventually identify "optimum" geometries. These findings seem to announce a brighter future for large-section solar stills, if their slope is so chosen as to facilitate the development of fast-moving natural convection vortices. Natural convection can compensate substantially, although apparently not completely, the loss in heat and mass transfer by diffusion. This means that large-section solar stills can be built to an advantage if their initial cost can be made about 30% lower than the cost of small-section solar stills, since that figure is the relative loss in overall efficiency incurred when diffusion-dominated heat and mass transfer is inhibited by the larger amounts of air volume contained in tall stills. The results of this study reveal that large-section stills can be advantageous in several aspects: manu-
350 o,I E
~--~ PRANDTL=O.? ~
200 "r 150
Heat transfer m e c h a n i s m s
Fig. 9. Heat flux through the glass, consideringthe effect of the heat transfer mechanisms.
A. PALACIOand J. L. FERNANDEZ
•+ G,!bo, Diffusive
ProndfNumber l 0.07
Fig. 11. Effect of Prandtl number on McCracken still.
facturing processes that depart from traditional artisan practices can be devised; stills capable of a c c o m m o dating an operator to inspect and clean the inside when needed can be tested in the future; a n d industrial-scale stills can now be explored. T h e added b o n u s is that with a larger slope on the roof sections, desert dust accumulates less, a factor that must be evaluated in each e n v i r o n m e n t from the e c o n o m i c viewpoint.
Acknowledgments--The experimental work reported herein was performed by Mr. M. A. Porta-Gandara, Director of Technological Design at the Centro de Investigaciones Biologicas of Baja California Sur, and his collaborators. The authors acknowledge the able help of Mr. Norberto Chargoy of the Institute of Engineering, UNAM, and the facilities given by CHAM de Mexico to perform the numerical simulations. Many valuable comments were received from Mr. Horace McCracken, President of McCracken Solar Company in AIturas, California. NOMENCLATURE A Cp C1, C2, C~ Ff fg k p Rp S, T t U u
coefficient in the algebraic equations specific heat capacity at constant pressure constant for the turbulence model constant for the turbulence model constant for the turbulence model body force gravitational vector turbulent kinetic energy static pressure residuals source term per unit volume static temperature time global heat transfer coefficient x-direction velocity component velocity vector , v y-direction velocity component x horizontal coordinate direction y vertical coordinate direction
Greek symbols a linear relaxation factor 3 thermal expansion coefficient dissipation rate of turbulent kinetic energy
q~ ,I, F n K u p V
dependent variable dissipation function exchange coefficient inertia term for false time relaxation thermal conductivity dynamic viscosity density vector differential operator
Subscripts 0 reference value i neighboring nodes p cell node t turbulent
1. P. I. Cooper, Solar distillation--State of the art and future prospects, In: Solar energy and the Arab world (also see SunWorld 1413]), Pergamon Press, New York (1983). 2. H. McCracken, Distillate production of a simple direct solar still, SunWorld 14(3), 83-87 (1990). 3. R. N. Morse and W. R. W. Read, A rational basis for the engineering development of a solar still, Solar Energy 12, 5-17 (1968). 4. H. McCracken and M. A. Porta, Solar still performance, La Paz and Alturas, Proceedings of the XIV National Meeting on Solar Energy (Mexico), Mexican Solar Energy Society (ANES), La Paz, B. Cfa. Sur (1990). 5. G. O. G. L6f, Solar distillation, In: Principles of desalination, 2nd ed., chapter 11, Part B, Academic Press, New York (1973). 6. P. I. Cooper, Digital simulation of transient solar still processes, Solar Energy 12, 313-331 (1969). 7. R. B. Bird, W. E. Stewart, and E. N. Lightfoot, Transport phenomena, John Wiley & Sons, New York (1960). 8. W. Rodi, Turbulence models and their application in hydraulics. A state of the art review, Book Publication, Delft, The Netherlands (1980). 9. J. C. Ludwig, H. Q. Qin, and D. B. Spalding, The PHOENICS reference manual for version 1.5, CHAM TR/200, CHAM Ltd, London, (1990). 10. R. Courant, E. Isaacson, and M. Rees, On the solution of non-linear hyperbolic differential equations by finite differences, Comm. Pure Appl. Math. 5, 243-245 ( 1952 ). 11. D. B. Spalding, A novel finite-difference formulation for differential equations involving first and second derivatives, Int. Z Num. Meth. Eng. 4, 551-559 (1972). 12. H. I. Rosten and D. B. Spalding, The PHOENICS equations, CHAM TR/99, CHAM Ltd, London (1987). 13. S. V. Patankar and D. B. Spalding, A calculation procedure for heat, mass and momentum transfer in threedimensional parabolic flows, Int. J. of Heat and Mass Transfer 15, 1787-1806 (1972). 14. D. B. Spalding, Four lectures on the PHOENICS computer code, CFD/82/5, CFDU, Imperial College, University of London, London (1982). 15. H. L. Stone, Iterative solution of implicit approximation of multi-dimensional partial differential equations, SlAM J. Num. Analysis 5, 530-533 (1968). 16. W. M. Roshenow, J. P. Hartnett, and E. N. Ganic, Handbook of heat transferfundamentals, 2nd ed., McGraw Hill, New York (1985). 17. T. H. Kuehn and R. J. Goldstein, "An experimental and theoretical study of natural convection in the annulus between horizontal concentric cylinders, J. Fluid Mech. 74(4), 695-719 (1976).