A novel multi-reactor system for thermochemical heat storage through detailed modeling of K2CO3 particles

Thermochemical energy storage (TCES) is becoming increasingly important in the energy transition, as it can effectively bridge the gap between renewable energy supply and demand. In this study, the reaction kinetics of K 2 CO 3 were characterized and validated. Based on this kinetic model, a numerical model of a packed bed of particles was developed using a coupled CFD-DEM approach. The results of the model were validated against experimental data of a particle bed, showing good agreement. The reaction rate of the system was found to be limited by the diffusion of water vapor into the material, which led to unsatisfactory performance on the bed scale due to significant temperature drop-offs. Although reducing particle size was found to be an effective way to improve system performance, practical concerns such as agglomeration and bed permeability limited its effectiveness. As an alternative, a multi-reactor system with adaptive flow rates was proposed, which improved system performance without the limitations of reducing particle size. The proposed modular system is capable of delivering 10 kW power at the temperature of 45 degrees for a duration of 19.5 h.


Introduction
With the world shifting towards sustainable and renewable energy sources to mitigate the impact of global warming, the demand for affordable and efficient energy storage solutions is on the rise.Until now, the focus has been largely on solid state batteries.However, their low energy density, high cost and reliance on rare elements are prohibitive to large scale adoption, resulting in a need for effective alternatives [1].Thermal energy storage (TES) is an attractive option for this, since more than 50% of the energy demand in buildings is projected to account for heating and cooling [2].However, current TES solutions for residential heating relying on sensible and latent heat have limited effectiveness due to their low energy density and high thermal losses [3].To make TES more accessible and affordable, there is a growing need for innovative solutions that offer high energy density, low cost and minimal thermal losses.Thermochemical energy storage (TCES) offers a promising solution.
TCES utilizes a thermochemical material (TCM) that stores heat through a reversible reaction.The process involves an endothermic reaction that stores heat, which can then be released on-demand through an exothermic reaction.The salt hydrate potassium carbonate (K 2 CO 3 ) is a TCM that has been identified as a strong candidate for residential heating applications due to its high energy density, suitable operating conditions and economic viability [4].The anhydrate form of K 2 CO 3 reacts with water vapor to form the sesquihydrate K 2 CO 3 .1.5H 2 O, One application of K 2 CO 3 is in a fixed bed open system.This application is specifically of interest for residential applications due to the simple design with few points of failure.In this system the heat transfer fluid (HTF), typically air or nitrogen gas, flows through a packed bed reactor.Temperature and humidity of the fluid at the inlet is controlled to facilitate the desired reaction, causing the HTF to heat up during the exothermic hydration reaction, and cool down during the endothermic dehydration reaction.
For this purpose, the use of K 2 CO 3 in its simplest form, which is a raw, powder-like state, is not recommended for two reasons.Firstly, permeability is lower due to the small particles [5,6].This causes a high pressure drop over the bed, resulting in a high compressive power requirement for the system.Secondly, the particles agglomerate as they get hydrated.This further obstructs mass flow through the bed, increasing the pressure drop and causing local variations in operating conditions and thus uneven bed conversion.
By processing the material into mm-sized particles agglomeration can be prevented and bed permeability can be increased.This causes issues however in terms of mechanical stability; cycling the material https://doi.org/10.1016/j.est.2023.110028Received 3 August 2023; Received in revised form 18 October 2023; Accepted 5 December 2023 causes repeated swelling and shrinking, forming cracks in the material and eventually breaking it down into smaller particles [7].Efforts have been made to increase the mechanical stability of K 2 CO 3 , as well as other salt hydrates which face similar difficulties, by using additives.Common additives are expanded graphite (EG) [8] and expanded vermiculite (EV) [9,10].These additives show improved cyclic stability at the cost of a reduction in energy density due to displacement of the TCM.
To optimize a system using such particles an accurate mathematical model of the packed bed is desired.The complex interplay between physical and chemical phenomena on the particle scale and how this translates to the bed scale makes the development of an accurate model a difficult undertaking.Spatially resolved numerical models are the preferred method of modeling the temporo-spatial evolution of state variables for this type of system [11].Although many such models have been made using various modeling methods, TCMs and reactor designs, these bed-scale models often make significant assumptions regarding the inner-particle physics, which are simplified into various constitutive relations.
A novel method of modeling the system is the extended discrete element method (XDEM), which was introduced and developed by Peters in 2003 [12], and first applied in TCES by Mahmoudi et al. in 2021 [13].In this method each salt particle is modeled separately using the discrete element method (DEM), allowing for internal gradients in the particle properties.This simulation is coupled with computational fluid dynamics (CFD) to model the continuous phase outside of the particles.This way a large amount of detail is considered on the particle scale, leading to more accurate results on the bed scale.
In this study the XDEM method will be applied to the hydration of an open system reactor.First, the hydration reaction kinetics of K 2 CO 3 will be characterized and validated using simultaneous thermal analysis (STA) measurements, for the first time accounting for cycling effects and close to equilibrium behavior.The particle model is then used to create a bed model which is validated with experimental measurements provided by an external party (Cellcius [14]).Furthermore, a set of criteria is constructed for residential heating applications which is used to analyze and optimize the system.Lastly, we introduce a novel approach involving a multi-reactor system with adaptive flow rates.This method aims to maintain a consistent outlet gas temperature over an extended duration while preserving the heating power.

