Shaolin Liu
email: shaolin.liu@u-bordeaux.fr

Azita Ahmadi-Senichault

Cyril Levet

Jean Lachaud

Development and validation of a local thermal non-equilibrium model for high-temperature thermal energy storage in packed beds

Keywords:

High-temperature thermal energy storage (TES) in packed beds is gaining interest for industrial energy recovery. The wide range of temperature distributions causes significant variations in thermophysical properties of the fluid and solid phases, leading to inaccuracies of classical TES models and heat transfer correlations. The objective of this work is to develop and validate a detailed but pragmatic model accounting for high-temperature effects. Based on a literature survey spanning over several communities interested in high-temperature porous media, we propose a generic local thermal non-equilibrium model for granulate porous media accounting for conservation of mass, momentum and energy (two-equation temperature model). The effective parameters needed to inform the model are the effective thermal conductivities of the different phases and the heat transfer coefficient. An experimental-numerical inverse analysis method is employed to determine these parameters. A dedicated experimental facility has been designed and built to study a model granulate made of glass bead of 16 mm diameter. Experiments are conducted using the Transient Single-Blow Technique (TSBT) by passing hot air (ranging from 293 K to 630 K) through cold particles at various mass flow rates, covering a Reynolds number range of 58 to 252. The new model was implemented in the Porous material Analysis Toolbox based on OpenFoam (PATO) used to compute the transient temperature fields. Two optimization algorithms were employed to determine the parameters by minimizing the error between experimental and simulated temperatures: a Latin Hypercube Sampling (LHS) method, and a local optimization method Adaptive nonlinear least-squares algorithm (NL2SOL). The results indicate that the value of heat transfer coefficient ℎ 𝑣 in the two-equation model falls in the range of 1.0 × 10 4 ∼ 1.93 × 10 4 W/(m 3 ⋅ K) under the given conditions. The axial dispersion gas thermal conductivity was found to be around 5.9 and 67.1 times higher than the gas thermal conductivity at Peclet numbers of around 55 and 165, respectively. Furthermore, two improved correlations of Nusselt number (𝑁𝑢 = 2 + 1.54𝑅𝑒(𝑇 ) 0.6 𝑃 𝑟(𝑇 ) 1∕3 ) and of axial dispersion gas thermal conductivity (𝑘 𝑑𝑖𝑠,|| = 0.00053𝑅𝑒(𝑇 ) 2.21 𝑃 𝑟(𝑇 ) ⋅ 𝑘 𝑔 ) are proposed and validated for a range of Reynolds number from 58 to 252. The overall approach is therefore validated for the model granulate of the study, opening new perspectives towards more precise design and monitoring of high-temperature TES systems.

Introduction

Thermal energy storage (TES) is receiving increased attention for the development of energy recovery technologies such as concentrated solar power (CSP) [START_REF] Esence | A review on experience feedback and numerical modeling of packed-bed thermal energy storage systems[END_REF][START_REF] Liu | Experimental study on heat storage and transfer characteristics of supercritical air in a rock bed[END_REF][START_REF] Anderson | Experimental results and modeling of energy storage and recovery in a packed bed of alumina particles[END_REF], advanced adiabatic compressed air energy storage (AACAES) [START_REF] Liu | Experimental study on heat storage and transfer characteristics of supercritical air in a rock bed[END_REF][START_REF] Yang | Theoretical evaluation on the impact of heat exchanger in advanced adiabatic compressed air energy storage system[END_REF][START_REF] Mozayeni | Performance study of an advanced adiabatic compressed air energy storage system[END_REF], and industrial waste heat recovery [START_REF] Alva | An overview of thermal energy storage systems[END_REF]. Amongst the major technologies (sensible, latent and chemical energy storage) [START_REF] Oró | Review on phase change materials (pcms) for cold thermal energy storage applications[END_REF][START_REF] Agyenim | A review of materials, heat transfer and phase change problem formulation for latent heat thermal energy storage systems (lhtess)[END_REF][START_REF] Gil | State of the art on high temperature thermal energy storage for power generation. part 1-concepts, materials and modellization[END_REF][START_REF] Yang | Thermal performance of a shell-and-tube latent heat thermal energy storage unit: Role of annular fins[END_REF], sensible heat storage is the most developed and used in the industry [START_REF] Esence | A review on experience feedback and numerical modeling of packed-bed thermal energy storage systems[END_REF][START_REF] Xu | Extending the validity of lumped capacitance method for large biot number in thermal storage application[END_REF][START_REF] Hänchen | High-temperature thermal storage using a packed bed of rocks -heat transfer analysis and experimental validation[END_REF]. Most industrial sensible heat storage systems rely on packed beds that offer advantages such as a large heat transfer area, low cost, and a wide operating range (293 K to 1000 K) [START_REF] Liu | Experimental study on heat storage and transfer characteristics of supercritical air in a rock bed[END_REF][START_REF] Qu | Heat transfer characteristics of mixed convection in packed beds[END_REF].

