Molten Silicon Storage of Concentrated Solar Power with Integrated Thermophotovoltaic Energy Conversion

. A new kind of systems combining latent heat energy storage in molten silicon and thermophotovoltaic (TPV) heat-to-power conversion are under development within the AMADEUS (www.amadeus-project.eu) project. The extremely high latent heat of silicon (1230 kWh/m 3 ) plus the very high electrical power density of TPV (several 10’s of kW/m 2 ) will eventually enable the fabrication of ultra-compact CSP systems that integrate thermal energy storage and power generation in the same unit. This work deals with the search of the optimal geometry of the PCM vessel, enabling the highest energy transfer to the TPV converter. Two kinds of geometries are explored: the inverted truncated pyramid (ITP) and the hollow cylinder (HC). A simplified quasi-1D semi-analytical model for the heat transfer in the PCM coupled to a TPV optical cavity model is used to simulate the system and determine some key figures of performance such as discharge efficiency, discharge time or electrical power. On top of that, the assumptions made in the 1-D model are verified against a well advanced 3-D Computational Fluid Dynamics (CFD) model, which takes into account buoyancy effects, dendrites formation and the PCM expansion during solidification, being thus in position to describe in more detail the induced complicated fluid structures governing both the phenomena of PCM melting and solidification. We find that both the analytical model and the advanced 3D CFD model coincide for most of the simulation time, enhancing the validity of the analytical model that can be implemented during the initial stages of the system design. Finally, the article discusses on the best geometries resulting in the maximum system efficiency and output power, as a result from the 1D analysis.

INTRODUCTION Solar thermal energy storage based on very high melting point PCMs, such as pure silicon and boron (melting points of 1414 ºC and 2076 ºC), plus thermophotovoltaic (TPV) energy conversion has been proposed theoretically in the past [1][2][3][4][5], the main motivation being the extremely high latent heat and thermal conductivity of these PCMs if compared with salt-based materials. In particular, silicon (1230 kWh/m 3 ) and boron (2680 kWh/m 3 ) have latent heats an order of magnitude greater than that of typical salts used in CSP such as NaNO3 (110 kWh/m 3 ) and KNO3 (156 kWh/m 3 ). Actually, silicon and/or boron PCMs have higher storage energy densities than most forms of energy storage, including electrochemical batteries (~ 250-500 Wh/m 3 ), or pressurized hydrogen (~ 1500 kWh/m 3 at 700 bars), and it is only clearly surpassed by fuels such as liquefied natural gas (~ 6000 kWh/m 3 ) or gasoline (~ 9500 kWh/m 3 ). The obvious technological challenge is the very high operation temperature, especially concerning the PCM vessel design and the heat-to-power conversion system. Maximum operation temperatures of conventional dynamic closed-cycle engines such as closed-Brayton, Stirling and Rankine, are typically well below 1000 ºC [6]. This is mostly because of the very serious concerns with the working fluid stability and structural mechanical strength of the engine steel parts at high temperatures. On the contrary, thermophotovoltaics (TPV) [7,8] is perfectly suited for such high operation temperatures. TPV directly convert radiative heat from incandescent bodies into electricity; thus, eliminating the need of working fluid and moving parts. Furthermore, TPV cells are planar devices that provide very high power densities of several 10's of kW/m 2 . This enables the development of very compact devices integrating storage and conversion in the very same unit. This modularity will eventually enable the development of decentralized CSP systems, or even new devices for energy storage in the housing and district sectors.
Among the several technological challenges of this technology, the geometrical configuration of the PCM container is particularly relevant, in order to maximize the efficiency and power output of the system. This geometry determines the specific location of the TPV converter, which will subsequently determine the heat flow through the PCM. In this article, we present a theoretical analysis of two possible geometrical configurations that are under consideration in the AMADEUS project ( Figure 1): the inverted truncated pyramid (ITP), and the hollow cylinder (HC). In both geometries, all the container walls are assumed adiabatic except for that one comprising the TPV emitter: in the ITP this is the smaller square surface at the bottom, and in the HC geometry this is the inner cylindrical surface ( Figure 1). The TPV cells face the emitter and transform its radiant heat directly into electricity. Consequently, heat flows in the direction of the TPV emitter during the solidification. At the beginning of the solidification process, the emitter is assumed to have the temperature of silicon's melting point (1680 K) and a solid-liquid interface is created very near this surface. During the solidification process, latent heat is released from that interface and flows towards the emitter. Thus, the solid-liquid interface moves away from the emitter, creating a solid silicon crust that will negatively impact on the heat transfer from the solid-liquid interface to the TPV emitter.
Selecting a proper PCM vessel geometry is essential to enhance the heat transfer and therefore, achieve high system efficiency. In this regard, notice that the two geometries considered in this work ( Figure 1) are characterized by having a cross-sectional area (perpendicular to the heat flux) that diminishes in the direction of heat flux. This makes heat flux density (in W/m 2 ) to increase towards the TPV emitter surface; thus, maximizing the power density at the TPV emitter surface, and subsequently maximizing the emitter temperature during the full solidification process. The purpose of the present work is to addresses the question of which of these two geometries (ITP or HC) leads to the higher efficiencies and power outputs during the solidification (discharge) of the system.