Theory
The hydration of K 2 CO 3 has been found to take place as a twostep process; (I) adsorption of water vapor on the particle surface, (II) the hydration reaction of Eq. ( 1) [15].When the conditions are close to the equilibrium conditions, adsorption is very slow due to low rates of nucleation.As the driving force increases nucleation rates increase rapidly, until adsorption is effectively instantaneous, making the hydration a one-step process.This adsorption limitation is known as metastability, and the region where this takes place is referred to as the metastable zone (MSZ).The p-T diagram of K 2 CO 3 is shown in Fig. 1, which includes the MSZ as found by Sögütoglu et al. [16], the deliquescence line and the saturation pressure.The equilibrium line is defined by the Van 't Hoff equation, which can be written as follows: With   [ ] the equilibrium pressure at temperature  [𝐾].  and   are constants that depend on the reaction enthalpy    and entropy    as follows:

Table 1
Relations for the pressure driving force of the hydration.
Driving force mechanism ℎ() ) The rate of the hydration reaction can be parameterized into the general kinetic equation: With [−] the conversion,   [1∕] an unknown constant,  () the reaction model, ( ) the temperature dependency, and ℎ () the pressure driving force based on the water vapor pressure [] which relates to the equilibrium conditions.The conversion is defined as follows: (6) where ,   and    are the masses of the hydrate during, before, and after hydration.

Temperature dependency
For the temperature dependency of the reaction the Arrhenius equation is used [13,15,17]: With   [J∕mol] the activation energy, [J∕mol K] the gas constant and  [K] the local temperature.

Pressure driving force
Various forms of the pressure driving force have been proposed and used in the literature.An overview is given in Table 1.These include linear [18], logarithmic, inverse [13] and exponential [19] (Table 1).Different relations should be used for dehydration.
Research by Sögütoglu et al. [16] however has shown that the vapor pressure dependency of the hydration rate is not well characterized by these relations.Instead they observed a near-zero reaction rate close to the equilibrium pressure (∕  < 2 at 40 • C), after which the reaction rate increased linearly with the driving force ∕  .The point at which this takes place will be referred to as  * .Furthermore they observed significant induction times in the MSZ (∕  < 3 at 40 • C).It should be noted that the same study showed that inducing the hydration at higher supersaturations and then lowering it into the MSZ bypasses the long induction time.
Gaeini et al. [17] observed a linear relation between the pressure driving force ∕  and the reaction rate, which is incompatible with the previous observations.They however only studied the reaction rate for high driving forces (i.e.∕  > 6), whereas the deviations from this only become obvious close to the equilibrium conditions.