To optimize the design of these systems, the coupled heat and mass transfer between the carrier fluid and the packed-bed material needs to be well understood and modeled. Many authors have contributed to this topic [START_REF] Esence | A review on experience feedback and numerical modeling of packed-bed thermal energy storage systems[END_REF][START_REF] Liu | Experimental study on heat storage and transfer characteristics of supercritical air in a rock bed[END_REF][START_REF] Anderson | Experimental results and modeling of energy storage and recovery in a packed bed of alumina particles[END_REF][START_REF] Hänchen | High-temperature thermal storage using a packed bed of rocks -heat transfer analysis and experimental validation[END_REF][START_REF] Qu | Heat transfer characteristics of mixed convection in packed beds[END_REF][START_REF] Schumann | Heat transfer: A liquid flowing through a porous prism[END_REF][START_REF] Pomeroy | Thermal energy storage in a packed bed of iron spheres with liquid sodium coolant[END_REF][START_REF] Saez | Dynamic response of a packed bed thermal storage system-a model for solar air heating[END_REF][START_REF] Wakao | Effect of fluid dispersion coefficients on particle-to-fluid heat transfer coefficients in packed beds: Correlation of nusselt numbers[END_REF][START_REF] Zanganeh | Packed-bed thermal storage for concentrated solar powerpilot-scale demonstration and industrial-scale design[END_REF][START_REF] Zanganeh | Design of packed bed thermal energy storage systems for high-temperature industrial process heat[END_REF][START_REF] Anderson | Packed bed thermal energy storage: A simplified experimentally validated model[END_REF]. Heat transfer in a packed bed can either be studied under the assumption of local thermal equilibrium (LTE) or local thermal non-equilibrium (LTNE) [START_REF] Esence | A review on experience feedback and numerical modeling of packed-bed thermal energy storage systems[END_REF]. These models are also referred to as one-equation and two-equation models, respectively. Pati et al. [START_REF] Pati | Critical review on local thermal equilibrium and local thermal nonequilibrium approaches for the analysis of forced convective flow through porous media[END_REF] reviewed the applicability of one-and two-equation models. The choice of model is influenced by parameters such as the gas-to-solid thermal diffusivity ratio (𝛼 𝑔 ∕𝛼 𝑠 ), Reynolds number (𝑅𝑒 = 𝜌 𝑔 𝑢𝑑 𝑝𝑎𝑟 ∕𝜇 𝑔 ), and Darcy number (𝐷𝑎 = 𝐾∕𝑑 2 𝑝𝑎𝑟 ). Specifically, when the Reynolds number is greater than 50 and the thermal diffusivity ratio between the gas and solid phases (𝛼 𝑔 ∕𝛼 𝑠 ) is higher than 30 [START_REF] Quintard | Introduction to heat and mass transport in porous media[END_REF], the two-equation model is required. Schumann [START_REF] Schumann | Heat transfer: A liquid flowing through a porous prism[END_REF] was the first to propose a two-equation energy model -one for carrier fluid and one for the packed-bed material; however, the heat capacity of the fluid and the thermal conductivities of both phases were neglected, and the model was one-dimensional. Saez et al. [START_REF] Saez | Dynamic response of a packed bed thermal storage system-a model for solar air heating[END_REF] flow direction (increase in conductivity due to the gas flow). The thermophysical properties were assumed to be constant because the inlet gas temperature was at most 70°C. Hanchen et al. [START_REF] Hänchen | High-temperature thermal storage using a packed bed of rocks -heat transfer analysis and experimental validation[END_REF] applied the model proposed by Saez et al. to high-temperature conditions (527°C) but neglected the dispersion term and assumed constant gas properties. Zanganeh et al. [START_REF] Zanganeh | Packed-bed thermal storage for concentrated solar powerpilot-scale demonstration and industrial-scale design[END_REF][START_REF] Zanganeh | Design of packed bed thermal energy storage systems for high-temperature industrial process heat[END_REF] proposed a more comprehensive model that considered the radiative term in the effective gas thermal conductivity and variable thermophysical properties for a packed bed of rocks in the temperature range of 20-650°C. Anderson et al. [START_REF] Anderson | Experimental results and modeling of energy storage and recovery in a packed bed of alumina particles[END_REF][START_REF] Anderson | Packed bed thermal energy storage: A simplified experimentally validated model[END_REF] further validated the two-equation one-dimensional model with varying thermophysical properties for a packed bed of alumina particles at both low (120°C) and high (700°C) temperatures. Ma et al. [START_REF] Ma | Numerical and experimental studies of packed bed thermal energy storage system based on a novel transient energy model[END_REF] developed a transient twoequation two-dimensional model that fully considers thermal dispersion and changes in thermophysical properties, expanding on the existing literature. In these works, the carrier fluid density and velocity are considered constant, which is a correct approximation for incompressible fluids (water, molten nitrate salts, ...). Using air as the carrier fluid, non-uniform temperature distributions inside packed beds leads to nonuniform gas density and viscosity, affecting the velocity and pressure distributions. We aim to fill this gap by accounting for mass and momentum conservation, namely by including the volume-averaged conservation of mass and momentum (Darcy-Forchheimer's law) into the model. Such models have been employed since the 1960's for very high-temperature applications such as space-vehicle heat shields [START_REF] Lachaud | A generic local thermal equilibrium model for porous reactive materials submitted to high temperatures[END_REF]. Also, in previous TES studies, the dimensionless numbers (Reynolds, Prandtl) used in the correlations were assumed to be constant. We would like to discuss this assumption when variations of velocity and thermophysical properties are not negligible. The two-equation energy model writes as follows [START_REF] Quintard | Two-medium treatment of heat transfer in porous media: numerical results for effective properties[END_REF][START_REF] Liu | Experimental investigation on the validity of the local thermal equilibrium assumption in ablative-material response models[END_REF]:

𝜕 𝜕𝑡 ( 𝜀 𝑠 𝜌 𝑠 𝑐 𝑝,𝑠 ⟨𝑇 𝑠 ⟩ 𝑠 ) = ∇ ⋅ (𝐤 𝐬,𝐞𝐟𝐟 ⋅ ∇⟨𝑇 𝑠 ⟩ 𝑠 ) + ℎ 𝑣 ( ⟨𝑇 𝑔 ⟩ 𝑔 -⟨𝑇 𝑠 ⟩ 𝑠 ) (1)
𝜕 𝜕𝑡

( 𝜀 𝑔 𝜌 𝑔 𝑐 𝑝,𝑔 ⟨𝑇 𝑔 ⟩ 𝑔 ) + ∇ ⋅ (𝜌 𝑔 𝑐 𝑝,𝑔 ⟨𝐯 𝐠 ⟩ 𝐠 ⟨𝑇 𝑔 ⟩ 𝑔 ) = ∇ ⋅ (𝐤 𝐠,𝐞𝐟𝐟 ⋅ ∇⟨𝑇 𝑔 ⟩ 𝑔 ) +ℎ 𝑣 ( ⟨𝑇 𝑠 ⟩ 𝑠 -⟨𝑇 𝑔 ⟩ 𝑔 ) (2)
where the gas volume fraction 𝜀 𝑔 is equal to the porosity of the packed bed, ℎ 𝑣 is the volumetric heat transfer coefficient, 𝑐 𝑝,𝑖 and 𝐤 𝐢,𝐞𝐟𝐟 respectively denote the heat capacity and effective thermal conductivity tensor of the 𝑖 phase, ⟨𝑇 𝑖 ⟩ 𝑖 and ⟨𝐯 𝐠 ⟩ 𝐠 denote the intrinsic phase average temperature for the 𝑖-phase and the intrinsic gas average velocity. For the development of TES models, three effective properties need to be evaluated: ℎ 𝑣 , 𝐤 𝐬,𝐞𝐟𝐟 and 𝐤 𝐠,𝐞𝐟𝐟 .

The heat transfer coefficient and corresponding correlation can be determined either experimentally [START_REF] Qu | Heat transfer characteristics of mixed convection in packed beds[END_REF][START_REF] Wakao | Effect of fluid dispersion coefficients on particle-to-fluid heat transfer coefficients in packed beds: Correlation of nusselt numbers[END_REF][START_REF] Wakao | Heat and mass transfer in packed bedsgordon and breach science publications[END_REF][START_REF] Ranz | Friction and transfer coefficients for single particles and packed beds[END_REF][START_REF] Galloway | A model of the mechanism of transport in packed, distended, and fluidized beds[END_REF][START_REF] Yang | Experimental analysis of forced convective heat transfer in novel structured packed beds of particles[END_REF][START_REF] Naghash | Measurement of convective heat transfer coefficients in a randomly packed bed of silica gel particles using ihtp analysis[END_REF][START_REF] Nazari | Experimental investigation of heat transfer and second law analysis in a pebble bed channel with internal heat generation[END_REF] or through pore-scale simulation in representative geometries [START_REF] Quintard | Two-medium treatment of heat transfer in porous media: numerical results for effective properties[END_REF][START_REF] Sassanis | Numerical determination of the volumetric heat transfer coefficient in fixed beds of wood chips[END_REF][START_REF] Augier | Numerical simulations of transfer and transport properties inside packed beds of spherical particles[END_REF]. The experimental methods include both steady state [START_REF] Ranz | Friction and transfer coefficients for single particles and packed beds[END_REF][START_REF] Galloway | A model of the mechanism of transport in packed, distended, and fluidized beds[END_REF][START_REF] Naghash | Measurement of convective heat transfer coefficients in a randomly packed bed of silica gel particles using ihtp analysis[END_REF][START_REF] Nazari | Experimental investigation of heat transfer and second law analysis in a pebble bed channel with internal heat generation[END_REF] and transient regime [START_REF] Qu | Heat transfer characteristics of mixed convection in packed beds[END_REF][START_REF] Yang | Experimental analysis of forced convective heat transfer in novel structured packed beds of particles[END_REF] approaches. With the steady state approach, particles within the packed beds are heated with a resistive heater. Compressed air, serving as the working fluid, extracts heat from these heated particles. After reaching steady-state, the temperatures of the gas flow at the inlet, outlet, within the packed bed, as well as the surface temperature of the particles are monitored. These measurements are used to estimate the heat exchange and infer the heat transfer coefficient. In the transient method, cold particles are heated using hot gas. The time evolution of average inlet, outlet, and inside gas flow temperatures, and particle temperatures, are recorded. The local heat transfer coefficient is derived from analyzing the spatial and time evolution of air and particle temperatures. The progress in optimization algorithms and computer resources has allowed the general usage of inverse methods to infer the parameters with more accuracy. Concerning the pore-scale simulation method, the volume averaging theory [START_REF] Quintard | Two-medium treatment of heat transfer in porous media: numerical results for effective properties[END_REF] have been applied to periodic structures, such are arrays of cylinders, to obtain ℎ 𝑣 . With the progress of the resolution of computed micro-tomography and computer resources, pore-scale simulations are becoming a realistic approach to determine ℎ 𝑣 [START_REF] Scandelli | Modélisation multi-échelle des milieux poreux réactifs à haute température appliquée aux protections thermique spatiale et à la pyrolyse de la biomasse[END_REF][START_REF] Petrasch | Tomography based determination of permeability, dupuit-forchheimer coefficient, and interfacial heat transfer coefficient in reticulate porous ceramics[END_REF]. Using these two methods, numerous correlations were obtained for the heat transfer coefficient ℎ 𝑣 in packed-bed under diverse conditions. A summary of the most widely used correlations and of their validity ranges is proposed in Table 1. However, the influence of temperature variations on these correlations is not specified. Taking the Wakao correlation as an example, it is based on published experimental data in the range of 3 < 𝑅𝑒 < 10000, mostly based on intermediate temperature conditions and without providing the temperature range of validity. Most simulation models for high-temperature storage systems used this correlation for temperatures up to 1173 K [START_REF] Liu | Experimental study on heat storage and transfer characteristics of supercritical air in a rock bed[END_REF][START_REF] Ma | Numerical and experimental studies of packed bed thermal energy storage system based on a novel transient energy model[END_REF][START_REF] Jalalzadeh-Azar | Heat Transfer in a High-Temperature Packed Bed Thermal Energy Storage System-Roles of Radiation and Intraparticle Conduction[END_REF][START_REF] Cascetta | A comparison between cfd simulation and experimental investigation of a packed-bed thermal energy storage system[END_REF]. At high temperatures, a wide temperature distribution can cause nonuniform gas density and viscosity inside the packed bed. These varying thermal physical properties affect 𝑅𝑒, which in turn affects 𝑁𝑢 in the Wakao correlation. Temperature-dependent parameters should replace constant parameters in such correlations; also the linear coefficients (1.1 in Wakao's correlation) should be validated or re-evaluated.

Concerning the effective gas and solid thermal conductivities, a large amount of theoretical work has been done by theoreticians and there is now a good agreement that they respectively write [START_REF] Quintard | Two-medium treatment of heat transfer in porous media: numerical results for effective properties[END_REF][START_REF] Degroot | Closure of non-equilibrium volume-averaged energy equations in high-conductivity porous media[END_REF][START_REF] Yang | A synthesis of tortuosity and dispersion in effective thermal conductivity of porous media[END_REF]

where 𝑘 𝑖 denotes the thermal conductivity of the 𝑖 phase, 𝑉 is a representative elementary volume of the porous medium, 𝐴 𝑔𝑠 is the area of the g-s interface contained in the averaging volume, 𝐧 𝐠𝐬 is the outwardly directed unit normal vector pointing from the gas phase toward the solid phase, 𝐧 𝐠𝐬 = -𝐧 𝐬𝐠 . 𝐯 𝑔 , ⟨𝐯 𝐠 ⟩ 𝐠 , and ṽ𝑔 =𝐯 𝑔 -⟨𝐯 𝑔 ⟩ 𝑔 are respectively the gas velocity within the pores, the intrinsic phase average velocity, and the deviation velocity. The closure variables 𝐛 𝐢𝐢 (𝑖 = 𝑔, 𝑠) are the vector fields that map ∇⟨𝑇 𝑖 ⟩ 𝑖 onto T𝑖 , where T𝑖 is the spatial deviation temperature in the 𝑖-phase, more precisely [START_REF] Quintard | Two-medium treatment of heat transfer in porous media: numerical results for effective properties[END_REF], T𝑔 = 𝐛 𝐠𝐠 ⋅ ∇⟨𝑇 𝑔 ⟩ 𝑔 -𝑠 𝑔 (⟨𝑇 𝑔 ⟩ 𝑔 -⟨𝑇 𝑠 ⟩ 𝑠 )+..., T𝑠 = 𝐛 𝐬𝐬 ⋅∇⟨𝑇 𝑠 ⟩ 𝑠 +𝑠 𝑠 (⟨𝑇 𝑠 ⟩ 𝑠 -⟨𝑇 𝑔 ⟩ 𝑔 )+.... Note that in these derivations the effective radiative conductivity may be added as a contribution the effective solid conductivity for optically thick media [START_REF] Van Eekelen | Numerical validation of an effective radiation heat transfer model for fiber preforms[END_REF]; this hypothesis is verified for packed beds as the grain size is small compared to the system. In many research works, regardless of the experimental conditions, simple expressions of 𝐤 𝐠,𝐞𝐟𝐟 and 𝐤 𝐬,𝐞𝐟𝐟 as a function of gas and solid phase conductivities and the porosity of the homogeneous and isotropic packed beds have been considered [START_REF] Hänchen | High-temperature thermal storage using a packed bed of rocks -heat transfer analysis and experimental validation[END_REF][START_REF] Anderson | Packed bed thermal energy storage: A simplified experimentally validated model[END_REF][START_REF] Abdulagatova | Effect of temperature and pressure on the thermal conductivity of sandstone[END_REF][START_REF] Lin Xia | Experiment on the convective heat transfer from airflow to skeleton in open-cell porous foams[END_REF]. In Eq.3, 𝐤 𝐠,𝐞𝐟𝐟 is composed of three parts: the conductivity term 𝐤 𝐠,𝐜𝐨𝐧 accounts for a simple average of the thermal conductivity, the tortuosity term 𝐤 𝐠,𝐭𝐨𝐫 represent the tortuosity of the phase repartition and the dispersion term 𝐤 𝐝𝐢𝐬 accounts for gas flow effects. 𝐤 𝐠,𝐭𝐨𝐫 and 𝐤 𝐝𝐢𝐬 could be numerically estimated by solving a closure problem on a periodic unit cell representative of the structure [START_REF] Quintard | Two-medium treatment of heat transfer in porous media: numerical results for effective properties[END_REF][START_REF] Degroot | Closure of non-equilibrium volume-averaged energy equations in high-conductivity porous media[END_REF]. 𝐤 𝐠,𝐭𝐨𝐫 is often neglected since 𝐤 𝐝𝐢𝐬 dominates [START_REF] Degroot | Closure of non-equilibrium volume-averaged energy equations in high-conductivity porous media[END_REF]. Wakao et al. [START_REF] Wakao | Heat and mass transfer in packed bedsgordon and breach science publications[END_REF] and DeGroot et al. [START_REF] Degroot | Closure of non-equilibrium volume-averaged energy equations in high-conductivity porous media[END_REF] described 𝐤 𝐝𝐢𝐬 in the axial and radial directions (respectively horizontal and transverse flow directions) as different coefficients multiplied with Reynolds 𝑅𝑒 and Prandtl 𝑃 𝑟 numbers. In Eq.4, 𝐤 𝐬,𝐞𝐟𝐟 consists of the conductivity of the solid 𝐤 𝐬,𝐜𝐨𝐧 corrected to account for its tortuosity 𝐤 𝐬,𝐭𝐨𝐫 [START_REF] Van Antwerpen | A review of correlations to model the packing structure and effective thermal conductivity in packed beds of mono-sized spherical particles[END_REF], whereas the closure variable 𝐛 𝐬𝐬 which is a function of the gas velocity 𝐯 𝐠 [START_REF] Quintard | Two-medium treatment of heat transfer in porous media: numerical results for effective properties[END_REF] integrates dispersion effects. Therefore, for high temperatures, the best approach is to use the models of Eq.3 and 4 with temperature-dependent thermal conductivities 𝑘 𝑔 and 𝑘 𝑠 .

In this context, the purpose of this article is to propose a generic two-equation model coupled with mass and momentum conservation for high-temperature conditions and to determine the effective parameters needed to inform the model, namely, the thermal conductivities and the heat transfer coefficient. For this purpose, we have designed and built a laboratory-scale heat storage facility where the gas temperatures at the inlet, at the outlet and within the packed bed are monitored. Glass beads with a diameter of 16 mm were selected as the heat storage medium. In Section 2, the experimental setup and test procedure are presented. In section 3, we present a numerical inverse analysis method that we have developed to infer the intrinsic parameters and that consists of a detailed numerical model of the packed bed and an advanced multi-objective optimization algorithm. In section 4, the experimental results and the estimated values of the two-equation parameters are presented and new correlations accounting for high-temperature effects are proposed and compared with the literature. Conclusions and outlook are discussed in section 5.

Experimental facility

The main characteristics of gas/solid packed-bed TES described in the literature have been reviewed by Esence et al [START_REF] Esence | A review on experience feedback and numerical modeling of packed-bed thermal energy storage systems[END_REF]. These encompass a mix of large-scale industrial systems [START_REF] Zanganeh | Packed-bed thermal storage for concentrated solar powerpilot-scale demonstration and industrial-scale design[END_REF][START_REF] Zanganeh | Design of packed bed thermal energy storage systems for high-temperature industrial process heat[END_REF] and a majority of reducedscale laboratory set-ups [START_REF] Qu | Heat transfer characteristics of mixed convection in packed beds[END_REF][START_REF] Ma | Numerical and experimental studies of packed bed thermal energy storage system based on a novel transient energy model[END_REF][START_REF] Yang | Experimental analysis of forced convective heat transfer in novel structured packed beds of particles[END_REF][START_REF] Okello | Experimental investigation of thermal de-stratification in rock bed tes systems for high temperature applications[END_REF][START_REF] Shitzer | Transient Behavior of a Rock-Bed Thermal Storage System Subjected to Variable Inlet Air Temperatures: Analysis and Experimentation[END_REF]. In studies conducted at the laboratory scale, the diameter ratios (tube diameter to particle size) range from 3 to 32, with gas flow rates varying between 0.035 to 3.23 m/s [START_REF] Esence | A review on experience feedback and numerical modeling of packed-bed thermal energy storage systems[END_REF]. To minimize the effect of diameter ratios on the distribution of fluid velocity, it is recommended that diameter ratio remains greater than 10 [START_REF] Halkarni | Estimation of volumetric heat transfer coefficient in randomly packed beds of uniform sized spheres with water as working medium[END_REF]. In most cases, the velocity of the gases is of the order 0.3 m/s [START_REF] Qu | Heat transfer characteristics of mixed convection in packed beds[END_REF][START_REF] Okello | Experimental investigation of thermal de-stratification in rock bed tes systems for high temperature applications[END_REF][START_REF] Shitzer | Transient Behavior of a Rock-Bed Thermal Storage System Subjected to Variable Inlet Air Temperatures: Analysis and Experimentation[END_REF]. The experimental facility presented in the following section was designed in alignment with these characteristics. It is based on a method called the transient single-blow technique (TSBT) [START_REF] Esence | A review on experience feedback and numerical modeling of packed-bed thermal energy storage systems[END_REF]. With this approach, hot gas is used to heat particles inside packed beds and the time evolution of the gas temperature is recorded at the inlet, at the outlet and also at different locations within the packed bed to bring robustness to the inverse analysis.

Experimental system and test procedure

A schematic drawing and photograph of the experimental setup that we have developed and an enlarged view of the test section, where thermocouple positions are shown, are displayed in Fig. 1. The setup consists of a hot air blower (BAK Thermoplastic) to supply and heat air, an insulated test section made of an iron tube and two metallic grids to maintain the particles, a differential pressure transmitter (EMER-SON FISHER ROSEMOUNT) to measure pressure drop and a thermoanemometer (VT 110-2014 THERMO-ANEMOMETER) to measure the gas velocity. The iron tube has dimensions of 940 mm in length, 194 mm in inner diameter, and a thickness of 3 mm. These dimensions are suitable to study glass beads with a diameter of 16 mm, while respecting the diameter aspect ratio of at least 10, recalled in the introduction.

In the experiment, the first step is to set the flow rate of the inlet gas. This parameter is controlled by the air blower and measured using the thermo-anemometer. The maximum flow rate is 800 liters/min, equivalent to a mass flow rate of 1.6 × 10 -2 kg/s. The maximum Darcy velocity 𝑢 (also known as superficial gas velocity, 𝑢 = 𝑞𝑚∕(𝜌 𝑔 𝜋𝑅 2 𝑏 )) reaches 0.451 m/s at room temperature (293 K), where 𝜌 𝑔 is the gas density of 1.2 kg/s and 𝑅 𝑏 is the radius of the tube at 0.097 m. Once the inlet air mass flow rate is determined at room temperature, the next step involves setting the heating power, which is also controlled by the air blower. To ensure data reproducibility, we repeated each experiment following a consistent procedure. After an initial experiment, with a twoday cooldown period to return particles to room temperature, subsequent runs exhibited only minor temperature differences (1-2 • C) initially due to the room temperature. After 15 minutes of heating, the temperature differences remained minimal, within 0.3 • C, confirming data reliability across trials.

The black points in Fig. 1 indicate the positions of the thermocouples (type K, 2 mm sheath diameter). Temperature measurements are taken at various locations: the inlet gas temperature (1) at 𝑧=0 mm, 𝑟=0 mm; the gas temperature inside the test section at seven central axis locations (from 2 to 8) at 𝑧=50, 150, 250, 350, 450, 650, and 850 mm, 𝑟=0 mm; and three radial locations (from 9 to 11) at 𝑧=50, 150, and 250 mm, 𝑟=-45 mm. The tube surface temperature is also measured at four axis locations (from 12 to 15) at 𝑥=50, 350, 650, and 900 mm, 𝑧=100 mm.

To ensure accurate gas temperature measurement, the thermocouple is securely fixed in the gap between the glass beads using ST-1000 thermal glue, as shown in Fig. 2. The thermal glue is capable of withstanding temperatures up to 950°C.

The thermocouples are connected to two display data loggers (Pico Technology TC-08) that record the temperature with a time step of 1.0 s. Experimental variations in the inlet air mass-flow rate and heating power were performed to verify the model's suitability for different Reynolds numbers (𝑅𝑒). The test section dimensions and the architectural properties of particles are listed in Table 2. The gas volume fraction (porosity, 𝜀 𝑔 ) in Table 2 is calculated using the Muller correlation [START_REF] Halkarni | Estimation of volumetric heat transfer coefficient in randomly packed beds of uniform sized spheres with water as working medium[END_REF][START_REF] Mueller | Radial void fraction distributions in randomly packed fixed beds of uniformly sized spheres in cylindrical containers[END_REF] as shown in Eq.5. The permeability and Forchheimer coefficient are calculated using the Ergun equation, as shown in Eq.6 [START_REF] Esence | A review on experience feedback and numerical modeling of packed-bed thermal energy storage systems[END_REF], where these two parameters are treated as scalars. These defined parameters will be used in the subsequent modeling section. Ergun equation is validated for modeling pressure variations in the experiments, as explained in Appendix A.

𝜀 𝑔 = 0.365 + 0.22∕(𝐷 𝑏 ∕𝑑 𝑝𝑎𝑟 ) (5)

Experimental uncertainty analysis

Experimental uncertainties can be categorized into two types: direct measurement errors resulting from the experimental device, and errors arising from the measurement process itself. In this work, the ranges and uncertainties associated with the experimental devices are presented in Table 3. Errors from the temperature measurement process are chiefly due to the positioning of the thermocouples in the tube. The thermocouples are fixed in the gap between the glass beads, which are positioned with a tolerance of ± 10 mm. The uncertainty in temperature measurement within the tube, attributed to the thermocouple positioning, is determined through repeated experiments, resulting in an uncertainty value of 1.19%.

Numerical inverse analysis method

To tackle the complexity of the problem, numerical inverse analysis was employed to estimate the desired quantities. Firstly, we present the model describing fluid flow and heat transfer within the test section. Secondly, we detail the multi-objective optimization algorithm used to infer the effective parameters.

Numerical model

To facilitate parameter estimation, we model the entire test section as shown in Fig. 1 (b). The corresponding numerical model, presented in Fig. 3, represents a two-dimensional axisymmetric geometry. The model consists of two regions, namely the packed bed (porous sample) and the tube (iron tube).

Mathematical model and numerical implementation

In the packed bed region, we consider the two-equation model (Eq.1 and 2) for energy conservation. For the tube region, transient heat conduction is considered (energy conservation in a solid). The Reynolds number based on the particle diameter (𝑑 𝑝𝑎𝑟 ) is defined as 𝑅𝑒 = 𝜌 𝑔 𝑢𝑑 𝑝𝑎𝑟 ∕𝜇 𝑔 . Throughout the entire experimental process, the minimum value of 𝑅𝑒 is 18.34 (corresponding to 𝑇 𝑔 =630 K, 𝜌 𝑔 =0.563 kg∕m 3 , 𝑢=0.067 m/s, 𝜇 𝑔 =3.29e-05 kg∕(m ⋅ s)), and the maximum value of 𝑅𝑒 is 453.36 (corresponding to 𝑇 𝑔 =293 K, 𝜌 𝑔 =1.2 kg∕m 3 , 𝑢=0.451 m/s, 𝜇 𝑔 =1.91e-05 kg∕(m ⋅ s)). Considering the range of applicability of the Forchheimer law inside a packed bed, which typically lies between 5 [START_REF] Sobieski | Darcy's and forchheimer's laws in practice. part 1. the experiment[END_REF] and 500 to 600 [START_REF] Montillet | Flow Through a Finite Packed Bed of Spheres: A Note on the Limit of Applicability of the Forchheimer-Type Equation[END_REF], the average gas velocity is determined by solving the Forchheimer law. The two regions are coupled by ensuring the conservation of heat fluxes at their interfaces. The mass and energy conservation equations for the two-equation model in the packed bed

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ 𝜕 𝜕𝑡 ( 𝜀 𝑔 𝑀⟨𝑝⟩ 𝑔 𝑅⟨𝑇 𝑔 ⟩ 𝑔 ) + ∇ ⋅ ( - 𝑀⟨𝑝⟩ 𝑔 𝐊𝐗 𝑅⟨𝑇 𝑔 ⟩ 𝑔 ⋅ ∇⟨𝑝⟩ 𝑔 ) = 0 (Mass equation) 𝜕 𝜕𝑡 ( 𝜀 𝑠 𝜌 𝑠 𝑐 𝑝,𝑠 ⟨𝑇 𝑠 ⟩ 𝑠 ) = ∇ ⋅ (𝐤 𝐬,𝐞𝐟 𝐟 ⋅ ∇⟨𝑇 𝑠 ⟩ 𝑠 ) + ℎ 𝑣 ( ⟨𝑇 𝑔 ⟩ 𝑔 -⟨𝑇 𝑠 ⟩ 𝑠 ) (Energy equation) 𝜕 𝜕𝑡 ( 𝜀 𝑔 𝜌 𝑔 𝑐 𝑝,𝑔 ⟨𝑇 𝑔 ⟩ 𝑔 ) + ∇ ⋅ (𝜌 𝑔 𝑐 𝑝,𝑔 ⟨𝐯 𝐠 ⟩ 𝐠 ⟨𝑇 𝑔 ⟩ 𝑔 ) = ∇ ⋅ (𝐤 𝐠,𝐞𝐟 𝐟 ⋅ ∇⟨𝑇 𝑔 ⟩ 𝑔 ) + ℎ 𝑣 ( ⟨𝑇 𝑠 ⟩ 𝑠 -⟨𝑇 𝑔 ⟩ 𝑔 ) (Energy equation) 𝜌 𝑡 𝑐 𝑝,𝑡 𝜕𝑇 𝑡 𝜕𝑡 = ∇ ⋅ (𝑘 𝑡 ⋅ ∇𝑇 𝑡 ) (Energy equation) (7)
where ⟨𝑝⟩ 𝑔 is the intrinsic average pressure, the mass equation is derived from the perfect gas law (Eq. 8), gas mass conservation (Eq. 9) [START_REF] Lachaud | A generic local thermal equilibrium model for porous reactive materials submitted to high temperatures[END_REF], and the Forchheimer equation (Eq. 10). The Forchheimer tensor, represented as 𝐗, is introduced to simplify the computation process [START_REF] Amao | Mathematical model for darcy forchheimer flow with applications to well performance analysis[END_REF].

𝜌 𝑔 = 𝑀⟨𝑝⟩ 𝑔 𝑅⟨𝑇 𝑔 ⟩ 𝑔 ( 8
)
𝜕 𝑡 ( 𝜀 𝑔 𝜌 𝑔 ) + ∇ ⋅ ( -𝜀 𝑔 𝜌 𝑔 ⟨𝐯 𝐠 ⟩ 𝐠 ) = 0 (9)
⟨𝐯 𝐠 ⟩ 𝐠 = - 1 𝜀 𝑔 (𝐊𝐗) ⋅ ∇⟨𝑝⟩ 𝑔 , 𝑋 𝑖𝑗 = 1 𝜇 𝑔 𝐾 𝑖𝑗 + 𝛽 𝑖𝑗 𝜌 𝑔 |⟨𝑣 𝑔 ⟩| (10)
It should be noted that glass beds exhibit significant opacity at a wavelength of 4.9 𝜇𝑚, which coincides with the peak intensity wavelength at 700 K. For 5 mm glass beds, radiation contributes to nearly 0% of the total heat transfer at temperatures below 700 K, while it increases to 50% at 1000 K [START_REF] Chen | Radiant heat transfer in packed beds[END_REF]. Therefore, it is reasonable to neglect radiation heat transfer effects when the temperature is below 700 K. The initial and boundary conditions for the systems of Eq.7 are described in detail in Table 4. In the packed bed region, 𝑇 𝑔 and 𝑇 𝑠 represent the gas and solid temperatures, respectively, while 𝑇 𝑡 represents the temperature in the tube region. The thermal conductivity of the tube is denoted as 𝑘 𝑡 . For simplicity, we use 𝑇 𝑔 , 𝑇 𝑠 , and 𝑝 to represent the average gas temperature ⟨𝑇 𝑔 ⟩ 𝑔 , average solid temperature ⟨𝑇 𝑠 ⟩ 𝑠 , and average pressure ⟨𝑝⟩ 𝑔 , respectively. The gas temperature 𝑇 𝑔 at the inlet side is determined from the experimental study. Numerically, we imposed the inlet mass flux using a pressure gradient inlet boundary condition based on Darcy-Forchheimer's law projected on the inlet patch normal 𝒏. The expression of the inlet pressure gradient imposed is provided in Table 4. In the tube region, the wall surface temperature 𝑇 𝑡 is measured throughout the experimental period. At the interface between the packed bed and the tube, the condition of 𝑇 𝑔 = 𝑇 𝑡 , 𝑘 𝑔,𝑒𝑓 𝑓 𝜕𝑇 𝑔 𝜕𝑟 = 𝑘 𝑡 𝜕𝑇 𝑡 𝜕𝑟 were applied for the gas phase. By implementing this boundary condition, which implies that the gas temperature 𝑇 𝑔 and the tube temperature 𝑇 𝑡 are equal at this interface, it also guarantees that the heat flux entering one region on one side of the interphase matches the heat flux leaving the other region on the opposite side of the domain. At the packed bed-tube interface, the condition of 𝜕𝑇 𝑠 ∕𝜕𝑟 = 0 was used for the solid phase. This boundary condition is in line with the practical situations and finds support in the existing literature [START_REF] Cascetta | A comparison between cfd simulation and experimental investigation of a packed-bed thermal energy storage system[END_REF]. A linear part variation of 𝑇 𝑡 within each interval is assumed in line with the findings reported by Cascetta et al. [START_REF] Cascetta | A comparison between cfd simulation and experimental investigation of a packed-bed thermal energy storage system[END_REF] in a similar experimental setup. The temperatures at the four measuring points (from 12 to 15) shown in Fig. 1 (b) are expressed as follows

⎧ ⎪ ⎨ ⎪ ⎩ 𝑇 𝑡 = 𝑇 12 + (𝑇 12 -𝑇 13)∕(𝑧13 -𝑧12) ⋅ (𝑧12 -𝑧) 0 < 𝑧 < 𝑧13 𝑇 𝑡 = 𝑇 13 + (𝑇 13 -𝑇 14)∕(𝑧14 -𝑧13) ⋅ (𝑧13 -𝑧) 𝑧13 ≤ 𝑧 < 𝑧14 𝑇 𝑡 = 𝑇 14 + (𝑇 14 -𝑇 15)∕(𝑧15 -𝑧14) ⋅ (𝑧14 -𝑧) 𝑧14 ≤ 𝑧 < 𝐿 (11)
Fig. 4 shows the temperature distribution of 𝑇 𝑡 on the wall surface, which is calculated using Eq.11. In this case, the inlet gas temperature 𝑇 𝑔 gradually increases from 292 K to 530K, and the inlet mass-flow rate 𝑞𝑚 is chosen as 5.97 × 10 -3 kg/s. As shown in Fig. 4, 𝑇 𝑡 exhibits a linear variation along the 𝑧-axis within each of the three zones as a function of time and space.

The numerical model was implemented in the Porous material Analysis Toolbox based on OpenFOAM (PATO) [START_REF] Lachaud | Porous-material analysis toolbox based on openfoam and applications[END_REF], using a finite volume method. A multi-block approach was employed, where different sets of equations were solved in each region at each time step, incorporating appropriate boundary conditions derived from the energy balances at the interfaces. A wedge-mesh for axisymmetrical simulations was generated using the blockMesh application of OpenFOAM. The mesh is refined at the near-tube region as shown in Fig. 5. In the packed bed region, the pressure equation was solved semi-implicitly using first-order schemes in time (Euler) and space (Gauss linear). Similarly, the energy equations were solved with an implicit treatment of the temperature terms and an explicit treatment of other quantities. In the tube region, only the energy equation needed to be solved, with the temperature terms treated implicitly and other quantities treated explicitly. The equations are solved in series, with appropriate mesh refinement and time steps to guarantee that the order of convergence is reached.

Table 4

Initial and boundary conditions of the two-equation model.

Physical properties of air and glass

This model takes into account the temperature-dependent properties of both the gas and solid phases. The density of air, 𝜌 𝑔 , is calculated using the perfect gas law, while the values for the heat capacity, 𝑐 𝑝,𝑔 , dynamic viscosity, 𝜇 𝑔 , and thermal conductivity, 𝑘 𝑔 , are obtained from the NASA-9 database [START_REF] Mcbride | Thermodynamic Data for Fifty Reference Elements[END_REF]. The density of glass, 𝜌 𝑠 , is assumed to be constant at 2500 kg/m 3 , while the heat capacity, 𝑐 𝑝,𝑠 , and thermal conductivity, 𝑘 𝑠 , are determined by fitting experimental data from Pertermann et al. [START_REF] Pertermann | Transport properties of low-sanidine single-crystals, glasses and melts at high temperature[END_REF] as described in Eq.12.

{

𝑘 𝑠 = 0.59206 + 0.00062𝑇 𝑠 + 1.0013 × 10 5 provides the temperature-dependent properties for both the gas and solid phases used in the numerical model. The values for each property are given at specific temperatures, and linear interpolation is employed between two adjacent values to estimate the property at intermediate temperatures.

As discussed in Section 1, the determination of three unknown effective parameters, namely the effective solid thermal conductivity 𝐤 𝐬,𝐞𝐟 𝐟 , effective gas thermal conductivity 𝐤 𝐠,𝐞𝐟𝐟 , and volumetric heat transfer coefficient ℎ 𝑣 , is required in the two-equation model (Eq.1 and 2). To facilitate the optimization algorithm, several dimensionless coefficients, including 𝑐 1 , 𝑐 2 , 𝑐 3 , and 𝑓 , are defined as presented in Eq.13.

⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ 𝐤 𝐬,𝐞𝐟𝐟 = 𝐤 𝐬,𝐜𝐨𝐧 + 𝐤 𝐬,𝐭𝐨𝐫 = 𝑐 1 ⋅ 𝑘 𝑠 𝐼 𝑘 𝑔,𝑒𝑓 𝑓 ,|| = 𝑘 𝑔,𝑐𝑜𝑛 + 𝑘 𝑔,𝑡𝑜𝑟,|| + 𝑘 𝑑𝑖𝑠,|| = 𝜀 𝑔 𝑘 𝑔 + 0 + 𝑐 2 ⋅ 𝑅𝑒 ⋅ 𝑃 𝑟 ⋅ 𝑘 𝑔 = 𝜀 𝑔 𝑘 𝑔 + 𝑐 2 ⋅ 𝑑 𝑝𝑎𝑟 ⋅ 𝜀 𝑔 ⋅ |⟨𝑣 𝑔 ⟩ 𝑔 | ⋅ 𝜌 𝑔 ⋅ 𝑐 𝑝,𝑔 𝑘 𝑔,𝑒𝑓 𝑓 ,⟂ = 𝑘 𝑔,𝑐𝑜𝑛 + 𝑘 𝑔,𝑡𝑜𝑟,⟂ + 𝑘 𝑑𝑖𝑠,⟂ = 𝜀 𝑔 𝑘 𝑔 + 0 + 𝑐 3 ⋅ 𝑅𝑒 ⋅ 𝑃 𝑟 ⋅ 𝑘 𝑔 = 𝜀 𝑔 𝑘 𝑔 + 𝑐 3 ⋅ 𝑑 𝑝𝑎𝑟 ⋅ 𝜀 𝑔 ⋅ |⟨𝑣 𝑔 ⟩ 𝑔 | ⋅ 𝜌 𝑔 ⋅ 𝑐 𝑝,𝑔 ℎ 𝑣 = (6𝜀 𝑠 ∕𝑑 𝑝𝑎𝑟 ) ⋅ (2 + 𝑓 𝑅𝑒 0.6 𝑃 𝑟 1∕3 )𝑘 𝑔 ∕𝑑 𝑝𝑎𝑟 (13)
The effective gas thermal conductivity, 𝐤 𝐠,𝐞𝐟𝐟 , is influenced by the tortuosity and dispersion terms, which have different values in the horizontal (||) and transverse (⟂) flow directions [START_REF] Esence | A review on experience feedback and numerical modeling of packed-bed thermal energy storage systems[END_REF][START_REF] Wakao | Heat and mass transfer in packed bedsgordon and breach science publications[END_REF]. Typically, the tortuosity term is neglected compared to the dispersion term [START_REF] Degroot | Closure of non-equilibrium volume-averaged energy equations in high-conductivity porous media[END_REF]. The relationship between 𝐤 𝐝𝐢𝐬 and the thermal properties of the gas is modeled similarly to the formulation proposed by Wakao et al. [START_REF] Wakao | Heat and mass transfer in packed bedsgordon and breach science publications[END_REF]. The volumetric heat transfer coefficient, ℎ 𝑣 , is modeled using the functional form proposed by Wakao et al. [START_REF] Wakao | Effect of fluid dispersion coefficients on particle-to-fluid heat transfer coefficients in packed beds: Correlation of nusselt numbers[END_REF]. The parameters 𝑐 1 , 𝑐 2 , 𝑐 3 , and 𝑓 will be optimized as described in section 3.2.

Optimization process

A numerical inverse analysis method was employed to infer the intrinsic parameters from experimental measurements in order to minimize errors between measured and predicted data. We have eleven data sequences

{𝑇 1 𝑖 } 𝑛 𝑖=0 , {𝑇 2 𝑖 } 𝑛 𝑖=0 , {𝑇 3 𝑖 } 𝑛 𝑖=0 , {𝑇 4 𝑖 } 𝑛 𝑖=0 , {𝑇 5 𝑖 } 𝑛 𝑖=0 , {𝑇 6 𝑖 } 𝑛 𝑖=0 ,

Table 5

Thermal properties of gas and solid as a function of temperature.

gas solid , which correspond to the measured gas temperatures collected at each time step indicated by the index 𝑖 (see Fig. 1 (b)). The transient inlet gas temperature data 𝑇 1 is fitted into a second-order polynomial using least squares and is used as a boundary condition. The objective function 𝑆 for the optimization processes is defined as the average of the root mean squared relative error between the measured and predicted temperatures at the ten positions mentioned (𝑇 2, 𝑇 3, 𝑇 4, 𝑇 5, 𝑇 6, 𝑇 7, 𝑇 8, 𝑇 9, 𝑇 10, 𝑇 11):

𝑇 (K) 𝑐 𝑝,𝑔 (J kg -1 K -1 ) 𝜇 𝑔 (kg m -1 s -1 ) 𝑘 𝑔 (W m -1 K -1 ) 𝑐 𝑝,𝑠 (J kg -1 K -1 ) 𝑘 𝑠 (W m -1 K -1 ) 273
𝑆 = 1 10 ⎛ ⎜ ⎜ ⎝ √ √ √ √ 1 𝑛 𝑛 ∑ 𝑖=1 ( 𝑇 2 𝑖 𝑛𝑢𝑚 -𝑇 2 𝑖 𝑇 2 𝑖 )2 + √ √ √ √ 1 𝑛 𝑛 ∑ 𝑖=1 ( 𝑇 3 𝑖 𝑛𝑢𝑚 -𝑇 3 𝑖 𝑇 3 𝑖 )2 +...... + √ √ √ √ 1 𝑛 𝑛 ∑ 𝑖=1 ( 𝑇 11 𝑖 𝑛𝑢𝑚 -𝑇 11 𝑖 𝑇 11 𝑖 )2 ⎞ ⎟ ⎟ ⎠ (14)
The minimization of 𝑆 was performed using the Design Analysis Kit for Optimization and Terascale Applications (Dakota) [START_REF] Adams | Dakota, a multilevel parallel object-oriented framework for design optimization[END_REF]. The optimization process involved a combination of the Latin hypercube sampling (LHS) method [START_REF] Olsson | On latin hypercube sampling for structural reliability analysis[END_REF] for sensitivity analysis and the Adaptive nonlinear least-squares algorithm (NL2SOL) [START_REF] Dennis | An adaptive nonlinear least-squares algorithm[END_REF] for local optimization. LHS was employed initially for sensitivity analysis to identify the most important input variables and their interactions. This stratified sampling technique divides the range of each uncertain variable into 𝑁 𝑠 segments of equal probability, resulting in a more accurate estimate of the mean value compared to random sampling. The next step involved local optimization using NL2SOL. This algorithm utilizes a quasi-Newton update to quickly converge to an optimal solution. By employing a simplification scheme to approximate the Hessian, NL2SOL achieves faster convergence compared to global optimization methods when performing least square calculations. An overview of this optimization strategy is presented in Fig. 6, which will be further explained in the subsequent section through examples of applications to two-equation model optimizations.

Results and discussions

In this section, we present the results of both measured and predicted temperatures. In the first subsection, we determine the volumetric heat transfer coefficient ℎ 𝑣 , the effective solid thermal conductivity 𝐤 𝐬,𝐞𝐟𝐟 , and the effective gas thermal conductivity 𝐤 𝐠,𝐞𝐟𝐟 used in the two-equation model. Additionally, we provide the predicted fields of gas properties, pressure, velocity, temperature, and volumetric heat transfer coefficient using the two-equation model. In the second subsection, we focus on evaluating the validity of the obtained effective parameters. Fig. 7 presents a comparison of the measured and predicted results for four mass-flow rates. The relative errors, represented by the objective function S, are below 10 -2 for all cases. Detailed data can be found in Table 7. A discussion of these results can be found in the following subsections. The measured results are represented by dots, while the predicted results obtained by solving the two-equation model are shown as solid and dashed lines. The recorded inlet gas temperature (𝑇 1), and gas temperature at ten locations (𝑇 2 to 𝑇 11) are plotted as a function of time. The presented temperature data represents the subset that was used in the optimization process. In the predicted results, gas temperatures at seven points along the central axis (𝑇 2 𝑛𝑢𝑚 to 𝑇 8 𝑛𝑢𝑚 ) are shown as solid lines, and the other three points (𝑇 9 𝑛𝑢𝑚 to 𝑇 11 𝑛𝑢𝑚 )) with dashed lines.

The experimental parameters at room temperature are shown in Table 6. In the first case, the inlet gas temperature gradually increases from room temperature to 630 K, the inlet gas mass-flow rate is 3.65 × 10 -3 kg/s, and the superficial gas velocity at room temperature is 0.10 m/s. The Reynolds number at room temperature is 101. The fluctuation in the inlet gas temperature might be caused by the tolerance of the thermocouple. Error bars indicate the uncertainty in the measured temperatures, as described in Section 2, where it was noted that the error due to the uncertain position of the thermocouples was estimated to be 1.19%. The major parameters that can affect the predicted temperatures are 𝐤 𝐬,𝐞𝐟𝐟 , 𝐤 𝐠,𝐞𝐟𝐟 , and ℎ 𝑣 . In the following, we present the procedure for determining these parameters.

Optimization and numerical results

In this subsection, we focus on determining the three effective parameters used in the two-equation model. Next, we present the numerical results, including the pressure, velocity, and temperature fields throughout the whole test section.

Determination of three effective parameters used in the two-equation model

As discussed in Subsection 3.2, the estimation of these effective parameters involves two modules: sensitivity analysis and the optimization process. To begin, a sensitivity analysis is performed to understand the relationship between the output result (𝑆) and the input variables (𝑓 , 𝑐 1 , 𝑐 2 , 𝑐 3 ). In this analysis, we use a high mass-flow rate case (𝑞𝑚 = 1.02 × 10 -2 kg/s) as an example to demonstrate the sensitivity analysis and optimization method. Additional case results are presented in Appendix B.

During the sensitivity analysis, it is important to select appropriate minimum and maximum values for the parameters to be optimized. To determine these values, we refer to the correlations provided by Esence et al. [START_REF] Esence | A review on experience feedback and numerical modeling of packed-bed thermal energy storage systems[END_REF] for packed beds. Based on their work, the coefficients 𝑓 , 𝑐 1 , 𝑐 2 , and 𝑐 3 are estimated to be 1.1, 0.617, 0.1, and 0.001, respectively, for Reynolds numbers less than 8000. In order to assess the sensitivity of the model to these parameters, we employ the Latin hypercube sampling (LHS) method. The ranges for 𝑓 , 𝑐 1 , 𝑐 2 , and 𝑐 3 are defined as follows: 𝑓 ranges from 0.7 to 2.5, 𝑐 1 ranges from 0.4 to 1.8, 𝑐 2 ranges from 0.01 to 0.5, and 𝑐 3 ranges from 0 to 0.01. The LHS method ensures that the input variables are evenly distributed within these minimum and maximum boundaries.

Fig. 8 presents a typical result of the sensitivity analysis. The red line represents the strength and direction of the relationship between the independent variables (𝑓 , 𝑐 1 , 𝑐 2 , and 𝑐 3 ) and the dependent variable (𝑆). The slope of the red line corresponds to the partial correlation coefficient (PCC), which quantifies the linear relationship between two variables (𝑥: 𝑓 or 𝑐 1 or 𝑐 2 ) and 𝑦 (𝑆) while accounting for the influence of one or more additional variables [START_REF] Adams | Dakota, a multilevel parallel object-oriented framework for design optimization[END_REF][START_REF] Helton | Sampling-based methods for uncertainty and sensitivity analysis[END_REF]. The PCC is calculated using the following formula:

PCC(𝑥, 𝑦) = ∑ 𝑖 ( 𝑥 𝑖 - x) ( 𝑦 𝑖 - ȳ) √ ∑ 𝑖 ( 𝑥 𝑖 - x) 2 ∑ 𝑖 ( 𝑦 𝑖 - ȳ) 2 x𝑖 = 𝑐 0 + 𝑛 ∑ 𝑝=1,𝑝≠𝑖 𝑐 𝑝 𝑥 𝑝 , ŷ = 𝑏 0 + 𝑛 ∑ 𝑝=1,𝑝≠𝑖
𝑏 𝑝 𝑥 𝑝 [START_REF] Pomeroy | Thermal energy storage in a packed bed of iron spheres with liquid sodium coolant[END_REF] where x and ȳ represent the mean values of 𝑥 and 𝑦. 𝑥 𝑖 and 𝑦 𝑖 denote individual values (samples) of the two variables. The regression equations x𝑖 and ŷ are defined by real number coefficients 𝑏 and 𝑐, and 𝑥 𝑝 represents the other variables except for 𝑥 𝑖 . PCC can help to identify the variables that have the most significant impact on the output results and should be prioritized in the optimization process. In this case, the PCC values between the input variables (𝑓 , 𝑐 1 , 𝑐 2 , and 𝑐 3 ) and the output variable (𝑆) are -0.12, -0.87, -0.22, and -0.04, respectively. A coefficient closer to ±1 indicates a strong correlation between the input and output variables. The results show a very low correlation between 𝑐 3 and 𝑆. Therefore, we shifted our focus to investigating the effects of the other three variables, 𝑓 , 𝑐 1 , and 𝑐 2 on the output result 𝑆.

During the optimization process, the boundary values for the three variables (𝑓 , 𝑐 1 , 𝑐 2 ) were set the same as in the sensitivity analysis: 𝑓 ranging from 0.7 to 2.5, 𝑐 1 ranging from 0.4 to 1.8, and 𝑐 2 ranging from 0.01 to 0.5. Convergence (𝑆 < 0.01) was reached after 50 iterations of the local optimization algorithm NL2SOL. Fig. 9 shows the influence of the three variables (𝑓 , 𝑐 1 , 𝑐 2 ) on the error 𝑆. The size of the points in Fig. 9 (a) represents the value of 𝑆, while the grey scale indicates the value of 𝑓 . As shown in Fig. 9 (a), when 𝑆 is below 0.01, the values of 𝑐 1 are concentrated around 1.29 and 𝑐 2 are concentrated around 0.42, while the values of 𝑓 range from 1.0 to 2.0. To further analyze the behavior of the error when varying 𝑓 , we set 𝑐 1 to 1.29 and 𝑐 2 to 0.42 and only varied 𝑓 as the input parameter. As shown in Fig. 9 (b), 𝑆 displays a decreasing and then increasing behavior. The minimum error 𝑆 is reached when the value of 𝑓 is 1.31 with an uncertainty of 5%. With the convergence of the error 𝑆 to the level of 10 -2 , we have determined the values of 𝑐 1 , 𝑐 2 , and 𝑓 to be 1.29, 0.42, and 1.31, respectively. Consequently, the corresponding values of ℎ 𝑣 , 𝑘 𝑠,𝑒𝑓 𝑓 , and 𝑘 𝑔,𝑒𝑓 𝑓 ,|| are established as well and summarized in Table 7.

As shown in Table 7, when considering the two-equation model coupled with the mass and momentum equations (Eq.7) and accounting for the temperature dependency of the gas and solid thermal properties, some coefficients exhibit different values compared to their original correlations. The coefficient 𝑐 1 , which is related to the solid phase volume fraction 𝜀 𝑠 , is found to increase from 0.617 to around 1.22. This value remains relatively constant across different inlet velocities. On the other hand, the coefficient 𝑐 2 , associated with gas dispersion, shows an increasing trend with higher inlet velocities and is not a constant value. This is in line with the fact that the dispersion effect increases with the 9.14 × 10 -3 , which is the average error across ten different locations (𝑇 2 to 𝑇 11). The maximum error among these locations is 1.41 × 10 -2 , For reference, the maximum temperature difference between the predicted and measured values is approximately 6.2 K. This indicates that the numerical model can accurately capture the transient temperature profiles in a packed bed with a high level of accuracy. To illustrate further the overall behavior of the packed bed, the upcoming subsection will present color maps of the simulation results.

Color maps of the simulation results in the full domain

In this subsection, we present a simulation result using the experimental conditions from case 1 at 4000 s. The inlet gas temperature 𝑇 𝑔 gradually increases from 292 K to 630 K within 4000 s. The inlet mass-flow rate 𝑞𝑚 is 3.56×10 -3 kg/s, resulting in a Darcy velocity 𝑢 of 0.1 m/s at room temperature. The values of 𝑐 1 , 𝑐 2 , 𝑐 3 , and 𝑓 are set to 1.21, 0.11, 0, and 1.61, respectively. These values were obtained from the optimization process in the previous subsubsection. Fig. 10 shows the distribution of dynamic viscosity 𝜇 𝑔 and air density 𝜌 𝑔 within the packed bed region at 4000 s. The thermodynamic properties of air are updated based on the pressure and gas temperature using a linear interpolation method, utilizing the data provided in Table 5. The coupling between air density and temperature is accounted for in the mass equation in Eq. 7. At the inlet of the packed bed, the air temperature is 630 K, accompanied by a dynamic viscosity of 3.3e-5 kg∕(m ⋅ s) and a density of 0.56 kg/m 3 . At the outlet, the air temperature remains at room temperature, resulting in a dynamic viscosity of 1.9e-5 kg∕(m ⋅ s) and a density of 1.20 kg/m 3 . The dynamic viscosity of air increases by 74% over the modeled temperature range of 290 K to 630 K, while the air density experiences a decrease of 53%. Fig. 12 shows the temperature fields in the packed bed region at 4000 s. The gas temperature 𝑇 𝑔 at the inlet is 630 K, while at the outlet, it is 300 K. The temperature of the gas and solid phases is updated by solving the energy equations in Eq.7. In the packed bed region, the absolute difference between the local average gas and solid temperatures (𝑑𝑇 = |𝑇 𝑔 -𝑇 𝑠 |) is around 4 K. The distribution of 𝑑𝑇 is influenced by the heat transfer coefficients ℎ 𝑣 . It is worth noting that in Fig. 12 (b) there is a significant temperature difference in the upper left region. This is primarily due to the strong coupling at the interface between the gas temperature inside the packed bed and the solid temperature at the tube surface, which is much lower than the solid temperature inside the packed bed near the interface. Fig. 13 shows the distribution of the volumetric heat transfer coefficients, ℎ 𝑣 , within the packed bed region at 4000 s. It is worth noting that ℎ 𝑣 is a function of space and time within the packed bed, as described by Eq. 13. The value of ℎ 𝑣 at the inlet is 13400 (W/m 3 ⋅ K) and at the outlet is 11400 (W/m 3 ⋅ K). Throughout the modeled temperature range of 290-630 K, ℎ 𝑣 experiences a decrease of 15 %, and its calculation incorporates the Wakao correlation with the factor 𝑓 . The distribution of ℎ 𝑣 aligns with the distribution of gas velocity, meaning that higher gas velocities correspond to higher ℎ 𝑣 values. To further investigate the transient distribution of ℎ 𝑣 within the packed bed, we plotted the values of ℎ 𝑣 along the central axis at different times for four massflow rates. ℎ 𝑣 is a function of the air density, viscosity, velocity, heat capacity, and thermal conductivity. Fig. 14 demonstrates that as the transient temperature increases, the heat transfer coefficient gradually increases as well. . Within the range of temperatures simulated (300-630 K), the calculated 𝐵𝑖 values range from 0.096 to 0.11. Given that 𝐵𝑖 ≤ 0.1, it is justifiable to neglect the intra-particle temperature gradient. Consequently, in this study, we can reasonably assume that internal bead resistance is negligible. It should be noted that even in the presence of a thermal gradient inside the solids, Eq.7 remains applicable. However, in that situation, ℎ 𝑣 also considers the conduction resistance of the solid.

Comparing the heat transfer coefficient and the effective gas thermal conductivity

Heat and mass transfer in packed beds with randomized packing of particles at high temperatures is a complex phenomenon. In this study, a new correlation is developed based on the numerical inverse analysis method, taking the glass bead diameter 𝑑 𝑝𝑎𝑟 as the characteristic length. The dimensionless heat transfer correlation is given by:

𝑁𝑢 = 2 + 1.54𝑅𝑒(𝑇 ) 0.6 𝑃 𝑟(𝑇 ) 1∕3 (16)
where the value of 𝑓 (1.54) is obtained as the average of the four 𝑓 values listed in Table 7. 𝑅𝑒(𝑇 ) and 𝑃 𝑟(𝑇 ) indicate that the Reynolds number (𝑅𝑒) and Prandtl number (𝑃 𝑟) are functions of temperature. Using the coefficient 𝑓 and its impact on the relative error (objective function 𝑆) as an example, as shown in Fig . 9 (b), we observed that employing the new coefficients (𝑓 = 1.54) resulted in 𝑆 = 0.0093. In contrast, when using the coefficients from the Wakao model (𝑓 = 1.1), we obtained 𝑆 = 0.01.

In Eq. 14, 𝑆 is defined as the average root mean squared relative error between measured and simulated temperatures at ten specific locations.

Fig. 16 (a) presents a comparison between heat transfer correlations reported in the work of Wakao [START_REF] Wakao | Effect of fluid dispersion coefficients on particle-to-fluid heat transfer coefficients in packed beds: Correlation of nusselt numbers[END_REF] and the results obtained from the inverse analysis in this study. The correlation proposed by Wakao et al is based on a large amount of experimental data. It is obvious that the Nusselt numbers increase with increasing Reynolds numbers. As presented in Fig. 16 (a), for lower Reynolds number values (𝑅𝑒 = 58 ∼ 252), the data from Wilke and Hougen [START_REF] Wilke | Mass transfer in the flow of gases through granular solids extended to low modified reynolds numbers[END_REF], Malling and Thodos [START_REF] Malling | Analogy between mass and heat transfer in beds of spheres: Contributions due to end effects[END_REF], and De Acetis and Thodos [START_REF] Acetis | Mass and heat transfer in flow of gases through spherical packings[END_REF] are slightly above the solid blue line (𝑁𝑢 = 2+1.1 𝑅𝑒 0.6 𝑃 𝑟 1∕3 ). The reason for the Galloway [START_REF] Galloway | Effect of configuration on mass and heat transfer in beds stacked spheres[END_REF] data being lower than predicted for all correlations is due to unique aspects of his study. Galloway investigated water evaporation into air from 17 mm porous spheres, arranged differently, and with solid particles maintaining constant surface temperatures. The newly proposed correlation (𝑁𝑢 = 2+1.54 𝑅𝑒(𝑇 ) 0.6 𝑃 𝑟(𝑇 ) 1∕3 ) demonstrates better suitability for the low Reynolds number range (𝑅𝑒 = 58 ∼ 252). It should be noted that the Reynolds number (𝑅𝑒) varies within the packed bed during each test, and thus the red dots in Fig. 16 (a) represent the average values of 𝑅𝑒.

The axial gas effective thermal conductivity, 𝑘 𝑔,𝑒𝑓 𝑓 ,|| , is expressed by Eq.13, which is a function of both Reynolds number (𝑅𝑒) and Prandtl number (𝑃 𝑟). Table 7 demonstrates that the value of the factor 𝑐 2 is not constant but increases with increasing 𝑅𝑒. To account for the dispersion effect resulting from the gas flow, a new correlation is developed based on the inverse analysis results obtained in this study. A previous work by DeGroot et al. [START_REF] Degroot | Closure of non-equilibrium volume-averaged energy equations in high-conductivity porous media[END_REF] suggests that the commonly used characterization of dispersion based solely on the Peclet number (𝑃 𝑒 = 𝑅𝑒 ⋅ 𝑃 𝑟) is not entirely accurate. In this study, we adopted a proposed function of the form 𝑛 1 𝑅𝑒 𝑛 2 𝑃 𝑟 to account for the dispersion effect caused by the gas flow. A least-squares analysis reveals that a suitable function for the dimensionless axial dispersion conductivity can be expressed as: 𝑘 𝑑𝑖𝑠,|| ∕𝑘 𝑔 = 0.00053𝑅𝑒(𝑇 ) 2.21 𝑃 𝑟(𝑇 )

Fig. 16 (b) presents a comparison between the dimensionless axial dispersion conductivity 𝑘 𝑑𝑖𝑠,|| ∕𝑘 𝑔 obtained in this study and the results of DeGroot et al. [START_REF] Degroot | Closure of non-equilibrium volume-averaged energy equations in high-conductivity porous media[END_REF]. The error bars indicate that the maximum difference between the two is within 20%. It is observed that the dispersion gas thermal conductivity increases with the Peclet number, as expected from the theory presented in the introduction. In this study, the axial dispersion gas thermal conductivity is about 5.9 times and 67.1 times higher than the gas thermal conductivity when the Peclet numbers are around 55 and 165, respectively.

Conclusions

The objective of this work was to develop and validate a generic two-equation model coupled with mass and momentum conservation for high-temperature TES and determine the effective conductivities and the heat transfer coefficients needed to inform the two-equation model. Packed beds undergoes a transient heat storage process, resulting in a broad range of temperatures for high-temperature systems. As a consequence, the thermophysical properties of both the fluid and solid exhibit significant variations. In this work, we developed a novel packed bed TES device, based on the transient single blow technique, to investigate the effective thermal conductivity of the gas and solid phases, and the heat transfer coefficient under high-temperature operating conditions. Experiments were conducted by passing hot air through cold particles at four different mass flow rates, within a Reynolds number range of 58 to 252. Unknown parameters are simultaneously determined using numerical inverse analysis, consisting in detailed simulations and state-of-theart multi-objective optimization. Flow inside packed beds is considered compressible and laminar; its thermal expansion is directly implemented in the mass equation by using the perfect gas law. The velocity field is modeled by the Forchheimer equation based on the Reynolds number. Heat transfer equations are described using the two-equation model in both phases inside packed beds. The variable gas and solid thermophysical properties are used over a wider range of temperatures (293K to 630K). The numerical model was implemented in finite volumes in PATO. The optimization part, which included sensitivity analysis using LHS method and the optimization process using a local optimization method (NL2SOL), was performed using the open-source optimization software Dakota. The results indicate that the value of heat transfer coefficient ℎ 𝑣 in the two-equation model falls in the range of 1.0 × 10 4 ∼ 1.93 × 10 4 W/(m 3 ⋅ K) under the given conditions. It was observed that the gas and solid temperatures exhibit local thermal non equilibrium within the packed bed. A new correlation has been proposed based on the Wakao correlation [START_REF] Wakao | Effect of fluid dispersion coefficients on particle-to-fluid heat transfer coefficients in packed beds: Correlation of nusselt numbers[END_REF] as follows: 𝑁𝑢 = 2 + 1.54𝑅𝑒(𝑇 ) 0.6 𝑃 𝑟(𝑇 ) 1∕3 . This correlation that accounts for the first time for the temperature variations of the dimensionless numbers was validated in the Reynolds number range of 𝑅𝑒 = 58 to 252. It was observed in this study that the dispersion gas thermal conductivity increased with increasing Peclet number. Specifically, the axial dispersion gas thermal conductivity was found to be around 5.9 and 67.1 times higher than the gas thermal conductivity at Peclet numbers of around 55 and 165, respectively. Results for the axial dispersion gas thermal conductivity have been fitted to a new correlation: 𝑘 𝑑𝑖𝑠,|| ∕𝑘 𝑔 = 0.00053𝑅𝑒(𝑇 ) 2.21 𝑃 𝑟(𝑇 ). The detailed high-temperature heat and mass transfer model proposed in this work was therefore validated for a large range of conditions, and new correlations were obtained to account for high-temperature effects for glass beads and air. This open perspectives towards more accurate sizing and monitoring models for high-temperature TES systems.

  (a) experimental setup (b) photograph of the experimental setup

Figure 1 :

 1 Figure 1: Schematic diagram and photograph of the experimental system and enlarged view of the test section showing the thermocouple identification numbers.

Figure 2 :

 2 Figure 2: A schematic of the thermocouple position.

Figure 3 :

 3 Figure 3: Two-dimensional axisymmetric numerical model of the test section.

  region and the energy equation in the tube region write,

  𝑇 𝑔 = 𝑇 𝑡 , 𝑘 𝑔,bed 𝑇 𝑠 = 𝑇 0 𝜕 𝑤 𝑝 = -( 𝜇 𝑔 𝐊⋅𝒏⋅𝒏 + 𝜌 𝑔 𝛽 ⋅ 𝒏 ⋅ 𝒏 |𝜖𝒗 𝒈 |

Figure 4 :

 4 Figure 4: Wall surface temperatures input in the simulation (𝑞𝑚=5.97×10 -3 kg/s).

Figure 5 :

 5 Figure 5: PATO multi-block mesh. (a) Zoom at the interface. (b) Zoom at the packed bed region (c) Side view of wedge-mesh.

Figure 6 :

 6 Figure 6: The steps of sensitivity analysis and optimization process.

Figure 7 :

 7 Figure 7: Comparison of predicted and measured temperatures with different mass-flow rates.

3 Figure 8 :

 38 Figure 8: Sensitivity analysis of the error 𝑆 to input parameters (𝑓 , 𝑐 1 , 𝑐 2 , 𝑐 3 ), 𝑞𝑚= 1.02 × 10 -2 kg/s.

Figure 9 :

 9 Figure 9: The influence of three variables (𝑐 1 , 𝑐 2 , 𝑓 )on the error 𝑆 .

Fig. 11

 11 Fig.11 presents the pressure 𝑝 and Darcy velocity magnitude 𝑢 fields within the packed bed region at 4000 s.The pressure at the inlet is 101358 Pa, while at the outlet, it is 101325 Pa. The significant pressure gradient observed in the packed bed region can be attributed to the inclusion of the Forchheimer equation in the model. The Darcy velocity magnitude 𝑢 is 0.22 m/s at the inlet and 0.1 m/s at the outlet. The velocity is updated based on the pressure and gas temperature. Inside the packed bed, the gas velocity ⟨𝐯 𝐠 ⟩ 𝐠 is determined by dividing the Darcy velocity 𝐮 by the porosity 𝜀 𝑔 .

Figure 10 :

 10 Figure 10: Dynamic viscosity and air density distribution predicted with the two-equation model.

  (a) pressure 𝑝 (Pa) (b) velocity magnitude 𝑢 (m/s)

Figure 11 :

 11 Figure 11: Pressure and velocity magnitude predicted with the two-equation model.

  (a) gas temperature 𝑇 𝑔 (K) (b) absolute temperature difference between gas and solid |𝑇 𝑔 -𝑇 𝑠 | (K)

Figure 12 :

 12 Figure 12: Gas temperature and absolute temperature difference in the porous materials region predicted with the twoequation model.

Figure 13 :

 13 Figure 13: Volumetric heat transfer coefficient ℎ 𝑣 (W/m 3 ⋅ K) distribution predicted with the two-equation model.

Figure 14 :

 14 Figure 14: The distribution of ℎ 𝑣 along the central axis with different mass flow rates.

Figure 15 :

 15 Figure 15: Biot number 𝐵𝑖 distribution predicted with the two-equation model.

  (a) Nusselt number vs. Reynolds number (b) dispersion gas thermal conductivity vs. Peclet number

Figure 16 :

 16 Figure 16: Comparison between the experimental data and data obtained using the numerical inverse analysis.

Fig. 15 presents

 15 Fig.15presents Biot number 𝐵𝑖 fields within the packed bed region at 4000 s. The 𝐵𝑖 was determined as follows: 𝐵𝑖 = ℎ 𝑠 𝐿 𝑐 ∕𝑘 𝑠 , where ℎ 𝑠 = ℎ 𝑣 ⋅𝑑 𝑝𝑎𝑟 ∕6(1-𝜀 𝑔 ), and 𝐿 𝑐 = 𝑉 𝑠 ∕𝐴 𝑠 = 𝑑 𝑝𝑎𝑟 ∕6[START_REF] Esence | A review on experience feedback and numerical modeling of packed-bed thermal energy storage systems[END_REF]. Within the range of temperatures simulated (300-630 K), the calculated 𝐵𝑖 values range from 0.096 to 0.11. Given that 𝐵𝑖 ≤ 0.1, it is justifiable to neglect the intra-particle temperature gradient. Consequently, in this study, we can reasonably assume that internal bead resistance is negligible. It should be noted that even in the presence of a thermal gradient inside the solids, Eq.7 remains applicable. However, in that situation, ℎ 𝑣 also considers the conduction resistance of the solid.

  improved Schumann's one-dimensional model by including heat capacity for the fluid and thermal conductivities and adding thermal dispersion in the 𝑁𝑢 Nusselt number based on ℎ 𝑠 𝐧 𝐠𝐬 outwardly directed unit normal vector pointing from the gas phase toward the solid phase, 𝐧 𝐠𝐬 =-𝐧 𝐬𝐠 𝑝 gas pressure, kg m -1 s -2

	Nomenclature
	Greek symbols
	𝛼	thermal diffusivity, 𝑘∕𝜌𝑐 𝑝 , m 2 s -1	⟨𝑝⟩ 𝑔	Intrinsic average pressure, kg m -1 s -2
	𝛽	Forchheimer coefficient, m -1	𝑃 𝑟	Prandtl number
	𝜀 𝑖	volume fraction of the 𝑖-phase	𝑞𝑚	gas mass-flow rate, kg s -1
	𝜇 𝑔	dynamic viscosity of gas, kg m -1 s -1	𝑅	gas constant, J K -1 mol -1
	𝜌 𝑖	density of the 𝑖-phase, kg m -3	𝑅 𝑏	tube radius, m
	Latin symbols	𝑅 𝑤	tube radius with wall thickness, m
	𝑎	specific gas-solid surface area, m -1	𝑅𝑒	Reynolds number based on the particle diameter
	𝐴 𝑠𝑔	area of the s-g interface contained in the averaging volume, 𝑉 ,	𝑆	the relative error
		m 2	𝑠 𝑖	scalar field used in LTNE model
	𝐛 𝐷 𝑏 𝑑 𝑝𝑎𝑟 𝐷𝑎 𝑓 ℎ 𝑠 ℎ 𝑣	tube diameter, m particle diameter, m Darcy number coefficients associated with Nusselt number interstitial heat transfer coefficient, W m -2 K -1 volumetric heat transfer coefficient, W m -3 K -1	𝑇 𝑎𝑣 ⟨𝑇 𝐗	average air temperature, K Forchheimer tensor
	𝐈	Identity tensor	Subscripts
	𝐤 𝐝𝐢𝐬	dispersion term of the thermal conductivity, W m -1 K -1	∥	horizontal flow directions
	𝐤 𝐢,𝐞𝐟 𝐟	effective thermal conductivity tensor of the 𝑖-phase, W m -1	⟂	transverse flow directions
		K -1	atm	atmosphere
	𝐤 𝐢,𝐬𝐭𝐚	static term of the thermal conductivity of the 𝑖 phase, W m -1 K -1	con	conductivity
			dis	dispersion
	𝑘 𝑖	thermal conductivity of the 𝑖-phase, W m -1 K -1	g	gas
	𝑘 𝑡	tube thermal conductivity, W m -1 K -1	in	inlet
	𝐤 𝐢,𝐜𝐨𝐧	conductivity term of the thermal conductivity of the 𝑖 phase, W m -1 K -1	inf	ambinent
			num	numerical
			s	solid
			sta	static
			t	tube

𝐢𝐢 vector field that maps ∇ ⟨𝑇 𝑖 ⟩ 𝑖 onto T𝑖 in LTNE model 𝐵𝑖 Biot number 𝑐 1 ,𝑐 2 , 𝑐 3 coefficients associated with thermal conductivity 𝑐 𝑝,𝑖 heat capacity of the 𝑖-phase, J kg -1 K -1 𝐤 𝐢,𝐭𝐨𝐫 tortuosity term of the thermal conductivity of the 𝑖 phase, W m -1 K -1 𝐊 permeability tensor, m 2 𝐿 tube length, m 𝑀 gas molar mass, kg mol -1 𝑖 ⟩ 𝑖 initrinsic phase average temperature in the 𝑖-phase, K 𝑇 𝑔 gas temperature, K 𝐮 Darcy velocity, 𝜀 𝑔 ⟨𝐯 𝐠 ⟩ 𝐠 , m s -1 𝑢 magnitude of the Darcy velocity, m s -1 𝐯 𝐠 velocity of the gas phase, m s -1 ⟨𝐯 𝐠 ⟩ 𝐠 intrinsic phase average velocity of the gas, m s -1 ⟨𝐯 𝐠 ⟩ superficial average velocity, m s -1 ṽ𝐠 the deviation of gas velocity, ṽ𝐠 = 𝐯 𝐠 -⟨ 𝐯 𝐠 ⟩ 𝐠 , m s -1

] 𝐤 𝐠,𝐞𝐟 𝐟 = 𝜀 𝑔 𝑘 𝑔 𝐈 + 𝑘 𝑔 𝑉 ∫ 𝐴 𝑔𝑠 𝐧 𝑔𝑠 𝐛 𝐠𝐠 𝑑𝐴 -𝜌 𝑔 𝑐 𝑝,𝑔 ⟨ ṽ𝑔 𝐛 𝐠𝐠 ⟩ = 𝐤 𝐠,𝐜𝐨𝐧 + 𝐤 𝐠,𝐭𝐨𝐫 + 𝐤 𝐝𝐢𝐬 (3) 𝐤 𝐬,𝐞𝐟 𝐟 = 𝜀 𝑠 𝑘 𝑠 𝐈 + 𝑘 𝑠 𝑉 ∫ 𝐴 𝑔𝑠 𝐧 𝑠𝑔 𝐛 𝐬𝐬 𝑑𝐴 = 𝐤 𝐬,𝐜𝐨𝐧 + 𝐤 𝐬,𝐭𝐨𝐫

Table 1

 1 Empirical correlations for the heat transfer coefficient in packed beds.

	Year	Investigators	Correlation		Range	Remarks
	1952 W. Ranz [28]	𝑁𝑢 = 2 + 1.8𝑅𝑒 0.5 𝑃 𝑟 1∕3		10 < 𝑅𝑒∕𝜀 𝑔 < 1000	experimental method, steady state
					not mentioned for 𝑇	constant properties
	1970 Galloway et al. [29] 𝑁𝑢 = 2 + 1.354𝑅𝑒 1∕2 𝑃 𝑟 1∕3 + 0.0326𝑅𝑒𝑃 𝑟 1∕2	𝑅𝑒∕𝜀 𝑔 < 5000	experimental method, steady state
					𝑇 =310 K	constant properties
	1979 Wakao et al. [17]	𝑁𝑢 = 2 + 1.1𝑅𝑒 0.6 𝑃 𝑟 1∕3		3 < 𝑅𝑒 < 10000	analysis of experimental data
					not mentioned for 𝑇	from the literature
	2012 Yang et al. [30]	𝑁𝑢 = 2.1 + 0.465𝑅𝑒 0.63 𝑃 𝑟 1∕3		100 < 𝑅𝑒∕𝜀 𝑔 < 5000 experimental method, transient
					298 K < 𝑇 < 358 K	constant properties
	2016 Naghash et al. [31] 𝑁𝑢 = (0.0012 ± 0.00273)𝑅𝑒 𝑑 ℎ	(1.647±0.501) 𝑃 𝑟 1∕3	40 < 𝑅𝑒 < 120	experimental method, steady state
			𝑅𝑒 𝑑ℎ = 𝑅𝑒∕1000(2𝜀 𝑔 ∕3𝜀 𝑠 )𝑑 𝑝𝑎𝑟		299 K < 𝑇 < 338 K	variable properties
	2022 Qu et al. [13]	𝑁𝑢 = 0.345∕𝜀 𝑔 ⋅ (2 + 1.033(𝑅𝑒∕𝜀 𝑔 ) 0.6 𝑃 𝑟 1∕3	200 < 𝑅𝑒∕𝜀 𝑔 < 1000 experimental method, transient
					283 K < 𝑇 < 345 K	constant properties
	2023 this work	𝑁𝑢 = 2 + 1.54𝑅𝑒(𝑇 ) 0.6 𝑃 𝑟(𝑇 ) 1∕3	58 < 𝑅𝑒 < 252	experimental method, transient
					293 K < 𝑇 < 630 K	variable properties

Table 2

 2 Dimension and structural properties of the packed bed and glass beds.

	Test section dimensions	Porosity Particle diameter	Particle density Permeability Forchheimer coefficient
	diameter 𝐷 𝑏 × length 𝐿 (mm)	𝜀 𝑔	𝑑 𝑝𝑎𝑟 (mm)	𝜌 𝑠 (kg m -3 )	K (m 2 )	𝛽 (m -1 )
	194 × 940	0.383	16	2500	2.52 × 10 -7	1201.18

Table 3

 3 The uncertainty analysis in the experimental measurements.

					Range	Uncertainty
	Temperature sensor (Pico Technology)	-270 • C to +1820 • C ± 0.2% of the reading
	K-type thermocouple	-100 • C to + 800 • C	± 0.5 %
	Thermo-anemometer	0.15 to 30 m/s	± 3 % of the reading
	Temperature measurement in the tube	293 K to 630 k	1.19 % of the reading
	(due to the position of the thermocouple)
	𝐾 =	𝑑 𝑝𝑎𝑟 150𝜀 𝑠 2 𝜀 𝑔 2 3	, 𝛽 =	1.75𝜀 𝑠 𝑑 𝑝𝑎𝑟 𝜀 𝑔 3	(6)

Table 6

 6 Experimental parameters at room temperature.

Table 7

 7 The value of 𝑓 , 𝑐 1 and 𝑐 2 obtained in the optimization process for different mass-flow rates.

	𝑞𝑚 (kg/s)	3.65 × 10 -3	5.97 × 10 -3	7.34 × 10 -3	1.02 × 10 -2
	𝑅𝑒	58∼100	100∼151	150∼190	223∼252
	𝑐 1	1.21 ± 0.03	1.22 ± 0.03	1.27 ± 0.03	1.29 ± 0.03
	𝑐 2	0.11 ± 0.005	0.19 ± 0.008	0.28 ± 0.02	0.42 ± 0.03
	𝑓	1.61 ± 0.09	1.61 ± 0.07	1.62 ± 0.07	1.31 ± 0.09
	𝑘 𝑠,𝑒𝑓 𝑓 (W∕(m ⋅ K))	1.1∼1.6	1.1∼1.5	1.1∼1.4	1.2∼1.4
	𝑘 𝑔,𝑒𝑓 𝑓 ,|| (W∕(m ⋅ K)) 0.132∼0.149	0.37∼0.39	1.14∼1.26	2.55∼2.64
	ℎ 𝑣 (W∕(m 3 ⋅K))	1.0 × 10 4 ∼1.3 × 10 4	1.16 × 10 4 ∼1.47 × 10 4	1.51 × 10 4 ∼1.80 × 10 4	1.68 × 10 4 ∼1.93 × 10 4
	𝑆	9.41 × 10 -3	9.84 × 10 -3	9.32 × 10 -3	9.14 × 10 -3

Acknowledgements

This work was supported by the China Scholarship Council (CSC) program (No. 201907040064).

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

A. Appendix. Validation of Ergun equation

To capture the pressure drop across the packed bed, a differential pressure transmitter from EMERSON FISHER ROSEMOUNT was used, which has an uncertainty of ±0.055% F.S. (Full Span is 0-620 mbar). The transmitter's electrical signals were relayed to a data acquisition system for real-time observation and subsequent analysis.

The pressure drop and superficial velocity were measured for the ten mass-flow rates (𝑞𝑚) considered in this work. Fig. 17 compares measured results (represented by dots) with calculated results using Ergun's equation (shown as solid lines). Ergun's equation is found to accurately predict the pressure drops across the packed beds, with an 𝑅 2 value of 0.998.

B. Appendix.Optimization results at different mass flow rates

In section 4.1.1, the determination of the effective thermal conductivity of gas and solid, as well as the heat transfer coefficient, was described using a case with a high mass-flow rate (𝑞𝑚=1.02×10 -2 kg/s). A sensitivity analysis was conducted, revealing that the coefficient 𝑐 3 , which is related to 𝑘 𝑔,𝑒𝑓 𝑓 ,⟂ , had minimal impact on the output results, even at high flow rates. Next, an optimization process was performed to obtain the values of 𝑐 1 , 𝑐 2 , and 𝑓 . In this appendix, we present the optimization results for three additional cases. The results are shown in Fig. 18, and the optimization process was conducted using the adaptive nonlinear least-squares algorithm (NL2SOL) in the Dakota software.

In the left side of Fig. 18, we can observe the effects of three variables, 𝑐 1 , 𝑐 2 , and 𝑓 , on the error 𝑆. On the right side, the impact of a single variable, 𝑓 , on the error 𝑆 is presented. The optimization process initially involves determining the values of 𝑐 1 and 𝑐 2 by including all


Pierre Carrive

Jfm Conseils

Exploiting cylinder focalisation of GPR signal for the detection of buried pipes and water leakages

Nowadays, the precise pipe localisation and water leakage detection have important societal impact such as avoiding explosions due to gas pipe degradation or large volume of water lost along its transportation. Ground-Penetrating Radar (GPR) is a common tool used for utility detection. A unique form of signal enhancement appears in GPR profiles recorded over spheres, ovoids and cylinders when presenting some specific dielectric contrasts with the surrounding medium. We study the potential applications of this peculiar phenomenon to detect buried pipes and water leakages. First we show field and laboratory measurements illustrating the signal enhancement effect. Then, numerical experiments explore the range of dielectric contrasts, depth and size of cylinders for which the effect appears, i) in homogeneous surrounding media, ii) for a cylinder embedded in a trench, iii) for a cylinder surrounded by a water infiltration bulb. In the homogeneous case, we show that the amplitude of the wave going around and through the cylinder can be 4.6 times higher that the one reflected on the top boundary of the cylinder. We investigate the interest of this effect to detect an anomalous water leakage along the pipe.

