Von der Fakultät Energie-, Verfahrens- und Biotechnik der Universität Stuttgart zur Erlangung des akademischen Grades eines Doktors der Ingenieurwissenschaften (Dr.-Ing) genehmigte Abhandlung Vorgelegt von Philipp Bartsch aus Würzburg s aus Würzbslleeecurg Entwicklung neuer Modelle zur Beschreibung von Schüttgutströmungen in Wanderbettwärmeübertragern Hauptberichter: Prof. Dr. André Thess Mitberichter: Prof. Dr.-Ing. Prof. E.h. Peter Eberhard Tag der mündlichen Prüfung: 09/06/2020 Institut für Gebäudeenergetik, Thermotechnik und Energiespeicherung der Universität Stuttgart 2020 II Eingeständigkeitserklärung Eingeständigkeitserklärung Ich versichere, dass ich die vorliegende Arbeit selbständig verfasst und keine anderen als die angegebenen Quellen und Hilfsmittel benutzt habe; aus fremden Quellen ent- nommene Passagen und Gedanken sind als solche kenntlich gemacht. Name: _______________________________ Unterschrift: ________________________________ Datum: ________________________________ Kurzfassung III Kurzfassung Granulare Medien weisen viele Eigenschaften auf, die ihre Verwendung als Wärmeträ- ger- und Wärmespeicher in Solarthermischen Kraftwerken oder Industrieprozessen nahe legen. Wanderbettwärmeübertrager (WBWÜ) mit horizontalen Rohren bieten sich dabei an, um thermische Energie aus heißen Schüttgütern auszukoppeln. Ihre thermische Leis- tungsfähigkeit wird durch das komplexe, granulare Strömungsprofil im Apparat be- stimmt. Daher werden präzise Simulationsmodelle benötigt, um eine Grundlage für neuartige Entwurfswerkzeuge zu bilden. In dieser Arbeit wird ein Kontinuum-Modell entwickelt, welches die Berechnung des Strömungsfeldes und Wärmetransports in einem WBWÜ ermöglicht. Ergänzend werden Simulationen mit Hilfe eines partikel-diskreten Modells durchgeführt, welches Einblick in Strömungsphänomene gewährt, welche sich auf Ebene der Einzelpartikel abspielen. Ein Fokus der Modellierungsarbeiten ist die genaue Abbildung der leistungsbestimmen- den Phänomene. Beide Modelle sagen übereinstimmend die Bildung einer Stauzone auf dem Rohrschei- tel voraus. Unterhalb des Rohres beginnen die Partikel sich von der Rohroberfläche abzulösen und bilden eine Leerzone. Dieser Effekt wird nur durch das partikeldiskrete Modell, nicht aber durch das Kontinuum-Modell erfasst. Um die daraus resultierende Verringerung des Wärmeübergangs zu berücksichtigen wird das Kontinuum-Modell modifiziert. Die Simulationsergebnisse zum Strömungsfeld werden anhand von optischen Messun- gen an einem transparenten Wärmeübertrager-Mockup validiert. Die Simulationsergeb- nisse zum lokalen Wärmeübergang am Einzelrohr werden ebenfalls mit Messdaten ver- glichen, wobei die Massenstromdichte, die Korngröße und die Rohranordnung variiert werden. Der Einfluss der Partikelgröße und Massenstromdichte auf den Wärmeübergang wird vom Kontinuum-Modell gut erfasst. Die Abweichungen zwischen Modell und Experi- ment hinsichtlich des mittleren Wärmeübergangskoeffizienten am Rohr liegen bei we- niger als 20 %. Anhand des partikeldiskreten Modells wird die Abhängigkeit des granularen Strö- mungsfeldes von unterschiedlichen Einflussgrößen untersucht. Es zeigt sich, dass sich IV Abstract die Größe und Ausdehnung der Stauzone signifikant durch die Rohranordnung beein- flussen lässt. Gleiches gilt für die Position in der unteren Rohrhälfte, an der die Partikel beginnen, sich von der Rohroberfläche abzulösen. Die Ausdehnung der Stauzone hängt zusätzlich vom Oberflächenreibungskoeffizient der Rohrwand ab. Andere Größen wie die Massenstromdichte und die innere Reibung der Schüttung zeigen einen geringen Einfluss auf das Strömungsmuster. Anhand dieser Erkenntnisse wird eine vereinfachte Beschreibung des Wärmeübergangs zwischen einem horizontalen Rohr und einer Schüttung entwickelt. Das vereinfachte Modell wird den Messergebnissen des mittleren Wärmeübergangskoeffizienten am Rohr gegenübergestellt wobei Abweichungen von weniger als 13 % auftreten. Die Ergebnisse der Arbeit zeigen damit eine deutliche Verbesserung gegenüber beste- henden Modellen, wie dem von Baumann [1], Niegsch [2] oder Schlünder [3], die je- weils eine Abweichung von 40 bis 50 % bezogen auf die experimentellen Daten aufwei- sen. Diese Arbeit behandelt somit die im WBWÜ ablaufenden thermophysikalischen Vor- gänge auf unterschiedlichsten Detaillierungsgraden und Modellierungsebenen. Sie lie- fert damit einen wichtigen Beitrag zum kosteneffizienten Einsatz granularer Wärmeträ- germedien in solarthermischen Kraftwerken als Teil einer nachhaltigen Energieversor- gung. Teile dieser Arbeit wurden in referierten Fachzeitschriften veröffentlicht [4] [5]. Abstract Granular materials offer many advantages which qualify them for use as heat transfer medium and heat storage material in solar thermal power plants and industrial proces- ses. Moving bed heat exchangers (MBHE) with horizontal tubes are favorable to extract thermal energy from hot granular materials. Their thermal performance is determined by a complex granular flow field in the device. This calls for accurate simulation models on which novel design tools can be based. In this work a continuum model is developed to calculate the granular flow field and heat transport in a MBHE. The model is supplemented by simulations using a discrete- particle model which provides insight into flow phenomena occuring on the particle level. Abstract V Both models agree in predicting the formation of a stagnant area at the top of the tubes. Below the tubes, the particles start to separate from the tube surface and form a void area. This phenomenon is captured only by the discrete particle model but not by the continuum model. To account for the resulting decrease in heat transfer in this area, the continuum model is modified. The simulation results of the granular flow field are validated based on optical meas- urements in a transparent heat exchanger mockup. Furthermore, the local heat transfer coefficient at a single tube is examined in detail, varying the mass flow rate, grain size and the tube arrangement. The continuum model well captures the influence of the mass flow rate and the grain size on the heat transfer rate. The averaged heat transfer coefficient per tube predicted by the continuum model deviates by less than 20 % from the measured results. Using the discrete particle model, the dependence of the granular flow field on different influencing parameters is studied. It is determined that the size and extent of the stag- nant zone can be significantly influenced by the tube arrangement. The same goes for the position in the lower half of the tube where the particles start to detach themselves from the tube surface. In addition the extent of the stagnant zone depends on the surface friction of the tube wall. Other parameters, like the inner friction of the bulk and the mass flow rate show little influence on the flow pattern. Based on these findings, a simplified description of the heat transfer between a horizon- tal tube and a moving bulk is developed. A comparison of the results of the simplified model with the measured data of the averaged heat transfer coefficient per tube shows the deviation to be less than 13 %. The results of this work show significant enhancement compared to existing models such as the one of Baumann [1], Niegsch [2] or Schlünder [3], who exhibit deviations of 40 to 50 % from the experimental data. This work hence analyzes the thermo-physical phenomena of a MBHE on different lev- els of detail and length scale. Its results contribute to a sustainable future energy supply by promoting the cost-efficiency of solar thermal power plants with granular heat trans- fer fluids. Parts of this work have been presented in journal publications [4] [5]. VI Table of contents Table of contents Eingeständigkeitserklärung ...........................................................................................II Kurzfassung .................................................................................................................. III Abstract ......................................................................................................................... IV Table of contents ........................................................................................................... VI List of abbreviations ..................................................................................................... IX 1 Introduction .......................................................................................................... 1 2 State of the art ...................................................................................................... 4 2.1 ...... Modelling of granular flows ................................................................................... 4 2.1.1 ... Discrete particle models ......................................................................................... 4 2.1.2 ... Continuum models ................................................................................................. 6 2.2 ...... Prior work on heat transfer and granular flow in moving packed beds ................. 9 2.3 ...... Scope of investigation .......................................................................................... 13 3 Experimental setup ............................................................................................ 16 3.1 ...... Measurement of flow field ................................................................................... 17 3.2 ...... Measurement of local heat transfer coefficient .................................................... 21 4 Modelling............................................................................................................. 25 4.1 ...... Continuum model equations................................................................................. 28 4.1.1 ... Frictional stress model ......................................................................................... 30 4.1.2 ... Boundary conditions at the tube walls ................................................................. 32 4.2 ...... Discrete particle model......................................................................................... 35 4.2.1 ... Contact models ..................................................................................................... 36 4.2.2 ... Local averaging .................................................................................................... 38 4.3 ...... Model parameters ................................................................................................. 39 4.3.1 ... Parameters specific to the continuum model........................................................ 40 4.3.2 ... Parameters specific to the discrete particle model ............................................... 41 5 Validation ............................................................................................................ 43 5.1 ...... Discrete particle model......................................................................................... 43 5.1.1 ... Velocity profile at the tube surface ...................................................................... 47 5.1.2 ... Influence of the front wall friction on the measured flow field ........................... 49 5.2 ...... Continuum model ................................................................................................. 50 5.2.1 ... Flow field and packing fraction ........................................................................... 51 Table of contents VII 5.2.2 ... Modification of the continuum model .................................................................. 58 5.2.3 ... Heat transfer ......................................................................................................... 62 5.2.4 ... Summarizing conclusions on the model quality .................................................. 71 6 Simplified model of heat transfer at a single tube ........................................... 72 6.1 ...... Local extent of the characteristic flow sections ................................................... 72 6.1.1 ... Stagnant area ........................................................................................................ 74 6.1.2 ... Void area .............................................................................................................. 79 6.1.3 ... Summarizing conclusions about the extent of the characteristic flow sections ... 82 6.2 ...... The model of Niegsch .......................................................................................... 84 6.3 ...... Modifications to Niegsch’ model ......................................................................... 87 6.4 ...... Enhanced model ................................................................................................... 89 6.4.1 ... Evaluation of the effective thermal resistances .................................................... 89 6.4.2 ... Comparison with measurement data .................................................................... 92 7 Conclusions and outlook .................................................................................... 95 Appendix A: Selected details on experimental procedures ....................................... 98 A.1 Grain size distribution of quartz sand ....................................................................... 98 A.2 Error estimate of measurement of local heat transfer ............................................... 98 A.3 Error estimate of determination of internal granular flow speed ............................ 100 A.4 Measurement of the wall friction coefficient between quartz sand and an acrylic glass surface ........................................................................................... 102 Appendix B: Selected details on the continuum model ............................................ 104 B.1 Heat transfer through the void area below the tubes............................................... 104 B.2 Thermal contact resistance ...................................................................................... 105 B.3 Thermal conductivity in the continuum model ....................................................... 107 B.4 Interphase heat transfer coefficient between phases ............................................... 108 B.5 Plasticity model without “critical state approximation” ......................................... 109 B.6 Grid independence of simulation results ................................................................ 111 Appendix C: Simplified model - Analytical approximation of effective heat transfer coefficient............................................................................................ 113 C.1 Section I .................................................................................................................. 113 C.2 Section II ................................................................................................................. 115 C.3 Section III ............................................................................................................... 117 Appendix D: DEM model – sensitivity on different model parameters ................. 119 D.1 Sensitivity to variation of Young’s modulus .......................................................... 119 D.2 Sensitivity to variation of tube arrangement (supplement to section 6.1) .............. 120 D.3 Sensitivity to variation of inner friction angle (supplement to section 6.1) ........... 125 VIII Table of contents D.4 Sensitivity to variation of internal friction parameters ........................................... 126 D.5 Sensitivity to variation of the reference velocity .................................................... 128 References..................................................................................................................... 130 List of abbreviations IX List of abbreviations Abbreviations CAES Compressed air energy storage CFD Computational fluid dynamics CSP Concentrating solar power CV control volume D dimension DEM Discrete element method ED Event-driven method EOS Equation of State Exp experiment, experimental KTGF Kinetic theory of granular fluids MBHE Moving bed heat exchanger MD Molecular dynamics NSCD/CD (Non-smooth) contact dynamics method PE Polyethylene PIV Particle image velocimetry PT Particle tracking (velocimetry) px Pixel Sim simulation TA Tube arrangement TC Thermopcouple TES Thermal energy storage Symbols a ’shape parameter’ of stagnant zone in equation (40) A area b width/breadth C12 radiation coefficient CD Coefficient in equation (10) cp specific heat Cs Stefan-Boltzmann constant d particle diameter D rate of deformation tensor X List of abbreviations DT outer tube diameter e coefficient of restitution E Youngs modulus F, F force vector, magnitude of force g acceleration due to gravity Gr Grashof number h specific enthalpy H height hR (total) surface roughness (see equation (20)) I electric current, unit tensor Is inertial number J inertia troque k contact stiffness Kn Knudsen number l mean free path of gas molecules, length (see equation (67) and equation (60)) m mass M torque, moment ṁ mass flow rate Mf molecular mass of fluid molecules n unit normal vector N number Nu Nusselt Number p pressure P Power pp particle-particle Q particulate quantity averaged per control volume (see equation (34)) q̇ heat flux, heat transfer rate R ohmic/thermal resistance, radius r,s exponents in equation (14) Ra Rayleigh number Re Reynolds number Rep Reynolds Number per particle diameter Rg universal gas constant s spacing between tubes, overlap between particles (see section 4.2.1) S deviatoric part of the rate of deformation tensor List of abbreviations XI t residence time, time step ‘t’ (superscript) T Temperature T0 ’reduced’ temperature at wall due to thermal contact resistance (see equation (17)) u, u velocity vector, magnitude of velocity vector V volume α heat transfer coefficient β interphase drag coefficient γ accomodation coefficient (see equation (75) and (76)) δ gap width between wall and contacting particle surface Δ Delta (difference) ε volume fraction ε1, ε2 emission coefficients of wall and bulk εmax closest random packing (maximum packing fraction in frictional regime) εmin loosest random packing (minimum packing fraction in frictional regime) λ thermal conductivity μ friction coefficient, shear viscosity ν kinematic viscosity ρ density σ stress tensor τ, τ shear stress tensor, shear stress (scalar) φ fraction of (tube) surface covered by particles (see equation (18)) ϕi internal friction angle ϕr angle of repose ω angular position (at tube surface) Indices abs absolute deviation/error av average c contact con conduction d dynamic e elastic el electric fr frictional regime XII List of abbreviations front front side of simulation geometry (see Figure 11) fs interfacial exchange quantity between phases (equation (12)) g gas phase H horizontal i internal quantity of the bulk i,j substitutes for gas and solid phase l, ul loading,unloading case lim limiting value loc local m mean min minimum n normal direction p particle rad radiativ, radial rear rear side of simulation geometry (see Figure 11) ref reference quantity rel relativ roll rolling s solid phase, static S0 bulk property sep ’separate’ – point at the tube surface where particles start to detach sl ’sliding’ SZ stagnant zone t tangential V Vertical VZ void zone w weighting factor (see equation (34)) W (tube) wall WP ’wall-particle’ 1 1 Introduction 1 Introduction The steadily increasing demand for electrical and thermal energy by the world popula- tion and the associated consequences for our environment and climate are one of the greatest challenges of our time. Although the energy demand is rising more slowly than in the past (2017), it is expected to expand by 30 % between now and 2040 [6]. On the other hand there is consensus both in politics and science that the rise of the average global temperature must be limited to 2 °C (“two-degree target”) to prevent irreversible damage to the environment and society [7]. The energy sector accounts for the major part of CO2-emissions and thus decarboniza- tion of this sector is a vital requirement to meet the above target. To this end, renewable energy technologies need to be further developed. However, the fluctuating availability of wind and sun constitutes a challenge to a flexible and adjustable power supply. Therefore, these technologies need to be combined with energy storage systems to en- sure a secure electricity supply from these renewable sources. Appropriate energy stor- age systems need to offer high capacities and high powers, long storage times at small losses and short reaction times. In this situation, thermal energy storage (TES) systems can be a key technology and a central element in several emerging applications: They are the main component in pumped heat electricity storage (PHES) systems. They can also be used to increase the efficiency and the flexibility of other storage technologies (e.g. CAES systems [8]), of industrial processes or of power plants. Common examples can be found in steel indus- tries, where high temperature heat storage systems recuperate excess heat from blast furnaces thus reducing the overall energy consumption. With respect to power genera- tion, TES systems are considered for increasing the flexibility of conventional thermal power plants [9]. TES systems are an integral part of concentrating solar power (CSP) plants, ensuring a demand-oriented electricity production. This type of power plant converts the solar irra- diance into thermal energy at a high temperature level which facilitates the integration of heat storage. Commonly used TES-concepts for CSP plants are based on molten salts or solid materi- als [10]. Solid materials allow high operating temperatures but need a secondary heat 2 1 Introduction transfer fluid to charge/discharge thermal energy to/from the storage. Molten salts, on the other hand, can be used both as heat transfer fluid and storage material. However, currently available and cost-effective materials like the commonly used “solar salt” are chemically stable only up to about 550 °C and exhibit a melting point around 200 °C [11]. A further promising alternative is to use flowable, fine grained solid materials which can serve both as a heat transfer fluid and a storage material in the power plant (see Fig- ure 1, left). Such concepts were already proposed in the 1980s [12] [13] [14]. Thus, high operating temperatures are feasible using cost-effective materials, such as natural stones like quartz sand or sintered bauxite particles [15], while no (gaseous) secondary fluid is needed. Besides the receiver, a key component of such a CSP plant is a particle heat exchanger which is used to discharge the thermal energy from the bulk and transfer it to the work- ing fluid of the power cycle. To this end, different alternatives like rotary kilns, fluid- ized beds and moving beds have been assessed [1]. Moving bed heat exchangers (MBHE, Figure 1, right) are considered a most promising option due to their low para- sitic loads, compact design, low investment cost and little need for maintenance and adjustment controls [16]. In a MBHE the granular material slowly moves in a vertical direction, driven by gravity. Immersed in the moving bed are heat-transferring surfaces such as tubes or plates. The design of a robust and efficient MBHE relies on detailed knowledge of the fluid- and thermodynamic processes inside the unit. In particular, the effects of the particle flow on the device’s thermal performance are determining. Therefore the aim of the present work is to provide insight into the governing phenomena and to constitute a sound basis for the design and optimization of MBHE. 1 Introduction 3 Figure 1. Left: Schematic of solar thermal power plant with granular heat transfer fluid and particle steam generator (HTS = High Temperature Storage, LTS = Low Tempera- ture Storage) [1]. Right: Schematic of moving bed heat exchanger. 4 2 State of the art 2 State of the art Granular materials display a variety of intriguing features. Most notably, they can ex- hibit different states of matter, i.e. they can behave like a liquid, a gas, or like a rigid body depending on the way they are handled [17]. This behavior as well as other com- plex phenomena leads to challenges regarding the description and modelling of granular materials. Furthermore, granular materials usually form a two-phase system. The inter- stitial fluid between the particles may have considerable impact on the motion, the heat transport and the thermodynamic properties of the bulk. The following section gives an overview of modelling approaches for granular materials in general. Subsequent to this, prior work specifically related to the heat transfer and flow in moving packed beds is reviewed. 2.1 Modelling of granular flows Models for the description of granular materials fall into two general groups: Discrete particle models and continuum models. While the former treat the material as collection of distinct particles, continuum models treat it as a continuous medium. Due to high computational effort, discrete particle models are restricted to relatively small geometries and particle numbers. However, they require no macroscopic input parameters such as the angle of repose. Instead, these macroscopic parameters can be extracted from the simulation as a result of the particle properties like friction and shape [18] [19]. Therefore, discrete particle simulations are a powerful tool to gain insight in the “microscopic” processes on the grain level and can be well combined with continu- um models which are applicable to engineering scale geometries but require further in- formation about the macroscopic behavior of the bulk material. The following sections give a brief overview of the two types of models with a focus on their application in moving bed heat exchangers. 2.1.1 Discrete particle models Cundall and Strack [20] were the first who introduced a discrete particle approach to model the movement and interaction forces of rock masses and granular materials. Their 2 State of the art 5 approach became known as the discrete element method (DEM) and was very similar in spirit to the so-called molecular dynamics (MD) approach which had been introduced earlier (e.g. [21]). The DEM differs from MD in that the particles are no point masses but possess a finite volume. Both methods evaluate the trajectories of the particles based on the Newtonian equations of motion. The interaction forces are calculated based on a small (virtual) overlap between particles which is interpreted as an elastic deformation of the particles. The method therefore is also called ‘soft particle’ method. To translate the overlap into interaction forces, sever- al contact force models were proposed such as continuous potential models (e.g. [22]), linear (e.g. [23]) and non-linear (e.g. [24]) viscoelastic models and hysteretic models (e.g. [25]) for the normal contact forces and linear and non-linear tangential force mod- els (see e.g. [26] for an overview of available models). The geometrical relations of two interacting particles according to the DEM method are illustrated in Figure 1. The evaluation of the contact forces requires a very fine time step to resolve the small time and length scales involved in the contact interactions. Figure 2. Schematic of DEM method: Geometrical relations between two interacting particles [27]. The DEM has become very popular in the last decades due to the rapid increase of available computational power. It has been used to investigate the behavior of particu- 6 2 State of the art late solids in different fields such as agricultural [28], chemical, pharmaceutical [29] [30] and heavy industries [31]. It has also been applied in thermal energy storage tech- nology to investigate thermally induced stresses in packed-bed heat storage systems [32]. Besides the soft sphere methods, also hard sphere methods have been developed treating the particles as perfectly rigid bodies. Representatives of this class of models are the event-driven method (ED) [33] [34] which is primarily applicable to dilute and agitated “granular gases” (e.g. [35] [36]) and the (non-smooth) contact dynamics method (CD or NSCD) [37] [38] [39] which was developed for dense particle assemblies with multiple contact partners (e.g. [40] [41]). Compared to DEM (MD), the CD method offers a much larger time step as the very small length and time scales of the (viscoelastic) particle interactions are not resolved. Instead, the contact forces between particles are determined such that no mutual inter- penetration of the particles occurs and that the contact laws (i.e. frictional stick-slip con- straints) are satisfied [42]. In a network of particles the set of contact forces which con- sistently satisfies all constraints usually has to be determined iteratively [43]. Hence, the computational effort per time step and particle is higher for CD than for DEM. CD is advantageous for quasi-static systems of rigid particles, especially for limited particle numbers. For large systems with finite rigidity of the particles, DEM is favora- ble [43]. All discrete particle methods may also incorporate heat transport inside the granular matter. However, in this case the interstitial gas between the particles usually has to be taken into account even if the motion of the fluid-phase is not explicitly modelled. For densely packed systems heat transport is commonly modelled by a network of thermal resistances between particles. Different approaches have been proposed for these sys- tems [44] [45] [46] [47]. 2.1.2 Continuum models A continuum model relies on the assumption of a continuous medium with locally aver- aged properties like density or velocity. The flow field and the density distribution of the continuum are obtained by solving the balance equations for mass and momentum. However, to close the system of equations, constitutive relations for the stresses in the 2 State of the art 7 continuum are required which are provided by rheological models of the granular mate- rial. Rheological models of granular flows fall into two groups depending on the flow re- gime that is investigated: Models for slow, dense granular flows and models for rapid, strongly agitated flows [17] [48]. In case of strongly agitated flows, such as fluidized beds, particles interact by binary collisions. Models for such rapidly moving granular flows have been developed by analogy with the kinetic theory of gases (“Kinetic theory of granular fluids” – KTGF) [17]. In contrast to molecules, though, collisions between particles are highly inelastic, resulting in a strongly dissipative nature of “granular gases”. The KTGF has been used to investigate various types of the granular flow such as vi- brated systems and chutes [49] and fluidized beds [35] [50]. The kinetic theory is based on the assumption of binary collisions between particles. However, if the packing fraction, εs, increases (εs > 0.4), as is the case in moving packed beds, particles interact by enduring, frictional contacts with multiple neighbors and hence the kinetic theory is no more valid. The behavior of such slow, dense granular flows is traditionally subject to soil mechan- ics and the so-called critical state theory [51] [52]. Central phenomena occurring in this type of flows is the transition from a solid-like to a yielding state as well as the dilation and the compaction of the material under shear deformation. Such (elasto-)plastic deformations of a granular material are modelled based on plastici- ty theory ( [48] [17]). These models, which are mostly based on metal plasticity ( [53] in [54]), rely on the concept of a pressure-dependent yield function [54]. The yield func- tion defines whether the material remains rigid or whether it flows under a given state of internal stresses. A simple yield function is the Drucker-Prager yield function [55] which is based on the Coulomb yield criterion according to which the shear stress at failure is proportional to the normal stress and a coefficient of internal friction. In three dimensions the Drucker-Prager yield surface (yield function) has the shape of cone (see Figure 3). 8 2 State of the art Figure 3. Drucker-Prager yield surface (without cohesion) in 3-D pricipal stress (σ1, σ2, σ3) space. For stress states inside the cone, the material behaves rigid/elastic. For stress states on the surface of the cone, yielding occurs. Stress states outside the cone are not allowed. Once yielding occurs, a flow rule defines how the plastic deformation takes place by relating the components of the stress tensor, σ, to those of the rate of deformation ten- sor, D. To this end an additional function, G, called the plastic potential is introduced. The gradient of the plastic potential gives the strain direction: 𝐷𝑖𝑗 = 𝛾 𝜕𝐺 𝜕𝜎𝑖𝑗 , (1) where γ is a scalar coefficient, often called the granular fluidity. A central feature of equation (1) is that the major principal axes of the two tensors are aligned (‘coaxiality’), i.e. that stresses cause deformations preferentially in the same direction [56] [57]. In many studies the plastic potential function is chosen identical with the yield function which is called an ‘associated flow rule’. Several extension to this basic approach have been proposed, such as non-associated flow rules (e.g. [58]) or curved yield surfaces (e.g. [59]) to include effects like dilatancy and contractancy, hardening and improved energy dissipation. Furthermore, there are models to include also angular momentum (Cosserat models [60]) or models which do not involve a yield condition but assume that stress rate can be expressed as a function of the stress tensor, rate of deformation tensor and the solids fraction (“Hypoplastic models” [61] [62]). An overview of different model extensions can be found in [54]. 2 State of the art 9 In many cases of continuously flowing granular materials it has turned out that the elas- tic part of the deformation can be neglected [54] and that the material deforms at ap- proximately constant volume (“critical state approximation”) [58] [48] [63] which leads to a simplification of the constitutive relations. Plasticity models have been used to calculate the granular flow in hoppers [59], bunkers and bins [64] [58], geophysics [57] [65], blast furnaces [66] and snow avalanches [67]. The bridging of both, the frictional flow regime described by plasticity models and the kinetic-collisional flow regime described by the KTGF is still a matter of research. The most common approach is to assume that the net stress in the granular medium is simply the sum of the frictional and the kinetic stresses, each determined as it acts alone [48] [68]. This approach however is only a preliminary step towards a comprehensive theory of the intermediate flow regime. Three other approaches have been proposed [69] [70] [71] which, according to Rao and Nott [48], have ingredients useful for further devel- opments but so far are quite speculative and without a firm basis. To include heat transport in a continuum model, the corresponding balance equation requires a formulation for the effective thermal conductivity of the bulk. If the motion of the fluid phase is considered, a heat transfer coefficient between the two phases is to be incorporated also. An overview of possible heat transfer correlations between parti- cles and a fluid is given by Syamlal et al. [72]. Models for the effective thermal conductivity of the granular bulk usually assume a stagnant interstitial fluid [73] [3] [74] [75]. Still, they are also commonly used for mov- ing packed beds (e.g. [76] [77]) and also in fluidized beds [78] [79]. 2.2 Prior work on heat transfer and granular flow in moving packed beds Heat transfer between granular media and immersed surfaces has been investigated for a long time. In the 1950s Mickley and Fairbanks [80] introduced the so-called “packet- theory” to account for the heat transfer in fluidized beds. They assumed a “packet” of particles to be in contact with the heat transferring surface for a short residence time t after which it is replaced by a fresh packet from the bulk. Based on the assumption of unsteady heat conduction into a homogeneous packet, they derived an expression of the instantaneous heat transfer coefficient, αt,S0, between the packet and the surface. Ac- 10 2 State of the art cording to their formulation, αt,S0 is proportional to the inverse square root of the resi- dence time (αt,S0 ~ t -0.5). However, succeeding researchers (e.g. [81]) observed that at very short residence times the heat transfer coefficient approaches a (finite) maximum which primarily depends on the particle size. To account for this result, Baskakov [82] introduced a contact re- sistance in series with the thermal resistance of the packets introduced by Mickley and Fairbanks. The physical origin of this contact resistance initially remained unclear. Later, Schlünder [83] suggested that the thermal conductivity of the interstitial gas de- creases around the contact point where the gap width s undercuts the mean free path l of the gas molecules (Kn = l/δ > 1). Hence, the heat transfer coefficient remains finite even around the contact point where the gap width approaches zero. The concept of the series connection of the two thermal resistances (contact resistance Rc = 1/αc at the surface and penetration resistance RS0 = 1/αS0 of the packets) is illustrat- ed in Figure 4. Figure 4. Heat transfer between a surface and a contacting granular bulk material. Con- cept of a series connection of two thermal resistances, Rc and RS0. Denloye and Botteril [84] investigated the heat transfer in flowing packed beds using different materials, grain sizes and interstitial gases. They found the heat transfer coeffi- cient to increase with increasing thermal conductivity of the gas, decreasing grain size and decreasing residence time. 2 State of the art 11 Other researchers published similar works (e.g. [85]) in various experimental set-ups and generally confirmed these findings. An overview is given in the works of Obus- kovic [77] and Niegsch [86]. Of particular interest with respect to the application in solar thermal power plants are tubular heat exchanger designs. This is due to the fact that the pressure of the secondary fluid – usually water/steam – is very high in this case. The granular material may flow through the shell of the heat exchanger while the pressurized fluid flows inside the tubes. A comparison between a vertical and a horizontal tube arrangement shows that the latter has the advantage that the heat transferring surface can be distributed over the entire width of the heat exchanger – especially in a staggered tube arrangement. Hence, a greater share of the particles comes in contact with the surface. This is of major im- portance as the thermal conductivity typically is very low in moving beds. In addition, the flow is laminar and almost no mixing lateral to the main flow direction occurs [86]. The granular flow pattern in a horizontal arrangement is more complex but leaves more potential for optimization as will be shown in the following paragraphs. Several works addressed the heat transfer between horizontally arranged tubes and a moving packed bed [76] [87] [88] [89] [86] [2] [1] of which three are to be emphasized here. These three studies also investigated the granular flow field around the tubes. Niegsch [86] [2] was the first who studied the heat transfer rate between a moving bed and a bundle of horizontally arranged tubes. He investigated different staggered and aligned tube arrangements and several materials like glass beads, corundum particles, ash and quartz sand. The heat transfer rate in the bundles increased with the mass flow rate but approached an upper limit at high flow rates. Besides the overall heat transfer of the tube bank, Niegsch also measured the local heat transfer coefficient along the circumference of a single tube of the bundle. He observed small heat transfer rates at the top and below the tube and high heat transfer rates at the lateral sides. He associated the measured profile with the granular flow pattern around the tubes. To visualize the granular flow, he used layers of colored quartz sand in a transparent acrylic glass tank and observed a stagnant area on top of the tubes, high flow speeds at the side and a cavity below the tubes. Therefore, he concluded that the stagnant zone 12 2 State of the art and the cavity constitute additional thermal resistances which hamper the heat transfer in the corresponding areas. Niegsch developed an empirical model to describe the bulk solids movement and the heat transfer at the tube surface. However, his model was based on certain assumption regarding the size, shape and extent of the stagnant zone and the void zone which have not been verified. For example he assumed that the stagnant area is solely dependent on the friction parameters of the bulk and the tube surface. Takeuchi [89] performed similar experiments with a bank of one to three tube rows in a moving bed of polypropylene particles. He measured the local heat transfer coefficient at a single tube and found a profile similar to the one of Niegsch. Furthermore, he varied the relative position of the tubes and showed that a staggered arrangement yields by far higher heat transfer rates than an aligned arrangement. To measure the granular flow field he tracked a horizontal layer of tracer particles using X-ray video technique. He noticed that the stagnant zone at the tube vertex remained independent of the mass flow rate but changed with the tube arrangement. Especially at the first row of tubes and in case of large horizontal pitches he reported that the stagnant zone disappeared. In these cases significantly more heat was transfer at the top of the tube. Baumann [1] performed heat transfer measurements at elevated temperatures (> 500 °C) for the specific application of MBHE in CSP plants using quartz sand and sintered bauxite proppants. He examined two staggered tube arrangements with different hori- zontal spacing. The narrower arrangement showed distinctly higher heat transfer rates which was attributed to smaller residence times due to elevated flow speed in the con- stricted flow cross-section between the tubes. However, similar to the findings of Nieg- sch the heat transfer rate approached an upper limit with increasing mass flow rate. Baumann used particle image velocimetry (PIV) to analyze the granular flow in a trans- parent acrylic glass heat exchanger mockup. From the results he estimated the size of the stagnant zone, which primarily depended on the tube arrangement while the flow speed and the type of granular material showed little effect. Furthermore, Baumann measured the granular flow speed along the tube surface. De- pending on the tube arrangement he observed a significant acceleration of the granular flow in the lower half of the tube. He inferred that the stagnant zones on top of the tubes grew into the space between the upstream tubes and thus reduced the effective flow 2 State of the art 13 cross section. Due to the reduced cross section, he argued, the flow speed in the lower part of the tubes increased. He found the void area below the tubes to remain independent of the tube arrangement and the mass flow rate. The inclination of the flanks of the void area was of the order of the repose angle of the material. Baumann proposed a 2D CFD1-model to calculate the flow field around the tubes which also included heat transport inside the bulk. However, his model didn’t capture the for- mation of the stagnant areas at the tube vertex. Furthermore, at elevated mass flow rates (flow speeds) the model tended to considerably overestimate the heat transfer rate. Nev- ertheless, Baumann considers the continuum model a promising approach which should be followed up in future works. The three works summarized in the preceding paragraphs primarily focused on the ex- perimental investigation of moving bed heat exchangers. They revealed the overall flow pattern around the tubes and the general heat transfer characteristics. However, besides the well-known dependence of the heat transfer on the residence time, no systematic studies exist on how the thermal performance of the device can be optimized. For ex- ample, it is not known whether the size and extent of the stagnant areas, which are ham- pering the heat transfer at the top of the tubes, can be influenced by certain measures. Especially, no comprehensive model of the granular flow around the tubes exists. The model of Niegsch relies on simplifying assumptions regarding the flow pattern which is also the case for the model of Lee et al. [90]. In the model of Baumann less ‘a priori’ information about the flow pattern is included but the stagnant area is not captured and the results regarding the heat transfer rate obviously are erroneous, especially at high mass flow rates. So, obviously there is still a need for further work regarding the model- ling and the understanding of moving bed heat exchangers. 2.3 Scope of investigation As shown in the previous section, a reliable and comprehensive tool for the design and optimization of moving bed heat exchangers with staggered, horizontal tubes is still missing. This work contributes to filling this gap by the accurate modelling of the hy- 1 CFD = Computational Fluid Dynamics 14 2 State of the art dro- and thermodynamic processes in the device. To this end, a two-step proceeding is chosen: Firstly, a comprehensive continuum model is developed which captures both the granu- lar flow field around the tubes as well as the heat transport inside the granular material. In particular, it captures the influence of granular-specific flow phenomena like stagnant and void zones. As it is expected that effects occurring on particle scale will also be relevant for the heat transfer process, the continuum model is to be supplemented by a discrete particle model providing insight into these phenomena. The results of the dis- crete particle model are used to modify the continuum model and to increase its preci- sion. The simulation results of the different models are validated by means of appropriate experiments. This includes the measurement of the flow pattern of the moving bed around the tubes. In addition, the heat transfer at a single heat exchanger tube is to be assessed in detail. Secondly, a simplified model of the heat transfer between the tubes of an MBHE and the moving bulk is to be developed. This model is based on the findings of both the dis- crete particle and the continuum model and provides a computationally efficient way to assess different heat exchanger designs. Furthermore, such a simplified model can be used in system simulations of industrial or power plant processes using a moving bed heat exchanger. The validated models are used to assess the influence of different operating, design and material parameters on the granular flow pattern with respect to the thermal perfor- mance of the heat exchanger. The work is structured as follows: In chapter 3 the experimental infrastructure and measurement setup used for the valida- tion of the numerical models is introduced. In chapter 4 the continuum model and the discrete particle model are introduced as well as the considered boundary conditions and calibration procedures of the two models. In chapter 5 the validation of the two models is addressed. In a first step the flow pattern is validated which leads to a modification of the continuum model to account for particle-scale effects. The modified continuum model is subsequently validated based on the measurement of the heat transfer at a sin- gle tube. 2 State of the art 15 Chapter 6 focuses on a simplified model for the heat transfer between a single tube and the bulk. To this end, the dependence of the flow pattern on several influencing parame- ters is assessed using the discrete particle model. The results of these parameter studies are subsequently used to develop the simplified model. The results of the work are summarized in chapter 7. 16 3 Experimental setup 3 Experimental setup The moving bed heat exchanger test infrastructure used for this study is shown in Figure 5 (left) and has been described in detail in [91]. It basically consists of a storage con- tainer for the granular material, a particle heat exchanger and a conveyor device. From the storage tank the material flows, driven by gravity, through the test section where arbitrary heat exchanger geometries can be inserted. The chain conveyor returns the material back from the heat exchanger outlet into the storage container. Figure 5. Left: Moving bed heat exchanger test rig. Right: Heat exchanger mockup. The rectangularly framed area marks a representative section of the tube arrangement, which is investigated numerically in later chapters. The heat exchanger test rig was filled with quartz sand of an average grain size of dp = 0.6 mm and dp = 1.2 mm, respectively, depending on the experiment to be conduct- ed. The average flow speed in the free cross-section of the heat exchanger was evaluat- 3 Experimental setup 17 ed from the mass flow rate of the tube chain conveyor and the bulk density (ρb ≈ 1600 kg/m³). The mass flow rate depended linearly on the rotational speed of the chain conveyor and was calibrated beforehand. The flow speed was varied in the range of uref = 1.5…6.0 mm/s. The measurement uncertainty of the reference velocity is ap- proximately 5 % (see appendix A.3). To allow a visual inspection of the bulk flow, the heat exchanger was replaced by an acrylic glass mockup which is shown in Figure 5 (right). When operated, the bulk mate- rial entered at the top and flowed around the tubes. Two types of experiments were conducted: On the one hand the granular flow field in- side the heat exchanger mock-up was examined using PIV analysis; on the other hand, the local heat transfer rate from a single tube to the bulk was examined. 3.1 Measurement of flow field Regarding the experimental examination of granular flows, different measurement tech- niques exist [17]. Though restricted to the inspection of the visible surface of the bulk flow, particle image velocimetry (PIV) and particle tracking (PT) are commonly used. They are easy to handle and are applicable to relatively large objects and measurement setups. In this work, PIV is used to examine the granular flow inside the mockup. PIV is based on the analysis of successive images of the flow. The images are subdivided by a grid of stationary control volumes. By detecting recurring patterns and determining their dis- placement in a pair of images the local velocity of the flow can be calculated. Unlike for fluids, in granular flows no extra tracer particles are needed as the patterns are formed by the grains themselves. A more detailed description of PIV in the context of granular flows can be found elsewhere [1] [92]. 18 3 Experimental setup Figure 6. Schematic of the experimental setup for PIV analysis. Figure 6 shows a schematic of the experimental setup. Pictures were taken at a rate of 25 Hz through the transparent acrylic glass front of the mockup. To reduce temporal fluctuations of the local flow speed the results were averaged over a sequence of 40 pairs of images. The grid discretization was determined based on the expectant flow speed. The dis- placement between two images should be about a quarter of the size of the control vol- umes [93]. As the flow is moving into vertical direction the flow speed in y-direction is generally higher than in x-direction. At a frequency of 25 Hz and an expectant flow speed of u ≈ 10 mm/s between the tubes, this leads to a vertical discretization of Δy = 1.6 mm. As the flow speed in x-direction is much smaller, a finer discretization of Δx = 0.4 mm is used. This means that a control volume contains only two or three parti- cles which is far below the recommended amount of 10…25 particles [93] for conven- tional PIV measurements. However, in the case at hand the images are rich in contrast and the PIV-algorithm doesn’t necessarily track the particles but rather an optical struc- ture in the image which might also be of sub-particle size. An exemplary vector plot obtained from the PIV analysis is shown in Figure 7. 3 Experimental setup 19 Figure 7. Exemplary vector plot obtained from PIV analysis. In Figure 8 the velocity profile along the horizontal plane y = 0 (see Figure 7) is depict- ed. Four peaks of the flow speed are found in Figure 8, which belong to the four gaps between the tubes. Due to continuity the mean flow speed, �̅�, between the tubes (over the entire cross section of the mockup) must be �̅�(𝑦 = 0) 𝑢ref = 𝑠H 𝑠H − 𝐷T = 3.7, (2) with sH being the horizontal spacing of the tubes and DT being the outer tube diameter. As visible from Figure 8, the mean flow speed between the tubes in the measurement is u/uref ≈ 1.7, which is attributed to the fact that the flow is slowed down due to surface friction of the acrylic glass front wall. Thus, for comparison of the measurement and simulation results (see section 5.1) the surface friction of the front wall must be taken into consideration in the numerical model. The flow inside the mockup cannot be meas- ured directly. The flow speed between the tubes deviates by up to 21 % from the average value of u/uref ≈ 1.7. Possible origins for the uneven velocity distribution are disturbing influ- ences of the front wall, the preceding row of tubes, the change of the free cross-section between two subsequent tube rows due to the guiding plates at the lateral sides of the 20 3 Experimental setup mockup (see Figure 5 (right)) and the tube chain conveyor (for further discussion of these influences see Appendix A.3). Figure 8. Horizontal velocity profile (velocity magnitude) at y = 0 (see Figure 7). Inaccuracies of the PIV measurement results are mainly caused by two sources of er- rors: Firstly, by the calibration procedure to find the conversion factor from image to physical space. Here, a conversion factor of 16.6 px/mm was determined based on a known length of 76 mm (1262 px). Assuming that the measurement accuracy is 10 px, the resulting error is approx. 0.8 %. Secondly, the PIV evaluation algorithm causes an inherent evaluation error of about 0.2 px [93]. Averaging over a series of 39 pairs of images reduces the inaccuracy to 0.03 px. At a frequency of 25 Hz this corresponds to an uncertainty of the calculated flow speed of ± 0.05 mm/s. This uncertainty is especially relevant in areas where the flow speed is low such as in the stagnant area. Based on the velocity in the free cross section of the mockup (uref = 4.0 mm/s) during the experiment, the resulting uncertainty is 1.25 %. These two sources of measurement errors thus lead to an overall accuracy of approx. 2 % of the PIV measurement results. 3 Experimental setup 21 3.2 Measurement of local heat transfer coefficient The local heat transfer coefficient between a single tube and the granular flow was eval- uated, following the approach of Takeuchi et al. [89]. It is based on measuring the local surface temperature of a heated test tube under operation. The local heat transfer coefficient is calculated from the heat transfer rate and the local temperature difference, as shown below. The test tube, made of polyethylene (PE), substituted one of the stainless steel tubes at a central position in the heat exchanger mockup. The PE-tube was covered by a stainless steel foil which was electrically contacted at both sides of the tube. By applying an elec- tric current, the foil was heated up due to its ohmic resistance. The surface temperature TW of the foil was measured by three thermocouples. The thermocouples were located under the foil in three notches carved into the PE tube with the notches being evenly distributed around the tube circumference (120°-displacement between the thermocou- ples). A drawing of the test tube is displayed in Figure 9. The heat generation rate q̇ is calculated according to equation (3), assuming a homoge- neous current density. Further below it is shown that this assumption is justified. �̇� = 𝑃el 𝐴 = 𝐼2𝑅 𝐴 (3) In equation (3), A is the area of the foil, I is the current and R is the ohmic resistance of the foil. Assuming adiabatic conditions at the inner surface of the tube (stagnant air), the total of the thermal energy of the foil is transferred to the surrounding granular medium after reaching a stationary state. With respect to the reference temperature Tref the local heat transfer coefficient αloc(ω) – at the position angle ω of the thermocouples – was calculated: 𝛼loc(𝜔) = �̇� 𝑇W(𝜔) − 𝑇ref (4) The reference temperature Tref was defined as the inflow temperature of the heat exchanger mock up which was measured in the center of the heat exchanger inlet. By rotating the test tube around its axis, a 360°-profile of the local heat transfer coefficient around the tube was derived. 22 3 Experimental setup Figure 9. Sketch (left) and photo (right) of the test tube for the measurement of the local heat transfer between tube surface and granular flow. To evaluate the local heat transfer rate q̇ according to equation (3) it was vital that the current density was distributed evenly around the circumference of the tube at the loca- tion where the surface temperature of the foil was measured. This prerequisite was test- ed beforehand: The tube was heated in the environment of stagnating air and the distri- bution of the surface temperature was measured using an infrared camera. The tempera- ture varied by ΔT ≈ 0.5 K around the circumference of the tube which corresponds to an uncertainty of the current density of less than 2.5 %. Together with the measurement uncertainty of the quantities R, I and A in equation (3) the total uncertainty of the heat flux density q̇ is approximately 4.5 % (see appendix A.2). The test tube was inserted into the heat exchanger mockup at a central position. At the beginning of every measurement series the flow rate was adjusted to the desired value and the foil was heated for 30 minutes to reach a stationary state. During the measure- ment the tube was rotated clockwise in steps of Δω = 10° beginning at the top of the tube (ω = 0°) and waiting about five minutes at every step before taking the local sur- face temperature. The local surface temperature appeared relatively constant during the measurement (± 0.1 K) except at ω ≈ 150° and ω ≈ 210° where fluctuations are ob- served of about ± 1 K were observed. These fluctuations are attributed to the formation of a void area below the tube which will be discussed in more detail in chapter 5. Figure 10 (left) shows exemplary profiles of the reference temperature Tref at the inlet of the mockup and of the surface temperature TW at the tube wall (of the stainless steel foil). The reference temperature didn´t stay at a constant level but slightly increased during the experiment. The angle-ranges “TC1” to “TC3” denote the ranges recorded by the three different thermocouples. 3 Experimental setup 23 Figure 10 (right) shows the corresponding profile of the local heat transfer coefficient. The profile is approximately symmetric with maxima at the lateral sides of the tube (ω ≈ 90°, ω ≈ 270°), a minimum below the tube (ω = 180°) and intermediate values at the top of the tube (ω ≈ 0°, ω ≈ 360°). The results are discussed in detail in section 5.2.3 where they are compared them to corresponding simulation results. Figure 10. Exemplary measured temperature profiles (left). Corresponding profile of local heat transfer coefficient (right). The accuracy of the measurement results according to equation (4) depends on the one hand on the accuracy of the heat flux q̇ discussed above, on the other hand on the accu- racy of the temperature difference (ΔT = TW(ω) - Tref) which is approximately 0.2 K. At the lateral sides of the tube with (ΔT)min ≈ 19 K this corresponds to a relative uncertain- ty of approximately 2 % (see also appendix A.2). Hence, the total uncertainty of the measured local heat transfer coefficient αloc is esti- mated to approximately 7 %. Three different tube arrangements and two different particle diameters were investigat- ed. The examined configurations are given in Table 1. For each configuration four dif- ferent mass flow rates/reference flow speeds were tested. The measurement results of the four different configurations are presented in section 5.2.3 together with the corre- sponding simulation results of the continuum model. 24 3 Experimental setup Table 1. Investigated configurations. Tube arrangement Horizontal tube spacing sH Vertical tube spacing sV Mean particle diameter dp Reference flow speed uref mm mm mm mm/s TA1 57 40 0.6 1.5 / 3.0 / 4.5 / 6.0 TA2 57 60 0.6 1.5 / 3.0 / 4.5 / 6.0 TA3 37 60 0.6 1.5 / 3.0 / 4.5 / 6.0 TA3 37 60 1.2 1.5 / 3.0 / 4.5 / 6.0 4 Modelling 25 4 Modelling In this work two different modelling approaches are used – a continuum model (CFD) and a discrete particle model (DEM) – to calculate the granular flow around the horizontally arranged tubes in a moving bed heat exchanger. The investigated geometry is depicted in Figure 11. It consists of a representative section of three rows of tubes arranged in a staggered manner (see framed area in Figure 5 (right), p.16). The granular material enters the geometry at the top, flows around the tubes and leaves at the bottom. The continuum model simulations are in 2D, transient. The motion of the interstitial fluid (air) and the heat transport inside the bulk are included. The model equations are given in the following section. The geometry is split along the vertical center line (see Figure 11, left), and only one half is simulated using symmetry boundary conditions in x-direction. The boundary conditions for the continuum model are as follows: At the inlet the pres- sure of the fluid phase is set, the volume fraction for both phases and the inlet tempera- ture (same for both phases). The boundary conditions at the tube walls deserve further explanation and are discussed in detail in section 4.1.2. At the outlet the gradient of the volume fraction normal to the boundary and the temperature for both phases is set. Fur- thermore, a fixed outlet velocity is defined which is the same for both phases. The mesh of the continuum model was refined until the simulation results converged. At the tube surface a local refinement was introduced to resolve the temperature gradient at the wall. This lead to a mesh size of about one particle diameter (0.6 mm) at the tube surface while in the rest of the geometry the mesh size was about three particle diame- ters in each spatial direction. The model equations of the continuum model follow in section 4.1. In the discrete particle simulation are in 3D as quantities like the porosity and the coor- dination number2 of the moving bed are three-dimensional phenomena. In z-direction the geometry is 13 to 20 particle diameters thick (depending on the studied case) and is confined by walls. The rear wall is friction-less in all cases. The friction coefficient at the front wall is non-zero for the comparison with the validation experiment in sec- 2 The coordination number is the number of contacts of a particle with its neighbors. It is of major im- portance for the heat transport in a packed bed [44] [74]. 26 4 Modelling tion 5.1. For all remaining simulations it is set to zero. Further explanation is given in section 4.3.2. In x-direction periodic boundaries are used in the discrete particle model. At the tube walls a wall friction coefficient is assigned which will be discussed in section 4.3. Heat transport as well as the interstitial fluid is not considered in the discrete particle model. After the simulation geometry was filled randomly with particles a fixed outlet velocity was specified by means of a moving horizontal plate. The model equations for the dis- crete particle model follow in section 4.2. Figure 11. Simulation geometry with boundary conditions. Side view only for DEM model (3D). Subscripts: “s” = solid phase, “f” = fluid phase, “i” is substitute for both phases in the continuum model. Remaining subscripts and symbols are explained in the following sections. 4 Modelling 27 The following Figure 12 shows exemplary simulation results (velocity magnitude of the particles/granular phase) of both models. Figure 12. Exemplary simulation results of both models. Left: Velocity magnitude of granular phase (Continuum model). Right: Magnitude of translational velocity of parti- cles (Discrete particle model). In chapter 5 the two models are compared to each other and to experimental flow meas- urement results. The comparison is performed based on the unit cell in the center of the geometry (see Figure 12, right). The areas above and below the unit cell, which are dominated by boundary effects, are not considered in the comparison. At this point, some differences regarding the boundary conditions between the two models and the experiment have to be mentioned. In the DEM model the particles are generated above the geometry (with a low packing fraction) and fall into the geometry. In contrast to that, in the continuum model the inlet packing fraction is high while the flow speed is low. This leads to deviations regarding the simulation results outside the unit cell. However, simulations with more tube rows arranged in vertical direction show a recur- ring, constant flow pattern for multiple unit cells. Therefore, in this work the geometry in Figure 12 with only one representative unit cell is used to analyze the flow pattern. 28 4 Modelling Furthermore, the fact that the flow pattern recurs constantly after the first row of tubes justifies the comparison of the simulation results to experimental data from a unit cell in the center of the heat exchanger geometry (see Figure 5). 4.1 Continuum model equations The continuum model mainly follows the works of Srivastava et al. [63] and Schnei- derbauer et al. [94]. The moving bed is modelled as two interpenetrating continua (Eu- ler-Euler approach), one representing the granular phase and the other one representing the gas phase. Following the approach of Ishii [95], for both phases averaged balance equations are solved as shown in equation (5)-(7). 𝜕 𝜕𝑡 𝜀i𝜌i + 𝛁 ∙ (𝜀i𝜌i𝒖i) = 0 (5) 𝜕 𝜕𝑡 (𝜀f𝜌f𝒖f) + 𝛁 ∙ (𝜀f𝜌f𝒖f𝒖f) = −𝜀f𝛁𝑝 + 𝛁 ∙ 𝜀f𝝉f − 𝛽(𝒖f − 𝒖s) + 𝜀f𝜌f𝒈 (6) 𝜕 𝜕𝑡 (𝜀s𝜌s𝒖s) + 𝛁 ∙ (𝜀s𝜌s𝒖s𝒖s) = −𝜀s𝛁𝑝 − 𝛁 ∙ 𝝈s + 𝛽(𝒖f − 𝒖s) + 𝜀s𝜌s𝒈 (7) 𝜕 𝜕𝑡 (𝜀i𝜌iℎi) + 𝛁 ∙ (𝜀i𝜌i𝒖iℎi) = −𝜀i 𝜕𝑝 𝜕𝑡 − 𝝉i: 𝛻𝒖i + 𝛻 ∙ 𝜀i𝜆i∇𝑇i ± 𝑄ij (8) Here, the indices “s” and “f” denote the solid phase and the fluid phase, respectively, and “i” and “j” are substitutes for both phases. ε is the volume fraction, u the velocity, ρ the density and h the enthalpy of the corresponding phase. For the interphase drag coef- ficient β a correlation of Gidaspow et al. [50] is used which applies the Ergun equation for high volume fractions: 𝛽 = { 150 𝜀s 2𝜇f (1 − 𝜀s)𝑑p2 + 1,75 𝜀s𝜌f|𝒖f − 𝒖s| 𝑑p if 𝜀s > 0.2 0,75𝐶D 𝜀s(1 − 𝜀s)𝜌f|𝒖f − 𝒖s| 𝑑p (1 − 𝜀s) −2.65 if 𝜀s ≤ 0.2 (9) 4 Modelling 29 𝐶D = { 24 𝑅𝑒f (1 + 0.15𝑅𝑒f 0.687) if (1 − 𝜀s)𝑅𝑒f < 1000 0.44 if (1 − 𝜀s)𝑅𝑒f ≥ 1000 (10) For the gas phase shear stress τg a Newtonian closure is used: 𝝉f = 2𝜇f𝑫f , (11) where 𝑫 = 1 2⁄ (𝛁𝒖 + (𝛁𝒖)T) is the rate-of-deformation tensor. In the energy conservation equation equation (8) two quantities remain to be defined: The thermal conductivity λi of each phase and the interphase heat-exchange term Qij. The heat exchange term depends on the temperature difference between the two phases and the interfacial area Afs: 𝑄ij = 𝛼fs𝐴fs(𝑇i − 𝑇j) (12) The interphase heat transfer coefficient αij is calculated using a correlation proposed by Gunn [96] which is valid for a wide range of packing fractions and Reynolds numbers (see appendix B.4). The thermal conductivities of the two phases, λs and λf, are effective transport properties and are not to be confused with the microscopic properties of the pure substances. They are defined based on the effective thermal conductivity λS0 of the bed which is calculat- ed according to a correlation of Zehner and Schlünder [3]. Subsequently, λS0 is split up into the two phase conductivities following the approach of Kuipers et al. [78] (see ap- pendix B.3). It is important to note that this approach assumes an isotropic thermal conductivity. In the direct vicinity of a wall, however, the packing structure becomes anisotropic, and due to the increased voidage in the near wall region the thermal conductivity decreases. In particular, the gas gap between the wall and the first layer of particles constitutes a significant thermal resistance. This thermal contact resistance is incorporated into the thermal boundary condition at the wall as will be described in section 4.1.2.1. 30 4 Modelling 4.1.1 Frictional stress model What remains to be defined is the stress tensor of the solid phase σs in equation (7). As shown in section 2.1.2, the constitutive relations for granular flow models depend on the investigated flow regime. As the packing fraction in a moving bed is close to the maxi- mum packing fraction a plasticity model is applied to account for the frictional interac- tion between particles. It is assumed that the granular flow is completely governed by friction so that kinetic and collisional influences can be neglected. The applied frictional closure follows the works of Jackson [97] and Tardos et al. [56] and has been used to model the discharging process of granular material from bins [63] [94] and hoppers [56]. The frictional stress is written in a compressible Newtonian form: 𝝈s = 𝝈fr = 𝑝fr𝑰 + 𝝉fr = 𝑝fr𝑰 + 2𝜇fr𝑺s (13) Here, 𝑺s = 𝑫s − 𝑡𝑟(𝑫𝐬) is the deviator of the strain rate tensor. The two parameters pfr and µfr are called frictional pressure and frictional viscosity, respectively, and are ex- plained in more detail in the following paragraphs. The frictional pressure accounts for the repulsive forces between grains and prevents the particle assembly from being compressed beyond a maximum packing fraction εmax. In soil mechanics pfr is often combined with a yield function to describe the effect that a particle assembly may dilate or compact under shear movement. The theoretical frame- work in this context is called the critical state theory. A simplifying assumption of this theory is that the dilation and compaction effects are small and that the state of the ma- terial is close to the so-called critical state where it deforms without any volume change (“Critical state approximation” [48]). This assumption has been proven to be justified for the flow in discharging bins [63] and hoppers [56]. Compared to this, the flow around a horizontal tube as it is being investigated in this work is more complex. Therefore, the critical-state-assumption has been tested before- hand and it was found that the influence on the simulation results is small. For example, the heat transfer rate changes by less than two percent. Hence, the cirtical-state- assumption is considered to be justified in the case at hand. The extended model (with- out critical-state approximation) is given in the appendix B.5. As a consequence pfr is a function only of the packing fraction εs. Several formulations for the critical state pres- sure have been proposed [98]. In this work the form of Johnson and Jackson [99] is used: 4 Modelling 31 𝑝𝑓𝑟 = 0.1𝜀𝑠 (𝜀𝑠 − 𝜀𝑚𝑖𝑛) 𝑟 (𝜀𝑚𝑎𝑥 − 𝜀𝑠)𝑠 (14) The parameter εmin is the minimum packing fraction above which frictional interaction between particles occurs (“loosest random packing”). If the packing fraction εs under- cuts the minimum packing fraction εmin, the frictional pressure pfr is set to zero. The quantities in equation (14) are set the same as in the work of Srivastava et al. [63] and are listed in Table 2. The second parameter in equation (13), the frictional viscosity μfr, accounts for the shear stresses inside the granular material. One of the major challenges regarding the viscosity is that it should model the transition of a granular material from a “solid” (static) to a “flowing” state. This behavior manifests in granular flows by the formation of static areas where the particle aggregation behaves like one rigid body. The transition from a static to a flowing state takes place when the shear stress reaches a certain threshold. The simplest assumption for the threshold at which yielding occurs is a coulomb friction correlation (τfr = μi ∙ pfr), where μi is an inner friction coefficient of the material. Based on this assumption Schaeffer [100] proposed a formulation for the granular viscosity: 𝜇fr = 𝑝fr𝜇i 2|𝑺s| (15) The denominator of equation (15) contains the Euclidian norm of the deviatoric part of the strain-rate tensor |𝑺𝐬| = √𝑆ij𝑆ij. If the strain-rate (shear) approaches zero the viscos- ity diverges, which ensures the existence of a yield stress: the stress does not vanish when the flow stops [17]. The internal friction coefficient µi of the granular material is related to the angle of in- ternal friction ϕi: 𝜇i = √2sin (𝜙i) (16) Jop et al. [70] suggested an empirical correlation for a dynamic friction coefficient μi(Is) depending on the dimensionless number Is. Is is called the inertial number and is a measure of whether the material exhibits a more solid-like or fluidized behavior. One could say, it is an indicator of the flow regime where the process takes place. The model of Jop et al. was primarily developed to describe granular flows with a free surface 32 4 Modelling down an inclined plane where a wide range of inertial numbers occur. In contrast to that, the present study investigates a confined flow. Furthermore, the inertial number of the investigated flow is very low (Is ≲ 0.01), which indicates that the process takes place in the quasi static flow regime where μi can be regarded as constant [17]. There- fore, the formulation of Schaeffer et al. using a constant value of μi according to equa- tion (16) is applied. 4.1.2 Boundary conditions at the tube walls At the tube surfaces the boundary conditions for the energy equation (8) and the mo- mentum equations (6) and (7) have to be defined. 4.1.2.1 Temperature boundary condition As mentioned in section 2.2, heat transfer between a surface and a contacting bulk in- volves a thermal contact resistance due to the gas gap between the surface and the first layer of particles. This contact resistance leads to a temperature drop ΔT = TW - T0, di- rectly at the surface as illustrated in Figure 13. Figure 13. Schematic of the heat transfer between a bulk material and an immersed sur- face. To avoid resolving this temperature drop, which would require a spatial resolution far below the particle size, it is directly incorporated into the boundary condition. In other 4 Modelling 33 words, the temperature at the boundary is not the temperature TW of the wall itself but a temperature T0 which is reduced in accordance with the contact resistance (αc) -1 and the heat transfer rate q̇: 𝑇0 = 𝑇𝑊 − �̇� 𝛼𝑐 (17) The contact heat transfer coefficient, αc, is calculated according to a formulation of Schlünder [3] which has been developed for resting bulks but has also been applied to moving beds [86]: 𝛼c = 𝜑𝛼WP − (1 − 𝜑)𝛼con + 𝛼rad (18) The contact resistance consists of three parts: • Heat transfer to contacting particles (φ∙αWP). • Heat transfer by conduction through the gas phase to the second layer of parti- cles (1-φ)∙αcon. • Heat transfer by radiation αrad. αWP denotes the heat transfer coefficient from the wall to a contacting particle which is being weighted by the factor φ denoting the fraction of the wall covered by contacting particles. 𝜑 = 𝑁p 𝜋𝑑p 2 4 𝐴W = 𝑁𝑝𝐴p 𝐴W (19) Here, Np is the number of contacting particles, Ap is the projected area of a particle and AW is the total surface area of the wall. Schlünder states an empirical value of φ ≈ 0.8 for random pebble beds. The second addend in equation (18), αcon, denoting the heat transfer coefficient to the second layer of particles, is accordingly weighted by the complement of the surface fraction (1 - φ). The magnitude of the contact heat transfer coefficient αc is usually dominated by the heat transfer coefficient to a contacting particle αWP. Schlünder evaluates αWP based on the thermal conduction through the gas gap between the wall and the contacting particle and gives the following equation: 34 4 Modelling 𝛼WP = 4𝜆f 𝑑p {[1 + 2(𝑙 + ℎR) 𝑑p ] ln [1 + 𝑑p 2(𝑙 + ℎR) ] − 1} (20) Equation (20) on the one hand includes two properties of the fluid, namely the mean free path way l and the thermal conductivity λf, on the other hand, two properties of the particle and the wall, namely the particle diameter dp and the total surface roughness of particles and wall hR. Maximum values of αWP are achieved at high conductivities λf, small particle diameters dp, small surface roughness hR, and small mean free path ways l. Of the last two quanti- ties the larger one dominates the process. So, at ambient pressure (l < 100 nm) the sur- face roughness (hR ≈ 1…10 µm) is usually the dominating quantity. The equations to calculate the remaining two heat transfer coefficients in equation (18), αcon and αrad as well as the mean free path of the gas can be found in the appendix B.2. The thermal boundary condition according to equation (17) defines the temperature at the boundary while the heat flux adjusts accordingly. In contrast to this, in the experi- mental setup described in section 3.2, the heat flux at the tube surface is defined and the corresponding temperature profile is measured. To compare the experiment with the simulation results (section 5.2.3), the wall temperature is set according to the measured temperature profile (see Figure 10, left). 4.1.2.2 Velocity boundary condition The velocity boundary condition for the granular phase must take into account that the granular material may either slip along the wall or stick to the wall. This is of particular interest with regard to the transition from a static to a flowing regime in the upper part of the tube. At the top of the tube the particles are at rest but start to slide along the cir- cumference when a certain angle of inclination is reached. In order to model this process the Coulomb friction law is used to define the granular velocity at a solid wall. The wall shear stress in slip case is 𝜏sl = 𝜇W𝑝fr , (21) where μw is a wall friction coefficient, τsl is the (scalar) shear stress at the wall and pfr is the normal stress. The wall shear stress is then compared to the viscose stress in tangen- tial direction τfr,t inside the particle assembly close to the wall [94] [99]: 4 Modelling 35 𝜏fr,t = 𝜇fr [ 𝑑𝒖t 𝑑𝒏 ] 𝑊 (22) The term [ 𝑑𝒖𝑡 𝑑𝒏 ] 𝑊 is the gradient of the tangential part of the velocity ut in the direction of the unit normal vector n of the wall. By comparing equation (21) and equation (22) the (vectorial) shear stress at the wall is determined: 𝝉W = − 𝒖s |𝒖s| { 𝜏fr,t if 𝜏sl > 𝜏fr,t 𝜏sl if 𝜏sl < 𝜏fr,t (23) If the viscose stress τfr,t falls below the shear stress in slip case τsl (τsl > τfr,t) no sliding occurs at the wall. In case of τsl < τfr,tan sliding occurs and the velocity at the boundary is calculated from τW. For the gas phase a no-slip boundary condition is used at the tube wall. 4.2 Discrete particle model The modelling approach used in this work was first introduced by Cundall and Strack [20] and is known in literature as Discrete Element Method (DEM). In contrast to con- tinuous model approaches, DEM tracks the motion of every single particle of the bulk. The method has become very popular in the last decades due to the rapid increase of available computational power and has been applied in various fields. The method is based on solving the Newtonian equations of motion for each particle for translational and rotational motion. ∑𝑭i,k k = 𝑚i�̈�i (24) ∑𝑴i,k k = 𝐽i�̇�i . (25) The acting forces and moments include body forces such as the gravitational force as well as external forces and moments which, for example, may originate from contacts with other particles and walls. When solving the equations of motion, two grains may turn out to overlap at the end of the time step. This overlap is interpreted as the elastic deformation which occurs for particles under stress [23]. The overlap is subsequently translated into tangential and normal interaction forces using certain contact models. 36 4 Modelling 4.2.1 Contact models Several contact models have been developed for the normal interaction force between particles [101]. Here, a linear hysteresis model is used which first has been introduced by Walton and Braun [25]. The instantaneous normal interaction force is calculated as follows: 𝐹n t = { min(𝐹n t−Δt + 𝑘ulΔ𝑠n, 𝑘l𝑠n t) , Δ𝑠n ≥ 0 max(𝐹n t−Δt + 𝑘ulΔ𝑠n, 0.001𝑘l𝑠n t) , Δ𝑠n < 0 (26) The model distinguishes between loading case (Δsn ≥ 0) and unloading case (Δsn < 0) where Δ𝑠n = (𝑠n t − 𝑠n t−Δt) is the change in normal overlap between two time steps. In the unloading case the normal force is limited to the value of 0.001𝑘l𝑠n t to ensure that no negative (attractive) normal force occurs. The inelastic nature of particle contacts is modelled by two different contact stiffnesses kl and kul for the loading and the unloading case, respectively. kl and kul are related to each other by the coefficient of restitution e: 𝑒 = −√𝑘l 𝑘ul⁄ (27) The contact stiffness for the loading case kl is calculated by material parameters: 𝑘l = 𝐸1𝑑p,1𝐸2𝑑p,2 𝐸1𝑑p,1 + 𝐸2𝑑p,2 (28) Here, indices “1,2” denote the two elements forming the considered contact, which might either be a particle-particle contact or a particle-wall contact. dp is the particle diameter, and E is the Young´s Modulus of the contacting elements [102]. The tangential interaction forces are calculated according to the following elastic- frictional force model: 𝐹t t = { min(𝐹t t−Δt + 𝑘lΔ𝑠t, 𝜇s𝐹n t) , if no sliding occurs min(𝐹t t−Δt + 𝑘lΔ𝑠t, 𝜇d𝐹n t) , if sliding occurs (29) As one can see from equation (29) the tangential contact force 𝐹t t between two contact- ing elements evolves with the relative tangential displacement at the current time Δst. Sliding occurs, if 𝐹t t exceeds the limit of (𝜇s ∙ 𝐹n t), with μs being the static friction coef- 4 Modelling 37 ficient. Once 𝐹t t falls below the value of (𝜇d ∙ 𝐹n t), with μd being the dynamic friction coefficient, the contact is considered non-sliding again. For a reduced computational expense spherical particles are used in the simulations. Unfortunately, with this simplification the model loses its ability to predict the high re- sistance of non-spherical particles against a rolling motion. A common way to remedy this deficit is to introduce a resistive torque applied to contacting particles. The concept is known as “rolling friction” [103], and different classes of rolling friction models have been introduced. In this work an elastic-plastic spring dashpot model is used, commonly referred to as “Model C” [103], using a “rolling friction coefficient” μroll as an input parameter. The resistive torque 𝑴roll t (at time t) is defined as follows [102]: 𝑴roll t = min(|𝑴roll,e t |,𝑀roll,lim) 𝑴roll,e t |𝑴roll,e t | (30) Mroll,lim is a limiting value of the torque, depending on the normal force Fn, the rolling resistance (rolling friction) coefficient µroll, and the rolling radius Rroll: 𝑀roll,lim = 𝜇roll𝑅roll𝐹n (31) In case of mono-sized, spherical particles the rolling radius is Rroll = dp/2. The rolling friction parameter μroll can be interpreted as the tangent of the maximum angle of a slope on which the rolling resistance moment counterbalances the moment produced by gravity in the particle. µroll is usually calibrated from experimental data as shown in the following section. Below the limiting value Mroll,lim, the resistive torque is allowed to vary continuously according to a linear elastic model: 𝑴roll,e t = 𝑴roll t−Δt − 𝑘roll𝝎relΔ𝑡 (32) 𝑴roll t−Δt is the resistive torque at the previous time step, ωrel is the relative angular veloci- ty between the two contacting particles (wall and particle) and kroll is the ‘rolling stiff- ness’: 𝑘roll = 𝑅roll 2 𝑘l (33) 38 4 Modelling Due to the elastic part 𝑴roll,e t of the resistive torque, discontinuities as they occur in other rolling friction models are avoided [103]. 4.2.2 Local averaging The DEM basically yields location and velocity of every single particle at a given time. To compare the DEM-simulation results to those of the continuum-model and of the PIV measurements, they have to be averaged in space and time. For this purpose locally fixed control volumes (CVs) have to be defined wherein the averaging is carried out. Lätzel [23] investigated two different averaging methods to compare his DEM simula- tions to experimental data of an annular shear cell. His basic averaging formalism for obtaining an averaged quantity Q inside a control volume V is 𝑄 = 1 𝑉 ∑𝑤i V𝑉i𝑄i 𝑖 𝜖 𝑉 , (34) with Vi being the particle volume and Qi being the considered quantity attributed to par- ticle i. In case of the averaged velocity (Q = ux, Q = uy) the quantity is averaged based on the number of particles Ni inside the control volume: 𝑄 = 1 𝑁i ∑𝑤i V𝑄i 𝑖 𝜖 𝑉 , (35) The parameter 𝑤i V is the weight of the particles contribution to the average. In this work the simplest choice of 𝑤i V is used which is 𝑤i V = { 1, if the center of the particle lies inside the CV 0, otherwise . (36) Lätzel showed that the method is sufficiently precise as long as the diameter of the CV is greater than the particle diameter. Consequently, the size of the CVs is set in such a way that this requirement is met. After averaging in space, time averaging is straight- forward by taking the mean over multiple time steps. Depending on the quantities to be analyzed, the CVs may either by of rectangular shape, similar to the 2D-CVs in the PIV analysis (see section 3.1), or they may be arranged in circular layers, especially when analyzing quantities along a tube surface. 4 Modelling 39 4.3 Model parameters The reference parameters for both models are listed in Table 2. The material properties of the particles (the solid phase) are chosen to match those of quartz sand3. The surface friction coefficient, μW, of the tubes takes different values depending on the experimental set-up with which the simulation results are to be compared. In the exper- imental setup for measuring the heat transfer coefficient between tube and bulk, the tube surface is covered by a stainless steel foil (see section 3.2). For the surface friction coef- ficient of the foil a value of μW,1 = 0.25 is used. This value is based on the data of Bau- mann et al. [1] who measured the surface friction for various granular materials and obtained values of μW = 0.2 – 0.5. For quartz sand on polished stainless steel he gives a value of μW = 0.25. The surface friction coefficient of the bare steel tubes (without the foil) used in the PIV measurement could not be measured directly. It is expected to be higher than that of polished stainless steel but to be within the bounds given by Baumann et al., i.e. 0.25 < μW ≲ 0.5. An intermediate value of μW,2 = 0.4 is set for the surface friction coef- ficient of the bare steel tubes. In section 6.1 the influence of the wall friction coefficient on the simulation results is investigated revealing a moderate influence on the extent of the stagnant area whereas the rest of the flow field remains approximately unaffected. Table 2. Reference parameters for the simulations. Parameter Symbol Value Unit Particle diameter4 dp 0.6 mm Outlet/reference velocity uref 5.0 mm/s Density particles/solid phase ρs 2600 kg/m³ Tube diameter DT 27 mm Vertical tube spacing sV variable mm Horizontal tube spacing sH variable mm Tube wall friction coefficient (section 5.2.3) μW,1 0.25 - 3 Except for the Young´s modulus (see section 4.3.2). 4 For the geometry variations using the DEM simulations in section 6.1 a particle diameter of dp = 1.0 mm is used to reduce the computational cost. 40 4 Modelling Tube wall friction coefficient (section 5.1 and 5.2.1) μW,2 0.4 - Continuum model Density gas phase (air) ρf f(T) kg/m³ Specific heat gas phase (air) cp,f f(T) J/kgK Dynamic viscosity gas phase (air) μf 33 10-6 Pa s Total surface roughness (equation (20)) hR 1.0 µm Inner friction angle ϕi 34 ° Share of tube surface covered by particles (equation (20)) φ 0.8 - Loosest random packing (equation (14)) εmin 0.5 - Closest random packing (equation (14)) εmax 0.65 - Exponent in equation for pfr(εs) (equation (14)) r 2 - Exponent in equation for pfr(εs) (equation (14)) s 5 - Discrete particle model Youngs modulus of particles Ep 106 N/m² Friction coefficient between particles μpp 0.2 - Rolling friction coefficient μroll 0.3 - Coefficient of restitution e 0.1 - Surface friction coefficient at front wall (1) μfront,1 0.0 - Surface friction coefficient at front wall (2) μfront,2 0.5 - Surface friction coefficient at rear wall μrear 0.0 - 4.3.1 Parameters specific to the continuum model The material properties of the gas phase (air) in the continuum model are calculated according to tabular values in [104]. The surface roughness hR and the parameter φ5 are needed for the evaluation of the thermal contact resistance between bulk and tube surface (equation (18) and (20)). The 5 φ denotes the fraction of the surface covered by contacting particles (see section 4.1.2.1, equation (18) and (20)). 4 Modelling 41 value of φ = 0.8 is chosen according to Schlünder et al. [3]. The total surface roughness hR is the sum of the roughness of the particles and the surface. Senetakis measured hR,Sand ≈ 0.5 µm in experiments with Leighton Buzzard sand. The surface roughness of cold-rolled stainless steel as it is used in the validation experiment in section 5.2.3 is typically hR,surface ≈ 0.3…0.5 µm [105]. Therefore set the total surface roughness is set to hR = hR,Sand + hR,surface ≈ 1.0 µm. The values of the loosest and the closest random packing, εmin and εmax, as well as the two exponents, r and s, which are needed for the calculation of the granular pressure (equation (14)) between particles are set the same as in the work of Srivastava and Sun- deresan [63]. 4.3.2 Parameters specific to the discrete particle model The input parameters of the DEM-model include material parameters such as the Young´s modulus as well as material interaction parameters like friction coefficients. A general review on the calibration of DEM models is given by Coetzee [106]. The Young´s modulus of the particles is set to Ep = 106 N/m² which is about three or- ders of magnitude lower than values of typical materials such as for example sand stone [107]. Using a reduced Young’s modulus (reduced contact stiffness) is motivated by increas- ing the time step of the simulation which scales with the inverse square root of the con- tact stiffness (Δ𝑡 ~ 1/√𝑘𝑙 ) [108], leading to a substantial saving of computation time. This measure is very common in DEM models [106] and is justified as long as the nor- mal overlap between particles is less than 1 % of the particle radius [30]. In the current case the maximum normal forces between particles are FN,max ≈ 0.001 N which corre- sponds to a normal overlap of 0.6 % of the particle radius (see equation (26)). Further- more, to ensure independence of the simulation results on the contact stiffness, a varia- tion of the Young’s modulus was conducted (up to Ep = 108 N/m²) and virtually no im- pact was found (see appendix D.1). In contrast to the contact stiffness, the friction parameters are expected to be determin- ing for the simulation results. These are, on the one hand the static and dynamic friction coefficient according to which the tangential forces between particles are calculated (see equation (29), p. 36). On the other hand, there is the “rolling-friction” coefficient which is to account for the non-spherical shape of the real grains (see section 4.2.1). 42 4 Modelling Senetakis et al. [109] investigated the inter-particle coefficient of friction of Leighton Buzzard sand. They found that the coefficient of dynamic friction (μd) and the coeffi- cient of static friction (μs) are of very similar magnitude and they measured values of μs ≈ µd = μpp = 0.1…0.23. Given the inter-particle coefficient of friction μpp = 0.2, the rolling friction coefficient μroll was obtained from a dedicated laboratory experiment. In the experiment the static angle of repose ϕr ≈ 34° of the bulk was determined. Subsequently, the rolling friction coefficient was adjusted such that the measured value was attained in the simulation, which leads to μroll = 0.3. In appendix D.3 a brief sensitivity study is given of the influ- ence of the friction parameters on the simulation results. At the front and rear boundary of the geometry (see side view in Figure 11, p. 26) the surface friction coefficient is set to zero (μfront,1, μrear). This reflects an idealized flow section inside the bulk. However, the validation experiment (see section 3.1) captures the flow field visible at the acrylic glass front wall. Therefore, surface friction of the front wall needs to be considered. To this end, an annular shear cell is used to measure the wall friction coefficient be- tween quartz sand and an acrylic glass surface. The friction coefficient turns out to be not constant but to depend on the normal stress between surface and bulk. Hence, to determine the friction coefficient applicable to the PIV-measurement setup, the normal (horizontal) stress on the acrylic glass front wall of the mockup needs to be estimated. This is done by applying the well-known model developed by Janssen (e.g. in [17]), which yields a horizontal stress σH ≲ 2000 Pa. The measured friction coefficient for this normal stress is μfront,2 ≈ 0.5. Further details on the measurement of the surface friction coefficient are given in appendix A.4. 5 Validation 43 5 Validation A two-pronged strategy is adopted to validate the discrete particle model and the con- tinuum model. In a first step the granular flow field obtained from the discrete particle model is validated using PIV analysis. The flow field of the continuum model is subse- quently validated using the discrete particle model. In a second step the heat transfer model, which is only part of the continuum model, is validated based on the measurement data of the local heat transfer coefficient. 5.1 Discrete particle model For the validation of the discrete particle model the PIV measurement setup is used as described in section 3.1. The tube arrangement was TA3 (see Table 1) and an interme- diate flow speed was adjusted (uref = 4.0 mm/s). The specific focus is on the simulation of the flow field close to the acrylic glass front wall as this is the area which is captured by the experiment. Therefore, the surface friction between bulk and acrylic glass at the front wall is considered in the model (μfront,2 = 0.5, see section 4.3.2). Figure 14 shows contour plots of the granular flow field around a single tube obtained from the PIV measurement (left) and DEM simulation (right). The tube in the experi- mental plot is located in the center of the mockup (see framed area in Figure 5 (right), (p. 16)). Essentially, Figure 14 shows good agreement between simulation and experiment. However, while the simulated velocity profile is symmetric to the vertical center line of the tube, the experimental profile is slightly asymmetric. Higher velocities are observed on the left half of the plot. This asymmetry of the measured profile is due to disturbing effects in the experimental setup originating from the upstream row of tubes, from the side walls and from the tube chain conveyor (see also Appendix A.3). 44 5 Validation Figure 14. Contour plots