Reaction model
The reaction model  () is an empirical relation which is used to model the decrease in reaction rate as the material becomes saturated.An extensive overview of reaction models for solid-state reactions is given by Vyazovkin et al. [20].Sögütoglu et al. investigated these reaction models for the hydration of K 2 CO 3 and proposed the use of an Avrami-Erofeev model (Eq.( 8)) for conditions inside the metastable zone where the rate of hydration is limited by the nucleation, and a reaction-order model (Eq.( 9)) for when the chemical reaction is rate limiting [15].
Using the reaction order model, Mahmoudi et al. experimentally obtained an exponential constant  of 0.7 for the hydration [13] and Gaeini et al. found values of approximately 0.66 based on four different cases [17].These correspond roughly with the contracting sphere model.However, these were only tested for cases with high pressure driving forces, raising concerns as to whether these values hold true for lower pressure driving forces (i.e.close to and inside the MSZ).

Numerical modeling methods
In the XDEM method every particle is modeled individually and includes the options for particle motion through dynamics calculations, conversion through chemical reactions, heat and mass transfer with the environment, and heat conduction and gas diffusion through the porous particles.Interaction with the environment is done through coupling with CFD.Inner particle heat and mass transfer uses radially distributed cells.Inter particle heat transfer is also considered through conduction and radiation.
The energy balance within the particle follows from the conservation of energy in the radial direction: ω   (10) with  the density,   the specific heat,  the temperature,    the effective thermal conductivity, ω the reaction rate of reaction k and   the reaction enthalpy of reaction k.Advective heat transfer through gas flow within the pore volume is neglected.Heat transfer at the boundaries is given by: with qconv the rate of convective heat transfer as given by the Wakao relation [21], qrad the radiative heat transfer as calculated by the Stefan-Boltzmann law, and qcond the conductive heat transfer based on [22] using the following relation: With  the thermal conductivity and  the temperature of the particle  and the neighboring particle . , is the conduction length, which is half of the radius of the computational cell.
The species mass equation of a gaseous species within the particle pore volume follows from the conservation equation: with   the particle porosity, ⟨  , ⟩  the partial density of gaseous species , ⟨   ⟩ the gas velocity,   the diffusion coefficient and ω,, the change in density of species  through reaction .
The momentum conversation equation is given by Darcy's Law: Mass transfer at the boundaries is given by: with   the particle surface area, ṁ the rate of diffusion as given by the Wakao relation [21], and ṁ the advective mass flow rate which follows from the momentum equation.
To model the hydration reaction given by Eq. ( 5), the reaction model, temperature dependency and pressure driving force must be obtained, for which STA measurements will be conducted.The details of these measurements are given in Section 3.2.

STA setup
For the characterization of the K 2 CO 3 grains and the validation of the mm-sized particles, a combination of Thermo-Gravimetric Analysis (TGA) and Differential Scanning Calorimetry (DSC) is used: Simultaneous Thermal Analysis (STA).The STA device (Netzsch STA 449 F3 Jupiter) is connected to a humidifier and a dry air gas supply.This way the mass change and the heat of reaction can be measured during the hydration and dehydration.The direction of the reaction is controlled by setting the relative humidity and the temperature of the inlet gas.The accuracy of the humidifier is ±0.8% relative humidity.The TGA sensor has an accuracy of ±0.030 mg.
The sample is tested in an alumina crucible (ø6.8 mm, 85 μL) which is weighed before it is put on the TGA-DSC sensor in the copper STA furnace.
For the characterization of the K 2 CO 3 grains around 10mg of material is placed in the crucible.The sample is then subjected to consecutive hydration and dehydration cycles.Hydration is performed at a temperature of 35 • C and a relative humidity of 30% for 180 min.Dehydration is performed at 120 • C in pure nitrogen over a period of 80 min.Once hydration rates have stabilized additional hydration cycles are performed at humidities of 17.7, 10.5, 6.4 and 4.8%.A heating rate of 5K/min and a cooling rate of 2.5 K/min are used.A 30 min rest period is included after the dehydration stage to ensure the temperature has stabilized before hydration.In line with prior research [23] which found hydration rates to be independent of nitrogen flow rates, we used a constant flow rate of 200 mL/min.

