Numerical study of heat and mass flow in layered heat storage

An analysis of argon gas flow in exemplary configuration of layered bed thermal energy storage is presented. The analysis incorporates URANS model with conjugate heat transfer between gas and solid storage core. The aim of this analysis was to identify key factors limiting exergy efficiency of this storage type and determine some details about storage transient behaviour. Three full cycles of storage loading and unloading having 17 hours physical time in total are simulated, with calculation of exergetic efficiency for each of the cycles. Conclusions regarding this storage type feasibility for indirect power storage in pumped heat systems are made.


Introduction and motivation
This work refers to article provided by authors [1] and awaiting publication.
Many energy markets experience rapid increase in renewable energy sources share, motivated by need to reduce CO2 emissions and related climate changes [2]. This solution suffers from general mismatch between renewable (mainly solar and wind) production and energy demand of the grid. Overproduction should be stored, otherwise it is wasted. Stored surplus can be then used when demand rises above production. Modern EES (electrical energy storage) methods are seen as solution to face this problem.
Considered EES methods must meet many economical and operational often conflicting requirements such as: meeting peak energy demand, flattening RE generation fluctuations, support of energy system management [3] and provide flexibility for introduction of new forms of modern economy like smart grids, electrical vehicle grids and distributed power generation.
Economic viability determined by daily electricity price ratio is fundamental. There are many important parameters of ES methods characterizing their utility such as overall efficiency, energy density, power density, specific energy, specific power, power rating and rated energy capacity. For commercial feasibility, system electric efficiency must be at least as high, as the average ratio of low energy price (when surplus of energy production occurs) to the high price (when there is a shortage of energy).
EES based on indirect storage of energy in the form of heat (Thermal Energy Storage -TES), offers economically favorable investment costs with relatively high energy density, especially comparing to batteries, considered as potential best solution for modern ES approaches [4]. Overall efficiency of TES systems is impacted most by compression/expansion inefficiencies and thermal store exergy losses [5]. Resulting efficiencies are usually too low, when referring to daily electricity price fluctuations in most energy markets. Therefore a concept called pumped heat energy storage (PHES) is believed to meet this obstacle reaching more than 70% overall efficiency. PHES concept is based on heat pump thermodynamic cycle principle during loading and heat engine cycle during discharging of storage tank. The diagram of one of the possible PHES configurations is shown in figure 1. G.F Frate et al. proved that PHES is promising ES technology and founded that the main approach to boost PHES performance is thermal integration. They achieved best round-trip efficiency results for working fluid R1233zd at 110°C, but for large industry scale applications it is crucial to prove economical reasonable gas [6].
Roskosch and Atakan investigated PHES with Rankine cycle (ORC cycle). They conclude that overall efficiency decreases deeply with increasing storage temperature and achieved 70% of overall efficiency in most favorable case [7].
Most of performed work on PHES take advantage of expensive working fluid like ORC, ammonia or refrigerants. It involves high costs of their purchase and further recycling or disposal. It impacts highly operating costs and final profitability. Therefore it is necessary to investigate less sophisticated working fluid to face this obstacle. Because of its thermodynamic characteristics, availability and safety reasons, argon seems to be very promising working fluid for utilization in PHES systems.
Small number of papers describe existing prototypes of PHES in the context of heat transfer and storage parts of the system [8]. This paper concerns the CFD calculation of fluid and heat flow in fixed bed heat storage reservoir with argon as working fluid.

Heat storage model
The numerical model used in the analysis represents a slice of a 32 layer storage, with layers of solid material with holes arranged in hexagonal pattern. Every second layer has the hole pattern shifted in a way that n th layer hole axes are positioned exactly in the middle of next layer's holes. Such arrangement makes the jet created in each hole to impinge on the surface of next layer, intensifying the heat transfer between argon gas and heat storing solid. A scheme of the arrangement, with smallest possible part of the arrangement that can be extracted from it with pattern symmetries is shown on figure 2. Such cut of the 32 layer storage was investigated numerically with a transient conjugate heat transfer analysis.

Physical model
Parameters of the storage core solid, are generic values of 500 J*kg -1 K -1 specific heat capacity, 8000 kg*m -3 density, and 50 Wm -1 K -1 thermal conductivity. This parameters are close to ingot iron parameters, which is one of the considered core material candidates because of its low cost and relatively high thermal conductivity compared to ceramic materials. The gas used in the simulation was argon, modelled as incompressible ideal gas model where its density is modelled as: Where: , , -are density, individual gas constant and temperature respectively, but is a constant reference pressure. Such model is suitable for current analysis, because in the presented storage device, gas flow is very slow and therefore pressure changes insignificantly in the fluid domain (absolute pressure changes by less than 0.3% throughout the analysis domain). Therefore it would be unjustified to use fully compressible flow model. The storage operating pressure was set to 2 bar, and this was the value of static pressure used for outlet boundary condition. Inlet boundary was set to constant flow velocity of 0,611 m/s for hot (800K) side inlet during charging process and 0,23 m/s inlet at the cold (300K) side inlet during storage discharge. Values of the velocities are matching a mass flow of 3.2x10 -3 kg/s for both charging and discharging stages of storage operation.
A geometry of two layers was created (shown on figure 3) and meshed with conformal periodicity on its maximum and minimum z coordinate. This mesh was then copied 16 times and merged into single, two zone mesh. Because of the shift between solid layers, no practical topology was determined for all hexa structured mesh. Therefore, solid region and contacting with it region of gas, are meshed as a structured multiblock meshes with a layer of pyramid and tetra elements connecting structured regions (see figure 4).   4. Coarse mesh detail view with exposed unstructured region joining two multiblock structured fluid regions.