Introduction

Anomalous radar signal enhancement appears in groundpenetrating radar profiles recorded over electrically resistive ground. Examples include reflections from localized water infiltration bulbs in sands [START_REF] Saintenoy | Evaluating Ground Penetrating Radar use for water infiltration monitoring[END_REF][START_REF] Léger | Hydrodynamic parameters of a sandy soil determined by Ground-Penetrating Radar monitoring of Porchet infiltrations[END_REF] and sedimentary studies [START_REF] Arcone | Sedimentary architecture beneath lakes subjected to storms: Control by turbidity current bypass and turbidite armouring, interpreted from ground-penetrating radar images[END_REF]. Airborne radar reflected events from parts of the Greenland ice sheet resemble those of Jupiters icy Galilean satellites [START_REF] Rignot | Unusual radar echoes from the Greenland ice sheet[END_REF], both being intense with anomalous circular and linear polarization ratios. These phenomena have been attributed to refraction within bodies containing concentric incremental increases in optical index, although [START_REF] Hagfors | Refraction scattering as origin of the anomalous radar returns of Jupiter's satellites[END_REF][START_REF] Hagfors | Refraction scattering in the anomalous reflections from icy surfaces[END_REF] explain that there is usually an abrupt transition between solid ice and compacted snow that contradicts the existence of gradual transition. Another study from Legall et al. (2010) explains radar-bright channels on Titan as riverbeds with rounded pebbles acting as efficient natural retro-reflectors. In all cases the likely anomaly is cylindrical or spherical, such as a conduit, boulder or saturation bulb. In such cases, refractive focusing within the anomaly may be the cause, whereby surface waves along the curved surface continually refract along short cuts through the anomaly, which results in enhanced backscatter that emanates from near the anomaly base [START_REF] Rheinstein | Backscatter from spheres: A short pulse view, Antennas and Propagation[END_REF]. Although the propagation dynamics of this phenomenon is well understood in terms of Mie scattering and optical ray paths, the associated backscatter attributes of amplitude, dispersion and phase are not conceptually obvious from the mathematics. Here we aim to provide some understanding of these backscatter attributes with laboratory experiments and finite difference time domain modeling. Our laboratory experiment comprised GPR reflection profiles of a sphere with relatively high permittivity buried in a sand matrix of relative low permittivity, and with diameter comparable to the dominant wavelength of the broadband GPR pulse. Our numerical model simulates different experiments.