Reactor performance assessment
To optimize the design of a TCES reactor an effective way of assessing the performance is needed.For this purpose, a set of clear and quantifiable assessment criteria is constructed.The values for these criteria will depend on the use case of the system.These criteria will be judged for a residential heating case, the parameters of which are shown in Table 2.The criteria are as follows: 1.Energy density: the energy content per unit volume.Depends on the hydration enthalpy, material density, particle porosity and bed porosity.2. Specific power: the amount of thermal heating power which can be supplied per unit volume.3. Functional volume: measured in percentage of total volume, this is the bed conversion percentage over which the bed can deliver a specified outlet temperature.This is based on research by Zhang et al. [24], who investigated the outlet temperature profiles of various studies from literature, each of which shows a drop off as the system becomes saturated with water.4. Startup time: time for the bed to heat up before the specified outlet gas temperature is reached.5. Electrical efficiency: the amount of electrical energy is required to drive the flow compared to the thermal energy that is obtained.
The main criteria under investigation are the specific power and functional volume.Additionally, we take into account the startup time and electrical efficiency, although they are not the primary focus of this study.Notably, the energy density, being largely independent of the parameters under examination, has been excluded from our current investigation.By emphasizing specific power and functional volume as key criteria, our research aims to optimize TCES reactor design and operation, ensuring efficient heat delivery and utilization within space-constrained environments.

Characterization
In this section the results for the characterization of the kinetic equation for the hydration of K 2 CO 3 crystals are given, resulting in the kinetic equation with models for  () and ℎ(): Before characterizing the kinetics of the sample it is cycled until stable performance across cycles is observed.Fig. 2 shows the effects of this, with the reaction kinetics drastically improving over the first cycles before stabilizing at around 10 cycles.These results are consistent with research by Beving et al. who found a similar stabilization period of around 12 cycles [7].

Reaction model
The reaction model  () determines the change in reaction rate with the conversion state of the particle.When considering a numerical model with dynamic operating conditions, a generic model which is valid over the entire operating range is desired.For this purpose, the reaction-order model from Eq. ( 9) will be considered for all operating conditions.The reaction model is determined by hydrating the cycled sample at constant temperature and humidity conditions, which causes the reaction rate to depend solely on the conversion state (Eq.( 19)).Fig. 3 shows the resulting relation for three relative humidities.A reaction order model with an exponent factor of 0.3 represents the data well, regardless of the pressure driving force or sample used, thus resulting in Eq. ( 20).This is much lower than in previous studies [13,

𝑑𝛼 𝑑𝑡
During the initial part of the hydration, however, the provided model fails to provide an accurate prediction.This is due to the assumption of instant induction of the reaction which is not valid for the given material, especially when operating within the metastable region.This however is not a product of the conversion state of the particle but rather of the state of nucleation.This in turn depends on the driving force of the reaction and the time given for nucleation [16].As such, when the pressure driving force is low (i.e. <   ) an additional term (, ) can be included.which depends on the time since induction start  and the induction time :  The value of (, ) is 0 initially ( = 0) and gradually increases to 1 at full induction ( ≥ ).The shape of this function depends on the induction time and characteristics, which depend on the pressure driving force and the temperature [16].As of yet, no mathematical formulation exists in literature to define this.In this study instant induction is assumed ( = 0, (, ) = 1).

Pressure driving force
By altering the experimental reaction rate using the reaction model from Section 4.1.1 an adjusted reaction rate is obtained (Eq.( 22)), which is independent of the conversion state and thus constant during hydration.Table 3 shows the obtained values for different vapor pressure conditions at a temperature of 35 • C. A linear relationship between the pressure and the reaction rate is observed down to a pressure  * , below which no reaction was observed during the 180 min hydration cycle, with   <  * <   .This result is consistent with research by Sögütoglu et al. [16].The relation for the pressure driving force thus becomes Eq. ( 23), with  * =   + 0.75 (see Fig. 4).

( 𝑑𝛼 𝑑𝑡
) Based on these results we can identify three different zones for the hydration: I.   <  <  * ∶ hydration rate is near zero.II. * <  <   ∶ hydration rate is proportional to ∕ * − 1, but with significant induction times.III.  <  ∶ hydration rate is proportional to ∕ * − 1 and without induction times.

Kinetic parameter estimation
The activation energy of the Arrhenius equation is 34,828 J/mol [13].Using this value in combination with the reaction model and the  24)) with experimental data.
pressure driving force relation, the kinetic parameter of the kinetic equation is   = 0.0074 [1∕min].by replacing all the functions and constants in Eq. ( 18), Eq. ( 24) can be written as: This reaction model is compared with the experimental data in Fig. 5.For the experiments with relative humidities of 30.0% and 17.7% the model fits the data well.For a relative humidity of 10.5% the experimental data shows a delayed response due to the induction time which was not included in the kinetic model.