Analytical Model
The modeling of these systems requires the solution of the Stefan problem during the solidification of the PCM. In this article we use a quasi-1D semi-analytical model developed in previous works [3,5]. This model assumes a quasi-stationary approach to simulate the transient response of the system. This assumption is reasonably accurate as the movement of the interface is expected to be slow due to the high latent heat of fusion of silicon. Details on this model for the ITP and HC geometries can be found in [5] and [3] respectively. These models are coupled to a TPV optical cavity model describing the complex radiative exchange taking place between the TPV emitter and the TPV cells. This model has been thoroughly described in previous works [6]. According to this model, the TPV cell is described by its semiconductor bandgap energy ( ) and the internal photoluminescence efficiency ( ), which accounts for the non-radiative recombination losses in the cell. The TPV cells comprise a Back Surface Reflector  (BSR) of reflectivity on their rear surface to reflect back to the emitter the infrared radiation not absorbed in the semiconductor; thus maximizing the conversion efficiency. Finally, the emitter is described by its spectral emissivity emit . Except otherwise indicated, this work assumes the following parameters for the TPV optical cavity: emit = 0.9 (e.g. silicon carbide), =0.7 eV (e.g. Ge, InGaAs and GaSb semiconductors), = 90% (typical of III-V semiconductors) and = 0.85.