Numerical setup
The analysis was performed with Fluent 18.1 software, with the use of pressure based coupled solver (coupled momentum and pressure-based continuity equations, turbulence model and energy equations solved in segregated manner). Two equation k--SST turbulence model was used. Second order bounded implicit scheme was chosen for time integration.
Meshes of two resolutions were prepared. A fine mesh, with boundary layers first element y+  1 (based on average velocities in the holes) having in total 14.2 million and second having y+ of around 15, having in total 2.1 million elements. Complete simulation was done on the coarse mesh.
Convective terms were discretized using first order upwind schemes, this together with coarse mesh approach was aimed at obtaining a highly diffusive setup, which would damp out small scale flow structures and allow for practical computation time.
The need for highly diffusive convective schemes and coarse mesh approach, was discovered after first attempts on computing the flow on the fine mesh. For this case, flow field oscillations are characterised by many orders of magnitude smaller time scales, compared to heat accumulation process in the solid. Therefore allowing for those oscillations to be resolved, requires impractically small time steps to be used when compared to the simulated process time of 17 hours. On the fine mesh, time step convergence was possible with given setup for time steps of the order of 10-20 milliseconds. After simulating first minute of the flow, it was estimated that the analysis would take around 32 years on 120 cores (almost 34 million CPU hours). Switching to coarse mesh resulted in numerical damping of the velocity field oscillations and allowed for two orders of magnitude higher timesteps. At the beginning of each storage operation stage, the time step was set to 1e-04 seconds and then gradually increased to 2.5 s which was the timestep size for major part of each charge/discharge stage. The simulation with that setup was finished in 13 days on 48 cores. To assess the inaccuracy of that approach, average heat fluxes absorbed by first five layers of the storage in first 50 seconds of the charging cycle (so in the part of cycle when highest temperature gradients are present so highest numerical errors are expected) where calculated on both meshes. Those average heat fluxes where 526.8 W and 479.2 W for fine and coarse meshes respectively. This 9% heat flux reduction, was found small enough to accept this approach because of its preliminary character aimed to assess how the storage core behaves and not to perform a detailed design.

Results and discussion
During analysis runtime, for each time step area averaged pressure and mass flow averaged temperature was recorded on inlet and outlet boundaries. Those data were later used to calculate inlet and outlet exergy fluxes as: Where: ̇ -exergy flux, Tref -reference temperature of 300K, pref -reference (operating) pressure of 2 bar, Cpspecific heat capacity at constant pressure, R -individual gas constant.
Exergy fluxes integrated in time, where then used to calculate amounts of exergy absorbed and recovered from the storage in each charge/discharge cycle to assess exergetic efficiency.
The time dependence of mass flow averaged temperatures at both storage ends is presented on figure 5. Constant temperature parts, represent prescribed values when given side acts as an inlet. Constant temperature of 800 K was assigned for hot side inlet for charging processes, and 300 K was assigned for cold side inlet for discharging processes. It can be seen, that except first cycle when hot gas is ingested into completely cold storage, outlet temperatures behave almost linearly, hence the heat flux exchanged with the storage drops in that almost linear manner throughout entire process. This observation suggests that the storage accumulates heat in its whole volume, instead of having a moving temperature front. This conclusion was also supported by recorded temperature profiles which are detailed in related article [1]. As can be also seen from figure 5, times of the subsequent processes are not equal. The storage heat capacity and mass flow, where designed for a 3.5 hours charge/discharge processes, that was the time of the first storage charging. During first discharge however, the outlet temperature was decreasing faster than expected and the discharge was stopped after 3 hours when outlet temperature dropped to almost 500 K. Second charge was also performed for 3 hours. After that it was decided (because first cycle exergy efficiency was calculated as 52.2% which is quite low) that the condition for process end should be the outlet temperature reaching 550K, which is halfway between charge/discharge inlet range instead of using prescribed stage time. This adaptive process allowed to increase the exergy efficiencies of subsequent cycles from 52.2% for first cycle, through 60.7% for second to 66.8% for the last operating cycle.
Main purpose of this study was to find qualitative results for gas-solid thermal interaction in conjugate heat transfer of the generic layered bed thermal storage. Moreover, limiting factors for the CFD pipeline of such devices.
The study in question combines two processes with large discrepancies in time scale: flow fluctuations occurring in order of milliseconds and heating of solids measured in hours. Numerical approach capable of time accurate capturing of flow fluctuations requires high density meshes with proper boundary layer resolution and leads to extremely slow computation and impractical walltime in order of magnitude of years. Coarsening the mesh resolution both in fluid and solid leads to filtering out the small time scale fluid flow fluctuations which leads to stabilisation of the numerical calculation and reducing the computation walltime to order of magnitude of hours. This is achieved at the cost of underestimation of conjugate heat flux.
It was also found, that the lower part of storage capacity is utilized during the thermal cycling, the higher exergy efficiency is obtained. The most useful conclusion however, is that the biggest temperature drop on the way of heat transfer between gas and solid is in the inner part of the gas boundary layer (closest to the solid wall, where the heat transfer coefficient is limited by molecular thermal conductivity of the gas). This leads to a conclusion, that in order to increase the thermal storage exergetic efficiency, layers (especially those nearest to its ends) should be of a much higher specific surface (gassolid interface surface per solid mass). It is further suggested that the use of granular beds, or porous materials in applications requiring high exergy efficiency might be more feasible.