Particle validation
In order to validate the hydration equation derived in Section 4.1.3on the particle scale, a DEM model is employed to simulate the hydration process of a mm-sized porous K 2 CO 3 particle.The model is then compared with STA experiments conducted on mm-sized particles provided by Cellcius [14].The used particle has a slightly non-spherical shape with an effective diameter of 4.43 mm and an initial mass of 96.28 mg.The porosity of the particle is assumed to be variable during cycling.The particle is subjected to six consecutive hydration cycles under a relative humidity of 30%.To facilitate comparison, a 4.43 mm spherical particle is modeled using the DEM simulation, employing 100 radial cells to ensure high radial resolution.The particle is initially composed of 97w% K 2 CO 3 and 3w% inert additive, and the particle size and porosity remain constant throughout the simulation, neglecting the effects of volume change.
The comparison between experimental data and the model is presented in Fig. 6.Since our model does not include the effect of volume change, the measurement data after 10 cycles were considered, from which the change in the volume is negligible.Very good agreement is achieved in subsequent cycles with the 13% porosity in the DEM model.
The conversion cross-section of a DEM particle (  = 4 mm,   = 13%) during a hydration cycle is shown in Fig. 7(a).A thin reaction front can be observed, which steadily moves inwards towards the center of the particle.The water vapor pressure shown in Fig. 7(b) exhibits a gradient in the vapor pressure going from the boundary to the reaction front.Based on this, the hydration rate of the mm-sized particle is largely dependent on the rate of diffusion into the particle, rather than the reaction rate of the material.This observation is consistent with previous experimental findings by Aarts et al. [25].

Bed model
In this section the particle model from the previous section will be applied to a packed bed system.This model then will be validated by the data provided by Cellcius [14], the details from which are outlined in a previous work (Kieskamp et al. 2022 [26]), and the bed's performance will be assessed for different configurations.
The presented system involves a cylindrical packed bed reactor operating within an open system, shown in Fig. 8.In this system the air flow is used for both vapor and heat transport.Humid air is introduced at the inlet and subsequently flows through the bed, leading to the consumption of water vapor by the salt.This process results in the heating of both the salt and the air.The heated air then exits the system, after which the heat is extracted through a heat exchanger.Afterwards the air is rehumidified in an evaporator/condenser before reaching the inlet again.Flow is facilitated by a ventilator to ensure efficient operation of the entire system.dehydration the process is reversed by supplying heat instead of extracting it.
The data consists of 13 cycles of hydration and dehydration, performed on a cylindrical packed bed with a diameter of 68 mm and a height of 120 mm.The initial bulk porosity is 42%.During the 10hour hydration process, humid air with a flow velocity of 0.42 m∕s and a vapor pressure of 14 mbar is used.The inlet air temperature of the experiment starts at 40 • C and drops to 32 • C over the first few hours.
In this paper, we have replicated this setup with a CFD-DEM model.The particles are initially fully dehydrated (97w% K 2 CO 3 , 3w% inert additive) at a temperature of 40 • C. Air with a water vapor pressure of  H 2 O = 14 mbar is blown in through the top at  N 2 = 0.42 m∕s, with an inlet temperature identical to that of the experiment.The pressure at the bottom of the bed (the outlet) is   = 1 bar, and no transport of mass, energy, or momentum through the walls is included in the model.By omitting the consideration of heat transport through the wall, radial variations in system properties such as air temperature and vapor pressure are effectively eliminated, leaving changes in the axial direction as the sole focus.
Fig. 9 compares the resulting outlet air temperature and vapor pressure for the model and the last experiment cycle versus the overall bed conversion (Eq.( 6)).The model and experiment show good agreement, with the outlet vapor pressure of the model closely following the experiment.However, the experiment's outlet temperature is lower than that of the model, indicating a significant influence of thermal losses not accounted for in the model.
Both the experimental data and the computational model show that the rate of vapor consumption has a gradual increase over time during the initial phase of hydration.This can be attributed to a decreasing inlet air temperature in the early stages of the process, which brings the system closer to equilibrium and thereby results in reduced reaction rates.
Based on the results of the bed validation, it can be concluded that the CFD-DEM model accurately captures the behavior of the packed bed system.However, the discrepancy in outlet air temperature between the experiment and the model highlights the need to include thermal losses in the model.This is particularly important when designing a larger scale system, as the impact of thermal losses would be more significant.