CFD Model
Apart from the 1-D model, an advanced 3-D CFD model -designated in this work as Fluent2Phase model-is developed on the ANSYS Fluent platform (v17.1) [9] for the simulation of the silicon solidification inside the ITP geometry. This model, based on the enthalpy-porosity method [10], can describe the unsteady heat transfer mechanisms occurring inside the PCM. Additionally, a simplified version of this model, named as Fluent1Phase model, with settings retrieved from the work of Veeraragavan et al. [5] is developed, as well. Results of the Fluent1Phase model are compared against the ones obtained from the respective 1Phase model developed by [5] in OpenFoam (meltFoam) solver [11] -called as OpenFoam model in the present work-in order to enhance the solidification/melting CFD model validity at different solvers. Finally, the 1-D model is compared against all 3-D models, i.e. Fluent2Phase, Fluent1Phase and OpenFoam model, in terms of temperature profiles and the dimensionless solid-liquid PCM interface at different time instants, in order to prove its accuracy.
Even though a substantial amount of published works in the open literature is dedicated mainly on the melting process of PCMs enclosed inside rigid casings of different shapes, mostly cylinder or sphere [12,13], limited attention has been paid on the PCM solidification. The literature is even more limited as concerns the study of phenomena, such as the PCM volumetric change during its solidification and floating or sinking of the solid material inside the liquid PCM due to their density difference. Veeraragavan et al. [5] in his work has studied the silicon solidification inside the closed truncated-cone geometry neglecting though the density difference. Assis et al. [14] in his work takes into account such an effect, by studying the solidification of paraffin wax at low temperatures inside spherical shell.
In the specific work, the developed Fluent2Phase model, follows the solidification process of molten silicon from beginning to end, occurring at high temperatures inside the ITP geometrical configuration. This model, compared to current state-of-art incorporates the effect of phenomena, as those of natural convection in the liquid phase and volumetric expansion due to solidification. The multiphase volume of fluid (VOF) approach is accompanied by the inclusion of compressibility effects of the inert gas, enclosed inside the PCM casing. The latter can be either compressed inside the closed rigid shell, or exit the container through a release valve, during the silicon expansion during its freezing, in order to avoid the exertion of high stresses in the PCM casing. In this work, the case of the closed shell, without a release valve, is being studied. For a more accurate representation of the gas-PCM and solidliquid interphases, in-house user-defined functions (UDFs) are implemented so that a grid adaptive local refinement [15], both in the PCM-gas and solid-liquid interphases, is applied to enhance the numerical results accuracy without substantial increase of the associated computational cost.
In the enthalpy-porosity model an artificial region, called as mushy zone, is assigned to encounter the phase change through a finite temperature range. Such an assumption is necessary to achieve a smooth transition from solid to liquid PCM properties around the melting point. In the particular model the temperature transition is equal to ΔΤ=2 K, i.e. Tliquidus=Tsolidus+2 K. Within this region the velocity varies from zero (solid PCM) to the liquid velocity (imposed mainly by the liquid natural-convection). The mushy zone parameter, Amush, used in the model [10] as a sink term to take into account the damping effect within the mushy region is dependent on the dendrites formed during solidification process, the material density and the liquid viscosity. In the Fluent2Phase model, this parameter is set to a more realistic value, i.e. 3·10 9 , rather than the default value of Fluent and OpenFoam solvers, i.e.10 5 . Calculations are based on the equation Amush =180 v/DAS 2 , where a small value of dendrites arm spacing (DAS), around 0.1 μm, is considered. The DAS is assumed that small, because under undercooling conditions, the dendrites formation is almost negligible [16]. It should be stressed out that with such a high value of the Amush the solid phase velocities are set equal to zero artificially and thus, the floating of the solid in the liquid silicon cannot be represented, because the solid PCM does not move. In a future model formulation the natural phenomenon of floating or sinking of the solid material inside the liquid PCM due to buoyancy forces, should be also taken into account, for a more realistic representation. This simplified Fluent1Phase model follows several assumptions. First of all, the molten silicon is supposed to behave as incompressible and the density of the two PCM phases -solid and liquid-are considered equal; thus, the PCM volumetric change and the floating of the solid phase inside the liquid silicon are not taken into account. The PCM density is solved by using the Bousinesq approximation, by simply adding a buoyancy source term in the moment equation. By this, the model takes into account only the natural convection within the liquid PCM. Additionally, the vessel is filled only by the PCM, since the 1Phase model can only model one medium. Such a design is not realistic due to the high stresses that will be exerted onto the casing during the PCM expansion. Finally, the mushy zone parameter, Amush, used in solidification/melting model is set at a default value of Fluent and OpenFoam equal to Amush=10 5 . In this simplified CFD model this value does not have a significant effect on the results, for two main reasons. Firstly, the mushy zone is not considered, as the solid temperature is the same with that of the liquid phase. Secondly, any actual density difference between the solid and liquid phases is not taken into account. Thus, the solid phase does not have the tendency to move upwards due to buoyancy forces and its velocities are almost equal to zero regardless the mushy parameter. However, in the advanced model this value should be calculated cautiously as it affects significantly the PCM melting/solidification rate, as highlighted in [17].
The PCM-gas thermophysical properties are listed in Table 1. It is worth noticing that in the 1Phase model the silicon thermal conductivity has a sharp transition around its melting point, whilst in the Fluent 2Phase model there is a small transition owing to the presence of the mushy zone. Additionally, in the Fluent 2Phase model the air compressibility is taken into account by solving Tait equation, while the PCM viscosity in the cells, where the liquid fraction is equal to zero, is set to a high value, i.e. in the order of 10, so that the solver can virtually understand that in the particular cell the material is behaved as being solid. This is a similar approach followed in the work of [18].
The upper and side walls of the tank are considered adiabatic, Figure 2. The ITP geometry is placed vertically, in order to take into account the liquid natural convection. The emitter during the night operation radiates heat, which is dependent on its temperature. For this reason, a custom-UDF is applied to calculate iteratively for each time step the temperature dependent heat flux at the emitter surface. Initially, the temperature is set equal to 1680 K at the emitter surface and equal to 1960 K at the upper part of the tank. For the rest of the domain a linear interpolation of these two values is used to patch the temperature profile. Transient calculations are performed, with a fixed time step set equal to Δt=0.05 s in the Fluent 1Phase model and a variable time step (Courant number=0.2, maximum Δt=0.01 s) in the Fluent 2Phase model. In the Fluent 1Phase model a uniform 3D coarse grid of 4356 hexahedral cells is used, as in the OpenFoam model. In the Fluent 2Phase model, the problem is solved as axisymmetric and the adaptive local refinement technique is implemented, which results in a maximum grid size of almost 7000 cells.  Table 2 shows the four geometrical cases that are analyzed in this article. These cases have been constructed for the ITP geometry and they are intended to explore the dependence of the two main parameters of this configuration [2,5]: the length and the tapering ratio = / = ( / ) (Figure 1). Under the assumptions of this work (adiabatic walls) these two parameters fully determine the dynamics of this system. This means that for a given and , all the values of and satisfying / = will lead to identical results. Thus, (or ) simply become a scaling parameter leading to different volume capacities with identical system performance (e.g. efficiency, discharge time, etc.). Only extensive variables (e.g. total stored energy, output power, etc.) will depend on (or ). Table 2 also shows the four corresponding cases for the HC systems. These cases have been defined to have the same length ( * = ), emitter area ( * = ) and total volume ( * = ) than that of the corresponding ITP cases. This enables a fare comparison between the two geometries, both having the same energy storage (PCM volume) and output power (TPV cell area) capacities. Figure 3 shows the average discharge efficiency and discharge time as a function of the total PCM volume for the eight geometries described in Table 2. Solid-thick lines represent the four ITP cases, while the four thinner lines represent the four cases of HC geometries. Discharge time is defined as the elapsed time from all-liquid to all-solid situations. Average discharge efficiency is defined as the ratio of the total delivered electrical energy to the total released heat. Figure 4 shows the maximum and minimum power during the discharge of the system. Maximum power corresponds to the initial situation in which the emitter is at silicon's melting point of 1680 K. This is why both ITP   Table 2. Thick solid lines represent the truncated-pyramid geometry.
and HC system provide the same initial maximum power (Figure 4a). Thus, the difference between the four cases is only attributed to the different TPV cell areas. Finally, Figure 5 show the main geometrical dimensions of the ITP and HC systems.
For the ITP geometry we clearly observe the benefit of reducing the tapering ratio = / (case 2 respect to case 1, or case 4 respect to case 3). Smaller TR boost the efficiency (Figures 3a) and the discharge time (Figure 3b) at the expense of reducing the output power (Figures 4) due to the smaller emitter area = (Figure 5a). On the other hand, increasing the length L drastically deteriorates the ITP system efficiency (Figures 3a), due to the thicker solid crust that difficult the heat flux from the liquid PCM to the TPV emitter. In summary, building a high efficient ITP system with large storage capacity requires a short length L, large , and small . In the limit, this leads to very high aspect ratios.
In the HC configuration, we can increase the volume without scarifying the system efficiency by increasing the length L of the system. However, for a given length, we observe that the highest efficiencies are obtained for small volumes, corresponding to small and − (case 4). This situation corresponds to the lowest "tapering ratio" and the thinnest solid crust. In brief, building a high efficient ITP system with large storage capacity requires a large length (L), small and small − . In the limit, this lead to very low aspect ratios.  . Maximum power output at the beginning of the discharge (a) and minimum output power at the end of discharge (b) for the different system geometries described in Table 2. Thick solid lines represent the ITP geometry.  Table 2.
From the discussion above, we may conclude that the selection of the system geometry will probably depend on the desired final aspect ratio of the system and the requirements concerning average power output and discharge time. ITP (HC) geometries will find the highest efficiency for large (low) aspect ratios. Intermediate cases must be assessed in detail taking into account practical limitations such as the minimum size of the TPV cell, the maximum length of the container taking into account the allowable pressure on its walls, etc. Very importantly, heat losses through the container walls must be assessed in order to determine the validity of these idealistic results. For instance, very small tapering ratios resulting in very low output power and long discharge times will be more sensitive to heat losses through the container walls than those cases in which discharge time is very short. All these aspects will be evaluated in the future using more the complete CFD simulations tools.

