- Email: [email protected]
Numerical Modeling of Solar Ponds
Available online at www.sciencedirect.com
ScienceDirect Energy Procedia 57 (2014) 2416 – 2425
2013 ISES Solar World Congress
Numerical Modeling of Solar Ponds M. C. Giestasa, J. P. Milhazesb, H. L. Pinab,* a
LNEG, Estrada do Paço do Lumiar, 1649-038 Lisboa, Portugal IDMEC/IST, Avenida Rovisco Pais, 1049-001 Lisboa, Portugal
b
Abstract A SGSP is a basin of water where solar energy is trapped due to an artificially imposed salinity gradient. In a SGSP three zones can be identified: the surface and bottom zones that are both convective and an intermediate zone in between which is intended to be non-convective. This zone acts as a transparent insulation and allows the storage of solar energy at the bottom where it is available for use. A numerical model where the SGSP dynamics is described in terms of velocity ݑ, pressure , temperature ߠ and salt concentration ߪ is presented. It is based on the Navier-Stokes equations for an incompressible fluid coupled to one advection-diffusion equation for ߠ and one advection-diffusion equation for ߪ . The fluid density ߩ is taken to depend on ߠ and ߪ and the Boussinesq hypothesis is adopted: the fluid density appearing in the LHS of the NavierStokes equation is supposed constant and equal to some reference value whereas it is assumed to be variable in the RHS. The space discretization of the governing equations is based on the respective weak formulations and the discretization employs finite elements with a pressure correction method used to decouple velocity and pressure. Integration in time is accomplished by a BDF (Backward Differentiation Formula) method with the above PDEs treated sequentially within each time step. A computer code was developed employing the finite element class library deal.II. Comparisons with available experimental results are made to validate this numerical model. © 2014 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license © 2013 The Authors. Published by Elsevier Ltd. (http://creativecommons.org/licenses/by-nc-nd/3.0/). Selection and/or peer-review under responsibility of ISES. Selection and/or peer-review under responsibility of ISES Keywords: Solar ponds, numerical modeling, finite elements, double advection-diffusion.
* H. L. Pina. Tel.: +351-965-683-057; Fax: +351-218-419-634. E-mail address: [email protected]
1876-6102 © 2014 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/3.0/). Selection and/or peer-review under responsibility of ISES. doi:10.1016/j.egypro.2014.10.250
M.C. Giestas et al. / Energy Procedia 57 (2014) 2416 – 2425
1. Introduction Nomenclature A B b g h L1, L2 n p σ θ t
Thermal diffusivity [m2/s] Salt diffusivity [m2/s] Buoyancy Gravity acceleration [m/s2] Convection coefficient [W/m2K] Domain length and height Unit normal vector Pressure [Pa] Salt concentration Temperature [K] Time [s]
Subscripts 0 f r a 1, 2, 3, 4
Initial Final Reference values Ambient values Left, right, bottom, top
u u Ω, ∂Ω α β γ ρ . Δ x
Generic scalar field Velocity field [m/s] Domain, boundary Thermal expansion coeff. [K-1] Salt contraction coefficient Kinematic viscosity [m2/s] Fluid specific mass [Kg/m3] Gradient Divergence Laplacian Space coordinates
Superscripts ż n
Time derivative Time instants
A Salt Gradient Solar Pond (SGSP) is a basin of water where solar energy is stored. A temperature gradient (hotter at the bottom and cooler at the top) is established in a natural way due to solar radiation absorption at the surface and the salt concentration gradient (denser at the bottom and lighter at the top) acting to prevent convective motions that would otherwise promote the return of the stored energy to the outside ambient destroying the SGSP very purpose. Thus a double advection-diffusion process occurs with the temperature and salinity fields making opposite contributions to the fluid density. The experimental, analytical and numerical studies of SGSP can be traced back to Kalecksinsky [16] and was pursued by Stommel [8], Weinberger [17], Turner [26] and Tabor [27]. Recently there was a renewed interest on SGSPs as a device coupled to desalination units or as a low cost seasonal thermal energy storage reservoir. The numerical solution of the governing equations that describe a SGSP has been attempted. Hull [3], Hawlader and Brinkworth [15] and Rubin [22] have applied a finite difference method while [4] has used a finite element technique. The SGSP stability – one of the key factors governing a SGSP performance – has been studied by several researchers [9], [6], [24], [11] and [12]. Some of those authors resorted in most cases to a linear perturbation theory and the results obtained have provided important information regarding the onset of the instabilities as well as the existence of several possible stable or unstable states that may arise. Weinberger [17] gave a mathematical formulation of a SGSP behavior. Meyer [19] developed a numerical model to predict the time dependent behavior of the interface between the convecting and the non-convecting zones. Panahi et al. [20] simulated the dynamics of the SGSP with a one-dimensional model. Angeli and Leonardi [1] and [2] studied the evolution of salt concentration profiles and investigated the salt diffusion and stability of the density gradient. Kurt [18] modeled the non-convective zone as a series of horizontal layers with the upper and lower convective zones modeled as single homogeneous layers. Mansour et al. [21] solved numerically the problem of transient heat and mass transfer and long term stability of a SGSP through a two-dimensional finite volume method.
2417
2418
M.C. Giestas et al. / Energy Procedia 57 (2014) 2416 – 2425
A comprehensive numerical model to simulate the SGPS dynamics is still lacking and the present paper is intended as a contribution to the development of such a tool. Thus, a rectangular cavity filled with a mixture of water and salt and heated at the bottom is considered. The fluid is treated as newtonian incompressible, heat conducting according to Fourier’s law and the salt diffusion obeys Fick’s law and are subject to an uniform gravitational field. Also the Boussinesq hypothesis is adopted. Due to the complexity of this time-dependent coupled system of equations it is extremely difficult to devise a numerical method of solution that performs equally well with respect to the several physical phenomena involved. Therefore the problem is separated into a sequence of simpler subproblems in correspondence with the underlying physics thus allowing to apply the most adequate numerical method to each subproblem. The present study calls for the splitting of the temperature, salinity and NavierStokes equations the later employing a pressure correction approach leading itself also to a decoupling of velocity and pressure. Integration in time is accomplished by a BDF(2) (Backward Differentiation Formula, second order) method with the above PDEs treated sequentially within each time step. A computer code was developed using the finite element class library deal.II (see [5]). To validate this numerical model some comparisons with available experimental results are presented. 2. The governing equations The governing equations are the standard conservation laws for continuum media. The water-salt solution is assumed to behave as a newtonian fluid conducting heat according to Fourier’s law and the salt diffusion is governed by Fick’s law. Also the Boussinesq hypothesis is adopted: the specific mass appearing in the left hand sides of the relevant equations is supposed constant and equal to some average value whereas in the right hand sides it is assumed to vary with temperature and salt. We collect here the governing partial differential equations as well as the boundary and initial conditions adopted for the case studied. Heat equation: μ௧ ߠ ሺ ڄ ܝሻߠ ൌ ܣȟߠπ ൈ ሺͲǡ ܶ ሿ ܖ ڄ ߠܣ ݄ሺߠ െ ߠ ሻ ൌ Ͳ μπ ൈ ሺͲǡ ܶ ሿ ߠሺܠǡ Ͳሻ ൌ ߠ π
(1) (2) (3)
μ௧ ߪ ሺ ڄ ܝሻߪ ൌ ܤȟߪπ ൈ ሺͲǡ ܶ ሿ ܖ ڄ ߪܤൌ Ͳ μπ ൈ ሺͲǡ ܶ ሿ ߪሺܠǡ Ͳሻ ൌ ߪ ሺܠሻπ
(4) (5) (6)
Salt equation:
Navier-Stokes system: μ௧ ܝ ሺ ڄ ܝሻ ܝൌ ߥȟ ܝെ ܊π ൈ ሺͲǡ ܶ ሿ ܝ ڄ ൌ Ͳπ ൈ ሺͲǡ ܶ ሿ ܖ ڄ ܝൌ Ͳ μπ ൈ ሺͲǡ ܶ ሿ ܝሺܠǡ Ͳሻ ൌ π
(7) (8) (9) (10)
with the buoyancy right hand side ܊given by (11) ܊ൌ ܴሺߠǡ ߪሻ ؠሺͳ െ ߙሺߠ െ ߠ ሻ ߚሺߪ െ ߪ ሻሻ where ߙ and ߚ are the expansion coefficients for temperature and salt and ߠ and ߪ are some reference values.
2419
M.C. Giestas et al. / Energy Procedia 57 (2014) 2416 – 2425
As can be seen, our model consists of a system of coupled partial differential equations comprising two time dependent advection-diffusion equations and the time dependent Navier-Stokes equation. This amounts to a Rayleigh-Bénard problem with two components (temperature and salt concentration). 3 Discretization The governing equations are solved by a sequential splitting method described in Fig. 1. Therefore the problem is decomposed into a sequence of simpler subproblems in correspondence with the underlying physics thus allowing to apply to each subproblem the most adequate numerical method. On the other hand this procedure is first order meaning that the splitting truncation error is only ࣩሺȟݐሻ (a general introduction to operator splitting techniques can be seen in [14] or [7]). Initializations: ߠ ൌ ߠ Ǣߪ ൌ ߪ Ǣܝ ൌ ܰ ൌ ሾܶ Ȁȟݐሿ FOR ݊ ൌ Ͳ TO ܰ DO Solve temperature equation for ߠ ାଵ Solve salt equation for ߪ ାଵ Solve Navier-Stokes equations for ሺܝାଵ ǡ ାଵ ሻ END Fig. 1: Program structure. 3.1 Time integration Time derivatives are approximated by the BDF(2) (Backward Differentiation Formula, second order) method with constant time step ȟݐ, that is, for ݊ ͳ ଵ ଷ ଵ ߠሶሺݐାଵ ሻ ൎ ቀ ߠ ାଵ െ ʹߠ ߠ ିଵ ቁ ߪሶ ሺݐାଵ ሻ
ൎ
௧ ଶ ଵ ଷ
ቀ ߪ
ାଵ
௧ ଶ ଵ ଷ
ଶ ଵ
െ ʹߪ ߪ ିଵ ቁ ଶ ଵ
(12)
ܝሶ ሺݐାଵ ሻ ൎ ቀ ܝାଵ െ ʹܝ ܝିଵ ቁ ௧ ଶ ଶ The BDF(2) method is known to be second order accurate and stiff A-stable (see [13] for details). As BDFs methods are not self starting a backward Euler method was employed for the first time step. that is, To compute approximations to the advective (non-linear) terms, the following linearly extrapolated values will be used (13) ߠכ ൌ ʹܝ െ ܝିଵ ǡߪכ ൌ ʹߪ െ ߪ ିଵ ǡܝ כൌ ʹܝ െ ܝିଵ The heat and salt equations are thus discretized in time as ଵ ଷ ାଵ ቀ ߠ ௧ ଶ ܖା
ଵ
െ ʹߠ ߠ ିଵ ቁ ሺܝ ڄ כሻߠ ାଵ ൌ ܣȟߠ ାଵ
(14)
ߠܣ ܖ ڄ ݄ሺߠ ାଵ െ ߠ ሻ ൌ Ͳ ଵ ଷ ାଵ ଵ ቀ ߪ െ ʹߪ ߪ ିଵ ቁ ሺܝ ڄ כሻߪ ାଵ ൌ ܤȟߪ ାଵ
(15) (16)
ܖ ߪܤା ܖ ڄൌ Ͳ
(17)
௧ ଶ
ଶ
ଶ
The Navier-Stokes system is discretized in time using a sequential operator splitting method to decouple velocity and pressure thus making this later variable explicit (see [10] for a detailed description). In the first substep, an intermediate (non divergence–free) velocity ܝ ෝ ାଵ is obtained by solving
2420
M.C. Giestas et al. / Energy Procedia 57 (2014) 2416 – 2425 ଵ
ଷ
ଵ
ቀ ܝ ෝ ାଵ െ ʹܝ ܝିଵ ቁ ሺܝ ڄ כሻܝ ෝ ାଵ ൌ ߥȟܝ ෝ ାଵ െ ͓ ܊כ ଶ ܖڄൌͲ ܝ ෝ Here ͓is an extrapolated pressure given by ସ ଵ ͓ൌ ߶ െ ߶ ିଵ ଷ ଷ and the buoyancy is computed explicitly by ܊ כൌ ܴሺߠכ ǡ ߪכ ሻ ௧ ଶ ାଵ
(18) (19) (20) (21)
In the second substep we correct to obtain a divergence–free velocity ܝାଵ and a pressure ାଵ by solving ଷ ሺܝାଵ െ ܝ ෝ ାଵ ሻ ൌ ሺ ͓ െ ାଵ ሻ (22) ଶ௧ ାଵ ൌͲ (23) ܝڄ (24) ܝାଵ ܖ ڄൌ Ͳ Applying the divergence operator to equation (22) we get a Poisson equation for the pressure ାଵ with homogeneous Neumann boundary conditions: ଷ ܝڄ ෝ ܖା (25) ȟሺାଵ െ ͓ሻ ൌ ଶ௧ ାଵ ൌͲ (26) μ ܖ The new pressure ାଵ is updated by ାଵ ൌ ߶ ାଵ and the new velocity ܝାଵ by ଶ௧ ሺ ͓െ ାଵ ሻ ܝାଵ ൌ ܝାଵ כ ଷ
(27) (28)
3.2 Space discretization For the space discretization, a Galerkin variational formulation is employed. The following finite element approximations will be used in the examples: (i) Temperature ߠ and salt ߪ: biquadratic elements (ܳଶ ); (ii) Velocity ܝand pressure : bilinear elements (ܳଵ ). 3.3 Stabilization For SGPS problems, equations (1) and (4) happen to be advection dominated. It is now well known that in this case Galerkin type methods lack stability and can return meaningless solutions usually displaying non physical oscillations or wiggles. It is therefore necessary to provide some form of stabilization in order to recover from these difficulties which is generally achieved by suitably modifying the underlying weak form. The method adopted was the isotropic artificial diffusion SOLD stabilization as described in [25]. 3.4 The linear algebraic systems The linear algebraic systems resulting from the space and time discretizations were solved by the conjugate gradient method with ILU(0) preconditioning (see for instance [23]).
M.C. Giestas et al. / Energy Procedia 57 (2014) 2416 – 2425
4. Results 4.1 Constants The constants values (in SI units) used to obtain the results were: Δ ݐ ൈ ͳͲିଶ ͷ ൈ ͳͲିଶ ିଵ ݄ଵ ͳͲǤͲ ͳ ൈ ͳͲ ݄ଶ ݄ଵ ͳǤͲͳͶ ൈ ͳͲଷ ݄ଷ ͵ʹͲǤͲ ͵Ǥ͵ͷ ൈ ͳͲିସ ݄ସ ͳͲͲǤͲ Ǥͷ ൈ ͳͲିସ ʹʹǤͻͷ ߠଵ ʹʹǤͻͷ ͲǤͲ ߠଶ ߠଵ ߠଷ ͷͲǤͲ ͲǤͺͲͷ ൈ ͳͲି ߠସ ʹʹǤͻͷ ͳǤ ൈ ͳͲି ߠ ʹʹǤͻͷ ͳǤͶ ൈ ͳͲିଽ ͶǤͻͳ ൈ ͳͲଷ ݔଶ ݔଶ ͲǤͺͷ ͲǤͲͷͳ െ ͲǤʹ ൬ͳ െ ൰ ǡ Ͳ ߪ ܮଶ ܮଶ ൞ ݔଶ ͲǤͲͲͳʹǡͲǤͺͷ ͳ ܮଶ
ܮଵ ܮଶ ߩҧ ߙ ߚ ߠ ߪ ߥ ܣ ܤ ܶ
4.2 Geometry and mesh The domain for the examples in Section 4 is a rectangle π ൌ ሾͲǡ ܮଵ ሿ ൈ ሾͲǡ ܮଶ ሿ as depicted in Fig. 2. Its boundary is μπ ൌ ڂସୀଵ Ȟ , where the Ȟ are the rectangle sides.
Fig. 2: Geometry and notation. The mesh comprises Ͳ ൈ ͳͶͲ ൌ ͻͺͲͲ elements and ͳͲͲͳͳ degrees of freedom for each component of the velocity and for the pressure, ͵ͻʹͳ degrees of freedom for the temperature and ͵ͻʹͳ degrees of freedom for the salt concentration. 4.3 Model validation A comparison of experimental temperature profiles and shadowgraph images is made against the numerical results. In the numerical case (NC) the vertical component of the density gradient was calculated in order to create a Schlieren image for the whole domain. This in turn made possible the comparison between the shadowgraph image and the Schlieren type. Temperature profiles were calculated at ൎ ͷ ൈ ͳͲିଷ from the right wall (see shadowgraph images). The ruler seen on the right side of the shadowgraph image gives an indication on the distance from the bottom in centimeters.
2421
2422
M.C. Giestas et al. / Energy Procedia 57 (2014) 2416 – 2425
Fig. 3-5 depict the shadowgraph images for the experimental case (EC) together with the temperature profile for a given time in seconds. Bellow we show the results of temperature profile and Schlieren obtained by the numerical model for comparison. Fig. 3 corresponds to an early state of the experiment. At time ݐൌ ͷͶͲs, a small but unstable lower convective zone (LCZ) is being developed. The upper convective zone (UCZ) in both NC and EC can also be observed and a second convective zone (CZ2) if beginning to form. At ݐൌ ͳͺͲ (see Fig. 4), one can notice that the separation between the LCZ and the gradient zones (GZ) is now sharply defined at ൎ ͵ ൈ ͳͲିଶ from the bottom. There is still a CZ2 but very unstable. In Fig. 5 we can verify that both in the NC and EC the CZ2 has vanished. The interface between the LCZ and GZ has is now at Ͷ ൈ ͳͲିଶ . Also, small movements can be observed on the UCZ due to heat loss from the surface. In all the cases shown above there is a reasonable approximation between the EC and NC. Fig. 6 shows a side-by-side comparison of the velocity magnitude vertical profile at a given time. In the EC the velocity was obtained using a PIV (Particle Image Velocimetry) technique. In both cases the vertical profile was taken along a vertical line at ൎ ʹ ൈ ͳͲିଶ m apart from the left wall.
(a)
(b) Fig. 3: (a) Temperature profile and shadowgraph image for EC. (b) Temperature profile and Schlieren image for NC. Time t=540 s
M.C. Giestas et al. / Energy Procedia 57 (2014) 2416 – 2425
(a)
(b) Fig. 4: (a) Temperature profile and shadowgraph image for EC. (b) Temperature profile and Schlieren image for NC. Time t=1680 s
2423
2424
M.C. Giestas et al. / Energy Procedia 57 (2014) 2416 – 2425
(a)
(b) Fig. 5: (a) Temperature profile and shadowgraph image for EC; (b) Temperature profile and Schlieren image for NC. Time t=4691 s
(a) (b) (c) (d) Fig. 6: Velocity magnitude profile for (a) EC; (b) NC at t=1860 s; (c) for EC; (d) NC for t=3000 s 5. Conclusions The results presented – both experimental and numerical – comprise a subset which are considered representative of the results available. Despite the small scale experimental setup it is believed that the information gathered can be useful for the full scale SGPS dynamics.
M.C. Giestas et al. / Energy Procedia 57 (2014) 2416 – 2425
The numerical model developed seems to be able to predict the main features of the physics involved. However it is still two-dimensional and its capacity to simulate the long term dynamics has yet to be demonstrated. The move to a three-dimensional model requires a particular attention to the computer implementation, namely code parallelization issues, and the employment of a robust time-integrator. Acknowledgments The work reported was supported by IDMEC/IST. J. P. Milhazes was supported by a FCT studentship grant. The computations described in this article were made using the finite element library deal.II (see [5]). References [1] C. Angeli and E. Leonardi. A one-dimensional numerical study of the salt diffusion in a salinity-gradient solar pond. Int. J. Heat Mass Transfer, (47):1–10, 2004. [2] C. Angeli and E. Lombardi. The effect of thermodiffusion on the stability of a salinity gradient solar pond. Int. J. Heat Mass Transfer, 48(21-22), 2005. [3] J. R. Hull. Computer simulation of solar pond thermal behaviour. Solar Energy, 25(33):33–40, 1980. [4] T. S. Jayadev and J. Henderson. Salt concentration gradient solar ponds. Modeling and Optimization. Colorado, 1979. [5] W. Bangerth, R. Hartmann, and G. Kanschat. deal.II — a general-purpose object-oriented finite element library. ACM Trans. Math. Softw., 33(4), 2007. [6] L. N. Da Costa, E. Knobloch, and N. O. Weiss. Oscillations in double-diffusive convection. Journ. Fluid Mech., 109(257):25–43, 1981. [7] I. Farago and A. Havasiy. Operator Splittings and their Applications. Nova Science Publishers, 2009. [8] H. Stommel, B. Arons, and D. Blanchard. An oceanographical curiosity: the perpetual salt fountain. Deep Sea Res., 3:152–153, 1956. [9] G. Veronis. Effect of a stabilizing gradient of solute on thermal convection. Journ. Fluid Mech., 34:315–336, 1968. [10] J. L. Guermond, P. Minev, and Jie Shen. An overview of projection methods for incompressible flow. Comput. Methods Appl. Mech. Engrg., (195):6011–6045, 2006. [11] M. Giestas, A. Joyce, and H. Pina. The influence of radiation absorption on solar ponds stability. J. Heat Mass Transfer, 39,18:3873–3885, 1996. [12] M. Giestas, A. Joyce, and H. Pina. The influence of non-constant diffusivities on solar ponds stability. J. Heat Mass Transfer, 40, 18:4379–4391, 1997. [13] E. Hairer and G. Wanner. Solving Ordinary Differential Equations II – Stiff and Differential-Algebraic Problems. Springer-Verlag, 1991. [14] W. Hundsdorfer and J. G. Verwer. Numerical Solution of Time-Dependent Advection-Diffusion-Reaction Equations. Springer-Verlag, 2003. [15] M. N. A. Hawlader and B. J. Brinkworth. An analysis of the non convecting solar pond. Solar Energy, 27:195–204, 1981. [16] A. V. Kalecsinsky. Ueber die ungarishen warmen und heissen kochsalzseen als natuerlich waermeaccumulatoren. Ann. Physik IV, 7:408. [17] H. Weinberger. The Physics of the solar pond. Solar Energy, 8(2). [18] H. F. Kurt, Halicir, and A. Binark. Solar pond conception - experimental and theoretical studies. Energy Conv. Manag., 41(9), 2000. [19] K. A. Meyer, D. P. Grimmer, and G. F. Jones. An experimental and theoretical study of salt-gradient pond interface behaviour. Progress in Solar Energy, 1982. [20] Z. Panahi, J. C. Batty, and J. P. Riley. Numerical simulation of the performance of a salt gradient solar pond. Transactions of ASME, J. of Solar Energy Enginnering, 105:361–374, 1983. [21] R. Mansour, C. Nguyen, and N. Galanis. Transient heat and mass transfer and long-term stability of a salt-gradient solar pond. Mech. Research Communications, 33:233–249, 2006. [22] H. Rubin, B. A. Benedict, and S. Bachu. Modeling the performance of a solar pond as a source of thermal energy. Solar Energy, 32:771–778, 1984. [23] Y. Saad. Iterative Methods for Sparse Linear Systems. SIAM, 2003. [24] R. S. Schechter, M. G. Velarde, and J. K. Platten. The two-component Bénard problem. Wiley, 1981. [25] A. C. Galeão and E. G. D. do Carmo. A consistent approximate upwind Petrov-Galerkin method for convectiondominated problems. Computer Methods in Applied Mechanics and Engineering, 68(1):83 – 95, 1988. [26] J. S. Turner. Double-diffusive phenomena. Ann. Rev. Fluid Mech., 6:37–56, 1974. [27] H. Tabor. Solar Ponds. Solar Energy, 12:549–552, 1969.
2425
ScienceDirect Energy Procedia 57 (2014) 2416 – 2425
2013 ISES Solar World Congress
Numerical Modeling of Solar Ponds M. C. Giestasa, J. P. Milhazesb, H. L. Pinab,* a
LNEG, Estrada do Paço do Lumiar, 1649-038 Lisboa, Portugal IDMEC/IST, Avenida Rovisco Pais, 1049-001 Lisboa, Portugal
b
Abstract A SGSP is a basin of water where solar energy is trapped due to an artificially imposed salinity gradient. In a SGSP three zones can be identified: the surface and bottom zones that are both convective and an intermediate zone in between which is intended to be non-convective. This zone acts as a transparent insulation and allows the storage of solar energy at the bottom where it is available for use. A numerical model where the SGSP dynamics is described in terms of velocity ݑ, pressure , temperature ߠ and salt concentration ߪ is presented. It is based on the Navier-Stokes equations for an incompressible fluid coupled to one advection-diffusion equation for ߠ and one advection-diffusion equation for ߪ . The fluid density ߩ is taken to depend on ߠ and ߪ and the Boussinesq hypothesis is adopted: the fluid density appearing in the LHS of the NavierStokes equation is supposed constant and equal to some reference value whereas it is assumed to be variable in the RHS. The space discretization of the governing equations is based on the respective weak formulations and the discretization employs finite elements with a pressure correction method used to decouple velocity and pressure. Integration in time is accomplished by a BDF (Backward Differentiation Formula) method with the above PDEs treated sequentially within each time step. A computer code was developed employing the finite element class library deal.II. Comparisons with available experimental results are made to validate this numerical model. © 2014 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license © 2013 The Authors. Published by Elsevier Ltd. (http://creativecommons.org/licenses/by-nc-nd/3.0/). Selection and/or peer-review under responsibility of ISES. Selection and/or peer-review under responsibility of ISES Keywords: Solar ponds, numerical modeling, finite elements, double advection-diffusion.
* H. L. Pina. Tel.: +351-965-683-057; Fax: +351-218-419-634. E-mail address: [email protected]
1876-6102 © 2014 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/3.0/). Selection and/or peer-review under responsibility of ISES. doi:10.1016/j.egypro.2014.10.250
M.C. Giestas et al. / Energy Procedia 57 (2014) 2416 – 2425
1. Introduction Nomenclature A B b g h L1, L2 n p σ θ t
Thermal diffusivity [m2/s] Salt diffusivity [m2/s] Buoyancy Gravity acceleration [m/s2] Convection coefficient [W/m2K] Domain length and height Unit normal vector Pressure [Pa] Salt concentration Temperature [K] Time [s]
Subscripts 0 f r a 1, 2, 3, 4
Initial Final Reference values Ambient values Left, right, bottom, top
u u Ω, ∂Ω α β γ ρ . Δ x
Generic scalar field Velocity field [m/s] Domain, boundary Thermal expansion coeff. [K-1] Salt contraction coefficient Kinematic viscosity [m2/s] Fluid specific mass [Kg/m3] Gradient Divergence Laplacian Space coordinates
Superscripts ż n
Time derivative Time instants
A Salt Gradient Solar Pond (SGSP) is a basin of water where solar energy is stored. A temperature gradient (hotter at the bottom and cooler at the top) is established in a natural way due to solar radiation absorption at the surface and the salt concentration gradient (denser at the bottom and lighter at the top) acting to prevent convective motions that would otherwise promote the return of the stored energy to the outside ambient destroying the SGSP very purpose. Thus a double advection-diffusion process occurs with the temperature and salinity fields making opposite contributions to the fluid density. The experimental, analytical and numerical studies of SGSP can be traced back to Kalecksinsky [16] and was pursued by Stommel [8], Weinberger [17], Turner [26] and Tabor [27]. Recently there was a renewed interest on SGSPs as a device coupled to desalination units or as a low cost seasonal thermal energy storage reservoir. The numerical solution of the governing equations that describe a SGSP has been attempted. Hull [3], Hawlader and Brinkworth [15] and Rubin [22] have applied a finite difference method while [4] has used a finite element technique. The SGSP stability – one of the key factors governing a SGSP performance – has been studied by several researchers [9], [6], [24], [11] and [12]. Some of those authors resorted in most cases to a linear perturbation theory and the results obtained have provided important information regarding the onset of the instabilities as well as the existence of several possible stable or unstable states that may arise. Weinberger [17] gave a mathematical formulation of a SGSP behavior. Meyer [19] developed a numerical model to predict the time dependent behavior of the interface between the convecting and the non-convecting zones. Panahi et al. [20] simulated the dynamics of the SGSP with a one-dimensional model. Angeli and Leonardi [1] and [2] studied the evolution of salt concentration profiles and investigated the salt diffusion and stability of the density gradient. Kurt [18] modeled the non-convective zone as a series of horizontal layers with the upper and lower convective zones modeled as single homogeneous layers. Mansour et al. [21] solved numerically the problem of transient heat and mass transfer and long term stability of a SGSP through a two-dimensional finite volume method.
2417
2418
M.C. Giestas et al. / Energy Procedia 57 (2014) 2416 – 2425
A comprehensive numerical model to simulate the SGPS dynamics is still lacking and the present paper is intended as a contribution to the development of such a tool. Thus, a rectangular cavity filled with a mixture of water and salt and heated at the bottom is considered. The fluid is treated as newtonian incompressible, heat conducting according to Fourier’s law and the salt diffusion obeys Fick’s law and are subject to an uniform gravitational field. Also the Boussinesq hypothesis is adopted. Due to the complexity of this time-dependent coupled system of equations it is extremely difficult to devise a numerical method of solution that performs equally well with respect to the several physical phenomena involved. Therefore the problem is separated into a sequence of simpler subproblems in correspondence with the underlying physics thus allowing to apply the most adequate numerical method to each subproblem. The present study calls for the splitting of the temperature, salinity and NavierStokes equations the later employing a pressure correction approach leading itself also to a decoupling of velocity and pressure. Integration in time is accomplished by a BDF(2) (Backward Differentiation Formula, second order) method with the above PDEs treated sequentially within each time step. A computer code was developed using the finite element class library deal.II (see [5]). To validate this numerical model some comparisons with available experimental results are presented. 2. The governing equations The governing equations are the standard conservation laws for continuum media. The water-salt solution is assumed to behave as a newtonian fluid conducting heat according to Fourier’s law and the salt diffusion is governed by Fick’s law. Also the Boussinesq hypothesis is adopted: the specific mass appearing in the left hand sides of the relevant equations is supposed constant and equal to some average value whereas in the right hand sides it is assumed to vary with temperature and salt. We collect here the governing partial differential equations as well as the boundary and initial conditions adopted for the case studied. Heat equation: μ௧ ߠ ሺ ڄ ܝሻߠ ൌ ܣȟߠπ ൈ ሺͲǡ ܶ ሿ ܖ ڄ ߠܣ ݄ሺߠ െ ߠ ሻ ൌ Ͳ μπ ൈ ሺͲǡ ܶ ሿ ߠሺܠǡ Ͳሻ ൌ ߠ π
(1) (2) (3)
μ௧ ߪ ሺ ڄ ܝሻߪ ൌ ܤȟߪπ ൈ ሺͲǡ ܶ ሿ ܖ ڄ ߪܤൌ Ͳ μπ ൈ ሺͲǡ ܶ ሿ ߪሺܠǡ Ͳሻ ൌ ߪ ሺܠሻπ
(4) (5) (6)
Salt equation:
Navier-Stokes system: μ௧ ܝ ሺ ڄ ܝሻ ܝൌ ߥȟ ܝെ ܊π ൈ ሺͲǡ ܶ ሿ ܝ ڄ ൌ Ͳπ ൈ ሺͲǡ ܶ ሿ ܖ ڄ ܝൌ Ͳ μπ ൈ ሺͲǡ ܶ ሿ ܝሺܠǡ Ͳሻ ൌ π
(7) (8) (9) (10)
with the buoyancy right hand side ܊given by (11) ܊ൌ ܴሺߠǡ ߪሻ ؠሺͳ െ ߙሺߠ െ ߠ ሻ ߚሺߪ െ ߪ ሻሻ where ߙ and ߚ are the expansion coefficients for temperature and salt and ߠ and ߪ are some reference values.
2419
M.C. Giestas et al. / Energy Procedia 57 (2014) 2416 – 2425
As can be seen, our model consists of a system of coupled partial differential equations comprising two time dependent advection-diffusion equations and the time dependent Navier-Stokes equation. This amounts to a Rayleigh-Bénard problem with two components (temperature and salt concentration). 3 Discretization The governing equations are solved by a sequential splitting method described in Fig. 1. Therefore the problem is decomposed into a sequence of simpler subproblems in correspondence with the underlying physics thus allowing to apply to each subproblem the most adequate numerical method. On the other hand this procedure is first order meaning that the splitting truncation error is only ࣩሺȟݐሻ (a general introduction to operator splitting techniques can be seen in [14] or [7]). Initializations: ߠ ൌ ߠ Ǣߪ ൌ ߪ Ǣܝ ൌ ܰ ൌ ሾܶ Ȁȟݐሿ FOR ݊ ൌ Ͳ TO ܰ DO Solve temperature equation for ߠ ାଵ Solve salt equation for ߪ ାଵ Solve Navier-Stokes equations for ሺܝାଵ ǡ ାଵ ሻ END Fig. 1: Program structure. 3.1 Time integration Time derivatives are approximated by the BDF(2) (Backward Differentiation Formula, second order) method with constant time step ȟݐ, that is, for ݊ ͳ ଵ ଷ ଵ ߠሶሺݐାଵ ሻ ൎ ቀ ߠ ାଵ െ ʹߠ ߠ ିଵ ቁ ߪሶ ሺݐାଵ ሻ
ൎ
௧ ଶ ଵ ଷ
ቀ ߪ
ାଵ
௧ ଶ ଵ ଷ
ଶ ଵ
െ ʹߪ ߪ ିଵ ቁ ଶ ଵ
(12)
ܝሶ ሺݐାଵ ሻ ൎ ቀ ܝାଵ െ ʹܝ ܝିଵ ቁ ௧ ଶ ଶ The BDF(2) method is known to be second order accurate and stiff A-stable (see [13] for details). As BDFs methods are not self starting a backward Euler method was employed for the first time step. that is, To compute approximations to the advective (non-linear) terms, the following linearly extrapolated values will be used (13) ߠכ ൌ ʹܝ െ ܝିଵ ǡߪכ ൌ ʹߪ െ ߪ ିଵ ǡܝ כൌ ʹܝ െ ܝିଵ The heat and salt equations are thus discretized in time as ଵ ଷ ାଵ ቀ ߠ ௧ ଶ ܖା
ଵ
െ ʹߠ ߠ ିଵ ቁ ሺܝ ڄ כሻߠ ାଵ ൌ ܣȟߠ ାଵ
(14)
ߠܣ ܖ ڄ ݄ሺߠ ାଵ െ ߠ ሻ ൌ Ͳ ଵ ଷ ାଵ ଵ ቀ ߪ െ ʹߪ ߪ ିଵ ቁ ሺܝ ڄ כሻߪ ାଵ ൌ ܤȟߪ ାଵ
(15) (16)
ܖ ߪܤା ܖ ڄൌ Ͳ
(17)
௧ ଶ
ଶ
ଶ
The Navier-Stokes system is discretized in time using a sequential operator splitting method to decouple velocity and pressure thus making this later variable explicit (see [10] for a detailed description). In the first substep, an intermediate (non divergence–free) velocity ܝ ෝ ାଵ is obtained by solving
2420
M.C. Giestas et al. / Energy Procedia 57 (2014) 2416 – 2425 ଵ
ଷ
ଵ
ቀ ܝ ෝ ାଵ െ ʹܝ ܝିଵ ቁ ሺܝ ڄ כሻܝ ෝ ାଵ ൌ ߥȟܝ ෝ ାଵ െ ͓ ܊כ ଶ ܖڄൌͲ ܝ ෝ Here ͓is an extrapolated pressure given by ସ ଵ ͓ൌ ߶ െ ߶ ିଵ ଷ ଷ and the buoyancy is computed explicitly by ܊ כൌ ܴሺߠכ ǡ ߪכ ሻ ௧ ଶ ାଵ
(18) (19) (20) (21)
In the second substep we correct to obtain a divergence–free velocity ܝାଵ and a pressure ାଵ by solving ଷ ሺܝାଵ െ ܝ ෝ ାଵ ሻ ൌ ሺ ͓ െ ାଵ ሻ (22) ଶ௧ ାଵ ൌͲ (23) ܝڄ (24) ܝାଵ ܖ ڄൌ Ͳ Applying the divergence operator to equation (22) we get a Poisson equation for the pressure ାଵ with homogeneous Neumann boundary conditions: ଷ ܝڄ ෝ ܖା (25) ȟሺାଵ െ ͓ሻ ൌ ଶ௧ ାଵ ൌͲ (26) μ ܖ The new pressure ାଵ is updated by ାଵ ൌ ߶ ାଵ and the new velocity ܝାଵ by ଶ௧ ሺ ͓െ ାଵ ሻ ܝାଵ ൌ ܝାଵ כ ଷ
(27) (28)
3.2 Space discretization For the space discretization, a Galerkin variational formulation is employed. The following finite element approximations will be used in the examples: (i) Temperature ߠ and salt ߪ: biquadratic elements (ܳଶ ); (ii) Velocity ܝand pressure : bilinear elements (ܳଵ ). 3.3 Stabilization For SGPS problems, equations (1) and (4) happen to be advection dominated. It is now well known that in this case Galerkin type methods lack stability and can return meaningless solutions usually displaying non physical oscillations or wiggles. It is therefore necessary to provide some form of stabilization in order to recover from these difficulties which is generally achieved by suitably modifying the underlying weak form. The method adopted was the isotropic artificial diffusion SOLD stabilization as described in [25]. 3.4 The linear algebraic systems The linear algebraic systems resulting from the space and time discretizations were solved by the conjugate gradient method with ILU(0) preconditioning (see for instance [23]).
M.C. Giestas et al. / Energy Procedia 57 (2014) 2416 – 2425
4. Results 4.1 Constants The constants values (in SI units) used to obtain the results were: Δ ݐ ൈ ͳͲିଶ ͷ ൈ ͳͲିଶ ିଵ ݄ଵ ͳͲǤͲ ͳ ൈ ͳͲ ݄ଶ ݄ଵ ͳǤͲͳͶ ൈ ͳͲଷ ݄ଷ ͵ʹͲǤͲ ͵Ǥ͵ͷ ൈ ͳͲିସ ݄ସ ͳͲͲǤͲ Ǥͷ ൈ ͳͲିସ ʹʹǤͻͷ ߠଵ ʹʹǤͻͷ ͲǤͲ ߠଶ ߠଵ ߠଷ ͷͲǤͲ ͲǤͺͲͷ ൈ ͳͲି ߠସ ʹʹǤͻͷ ͳǤ ൈ ͳͲି ߠ ʹʹǤͻͷ ͳǤͶ ൈ ͳͲିଽ ͶǤͻͳ ൈ ͳͲଷ ݔଶ ݔଶ ͲǤͺͷ ͲǤͲͷͳ െ ͲǤʹ ൬ͳ െ ൰ ǡ Ͳ ߪ ܮଶ ܮଶ ൞ ݔଶ ͲǤͲͲͳʹǡͲǤͺͷ ͳ ܮଶ
ܮଵ ܮଶ ߩҧ ߙ ߚ ߠ ߪ ߥ ܣ ܤ ܶ
4.2 Geometry and mesh The domain for the examples in Section 4 is a rectangle π ൌ ሾͲǡ ܮଵ ሿ ൈ ሾͲǡ ܮଶ ሿ as depicted in Fig. 2. Its boundary is μπ ൌ ڂସୀଵ Ȟ , where the Ȟ are the rectangle sides.
Fig. 2: Geometry and notation. The mesh comprises Ͳ ൈ ͳͶͲ ൌ ͻͺͲͲ elements and ͳͲͲͳͳ degrees of freedom for each component of the velocity and for the pressure, ͵ͻʹͳ degrees of freedom for the temperature and ͵ͻʹͳ degrees of freedom for the salt concentration. 4.3 Model validation A comparison of experimental temperature profiles and shadowgraph images is made against the numerical results. In the numerical case (NC) the vertical component of the density gradient was calculated in order to create a Schlieren image for the whole domain. This in turn made possible the comparison between the shadowgraph image and the Schlieren type. Temperature profiles were calculated at ൎ ͷ ൈ ͳͲିଷ from the right wall (see shadowgraph images). The ruler seen on the right side of the shadowgraph image gives an indication on the distance from the bottom in centimeters.
2421
2422
M.C. Giestas et al. / Energy Procedia 57 (2014) 2416 – 2425
Fig. 3-5 depict the shadowgraph images for the experimental case (EC) together with the temperature profile for a given time in seconds. Bellow we show the results of temperature profile and Schlieren obtained by the numerical model for comparison. Fig. 3 corresponds to an early state of the experiment. At time ݐൌ ͷͶͲs, a small but unstable lower convective zone (LCZ) is being developed. The upper convective zone (UCZ) in both NC and EC can also be observed and a second convective zone (CZ2) if beginning to form. At ݐൌ ͳͺͲ (see Fig. 4), one can notice that the separation between the LCZ and the gradient zones (GZ) is now sharply defined at ൎ ͵ ൈ ͳͲିଶ from the bottom. There is still a CZ2 but very unstable. In Fig. 5 we can verify that both in the NC and EC the CZ2 has vanished. The interface between the LCZ and GZ has is now at Ͷ ൈ ͳͲିଶ . Also, small movements can be observed on the UCZ due to heat loss from the surface. In all the cases shown above there is a reasonable approximation between the EC and NC. Fig. 6 shows a side-by-side comparison of the velocity magnitude vertical profile at a given time. In the EC the velocity was obtained using a PIV (Particle Image Velocimetry) technique. In both cases the vertical profile was taken along a vertical line at ൎ ʹ ൈ ͳͲିଶ m apart from the left wall.
(a)
(b) Fig. 3: (a) Temperature profile and shadowgraph image for EC. (b) Temperature profile and Schlieren image for NC. Time t=540 s
M.C. Giestas et al. / Energy Procedia 57 (2014) 2416 – 2425
(a)
(b) Fig. 4: (a) Temperature profile and shadowgraph image for EC. (b) Temperature profile and Schlieren image for NC. Time t=1680 s
2423
2424
M.C. Giestas et al. / Energy Procedia 57 (2014) 2416 – 2425
(a)
(b) Fig. 5: (a) Temperature profile and shadowgraph image for EC; (b) Temperature profile and Schlieren image for NC. Time t=4691 s
(a) (b) (c) (d) Fig. 6: Velocity magnitude profile for (a) EC; (b) NC at t=1860 s; (c) for EC; (d) NC for t=3000 s 5. Conclusions The results presented – both experimental and numerical – comprise a subset which are considered representative of the results available. Despite the small scale experimental setup it is believed that the information gathered can be useful for the full scale SGPS dynamics.
M.C. Giestas et al. / Energy Procedia 57 (2014) 2416 – 2425
The numerical model developed seems to be able to predict the main features of the physics involved. However it is still two-dimensional and its capacity to simulate the long term dynamics has yet to be demonstrated. The move to a three-dimensional model requires a particular attention to the computer implementation, namely code parallelization issues, and the employment of a robust time-integrator. Acknowledgments The work reported was supported by IDMEC/IST. J. P. Milhazes was supported by a FCT studentship grant. The computations described in this article were made using the finite element library deal.II (see [5]). References [1] C. Angeli and E. Leonardi. A one-dimensional numerical study of the salt diffusion in a salinity-gradient solar pond. Int. J. Heat Mass Transfer, (47):1–10, 2004. [2] C. Angeli and E. Lombardi. The effect of thermodiffusion on the stability of a salinity gradient solar pond. Int. J. Heat Mass Transfer, 48(21-22), 2005. [3] J. R. Hull. Computer simulation of solar pond thermal behaviour. Solar Energy, 25(33):33–40, 1980. [4] T. S. Jayadev and J. Henderson. Salt concentration gradient solar ponds. Modeling and Optimization. Colorado, 1979. [5] W. Bangerth, R. Hartmann, and G. Kanschat. deal.II — a general-purpose object-oriented finite element library. ACM Trans. Math. Softw., 33(4), 2007. [6] L. N. Da Costa, E. Knobloch, and N. O. Weiss. Oscillations in double-diffusive convection. Journ. Fluid Mech., 109(257):25–43, 1981. [7] I. Farago and A. Havasiy. Operator Splittings and their Applications. Nova Science Publishers, 2009. [8] H. Stommel, B. Arons, and D. Blanchard. An oceanographical curiosity: the perpetual salt fountain. Deep Sea Res., 3:152–153, 1956. [9] G. Veronis. Effect of a stabilizing gradient of solute on thermal convection. Journ. Fluid Mech., 34:315–336, 1968. [10] J. L. Guermond, P. Minev, and Jie Shen. An overview of projection methods for incompressible flow. Comput. Methods Appl. Mech. Engrg., (195):6011–6045, 2006. [11] M. Giestas, A. Joyce, and H. Pina. The influence of radiation absorption on solar ponds stability. J. Heat Mass Transfer, 39,18:3873–3885, 1996. [12] M. Giestas, A. Joyce, and H. Pina. The influence of non-constant diffusivities on solar ponds stability. J. Heat Mass Transfer, 40, 18:4379–4391, 1997. [13] E. Hairer and G. Wanner. Solving Ordinary Differential Equations II – Stiff and Differential-Algebraic Problems. Springer-Verlag, 1991. [14] W. Hundsdorfer and J. G. Verwer. Numerical Solution of Time-Dependent Advection-Diffusion-Reaction Equations. Springer-Verlag, 2003. [15] M. N. A. Hawlader and B. J. Brinkworth. An analysis of the non convecting solar pond. Solar Energy, 27:195–204, 1981. [16] A. V. Kalecsinsky. Ueber die ungarishen warmen und heissen kochsalzseen als natuerlich waermeaccumulatoren. Ann. Physik IV, 7:408. [17] H. Weinberger. The Physics of the solar pond. Solar Energy, 8(2). [18] H. F. Kurt, Halicir, and A. Binark. Solar pond conception - experimental and theoretical studies. Energy Conv. Manag., 41(9), 2000. [19] K. A. Meyer, D. P. Grimmer, and G. F. Jones. An experimental and theoretical study of salt-gradient pond interface behaviour. Progress in Solar Energy, 1982. [20] Z. Panahi, J. C. Batty, and J. P. Riley. Numerical simulation of the performance of a salt gradient solar pond. Transactions of ASME, J. of Solar Energy Enginnering, 105:361–374, 1983. [21] R. Mansour, C. Nguyen, and N. Galanis. Transient heat and mass transfer and long-term stability of a salt-gradient solar pond. Mech. Research Communications, 33:233–249, 2006. [22] H. Rubin, B. A. Benedict, and S. Bachu. Modeling the performance of a solar pond as a source of thermal energy. Solar Energy, 32:771–778, 1984. [23] Y. Saad. Iterative Methods for Sparse Linear Systems. SIAM, 2003. [24] R. S. Schechter, M. G. Velarde, and J. K. Platten. The two-component Bénard problem. Wiley, 1981. [25] A. C. Galeão and E. G. D. do Carmo. A consistent approximate upwind Petrov-Galerkin method for convectiondominated problems. Computer Methods in Applied Mechanics and Engineering, 68(1):83 – 95, 1988. [26] J. S. Turner. Double-diffusive phenomena. Ann. Rev. Fluid Mech., 6:37–56, 1974. [27] H. Tabor. Solar Ponds. Solar Energy, 12:549–552, 1969.
2425