Performance investigation
The drop in the outlet air temperature observed in both the experiments and the model is a considerable drawback of the system.Based   on the requirements set in Table 2, the system should be able to provide a temperature increase of at least 10 • C. In the experiments however, the observed increase drops below this at a bed conversion of 50%, thereby drastically reducing the functional volume and thus the energy density of the system.To investigate the cause of this phenomenon it is useful to look at the reaction rate profile.Fig. 10 shows the reaction rate as a function of the bed height and the time.A reaction wave profile similar to those described by Lin et al. [27] and Huinink et al. [28] is found.Initially the reaction mostly takes place at the inlet.As the material becomes saturated over time, the wave travels along the length of the reactor.Slow diffusion of water vapor into the particles however results in a wide reaction profile that does not reach equilibrium at the outlet of the 13 cm reactor.
To improve the performance of the system the diameter of the particles in the model is reduced to 1 mm.Fig. 11 compares the power densities as a function of the overall bed conversion, which is calculated by dividing the heat output with the bed volume (Eq.( 25)).A much more consistent power density profile is observed when using the smaller particles.This happens because the individual particles react much faster compared to 4 mm particles.This also corresponds to an increased volumetric power density of the system.
The reaction rate profile for the system using   = 1 mm is presented in Fig. 12.The profile exhibits a consistent reaction wave moving through the system after an initial startup period.In contrast to the   = 4 mm particle bed shown in Fig. 10, where the reaction takes place over the entire bed, the thin profile in Fig. 12 indicates a more localized reaction.This behavior is attributed to the increased reaction rate resulting from using smaller particles.Moreover, the full hydration cycle time is reduced by half compared to the original system due to the enhanced reaction rate.
Fig. 13 illustrates the reaction rate profile of the   = 1 mm system at  = 4000 s, revealing a distinct bell-shaped curve.On the inlet side, where particles are reaching complete conversion, the reaction rates are observed to be low.This phenomenon can be attributed to the limitation in vapor diffusion rates caused by the approaching full conversion, consequently decelerating the overall reaction kinetics.Conversely, as we progress towards the outlet, the conversion state of the particles decreases, leading to an increase in reaction rates.However, the gradual reduction in water vapor concentration along this trajectory impedes the reaction kinetics, resulting in the observed decline in reaction rates.At the end of the profile the reaction has mostly stopped, having achieved equilibrium conditions.
Based on these results, it can be concluded that the reactor outlet conditions will be at the equilibrium state as long as the reaction wave remains fully within the bed.Only once the reaction wave has reached the outlet will the outlet air temperature start to drop.
In a recent study conducted by Huinink et al. [28], a relationship was proposed for predicting the width of a reaction front:  With  the width of the reaction front,   a model dependant constant,   the bed bulk porosity,  2 the particle radius,  and  the velocities of the air and the hydration front respectively, and   the water vapor diffusivity inside the porous particle.
Applying this equation, the width of the reaction profiles for the   = 4 mm and   = 1 mm systems were determined to be 724 mm and 45 mm, respectively.It is noteworthy that for the numerical results of the   = 1 mm system approximately 78.9% of the reaction occurs within this width.However, it should be acknowledged that due to the asymptotic nature of the obtained results an exact measurement of the reaction profile's width cannot be provided.