Laboratory experiments

Figure 1 : a) B-scan profile acquired across the 12-cm diameter glass-sphere with 2.6 GHz GSSI antennas. b) The A-scan acquired at a 2-m distance, above the sphere.

The experiment took place in a sandbox, 4 m long and 1.2 m wide. A metal sheet was placed at the bottom of the sand layer at approximately 1 m depth. A glass sphere of 12 cm diameter was buried 0.3 m deep (relative to the top of the sphere) in dry sand with a known relative dielectric permittivity of εm = 2.6. The glass relative dielectric permittivity value is reported to be εs = 6.9. A profile was acquired on the surface of the sand layer using 2.6 GHz GSSI antennas. The zero time reference was estimated from the start of the direct coupling in the raw data minus 0.13 ns for the 4 cm antenna separation. No range gain was applied. The four perfect hyperbolas indicate nearly-point scatterers. The two linear reflection horizons are from the metal sheet and from the box bottom edges, the latter of which occurs at about the correct time delay of 1.5 ns behind the bottom reflection (which is at 88 cm). Several weak events indicate that the sand is not perfectly homogeneous (also evidenced by the slightly uneven bottom reflection). The hyperbola (a) is assumed to be from the top of the sphere, (b) from its bottom, and (b') and (c) origins will be explained in the rest of this paper. Hyperbola (d) is slightly delayed from the metal sheet reflection and originated from waves that diffracted either around the sphere or went through the sphere. Note that there is almost no discontinuity in the metal sheet horizon especially directly under the sphere. The theoretical delay caused by signals diffracting around the sphere relative to the normal (no-sphere) sheet reflection is only 0.74 ns. Hyperbola (a) was perfectly migrated using ε = 2.6, and the Kirchhoff migration routine with a 127 trace window, which verifies εm = 2.6 for the sand medium. Its half-cycle amplitude polarity sequence is the inverse of the transmitted pulse as expected for a reflection from an interface between an higher permittivity material below (glass) and lower one above (dry porous sand)

Numerical experiments with a single cylinder

To better understand the contrast of amplitude of the different events resulting from a spherical anomaly, we carried out several numerical simulations with GprMax [START_REF] Giannopoulos | Modelling Ground Penetrating Radar by gprMax[END_REF]. A point source Ricker wavelet centered on 2.6 GHz, is placed in the same place as the receiver 0.3 m above the top of the cylinder. The cells are millimeter sized. Three 0.80 m by 0.80 m domain were considered. The first one consists of a homogeneous medium with εm = 2.6 in which a cylinder of relative dielectric permittivity εs = 6.9 of radius R = 0.06 m is embedded at the depth h = 0.30 m (Fig. 2). On the obtained trace (Fig. 3), we recognise reflections a and b by their arrival times and clearly the amplification of reflection b compared to reflection a.

Keeping the geometry of Fig. 2, with a relative permittivity value of the surrounding medium constant to εm=9, we run different simulations varying the relative permittivity inside the cylinder from 12 to 78 by step of 3. The traces obtained for each εs value is displayed Fig. 4 and the ratio between amplitudes of reflections Ab /Aa is on Fig. 5. On this figure we see that for εs=33, we obtain the highest amplitude ratio of 4.6. Reflection b, which propagates a longer way than reflection a has an amplitude 4.6 times higher than the one of reflection a. For all considered relative permittivity value the amplification is happening as the ratio Ab /Aa stays above 1. A second analysis consists in keeping geometry of Fig. 2, with relative permittivities εm=9 and εs=33 (highest amplification observed on Fig. 5) and increasing the radius of the cylinder from 0.05 m to 0.29 m by step of 0.02 m.

The resulting radargram is presented on Fig. 6 and the amplitude ratio (Fig. 7) shows that the highest amplitude amplification is obtained for a radius of 0.07 m.

Pipe in a trench simulation

One of the most used application of GPR is locating underground pipes. Most of these pipes are situated down at the bottom of a trench to set them in. We constructed a theoretical geometry in Fig. 7 where the pipe is a cylinder of radius 0.06 m with a relative permittivity εs=25, embedded at the depth h=0.3 m in a parallelepiped with εt=12, in a surrounding homogeneous medium with a dielectric constant εm=9. The trace obtained with antennas positioned on top middle of the pipe (Fig. 8) show again a reflection b with higher amplitude than reflection a.

Interferences with corners and bottom of the trench is happening but do not prevent the enhancement phenomena.

Figure 8 : Geometry of the second domain: a pipe in a trench. The yellow media is at εt=12, the white media is at εm=9 and the pipe is at εs=25.

Leaking water pipe simulation

One big challenge is to locate leaks in city's drinking water network. A third conceptual model (Fig. 9) consists of an homogeneous medium with εm=9, an outside cylinder representing a medium with high water content (εl=35, R=0.12 m) and a small cylinder representative of a plastic water pipe (εs=81, R=0.06). In the case of a plastic water pipe, reflection b is higher amplitude than reflection a. If we consider a surrounding zone simulating the medium with higher water content due to some water leak (bottom of Fig. 11), the reflection from the bottom of the pipe and its leak is even higher amplitude that the reflections on top of the infiltration bulb. It looks like the zone with higher water content is amplifying the enhancement effect. This effect should help to recognize water leak along a water pipe.

Conclusion

In this study we showed laboratory experiment exacerbating amplitude enhancement due to the geometry of the buried object. Because of its forms (cylindrical and spherical), the backscattered electromagnetic wave. We showed that we could retrieve the same behavior using numerical simulations, and found relationship between the cylinder radius, the permittivity contrasts, and the amplitude enhancement. In addition, we showed that this phenomena could be a help for detecting water leak along a plastic water pipe. In general, we want to point out the risk of miscalculating the depth of a pipe when enhancement conditions are meet. Indeed, the largest amplitude is usually expected to be related to the top of the pipe, but we show that it is actually related to a reflection coming from the bottom of the pipe.

Figure 2 :

 2 Figure 2 : Geometry of the numerical domain: a pipe is surrounded by an homogeneous medium

Figure 4 :

 4 Figure 4 : Traces obtained when varying εs from 12 to 78 using the geometry of Fig. 2 and εm=9.

Figure 5 :

 5 Figure 5 : Ratio of the amplitude of reflection b over the amplitude of reflection a when εs varies from 12 to 78 (red: maximum amplitudes, blue: absolute value of minimum amplitudes).

Figure 6 :

 6 Figure 6 : Traces obtained with the geometry of Fig. 2 when the cylinder radius varies from 0.05 to 0.29 m.

Figure 7 :

 7 Figure 7 : Ratio of the amplitude of reflection b over the amplitude of reflection a when the radius of the cylinder varies from 0.5 m to 0.29 m (red: maximum amplitudes, blue: absolute value of minimum amplitudes).

Figure 9 :

 9 Figure 9 : Numerical A-scan calculated for the geometry of the Figure 7.

Figure 10 :

 10 Figure 10: Geometry of the third domain: leaking water pipe