Analytical Model Validation against Advanced CFD Modelling
In order to facilitate the comparison with previous works, the different CFD models have been applied to the particular case of ITP geometry described in [5], in which TR = 0.45 and the TPV optical cavity is described by the following parameters: emit = 1 (black body), = 0.51 eV, = 100% and = 0.9. From the CFD analysis it was revealed that for most of the simulation time results of the advanced Fluent 2Phase model virtually coincide with that of the Fluent 1Phase, OpenFoam and Analytical 1-D model (Figure 6). At approximately 2 hours of simulation time, the temperature results indicate that the solidification is terminated and afterwards the temperature decreases rapidly. Such a good agreement enhances the applied models validity. More specifically, it is proven that during the silicon solidification inside the closed vertical shell, the dominant heat transfer mechanism is thermal conduction. Thus, neglecting in the 1D model of phenomena, such as the liquid natural convection and PCM expansion, is expected to play a minor effect on the temperature profile prediction inside the casing for most of the solidification process. A small deviation, of almost 4 %, between the 1Phase and 2Phase model is tracked as concerns the prediction of the overall discharge time. Such discrepancy noticed is mainly attributed to the fact that the solidification process in the 2Phase model lasts more than in the 1Phase one, due to PCM expansion. Actually, the PCM solidification lasts for 2 hours based on the 2Phase model results, whilst in the 1Phase model lasts around to 1.9 hours. An additional small effect on the solidification process plays the high Amush value applied, which assumes that the dendrites formation during silicon solidification is negligible. This approximation is in accordance with the silicon behavior for such low undercooling conditions, as observed as well by [16]. Even so, the silicon  It is worth noticing that with the 1D and 1Phase CFD models any PCM volume changes and dendrites formation effect on the phenomenon temporal evolution are not taken into account; the PCM volume change during solidification is approximately 7%. Contrary to that the Fluent 2Phase model predicts the PCM solidification process in a more realistic manner, since volume expansion due to density variation is considered. Nevertheless, for the specific PCM material, this phenomenon is of minor importance; though it should be included especially when density variations between solid and liquid phases are high enough (>15%). In the specific problem modelled, the silicon is contained in a rigid casing in a shape of truncated cone, where the direction of the smallest resistance to expansion is the axial. During the PCM solidification a slightly higher expansion can be observed near the PCM centerline than the one noticed near the walls (Figure 7). This is reasonable since the distance between the lateral walls of the vessel increase with height and therefore in those areas both axial and radial expansion is observed. Such an observation cannot be predicted when the 1Phase model is applied, where the PCM volume change is not taken into account.
Concluding, many aspects of the real phenomenon are neglected in the 1D and 1Phase CFD models. However, the analytical 1-D model can be considered a valuable tool for the PCM casing design, especially during the earlier design stages, as it can quantify very fast important design parameters, without high error. On the other hand, the advanced CFD 2Phase model should be applied for more complicated and with more phases inter-reacting applications, characterized by high gradients of different properties, and for the more realistic representation of the solidification process in the case of more complex and promising designs. In addition, influential operating parameters that cannot be easily predicted from 1D models, such as pressure distribution, thermal stresses, thermal losses etc., can be described in detail by more sophisticated than 1-D models.

SUMMARY AND CONCLUSIONS
This article presents the simulation of a new kind of energy storage devices that are under development in the AMADEUS project, which comprises pure silicon PCM and thermophotovoltaic converters. The systems have been simulated using a quasi-1D semi-analytical model that enabled the exploration of a set of different geometrical configurations belonging to two main families: the inverted truncated pyramid (ITP) and the hollow cylinder (HC). In the ITP (HC) system the highest efficiency requires high (low) aspect ratios, i.e. low and wide (high and narrow) designs. Thus, the final geometrical choice will depend on concerns such as preferred aspect ratio, maximum allowable length of the container, minimum size of the TPV cell, etc. An advanced 3D CFD solidification model has been also developed to conduct a more profound analysis of the systems. This model considers the effects of a) PCM expansion during solidification, and b) liquid natural convection and dendrites formation during such type of processes.
Simulation results indicate that both 1D-analytical and CFD models agree very well for most of the simulation time.
In future simulations, sophisticated CFD models should be applied for the needs of a more thorough investigation concerning the prediction of PCM solidification/melting process, such as the analysis of heat losses through the container walls.