Reactor optimization
Based on the findings from the performance investigation in the previous section, it is evident that the investigated particle bed with large (  = 4 mm) particles is unable to maintain the outlet conditions required for residential use for an extended period.This is because the outlet air temperature drops below the desired 40 • C at a bed conversion of only 40%.Here we first provide an overview of parameters that can be altered to achieve a better performing system.
• Particle size: Lowering the particle size to   = 1 mm was found to double the duration of maintaining the desired outlet air temperature, but it also poses significant risks of agglomeration and results in an increased pressure drop over the bed.• Reactor length: The length of the reactor can be increased such that the aforementioned reaction wave profile stays fully within the bed, ensuring equilibrium conditions at the outlet.This however reduces the power density of the system, and also increases the pressure drop over the bed.Furthermore there may also be spacial constraints limiting the reactor length.• Particle porosity: Increasing the particle porosity would allow for higher rates of diffusion into the particle, improving reaction kinetics and thus resulting in an overall thinner reaction profile that maintains outlet conditions for a longer time.However, mechanical stability of the particles is a major concern, as this change could lead to the particles breaking down into smaller particles, which would result in other complications as discussed previously • Flow rate: Another possibility is to operate the system at lower flow rates.The lower flow rates extract less heat from the particles and allow for more time for the diffusion of water vapor into the particle, narrowing the reaction profile and thus providing the desired outlet air temperature for a longer period.However, this also results in a reduction in heating power, which is directly proportional to the flow velocity.In addition, the startup time of the system increases proportionally, which may not be ideal for residential heating applications where heat should be readily available in a limited space.
To combat the performance loss without balancing the trade-offs that come with the aforementioned parameter variations, we present a novel approach involving a multi-reactor system with adaptive flow rates.This system aims to sustain a consistent outlet air temperature over an extended duration, while simultaneously preserving the heating power.The proposed method is visually depicted in Fig. 15, and its operational principles can be outlined as follows: 1. Multiple reactors are arranged in parallel configuration.2. The flow rate through each reactor is independently controlled, and only the minimum number of reactors necessary to meet the required demands are activated at any given time.3.As the outlet air temperature of an active reactor approaches a predefined threshold, the flow rate through that specific reactor is reduced to sustain the desired temperature level.4. To compensate for the decrease in heating power resulting from the reduced flow rate, another reactor is simultaneously activated in parallel.5.This process is iterated until no reactors remain inactive, indicating that the maximum system conversion has been achieved.
This approach offers the advantage of maintaining a constant outlet air temperature and high heating power, while also enabling efficient operation as the bed becomes saturated with water.Thus, the proposed multi-reactor system with adaptive flow rates presents a promising solution for achieving sustained and efficient operation.In addition to maintaining a steady outlet air temperature for an extended period, this approach offers modularity and reduces sensible heat losses during frequent short hydration-dehydration cycles by incorporating multiple reactors.
This approach is applied to the numerical model from Section 4.2.As shown in Fig. 14, the local air temperature of the   = 4 mm reactor with adaptive flow rate remains steady at 45 • C for the bulk of the hydration.However, a full hydration cycle for the system with adaptive flow rate takes twice as long.By extension, Fig. 11 illustrates a faster decrease in volumetric power density with the overall bed conversion.
For the case study described in Section 3.3, the TCM is distributed over five equal-sized reactor modules, resulting in a TCM content of 100 L per bed.To achieve the 10 kW power requirement using a single reactor, a power density of 0.1 kW∕L is required.This is reached at 50.8% conversion of the bed, or 130 min of continuous operation.Once this point is reached, the next module is turned on to cover the remaining heating load.This process is repeated for all modules until the system is exhausted.No thermal losses are taken into account.The resulting load distribution for the studied system over 25 h is presented in Fig. 16.The system's performance was evaluated based on its ability to provide a continuous heat output of 10 kW at a temperature of 45 • C for a prolonged period.The system achieved an operational time of 1171 min, during which 90.2% of the TCM has converted.This marks a significant enhancement over a similar modular system without adaptive flow rate reduction, in which the outlet air temperature and power were found to drop much faster

Conclusion
In this work, we have investigated the potential of thermochemical energy storage using K 2 CO 3 as the storage material.We have characterized the reaction kinetics of the material on the particle scale, and developed a numerical model of a packed bed of particles using a coupled CFD-DEM approach.We have validated the model against experimental data, and used it to assess the performance of the bed for different bed configurations.
Our results indicate that the reaction rate of the system is limited by the diffusion of water vapor into the material, and that reducing particle size is an effective way of increasing the performance of the system.However, practical concerns of agglomeration and bed permeability limit the effectiveness of this approach.As an alternative, we propose a multi-reactor system with adaptive flow rates, which improves the performance without the limitations of particle size reduction.
In future work, we would like to extend the characterization study in order to include the induction period in the kinetic equation, as well as incorporate a larger variety of conditions.Furthermore, further studies should also include dehydration, which was excluded in this study.

Fig. 2 .
Fig. 2. Conversion of a 14.96 mg K 2 CO 3 sample over the first 10 hydration cycles.

Fig. 3 .
Fig. 3. Normalized hydration rates compared with the obtained reaction order model at  = 35 • C.
17], likely due to the precycling of the samples improving the fluid pathways within the crystalline structure.

Fig. 4 .
Fig. 4. p-T diagram of K 2 CO 3 featuring MSZ boundaries as found by Sögütoglu et al. [16], with the inclusion of the induction transition region  * .

Fig. 6 .
Fig. 6.Comparison with experimental data of the conversion for a 4.43 mm particle.

Fig. 7 .
Fig. 7. Reaction front for a spherical DEM particle with   = 4 mm,   = 13% showing (a) the local molar H 2 O uptake and (b) the local H 2 O vapor pressure.

Fig. 8 .
Fig. 8. Schematic of the packed bed reactor filled with K 2 CO 3 particles.The bed height and diameter are 120 mm and 68 mm respectively.The humid air has a superficial velocity of 0.42 m∕s.

B
.Kieskamp et al.

Fig. 9 .
Fig. 9. Outlet air temperature and vapor pressure comparison model and experiment.

Fig. 10 .
Fig. 10.Reaction rate profile of the original system with   = 4 mm.The plot shows the changes of the average reaction rate in the axial direction as a function of time.

Fig. 11 .
Fig. 11.Volumetric power density as a function of bed conversion for particle beds consisting of 1 mm and 4 mm particles, as well as an adaptive flow velocity discussed in Section 4.3.

Fig. 12 .
Fig. 12. Reaction rate profile of the system with   = 1 mm.The plot shows the changes of the average reaction rate in the axial direction as a function of time.

Fig. 13 .
Fig. 13.Reaction rate, conversion, and vapor pressure along the length of the packed bed reactor with   = 1 mm, evaluated at  = 4000 s.

Fig. 14 .
Fig. 14.Air temperature profile of the   = 4 mm system with an adaptive flow rate in the axial direction as a function of time.Humid air is supplied at the bottom and hot, dry air flows out at the top.

Fig. 15 .
Fig. 15.Schematic overview of the multiple module TCES system using adaptive flow velocities.

Fig. 16 .
Fig. 16.Load distribution of five 100L variable flow rate reactor modules for a continuous heating load of 10 kW.

Table 2
Reference values for the system performance assessment, based on a typical Dutch home with floor heating.

Table 3
Conversion rate ∕ as a function of the pressure.